- Email: [email protected]
Quantum flux through a conical intersection
Chemical Physics Letters 372 (2003) 128–138 www.elsevier.com/locate/cplett
Quantum flux through a conical intersection Jingrui Li, Clemens Woywod
*
Institute of Physical and Theoretical Chemistry, Technical University of Munich, D-85747 Garching, Germany Received 20 December 2002; in final form 19 February 2003
Abstract Nonadiabatic dynamics at a conical intersection is investigated by solving the time-dependent Schr€ odinger equation numerically. The nuclear coordinate space is spanned by two degrees of freedom, a ÔtuningÕ and a ÔcouplingÕ coordinate. The analogy of a wavepacket with a quantum liquid is used to illustrate the nuclear motion in a strong vibronic coupling regime. An analysis of the time of evolution of the probability densities and velocity fields in the diabatic and adiabatic electronic picture reveals significant differences of the quantum flux between the two representations. The complexity of quantum currents in both electronic states increases rapidly with propagation time. Ó 2003 Elsevier Science B.V. All rights reserved.
1. Introduction The description of nuclear motion in molecules by quantum mechanical wavepacket simulations is in many ways similar to classical fluid dynamics. A moving wavepacket can therefore be referred to as a quantum liquid and the time-dependent Schr€ odinger equation can be transformed into two coupled differential equations, a continuity equation for the probability density and the so-called quantum Hamilton–Jacobi equation for the phase of the wavefunction [1–3]. This system of differential equations defines the framework for quantum fluid dynamics (QFD). The hydrodynamical equations can be solved either in the Eulerian or Lagrangian representation, i.e., with respect to a
*
Corresponding author. E-mail address: [email protected] (C. Woywod).
space-fixed or convected coordinate system. In the Eulerian picture, flux maps can be generated that reveal currents in the liquid while the Lagrangian point of view allows to track down the time evolution of individual fluid elements [4]. In the context of the treatment of dynamical processes in molecules, the ideas of quantum stream lines and trajectories have recently attracted significant interest because they provide a pictorial representation of quantum phenomena and may lead to new computational schemes for the solution of the Schr€ odinger equation [5–8]. Quantum flux maps have for example been employed to visualize the tunneling currents that accompany electron transfer in peptides [9] or local wavepacket motion in reactive scattering [10,11]. Studies involving the development of Lagrangian quantum trajectories have been carried out for problems such as transmission at an Eckart barrier [2,12] or rapid photodissociation [3,13].
0009-2614/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0009-2614(03)00378-6
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Solution of the QFD equations for nuclear motion on one or multiple potential energy surfaces requires a discretization of the wavepacket in coordinate space. The limitation imposed by representing the wavefunction in terms of a finite number of fluid elements, respectively, on a finite number of grid points constitutes the central approximation of QFD approaches [14]. A description of the motion of the quantum fluid in the Eulerian representation is complicated by the appearance of advective terms in the quantum Hamilton–Jacobi equation. Most theoretical studies have therefore been carried out in the Lagrangian framework and have been primarily devoted to dynamics in a single electronic state. In the Lagrangian view, trajectories are developed in analogy to classical fluid mechanics, but with the important difference that the quantum particles are correlated via the quantum potential and do not individually represent physical particles. The main disadvantage of the Lagrangian approach to QFD is the irregular distribution of trajectory endpoints that can result already after short propagation times even if the starting grid is regular and makes the numerical differentiation of the velocity field and of the wavefunction amplitude for an update of the density and of the particle positions, respectively, difficult [14–16]. Quantum trajectory studies of electronically nonadiabatic processes have been restricted to problems including only one nuclear degree of freedom so far. Several approximations to the full QFD equations have been examined in [17] by application to two-state avoided crossing problems. Unlike in the present work, the QFD formalism is developed by employing adiabatic electronic states in [17]. A so-called classical limit Schr€ odinger equation is derived by neglecting the quantum potential term in the quantum Hamilton–Jacobi equation. A further simplification is obtained by dropping two nonadiabatic coupling terms from the continuity equation for the density, leading to the so-called velocity coupling approximation. The approximate QFD models are shown to describe the population dynamics significantly better than the fewest-switches surface hopping scheme [18] for a dual avoided crossing example.
129
In [19,20], the QFD equations are formulated in the diabatic electronic picture (i.e., the working equations are identical to Eqs. (2.10) and (2.11)) and solved without making approximations beyond discretization for a nonadiabatic collision model involving two electronic states that are locally coupled. A comparison with the results of a conventional quantum mechanical grid calculation demonstrates the high numerical accuracy achieved by the Lagrangian quantum trajectory method. Burghardt and Cederbaum [21] further present a hydrodynamic formulation of nonadiabatic dynamics derived from the Liouville–von-Neumann equation for the density operator, considering the diabatic as well as the adiabatic electronic representation. In the present study, we will investigate the time dependence of the population transfer between two electronic states and of the quantum currents on both potential energy surfaces. The dynamical model we use is related to the S1 –S2 vibronic coupling problem in pyrazine [22,23] and can be described as a conical intersection spanned by one tuning and one coupling mode [24]. By analyzing the probability densities and velocity fields in the adiabatic and diabatic electronic pictures the nuclear motion is illuminated from different perspectives. Detailed information about the quantum flux through a conical intersection is of interest since electronic degeneracies are well known to play a key role for the origin of nonadiabatic dynamics in molecules [25].
2. Definition of model Hamiltonian The model Hamiltonian defined below describes wavepacket dynamics on two coupled electronic potential energy surfaces in a nuclear coordinate space spanned by two degrees of freedom. The electronic energies exhibit a symmetry induced conical intersection in this two-dimensional space, that is, the coordinates correspond to a totally symmetric or tuning mode (Qx ) and to a nontotally symmetric or coupling vibration (Qy ). Nuclear motion along Qx modulates the energy difference between the noninteracting electronic
130
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
states and displacement along Qy reduces the symmetry of the molecule and leads to a mixing of the two diabatic configurations in the adiabatic wavefunctions. The model is appropriate for the simulation of ultrafast internal conversion dynamics triggered by a strong vibronic interaction between the electronic states. Since we are interested mainly in the first 40 fs following preparation of the wavepacket in the S2 state, we will represent the vibrational motion in terms of normal coordinates, that is, we will neglect mode–mode couplings in the kinetic energy operator. The model Hamiltonian is constructed in the diabatic electronic picture to facilitate the numerical treatment of nonadiabatic transitions [24] and reads in atomic units ( h ¼ 1) H ¼ H^0 1 þ Vdia 0 E1 þ j1 Qx ¼ H^0 1 þ kQy
kQy ; E20 þ j2 Qx
where H^0 is defined as 2 1 o 1 2 H^0 ¼ xx þ Q þ xy x 2 2 oQ2x
ð2:1Þ
! o2 2 þ Qy : oQ2y ð2:2Þ
The system described by H can accordingly be characterized as a two-particle ensemble and each particle is limited to one-dimensional motion. The parameters used for the calculations are taken from a model developed for the S1 –S2 vibronic interaction problem in pyrazine [22] and are given in Table 1. The time-dependent Schr€ odinger equation dia o jWdia jW1 ðtÞi 1 ðtÞi i ¼ H ð2:3Þ ot jWdia jWdia 2 ðtÞi 2 ðtÞi for H is solved on a rectangular grid spanned by Qx and Qy . The vibronic wavefunction ðjWdia 1 ðtÞi; jWdia ðtÞiÞ at time t is obtained by applying the time 2
propagator UðtÞ ¼ eiHt to the initial wavepacket dia ðjWdia 1 ð0Þi; jW2 ð0ÞiÞ: dia dia jW1 ðtÞi jW1 ð0Þi ¼ UðtÞ : ð2:4Þ jWdia jWdia 2 ðtÞi 2 ð0Þi For short propagation times, the linear approximation UðtÞ 1 iHt is sufficiently accurate and has been chosen for the present study. The wavefunction at time t þ Dt is calculated by a finite difference scheme according to dia dia jW1 ðt þ DtÞi jW1 ðtÞi ð1 iHDtÞ : jWdia jWdia 2 ðt þ DtÞi 2 ðtÞi ð2:5Þ The action of the kinetic energy operator on the vibronic wavefunction is evaluated by the FFT technique [26]. The two-component wavefunction (jWdia 1 ðtÞi; dia jW2 ðtÞi) is defined with respect to a diabatic electronic basis and can be transformed into the adiabatic picture via [24] ! cos½hðQx ; Qy Þ sin½hðQx ; Qy Þ
jWadia 1 ðtÞi ¼ sin½hðQx ; Qy Þ cos½hðQx ; Qy Þ
jWadia 2 ðtÞi ! jWdia 1 ðtÞi : ð2:6Þ jWdia 2 ðtÞi The diabatic-to-adiabatic mixing angle hðQx ; Qy Þ is not time dependent since it is entirely determined by Vdia and is given by hðQx ; Qy Þ ¼
1 tan1 ½2V12dia =ðV22dia V11dia Þ : 2
For a graphical representation of hðQx ; Qy Þ, see [23]. We will assume in the present example that the oscillator strength is exclusively carried by jWdia 2 ðtÞi so that the starting wavepacket for the time propagation can be written as (0; jWdia 2 ð0Þi). However, due to the vibronic interaction, nonzero wavefunction amplitudes are obtained in both
Table 1 The parameters for the model Hamiltonian (Eq. (2.1)) are given belowa E10
E20
xx
xy
j1
j2
k
3.94
4.84
0.074
0.118
)0.105
0.149
0.262
a
ð2:7Þ
The values have been taken from a model for the S1 –S2 vibronic coupling in pyrazine [22]. Units are eV.
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
diabatic electronic states for t > 0. The adiabatic populations are different from zero for all times. At each time step, we can now define the velocity adia fields vdia n;x=y and vn;x=y describing the motion of the wavepackets on the diabatic and adiabatic potential energy surfaces, respectively: dia=adia
dia=adia
vn;x=y
ðtÞ ¼
jn;x=y
ðtÞ
dia=adia qn ðtÞ
! rx=y jWdia=adia ðtÞi 1 n ¼ Im ; mx=y jWdia=adia ðtÞi n n ¼ 1; 2;
ð2:8Þ
where density corresponding electronic state and representation. The reduced masses along the x- and y-directions are labeled as mx=y . Another expression for the velocity fields can be derived within the framework of quantum fluid dynamics (QFD). Inserting the polar ansatz [1] dia R1 ðtÞeiS1 ðtÞ jW1 ðtÞi ¼ ð2:9Þ R2 ðtÞeiS2 ðtÞ jWdia 2 ðtÞi with real valued amplitude (Rn ðtÞ) and phase functions (Sn ðtÞ) into Eq. (2.3) we obtain four coupled equations describing the time dependence of Sn ðtÞ and of the diabatic probability densities 2 qdia [19]. Wavepacket motion in n ðtÞ ¼ ðRn ðtÞÞ electronic state jWdia 1 ðtÞi with respect to a spacefixed coordinate frame (Euler picture) is described by: oS1 1 X dia 2 ¼ mk v1;k V11dia Q1 2 k¼x;y ot R2 dia V cos½S2 S1 ; R1 12
X oqdia 1 ¼ rk ðqvdia 1;k Þ ot k¼x;y pffiffiffiffiffiffiffiffiffi þ 2V12dia q1 q2 sin½S2 S1 ;
Nuclear dynamics in electronic state jWdia 2 ðtÞi follows by interchanging the indices 1 and 2 in Eqs. (2.10) and (2.11). The force Fn;x=y acting on a fluid element in electronic state n at a certain time and point in coordinate space can be evaluated by acting with the derivative operator rx=y on the quantum Hamilton–Jacobi equation (2.10) and subsequently switching to convected coordinates (Lagrange picture) d dia v dt 1;x=y dia ¼ rx=y V11 þ Q1 R2 dia rx=y V cos½S2 S1 : R1 12
F1;x=y ¼ mx=y
dia=adia jn;x=y ðtÞ denotes the quantum current and qdia=adia ðtÞ the probability density in the n
131
ð2:13Þ
Time integration of Eq. (2.13) yields the velocity dia field vdia 1;x=y ðtÞ. F2;x=y , and consequently v2;x=y ðtÞ, can again be obtained from Eq. (2.13) by permutation of the indices 1 and 2. 3. Results and discussion In Fig. 1, the evolution of the diabatic (solid line) and adiabatic (dotted line) electronic populations in the upper state is shown for the first 200 fs following preparation of the wavepacket in adia state jWdia 2 ðtÞi. The jW2 ðtÞi population is slightly smaller than 1 at time t ¼ 0 because of the configurational mixing induced by the zero point
ð2:10Þ
ð2:11Þ
where the quantum potential Q1 is defined according to Q1 ¼
X 1 r 2 R1 k : 2m R k 1 k¼x;y
ð2:12Þ
Fig. 1. Diabatic (solid curve) and adiabatic (dotted curve) population probability of the S2 state for the two-mode model for pyrazine.
132
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
amplitude of the coupling coordinate Qy . The activity of the coupling mode is also the reason for the larger asymptotic value of the diabatic S2 population as compared to the corresponding adiabatic limit because the jWdia 2 ðtÞi configuration contributes to the jWadia 1 ðtÞi wavefunction although the relaxation process is essentially complete. The transfer of adiabatic population occurs mainly within the first vibrational period, followed by a strong recurrence after 100 fs. In the diabatic picture, the oscillations of the population are a lot more pronounced, even leading to population inversion from 100 to 120 fs and again from 150 to 160 fs. For a more detailed discussion of population dynamics in pyrazine see [22]. We will now analyze the evolution of the vibrational wavepackets on both electronic potential energy surfaces in the diabatic and adiabatic electronic representations. Fig. 2 shows snapshots separated by time steps of 5, 10 and 20 fs of 2 2 dia jjWdia 2 ðtÞij (a)–(d) and of jjW1 ðtÞij (e)–(h) as a function of Qx and Qy on a grid ranging from )8 to 8 dimensionless normal coordinate units. In the coordinate system of Fig. 2, the ground state equilibrium geometry and the S1 –S2 intersection point are located at (Qx ¼ 0, Qy ¼ 0) and (Qx ¼ 3:54, Qy ¼ 0:0), respectively. The initial wavepacket representing the vibrational ground states of the tuning and coupling mode in the diabatic S2 state is centered at the reference geometry, no population probability is present on the S1 surface. Five femtoseconds later, the wavepacket (Fig. 2a) is accelerated along the negative Qx axis towards the minimum of the diabatic S2 state with coordinates (Qx ¼ 1:0, Qy ¼ 0:0) and a small amount of population has already been transferred onto the S1 surface (Fig. 2e). The wavepacket in the S1 state splits symmetrically with respect to the plane defined by Qy ¼ 0 because of the double well topology of the adiabatic S1 potential energy as a function of Qy (Fig. 2e). The local motion of the quantum liquid after 5 fs is shown in Figs. 3a and e. The orientation of the vector field in S2 along the Qx axis is quite uniform in the interval 6 < Qx < 3 (Fig. 3a). The acceleration of the main portion of the wavepacket in the negative Qx direction leads to the accruement of a boundary at which the orientation of the
velocity field rapidly changes. At the boundary, nodes of the wavefunction induce vortexes at which the vectors switch from pointing in the direction of the minimum energy path on the diabatic S2 surface to a distribution that is normal with respect to concentric circles of a radius up to 6 around the center of the wavepacket. For larger radii, at which the population is close to zero, the velocity field is essentially unstructured. In this region, the wavefunction amplitudes have become too small for an accurate numerical differentiation so that the orientation of the flux vectors cannot be unambiguously determined. The transition of the vector field from ordered to random occurs at about the same circular border in the S1 state (Fig. 3e). The velocity vectors point largely radially outwards except for two regions. First, we note the development of vortexes around two nodes located adjacent to the wavepacket maxima. Second, a triangular domain of quantum flux defined by two straight lines connecting the origin with the points ðQx ¼ 8:0; Qy 2:5Þ is formed. The two lines indicate a sharp reorientation of the quantum trajectories since they are aligned strictly antiparallel to the Qx axis within the triangle. The uphill motion of the S1 wavepacket in this triangular region demonstrates that the momentum of the quantum fluid on the S2 surface has been transferred to the S1 state. After 10 fs, the gradient on the S2 surface has driven the center of the wavepacket towards the crossing region and the width of the wavepacket along Qx has increased (Fig. 2b). More probability density has further accumulated in the S1 state (Fig. 2f). The velocity maps show that the structured regions have increased and taken an elliptic shape nearly identical in both electronic states, showing the expansion of the wavepacket. The motion of the quantum fluid on the S2 (Fig. 3b) surface is characterized by a continuing stream towards the S2 minimum and beyond to the S1 –S2 degeneracy seam. The increasing potential energy for displacements Qx < 1:0 leads to a deflection of the quantum trajectories away from the Qx axis and to the appearance of additional vortexes at the border separating velocity vectors that are oriented parallel to the Qx axis from rays that are directed
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138 Fig. 2. Diabatic probability densities for the S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
133
134
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Fig. 3. Velocity fields for the wavepacket in the diabatic S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
outwards. In the S1 (Fig. 3f) state, the motion along a positive energy gradient parallel to the negative Qx axis is still present, but the Qy components of the flow vectors are increasing. Fluid elements in the rear portion of the wavepacket are accelerated towards the minimum of the diabatic S1 surface located at ðQx ¼ 0:71; Qy ¼ 0:0Þ. This antidromic pattern induces angular momentum in the quantum fluid, leading to a reinforcement of the already present vortexes. Five femtoseconds later, the S2 wavepacket is passing over the conical intersection. The vibronic interaction reaches a maximum during this passage and a significant amount of population is built up in the S1 state accordingly. Evolution of the velocity fields between 10 fs and 15 fs leads to the simultaneous formation of two satellite vortexes at ðQx 0:5; Qy 2:5Þ in both electronic states. The turning point of the S2 wavepacket is reached after a total time of 20 fs has elapsed as is evident from the reorientation of the vector field along the negative Qx axis close to the edge of the
grid (Fig. 3c). The entire pattern of motion in both the S1 and S2 states is undergoing a pronounced change from 15 to 20 fs, manifest for example in the increasing number of points on the Qx axis at which the quantum flux alternates direction. On both diabatic potential surfaces, the quantum particles representing the main portion of the wavepacket are moving outwards in both directions along the Qy axis and the outermost fluid elements begin to bend towards the minimum of the respective electronic state (Figs. 3g,c). The snapshots taken at t ¼ 40 fs demonstrate that on each potential surface the two wavepacket fragments describe rotational motion and are in the process of converging back to the Qx axis (Figs. 3d,h). A comparison of Figs. 2b,f with 2c,g shows, that the wavepackets in the S1 and in particular in the S2 state are stretched along the Qx axis after a propagation time of 20 fs. During the next 20 fs, the probability density on both potential energy surfaces migrates in the positive Qx direction and, mainly in the S1 state, moves outwards along the
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138 Fig. 4. Adiabatic probability densities for the S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
135
136
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Qy axis. A total propagation time of 40 fs therefore renders both wavepackets quite delocalized (Figs. 2d,h). The adiabatic picture will be considered next. Fig. 4 illustrates the time evolution of the popu2 2 adia lations jjWadia 1 ðtÞij (e)–(h) and of jjW2 ðtÞij (a)– (d). At t ¼ 0 (not shown here), the diabatic and adiabatic wavepackets in S2 are nearly indistinguishable, but a very small probability amplitude is already present in the adiabatic S1 state. 2 jjWdia 1 ð0Þij is strictly zero, on the other hand, as discussed above. Within the next 5 fs, the S2 wavepacket begins to move towards the conical intersection (Fig. 4a) and more population is transferred to S1 (Fig. 4e). The adiabatic and diabatic quantum fluxes in S1 and S2 after 5 fs differ significantly as can be seen from a comparison of Figs. 5a,e and 3a,e. The adia velocity fields for jWadia 1 ðtÞi and jW2 ðtÞi are mainly oriented parallel to the energy gradients on the adiabatic double cone potential surfaces and no vortexes are present in contrast to the di-
abatic representation. Accordingly, the vibrational wavefunctions do not exhibit any nodes (see below). This changes over the next 5 fs, which lead to the formation of vortexes in both electronic states so that the appearance of the adiabatic flux maps largely resembles the diabatic quantum currents (Figs. 5b,f and 3b,f). The acceleration of the center of the S2 wavepacket to the intersection point becomes evident from Fig. 5b. In the S1 state, the nuclear motion is generally directed away from the conical intersection which represents a local maximum of the S1 potential energy function (Fig. 5f). Interestingly, no quantum trajectories are aligned parallel to the Qy axis. One fraction of the velocity vectors is characterized by a significant component along the negative Qx axis. This orientation is imposed by momentum transfer from the S2 state. The other part of the wavepacket begins to move along the two relaxation channels formed on the adiabatic S1 potential surface. At t ¼ 20 fs, the nuclei in the adiabatic S1 state are still moving away from the conical intersection
Fig. 5. Velocity fields for the wavepacket in the adiabatic S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
(Fig. 5g). The velocity map closely resembles the picture obtained in the diabatic representation (Fig. 3g). The driving force for the splitting of the wavepacket is in both cases energy relaxation as the positions of the minima of the adiabatic and diabatic S1 states are quite similar (ðQx ¼ 1:4; Qy ¼ 0:0Þ vs. ðQx ¼ 0:71; Qy ¼ 0:0Þ). The topology of the adiabatic and diabatic energy landscapes in the S2 state, however, can be characterized as funnel-like and harmonic, respectively, and this qualitative difference leads to the greatly deviating flux maps seen in Figs. 3c and 5c. The conical structure of the adiabatic S2 well directs the nuclei to the intersection point, while the wavepacket on the diabatic S2 surface is slowed down after passing through the minimum of the S2 state and bifurcates to eventually move backwards, describing two half circles as can be nicely seen after 40 fs (Fig. 3d). Fig. 5d demonstrates that only small amplitude motion is possible in the S2 state adiabatically. The S2 wavepacket approaches (10 fs, Fig. 4b) and reaches the conical intersection after 15 fs, leading to a rapid depletion and buildup of population in the upper and lower state, respectively. After a total propagation time of 15 fs, the adiabatic S1 population is localized at the degeneracy point and not split symmetrically along Qy as in the diabatic case. This reflects the pronounced mixing of the diabatic configurations in the adiabatic electronic wavefunction in the vicinity of the crossing point. Population is continuously leaking from S2 to S1 until the turning point along Qx is reached at t ¼ 20 fs (Fig. 4c). The adiabatic wavepacket on the S1 surface (20 fs, Fig. 4g) splits shortly afterwards and is found in a delocalized distribution at t ¼ 40 fs (Fig. 4h) very close to the diabatic picture (Fig. 2g). The adiabatic S2 population is still localized (Fig. 4d), however, unlike in the diabatic scenario (Fig. 2d). The trapping of the adiabatic upper state wavepacket at the bottom of the inverted cone is typical for internal conversion dynamics and persists as the propagation proceeds. For longer propagation times, the S2 population spreads out and decreases in the diabatic representation until the second strong recurrence at t ¼ 120 fs, a period corresponding to 2xx , (cf. Fig. 1) recreates a localized wavepacket in the
137
upper state. The diabatic S1 state becomes increasingly delocalized over the same period. The degree of delocalization of the S1 state in the diabatic and adiabatic pictures is of comparable magnitude for a time evolution beyond 20 fs, except for the localization of a part of the adiabatic wavepacket at the conical intersection after 120 fs that is associated with the repopulation of the S2 state (cf. Fig. 1).
4. Conclusions The quantum flux at the conical intersection is found to be characterized by a rapid development of complicated patterns of local motion in the diabatic as well as adiabatic electronic representation. The formation of an increasing number of vortexes in both electronic states as the propagation proceeds reflects the generation of nodes of the vibrational wavefunctions [9,27]. The quantum trajectories rotate around these nodes which leads to a quantization of the velocity field at each vortex. The location of the vortexes actually reveals the position of the nodes which are difficult to see from the plots of the probability density. The appearance of vortexes in quantum stream lines is a phenomenon also known from other systems such as superfluid liquids [28,29], superconductors, probability flux maps for chemical reactions [2,10,11], electron tunneling flows in polypeptide chains [9] and others [30]. At places of zero probability amplitude, the phase S of the wavefunction and accordingly the velocity field v are not defined since v ¼ rS=m. In the vicinity of a vortex, the orientation of the velocity field is tangential to the nodal lines and the absolute magnitude of the velocity is inversely proportional to the distance from the center of the vortex [27]. The radial and even more the angular velocity gradients are consequently large in the vicinity of a node. The eigenvalues of the local stability matrix rv will be large near a wavefunction node so that the orbits have a tendency to diverge when approaching a vortex. Since extreme sensitivity of trajectories to the initial conditions is the hallmark of chaos, vortexes can be considered as the most important sources of chaotic behavior in QFD [31,32].
138
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
We have considered here the time dependence of the quantum flux in the Eulerian picture, i.e., the motion of the fluid elements has been observed by stationary monitors. The Eulerian point of view does not reveal, however, if the trajectories described by individual quantum particles are periodic or aperiodic and possibly chaotic. Parmenter and Valentine have derived minimal conditions for the existence of chaotic trajectories in a quantum system. According to [33], the Hamiltonian must include at least two degrees of freedom, the timedependent wavefunction must be a linear combination of at least three stationary states and at least one pair of these states must have mutually incommensurate eigenvalues. According to these requirements, the present vibronic coupling model could in principle lead to quantum trajectories that follow chaotic orbits in both electronic states for certain choices of the parameters. This does not mean, however, that deterministic chaos would occur in the wavefunction W itself since the Schr€ odinger equation (2.3) is linear in W. The dynamics of W can only become chaotic if the Hamiltonian is a functional of the state of the system [34]. In order to trace the paths taken by the components of the discretized wavepacket, the motion of the quantum fluid must be described in the Lagrangian point of view. An investigation of these orbits at a conical intersection, in particular with respect to the possible emergence of quantum chaos and geometric phase effects on quantum trajectories [35], will be deferred to a future study. Acknowledgements We gratefully acknowledge support from the DFG. References [1] P.R. Holland, The Quantum Theory of Motion, Cambridge, New York, 1993.
[2] C.L. Lopreore, R.E. Wyatt, Phys. Rev. Lett. 82 (1999) 5190. [3] B.K. Dey, A. Askar, H. Rabitz, J. Chem. Phys. 109 (1998) 8770. [4] N.P. Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, Gulf Publ., Houston, 1990. [5] R.E. Wyatt, J. Chem. Phys. 117 (2002) 9569. [6] E. Gindensperger, C. Meier, J.A. Beswick, J. Chem. Phys. 113 (2000) 9369. [7] R.G. Brook, P.E. Oppenheimer, C.A. Weatherford, I. Banicescu, J. Zhu, J. Chem. Phys. 85 (2001) 263. [8] O.V. Prezhdo, C. Brooksby, Phys. Rev. Lett. 86 (2001) 3215. [9] I. Daizadeh, J. Guo, A. Stuchebrukhov, J. Chem. Phys. 110 (1999) 8865. [10] E.A. McCullough, R.E. Wyatt, J. Chem. Phys. 54 (1971) 3578. [11] R.E. Wyatt, J. Chem. Phys. 111 (1999) 4406. [12] D. Nerukh, J.H. Frederick, Chem. Phys. Lett. 332 (2000) 145. [13] F. Sales-Mayor, A. Askar, H. Rabitz, J. Chem. Phys. 111 (1999) 2423. [14] R.E. Wyatt, E.R. Bittner, J. Chem. Phys. 113 (2000) 8898. [15] X.-G. Hu, T.-S. Ho, H. Rabitz, Phys. Rev. E 61 (2000) 5967. [16] E.R. Bittner, R.E. Wyatt, J. Chem. Phys. 113 (2000) 8888. [17] J.C. Burant, J.C. Tully, J. Chem. Phys. 112 (2000) 6097. [18] J.C. Tully, J. Chem. Phys. 93 (1990) 1061. [19] R.E. Wyatt, C.L. Lopreore, G. Parlant, J. Chem. Phys. 114 (2001) 5113. [20] C.L. Lopreore, R.E. Wyatt, J. Chem. Phys. 116 (2002) 1228. [21] I. Burghardt, L.S. Cederbaum, J. Chem. Phys. 115 (2001) 10312. [22] R. Schneider, W. Domcke, H. K€ oppel, J. Chem. Phys. 92 (1990) 1045. [23] C. Woywod, W. Domcke, A.L. Sobolewski, H.-J. Werner, J. Chem. Phys. 100 (1994) 1400. [24] H. K€ oppel, W. Domcke, L.S. Cederbaum, Adv. Chem. Phys. 57 (1984) 59. [25] M. Klessinger, J. Michl, Excited States and Photochemistry of Organic Molecules, VCH, New York, 1995. [26] D. Kosloff, R. Kosloff, J. Comput. Phys. 52 (1983) 35. [27] H. Wu, D.W.L. Sprung, Phys. Lett. A 261 (1999) 150. [28] L. Onsager, Nuovo Cimento 6 (Suppl. 2) (1949) 249. [29] R.P. Feynman, Phys. Rev. 94 (1954) 262. [30] J.O. Hirschfelder, A.C. Christoph, W.E. Palke, J. Chem. Phys. 61 (1974) 5435. [31] H. Frisk, Phys. Lett. A 227 (1997) 139. [32] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 227 (1997) 5. [33] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 201 (1995) 1. [34] R.H. Parmenter, L.Y. Yu, Physica D 80 (1995) 289. [35] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 219 (1996) 7.
Quantum flux through a conical intersection Jingrui Li, Clemens Woywod
*
Institute of Physical and Theoretical Chemistry, Technical University of Munich, D-85747 Garching, Germany Received 20 December 2002; in final form 19 February 2003
Abstract Nonadiabatic dynamics at a conical intersection is investigated by solving the time-dependent Schr€ odinger equation numerically. The nuclear coordinate space is spanned by two degrees of freedom, a ÔtuningÕ and a ÔcouplingÕ coordinate. The analogy of a wavepacket with a quantum liquid is used to illustrate the nuclear motion in a strong vibronic coupling regime. An analysis of the time of evolution of the probability densities and velocity fields in the diabatic and adiabatic electronic picture reveals significant differences of the quantum flux between the two representations. The complexity of quantum currents in both electronic states increases rapidly with propagation time. Ó 2003 Elsevier Science B.V. All rights reserved.
1. Introduction The description of nuclear motion in molecules by quantum mechanical wavepacket simulations is in many ways similar to classical fluid dynamics. A moving wavepacket can therefore be referred to as a quantum liquid and the time-dependent Schr€ odinger equation can be transformed into two coupled differential equations, a continuity equation for the probability density and the so-called quantum Hamilton–Jacobi equation for the phase of the wavefunction [1–3]. This system of differential equations defines the framework for quantum fluid dynamics (QFD). The hydrodynamical equations can be solved either in the Eulerian or Lagrangian representation, i.e., with respect to a
*
Corresponding author. E-mail address: [email protected] (C. Woywod).
space-fixed or convected coordinate system. In the Eulerian picture, flux maps can be generated that reveal currents in the liquid while the Lagrangian point of view allows to track down the time evolution of individual fluid elements [4]. In the context of the treatment of dynamical processes in molecules, the ideas of quantum stream lines and trajectories have recently attracted significant interest because they provide a pictorial representation of quantum phenomena and may lead to new computational schemes for the solution of the Schr€ odinger equation [5–8]. Quantum flux maps have for example been employed to visualize the tunneling currents that accompany electron transfer in peptides [9] or local wavepacket motion in reactive scattering [10,11]. Studies involving the development of Lagrangian quantum trajectories have been carried out for problems such as transmission at an Eckart barrier [2,12] or rapid photodissociation [3,13].
0009-2614/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0009-2614(03)00378-6
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Solution of the QFD equations for nuclear motion on one or multiple potential energy surfaces requires a discretization of the wavepacket in coordinate space. The limitation imposed by representing the wavefunction in terms of a finite number of fluid elements, respectively, on a finite number of grid points constitutes the central approximation of QFD approaches [14]. A description of the motion of the quantum fluid in the Eulerian representation is complicated by the appearance of advective terms in the quantum Hamilton–Jacobi equation. Most theoretical studies have therefore been carried out in the Lagrangian framework and have been primarily devoted to dynamics in a single electronic state. In the Lagrangian view, trajectories are developed in analogy to classical fluid mechanics, but with the important difference that the quantum particles are correlated via the quantum potential and do not individually represent physical particles. The main disadvantage of the Lagrangian approach to QFD is the irregular distribution of trajectory endpoints that can result already after short propagation times even if the starting grid is regular and makes the numerical differentiation of the velocity field and of the wavefunction amplitude for an update of the density and of the particle positions, respectively, difficult [14–16]. Quantum trajectory studies of electronically nonadiabatic processes have been restricted to problems including only one nuclear degree of freedom so far. Several approximations to the full QFD equations have been examined in [17] by application to two-state avoided crossing problems. Unlike in the present work, the QFD formalism is developed by employing adiabatic electronic states in [17]. A so-called classical limit Schr€ odinger equation is derived by neglecting the quantum potential term in the quantum Hamilton–Jacobi equation. A further simplification is obtained by dropping two nonadiabatic coupling terms from the continuity equation for the density, leading to the so-called velocity coupling approximation. The approximate QFD models are shown to describe the population dynamics significantly better than the fewest-switches surface hopping scheme [18] for a dual avoided crossing example.
129
In [19,20], the QFD equations are formulated in the diabatic electronic picture (i.e., the working equations are identical to Eqs. (2.10) and (2.11)) and solved without making approximations beyond discretization for a nonadiabatic collision model involving two electronic states that are locally coupled. A comparison with the results of a conventional quantum mechanical grid calculation demonstrates the high numerical accuracy achieved by the Lagrangian quantum trajectory method. Burghardt and Cederbaum [21] further present a hydrodynamic formulation of nonadiabatic dynamics derived from the Liouville–von-Neumann equation for the density operator, considering the diabatic as well as the adiabatic electronic representation. In the present study, we will investigate the time dependence of the population transfer between two electronic states and of the quantum currents on both potential energy surfaces. The dynamical model we use is related to the S1 –S2 vibronic coupling problem in pyrazine [22,23] and can be described as a conical intersection spanned by one tuning and one coupling mode [24]. By analyzing the probability densities and velocity fields in the adiabatic and diabatic electronic pictures the nuclear motion is illuminated from different perspectives. Detailed information about the quantum flux through a conical intersection is of interest since electronic degeneracies are well known to play a key role for the origin of nonadiabatic dynamics in molecules [25].
2. Definition of model Hamiltonian The model Hamiltonian defined below describes wavepacket dynamics on two coupled electronic potential energy surfaces in a nuclear coordinate space spanned by two degrees of freedom. The electronic energies exhibit a symmetry induced conical intersection in this two-dimensional space, that is, the coordinates correspond to a totally symmetric or tuning mode (Qx ) and to a nontotally symmetric or coupling vibration (Qy ). Nuclear motion along Qx modulates the energy difference between the noninteracting electronic
130
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
states and displacement along Qy reduces the symmetry of the molecule and leads to a mixing of the two diabatic configurations in the adiabatic wavefunctions. The model is appropriate for the simulation of ultrafast internal conversion dynamics triggered by a strong vibronic interaction between the electronic states. Since we are interested mainly in the first 40 fs following preparation of the wavepacket in the S2 state, we will represent the vibrational motion in terms of normal coordinates, that is, we will neglect mode–mode couplings in the kinetic energy operator. The model Hamiltonian is constructed in the diabatic electronic picture to facilitate the numerical treatment of nonadiabatic transitions [24] and reads in atomic units ( h ¼ 1) H ¼ H^0 1 þ Vdia 0 E1 þ j1 Qx ¼ H^0 1 þ kQy
kQy ; E20 þ j2 Qx
where H^0 is defined as 2 1 o 1 2 H^0 ¼ xx þ Q þ xy x 2 2 oQ2x
ð2:1Þ
! o2 2 þ Qy : oQ2y ð2:2Þ
The system described by H can accordingly be characterized as a two-particle ensemble and each particle is limited to one-dimensional motion. The parameters used for the calculations are taken from a model developed for the S1 –S2 vibronic interaction problem in pyrazine [22] and are given in Table 1. The time-dependent Schr€ odinger equation dia o jWdia jW1 ðtÞi 1 ðtÞi i ¼ H ð2:3Þ ot jWdia jWdia 2 ðtÞi 2 ðtÞi for H is solved on a rectangular grid spanned by Qx and Qy . The vibronic wavefunction ðjWdia 1 ðtÞi; jWdia ðtÞiÞ at time t is obtained by applying the time 2
propagator UðtÞ ¼ eiHt to the initial wavepacket dia ðjWdia 1 ð0Þi; jW2 ð0ÞiÞ: dia dia jW1 ðtÞi jW1 ð0Þi ¼ UðtÞ : ð2:4Þ jWdia jWdia 2 ðtÞi 2 ð0Þi For short propagation times, the linear approximation UðtÞ 1 iHt is sufficiently accurate and has been chosen for the present study. The wavefunction at time t þ Dt is calculated by a finite difference scheme according to dia dia jW1 ðt þ DtÞi jW1 ðtÞi ð1 iHDtÞ : jWdia jWdia 2 ðt þ DtÞi 2 ðtÞi ð2:5Þ The action of the kinetic energy operator on the vibronic wavefunction is evaluated by the FFT technique [26]. The two-component wavefunction (jWdia 1 ðtÞi; dia jW2 ðtÞi) is defined with respect to a diabatic electronic basis and can be transformed into the adiabatic picture via [24] ! cos½hðQx ; Qy Þ sin½hðQx ; Qy Þ
jWadia 1 ðtÞi ¼ sin½hðQx ; Qy Þ cos½hðQx ; Qy Þ
jWadia 2 ðtÞi ! jWdia 1 ðtÞi : ð2:6Þ jWdia 2 ðtÞi The diabatic-to-adiabatic mixing angle hðQx ; Qy Þ is not time dependent since it is entirely determined by Vdia and is given by hðQx ; Qy Þ ¼
1 tan1 ½2V12dia =ðV22dia V11dia Þ : 2
For a graphical representation of hðQx ; Qy Þ, see [23]. We will assume in the present example that the oscillator strength is exclusively carried by jWdia 2 ðtÞi so that the starting wavepacket for the time propagation can be written as (0; jWdia 2 ð0Þi). However, due to the vibronic interaction, nonzero wavefunction amplitudes are obtained in both
Table 1 The parameters for the model Hamiltonian (Eq. (2.1)) are given belowa E10
E20
xx
xy
j1
j2
k
3.94
4.84
0.074
0.118
)0.105
0.149
0.262
a
ð2:7Þ
The values have been taken from a model for the S1 –S2 vibronic coupling in pyrazine [22]. Units are eV.
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
diabatic electronic states for t > 0. The adiabatic populations are different from zero for all times. At each time step, we can now define the velocity adia fields vdia n;x=y and vn;x=y describing the motion of the wavepackets on the diabatic and adiabatic potential energy surfaces, respectively: dia=adia
dia=adia
vn;x=y
ðtÞ ¼
jn;x=y
ðtÞ
dia=adia qn ðtÞ
! rx=y jWdia=adia ðtÞi 1 n ¼ Im ; mx=y jWdia=adia ðtÞi n n ¼ 1; 2;
ð2:8Þ
where density corresponding electronic state and representation. The reduced masses along the x- and y-directions are labeled as mx=y . Another expression for the velocity fields can be derived within the framework of quantum fluid dynamics (QFD). Inserting the polar ansatz [1] dia R1 ðtÞeiS1 ðtÞ jW1 ðtÞi ¼ ð2:9Þ R2 ðtÞeiS2 ðtÞ jWdia 2 ðtÞi with real valued amplitude (Rn ðtÞ) and phase functions (Sn ðtÞ) into Eq. (2.3) we obtain four coupled equations describing the time dependence of Sn ðtÞ and of the diabatic probability densities 2 qdia [19]. Wavepacket motion in n ðtÞ ¼ ðRn ðtÞÞ electronic state jWdia 1 ðtÞi with respect to a spacefixed coordinate frame (Euler picture) is described by: oS1 1 X dia 2 ¼ mk v1;k V11dia Q1 2 k¼x;y ot R2 dia V cos½S2 S1 ; R1 12
X oqdia 1 ¼ rk ðqvdia 1;k Þ ot k¼x;y pffiffiffiffiffiffiffiffiffi þ 2V12dia q1 q2 sin½S2 S1 ;
Nuclear dynamics in electronic state jWdia 2 ðtÞi follows by interchanging the indices 1 and 2 in Eqs. (2.10) and (2.11). The force Fn;x=y acting on a fluid element in electronic state n at a certain time and point in coordinate space can be evaluated by acting with the derivative operator rx=y on the quantum Hamilton–Jacobi equation (2.10) and subsequently switching to convected coordinates (Lagrange picture) d dia v dt 1;x=y dia ¼ rx=y V11 þ Q1 R2 dia rx=y V cos½S2 S1 : R1 12
F1;x=y ¼ mx=y
dia=adia jn;x=y ðtÞ denotes the quantum current and qdia=adia ðtÞ the probability density in the n
131
ð2:13Þ
Time integration of Eq. (2.13) yields the velocity dia field vdia 1;x=y ðtÞ. F2;x=y , and consequently v2;x=y ðtÞ, can again be obtained from Eq. (2.13) by permutation of the indices 1 and 2. 3. Results and discussion In Fig. 1, the evolution of the diabatic (solid line) and adiabatic (dotted line) electronic populations in the upper state is shown for the first 200 fs following preparation of the wavepacket in adia state jWdia 2 ðtÞi. The jW2 ðtÞi population is slightly smaller than 1 at time t ¼ 0 because of the configurational mixing induced by the zero point
ð2:10Þ
ð2:11Þ
where the quantum potential Q1 is defined according to Q1 ¼
X 1 r 2 R1 k : 2m R k 1 k¼x;y
ð2:12Þ
Fig. 1. Diabatic (solid curve) and adiabatic (dotted curve) population probability of the S2 state for the two-mode model for pyrazine.
132
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
amplitude of the coupling coordinate Qy . The activity of the coupling mode is also the reason for the larger asymptotic value of the diabatic S2 population as compared to the corresponding adiabatic limit because the jWdia 2 ðtÞi configuration contributes to the jWadia 1 ðtÞi wavefunction although the relaxation process is essentially complete. The transfer of adiabatic population occurs mainly within the first vibrational period, followed by a strong recurrence after 100 fs. In the diabatic picture, the oscillations of the population are a lot more pronounced, even leading to population inversion from 100 to 120 fs and again from 150 to 160 fs. For a more detailed discussion of population dynamics in pyrazine see [22]. We will now analyze the evolution of the vibrational wavepackets on both electronic potential energy surfaces in the diabatic and adiabatic electronic representations. Fig. 2 shows snapshots separated by time steps of 5, 10 and 20 fs of 2 2 dia jjWdia 2 ðtÞij (a)–(d) and of jjW1 ðtÞij (e)–(h) as a function of Qx and Qy on a grid ranging from )8 to 8 dimensionless normal coordinate units. In the coordinate system of Fig. 2, the ground state equilibrium geometry and the S1 –S2 intersection point are located at (Qx ¼ 0, Qy ¼ 0) and (Qx ¼ 3:54, Qy ¼ 0:0), respectively. The initial wavepacket representing the vibrational ground states of the tuning and coupling mode in the diabatic S2 state is centered at the reference geometry, no population probability is present on the S1 surface. Five femtoseconds later, the wavepacket (Fig. 2a) is accelerated along the negative Qx axis towards the minimum of the diabatic S2 state with coordinates (Qx ¼ 1:0, Qy ¼ 0:0) and a small amount of population has already been transferred onto the S1 surface (Fig. 2e). The wavepacket in the S1 state splits symmetrically with respect to the plane defined by Qy ¼ 0 because of the double well topology of the adiabatic S1 potential energy as a function of Qy (Fig. 2e). The local motion of the quantum liquid after 5 fs is shown in Figs. 3a and e. The orientation of the vector field in S2 along the Qx axis is quite uniform in the interval 6 < Qx < 3 (Fig. 3a). The acceleration of the main portion of the wavepacket in the negative Qx direction leads to the accruement of a boundary at which the orientation of the
velocity field rapidly changes. At the boundary, nodes of the wavefunction induce vortexes at which the vectors switch from pointing in the direction of the minimum energy path on the diabatic S2 surface to a distribution that is normal with respect to concentric circles of a radius up to 6 around the center of the wavepacket. For larger radii, at which the population is close to zero, the velocity field is essentially unstructured. In this region, the wavefunction amplitudes have become too small for an accurate numerical differentiation so that the orientation of the flux vectors cannot be unambiguously determined. The transition of the vector field from ordered to random occurs at about the same circular border in the S1 state (Fig. 3e). The velocity vectors point largely radially outwards except for two regions. First, we note the development of vortexes around two nodes located adjacent to the wavepacket maxima. Second, a triangular domain of quantum flux defined by two straight lines connecting the origin with the points ðQx ¼ 8:0; Qy 2:5Þ is formed. The two lines indicate a sharp reorientation of the quantum trajectories since they are aligned strictly antiparallel to the Qx axis within the triangle. The uphill motion of the S1 wavepacket in this triangular region demonstrates that the momentum of the quantum fluid on the S2 surface has been transferred to the S1 state. After 10 fs, the gradient on the S2 surface has driven the center of the wavepacket towards the crossing region and the width of the wavepacket along Qx has increased (Fig. 2b). More probability density has further accumulated in the S1 state (Fig. 2f). The velocity maps show that the structured regions have increased and taken an elliptic shape nearly identical in both electronic states, showing the expansion of the wavepacket. The motion of the quantum fluid on the S2 (Fig. 3b) surface is characterized by a continuing stream towards the S2 minimum and beyond to the S1 –S2 degeneracy seam. The increasing potential energy for displacements Qx < 1:0 leads to a deflection of the quantum trajectories away from the Qx axis and to the appearance of additional vortexes at the border separating velocity vectors that are oriented parallel to the Qx axis from rays that are directed
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138 Fig. 2. Diabatic probability densities for the S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
133
134
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Fig. 3. Velocity fields for the wavepacket in the diabatic S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
outwards. In the S1 (Fig. 3f) state, the motion along a positive energy gradient parallel to the negative Qx axis is still present, but the Qy components of the flow vectors are increasing. Fluid elements in the rear portion of the wavepacket are accelerated towards the minimum of the diabatic S1 surface located at ðQx ¼ 0:71; Qy ¼ 0:0Þ. This antidromic pattern induces angular momentum in the quantum fluid, leading to a reinforcement of the already present vortexes. Five femtoseconds later, the S2 wavepacket is passing over the conical intersection. The vibronic interaction reaches a maximum during this passage and a significant amount of population is built up in the S1 state accordingly. Evolution of the velocity fields between 10 fs and 15 fs leads to the simultaneous formation of two satellite vortexes at ðQx 0:5; Qy 2:5Þ in both electronic states. The turning point of the S2 wavepacket is reached after a total time of 20 fs has elapsed as is evident from the reorientation of the vector field along the negative Qx axis close to the edge of the
grid (Fig. 3c). The entire pattern of motion in both the S1 and S2 states is undergoing a pronounced change from 15 to 20 fs, manifest for example in the increasing number of points on the Qx axis at which the quantum flux alternates direction. On both diabatic potential surfaces, the quantum particles representing the main portion of the wavepacket are moving outwards in both directions along the Qy axis and the outermost fluid elements begin to bend towards the minimum of the respective electronic state (Figs. 3g,c). The snapshots taken at t ¼ 40 fs demonstrate that on each potential surface the two wavepacket fragments describe rotational motion and are in the process of converging back to the Qx axis (Figs. 3d,h). A comparison of Figs. 2b,f with 2c,g shows, that the wavepackets in the S1 and in particular in the S2 state are stretched along the Qx axis after a propagation time of 20 fs. During the next 20 fs, the probability density on both potential energy surfaces migrates in the positive Qx direction and, mainly in the S1 state, moves outwards along the
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138 Fig. 4. Adiabatic probability densities for the S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
135
136
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
Qy axis. A total propagation time of 40 fs therefore renders both wavepackets quite delocalized (Figs. 2d,h). The adiabatic picture will be considered next. Fig. 4 illustrates the time evolution of the popu2 2 adia lations jjWadia 1 ðtÞij (e)–(h) and of jjW2 ðtÞij (a)– (d). At t ¼ 0 (not shown here), the diabatic and adiabatic wavepackets in S2 are nearly indistinguishable, but a very small probability amplitude is already present in the adiabatic S1 state. 2 jjWdia 1 ð0Þij is strictly zero, on the other hand, as discussed above. Within the next 5 fs, the S2 wavepacket begins to move towards the conical intersection (Fig. 4a) and more population is transferred to S1 (Fig. 4e). The adiabatic and diabatic quantum fluxes in S1 and S2 after 5 fs differ significantly as can be seen from a comparison of Figs. 5a,e and 3a,e. The adia velocity fields for jWadia 1 ðtÞi and jW2 ðtÞi are mainly oriented parallel to the energy gradients on the adiabatic double cone potential surfaces and no vortexes are present in contrast to the di-
abatic representation. Accordingly, the vibrational wavefunctions do not exhibit any nodes (see below). This changes over the next 5 fs, which lead to the formation of vortexes in both electronic states so that the appearance of the adiabatic flux maps largely resembles the diabatic quantum currents (Figs. 5b,f and 3b,f). The acceleration of the center of the S2 wavepacket to the intersection point becomes evident from Fig. 5b. In the S1 state, the nuclear motion is generally directed away from the conical intersection which represents a local maximum of the S1 potential energy function (Fig. 5f). Interestingly, no quantum trajectories are aligned parallel to the Qy axis. One fraction of the velocity vectors is characterized by a significant component along the negative Qx axis. This orientation is imposed by momentum transfer from the S2 state. The other part of the wavepacket begins to move along the two relaxation channels formed on the adiabatic S1 potential surface. At t ¼ 20 fs, the nuclei in the adiabatic S1 state are still moving away from the conical intersection
Fig. 5. Velocity fields for the wavepacket in the adiabatic S2 (a)–(d) and S1 (e)–(h) states are shown on a grid spanned by the normal coordinates Qx and Qy . Propagation times: 5, 10, 20 and 40 fs.
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
(Fig. 5g). The velocity map closely resembles the picture obtained in the diabatic representation (Fig. 3g). The driving force for the splitting of the wavepacket is in both cases energy relaxation as the positions of the minima of the adiabatic and diabatic S1 states are quite similar (ðQx ¼ 1:4; Qy ¼ 0:0Þ vs. ðQx ¼ 0:71; Qy ¼ 0:0Þ). The topology of the adiabatic and diabatic energy landscapes in the S2 state, however, can be characterized as funnel-like and harmonic, respectively, and this qualitative difference leads to the greatly deviating flux maps seen in Figs. 3c and 5c. The conical structure of the adiabatic S2 well directs the nuclei to the intersection point, while the wavepacket on the diabatic S2 surface is slowed down after passing through the minimum of the S2 state and bifurcates to eventually move backwards, describing two half circles as can be nicely seen after 40 fs (Fig. 3d). Fig. 5d demonstrates that only small amplitude motion is possible in the S2 state adiabatically. The S2 wavepacket approaches (10 fs, Fig. 4b) and reaches the conical intersection after 15 fs, leading to a rapid depletion and buildup of population in the upper and lower state, respectively. After a total propagation time of 15 fs, the adiabatic S1 population is localized at the degeneracy point and not split symmetrically along Qy as in the diabatic case. This reflects the pronounced mixing of the diabatic configurations in the adiabatic electronic wavefunction in the vicinity of the crossing point. Population is continuously leaking from S2 to S1 until the turning point along Qx is reached at t ¼ 20 fs (Fig. 4c). The adiabatic wavepacket on the S1 surface (20 fs, Fig. 4g) splits shortly afterwards and is found in a delocalized distribution at t ¼ 40 fs (Fig. 4h) very close to the diabatic picture (Fig. 2g). The adiabatic S2 population is still localized (Fig. 4d), however, unlike in the diabatic scenario (Fig. 2d). The trapping of the adiabatic upper state wavepacket at the bottom of the inverted cone is typical for internal conversion dynamics and persists as the propagation proceeds. For longer propagation times, the S2 population spreads out and decreases in the diabatic representation until the second strong recurrence at t ¼ 120 fs, a period corresponding to 2xx , (cf. Fig. 1) recreates a localized wavepacket in the
137
upper state. The diabatic S1 state becomes increasingly delocalized over the same period. The degree of delocalization of the S1 state in the diabatic and adiabatic pictures is of comparable magnitude for a time evolution beyond 20 fs, except for the localization of a part of the adiabatic wavepacket at the conical intersection after 120 fs that is associated with the repopulation of the S2 state (cf. Fig. 1).
4. Conclusions The quantum flux at the conical intersection is found to be characterized by a rapid development of complicated patterns of local motion in the diabatic as well as adiabatic electronic representation. The formation of an increasing number of vortexes in both electronic states as the propagation proceeds reflects the generation of nodes of the vibrational wavefunctions [9,27]. The quantum trajectories rotate around these nodes which leads to a quantization of the velocity field at each vortex. The location of the vortexes actually reveals the position of the nodes which are difficult to see from the plots of the probability density. The appearance of vortexes in quantum stream lines is a phenomenon also known from other systems such as superfluid liquids [28,29], superconductors, probability flux maps for chemical reactions [2,10,11], electron tunneling flows in polypeptide chains [9] and others [30]. At places of zero probability amplitude, the phase S of the wavefunction and accordingly the velocity field v are not defined since v ¼ rS=m. In the vicinity of a vortex, the orientation of the velocity field is tangential to the nodal lines and the absolute magnitude of the velocity is inversely proportional to the distance from the center of the vortex [27]. The radial and even more the angular velocity gradients are consequently large in the vicinity of a node. The eigenvalues of the local stability matrix rv will be large near a wavefunction node so that the orbits have a tendency to diverge when approaching a vortex. Since extreme sensitivity of trajectories to the initial conditions is the hallmark of chaos, vortexes can be considered as the most important sources of chaotic behavior in QFD [31,32].
138
J. Li, C. Woywod / Chemical Physics Letters 372 (2003) 128–138
We have considered here the time dependence of the quantum flux in the Eulerian picture, i.e., the motion of the fluid elements has been observed by stationary monitors. The Eulerian point of view does not reveal, however, if the trajectories described by individual quantum particles are periodic or aperiodic and possibly chaotic. Parmenter and Valentine have derived minimal conditions for the existence of chaotic trajectories in a quantum system. According to [33], the Hamiltonian must include at least two degrees of freedom, the timedependent wavefunction must be a linear combination of at least three stationary states and at least one pair of these states must have mutually incommensurate eigenvalues. According to these requirements, the present vibronic coupling model could in principle lead to quantum trajectories that follow chaotic orbits in both electronic states for certain choices of the parameters. This does not mean, however, that deterministic chaos would occur in the wavefunction W itself since the Schr€ odinger equation (2.3) is linear in W. The dynamics of W can only become chaotic if the Hamiltonian is a functional of the state of the system [34]. In order to trace the paths taken by the components of the discretized wavepacket, the motion of the quantum fluid must be described in the Lagrangian point of view. An investigation of these orbits at a conical intersection, in particular with respect to the possible emergence of quantum chaos and geometric phase effects on quantum trajectories [35], will be deferred to a future study. Acknowledgements We gratefully acknowledge support from the DFG. References [1] P.R. Holland, The Quantum Theory of Motion, Cambridge, New York, 1993.
[2] C.L. Lopreore, R.E. Wyatt, Phys. Rev. Lett. 82 (1999) 5190. [3] B.K. Dey, A. Askar, H. Rabitz, J. Chem. Phys. 109 (1998) 8770. [4] N.P. Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, Gulf Publ., Houston, 1990. [5] R.E. Wyatt, J. Chem. Phys. 117 (2002) 9569. [6] E. Gindensperger, C. Meier, J.A. Beswick, J. Chem. Phys. 113 (2000) 9369. [7] R.G. Brook, P.E. Oppenheimer, C.A. Weatherford, I. Banicescu, J. Zhu, J. Chem. Phys. 85 (2001) 263. [8] O.V. Prezhdo, C. Brooksby, Phys. Rev. Lett. 86 (2001) 3215. [9] I. Daizadeh, J. Guo, A. Stuchebrukhov, J. Chem. Phys. 110 (1999) 8865. [10] E.A. McCullough, R.E. Wyatt, J. Chem. Phys. 54 (1971) 3578. [11] R.E. Wyatt, J. Chem. Phys. 111 (1999) 4406. [12] D. Nerukh, J.H. Frederick, Chem. Phys. Lett. 332 (2000) 145. [13] F. Sales-Mayor, A. Askar, H. Rabitz, J. Chem. Phys. 111 (1999) 2423. [14] R.E. Wyatt, E.R. Bittner, J. Chem. Phys. 113 (2000) 8898. [15] X.-G. Hu, T.-S. Ho, H. Rabitz, Phys. Rev. E 61 (2000) 5967. [16] E.R. Bittner, R.E. Wyatt, J. Chem. Phys. 113 (2000) 8888. [17] J.C. Burant, J.C. Tully, J. Chem. Phys. 112 (2000) 6097. [18] J.C. Tully, J. Chem. Phys. 93 (1990) 1061. [19] R.E. Wyatt, C.L. Lopreore, G. Parlant, J. Chem. Phys. 114 (2001) 5113. [20] C.L. Lopreore, R.E. Wyatt, J. Chem. Phys. 116 (2002) 1228. [21] I. Burghardt, L.S. Cederbaum, J. Chem. Phys. 115 (2001) 10312. [22] R. Schneider, W. Domcke, H. K€ oppel, J. Chem. Phys. 92 (1990) 1045. [23] C. Woywod, W. Domcke, A.L. Sobolewski, H.-J. Werner, J. Chem. Phys. 100 (1994) 1400. [24] H. K€ oppel, W. Domcke, L.S. Cederbaum, Adv. Chem. Phys. 57 (1984) 59. [25] M. Klessinger, J. Michl, Excited States and Photochemistry of Organic Molecules, VCH, New York, 1995. [26] D. Kosloff, R. Kosloff, J. Comput. Phys. 52 (1983) 35. [27] H. Wu, D.W.L. Sprung, Phys. Lett. A 261 (1999) 150. [28] L. Onsager, Nuovo Cimento 6 (Suppl. 2) (1949) 249. [29] R.P. Feynman, Phys. Rev. 94 (1954) 262. [30] J.O. Hirschfelder, A.C. Christoph, W.E. Palke, J. Chem. Phys. 61 (1974) 5435. [31] H. Frisk, Phys. Lett. A 227 (1997) 139. [32] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 227 (1997) 5. [33] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 201 (1995) 1. [34] R.H. Parmenter, L.Y. Yu, Physica D 80 (1995) 289. [35] R.H. Parmenter, R.W. Valentine, Phys. Lett. A 219 (1996) 7.