7896
J. Phys. Chem. 1996, 100, 7896-7902
A Wavepacket-Path Integral Method for Curve-Crossing Dynamics Rob D. Coalson Department of Chemistry, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15260 ReceiVed: October 30, 1995; In Final Form: December 14, 1995X
A combined wavepacket-path integral method is proposed for computing the time evolution of nuclear coordinate wavepackets on several nonradiatively coupled potential energy surfaces. In order to make the technique useful for multidimensional problems, the method uses the Gaussian wavepacket dynamics (GWD) algorithm (Heller, E. J.; J. Chem. Phys. 1975, 62, 1544) to perform nuclear coordinate wavepacket evolution on successive individual surfaces as prescribed by the path integral formalism. Since the GWD algorithm is not exact for anharmonically coupled potential energy surfaces, this introduces some error into the calculation. However, as long as the GWD algorithm is sufficiently accurate to treat single-surface propagation on all of the relevant potential energy surfaces, the proposed method is expected to yield useful results. A simple illustration is given using a one-dimensional model of a direct photodissociation process involving two nonradiatively coupled electronic states. Good agreement with numerically exact results is obtained, even in the case of strong, and position-dependent, nonradiative coupling.
1. Introduction Processes which involve transitions between nonradiatively coupled Born-Oppenheimer electronic states play a central role in many quantum dynamical processes ranging from excitedstate photochemistry1 to electron transfer.2 Of particular interest in many chemical and condensed matter physics applications are effects on the rates of these transitions (often referred to as “nonadiabatic transitions”) due to the presence of a condensedphase environment. In order to study quantum dynamics in many-body systems (a classification which includes most isolated polyatomic molecules as well), special numerical techniques are needed to solve the coupled multidimensional Schro¨dinger equations which arise in problems involving crossings between several potential energy curves. The prognosis for an exact solution of multidimensional Schro¨dinger equations is not encouraging. Due to the necessity of expanding the unknown wavefunction or wavepacket in a direct product basis (a product of single-coordinate basis functions, with all possible number of independent excitations in each coordinate included in the expansion), basis set and grid wavepacket techniques utilize matrices which scale exponentially with the number of nuclear coordinate degrees of freedom in the system. Thus, the state of the art in propagating nuclear coordinate wavepackets on general potential energy surfaces is 3-4 degrees of freedom,3 far from the capacity needed to treat condensed-phase systems. For propagation on coupled potential surfaces, the situation is similar to the single-surface case: for a fixed number of spatial degrees of freedom, basis set size scales linearly with the number of potential surfaces that are coupled together. So it is not the curve-crossing aspect but simply the high spatial dimensionality of condensed-phase problems that severely restricts the application of exact integration schemes. This state of affairs has spawned a great deal of effort to develop reliable many-body techniques for solving coupled nuclear-coordinate Schro¨dinger equations appropriate to nonadiabatic transition phenomena. Most are, by virture of the difficulties cited above, approximate.4-6 One major exception is the use of path integral techniques to provide numerically exact solutions to the Spin-Boson model,8 which has proved X
Abstract published in AdVance ACS Abstracts, April 1, 1996.
S0022-3654(95)03182-0 CCC: $12.00
relevant to chemically interesting nonadiabatic transition processes (see below).9-13 Each of the resultant techniques has its strengths and weaknesses. For example, surface hopping techniques are simple to implement, even in realistically complex systems, because they require only ensembles of classical trajectories, plus a probabilistic rule that governs hopping between potential surfaces.4 These have provided extremely useful guidelines for many nonadiabatic transition processes. However, there is always some concern about accuracy. There is no way to successively improve these techniques toward an exact solution of the Schro¨dinger equation (SE). Indeed, by construction, they suppress certain wavemechanical features of the SE; hence, effects such as phase interference that can emerge from exact solutions are lost at the outset. Time-dependent Hartree grid wavepacket methods have also proved valuable in studying many-body quantum dynamics where curve crossing is involved.5 In this procedure, the multidimensional nuclear coordinate wavepacket is factorized into one-dimensional pieces, each of which is propagated using an effective one-dimensional effective potential (which depends in a mean-field way on the instantaneous probability density associated with all the other degrees of freedom). Timedependent Hartree methods are fully wavemechanical so that complex phase effects, wavepacket spreading, and nuclear tunneling are included. Again, there is a problem of knowing how accurate the solutions are. It is possible to develop systematic configuration interaction correction schemes, but these require much more effort.14,15 Also, the factorization appropriate to an amorphous condensed-phase environment such as a liquid solvent is not obvious. Finally, among approximate techniques for computing multidimensional curve-crossing dynamics, let us mention here the Golden Rule for nonadiabatic transitions.6,7 This can be successfully implemented for many-dimensional potential surfaces in certain useful regimes (generally, high temperature and high density of accepting states). It also provides a simple firstorder kinetics picture of the probability flow between electronic states. Unfortunately, it breaks down when the nonradiative coupling is strong. The Spin-Boson model may be thought of as a Hamiltonian which consists of two multidimensional, linearly displaced © 1996 American Chemical Society
Wavepacket-Path Integral Method for Dynamics
J. Phys. Chem., Vol. 100, No. 19, 1996 7897
harmonic oscillator potential surfaces coupled via a constant nonadiabatic hopping matrix element.8-11 Using the path integral formalism, this problem can be transformed into one requiring sequential propagation of nuclear coordinate wavepackets on the two harmonic potential surfaces. That is, we must propagate a wavepacket for some time on surface 1, then propagate it on surface 2 for another time interval, then “hop” it back to surface 1, and so on. All possible sequences of hopping times must be considered, and the amplitudes associated with each path through the (two-level) electronic-state space summed coherently. (This procedure can be made precise, as is done below.) An essential requirement for the success of such a scheme is the ability to accurately and efficiently perform the nuclear coordinate wavepacket dynamics on either of the potential surfaces, 1 or 2. The Spin-Boson model, where both potential surfaces are sums of separable harmonic oscillator potentials, is ideal in this regard. The appropriate nuclear coordinate propagation along each electronic state “path” (prescription of all temporal hopping points) can be done analytically, for any number of oscillator modes. Thus, in the Spin-Boson problem, the difficulty is reduced to summing the amplitude contributions obtained from all allowed spin paths. This number grows rapidly with the time interval of interest (see below), but even simple direct enumeration enables multidimensional dynamics to be computed over useful time periods.7a,b,9,10 Recent innovations which circumvent the direct spin-path enumeration bottleneck have been proposed and successfully implemented.11-13 The difficulty faced in extending the Spin-Boson path integral methodology to general coupled anharmonic potential surfaces is the propagation of the appropriate wavepackets efficiently and accurately on either of the indiVidual potential energy surfaces. For reasons discussed above, exact numerical propagation cannot be performed with current computers. Hence, we propose here to use an approximate single-surface propagation technique which is rapid and multidimensional, namely, Gaussian Wavepacket Dynamics (GWD).16 In spite of the fact that GWD applied to motion in anharmonic potentials is not an exact algorithm by any means, it has been very successful in describing a number of short-time (subpicosecond) processes including molecular and atom-surface scattering,16,17 excited-state photophysics,18 and even evaluation of the Golden Rule formula for nonadiabatic transition rate constants.19 Provided that we do not abuse GWD by applying it in situations where it is obviously not suitable, we can use the path integral formulation of molecular curve crossing to extract curvecrossing dynamics from simple GWD wavepacket dynamics on individual potential surfaces. This paper is organized as follows. In Section 2, we review the path integral theory of nuclear coordinate wavepacket dynamics on coupled potential surfaces, with emphasis on approximate implementation using GWD to perform the required single-surface wavepacket propagations. In Section 3, we present a simple example involving photodissociation in a 1-D model as a practical illustration of the method. In Section 4, strengths and weaknesses of the method are discussed, as well as potential applications to more complex systems. 2. Formulation of Curve Crossing as a Wavepacket Hopping Process Consider the standard curve-crossing Hamiltonian:
H ) |1〉〈1|hˆ 1 + |2〉〈2|hˆ 2 + gˆ (|1〉〈2| + |2〉〈1|)
(1)
Here |1,2〉 are (structureless) electronic states, hˆ 1,2 are nuclear
coordinate Hamiltonians corresponding to single potential energy surfaces, and gˆ is the nonradiative or “nonadiabatic” coupling operator.20 The nuclear coordinate Hamiltonians have the form hˆ j ) Tˆ + Vˆ j, where Tˆ is the kinetic energy operator and Vˆ j the potential surface corresponding to electronic state j ) 1, 2 (which is a known function of the nuclear coordinate configuration). The dimensionality of the nuclear coordinate space and the choice of coordinate system (Cartesian, curvilinear, ...) are arbitrary. The nuclear coordinate operator gˆ can also be an arbitrary function of nuclear coordinates and momenta. The corresponding wavepacket states associated with the Hamiltonian in eq 1 are
|Ψ(t)〉 ) φ1(x,t)|1〉 + φ2(x,t)|2〉
(2)
Here φ1,2(x,t) are the nuclear coordinate wavepackets associated with electronic states 1 and 2. We note that they depend on nuclear coordinate configuration (which will be multidimensional in general) and time. The time-dependent Schro¨dinger equation which governs the dynamics of two coupled wavepackets is (with p ) 1)
∂ ˆ |Ψ(t)〉 i |Ψ(t)〉 ) H ∂t
(3)
For clarity, it is worthwhile to express the content of eqs 1-3 in an explicit matrix/vector form:
[ ] [
i
hˆ 1 φ˙ 1(x,t) ) φ˙ 2(x,t) g(xˆ )
][ ]
g(xˆ ) φ1(x,t) hˆ 2 φ2(x,t)
(4)
The elements of the Hamiltonian matrix on the right-hand side of eq 4 are all operators on nuclear coordinates. We have notated the nonadiabatic coupling as a function of the nuclear coordinate configuration. This is the case we will consider in detail in this paper. In this form, it is clear from eq 4 how the evolution of the two nuclear wavepacket components is prescribed. In fact, g could be a function of the conjugate nuclear momenta as well. A path integral prescription for solving eq 3 or eq 4 begins with the formal solution |Ψ(t)〉 ) exp(-iH ˆ t)|Ψ(0)〉. One then “shreds” the propagator into N short-time propagators21,22 corresponding to time interval (such that N ) t); i.e., exp(-iH ˆ t) ) [exp(-iH ˆ )]N. Progress is made by inserting complete sets of electronic states between each short-time propagator.8-10 We can then take advantage of explicit expressions for the matrix elements of the short-time propagator with respect to the electronic basis states. That is,
|Ψ(t)〉 )
∑
eN-1)1,2
...
∑ exp(-iH)|eN-1〉〈eN-1|exp(-iH)|eN-2〉...
e1)1,2
〈e2|exp(-iH)|e1〉〈e1|exp(-iH)|Ψ(0)〉 (5) Each factor 〈ej|exp(-iH)|ek〉 is an operator on nuclear coordinates which is specified below. Given explicit forms for these operators, each term in the electronic-state sum propagates the initial state |Ψ(0)〉 into some final two-component wavepacket state. To get the solution of the time-dependent SE (3) at output time t, we have to sum the contributions from all intermediate electronic-state configurations, as indicated in eq 5. [Note that carets on the nuclear coordinate operators have been suppressed in eq 5 for notational simplicity, a policy which will be followed for the remainder of the paper.]
7898 J. Phys. Chem., Vol. 100, No. 19, 1996
Coalson
The matrix elements 〈ej|exp(-iH)|ek〉 play a critical role in the procedure just outlined. Since we take f 0 to obtain convergence, the short-time propagator can be disentangled in any desired way.21,22 For our purposes, it proves convenient to use the decomposition
exp(-iH) ) exp(-ih0/2) exp(-ig(x)σx) exp(-ih0/2) (6) where h0 ≡ |1〉〈1|h1 + |2〉〈2|h2 is the electronically diagonal part of H (i.e., in the absence of nonadiabatic coupling) and σx is the Pauli matrix. The electronic matrix elements of the righthand side can easily be taken, yielding
〈ej|exp(-iH)|ek〉 ) exp(-ihej/2)Fej,ek(x) exp(-ihek/2) (7) where Fej,ek(x) is the function
Fej,ek(x) )
{
ej ) ek ej * ek
cos(g(x)); -i sin(g(x));
(8)
Putting together these short time propagators leads to an expression for the desired coupled surface wavepacket state |Ψ(t)〉. For simplicity, let us consider the case that at t ) 0 there is amplitude only on one electronic state |e0〉 (where e0 may be either 1 or 2). That is, |Ψ(0)〉 ) φe0(x)|e0〉.23 This will generate an output state |Ψ(t)〉 which has components in both electronic states. The output nuclear coordinate state on electronic state |eN〉, eN ) 1, 2, is then given by
φeN(x,t) ) ∑ exp(-iheN/2)FeN,eN-1(x) C
exp(-iheN-1)FeN-1,eN-2(x) exp(-iheN-2)...Fe1,e0(x) exp(-ihe0/2)φe0(x) (9) with ∑C ≡ ∑e1)1,2‚‚‚∑eN-1)1,2 representing the sum over all possible configurations (“C”) of the N - 1 intermediate electronic-state indices. Equation 9 gives a formally exact prescription for calculating the quantum dynamics between nuclear coordinate wavepackets on two coupled potential energy surfaces in terms of sums of amplitudes obtained from sequential single-surface propagations. There are two practical bottlenecks. First, it is impossible to do exact single-surface wavepacket propagations on general potential surfaces for systems with more than 3-4 degrees of freedom. Second, it is necessary to compute many singlesurface wavepacket trajectories to evaluate eq 9. The number of configurations of electronic intermediate states is equal to 2N-1. The first limitation can be overcome with no loss in accuracy in special cases where both potential surfaces are at most quadratic in all degrees of freedom, and the nonadiabatic coupling function is a constant. In this case, an initially Gaussian wavepacket remains Gaussian for all time whether it is propagated on V1 or V2. Each sequence of nuclear coordinate wavepacket propagations indicated in eq 9 then generates an output wavepacket which is Gaussian and whose parameters can be determined by evolving an appropriate finite parameter set using the GWD algorithm.16 Each output Gaussian is weighted by a complex number, independent of x (given by the product of Fek,ej appropriate to the selected intermediatestate configuration C). In fact, in the case where V1 and V2 both consist of separable harmonic oscillator potentials for each nuclear degree of freedom (“mode”) in the system and, furthermore, where these potentials vary only with respect to the equilibrium position of each mode
(and an overall vertical displacement in energy), the propagation of the multimode wavepacket factorizes into one-dimensional and analytical pieces. This is the Spin-Boson Hamiltonian, and the ease of nuclear coordinate wavepacket propagation for this problem is an important feature which has enabled computation of real-time many-body quantum dynamics using the path integral formalism.8-13 Beyond the Spin-Boson model, but still within the classification of “harmonic potential surface, constant g” are general harmonic potentials, which may involve frequency shifts and Duschinsky rotations of the normal modes of one surface relative to the other. Although the relevant nuclear coordinate wavepacket dynamics cannot be done analytically in this case, numerical solution by integrating a finite number of first-order coupled ordinary differential equations (that specify the parameters in the multidimensional Gaussian) can be performed according to the GWD equations of motion. Thus, the method outlined above provides a route to numerically exact solution of curve-crossing dynamics involving arbitrary quadratic surfaces coupled by a constant nonadiabatic coupling matrix element. If the potential surfaces are multidimensional, coupled, and anharmonic or if g depends on the (multidimensional) coordinate x, it is not obvious how to achieve exact solution by the “sequential wavepacket propagation” strategy outlined above.24 This leads us to suggest the use of the GWD approximation for the nuclear coordinate wavepacket propagations. The advantage is that nuclear coordinate wavepackets for a large number of coupled degrees of freedom can be efficiently evolved with this method. Also, since the wavepackets remain relatively welllocalized, they make possible a further approximation that enables treatment of problems where g depends on x. Specifically, at each “hopping point” [t ) (j + 1/2), j ) 0, N - 1], we are required to multiply the wavepacket obtained upon propagation up to this time by either cos(g(x)) or -i sin(g(x)) depending on whether the electronic state remains the same or changes in the next interval. This operation will distort a Gaussian wavepacket into a non-Gaussian one, which imperils the proposal to use GWD to propagate the wavepacket all the way through any sequence of electronic-state hops. To prevent this complication, we suggest the replacement of cos(g(x)) f cos(g(xt)), where xt is a parameter in the Gaussian wavepacket that specifies its center in coordinate space at the time where multiplication by cos(g(x)) is called for in eq 8 above (and analogously if the intermediate electronic-state configuration selects -i sin(g(x))). Thus, the Gaussian remains Gaussian as it passes through each hopping point, and we simply collect an overall multiplying factor that reflects the trajectory of the nuclear coordinate wavepacket in the way just indicated. Again, the success of this last approximation requires that the propagating Gaussian remain spatially well localized for all times. This condition is generally satisfied to a reasonable approximation when the GWD approximation is valid. (If the wavepacket spreads too much, any quadratic approximation to the anharmonic potential energy function will fail.) Let us summarize the proposed Gaussian Wavepacket Dynamics-path integral (GWD-PI) procedure for computing the nuclear coordinate wavepacket φeN(x,t) which emerges at time t in electronic state eN, given initial preparation in electronic state e0 of a Gaussian nuclear coordinate wavepacket φe0(x): (i) Select an intermediate state configuration (e1, ..., eN-1). (ii) Propagate on potential surface Ve0 for time /2 using GWD. This generates an output Gaussian wavepacket from the initial wavepacket; call it φ(x,/2). The coordinate space
Wavepacket-Path Integral Method for Dynamics
Figure 1. Intermediate electronic-state configuration (e1, e2, e3) ) (1, 1, 2) for the case where e0 ) eN ) 2 and N ) 4. An equivalent linear spin chain is also indicated.
center of the output wavepacket is given directly by the GWD algorithm; call it xt. (iii) Multiply φ(x,/2) by a (coordinate independent) factor Fe1,e0(xt). Thus, the output wavepacket at this stage is still Gaussian. (iv) Propagate φ(x,/2) forward by a time on potential surface Ve1 using GWD. (The multiplicative factor introduced in step iii is simply “carried along”.) The output wavepacket, φ(x,3/2), remains Gaussian. (v) Repeat steps iii and iv until the desired final time is reached. [After the last hopping point, namely (N - 1/2), propagation to N ) t is performed on VeN.] (vi) Repeat steps ii-iv for all 2N-1 intermediate electronicstate configurations. Sum (coherently) all output Gaussians obtained by this procedure. This gives φeN(x,t). Note that the procedure for calculating the output wavepacket associated with a given intermediate-state configuration can be summarized in a pictorial way by the kind of diagram or “roadmap” in Figure 1. The sequence of propagations on V1 and V2 is clearly indicated. At each potential hopping time (demarcated via dashed vertical lines), the wavepacket is multiplied by one of two factors, depending on whether it hops or does not hop at this juncture. (In the prototypical twoelectronic-state problem under direct study here, the electronicstate configuration can be thought of as a linear spin chain, with spin up, down signifying propagation on V1, V2, respectively. This correspondence is indicated in Figure 1. Since we do not utilize the analogy in any essential way, we will not pursue it further in this paper.) A number of variations on the theme just elaborated are possible. For example, the approximation cos(g(x)) f cos(g(xt)) could be replaced by cos(g(x)) f 〈cos(g(x))〉, where the average is over the probability density associated with the Gaussian wavepacket just before the appropriate hopping time (and analogously if the “sin” alternative is selected). The latter option may be more accurate but will require more computational labor. Furthermore, the case where g(x) is explicitly time dependent, as it is in related problems of excited-state optical spectroscopy, can easily be treated by the proposed method. Assuming the time slice is short enough that g(x,t) is constant over this interval, we can simply replace cos(g(x,t)) f cos(g(xthop,thop)), where thop is the hopping time in question, and xthop is the center of the Gaussian wavepacket at this time (and likewise if the “sin” alternative is selected.)25 It is also of interest to note briefly the connection between the electronic-state path integral strategy considered here and a perturbative strategy termed the wavepacket perturbation theory (WPT) introduced in earlier work.26-28 The two techniques set out by different routes to achieve the same goal, namely, to sum wavepacket amplitudes over “all possible” sequences of intermediate electronic-state configurations. In WPT, all n-hop paths are counted at the nth order of perturbation theory (which
J. Phys. Chem., Vol. 100, No. 19, 1996 7899 involves an n-dimensional quadrature over hopping times). In the path integral method, at a given level of temporal discretization (N), many paths of many orders are considered (but not all paths at any order). In practice, it proves difficult to evaluate the nth-order time quadratures associated with WPT beyond about fourth order, while convergence of the discretized electronic-state path integral can often be achieved using a fairly coarse time slice. Therein lies part of the utility of the nonperturbative path integral strategy. Finally, we return to the issue of the growth of the number of intermediate electronic-state configurations with the number of time slices. Currently, one can perform explicit path summations with 20-30 time slices (requiring enumeration of 220-230 intermediate electronic-state configurations). This certainly restricts the interval of time over which these calculations can be performed. Fortunately, for many processes, particularly in the condensed phase, temporal signatures of experimental interest decay rapidly (on a subpicosecond time scale). Direct enumeration over electronic-state paths has proven useful in previous condensed-phase studies of electronicstate population evolution and excited-state spectroscopy within the Spin-Boson model.7,9,10 We do not expect the additional effort of propagating Gaussian wavepackets by integration of the necessary first-order differential equations to further encumber the computational effort to a significant degree. 3. Numerical Illustration We consider a scenario in which a molecule, initially in the vibrational ground state φG(x) associated with ground electronic state |G〉, is excited by a weak monochromatic electric field to a pair of nonadiabatically coupled electronic excited states. For concreteness, we focus on the common experimental situation where one member of this excited-state pair, 1, is radiatively coupled to G, i.e., “bright”, while the other excited state, 2, is radiatively “dark”. In this scenario, σ(ωL), the absorption cross section at laser frequency ωL, may be calculated as
σ(ωL) )
Re ∞ ∫ dt ei(ωL+EG)t〈φG(x)|φ1(x,t)〉 π 0
(10)
where φ1(x,t) is the component of the initial nuclear coordinate wavepacket φG(x) which evolves on the excited-state potential surface V1(x), 〈φG(x)|φ1(x,t)〉 signifies the inner product of the indicated nuclear coordinate states (computed by integration over all coordinate space), and EG is the energy eigenvalue associated with vibrational eigenstate φG(x) on the ground-state potential surface.26,29 The key ingredient is φ1(x,t), which includes effects due to nonadiabatic coupling between excited states 1 and 2. Evaluation of φ1(x,t) requires solution of the coupled surface Schro¨dinger equation (4). Given such a solution, a variety of interesting experimentally observable signatures can be extracted. In addition to the total absorption cross section, they include resonance Raman excitation profiles (RREP’s) and partial cross sections for photodissociation. To compute RREP’s, we would simply overlap φ1(x,t) with the appropriate final vibrational state of the ground potential surface to which the molecule returns after the Raman scattering event. HalfFourier transformation as indicated in eq 10 (without taking the real part) gives a complex amplitude whose square is essentially the measured scattering intensity.26 Similarly, photodissociation cross sections are obtained by propagating the evolution of both components φ1,2(x,t) until these wavepackets are “asymptotic” (i.e., correspond to a configuration where the molecule is dissociated).27 Projection of these wavepackets onto appropriate final states of the dissociated complex yields the
7900 J. Phys. Chem., Vol. 100, No. 19, 1996
Coalson
cross section to wind up in electronic state 1 or 2 at any laser frequency (see below). We should also note that the evolution of φ1,2(x,t) is the central quantity involved in another type of experiment, namely, excitation of the system via an ultrashort laser pulse.31,32 A sufficiently short pulse creates at time t ) 0 a copy of φG(x) on potential surface V1(x). The subsequent fate of this initial condition is determined directly from the solution of eq 4. Interesting experimentally addressable questions include the probability that the system dissociates in electronic state 1 or 2 and the internal vibrational states accessed by the dissociated species. As a specific illustration of the use of the proposed method to study these dynamical properties, we consider a onedimensional system of mass m, which is bound in a harmonic oscillator potential in its electronic ground state. It is then photoexcited into a radiatively coupled electronic state 1, which is nonradiatively coupled to another excited electronic state, 2. The potential functions associated with electronic states 1 and 2 are chosen to represent direct photodissociation in both states. Specifically,
V1(x) ) A1e-R1x; V2(x) ) A2e-a2x + V0
(11)
The nonadiabatic coupling function g(x) is also chosen to be an exponentially decaying function of position; i.e., g(x) ) g0 exp(-agx). The numerical values in the calculations presented below are as follows. The vibrational quantum of VG is 670 cm-1. The parameters defining the excited-state potentials in eq 11 and the nonadiabatic coupling function are
A1 ) 0.41 eV; g0 ) 0.03 eV;
A2 ) 0.82 eV; -1
a1 ) 3 Å ;
Figure 2. Potential energy functions and initial probability density. The nonadiabatic coupling function, which decays exponentially with increasing nuclear coordinate, is not shown.
V0 ) -0.205 eV; a2 ) 3 Å-1;
ag ) 1 Å-1 (12)
The initial equilibrium point of the ground-state surface is taken as 0.1 Å in the coordinate system used to specify the excitedstate potentials in eq 11, which is 0.13 Å “inside” the V1/V2 crossing point (cf. Figure 2). The particle mass is 10 amu. Thus, at t ) 0, a Gaussian wavepacket corresponding to the ground vibrational eigenfunction of VG was placed on V1. It was subsequently evolved on the coupled excited surfaces V1,2 according to eq 4 using two methods. The first was the GWDPI method proposed in this paper. A variationally optimized time-dependent quadratic fit to the exponential potential functions V1 and V2 was utilized to guide the Gaussian wavepacket evolution. (This effective harmonic potential has previously been shown27,30 to be somewhat more accurate for this type of potential than the thawed Gaussian approximation,16 which entails a local Taylor series expansion around the instantaneous center of the wavepacket.) The appropriate sum over intermediate spin configurations was performed by direct enumeration. Since the dissociation is direct and rapid, the wavepackets pass through the crossing point and into the asymptotic regime (i.e., where all exponential terms in V1,2 and g have gone to zero) rapidly. Converged wavepackets φ1,2(x,t) in the asymptotic regime were obtained with N ) 4-6. The second method employed to propagate the coupled surface wavepackets φ1,2 was a standard two-surface numerical grid integrator based on the split-operator algorithm of Feit and Fleck.33 This method gives numerically exact results for the problem of interest here and hence is ideal for calibrating the accuracy of the GWD-PI method. Of course, the split-operator integrator is not a many-body technique. It becomes impractical for systems consisting of more than 3-4 degrees of freedom.
Figure 3. Re[φ1(x)] and |φ1(x)|2. GWD-PI approximation (straight line) is compared to the exact result (dashed line).
The ultimate goal of developing the GWD-PI method is to employ it on multidimensional (e.g., condensed phase) systems, where extant methods are not completely satisfactory. An essential test of the accuracy of the GWD-PI method for this problem is to examine the asymptotic wavefunctions on both surfaces. If these are accurate, then any observable property extracted from them, such as the quantum yield (total probability) to exit on either surface and partial cross sections for photodissociation as a function of laser frequency, will be accurate. Results are shown in Figure 3 for the real part and the squared magnitude of φ1(x,t) at time t ) 3.6 × 10-14 s, when the nuclear coordinate wavepackets emerging on both surfaces are essentially asymptotic; i.e., dissociation is complete. Analogous results for φ2(x,t) are shown in Figure 4. In all panels, the GWD-PI approximation (solid line) is compared to the exact numerical grid solution (dashed line). These results are encouraging. A high degree of accuracy in the underlying wavefunctions on both surfaces is evident. This guarantees the accuracy of any observable property
Wavepacket-Path Integral Method for Dynamics
J. Phys. Chem., Vol. 100, No. 19, 1996 7901 probability densities indicated in Figures 3 and 4 is that they are non-Gaussian (particularly the bright surface component |φ1(x,t)|2). The GWD-PI method has no difficulty reproducing this feature since it represents the output wavepackets φ1,2(x,t) as coherent superpositions of many complex-valued Gaussians. Complicated waveforms can be synthesized in this way. In particular, distortion of the wavepacket from a simple Gaussian as it goes through a crossing of the two potential surfaces and other wavemechanical interference effects associated with the curve-crossing process are naturally accounted for by the GWDPI method.
Figure 4. Re[φ2(x)] and |φ2(x)|2. GWD-PI approximation (straight line) is compared to the exact result (dashed line).
Figure 5. Partial photodissociation cross sections σ1,2 vs photon excitation energy. GWD-PI (straight line) is compared to the exact (dashed line) results. The excitation energies assume a value of EG ) 0.
extracted from them. To demonstrate this, we show in Figure 5 partial cross section for photodissociation into both electronic channels as a function of laser frequency. These are obtained by projecting the asymptotic outgoing wavepackets onto outgoing energy normalized plane waves corresponding to a wavevector whose dependence on ωL is determined by conservation of energy. Specifically, the cross section to absorb a photon of frequency ωL and dissociate on surface V2 is given by27
σ2(ωL) )
∞ m |∫-∞ dx e-ik2(ωL)xφ2(x,tf∞)|2 (13) 2πk2(ωL)
with k2(ωL) determined by the energy conservation constraint pωL + EG ) k22(ωL)/2m + V0 and analogously for σ1(ωL). (The spectra in Figure 5 correspond to a value of EG ) 0. Changing this value simply shifts all spectra by a constant value, namely EG.) As expected, the accuracy of the GWD-PI results is respectable. Examination of the probability densities indicated in the righthand panels of Figures 3 and 4 reveals two interesting features. First, the integrated area under |φ2(x,t)|2 is 0.60 (vs, of course, 0.40 for the area under |φ1(x,t)|2). Thus, there is a 60% probability that a system prepared in the wavepacket state φg(x) on electronic state 1 will dissociate into state 2.34 Clearly, this is a result of strong nonadiabatic coupling between the two states |e1,2〉. The GWD-PI technique is nonperturbatiVe in the nonadiabatic coupling strength, which makes it complementary to weak-coupling approximations (like the golden rule approach to nonadiabatic transition rates and wavepacket perturbation theory). The latter are easier to implement but break down in strong coupling situations. A second interesting aspect of the
4. Discussion and Conclusion In this paper, we have presented an approximate method for calculating multidimensional curve-crossing dynamics. The method consists of propagating wavepackets sequentially on the individual potential surfaces V1,2. A practical prescription for determining the hopping times (times at which transfer to the other surface is considered) can be specified. Coupling between the surfaces is obtained by coherently summing the set of output wavepackets obtained from all possible hopping paths through the electronic-state space. To treat systems consisting of coupled multidimensional nuclear coordinates, we have proposed using the well-known gaussian wavepacket dynamics (GWD) algorithm. Let us summarize again the weaknesses and strengths of this Gaussian wavepacket dynamics-path integral (GWD-PI) method. The major weakness of the method is the possibility that the GWD approximation (which forces the evolving nuclear coordinate wavepacket to remain Gaussian for all times) will break down. On sufficiently anharmonic potential surfaces and over sufficiently long time intervals, such a breakdown is inevitable. However, for more mildly anharmonic systems and in particular for condensed-phase systems where short time (subpicosecond) dynamics and correspondingly low spectral resolution are the rule rather than the exception, this method may prove useful, since the propagation of multidimensional Gaussians is not much harder than the integration of Newton’s equations of classical mechanics. Furthermore, one can hope for the continued development of improved semiclassical techniques, which may eventually allow their substitution for the GWD ingredient in the present scheme.35 A related weakness is the additional approximation invoked in order to treat position dependence in the nonadiabatic coupling function. As discussed in the text, the narrower our Gaussian wavepackets, the better the approximation. We expect that if GWD is satisfactory for performing the required nuclear coordinate wavepacket evolutions, the further approximation of replacing the full coordinate dependence of the terms which contain the nonadiabatic coupling by their value at the center of the Gaussian wavepacket that they directly multiply will be reasonable. A final problem is the geometric growth in intermediate electronic intermediate states with the number of time slices, N. Fortunately, this growth is independent of the spatial dimensionality of the nuclear coordinate space, so it does not represent a problem which becomes unmanageable as the spatial dimensionality grows. For a system consisting of M coupled excited states, the number of intermediate electronic configurations goes as MN. For the prototypical case where M ) 2, we can achieve N ) 20-30 on currently available workstations. Obviously, the maximum value of N accessible will decrease with increasing M. Perhaps tricks inspired by statistical mechanics, which have proven very successful in alleviating this problem in the standard Spin-Boson model,11-13 can be developed for the more general systems considered herein.
7902 J. Phys. Chem., Vol. 100, No. 19, 1996 Now let us summarize the strengths of the method. (i) In the special case that the nonadiabatic coupling strength is constant and both potential surfaces V1,2 are harmonic, the GWD-PI method gives exact results upon convergence of the temporal discretization procedure. Application of the method in cases where the quadratic potential functions are more elaborate than the commonly invoked “linearly displaced harmonic oscillator” form would represent some advance over the capabilities of the standard Spin-Boson model. In the general case, where the results are approximate, we can still point to the following potential advantages: (ii) The method can be used to study systems of high dimensionality (100’s or even 1000’s of degrees of freedom are not out of the question). (iii) It can treat nonadiabatic coupling functions with arbitrary functional dependence on the nuclear coordinate configuration and arbitrary externally imposed time dependence, too. (iv) It can be used in principle to treat systems with any number of nonadiabatically coupled surfaces. (Of course, in practice, if direct enumeration of intermediate electronic states is employed, then the larger the number of coupled surfaces, the smaller the time interval that can be studied.) (v) The method is not restricted to weak nonadiabatic coupling. Any coupling strength can be studied. (Of course, the higher the coupling, the smaller the value of necessary to obtain convergence. However, we demonstrated in the previous section that substantial nonadiabatic coupling effects can be treated in practice without difficulty.) (vi) The method is completely wavemechanical. Effects due to superposition of complex amplitudes, which can play an important role in curve-crossing problems, are naturally included. (vii) In principle, the method can be used to treat nonadiabatic coupling functions that depend on nuclear coordinate momenta, pˆ . This situation arises when effects due to breakdown of the Born-Oppenheimer approximation occur, and the problem is formulated in the “adiabatic” representation.29 Following the strategy employed here to treat coordinate dependence in the nonadiabatic coupling function, momentum dependence could be included in an approximate way by replacing the operator pˆ by the average momentum of the Gaussian wavepacket which it operates directly on. For all the reasons just cited, the GWD-PI method appears to be worth pursuing further, particularly in condensed-phase applications, where it could be benchmarked for small nonadiabatic coupling strengths by golden rule calculations and then, if it proves successful, utilized to further explore the strong nonadiabiatic coupling regime. Acknowledgment. This work was supported by NSF Grant CHE-9101432. I thank M. Kurnikova for valuable assistance with the figures and an anonymous referee for helpful suggestions. Finally, I am grateful to Jim Kinsey for introducing me to the fascinating subject of nonadiabatic transitions in molecules. References and Notes (1) See, for example: Scherer, N. F.; Zeigler, L. D.; Fleming, G. R. J. Chem. Phys. 1992, 96, 5544. Black, J. F.; Waldeck, J. R.; Zare, R. N. J. Chem. Phys. 1990, 92, 3519. Trulson, M. O.; Mathies, R. A. J. Phys. Chem. 1990, 94, 5741. (2) See, for example: Sumi, H.; Marcus, R. A. J. Chem. Phys. 1986, 84, 4894. Ulstrup, J. Charge Transfer Processes in Condensed Media; Lecture Notes in Chemistry 10; Springer-Verlag: Berlin, 1979. (3) Kosloff, R. J. Phys. Chem. 1988, 92, 2087. (4) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562. Tully, J. C. J. Chem. Phys. 1990, 93, 1061. Krylov, A. I.; Gerber, R. B.; Apkarian, V. A. J. Phys. Chem. 1994, 189, 261.
Coalson (5) Makri, N.; Miller, W. H. J. Chem. Phys. 1987, 87, 5781. Kotler, Z.; Nitzan, A.; Kosloff, R. Chem. Phys. Lett. 1988, 153, 483. Waldeck, J. R.; Campos-Martinez, J.; Coalson, R. D. J. Chem. Phys. 1991, 94, 2773. (6) Nitzan, A.; Jortner, J. Chem. Phys. Lett. 1972, 15, 350; J. Chem. Phys. 1973, 58, 2412. Todd, M. D.; Nitzan, A.; Ratner, M. A. J. Phys. Chem. 1993, 97, 29. (7) (a) Coalson, R. D.; Evans, D. G.; Nitzan, A. J. Chem. Phys. 1994, 101, 486. (b) Evans, D. G.; Coalson, R. D. J. Chem. Phys. 1995, 102, 5658. (c) Evans, D. G.; Coalson, R. D. J. Chem. Phys. 1996, 104, 3598. (8) Chakravarty, S.; Leggett, A. Phys. ReV. Lett. 1984, 52, 5. Leggett, A.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, W. ReV. Mod. Phys. 1987, 59, 1. (9) (a) Coalson, R. D. J. Chem. Phys. 1987, 86, 995. (b) Coalson, R. D. Phys. ReV. 1989, B39, 12502. (c) Coalson, R. D. J. Chem. Phys. 1991, 94, 1108. (10) Evans, D. G.; Coalson, R. D. J. Chem. Phys. 1993, 99, 6264; ibid. 1994, 100, 5605. (11) Makarov, D. E.; Makri, N. Chem. Phys. Lett. 1994, 221, 482. Makri, M.; Makarov, D. E. J. Chem. Phys. 1995, 102, 4600, 4611. (12) Egger, R.; Mak, C.; Weiss, U. J. Chem. Phys. 1993, 100, 2651. Mak, C. H.; Weiss, U. Phys. ReV. E 1994, 50, 655. (13) Winterstetter, M.; Domcke, W. Phys. ReV. A 1993, 47, 2838; 1993, 48, 4272. Krempl, S.; Winterstetter, M.; Plo¨hn, H.; Domcke, W. J. Chem. Phys. 1994, 100, 926. (14) Campos Martinez, J.; Waldeck, J. R.; Coalson, R. D. J. Chem. Phys. 1992, 96, 3613. (15) Jungwirth, P.; Gerber, R. B. J. Chem. Phys., in press. (16) Heller, E. J. J. Chem. Phys. 1975, 62, 1544. Heller, E. J. Acc. Chem. Res. 1981, 14, 368. (17) Drolshagen, G.; Heller, E. J. J. Chem. Phys. 1983, 79, 4005. Drolshagen, G.; Heller, E. J. Surf. Sci. 1984, 139, 260. (18) Kulander, K.; Heller, E. J. J. Chem. Phys. 1978, 69, 2439. (19) Neria, E.; Nitzan, A. J. Chem. Phys. 1993, 99, 1109; Chem. Phys. 1994, 183, 351. (20) The method developed in this work can be applied to a number of curve-crossing processes, ranging from intersystem crossing to electron transfer. We shall therefore use the terms “nonadiabatic” and “nonradiative” (somewhat loosely) interchangeably to refer to any process governed by a Hamiltonian of the form indicated in eq 1. (21) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (22) Schulman, L. S. Techniques and Applications of Path Integration; Wiley and Sons: New York, 1981. (23) Because the time-dependent SE is linear, no loss of generality is entailed by this “specialization”. If the initial state has nuclear wavepacket components on both surfaces, then we can do the temporal propagation for each component separately and coherently add the results of the two calculations. (24) An important special case of this situation is when both potential energy surfaces separate into a sum of one-dimensional contributions (in the same coordinate system). Then provided that g factorizes in the same coordinate system and that the initial wavepacket state also factorizes in these coordinates, all the evolved nuclear wavepackets remain factorized for all times and can be calculated by independent one-dimensional propagations. This simplifies the computation immensely. This feature has been nicely exploited by: Ilk, G.; Makri, N. J. Chem. Phys. 1994, 101, 6708. However, in the general case, these conditions are not satisfied and a different approach appears to be required. (25) These possibilities were pointed out by an anonymous referee. (26) Coalson, R. D.; Kinsey, J. L. J. Chem. Phys. 1986, 85, 4322. (27) Coalson, R. D. J. Chem. Phys. 1987, 86, 6823. Coalson, R. D. AdV. Chem. Phys. 1989, 73, 605. (28) Rama Krishna, M. V. J. Chem. Phys. 1990, 93, 3258. (29) Ko¨ppel, H.; Domcke, W.; Cederbaum, L. S. AdV. Chem. Phys. 1984, 57, 59. (30) Coalson, R. D.; Karplus, M. J. Chem. Phys. 1990, 93, 3919. (31) Grubele, M.; Roberts, G.; Dantus, M.; Bowman, R. M.; Zewail, A. H. Chem. Phys. Lett. 1990, 166, 459. (32) Pollard, W. T.; Lee, S.-Y.; Mathies, R. A. J. Chem. Phys. 1990, 92, 4012. Metiu, H.; Engel, V. J. Chem. Phys. 1990, 93, 5693. Kohler, B.; Yakovlev, V. V.; Che, J.; Krause, J. L.; Messina, M.; Wilson, K. R.; Schwentner, N.; Whitnell, R. M.; Yan, Y. J. Phys. ReV. Lett. 1995, 74, 3360. (33) Feit, M. D.; Fleck, J. A., Jr.; Steiger, A. J. Comput. Phys. 1982, 47, 412. (34) The integrated probability densities for electronic states 1 and 2 obtained by the GWD-PI method are each high by about 0.005 (roughly 1%). This method does not conserve probability automatically, so the deviation of the total probability to be in both electronic states from unity is a meaningful indicator of the accuracy of the approximation. (35) Heller, E. J. J. Chem. Phys. 1991, 94, 2723. Sepu´lveda, M. A.; Heller, E. J. J. Chem. Phys. 1994, 101, 8004; 1994, 101, 8016.
JP9531826
