Quantum Simulation of Conical Intersections

We explore the simulation of conical intersections (CIs) on quantum devices, setting the groundwork for potential applications in nonadiabatic quantum dynamics within molecular systems. The intersecting potential energy surfaces of H$_{3}^{+}$ are computed from a variance-based contracted quantum eigensolver. We show how the CIs can be correctly described on quantum devices using wavefunctions generated by the anti-Hermitian contracted Schr{\"o}dinger equation ansatz, which is a unitary transformation of wavefunctions that preserves the topography of CIs. A hybrid quantum-classical procedure is used to locate the seam of CIs. Additionally, we discuss the quantum implementation of the adiabatic to diabatic transformation and its relation to the geometric phase effect. Results on noisy intermediate-scale quantum devices showcase the potential of quantum computers in dealing with problems in nonadiabatic chemistry.


I. INTRODUCTION
Nonadiabatic processes involve nuclear motion on multiple potential energy surfaces (PESs).These processes are ubiquitous in nature and have been studied extensively in diverse areas such as spectroscopy, solar energy conversion, chemiluminescence, photosynthesis, and photostability of biomolecules.[1][2][3][4][5][6][7][8][9][10][11][12][13][14] Different potential energy surfaces can intersect at regions that exhibit a conical-shaped topography, known as conical intersections (CIs).[15][16][17][18] In the vicinity of CIs, the Born-Oppenheimer approximation which assumes adiabaticity breaks down.Systems with nonadiabaticity can undergo sudden changes in their dominant configurations at CIs, leading to the classical and well-known "hop" picture between different electronic states.[19] CIs act as highly efficient channels for converting external excitation energy, usually carried by a photon, to internal electronic energy.Their characterization is crucial to understand rich photochemistry and photobiology processes involving energy conversion.
CIs are in general difficult to treat with quantum mechanical methods for several reasons.First, from the perspective of electronic structure theory, the excited states are harder to compute than the ground state as they correspond to first-order critical points rather than the global minimum.Moreover, the nonunitary ansatz for the wavefunction employed in some methods, such as standard coupled cluster (CC) methods, gives an incorrect topography of CIs.[20][21][22][23] Second, since most electronic structure programs work under the Born-Oppenheimer approximation, results obtained from these programs are not readily applicable for subsequent chemical dynamics studies, especially near CIs.The process of converting the original adiabatic electronic structure data to a diabatic representation, referred as diabatization, is an active yet non-unified field due to the non-uniqueness of quasi-diabatic representations.[24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] Third, the dynamics of nonadiabatic systems typically require a more complex treatment than the dynamics on a single potential energy surface.For example, we need to expand wavefunctions in the basis of every diabatic state to account for effective state transition in quantum dynamics.[1] Quantum computers could be a natural solution for nonadiabatic chemistry.[39][40][41][42][43] To address some of the concerns in the last paragraph, we observe first that the gate operations are unitary, which makes it convenient to implement a unitary ansatz of wavefunctions (e.g., the unitary coupled cluster (UCC) ansatz [44,45] or the anti-Hermitian contracted Schrödinger equation (ACSE) ansatz [46][47][48][49][50][51][52][53][54]).They offer robust and accurate solutions to the electronic structure data near CIs.In fact, for the ACSE ansatz used in this paper, classical calculations of CIs are well established.[55][56][57][58] Second, quantum computers are ideal tools to perform unitary and even nonunitary propagation [59,60] with a possible polynomial scaling advantage over classical computers where the coupling potential term can be expressed as an entanglement of encoded qubits.[39] Third, the transformation from adiabatic wavefunctions to diabatic wavefunctions is unitary and can be easily implemented as parametric gates during state preparation on quantum computers.The geometric phase, [42,43,61] a global phase factor dressing the wavefunctions near CIs, can also be encoded with simple rotations in the Pauli basis, which is a natural advantage of quantum computers.
In this paper we evaluate the performance of quantum computers in describing CIs.Some key issues associated with CIs, such as seam curvature, optimization and geometric phase are discussed.We implement the electronic structure simulation of H + 3 with and without noise using the excited-state contracted quantum eigensolver (CQE) proposed in Ref. [62] in which the wavefunctions are generated by the ACSE ansatz.The theory and methodology of CI including an overview of CQE are presented in section II.Results and outlooks are further discussed in section III and IV, respectively.

II. THEORY A. Diabatic Hamiltonian matrix and CIs
The electronic Schrödinger equation can be written as in which H d is the quasi-diabatic Hamiltonian matrix.E J and d J are the adiabatic energies and wavefunctions of the Jth state respectively, corresponding to nuclear coordinate R. We consider a two-state example in which H d is a two by two matrix throughout this work, where most conclusions can be readily extended to additional electronic states.To obtain degenerate eigenvalues, we require where H IJ are matrix elements of H d .When wavefunctions are non-Hermitian as generated from nonunitary exponential ansatz (for example, standard coupled cluster ansatz), the above equation is the only constraint to form CIs. However, this creates a nonphysical (N −1) artifact that accompanies complex eigenvalues in the vicinity of true CIs, where N is the molecular degree of freedom.[20][21][22][23] The true CI is a submanifold of the potential energy surfaces where the following two constraint equations, one for diagonal and one for off-diagonal term, are simultaneously satisfied in the diabatic representation.
It is then easily recognized that the dimension of CIs is (N − 2).While the diagonal condition is easy to constrain, the off-diagonal condition is subtle because it is not directly available in electronic structure programs.For some molecules, as we will show in this paper, it is possible to find two states with high symmetry such that the couplings between them are strictly zero by symmetry, a situation known as symmetry-required CIs.The symmetry, however, does not serve as a necessary condition for the existence of CIs as some CIs occur in the more general category of "accidental" CIs.[16] B. Geometric phase effect Most modern quantum chemistry programs assume the Born-Oppenheimer approximation and thus, produce electronic structure data in the adiabatic representation.The adiabatic representation, however, fails to describe nonadiabatic dynamics because state couplings are always zero.The diabatic representation, on the contrary, produces well-described state couplings and smooth diabatic wavefunctions.There has been significant research effort directed towards determining the transformation from the adiabatic representation to the diabatic representation.One reason for the abundance of such diabatization techniques is that a strict diabatic representation does not exist for polyatomic molecules; hence, we refer to their states as "quasi-diabatic."[24] For two-state diabatization, the transformation is written as in which U is a unitary matrix.Some literature expresses the unitary as a rotation matrix parameterized by angle θ, We remind the reader that the expression might naturally lead to the assumption that θ is a continuous function of R, but this is not necessarily true in the presence of CIs due to the geometric phase effect.The geometric phase effect requires that wavefunctions that are transported around a path enclosing a CI acquire an additional phase factor.[61,63,64] A geometry-dependent and statedependent factor e iA K (R) (K = I, J) must be included in the adiabatic wavefunction.The natural advantage of using qubits to represent this two-state diabatization is that they are both isomorphic to the SU(2) group.Indeed on quantum computers, the phase factor can be implemented as a simple rotation gate parametrized by A K (R).One of the authors has shown in previous work [65] that A J (R) can be evaluated from the integral below, where the integrand is the well-known derivative coupling vector.[16] As this paper focuses on the topography of CIs, namely a more "static" description, a detailed analysis of the nonadiabatic quantum dynamics, the geometric phase factor and its implementation on quantum platforms is reserved for future work.

C. Variance-based contracted quantum eigensolver
The electronic structure calculation in this work is performed with a variance-based contracted quantum eigensolver (CQE) that has been proposed in a previous paper.[62] The algorithm is briefly reviewed here.
The variance (denoted as Var) of the system is defined as: We minimize the variance with respect to the parametric two-body anti-hermitian operator Fm , where the wave-function at the m th iteration is given by the unitary ansatz as where in which â † p and âp are the creation and annihilation operators, respectively.The key equation guiding the optimization is derived by taking the gradient of the variance with respect to Fm : in which Γpq st = â † p â † q ât âs and the equation defines the elements of the 2-RDM.Through a selfconsistent update of energy and 2 F m , we can converge the variance to a minima which corresponds to an excited or ground state.More details including an ancillary-assisted measurement of the variance has been reported in previous work.[62] We provide additional comments regarding why variance-based CQE is suitable to describe the CIs.The convergence depends on the choice of the initial guess, which can be generated from single Slater determinant or a linear combination of them.It will converge to the nearest minimum of the variance without knowledge of the lower states.Here by nearest, we mean the most similar in configuration composition.This state-specific feature can be beneficial in studying the CIs.It allows us to tackle a specific state during the slow variation of molecular geometry without concern that the adiabatic states will cross.Note this also coincides with the idea of configurational-uniformity-based diabatization as first proposed by Nakamura and Truhlar.[26]

III. RESULTS
We demonstrate the approach to computing the CI with the molecule H + 3 .The relative positions of the three co-planar hydrogen atoms are described in polar coordinates as (R, 0), (R, π) and (ρ, θ) where R ≥ 0, ρ ≥ 0, 0 ≤ θ < 2π, allowing us to represent the molecular geometry by the set of coordinates (R, ρ, θ).Calculations are performed with the IBMQ statevector simulator and FakeLagosV2 backend.The quantum simulation result is benchmarked with full configuration interaction calculations.All computations are performed in the minimal Slater-type orbital (STO-3G) basis set.Here and below we denote full configuration interaction as FCI to distin- guish it from the abbreviation for the conical intersection (CI).

A. Electronic structure of H + 3
The H + 3 molecule exhibits arguably the simplest CI.Nonetheless, despite its simplicity, the molecule is an important species in astrochemistry, providing a useful benchmark for the study of CIs.[66][67][68][69][70][71][72][73] We compute the first three states of H + 3 with S z = 0.A compact mapping is used to reduce the number of required qubits to three for the first and second excited states (denoted as E 1 and E 2 ) of H + 3 .The mapping is described here.We denote the configuration state function as |ij⟩, (1 ≤ i, j ≤ 3) with the S z = +1/2 electron occupying the i th molecular orbital and the S z = −1/2 electron occupying the j th molecular orbital.The dimension of FCI matrix is 9.A further reduction is performed by eliminating |11⟩ by observing that it has almost no coupling to the E 1 and E 2 states.Although |11⟩ can couple to other higher states and in principle affect the diagonalization result, the truncation has negligible effect on the energy of E 1 and E 2 (< 10 −8 hartree), resulting in a total qubit number of log 2 8 = 3.
We analyze the electronic structure property of H + 3 using the highly-symmetric D 3h point group.The first and second excited states correspond to the two components of an E ′ irreducible representation and thus form the symmetry-required CIs.We plot the potential energy curve in Fig. 1 obtained from the statevector simulator and a fake-backend simulator.A zoomed region of the degeneracy is given in the figure as well.The two excited states always overlap in the FCI scheme, which is consistent with our electronic structure knowledge of the system.On a noiseless statevector simulator we achieve an energy accuracy of 10 −6 hartree, where the only error comes from the trotterization, proving the exactness of the ACSE ansatz.After introducing device noise, the error for each individual state is around 12 mhartree.It is worth noting that since the error is quite uniform for both states, the error of their energy gap is significantly smaller, which is quite promising for predicting the energetic degeneracy.
We next plot the dissociation curve at a lower symmetry, namely C 2v in Fig. 2. The discontinuity of the E 2 curve is due to crossings with intruder states.We are particularly interested in the CIs between E 1 and E 2 , where the two states coincide at a D 3h geometry.Despite the relatively large error of individual states, the prediction of the location of the CIs is surprisingly accurate (<0.01 bohr).As mentioned before, if the error induced by noise is nearly uniform for both states and for all geometries, then the effect of noise is only to shift both potential energy surfaces by a similar amount, which should not significantly affect the topography of the CIs.To verify this, in Fig. 3 we plot the coupled potential energy surfaces as a function of the coordinates of the third hydrogen atom.It can be seen that the topography of the CIs is well reproduced.The expected cusps induced by random noise are barely discernible due to the uniformity of noise.We note, however, that although the potential energy surfaces and relative energy gap are well reproduced, the absolute error still remains challenging on noise intermediate-scale quantum (NISQ) devices and further error mitigation techniques are needed.

B. Locating the seam of CI
The search of the minimum energy CI (MEX) is done by minimizing the constrained Lagrangian where C k are certain geometry constraint equations and ∆E IJ = 0 constrains the CIs in the adiabatic representation.In previous work, we showed that the gradient of Lagrangian corresponding to the geometry parameter set R can be obtained with classical calculations.[74] Here we use a hybrid method, where single-point calculations are performed on quantum computers and the gradient is obtained numerically and classically by varying geometries.
The dimension of CI for the triatomic molecule is only one.By varying one molecular coordinate and fixing the rest, we should obtain a one-dimensional curve that corresponds to the seam of the CI.For the special case between the E 1 and E 2 of H + 3 , we know that such curve is unique and corresponds to the D 3h geometries in Fig. 1.We report the optimization results by setting R to 2.0 bohr and optimizing over the position of the third hydrogen (r, θ).To keep things simple, we used a gradient-like Newton-Raphson method with fixed step size of 0.1 where the Hessian of the Lagrangian is approximated by the identity matrix.A typical optimization process is shown in Table I.The performance is quite robust despite the simple setup, suggesting the gradient from quantum devices is resilient enough for geometry optimization purposes.The energy difference decreases during the optimization, except in the first iteration.This exception can occur because the Lagrangian includes contributions beyond the energy difference.An important observation is that, the noise on NISQ simulators introduces a small oscillation around our targeted D 3h CI.The errors in the bond distance and bond angle, however, are quite small, which helps to demonstrate the accuracy of our description of the CIs.

IV. CONCLUSIONS AND OUTLOOK
Current state-of-art nonadiabatic quantum dynamics are limited to small molecules due to their exponential scaling with respect to the active vibrational modes.Quantum computers may potentially provide a solution.In the future quantum devices with hundreds of qubits may be able to perform nonadiabatic quantum dynamics simulations that are either too expensive or intractable on high-performance classical computers.This paper provides a foundation for simulating the CIs on quan-tum computers, paving the way for advancements in quantum-based simulations for nonadiabatic molecular systems.Using a variance-based CQE, the energies of intersecting potential energy surfaces of molecular H + 3 are accurately computed.The study achieves a correct representation of CI topography through a Hermitian wavefunction generated by the ACSE ansatz.Future work includes realizing diabatization and quantum dynamics of complex nonadiabatic molecular system on quantum devices.

FIG. 1 .
FIG. 1. Potential energy curve calculated by FCI as well as noiseless and noise simulators.Molecule is treated in the D 3h symmetry, where the polar coordinates of the three hydrogen atoms are (R, 0),(R, π) and ( √ 3R, π/2).Note that E1 and E2 are degenerate in FCI result due to symmetry.

FIG. 3 .
FIG. 3. Intersecting adiabatic potential energy surfaces calculated from variance CQE.The grid size is 20×12.Additional points are placed in the vicinity of the CI.For better illustration, we use Cartesian coordinate with the coordinates of the three atoms being (-1,0,0), (1,0,0) and (x,y,0).

TABLE I .
Energy gap during a geometry optimization with respect to θ and ρ starting from random guesses for these variables.The remaining parameter R is fixed at 2.0 bohr.