Electronic Control of Initial Nuclear Dynamics Adjacent to a Conical

Article pubs.acs.org/JPCA

Electronic Control of Initial Nuclear Dynamics Adjacent to a Conical Intersection Morgane Vacher, Jan Meisner,† David Mendive-Tapia,‡ Michael J. Bearpark, and Michael A. Robb* Department of Chemistry, Imperial College London, London SW7 2AZ, United Kingdom S Supporting Information *

ABSTRACT: Photoionization can create a nonstationary electronic state and therefore initiates coupled electron−nuclear dynamics in molecules. Using a CASSCF implementation of the Ehrenfest method, we study the nuclear dynamics following vertical ionization of toluene, starting close to the conical intersection between ground and first excited states of its cation. The results show how the initial nuclear dynamics is controlled by the nonstationary electronic state character. In particular, ionization of this system leading to an equal superposition of the two lowest energy states can initiate nuclear dynamics in an orthogonal direction in the branching space to dynamics on the ground or first excited state potential energy surfaces alone.

■

INTRODUCTION Photoionization can create a coherent superposition of electronic states (i.e., a nonstationary state) and therefore initiates coupled electron−nuclear dynamics in molecules. Experiments suggest that these may play a chemical role: for example, Weinkauf et al. suggest that fragmentation sites in small polypeptide cations are controlled by the dynamics of a superposition of quasi-degenerate electronic states.1,2 Geometries where two electronic states (at least) become degenerate are called conical intersections (CIs). They are important in photochemistry as they allow efficient nonradiative electronic transitions3−5 and may lead to several products. In the vicinity of a CI, the Born−Oppenheimer approximation breaks down because of a potentially strong interaction between electronic and nuclear degrees of freedom. The control of electron dynamics by an external electric field during molecular dissociation has been demonstrated experimentally in diatomics.6 De Vivie-Riedle extended the concept theoretically to larger molecular systems and suggested to control the relative phases and amplitudes of a superposition of electronic states near a CI to “steer” chemical reactions.7,8 Studying how the nuclear motion is affected by a nonstationary electronic state is challenging. We have used our CASSCF implementation of the Ehrenfest method9−11 as a general method to simulate coupled electron−nuclear dynamics on-the-fly for any system. Here, the classical nuclear motion is determined by the gradient of the superposition of electronic states (i.e., the relative phases and amplitudes) rather than the gradient of a single adiabatic electronic state (eigenstate of the Born−Oppenheimer Hamiltonian). Differences in the gradients will lead to different motions of the nuclei in each case.12 In this article, we demonstrate explicitly how the initial coherent population of two electronic adiabatic states affects the subsequent nuclear motion in a real molecular system. We also study the effect of sampling on nuclear dynamics using a Wigner distribution for the initial geometries and velocities. © XXXX American Chemical Society

Note that in this work we focus on the nuclear dynamics following the creation of a superposition of electronic states, not on the detailed electron dynamics. Ehrenfest is an approximate method that could eventually lead to nonphysical asymptotic behavior in the case of bifurcating paths due to the single-configuration approach. In our simulations, we start with a coherent superposition of electronic states, and each nuclear trajectory follows the gradient of this superposition. In a full quantum treatment, each nuclear wave packet would evolve on its electronic state under coupling with the other electronic states. We expect the Ehrenfest method to be valid at short times before the nuclear wave packets move too far apart from each other. We choose the methyl-substituted benzene − toluene − as an example because vertical ionization takes place at a geometry slightly displaced from the CI between ground and first excited states of the cation. The resulting nonzero energy gap is needed to generate the nonstationary electronic states we study. Figure 1 represents the ionization process. Two nuclear coordinates lift the subsequent degeneracy: the gradient difference vector (GD) and the derivative coupling vector (DC). These are defined in the lower part of Figure 2: the GD corresponds to the shortening of the side carbon−carbon bonds, and the DC is a shearing motion where the top carbon atom of the ring moves to the left and the bottom carbon to the right. The two-dimensional space made of these two vectors is called the branching space. Figure 2 is a projection of the surrounding “moat” of the Jahn−Teller-like conical intersection13−15 in the branching space (also shown centrally in Figure 1). The red dot represents the equilibrium geometry of the neutral species, displaced from the crossing along GD. The Special Issue: Jacopo Tomasi Festschrift Received: September 27, 2014 Revised: November 26, 2014

A

dx.doi.org/10.1021/jp509774t | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 2. Diagram representing the structures in the moat around the conical intersection in toluene. The quinoid structure on the right is a minimum on the ground state of the cation, while the antiquinoid structure on the left is a transition state. The red circle represents the equilibrium geometry of the neutral species. The branching space consists of two vectors, the gradient difference (GD) and the derivative coupling (DC) (shown in the lower part of the figure): the GD is the shortening of the side carbon−carbon bonds, and the DC is a shearing motion where the top carbon atom of the ring moves toward the left and the bottom carbon toward the right.

Figure 1. Ionization process in toluene: the lower surface represents the ground state of the neutral species and the upper two represent the ground and first excited states of the cation in the branching space of the conical intersection (GD is the gradient difference vector and DC the derivative coupling vector). A projection of the cationic conical intersection moat expressed in terms of valence bond structures is shown centrally.

■

THEORETICAL METHODS Initial Conditions. Our aim is not to reproduce any particular photoionization experiment but rather to extract physical insights from some model calculations. For this reason, we compare the nuclear dynamics initiated with different electronic wave functions. A coherent superposition of the two lowest-energy eigenstates |ψ0⟩ and |ψ1⟩24 can be expressed using the Bloch representation:

moat contains two types of valence bond (VB) resonance structures: three quinoid and three antiquinoid. The quinoid structure along the GD is the minimum of the ground state of the cation with the side carbon−carbon bonds shortened to 1.37 Å; the antiquinoid structure along the GD is a transition state with the side carbon−carbon bonds lengthened to 1.44 Å. Since the branching space defines the conical intersection, we are interested in how motion in this space can be controlled by changing the initial superposition of states. Note that although the substitution of a hydrogen by a methyl group removes the orbital degeneracy,16,17 the toluene cation can be seen as a slightly perturbed benzene cation (a Jahn−Teller18,19 system). The latter has been studied using simple analytical models in two dimensions.20 More recently, quantum dynamics simulations in benzene cation were carried out21,22 using the MCTDH package,23 in which a reduced number of nuclear coordinates was selected and the potential energy surfaces in this subspace were computed and fitted. Our aim is rather to illustrate the principles of nonadianatic dynamics initiated with a superposition of two (or more in general) adiabatic states using a general method, a well-defined approximation, that could easily be applied to other systems in full dimensionality.

|Ψ⟩ = sin(θ )|ψ0⟩ + cos(θ ) eiϕ|ψ1⟩

(1)

Two angles are needed to define the wave function: θ determines the relative weight (0° ≤ θ ≤ 90°) and ϕ is the relative phase between the two states (0° ≤ ϕ ≤ 180°). The superposition is by definition normalized since the electronic eigenstates are orthonormal and (sin2(θ) + cos2(θ) = 1. Note that the absolute phase of the total wave function is meaningless. Also, if only one eigenstate is populated, then the relative phase ϕ becomes the absolute phase and is therefore meaningless. In our implementation, we use adiabatic states to define the electronic wave function.25 One could consider the two quasi-diabatic states involved in the CI: these are the “quinoid” and “antiquinoid” states. At the equilibrium geometry of the neutral species, the cationic adiabatic states coincide with the quasi-diabatic states; the ground adiabatic state has a quinoid character and the first excited adiabatic state has an antiquinoid character. The equilibrium geometry (minimum on the ground state) of the neutral species was optimized at the B3LYP/6-31G* level. Its position in the branching space is represented by a white circle in Figure 3: it is displaced from the cation CI along GD. The zero-point energy (ZPE) of all 30 degrees of freedom B



Article

Note that the (quantum) sampling carried out implies a harmonic approximation of the potential energy surface (PES) in rectilinear normal modes. It has recently been demonstrated in the case of a classical distribution that the use of curvilinear normal modes reduces the sampling error without going beyond the harmonic approximation.27 CASSCF Implementation of the Ehrenfest Method and Analysis. Our method and implementation are described in detail in our previous article.28 Here, we only mention the important points and provide technical details. In the Ehrenfest method, the nuclear dynamics is treated classically by solving Newton’s equation of motion. The Hessian-based predictor−corrector algorithm designed by Hase and Schlegel29 is used with a mass-weighted step size of 0.0025 amu1/2 bohr (corresponding to a time step of about 0.04 fs). The electron dynamics is treated quantum mechanically by solving the time-dependent Schrö dinger equation. The electronic structure is computed using the state-averaged CASSCF30 method. Using the standard 6-31G* basis set, we choose the six π orbitals as active. Nuclear motion is analyzed in the two dimensions of the branching space by projecting the geometry along the GD and DC vectors, taking the CI as the geometrical origin. Note that the lengths of the GD and DC vectors are scaled so that the minimum of the ground state of the cation is located at (1,0). Coupled Electron−Nuclear Dynamics. The initial electronic density and also the energy gradient will be different for the superposition of states compared to the single adiabatic states. Indeed, the electronic density of the superposition of two electronic states defined in eq 1 is

Figure 3. Initial distribution of the sampled nuclear geometries along the two directions of the branching space. The ZPE of the equilibrium geometry of the neutral species (represented with the white circle) is sampled by 500 points as a Wigner distribution. The white square is the average position of the ensemble of points. The lengths of the GD and DC vectors are scaled so that the minimum of the ground state of the cation is located at (1,0).

was sampled (500 points of the phase space) with a Wigner distribution using the Newton-X package.26 The resulting distribution along the two vectors of the branching space is shown in Figure 3 as a heat map. It is centered on the equilibrium geometry of the neutral species (white circle). The average position of the ensemble of points is represented by the white square: it is very close to the equilibrium geometry as expected. The slight displacement is due to the finite sampling. To study the effect of sampling the initial geometries and velocities, an “unsampled” trajectory is simulated starting at the equilibrium geometry of the neutral species with no initial kinetic energy.

|Ψ|2 = sin 2(θ )|ψ0|2 + cos2(θ )|ψ1|2 + sin(2θ ) cos(ϕ) Re(ψ0*ψ1)

(2)

The first two components are the average electronic density where each adiabatic electronic density is weighted by the occupation of the adiabatic state. The third term corresponds to the interference between the two electronic states.31

Figure 4. Dependence of the electronic density (eq 2) and nuclear gradient (eq 3) on the angles θ and ϕ. (a) Coefficients determining the amplitude of the first two terms in each equation. (b) Two coefficients determining the amplitude of the third term. Gray arrows highlight the values corresponding to some specific angles. C



Article

Figure 5. Time evolution of the ensemble of the nuclear trajectories in the branching space: along the GD on the left (a,c,e) and along the DC on the right (b,d,f). The initial electronic wave function is |ψ0⟩ for the upper panels (a,b), 1/√2(|ψ0⟩ + |ψ1⟩) for the middle panels (c,d), and |ψ1⟩ for the lower panels (e,f). The initial geometries and velocities are sampled from a Wigner distribution of the equilibrium geometry of the neutral species (see Figure 3). The dashed white line shows the average evolution for the ensemble of trajectories. The solid white line shows the evolution of the unsampled trajectory starting at the equilibrium geometry of the neutral species with no initial kinetic energy.

considered here to simplify the expression, the coupledperturbed equations are in practice solved approximately.28 The electronic density and nuclear gradient in the branching space are determined by the composition of the electronic wave function through the angles θ and ϕ: - the amplitudes of the first two terms are controlled by the populations of the adiabatic states sin2(θ) and cos2(θ), respectively. The evolution of these quantities is shown in Figure 4a. - the amplitude of the third term is given by the product of sin(2θ) and cos(ϕ). The evolution of those quantities is shown in Figure 4b.

The associated energy gradient with respect to the position of nucleus I is32 ∇⟨Ψ| /e|Ψ⟩ = sin 2(θ )⟨ψ0|∇I /e|ψ0⟩ + cos2(θ )⟨ψ1|∇I /e|ψ1⟩ I + sin(2θ ) cos(ϕ)⟨ψ0|∇I /e|ψ1⟩

(3)

where /e is the electronic Hamiltonian. The first two components are the average gradient weighted by the occupation of the adiabatic states: in the case of toluene, they are along the GD in opposite directions. The third component is due to the change in occupation of the adiabatic states because of nonadiabatic transitions: here, it is along the DC. Note that although only the Hellmann−Feynman terms are D



Article

character and all but the first of the gradient terms in eq 3 vanish. As a result, the ensemble evolves toward positive GD (Figure 5a) and barely moves along the DC (Figure 5b): it follows the ground state gradient ⟨ψ0|∇I /e|ψ0⟩ and evolves toward the quinoid-like structure, which is the minimum on the ground state surface. In other words, the ensemble of trajectories moves to the right of the diagram in Figure 2: the side carbon−carbon bonds are shortened by approximately 0.01 Å. The ensemble of trajectories spreads by ∼25−30% (see Table 1): the shape of the ground state PES does not drive the trajectories toward a specific structure but allows them to spread out around the minimum (structure Min in Figure 2). If only the first excited state of the cation is populated (θ = 0°), the amplitudes of the three terms in eqs 2 and 3 are sin2(θ) = 0, cos2(θ) = 1, and sin(2θ) cos(ϕ) = 0, respectively (Figure 4). The initial electronic density |ψ1|2 has an antiquinoid character. It again barely moves along the DC (Figure 5f), but this time evolves toward negative GD (Figure 5e): it follows the gradient of the first excited state ⟨ψ1|∇I /e|ψ1⟩ and goes toward the CI, in the opposite direction compared to the dynamics on the ground state. In other words, the ensemble of trajectories moves to the left of the diagram in Figure 2: the side carbon− carbon bonds are lengthened. The ensemble of trajectories only spreads by ∼5−10% (see Table 1) because the conical shape of the excited state PES constrains the trajectories toward the CI. (In the Supporting Information, we show the extent to which the population transfers to D0 during this time.) If the in-phase equal superposition of states is populated (θ = 45° and ϕ = 0°), the amplitudes of the three terms in eqs 2 and 3 are sin2(θ) = 0.5, cos2(θ) = 0.5, and sin(2θ) cos(ϕ) = 1, respectively (Figure 4). The initial electronic density 1/2|ψ0|2 + 1/2|ψ1|2 + Re(ψ*0 ψ1) is intermediate between quinoid and antiquinoid. The ensemble of trajectories does not move significantly along the GD (Figure 5c) and moves instead toward positive DC (Figure 5d). It follows the gradient of the 1 1 superposition 2 ⟨ψ0|∇I /e|ψ0⟩ + 2 ⟨ψ1|∇I /e|ψ1⟩ + ⟨ψ0|∇I /e|ψ1⟩: the first two terms of the gradient have similar magnitudes but opposite directions along the GD so they cancel out leaving the third term that is along the DC dominating the gradient. In other words, the ensemble of trajectories moves to the top of the diagram in Figure 2: the nuclear geometry undergoes a shear motion where the top carbon atom of the ring moves to the left and the bottom carbon to the right. The ensemble of trajectories therefore moves in an orthogonal manner to the two adiabatic dynamics. In this case, the ensemble of trajectories spreads by ∼5−11% (see Table 1). Note that in a full quantum calculation, the nuclear wave packets on each electronic state would have the freedom to move along the GD as well as along the DC. In the Ehrenfest method considered here, the motion of the nuclear trajectories along the DC may be exaggerated because of the cancellation of the adiabatic gradients. The dashed white lines in Figure 5 show the average evolution of the ensemble, while the solid white lines show the evolution of the unsampled trajectory, i.e., the trajectory initiated at the equilibrium geometry of the neutral species with no initial kinetic energy. (Note that looking at the average evolution is meaningful only because the ensemble does not split in the first 5 fs.) In Figure 6, the unsampled and averaged trajectories are plotted in the branching space for the three initial electronic wave functions. They are in good qualitative agreement for each of the three initial electronic states (ground,

We show below the result of choosing specific values of the angles θ and ϕ on the gradient and on the initial motion of the nuclei. The initial nuclear dynamics adjacent to a conical intersection can be controlled by choosing a superposition with a specific gradient. Before presenting the results, we should comment on their sensitivity to the method used in this study. One potential limitation of the Ehrenfest method is the lack of electronic decoherence: the electronic eigenstates populated share the same nuclear geometry by definition. This could lead to nonphysical asymptotic behavior, but we do not expect this to be a problem here as we are only interested in relatively short time scale dynamics. Although the initial motion along the DC (third term of eq 3) and its sensitivity to the relative phase ϕ between the two electronic states is a real effect, it could be magnified because of some cancellation between the two adiabatic gradients (first two terms of eq 3). One way to investigate this in the future would be to work with an exact factorization of the total wave function33−35 or to use a Hellertype Gaussian wave packet representation23,36,37 of the nuclear wave packet.

■

RESULTS AND DISCUSSION Nuclear Dynamics Initiated on the Ground State, First Excited State, and the In-Phase Equal Superposition of Ground and First Excited States of the Cation. In order to demonstrate the effect of a nonstationary electronic state on the nuclear dynamics, we first compare the nuclear dynamics initiated with three different electronic wave functions: (i) the adiabatic ground state of the cation |Ψt=0⟩ = |ψ0⟩ corresponding to the angle θ = 90° in eqs 1−3, (ii) an equal superposition of the two lowest-energy eigenstates of the cation in phase |Ψt=0⟩ = 1/√2(|ψ0⟩ + |ψ1⟩) corresponding to θ = 45° and ϕ = 0°, and (iii) the adiabatic first excited state of the cation |Ψt=0⟩ = |ψ1⟩ corresponding to θ = 0°. For each initial electronic state, 500 Ehrenfest trajectories were run with initial geometries and velocities sampled as described above. Figure 5 shows the time evolution of the ensemble of trajectories along the two directions of the branching space, the GD (on the left panels) and the DC (on the right panels), for the three initial electronic wave functions: the ground state of the cation (Figure 5a,b), the in-phase equal superposition of the ground and first excited states of the cation (Figure 5c,d), and the first excited state of the cation (Figure 5e,f). In general, all distributions spread as time evolves: the standard deviations of geometries at t = 0 and t = 5 fs are indicated in Table 1 for each ensemble of trajectories. If only the ground state of the cation is initially populated (θ = 90°), the amplitudes of the three terms in eqs 2 and 3 are sin2(θ) = 1, cos2(θ) = 0, and sin(2θ)cos(ϕ) = 0, respectively (Figure 4). The initial electronic density |ψ0|2 has a quinoid Table 1. Standard Deviations of Geometries along the GD and the DC for Each Ensemble of Nuclear Trajectories: Dynamics Initiated on the Ground State θ = 90°, the InPhase Equal Superposition of States θ = 45° and ϕ = 0°, and the First Excited State θ = 0° t = 5 fs coordinate

t=0

θ = 90°

θ = 45°

θ = 0°

GD DC

0.39 0.38

0.48 (+23%) 0.50 (+32%)

0.41 (+5%) 0.42 (+11%)

0.43 (+10%) 0.40 (+5%) E



Article

Figure 6. Nuclear trajectories in the branching space initiated with |ψ0⟩ in green, 1/√2(|ψ0⟩ + |ψ1⟩) in blue, and |ψ1⟩ in red. The dashed lines show the average evolution for the ensemble of trajectories. The solid lines show the evolution of the unsampled trajectory starting at the equilibrium geometry of the neutral species with no initial kinetic energy. The time of simulation is 5 fs.

excited, and equal superposition): the unsampled trajectory gives the same direction of the nuclear motion in the branching space as the average evolution of the ensemble, but it tends to overestimate the amplitude of the nuclear motion. In other words, even though sampling is necessary for quantitative details of the nuclear dynamics, the basic behavior is already provided by the unsampled trajectory. These results show that ionization of toluene leading to an equal superposition of the two lowest energy states initiates nuclear dynamics in an orthogonal direction to dynamics on the ground or first excited state potential energy surfaces alone. Note that, in this particular case, observing this change in direction experimentally would be challenging due to the width of the distribution as shown in Figure 5. General Investigation of the Influence of an Electronic Superposition onto the Nuclear Motion. Having studied limiting cases and understood the effect of geometry sampling, we now systematically investigate the effect of varying the angles θ and ϕ in the electronic superposition. In the following, we have simulated only the unsampled trajectory to give the qualitative trend; the approximation is justified in the previous section. The time evolution of the nuclear trajectories in the branching space is shown in Figure 7 for different initial electronic wave functions. The projections onto the three 2D planes are also shown: GD−time and DC−time planes (equivalent to the left and right panels of Figure 5, respectively) and GD−DC plane (equivalent to Figure 6). The different colors correspond to different trajectories initiated with different angles in the initial electronic superposition (eq 1). In the upper panel (Figure 7a), all the initial electronic wave functions are in-phase superpositions of the ground and first excited states of the cation: ϕ = 0°. The relative weight between the two states varies from a pure ground state to a pure first excited state with all intermediate mixing: 0° ≤ θ ≤ 90°. According to Figure 4a, the amplitudes of the first two terms in

Figure 7. Nuclear trajectories in the branching space initiated with different initial electronic wave functions. (a) Initial electronic wave functions correspond to in-phase superpositions of the ground and first excited states ϕ = 0°, their relative weight varying with 0° ≤ θ ≤ 90°. (b) Initial electronic wave functions correspond to equal superpositions of the ground and first excited states θ = 45°, their relative phase varying with 0° ≤ ϕ ≤ 180°. The initial geometry is the equilibrium geometry of the neutral species, and the initial kinetic energy is null.

the gradient (eq 3) vary from 0 to 1 in an opposite manner, determining the component of the nuclear trajectory along the GD. At the same time, the amplitude of the third term in the gradient varies from 0 to 1 with a maximum at θ = 45° (Figure 4b), determining the component of the nuclear trajectory along the DC. In Figure 7a, the projection onto the GD−DC plane shows that the different trajectories evolve toward different directions in the branching space at approximately the same velocity. For θ = 90° (and 0°), the nuclear motion is the shortening (and lengthening, respectively) of the side carbon− carbon bonds. For θ = 45°, the nuclear motion is a shearing motion. For intermediate mixing angles, the nuclear geometry undergoes a combination of the shortening (lengthening) of the side bonds and the shearing motion. In the lower panel (Figure 7b), all the initial electronic wave functions are equal superpositions of the ground and first excited states of the cation: θ = 45°. The relative phase between the two states varies from an in-phase superposition to an outof-phase superposition with all intermediate phases: 0° ≤ ϕ ≤ 180°. According to Figure 4a, the amplitudes of the first two F



Article

terms in the gradient (eq 3) are both 0.5: as the first two components of the gradient have similar magnitude and opposite direction along the GD, they cancel out, and the component of the nuclear trajectory along the GD is negligible. At the same time, the amplitude of the third term in the gradient varies from 1 to −1 (Figure 4b), determining the component of the nuclear trajectory along the DC. In Figure 7b, the projection into the GD−DC plane shows that all the different trajectories evolve along the DC direction. The projections in the GD−time plane show barely any component along the GD during the first 4 or 5 fs. However, the trajectories move to either positive or negative DC and at different velocities. For instance, the nuclear trajectory initiated with ϕ = 0° (in turquoise) goes toward positive DC, which means a shearing motion with the top carbon atom of the ring moving to the left and the bottom carbon to the right; the trajectory initiated with ϕ = 180° (in pink) goes toward negative DC, which means a shearing motion with the top carbon atom of the ring moving to the right and the bottom carbon to the left. Note that the nuclear trajectory with ϕ = 90° feels an initial gradient that is null in the branching space (it moves only in the intersection space14 initially). The results show how both the relative weight and phase in the electronic superposition determine the initial nuclear gradient. By choosing the initial combination of electronic states, one can direct the initial nuclear motion in the branching space, i.e., how much of a shortening/lengthening motion along the GD and of a shearing distortion along the DC (leading sometimes to no initial nuclear motion in the branching space) takes place.

states. This material is available free of charge via the Internet at http://pubs.acs.org.

■

Corresponding Author

*E-mail: [email protected]. Present Addresses †

Institute for Theoretical Chemistry, University of Stuttgart, Germany. ‡ Laboratoire CEISAM-UMR CNR 6230, Université de Nantes, 2 Rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was supported by UK-EPSRC Grant EP/I032517/1. J.M. thanks the SimTech Cluster of Excellence for funding. All calculations were run using the Imperial College High Performance Computing service. The authors thank João Pedro Malhado, Lee Steinberg, Misha Ivanov, and Serguei Patchovskii for helpful discussions.

(1) Weinkauf, R.; Schlag, E. W.; Martinez, T. J.; Levine, R. D. Nonstationary Electronic States and Site-Selective Reactivity. J. Phys. Chem. A 1997, 101, 7702−7710. (2) Remacle, F.; Levine, R. D.; Schlag, E. W.; Weinkauf, R. Electronic Control of Site Selective Reactivity: A Model Combining Charge Migration and Dissociation. J. Phys. Chem. A 1999, 103, 10149−10158. (3) Yarkony, D. R. Diabolical Conical Intersections. Rev. Mod. Phys. 1996, 68, 985−1013. (4) Robb, M. A. In Conical Intersections: Theory, Computation and Experiment; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; World Scientific: Singapore, 2011; Chapter 1, pp 3−50. (5) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325−352. (6) von den Hoff, P.; Znakovskaya, I.; Kling, M. F.; de Vivie-Riedle, R. Attosecond Control of the Dissociative Ionization via Electron Localization: A Comparison Between D2 and CO. Chem. Phys. 2009, 366, 139−147. (7) von den Hoff, P.; Thallmair, S.; Kowalewski, M.; Siemering, R.; de Vivie-Riedle, R. Optimal Control Theory: Closing the Gap Between Theory and Experiment. Phys. Chem. Chem. Phys. 2012, 14, 14460− 14485. (8) Kling, M. F.; von den Hoff, P.; Znakovskaya, I.; de Vivie-Riedle, R. (Sub-)femtosecond Control of Molecular Reactions via Tailoring the Electric Field of Light. Phys. Chem. Chem. Phys. 2013, 15, 9448− 9467. (9) Mendive-Tapia, D.; Vacher, M.; Bearpark, M. J.; Robb, M. A. Coupled Electron-Nuclear Dynamics: Charge Migration and Charge Transfer Initiated near a Conical Intersection. J. Chem. Phys. 2013, 139, 044110. (10) Vacher, M.; Bearpark, M. J.; Robb, M. A. Communication: Oscillating Charge Migration Between Lone Pairs Persists Without Significant Interaction With Nuclear Motion in the Glycine and GlyGly-NH-CH3 Radical Cations. J. Chem. Phys. 2014, 140, 201102. (11) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian Development Version, revision H.10; Gaussian, Inc.: Wallingford, CT, 2010. (12) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106.

■

CONCLUSIONS Coupled electron−nuclear dynamics following ionization was simulated in toluene starting with different electronic wave functions. The initial geometries and velocities were sampled with a Wigner distribution to reproduce the spread of the equilibrium geometry as a result of zero-point motion of the neutral species. We have shown in the case of toluene that the ensemble of trajectories does not split within the first 5 fs. The average evolution is qualitatively well described by an unsampled trajectory (starting at the equilibrium geometry of the neutral species with no initial kinetic energy). The latter was then used for a more systematic investigation. Within the Ehrenfest approximation, our results for toluene show how the initial electronic wave function determines the initial nuclear dynamics: both the relative weight and phase of an electronic superposition of states play a role in the nuclear gradient. By choosing the initial electronic state, one can control the initial nuclear motion in the branching space, i.e., how much certain bond lengths are shortened or lengthened (left or right distortion along the GD) and how much of a shearing motion takes place (up or down distortion along the DC). For instance, if vertical ionization of toluene prepares an equal superposition of the two lowest energy states, the nuclear relaxation takes place in a direction orthogonal to the adiabatic dynamics, following the gradient of the superposition.

■

AUTHOR INFORMATION

ASSOCIATED CONTENT

S Supporting Information *

Equilibrium nuclear geometry of the neutral species. Evolution of the populations on the ground and first excited states of the cation during the dynamics initiated on the ground state, the first excited state, and the in-phase equal superposition of both G



Article

electronic density.38−43 To recover this feature within our formalism, one has to realize that the time-dependent electronic superposition defined in eq 1 can be expressed as |Ψ(t)⟩ = sin(θ(t))|ψ0⟩ + cos(θ(t)) ei(ϕ(t))|ψ1⟩. There is no transfer of population between the two states so θ(t) = θ(t0) is time-independent. The relative phase ϕ(t) = ϕ(t0) + (E0 − E1)t/ℏ is the sum of the initial relative phase between the two states and the time-dependent component that evolves linearly with time. The latter leads to alternating constructive and destructive interference. The period of the oscillations in the electronic density is inversely proportional to the energy gap between the two eigenstates. (32) Tully, J. C. In Modern Methods for Multidimensional Dynamics Computations in Chemistry; Thompson, D. L., Ed.; World Scientific: Singapore, 1998; pp 34−72. (33) Hunter, G. Conditional Probability Amplitudes in Wave Mechanics. Int. J. Quantum Chem. 1975, 9, 237−242. (34) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. Phys. Rev. Lett. 2010, 105, 123002. (35) Cederbaum, L. S. The Exact Molecular Wavefunction as a Product of an Electronic and a Nuclear Wavefunction. J. Chem. Phys. 2013, 138, 224110. (36) Heller, E. J. Time-Dependent Approach to Semiclassical Dynamics. J. Chem. Phys. 1975, 62, 1544−1555. (37) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab Initio Multiple Cloning Algorithm for Quantum Nonadiabatic Molecular Dynamics. J. Chem. Phys. 2014, 141, 054110. (38) Breidbach, J.; Cederbaum, L. S. Migration of Holes: Formalism, Mechanisms, and Illustrative Applications. J. Chem. Phys. 2003, 118, 3983−3996. (39) Remacle, F.; Levine, R. Probing Ultrafast Purely Electronic Charge Migration in Small Peptides. Z. Phys. Chem. 2007, 221, 647− 661. (40) Kuleff, A. I.; Cederbaum, L. S. Charge Migration in Different Conformers of Glycine: The Role of Nuclear Geometry. Chem. Phys. 2007, 338, 320−328. (41) Lünnemann, S.; Kuleff, A. I.; Cederbaum, L. S. Ultrafast Charge Migration in 2-Phenylethyl-N,N-dimethylamine. Chem. Phys. Lett. 2008, 450, 232−235. (42) Periyasamy, G.; Levine, R.; Remacle, F. Electronic Wave Packet Motion in Water Dimer Cation: A Many Electron Description. Chem. Phys. 2009, 366, 129−138. (43) Kuleff, A. I.; Lünnemann, S.; Cederbaum, L. S. Ultrafast Charge Migration Following Valence Ionization of 4-Methylphenol: Jumping over the Aromatic Ring. J. Phys. Chem. A 2010, 114, 8676−8679.

(13) Applegate, B. E.; Barckholtz, T. A.; Miller, T. A. Explorations of Conical Intersections and their Ramifications for Chemistry Through the Jahn-Teller Effect. Chem. Soc. Rev. 2003, 32, 38−49. (14) Paterson, M. J.; Bearpark, M. J.; Robb, M. A.; Blancafort, L.; Worth, G. A. Conical Intersections: A Perspective on the Computation of Spectroscopic Jahn-Teller Parameters and the Degenerate ‘Intersection Space’. Phys. Chem. Chem. Phys. 2005, 7, 2100−2115. (15) Hall, K. F.; Tokmachev, A. M.; Bearpark, M. J.; Boggio-Pasqua, M.; Robb, M. A. Molecular Mechanics-Valence Bond Method for Planar Conjugated Hydrocarbon Cations. J. Chem. Phys. 2007, 127, 134111. (16) Miller, T. A. Light and Radical Ions. Annu. Rev. Phys. Chem. 1982, 33, 257−282. (17) Whetten, R. L.; Fu, K.-J.; Grant, E. R. High Rydberg States of Jet-Cooled Toluene Observed by Ultraviolet Two-Photon Absorption Spectroscopy: Ultrafast Radiationless Decay and Pseudo-Jahn−Teller Effects. Chem. Phys. 1984, 90, 155−165. (18) Teller, E. The Crossing of Potential Surfaces. J. Phys. Chem. 1937, 41, 109−116. (19) Bersuker, I. B. Modern Aspects of the Jahn−Teller Effect Theory and Applications To Molecular Problems. Chem. Rev. 2001, 101, 1067−1114. (20) Longuet-Higgins, H. C.; Opik, U.; Pryce, M. H. L.; Sack, R. A. Studies of the Jahn-Teller Effect. II. The Dynamical Problem. Proc. R. Soc. London Ser. A 1958, 244, 1−16. (21) Döscher, M.; Köppel, H.; Szalay, P. G. Multistate Vibronic Interactions in the Benzene Radical Cation. I. Electronic Structure Calculations. J. Chem. Phys. 2002, 117, 2645−2656. (22) Köppel, H.; Döscher, M.; Bâldea, I.; Meyer, H.-D.; Szalay, P. G. Multistate Vibronic Interactions in the Benzene Radical Cation. II. Quantum Dynamical Simulations. J. Chem. Phys. 2002, 117, 2657− 2671. (23) Worth, G. A.; Robb, M. A.; Lasorne, B. Solving the TimeDependent Schrödinger Equation for Nuclear Motion in One Step: Direct Dynamics of Non-Adiabatic Systems. Mol. Phys. 2008, 106, 2077−2091. (24) Several electronic states can be populated coherently using for instance a pulse with a relatively large bandwidth. A coherent superposition of only the two lowest-energy eigenstates |ψ0⟩ and |ψ1⟩ is reasonable in the toluene cation since at the equilibrium geometry of the neutral species the energy gap between the ground and first excited states is small (0.2 eV at the CASSCF(5,6)/6-31G* level) compared to the energy gap between the first and second excited states (more than 2 eV). (25) In practice, the eigenstates are obtained by solving the timeindependent Schrödinger equation at a fixed nuclear geometry: in most quantum chemistry packages, they are real but their absolute sign is not uniquely defined. Therefore, when constructing a superposition of eigenstates where the relative phase matters, one must ensure consistency of the individual phases of each eigenstate among different simulations. (26) Barbatti, M.; Granucci, G.; Lischka, H.; Persico, M.; Ruckenbauer, M. Newton-X: A Package for Newtonian Dynamics Close to the Crossing Seam, Version 0.11b. www.univie.ac.at/newtonx. (27) Rybkin, V. V.; Ekström, U. Sampling Microcanonical Ensembles of Trajectories using Harmonic Approximation in Internal Coordinates. J. Chem. Phys. 2014, 141, 064108. (28) Vacher, M.; Mendive-Tapia, D.; Bearpark, M. J.; Robb, M. A. The Second-Order Ehrenfest Method. Theor. Chem. Acc. 2014, 133, 1505. (29) Millam, J. M.; Bakken, V.; Chen, W.; Hase, W. L.; Schlegel, H. B. Ab Initio Classical Trajectories on the Born−Oppenheimer Surface: Hessian-Based Integrators using Fifth-Order Polynomial and Rational Function Fits. J. Chem. Phys. 1999, 111, 3800−3805. (30) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. A Complete Active Space SCF Method (CASSCF) using a Density Matrix Formulated Super-CI Approach. Chem. Phys. 1980, 48, 157−173. (31) Theoretical studies of electron dynamics at fixed nuclear geometries in molecules have demonstrated oscillations in the H


Electronic Control of Initial Nuclear Dynamics Adjacent to a Conical

Recommend Documents