Geometric Rotation of the Nuclear Gradient at a Conical Intersection: Extension to Complex Rotation of Diabatic States

: Nonadiabatic dynamics in the vicinity of conical intersections is of essential importance in photochemistry. It is well known that if the branching space is represented in polar coordinates, then for a geometry represented by angle θ , the corresponding adiabatic states are obtained from the diabatic states with the mixing angle θ /2. In an equivalent way, one can study the relation between the real rotation of diabatic states and the resulting nuclear gradient. In this work, we extend the concept to allow a complex rotation of diabatic states to form a nonstationary superposition of electronic states. Our main result is that this leads to an elliptical transformation of the e ﬀ ective potential energy surfaces; i.e., the magnitude of the initial nuclear gradient changes as well as its direction. We fully explore gradient changes that result from varying both θ and ϕ (the complex rotation angle) as a way of electronically controlling nuclear motion, through Ehrenfest dynamics simulations for benzene cation.


INTRODUCTION
Conical intersections are points of degeneracy between two (or more) electronic states which exist in molecules by nature 1−4 or can be induced by light. 5 They are important in photochemistry since they serve as funnels allowing radiationless electronic transitions. 6−8 The control of reactivity in the vicinity of conical intersections would allow one to control the possible products of photochemical reactions. de Vivie-Riedle et al. proposed to "steer" chemical reactions near a conical intersection by using an external electric field to control the amplitudes and relative phases of the populated electronic states. 9−11 We have recently shown how a superposition of electronic states influences the initial nuclear dynamics next to a conical intersection in toluene cation. 12 Here we generalize this idea, by systematically exploring the effect of the relative phase of the superposition on the initial nuclear motion, through changing the gradient of an effective potential energy surface.
The control of nuclear dynamics near conical intersections requires a clear understanding of the forces experienced by the nuclei due to the electrons. The electronic degeneracy is lifted along two nuclear coordinates: the gradient difference vector ⎯ → X 1 and the derivative coupling vector ⎯→ X 2 . These span a twodimensional subspace of the nuclear coordinates called the branching space, the remaining dimensions forming the intersection space. In the branching space, the adiabatic potential energy surfaces form a double cone with a vertex at the origin.
It is well-known that near a conical intersection the adiabatic states at a geometric polar angle of θ are constructed from the real rotation of two diabatic basis states with an angle θ/2. 13−16 It is equivalent to state that a real rotation of two diabatic basis states with an angle θ/2 results in a geometric rotation of the nuclear gradients by an angle θ. In this work, it is more convenient to speak in the latter terms. But one has to keep in mind that we are not dealing with the gradient of an adiabatic state energydiscontinuous at the exact point of degeneracy but rather the gradient of the energy of a superposition of two diabatic basis states which defines another diabatic state energy, the gradient of which is continuous everywhere.
In this work we extend this concept to complex rotation of two diabatic basis states defining a relative phase ϕ as well as the relative weight determined by θ/2. We will present an analytical expression for the resulting nuclear gradient, with a graphical interpretation using the Bloch sphere. The fact that the nuclear gradient depends on the composition of the superposition could be used for coherent electronic control as θ and ϕ can both be changed experimentally in principle. We will demonstrate this idea with nonadiabatic dynamics simulations in benzene radical cationa typical example of a Jahn−Teller system. 17−20 Quantum dynamical simulations 21,22 have been performed on benzene cation using the MCTDH package: 23,24 using the first five electronic states of the benzene radical cation, Koppel and co-workers were able to show that the deexcitation of these states to the electronic ground state occurs on a time scale of 100 fs. 23 Here, we are interested in short time scale nuclear dynamics, with the two lowest energy electronic states populated initially. We will use the Ehrenfest method, 25 able to describe the coupled electron−nuclear dynamics of small organic molecules as previously exhibited. 12,26−28 In the case of an untilted conical intersection (similar to Jahn−Teller systems), the relation ⟨ψ I |∇ ⃗ Ĥe|ψ I ⟩ = −⟨ψ II |∇ ⃗ Ĥe|ψ II ⟩ is valid to first order. Using it, we come to the expression

Journal of Chemical Theory and Computation
Article Note that the vectors ⎯ → X 1 and ⎯→ X 2 are determined by the diabatic states |ψ I ⟩ and |ψ II ⟩ and therefore are the same for all superpositions |χ⟩.
It is useful here to introduce polar coordinates (r,α) for the nuclear gradient instead of the Cartesian coordinates (x 1 ,x 2 ) = (cos(θ), sin(θ) cos(ϕ)) along the ⎯ → X 1 and ⎯→ X 2 axis: r is the magnitude from the origin, and α is the counterclockwise angle from the ⎯ → X 1 axis. The magnitude of the nuclear gradient of the energy of |χ⟩ in the branching space is expressed by and its direction by α θ ϕ θ = atan2(sin( ) cos( ), cos( )) (8) where the function atan2 is the arctangent function with two arguments.
When ϕ = 0°, we get back the results of a real rotation which were first discussed by Longuet-Higgins: 13−16 eq 6 becomes ∇ ⃗ ⟨χ|Ĥe |χ⟩ = cos(θ) ⎯ → X 1 + sin(θ) ⎯→ X 2 ; i.e., the rotation of diabatic states with an angle θ/2 leads to a geometric rotation of the nuclear gradient by an angle α = θ without affecting the magnitude of the nuclear gradient r = 1. Here, eqs 6−8 give the more general results of a complex rotation of diabatic states: the relative phase ϕ affects the component of the gradient along ⎯→ X 2 , but not the component of the gradient in the direction of ⎯ → X 1 . For a given relative phase ϕ, the magnitude of the gradient is an ellipse (eq 7) and the geometric rotation α now depends on both θ and ϕ in a less trivial way (eq 8).
One can use the Bloch sphere to represent graphically the nuclear gradient of the energy of a superposition of two basis states. Indeed, the projection of the Bloch vector representing the superposition |χ⟩ onto the real plane of the sphere has Cartesian coordinates (x = cos(θ), y = sin(θ) cos(ϕ)); these match the Cartesian coordinates of the nuclear gradient in the branching space (eq 6). One can therefore interpret the projected vector (blue arrow in Figure 2) as the corresponding nuclear gradient in the branching space ( ⎯ → X 1 , ⎯→ X 2 ). The nuclear gradient of the energy of the orthogonal superposition |χ′⟩ has the same magnitude but opposite direction: In polar coordinates, we have r′ = r and α′ = α + 180°. Note that an electronic wave function |χ′⟩ = cos(θ/2)|ψ I ⟩ + e iϕ sin(θ/ 2)|ψ II ⟩ and its complex conjugate wave function |χ*⟩ = cos(θ/ 2)|ψ I ⟩ + e −iϕ sin(θ/2)|ψ II ⟩ lead to Bloch vectors that are mirror images with respect to the real plane: they give the same projection onto the real plane and thus the same nuclear gradient in the branching space. The reason is that the sign of the relative phase ϕ is not relevant and only its absolute valuethe absolute dif ference of phasesdetermines the nuclear gradient.
2.3. Effective Potential Energy Surfaces. One can construct a representation of a potential energy surface to first order using the nuclear gradient. Here, we deal with superpositions of diabatic states. Although, at the exact point of degeneracy, a superposition of the two real diabatic basis states is also an eigenstate of the electronic Hamiltonian, a complex superposition of the diabatic states is in general not an electronic eigenstate outside the intersection seam. The nuclear gradients of the energies of the two orthogonal superpositions |χ⟩ and |χ′⟩ lead thus to ef fective potential energy surfaces. These represent the averaged forces felt by the nuclei due to the electrons.
We will use for reference the adiabatic potential energy surfaces obtained with ϕ = 0°(and 0°≤ θ ≤ 360°); they are schematically represented by the isotropic inner cones in Figure  3 and Figure 4. The well-known dependence of the nuclear

Journal of Chemical Theory and Computation
Article gradient with respect to the mixing angle θ/2 is indicated by the coloring of the potential energy surfaces: yellow for high population in |ψ I ⟩ and therefore a gradient close to ⟨ψ I |∇ ⃗ Ĥe|ψ I ⟩, red for high population in |ψ II ⟩ and therefore a gradient close to ⟨ψ II |∇ ⃗ Ĥe|ψ II ⟩ and orange for intermediate cases.
Increasing the absolute value of ϕ decreases the magnitude of the gradient along ⎯→ X 2 (eq 6). This leads to an elliptical distortion of the effective conical intersection (eq 7) as shown in Figure 3. In the case of ϕ = 90°, where there is no contribution of the nuclear gradient along ⎯→ X 2 , the effective potential hypersurfaces are two intersecting flat sheets (not shown).
The reference case of ϕ = 0°also corresponds to r = 1 (eq 7). One can modify both θ and ϕ to decrease the value of r, leading to a symmetric flattening of the effective potential energy surfaces, see Figure 4. In the limiting case of r = 0 (obtained with θ = 90°and ϕ = 90°or 270°), the two surfaces become completely flat (not shown) and the magnitude of the gradient in the branching space becomes zero: there would not be any nuclear motion in the branching space.
We have shown how the nuclear gradient depends on the composition of the superposition of electronic diabatic states. This idea can be used to control the nuclear motion via a coherent superposition of electronic states.

ILLUSTRATIVE EXAMPLE: BENZENE CATION
We now illustrate the effective potential that the nuclei feel due to the electronic superposition in benzene cation, a typical Jahn−Teller system. For this, we simulated the first 5 fs following the instantaneous valence ionization. If dynamics occurs on a single electronic eigenstate, one can use the Born− Oppenheimer molecular dynamics method to simulate the adiabatic dynamics. However, in general, the complex electronic wave functions are nonstationary electronic states: here, we used the Ehrenfest method to describe the nonadiabatic dynamics that occurs.
3.1. Computational Details. The equilibrium geometry of the neutral benzene has D 6h symmetry. It lies within the seam of a conical intersection between the two lowest electronic eigenstates of the cationic species, corresponding to ionization from the degenerate occupied π orbitals (HOMO and HOMO-1). 34 The geometry of the conical intersection corresponds to (x 1 ,x 2 ) = (0,0) in the branching space. The surrounding hypersurface of the adiabatic ground state of the cation is of C 3 symmetry: 25,27,35 there are three equivalent minima and three equivalent transition structures forming a moat around the conical intersection. At each of these minima, the adiabatic ground and first excited states have respectively a quinoid and antiquinoid character. 25,27,28 We have chosen the adiabatic ground and first excited states at one of these equivalent minima to define the diabatic basis states |ψ I ⟩ (quinoid) and |ψ II ⟩ (antiquinoid). These two diabatic basis states determine the directions of the diabatic gradient difference ⎯ → X 1 and diabatic derivative coupling ⎯→ X 2 ; see Figure 5. The selected minimum structure of the cationic ground state therefore lies in the direction of ⎯ → X 1 and corresponds to (x 1 ,x 2 ) = (1,0). The electronic structure was treated with the CASSCF method. To include the static correlation and part of the dynamic correlation, an active space consisting of all six π orbitals and five π electrons was chosen. The electronic structure was calculated with the 6-31G* basis set. The quality of the CASSCF[5,6]/6-31G* potential hypersurface around the conical intersection was validated against explicitly correlated multireference configuration interaction (MRCI-F12); see the Supporting Information. For the dynamics, we used the secondorder Ehrenfest method implemented in a development version of Gaussian. 36 The quantum electronic wave function is propagated as a solution of the electronic time-dependent Schrodinger equation. The nuclei are treated as classical particles by integration of Newton's equations of motion. Our implementation of the nuclear gradient includes the non-Hellmann−Feynman terms approximately. The term due to the molecular orbital derivative is neglected because the orbitals are assumed not to change much over a time step. The term due to the configuration interaction derivative is obtained solving the coupled perturbed equation but neglecting the phase of the time-dependent wave function. These approximations result in some small errors in the gradient, which are corrected by a fifthorder Hessian based predictor−corrector scheme. 37 Here, we used a step size of 0.006 amu 0.5 bohr, corresponding to a time step of approximately 0.1 fs. More details of our CASSCF implementation of the Ehrenfest method can be found in our previous work. 25 The simulations are started at the Franck−Condon point (a point on the seam of the conical intersection) with no initial kinetic energy. As our previous work on toluene shows, sampling of the phase space does not change the qualitative nature of the trajectories in the Ehrenfest approximation. 12 Here, our focus lies in the initial nuclear motion in the branching space, although motion in the intersection space also happens in the simulations. Note that the interaction between the photoelectron and the resulting cation is neglected.
3.2. Adiabatic Dynamics. Adiabatic states are obtained with ϕ = 0°, which also leads to r = 1 in eq 7. To show the effect of the relative weight θ, we now present the results of initial nuclear dynamics upon ionization for a reference set of initial electronic wave functions: we varied the mixing angle from θ = 0°to 350°in steps of 10°while fixing the relative phase to ϕ = 0°.
In Figure 6, bottom right, the initial electronic wave functions are represented by vectors on the Bloch sphere.

Article
Since ϕ = 0°, they all lie in the real plane. Above that, the projections of the Bloch vectors onto the real plane (here trivial) are shown. This gives the initial nuclear gradients in the branching space, all lying on the unit circle. In Figure 6, left, the trajectories up to 5 fs are shown. The projection onto the bottom plane shows the motion in the branching space. The projections onto the side walls show the time evolution of the position in the direction of gradient difference and derivative coupling.
All trajectories are straight lines in the branching space, indicating nuclear gradients with time-independent directions α. This is consistent with the fact that the electronic wave functions are adiabatic states, i.e., stationary states. All trajectories are approximately of the same length; thus, nuclear motion proceeds with the same velocity independent of the mixing angle θ/2. This makes sense because the magnitude of the initial nuclear gradients is constant r = 1 for all θ. Here we can clearly see, as expected, how the mixing angle θ/2 determines the direction of initial nuclear motion in the branching space by a geometric rotation of α = θ. Remember that ⎯→ X 2 is the diabatic derivative coupling and not the adiabatic derivative coupling. Therefore, adiabatic dynamics along the diabatic ⎯→ X 2 should not be surprising. To summarize, in adiabatic dynamics, by choosing the mixing angle θ/2 we can predict the initial direction of the nuclear motion from the conical intersection while keeping the velocity constant.
3.3. Nonadiabatic Dynamics and the Relative Phase ϕ. To study the effect of the relative phase ϕ, we present the results of initial nuclear dynamics upon ionization for two new sets of initial electronic wave functions and compare them to the reference set of the previous subsection. We again varied the mixing angle from θ = 0°to 350°in steps of 10°while this

Journal of Chemical Theory and Computation
Article time fixing the relative phase to ϕ = 45°and 90°, in the first and second sets, respectively.
The first set of initial electronic wave functions are represented by Bloch vectors in Figure 7, bottom right. The projected points onto the real plane all lie on an ellipse as predicted by eq 7: the magnitude of the initial nuclear gradients is preserved along the gradient difference direction but diminished along the derivative coupling, compared to the reference set with ϕ = 0°( Figure 6, top right). The reason for this is in eq 6: the second term becomes smaller as cos(45°) = 1/√2. The trajectories show an initial motion in the direction of the derivative coupling that is reduced while the initial motion in the direction of the gradient difference is similar to the reference set ( Figure 6, left). This corresponds to motion on an elliptical cone (Figure 3). The mixing angle θ/2 still determines the direction of the nuclear motion α. However, by fixing ϕ = 45°, we scale down the nuclear motion along ⎯→ X 2 .
Note that here the trajectories in the branching space are not straight lines but they are slightly curved. Why? The nuclear gradients break D 6h symmetry, and the trajectories move in the branching space away from the conical intersection. In the reference set, the real superpositions become the pure adiabatic ground state. Here, the complex electronic wave functions are still a superposition of the two lowest adiabatic states (not degenerate anymore); they are not stationary, i.e., the angles θ and ϕ in eq 3 evolve with time, leading to time-dependent directions of nuclear gradients (eq 8).
The second set of initial electronic wave functions are represented by Bloch vectors in Figure 8, bottom right. The projected points onto the real plane lie all on the ⎯→ X 2 axis: the magnitude of the initial nuclear gradients is again preserved along the gradient difference direction compared to the reference set with ϕ = 0°( Figure 6), but is now null along the derivative coupling. The reason for that is in eq 6 again: the

Journal of Chemical Theory and Computation
Article second term becomes zero as cos(90°) = 0. The trajectories show approximately no initial motion in the direction of the derivative coupling while the initial motion in the direction of the gradient difference is similar to the reference and last sets. This corresponds to a nuclear motion on a flat sheet, an infinitely elongated ellipse. The mixing angle θ/2 still determines the direction of the nuclear motion α. However, by fixing ϕ = 90°, we suppress the initial nuclear motion along ⎯→ X 2 . As a consequence, changing the mixing angle θ/2 allows one to choose the magnitude of the nuclear motion along ⎯ → X 1 and its sign. The slight deviations from x 2 = 0 are due to the time evolution of the electronic wave functions that are nonstationary.
3.4. Nonadiabatic Dynamics and the Gradient Magnitude r. In the reference set, the mixing angle θ/2 allows one to choose the initial nuclear motion but all trajectories show the same length after 5 fs; i.e., the nuclei have moved at the same velocity: all initial gradient magnitudes are r = 1. We now investigate symmetric distortion of the effective potential energy surfaces with two other sets of initial electronic wave functions. We varied the initial gradient direction from α = 0°to α = 350°in steps of 10°while this time fixing the gradient radius to r = 0.5 and 0, in the third and fourth sets, respectively. The angles θ and ϕ are chosen corresponding to eqs 7 and 8.
In Figure 9 we can see the third set of initial electronic wave functions represented by Bloch vectors. In this set, the projected points all lie on a circle with a radius of r = 0.5. In Figure 9, left, the trajectories of nuclear motion of this set are shown. They are very similar to those of the reference set but smaller in magnitude, indicating slower nuclear motion on a flatter cone. Also, they are curved lines because here, the electronic wave functions are nonstationary and the nuclear gradient directions evolve with time.
The fourth set is actually only made of two superpositions {θ = 90°, ϕ = 90°} and {θ = 90°, ϕ = 270°}; see Figure 10. For both, the diabatic states are equally weighted. This leads to initial nuclear gradients vanishing in the branching space (blue point at (x 1 ,x 2 ) = (0,0) in Figure 10, top right): the initial nuclear gradients have D 6h symmetry. The corresponding nuclear trajectories stay in the intersection space (x 1 ,x 2 ) = (0,0); the electronic diabatic states stay degenerate, and therefore the superposition is stationary. Note that even if there is no nuclear motion in the branching space, motion still occurs in the intersection space where the degeneracy is not lifted.

CONCLUSION
We have extended the idea of geometric rotation due to the real rotation of diabatic states and presented an analytical expression for the nuclear gradient of the energy of a complex superposition of two diabatic electronic states in the vicinity of a conical intersection. For that we have introduced the Bloch sphere to represent the superposition as a vector in this context. The two angles determining the position of the vector are the mixing angle θ/2 and the relative phase ϕ. We have shown that the projection of this vector onto the real plane can be interpreted as the nuclear gradient in the branching space. Different superpositions of electronic states lead to different effective potential energy surfaces.
We have illustrated the theory with nonadiabatic Ehrenfest dynamics in benzene cation. For this, trajectories with different initial electronic wave functions generated by systematically changing the two angles θ and ϕ were run. The results show that real-valued superpositions induce nuclear motion of the same magnitude with direction in the branching space determined by the mixing angle θ/2. Considering complexvalued superpositions allows one to also control the velocity as well as the direction of nuclear dynamics in the branching space. Although we have not considered a specific experiment, θ and ϕ can both be varied experimentally, providing a route toward coherent electronic control of initial nuclear motion following ionization.