# Non-adiabatic Quantum Wavepacket Dynamics Simulation Based on Electronic
Structure Calculations using the Variational Quantum Eigensolver

Hirotoshi Hirai

Toyota Central Research and Development Labs., Inc.,
41-1, Yokomichi, Nagakute, Aichi 480-1192, Japan

Sho Koh

QunaSys Inc., Aqua Hakusan Building 9F,
1-13-7 Hakusan, Bunkyo, Tokyo 113-0001, Japan e-mail:
<EMAIL_ADDRESS>
###### Abstract

A non-adiabatic nuclear wavepacket dynamics simulation of the H2O+ de-
excitation process is performed based on electronic structure calculations
using the variational quantum eigensolver. The adiabatic potential energy
surfaces and non-adiabatic coupling vectors are computed with algorithms for
noisy intermediate-scale quantum devices, and time propagation is simulated
with conventional methods for classical computers. The results of non-
adiabatic transition dynamics from the $\tilde{B}$ state to $\tilde{A}$ state
reproduce the trend reported in previous studies, which suggests that this
quantum-classical hybrid scheme may be a useful application for noisy
intermediate-scale quantum devices.

##  1 Introduction

Elucidation of the non-adiabatic processes [1, 2] that appear in photochemical
reactions is important to understand the function of photocatalysts [3] and to
determine a visual mechanism [4]. The quantum wavepacket simulation is a
useful tool to investigate the mechanism of non-adiabatic processes. In the
adiabatic representation, it is necessary to compute the adiabatic potential
energy surfaces (PESs) and the non-adiabatic couplings (NACs) between each
pair of states for nuclear wavepacket simulations. Obtaining these physical
quantities in high accuracy in vicinity of the intersection (degenerate)
region, where the system has strong non-adiabaticity, requires the use of a
high-level quantum chemistry method, such as the full configuration
interaction (FCI) method. However, the computational cost increases
exponentially with the size of the system on conventional classical computers,
so that application is limited to small molecular systems [5].

Quantum computers are considered to be a powerful tool for quantum chemistry
computation because the computational cost can be suppressed into the
polynomial scale [6]. For a noise-less quantum computer, a so-called fault
tolerant quantum computer (FTQC), the quantum phase estimation algorithm (QPE)
[7, 8] has been proposed to compute the eigenenergies of electronic states.
Aside from the algorithms based on the adiabatic representation, an algorithm
has been proposed to directly solve the time evolution of the wavefunction of
nuclei and electrons[9]. However, the FTQC requires a sufficient number of
qubits for error correction; therefore, it has been estimated that it will
take 20 to 30 years to realize this [10].

Against FTQC, the variational eigenvalue solver (VQE) algorithm [11, 12] has
been reported in recent years to utilize current noisy intermediate-scale
quantum (NISQ) devices [13] that consist of a small number of noisy qubits and
allow only shallow quantum circuits. The VQE algorithm, which was originally
focused on the electronic ground state, has been expanded to compute the
energies of electronically excited states [14, 15, 16, 17, 18, 19]. In
addition, an algorithm used to compute NACs has also been proposed [20], so
that the components for quantum wavepacket simulations, PESs and NACs, can be
obtained using the VQE.

In the present study, the possibility that the NISQ devices can be applied to
quantum wavepacket simulations is investigated by non-adiabatic quantum
dynamics simulation in the adiabatic representation using a quantum-classical
hybrid scheme. The adiabatic PESs and NACs required in the simulations are
computed with the algorithms for NISQ devices using the simulator, and a non-
adiabatic quantum dynamics simulation is performed with a conventional
algorithm for classical computers [2]. We demonstrate the relaxation process
of H2O+ from $\tilde{B}^{2}B_{2}$ (the second excited state) to
$\tilde{A}^{2}A_{1}$ (the first excited state) can reproduce the trend
reported in the previous studies.

##  2 Method

###  2.1 Hamiltonian in the adiabatic representation

The degrees of freedom of electrons and nuclei must be treated quantum
mechanically to accurately describe the dynamics of non-adiabatic systems. The
Hamiltonian of the nuclear-electron system under non-relativity can be written
as

$H(\\{\vec{R}\\},\\{\vec{r}\\})=-\sum_{i=1}^{N}\frac{\hbar^{2}}{2M_{i}}\frac{\partial^{2}}{\partial\vec{R}_{i}^{2}}-\sum_{j=1}^{n}\frac{\hbar^{2}}{2m_{j}}\frac{\partial^{2}}{\partial\vec{r}_{j}^{2}}+V(\\{\vec{R}\\},\\{\vec{r}\\}),$
(1)

where $M_{i}$ and $m_{j}$ represent the mass of each nucleus and each
electron, and $N$ and $n$ represent the numbers of each particle,
respectively. $\\{\vec{R}\\}=\\{\vec{R}_{1},\vec{R}_{2},\cdots\vec{R}_{N}\\}$
and $\\{\vec{r}\\}=\\{\vec{r}_{1},\vec{r}_{2},\cdots\vec{r}_{n}\\}$ are the
collective representations for the positional coordinates of each nucleus and
each electron, respectively. $V(\\{\vec{R}\\},\\{\vec{r}\\})$ contains
coulombic interactions between electrons, between nuclei, and between
electrons and nuclei. The dynamics of the quantum system are described by the
following time-dependent Schrödinger equation,

$i\hbar\frac{\partial}{\partial
t}\Theta(\\{\vec{R}\\},\\{\vec{r}\\},t)=H\Theta(\\{\vec{R}\\},\\{\vec{r}\\},t),$
(2)

where the total wavefunction $\Theta(\\{\vec{R}\\},\\{\vec{r}\\},t)$ is
expanded with the eigenstates of the adiabatic electron Hamiltonian, which is
defined later,

$\Theta(\\{\vec{R}\\},\\{\vec{r}\\},t)=\sum_{l}\chi_{l}(\\{\vec{R}\\},t)\psi_{l}(\\{\vec{r}\\};\\{\vec{R}\\}).$
(3)

Such expansion is possible because $\psi_{l}(\\{\vec{r}\\};\\{\vec{R}\\})$ is
a complete set for the system. Here, $\psi_{l}(\\{\vec{r}\\};\\{\vec{R}\\})$
is the eigenstates of the adiabatic electron Hamiltonian,

$H_{el}(\\{\vec{r}\\};\\{\vec{R}\\})\psi_{l}(\\{\vec{r}\\};\\{\vec{R}\\})=E_{l}(\\{\vec{R}\\})\psi_{l}(\\{\vec{r}\\};\\{\vec{R}\\})$
(4)

where $H_{el}(\\{\vec{r}\\};\\{\vec{R}\\})$ is defined as

$H_{el}(\\{\vec{r}\\};\\{\vec{R}\\})=-\sum_{j=1}^{n}\frac{\hbar^{2}}{2m_{j}}\frac{\partial^{2}}{\partial\vec{r}_{j}^{2}}+V(\\{\vec{r}\\};\\{\vec{R}\\}).$
(5)

It should be noted that the adiabatic electron Hamiltonian
$H_{el}(\\{\vec{r}\\};\\{\vec{R}\\})$ and its eigenvalues and eigenstates are
parametrically dependent on $\\{\vec{R}\\}$. Here, the expansion coefficient
$\chi_{l}(\\{\vec{R}\\},t)$ can be interpreted as a nuclear wavefunction that
represents the chemical reactions, and $E_{l}(\\{\vec{R}\\})$ gives the PES of
the $l$-th electronic excited state.

Using these notations, the previous time-dependent Schrödinger equation can be
rewritten as

$\begin{split}i\hbar\frac{\partial}{\partial
t}\chi_{p}(\\{\vec{R}\\},t)&=[-\sum_{i}^{N}\frac{\hbar^{2}}{2M_{i}}\frac{\partial^{2}}{\partial\vec{R}_{i}^{2}}+E_{p}(\\{\vec{R}\\})]\chi_{p}(\\{\vec{R}\\},t)\\\
&-\sum_{q}\sum_{i}\frac{\hbar^{2}}{M_{i}}\bra{\psi_{p}}\frac{\partial}{\partial\vec{R}_{i}}\ket{\psi_{q}}\frac{\partial}{\partial\vec{R}_{i}}\chi_{q}(\\{\vec{R}\\},t)\\\
&-\sum_{q}\sum_{i}\frac{\hbar^{2}}{2M_{i}}\bra{\psi_{p}}\frac{\partial^{2}}{\partial\vec{R}_{i}^{2}}\ket{\psi_{q}}\chi_{q}(\\{\vec{R}\\},t),\end{split}$
(6)

where $\bra{\psi_{p}}\frac{\partial}{\partial\vec{R}_{i}}\ket{\psi_{q}}$ is
the first-order NAC, and this term is now referred to as the NAC vector.
$\bra{\psi_{p}}\frac{\partial^{2}}{\partial\vec{R}_{i}^{2}}\ket{\psi_{q}}$ is
the second-order NAC coefficient. The diagonal components of the second-order
NACs give non-adiabatic corrections to the adiabatic PESs. On the other hand,
the off-diagonal components of the NAC vectors and second-order NAC
coefficients have the effect of mixing the adiabatic states. It should be
noted that the first-order NAC term is proportional to the inner product of
the NAC vector and the velocity of the nuclear wavepacket, and the non-
adiabatic transitions occur when the both are large. If these non-adiabatic
terms are negligibly small, then the time-dependent Schrödinger equation under
the Born-Oppenheimer approximation is obtained as

$i\hbar\frac{\partial}{\partial
t}\chi_{p}(\\{\vec{R}\\},t)=[-\sum_{i}^{N}\frac{\hbar^{2}}{2M_{i}}\frac{\partial^{2}}{\partial\vec{R}_{i}^{2}}+E_{p}(\\{\vec{R}\\})]\chi_{p}(\\{\vec{R}\\},t).$
(7)

This equation shows that the motions of nuclei and electrons can be separated,
and the wavepacket $\chi_{p}(\\{\vec{R}\\},t)$ moves on only one adiabatic
PES, $E_{p}(\\{\vec{R}\\})$. This approximation assumes that the non-adiabatic
terms are negligibly small because the nuclear mass $\\{M_{i}\\}$ is greater
than the electron mass, $m$. This approximation holds well near the stable
structure of the ground state and is the starting point for the theory of
first-principles calculations and quantum chemistry calculations for various
materials and molecules. However, this approximation breaks down in various
phenomena involving excited states. This can be understood as follows. The NAC
vectors can be rewritten as

$\bra{\psi_{p}}\frac{\partial}{\partial\vec{R}_{i}}\ket{\psi_{q}}=\frac{\bra{\psi_{p}}\frac{\partial
H_{el}}{\partial\vec{R}_{i}}\ket{\psi_{q}}}{E_{q}-E_{p}},$ (8)

which shows that the NAC vector is inversely proportional to the gap between
adiabatic potential energies. The adiabatic approximation (Born-Oppenheimer
approximation) thus breaks down near the conical intersection where the
adiabatic potentials become close to each other. For this reason, the NAC
vectors cannot be ignored in order to describe the non-adiabatic transition,
in which the order of electronic states changes, in the de-excitation process,
as often is the case in photochemical reactions. Against the first-order NAC,
the contribution of the second-order NAC can be considered negligible because
the contribution of the NAC vector dominates in the semiclassical limit,
$\hbar\rightarrow 0$ [1, 2]. Thereby the second-order NAC is ignored in this
study.

###  2.2 H2O+ molecule (symmetric extension model, light-heavy-light
approximation)

In this study, the relaxation process of the H2O+ molecule from
$\tilde{B}^{2}B_{2}$ (the second excited state) to $\tilde{A}^{2}A_{1}$ (the
first excited state) is studied to demonstrate the proposed method. The
process of ionization of water molecules and subsequent decomposition
(generation of H+, OH+, O+, etc.) is important for investigation of the
universe (interstellar medium, molecular clouds, planetary atmosphere, in
comets) and radiation damage (radiation therapy for cancer) in living
organisms [21]. Therefore, it has been investigated from both experimental and
theoretical perspectives [22, 23, 24, 25, 26], and is suitable as a subject
for verification of the proposed method.

There are three nuclei in H2O+, and the Hamiltonian of a three-body system is
given in the literature [27] for the case of H2O, where the degrees of freedom
for rotation and the coupling between rotation and vibration (Coriolis force)
are ignored, and only the intramolecular vibration is focused on.

For further simplification, we assume a symmetric extension model,
$r_{1}=r_{2}=r$, and use a light-heavy-light (LHL) approximation (i.e.,
$M_{O}\rightarrow\infty$ from $M_{H}\ll M_{O}$). The kinetic energy term of
the Hamiltonian for H2O+ can be written as follows.

$T=-\frac{1}{M_{H}}\frac{\partial^{2}}{\partial
r^{2}}-\frac{1}{M_{H}r^{2}}\frac{\partial^{2}}{\partial\theta^{2}}-\frac{1}{4M_{H}r^{2}}(1+csc^{2}(\theta)),$
(9)

where $r=r_{1}=r_{2}$ represents the length of OH, $\theta$ represents the
angle of H-O-H, and $csc(\theta)=1/sin(\theta)$. The third term is small
(0.0005 Ha or less) in the range dealt with in this study
($\theta>90^{\circ}$) and is thus ignored in this work. The above
approximations are justified in the relaxation process immediately after
excitation and before energy transfer to other vibrational modes occurs [28].

The transition to $\tilde{a}^{4}B_{1}$ state is also known to occur, either
from $\tilde{B}^{2}B_{2}$ or $\tilde{A}^{2}A_{1}$ states by spin-orbit
couplings. The previous theoretical study [21] and experimental result [29]
indicated that this had a minor effect on the de-excitation process due to the
small magnitude of the spin-orbit coupling term. Therefore, the effect of the
spin-orbit coupling is omitted in this study.

###  2.3 Computation of PESs and NAC vectors

The H2O+ molecule is set on the y-z plane where the O atom is placed on the
origin and the C2v axis is set along the z-axis, as shown in FIG. 1. The three
lowest electronic states of $\tilde{X}$, $\tilde{A}$ and $\tilde{B}$ are
calculated using the subspace-search variational quantum eigensolver (SSVQE)
method [14], which is an expansion of VQE for obtaining excited states, with
the quantum algorithm simulator Qulacs [30] and OpenFermion [31]. The three
minimum molecular orbitals $1b_{2}$, $3a_{1}$ and $1b_{1}$ are taken into the
active space to describe the conical intersection structure between the
$\tilde{A}$ and $\tilde{B}$ states. The molecular orbitals used for the second
quantized Hamiltonian are optimized using the complete active space self-
consistent field (CASSCF) method with the basis set of 6-31G(d,p) implemented
in the PySCF program [32].

We use Jordan-Wigner mapping [33] to map spin orbitals onto the qubits so that
different spins alternate on 6 qubits. In the SSVQE calculation, the three
initial states of $\ket{101111}$, $\ket{111011}$ and $\ket{1111110}$ on
quantum devices are prepared to correspond to the electron configurations with
one hole in $1b_{1}$, $3a_{1}$ and $1b_{1}$, respectively. Symmetry preserving
real ansatz [34] $U(\phi)$ is used to represent the electronic eigenstates of
the second quantized Hamiltonian as shown in FIG. 2. The weighted energy
summation of the three electronic states are given by

$E_{SSVQE}(\phi)=w_{0}E_{0}(\phi)+w_{1}E_{1}(\phi)+w_{2}E_{2}(\phi),$ (10)

where the weights for the electronic ground state, the first excited state and
the second excited state are given as $w_{0}=9.0,w_{1}=4.0,$ and $w_{2}=1.0$,
respectively, and the eigenenergies of each electronic state are calculated as

$\displaystyle E_{0}(\phi)$ $\displaystyle=$
$\displaystyle\bra{101111}U^{\dagger}(\phi)H_{el}U(\phi)\ket{101111},$ (11)
$\displaystyle E_{1}(\phi)$ $\displaystyle=$
$\displaystyle\bra{111011}U^{\dagger}(\phi)H_{el}U(\phi)\ket{111011},$ (12)
$\displaystyle E_{2}(\phi)$ $\displaystyle=$
$\displaystyle\bra{111110}U^{\dagger}(\phi)H_{el}U(\phi)\ket{111110}.$ (13)

The parameters $\phi$ are optimized using the BFGS optimizer in the SSVQE
calculation. Since only the concept of the method is presented in this study,
we do not consider the effect of shot noise. The obtained adiabatic PESs of
the $\tilde{B}^{2}B_{2}$ and $\tilde{A}^{2}A_{1}$ states are shown in FIG. 3.

After optimization of the electronic states, we compute the NAC matrix
elements given in Eq. (8). The operators $\frac{\partial
H_{el}}{\partial\vec{R}}$ are calculated using the finite numerical derivative
method according to

$\frac{\partial
H_{el}(\vec{R}_{i})}{\partial\vec{R}_{i}}=\frac{H_{el}(\vec{R}_{i}+\Delta
R\vec{e_{i}})-H_{el}(\vec{R}_{i}-\Delta R\vec{e_{i}})}{2\Delta R},$ (14)

where $\Delta R$ is set to be 0.001 bohr and $\vec{e_{i}}$ is the unit vector
along $\vec{R}_{i}$. The transition amplitude $\bra{\psi_{p}}\frac{\partial
H_{el}(R)}{\partial\vec{R}}\ket{\psi_{q}}$ can be evaluated on quantum
computers, so that the NAC vectors can be computed by dividing this value by
the energy gap, $E_{q}-E_{p}$, as reported in the previous study [20].

The conversion from the $Y$ and $Z$ NAC vector components to the $r$ and
$\theta$ components can be calculated as follows:

$\bra{\Psi_{p}}\frac{\partial}{\partial
r_{1}}\ket{\Psi_{q}}=-\bra{\Psi_{p}}\frac{\partial}{\partial
Y_{1}}\ket{\Psi_{q}}\sin\frac{\theta}{2}-\bra{\Psi_{p}}\frac{\partial}{\partial
Z_{1}}\ket{\Psi_{q}}\cos\frac{\theta}{2},$ (15)
$\bra{\Psi_{p}}\frac{\partial}{\partial
r_{2}}\ket{\Psi_{q}}=\bra{\Psi_{p}}\frac{\partial}{\partial
Y_{2}}\ket{\Psi_{q}}\sin\frac{\theta}{2}-\bra{\Psi_{p}}\frac{\partial}{\partial
Z_{2}}\ket{\Psi_{q}}\cos\frac{\theta}{2},$ (16)
$\bra{\Psi_{p}}\frac{\partial}{\partial\theta_{1}}\ket{\Psi_{q}}=-\bra{\Psi_{p}}\frac{\partial}{\partial
Y_{1}}\ket{\Psi_{q}}\cos\frac{\theta}{2}+\bra{\Psi_{p}}\frac{\partial}{\partial
Z_{1}}\ket{\Psi_{q}}\sin\frac{\theta}{2},$ (17)
$\bra{\Psi_{p}}\frac{\partial}{\partial\theta_{2}}\ket{\Psi}_{q}=\bra{\Psi_{p}}\frac{\partial}{\partial
Y_{2}}\ket{\Psi_{q}}\cos\frac{\theta}{2}+\bra{\Psi_{p}}\frac{\partial}{\partial
Z_{2}}\ket{\Psi_{q}}\sin\frac{\theta}{2}.$ (18)

And the following relationships hold:

$\bra{\Psi_{p}}\frac{\partial}{\partial
r}\ket{\Psi_{q}}=\bra{\Psi_{p}}\frac{\partial}{\partial
r_{1}}\ket{\Psi_{q}}=\bra{\Psi_{p}}\frac{\partial}{\partial
r_{2}}\ket{\Psi_{q}},$ (19)
$\bra{\Psi_{p}}\frac{\partial}{\partial\theta}\ket{\Psi_{q}}=\bra{\Psi_{p}}\frac{\partial}{\partial\theta_{1}}\ket{\Psi_{q}}=\bra{\Psi_{p}}\frac{\partial}{\partial\theta_{2}}\ket{\Psi_{q}}.$
(20)

Figure 1: Coordinates for adiabatic PESs and NACs computations. The H2O+
molecule is set on the y-z plane where the O atom is placed on the origin and
the C2v axis is set along the z-axis. Figure 2: Quantum circuit of the
symmetry preserving ansatz. $R_{y}(\phi)=\exp(-i\phi Y/2)$, denotes the y-axis
rotation gate and D $=5$ denotes the depth of the circuit. The rotation angles
implemented in two-qubit unitary gates $U(\phi)$ are optimized in VQE
calculation. Figure 3: Adiabatic PESs of the $\tilde{B}^{2}B_{2}$ (green) and
$\tilde{A}^{2}A_{1}$ (violet) states for the H2O+ molecule. The corresponding
contour plots are shown at bottom panel of the figure.

(a) NAC vector for the $r$-direction (b) NAC vector for the $\theta$-direction

(c) NAC vector for the $r$-direction (CASCF) (d) NAC vector for the
$\theta$-direction (CASSCF)

Figure 4: NAC vector’s components for the (a)$r$\- and (b)$\theta$-directions
between the $\tilde{B}^{2}B_{2}$ and $\tilde{A}^{2}A_{1}$ states for the H2O+
molecule. For comparison, NAC vector’s components computed by CASSCF are shown
in (c) and (d) for each direction, respectively.

Figure. 4 shows NAC vector components for the $r$\- and $\theta$-directions
between the $\tilde{B}^{2}B_{2}$ and $\tilde{A}^{2}A_{1}$ states for the H2O+
molecule. For comparison, the corresponding components computed by CASSCF are
shown. In the CASSCF computation for NAC vectors, OpenMolcas [35] is used and
the computational conditions such as the basis set and the active space are
set to the same ones as SSVQE computations except the derivative method (the
analytical derivative method is used for CASSCF). Figure. 4 shows that there
are large NAC vector components on the conical intersection seam between the
$\tilde{B}^{2}B_{2}$ and $\tilde{A}^{2}A_{1}$ states shown in FIG. 3 for the
both results of SSVQE and CASSCF methods. There is a slight difference between
the NAC vector components of SSVQE and CASSCF, which can be attributed to the
error in the finite numerical differentiation of SSVQE.

###  2.4 Grid notation of the Hamiltonian

A non-adiabatic quantum dynamics simulation was conducted on a grid in real
space. The adiabatic PESs and NAC vectors were computed with the method
explained in the previous section. Since these values were computed on the
coarse grid, these values should be complemented for a fine grid to conduct
the quantum dynamics. For completion, the three-dimensional spline method is
adopted. The second-order central difference method is used to represent the
kinetic energy term,

$\frac{\partial^{2}}{\partial
r^{2}}\chi(r,\theta,t)=\frac{\chi_{i}(r+dr,\theta,t)-2\chi(r,\theta,t)+\chi(r-dr,\theta,t)}{dr^{2}}+O(dr^{2}).$
(21)

The same applies to the $y$-direction. The second-order central difference
method is also used for the time derivative,

$\frac{\partial}{\partial
t}\chi(r,\theta,t)=\frac{\chi_{i}(r,\theta,t+dt)-\chi(r,\theta,t-dt)}{2dt}+O(dt^{2}).$
(22)

Therefore, the time evolution of the system can be computed by the following
formula:

$\chi(r,\theta,t+dt)=\chi(r,\theta,t-dt)-2dti\hbar H\chi(r,\theta,t),$ (23)

where $H$ is the non-adiabatic Hamiltonian with an adiabatic representation
under the symmetric extension model and the LHL approximation. In this study,
$dx=0.0443$ bohr and $dy=0.0409$ rad were used, and the ranges were set to be
0.9449-3.7352 bohr for $x$ and 0.5236-3.1007 rad for $y$. The grid size is
then $(N_{r},N_{\theta})=(64,64)$.

##  3 Results and discussion

###  3.1 Preparation of initial wavepacket

In this study, a vertical transition was assumed according to the Franck-
Condon principle [36]. This approximation is based on the instantaneous
occurrence of electronic transitions compared to the time scale of nuclear
motion, so that the ionization occurs without changing the position of the
nuclei. At first, the nuclear wavepacket of the zero-point vibrational state
in the electronic ground state PES was computed, and the initial wavepacket
was constructed by placing it on the adiabatic PES of the second excited state
($\tilde{B}^{2}B_{2}$). The electron wavepacket in the electronic ground state
was obtained by diagonalization of the adiabatic Hamiltonian. The sum of the
zero-point vibrational energies of the $r$ and $\theta$ degrees of freedom
known in the experiment should match the energy of this wavepacket. The energy
of the wavepacket and the energies obtained by the photoelectron spectroscopy
experiment [37] are shown in TABLE. 1, which also shows the results for D2O+.
The results obtained in this study are in good agreement with the experimental
values, including the isotopic effect.

Table 1: Comparison of zero-point vibrational energies [hartree]. System | This work | Experiment[37]
---|---|---
H2O+ | 0.0210 | 0.0214
D2O+ | 0.0149 | 0.0157

###  3.2 De-excitation process from $\tilde{B}^{2}B_{2}$ to
$\tilde{A}^{2}A_{1}$ in the H2O+ molecule

Figure. 5 shows the time-evolution of the populations for the
$\tilde{B}^{2}B_{2}$ and $\tilde{A}^{2}A_{1}$ states. For comparison, the
results simulated with the PESs and NAC vectors computed by CASSCF method are
also shown.

Figure 5: Time-evolution of the populations of the $\tilde{B}^{2}B_{2}$ and
$\tilde{A}^{2}A_{1}$ states and the time integrated value of the flux that
reached $\theta=\pi$ on $\tilde{A}^{2}A_{1}$ for the H2O+ molecule based on
the PESs and NAC vectors computed by SSVQE method (solid lines) and CASSCF
method (dashed lines).

The both results are qualitatively consistent, but there is a slight
difference in the curve of the population as a function of time, reflecting
the difference of NAC vector of SSVQE and CASSCF as shown in Fig. 4. The
transition from $\tilde{B}^{2}B_{2}$ to $\tilde{A}^{2}A_{1}$ can be confirmed
approximately 10 fs after the start of the simulation for the both results.
Figure. 6 shows a snapshot of the time-evolution of the wavepacket at each
time.

(a) $t=0$ fs on $\tilde{B}^{2}B_{2}$ (b) $t=0$ fs on $\tilde{A}^{2}A_{1}$ (c)
$t=2.4$ fs on $\tilde{B}^{2}B_{2}$ (d) $t=2.4$ fs on $\tilde{A}^{2}A_{1}$ (e)
$t=4.8$ fs on $\tilde{B}^{2}B_{2}$ (f) $t=4.8$ fs on $\tilde{A}^{2}A_{1}$ (g)
$t=8.4$ fs on $\tilde{B}^{2}B_{2}$ (h) $t=8.4$ fs on $\tilde{A}^{2}A_{1}$

Figure 6: Snapshots of the wavepackets on the $\tilde{B}^{2}B_{2}$ and
$\tilde{A}^{2}A_{1}$ states for H2O+ molecule.

The wavepacket moves in the $r$-positive direction from $t=0$ to $t=2.4$ fs,
then turns around and moves in the $r$-negative direction from $t=2.4$ to
$t=4.8$ fs. The wavepacket reaches the conical intersection near $t=4.8$ fs,
and a part of the wavepacket transits from $\tilde{B}^{2}B_{2}$ to
$\tilde{A}^{2}A_{1}$. However, the NAC vector in the $r$ direction is not so
large at this position; therefore, only a part of the transition occurs,
whereas most of the wavepacket transitioned to $\tilde{A}^{2}A_{1}$ in the
previous study [21]. The reason for this difference can be considered to be
the small active space used to describe the conical intersection, so that the
rate of transition from $\tilde{B}^{2}B_{2}$ to $\tilde{A}^{2}A_{1}$ is
smaller than that of the previous report. The $\theta$ positive motion is
subsequently promoted and almost all the wavepackets transition to
$\tilde{A}^{2}A_{1}$ when passing through the conical intersection. This is
due to the large NAC vector in the $\theta$ direction at this position. The
wavepacket that transitions to $\tilde{A}^{2}A_{1}$ has momentum in the
$\theta$-positive direction and reaches to $\theta=\pi$, where it is absorbed
at the boundary by the absorption potential. This wavepacket is considered to
make the transition to $\tilde{X}^{2}B_{1}$ (electronic ground state) by
Renner-Teller coupling, which has large coupling near $\theta=\pi$ (the linear
geometry). The motion of the nuclear wavepacket is substantially consistent
with the previous reports [21, 29].

##  4 Conclusion

In this study, we performed a non-adiabatic quantum dynamics simulation in the
adiabatic representation using a quantum-classical hybrid scheme. The non-
adiabatic quantum dynamics computations were conducted with the conventional
algorithm for classical computers using the Hamiltonian components in the
adiabatic representation, which consisted of PESs and NACs computed with NISQ
algorithms. The simulation results for the de-excitation process from
$\tilde{B}^{2}B_{2}$ (second excited state) to $\tilde{A}^{2}A_{1}$ (first
excited state) of the H2O+ molecule reproduced the trend in the previous study
and consistent with the the results based on the PESs and NAC vectors computed
with the CASSCF method, which indicates the effectiveness of this method. The
quantum-classical hybrid calculation scheme proposed in this study may be a
promising method for non-adiabatic quantum dynamics simulations in the NISQ
era.

## Acknowledgement

The authors thank Dr. Y. O. Nakagawa for his valuable comments.

## CRediT authorship contribution statement

Hirotoshi Hirai: Conceptualization, Methodology, Software, Writing - original
draft, Investigation, Visualization. Sho Koh: Methodology, Software, Writing -
review & editing, Investigation.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] Michael Baer.  Introduction to the theory of electronic non-adiabatic coupling terms in molecular systems.  Phys. Rep., 358(2):75–142, 2002.
  * [2] Graham A Worth and Lorenz S Cederbaum.  Beyond born-oppenheimer: molecular dynamics through a conical intersection.  Annu. Rev. Phys. Chem., 55:127–158, 2004.
  * [3] Nick Serpone.  Photocatalysis.  Kirk-Othmer Encyclopedia of Chemical Technology, 2000.
  * [4] Gunnar Svaetichin and Edward F MacNichol Jr.  Retinal mechanisms for chromatic and achromatic vision.  Ann. N. Y. Acad. Sci., 74(2):385–404, 1958.
  * [5] Jeppe Olsen, Poul Jørgensen, and Jack Simons.  Passing the one-billion limit in full configuration-interaction (fci) calculations.  Chemical Physics Letters, 169(6):463–472, 1990.
  * [6] Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, Simon C Benjamin, and Xiao Yuan.  Quantum computational chemistry.  Rev. Mod. Phys., 92(1):015003, 2020.
  * [7] A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv quant-ph/9511026, 1995.
  * [8] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca.  Quantum algorithms revisited.  Proc. R. Soc. A, 454(1969):339–354, 1998.
  * [9] Ivan Kassal, Stephen P Jordan, Peter J Love, Masoud Mohseni, and Alán Aspuru-Guzik.  Polynomial-time quantum algorithm for the simulation of chemical dynamics.  Proc. Natl. Acad. Sci., 105(48):18681–18686, 2008.
  * [10] Matt Langione, Corban Tillemann-Dick, Amit Kumar, and Vikas Taneja.  Where will quantum computers create value—and when.  Boston Consulting Group, November 2019, 2019.
  * [11] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Comm., 5(1):4213, 2014.
  * [12] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, 2017.
  * [13] John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, 2018.
  * [14] Ken M Nakanishi, Kosuke Mitarai, and Keisuke Fujii.  Subspace-search variational quantum eigensolver for excited states.  Phys. Rev. Res., 1(3):033062, 2019.
  * [15] Robert M. Parrish, Edward G. Hohenstein, Peter L. McMahon, and Todd J. Martínez.  Quantum computation of electronic transitions using a variational quantum eigensolver.  Phys. Rev. Lett., 122:230401, Jun 2019.
  * [16] Oscar Higgott, Daochen Wang, and Stephen Brierley.  Variational Quantum Computation of Excited States.  Quantum, 3:156, July 2019.
  * [17] Tyson Jones, Suguru Endo, Sam McArdle, Xiao Yuan, and Simon C. Benjamin.  Variational quantum algorithms for discovering hamiltonian spectra.  Phys. Rev. A, 99:062304, Jun 2019.
  * [18] Jarrod R. McClean, Mollie E. Kimchi-Schwartz, Jonathan Carter, and Wibe A. de Jong.  Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states.  Phys. Rev. A, 95:042308, Apr 2017.
  * [19] Pauline J. Ollitrault, Abhinav Kandala, Chun-Fu Chen, Panagiotis Kl. Barkoutsos, Antonio Mezzacapo, Marco Pistoia, Sarah Sheldon, Stefan Woerner, Jay M. Gambetta, and Ivano Tavernelli.  Quantum equation of motion for computing molecular excitation energies on a noisy quantum processor.  Phys. Rev. Res., 2:043140, Oct 2020.
  * [20] Shiro Tamiya, Sho Koh, and Yuya O. Nakagawa.  Calculating nonadiabatic couplings and berry’s phase by variational quantum eigensolvers.  Phys. Rev. Research, 3:023244, Jun 2021.
  * [21] Jaime Suárez, L Méndez, and I Rabadán.  Nonadiabatic quantum dynamics predissociation of h2o+ ($\tilde{B}^{2}$b2).  J. Phys. Chem. Lett., 6(1):72–76, 2015.
  * [22] KH Tan, CE Brion, Ph E Van der Leeuw, and MJ Van der Wiel.  Absolute oscillator strengths (10–60 ev) for the photoabsorption, photoionisation and fragmentation of h20s.  Chem. Phys., 29(3):299–309, 1978.
  * [23] Udo Werner, K Beckord, J Becker, and Hans Otto Lutz.  3d imaging of the collision-induced coulomb fragmentation of water molecules.  Phys. Rev. Lett., 74(11):1962, 1995.
  * [24] F Gobet, B Farizon, M Farizon, MJ Gaillard, M Carré, M Lezius, P Scheier, and TD Märk.  Total, partial, and electron-capture cross sections for ionization of water vapor by 20–150 kev protons.  Phys. Rev. Lett., 86(17):3751, 2001.
  * [25] H Luna and EC Montenegro.  Fragmentation of water by heavy ions.  Phys. Rev. Lett., 94(4):043201, 2005.
  * [26] H Luna, ALF De Barros, JA Wyer, SWJ Scully, Julien Lecointre, PMY Garcia, GM Sigaud, ACF Santos, V Senthil, MB Shah, et al.  Water-molecule dissociation by proton and hydrogen impact.  Phys. Rev. A, 75(4):042711, 2007.
  * [27] Bruce R Johnson and William P Reinhardt.  Adiabatic separations of stretching and bending vibrations: Application to h2o.  J. Chem. Phys., 85(8):4538–4556, 1986.
  * [28] Geert-Jan Kroes.  Photodissociation of ch2: A test case for the light–heavy–light approximation.  J. Chem. Phys., 101(7):5792–5803, 1994.
  * [29] CR Brundle and David Warren Turner.  High resolution molecular photoelectron spectroscopy ii. water and deuterium oxide.  Proc. R. Soc. A, 307(1488):27–36, 1968.
  * [30] Yasunari Suzuki, Yoshiaki Kawase, Yuya Masumura, Yuria Hiraga, Masahiro Nakadai, Jiabao Chen, Ken M. Nakanishi, Kosuke Mitarai, Ryosuke Imai, Shiro Tamiya, Takahiro Yamamoto, Tennin Yan, Toru Kawakubo, Yuya O. Nakagawa, Yohei Ibe, Youyuan Zhang, Hirotsugu Yamashita, Hikaru Yoshimura, Akihiro Hayashi, and Keisuke Fujii.  Qulacs: a fast and versatile quantum circuit simulator for research purpose.  arXiv:2011.13524, 2020.
  * [31] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby, Pranav Gokhale, Thomas Häner, Tarini Hardikar, Vojtěch Havlíček, Oscar Higgott, Cupjin Huang, Josh Izaac, Zhang Jiang, Xinle Liu, Sam McArdle, Matthew Neeley, Thomas O’Brien, Bryan O’Gorman, Isil Ozfidan, Maxwell D Radin, Jhonathan Romero, Nicolas P D Sawaya, Bruno Senjean, Kanav Setia, Sukin Sim, Damian S Steiger, Mark Steudtner, Qiming Sun, Wei Sun, Daochen Wang, Fang Zhang, and Ryan Babbush.  OpenFermion: the electronic structure package for quantum computers.  Quantum Sci. Technol., 5(3):034014, jun 2020.
  * [32] Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Marc Barbry, Nick S. Blunt, Nikolay A. Bogdanov, George H. Booth, Jia Chen, Zhi-Hao Cui, Janus J. Eriksen, Yang Gao, Sheng Guo, Jan Hermann, Matthew R. Hermes, Kevin Koh, Peter Koval, Susi Lehtola, Zhendong Li, Junzi Liu, Narbe Mardirossian, James D. McClain, Mario Motta, Bastien Mussard, Hung Q. Pham, Artem Pulkin, Wirawan Purwanto, Paul J. Robinson, Enrico Ronca, Elvira R. Sayfutyarova, Maximilian Scheurer, Henry F. Schurkus, James E. T. Smith, Chong Sun, Shi-Ning Sun, Shiv Upadhyay, Lucas K. Wagner, Xiao Wang, Alec White, James Daniel Whitfield, Mark J. Williamson, Sebastian Wouters, Jun Yang, Jason M. Yu, Tianyu Zhu, Timothy C. Berkelbach, Sandeep Sharma, Alexander Yu. Sokolov, and Garnet Kin-Lic Chan.  Recent developments in the pyscf program package.  J. Chem. Phys., 153(2):024109, 2020.
  * [33] Alán Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon.  Simulated quantum computation of molecular energies.  Science, 309(5741):1704–1707, 2005.
  * [34] Yohei Ibe, Yuya O Nakagawa, Nathan Earnest, Takahiro Yamamoto, Kosuke Mitarai, Qi Gao, and Takao Kobayashi.  Calculating transition amplitudes by variational quantum deflation.  arXiv:2002.11724, 2020.
  * [35] Francesco Aquilante, Jochen Autschbach, Alberto Baiardi, Stefano Battaglia, Veniamin A Borin, Liviu F Chibotaru, Irene Conti, Luca De Vico, Mickaël Delcey, Ignacio Fdez. Galván, et al.  Modern quantum chemistry with [open] molcas.  J. Chem. Phys., 152(21):214117, 2020.
  * [36] James Franck and EG Dymond.  Elementary processes of photochemical reactions.  Trans. Faraday Soc., 21(February):536–542, 1926.
  * [37] JE Reutt, LS Wang, YT Lee, and DA Shirley.  Molecular beam photoelectron spectroscopy and femtosecond intramolecular dynamics of h2o+ and d2o+.  J. Chem. Phys., 85(12):6928–6939, 1986.