Acc. Chem. Res. 2006, 39, 127-134
Guided Gaussian Wave Packets E. J. HELLER Departments of Physics and Chemistry, Harvard University, Cambridge, Massachusetts 02138 Received May 5, 2005 ABSTRACT We discuss several semiclassical Gaussian wave packet approaches with emphasis on one that is not very well-known, that is, the offcenter guiding approach. Off-center guiding of (thawed) Gaussian wave packets uses the same information as other Gaussian propagation schemes, that is, the full van Vleck determinant and multiple trajectories. It retains the well-known caustic smoothing property of Gaussians. The off-center guiding of Gaussians can handle hard chaos and other highly nonlinear situations, where the van Vleck propagator, for example, has over 30 000 separate branches.
1. Introduction Quantum dynamics is undergoing a rapid development, which parallels earlier progress in electronic structure theory. The latter can be thought of as the necessary “static” precursor to the nuclear dynamics problem occurring on and between Born-Oppenheimer potential energy surfaces: this is the essence of chemistry and chemical reactivity. (“On the fly” methods, which do the quantum dynamics in conjunction with electronic structure, compute structure information only where and when it is needed; these are very powerful and are rapidly gaining favor.1) The issue of how to deal with the potentially huge basis sets required to follow a quantum time developing wave function of say 6 or 60 or 600 degrees of freedom is the subject of many approximate methods, but one of the leading ideas is to use the convenient and physically motivated Gaussian wave packets as a basis for description. This brings mathematical convenience as well as semiclassical motivation and approximation into the mix, since Gaussian wave packets are venerable workhorses of heavy particle quantum dynamics, dating back to Schro¨dinger’s use of them for the harmonic oscillator, where he noticed they behaved classically. Semiclassical Gaussian wave packet methods are now seeing widespread use in a wide variety of applications to chemical systems. The range of problems that they have addressed is very large and growing. Hybrid methods that are partly semiclassical and partly ab initio are also receiving much attention. These methods go beyond merely solving the Schro¨dinger equation for a physical system by providing intuition and mechanism behind the results. It is impossible to do justice to all the
Eric Heller received his B.A. at the University of Minnesota in 1968 and a Ph.D. in Chemical Physics at Harvard in 1973. He was then a postdoctoral fellow at the University of Chicago. A major research focus since 1976 has been semiclassical and wave packet approaches to chemical dynamics. 10.1021/ar040196y CCC: $33.50 Published on Web 12/31/2005
2006 American Chemical Society
FIGURE 1. Anatomy of a Gaussian. A Gaussian wave packet is displayed in coordinate space with the role of the parameters at, bt, qt, pt, and φt shown. Within the TGA, the parameters are all controlled by classical trajectories with qt,pt following the classical trajectory leading from q0,p0 at time t ) 0. The phase φt is the classical action. significant ideas and papers without writing a large review article. This Account has the more modest goal of reviewing purely semiclassical guided wave packet methods and discussing one approach that has not received the attention it might deserve in chemical applications. This is the method of “off-center guiding” of wave packets,2-5 which may be thought of as a significant improvement on the older “thawed Gaussian approximation” or TGA.6 It is more informal and much simpler to implement than the full stationary phase approach to the coherent state semiclassical representation (CSSR).7,8 A Gaussian wave packet is associated with a classical point in phase space through its expectation value with the position and momentum operators. It is useful to define a few things about Gaussians. A general Gaussian wave packet in one dimension can be written as
ψt(q) ≡
( ) [ ( ) [
)
Re Rt πp 2at πp
1/4
exp -
1/4
exp -
]
Rt i i (q - qt)2 + pt(q - qt) + φt 2p p p
]
Rt ipt(q - qt) i (q - qt)2 + + φt 2p p p (1)
where Rt/2 ≡ (at + ibt) is a complex number with positive real part at and imaginary part bt of either sign. The Gaussian has average position qt ) ∫ψ/t (q)qψt(q) dq ) ∫|ψt(q)|2q dq and average momentum pt ) ∫ψ/t (q)(-ip(∂/∂q)ψt(q) dq. A little “anatomy lesson” is worthwhile. In Figure 1, we show a picture of a normalized general Gaussian in coordinate space and the aspects that each parameter in the Gaussian controls. The parameters at, bt, qt, pt, and φt are all real. Since the Gaussian wave packet has uncertainty in position and momentum, the classical analogy has to be extended to a distribution of trajectories with collectively similar uncertainties in position and momenVOL. 39, NO. 2, 2006 / ACCOUNTS OF CHEMICAL RESEARCH
127
Guided Gaussian Wave Packets Heller
FIGURE 2. The TGA is derived by locally expanding the potential quadratically around the moving center of the Gaussian. tum. The best way of doing this is through the Wigner transform,
FW(q,p) )
(πp1 )∫
∞
-∞
e-2ips/p ψ*(q - s)ψ(q + s) ds (2)
The Wigner function for the Gaussian eq 1 is obtained by carrying out the integral eq 2, which results in
FW 1 (q,p) )
(πp2 ) exp[-(2(a
2
1
+ b12)/(a1p))(q - q1)2 -
(1/(2a1p))(p - p1)2 + (2b1/(a1p))(q - q1)(p - p1)] (3) Phase space representations and pictures are the best way to understand the issues facing wave packet propagation. We will exploit the intuition suggested by them below.
2. Thawed Gaussians The original TGA idea is based on the notion that, initially, the expectation values of the position and momentum of a Gaussian obey the classical equations of motion (Ehrenfest’s theorem). This correspondence is exact for all times for a quadratic potential but still true at short times for an anharmonic potential. The TGA idea is to “cradle” the wave packet in a local harmonic expansion surrounding the center of the wave packet; this notion is clear from Figure 2. This leads to an effective time-dependent quadratic Hamiltonian as shown in Figure 2. If a Gaussian is propagated with a different initial q0,p0, then a different time-dependent effective Hamiltonian results, one that is optimized for the new Gaussian. Many more details about Gaussian wave packet propagation may be found in refs 9-11. The difficulties with the TGA increase with time; initial Gaussian wave packets do not remain Gaussian in anharmonic potentials. The problem arises from a global linearization of the classical dynamics around the guiding center of the wave packet. It is useful to make phase space pictures to see what happens when the dynamics is nonlinear. We take an ellipse (here a circle) to represent the initial Gaussian wave packet; the parameters of the ellipse are taken from the uncertainties of the inital Gaussian. Think of the ellipse as filled with δ function classical trajectories, and propagate them all in time. A typical result is shown in Figure 3. The TGA uses the center of the initial wave packet as initial conditions for the guiding trajectory and linearizes around that. The TGA approximation to the wave packet dynamics is correct in the “core” of the wave packet, but fails badly in other regions of phase space. The overlap with the final Gaussian near the top of the picture would be quite wrong. Note that the Gaussian propagated under 128 ACCOUNTS OF CHEMICAL RESEARCH / VOL. 39, NO. 2, 2006
FIGURE 3. (A) An initially Gaussian wave packet, represented by the phase space circle on the left, is propagated under a nonlinear Hamiltonian with the resulting nonlinearly distorted phase space distribution shown on the right. Also shown is the overlap of the propagated distribution with a second Gaussian. (B) TGA approximation shown in gray approximates the situation in the “core” of the propagated wave packet but does a very poor job in representing the overlap of the true distorted state with the Gaussian depicted near the top of the picture. the linear dynamics of the TGA remains Gaussian; this means that it’s phase space plot must be an ellipse, as shown in gray in Figure 3B. We should not have expected too much from the TGA; it propagates an entire, complete wave function with one trajectory (and its immediate environs, which is the information in the stability matrix and van Vleck determinant; see the next section). Nonetheless, the TGA has done remarkable service, in particular in deriving analytical approximations to such things as absorption and Raman spectroscopy, where many quantities depend only on very short time dynamics.12
3. Beyond the TGA The beginning of a “proper” and less heuristic theory of semiclassical wave packet propagation is to ask: what is the stationary phase limit of the Feynman propagator? The answer is the van Vleck-Gutzwiller semiclassical propagator, GVVG(q,q′;t):
GVVG(q,q′;t) )
( ) 1
2πip
d/2
∑j
(
|Det
) (
∂2Sj(q,q′;t) ∂q∂q′
1/2
exp iSj(q,q′;t)/p -
)
iπνj
2 (4)
The sum over j is over all the trajectories connecting q to q′ in time t; d is the number of degrees of freedom. The prefactor determinant is in fact the square root of the classical probability of starting at q′ at time t ) 0 and ending at q at time t. This determinant involves the
Guided Gaussian Wave Packets Heller
stability matrix and thus the linearized dynamics in the vicinity of each trajectory connecting q with q′ in time t: we can write, in one dimension,
∂2Sj(q,q′;t) 1/2 ∂p | ) | |1/2 ∂q∂q′ ∂q′
|
(5)
The phase is determined by the classical action Sj(q,q′;t) and the Maslov index νj (see below). The classical action is the time integral of the Lagrangian, L
Sj(q,q′;t) )
∫0t dt′ L ) ∫0t dt′ {p(t′)‚q3 (t′) H(p(t′),q′(t′))} (6)
along the jth classical path. H is the classical Hamiltonian, which is presumed to be the classical limit of H ˆ . Equation 4 was originally written down by van Vleck in 1928 without the summation or index ν.13 We note in passing that one can derive the TGA under the linearizing of the dynamics in the following way: quadratically expand the action S(q, q′, t) around q′ ) q0 and q ) qt, where qt ) qt(q0,p0) with (q0,p0) and (qt,pt) real. We will not carry this out here, but the result is the TGA when this linearized propagator is applied to the initial Gaussain centered at (q0,p0). This has all the problems already mentioned and is clearly more primitive than keeping the full nonlinear VVG propagator; it is also much less work. The TGA is not the full semiclassical result. Suppose we do want to semiclassically propagate the initial wave function, ψ0(q′):
ψt(q) )
∫dq′ GVVG(q,q′,t)ψ0(q′)
(7)
using the VVG propagator. There are two alternatives, namely, do the integral above by stationary phase or numerically. The latter is a “uniform” approximation; the former is a “primitive” semiclassical result. We consider both approaches, starting with the uniform or numerical approach. To carry this out, we need to determine Sj(q,q′;t) and do the integration over q,q′ in eq 4 numerically. One way of doing this, called cellular dynamics (CD), was introduced in ref 14; later, Miller correctly pointed out15 that cellular dynamics had an initial value representation (IVR) as its modus operandi.16 What was new about cellular dynamics however was the direct calculation (by whatever means) of the VVG time-dependent propagator applied to smooth wave functions. For some reason, despite the prime importance of the VVG propagator, no one had previously evaluated it’s performance as a semiclassical tool, except to note that the propagator itself is badly singular.17 This is not to say that the result is badly singular when the VVG propagator is applied to smooth wave functions; indeed it is not. The CD-IVR approach makes the integral easier by regularizing the integrand and eliminating “root searches”. The end result was that even when the dynamics is strongly nonlinear, the VVG propagator acting on a smooth wave function may be remarkably accurate.14 Cellular dynamics was used in refs 18 and 19 to successfully treat systems with mixed chaotic and
FIGURE 4. The initial Gaussian wave packet in the Morse potential at the top corresponds to a vertical manifold of trajectories, since it is very narrow. This wave packet is propagated according to the VVG semiclassical cellular dynamics and exact fast Fourier transform. There is almost no difference, even after the classical manifold (initially a vertical line, later escaping to the right and twisting) begins to have many branches. Both the semiclassical (dashed) and the exact quantum results for the last frame are shown enlarged at the bottom. integrable dynamics. Figure 4 shows a successful implementation of cellular dynamics to the highly nonlinear long time dynamics of a Morse oscillator from ref 14. The CD method is not specifically a guided wave packet theory, although it can be used to propagate wave packets. It is an example of what is called a “uniform” approximation to semiclassical amplitudes, in which an explicit integration is introduced that is not itself Gaussian. If the integration is done instead by stationary phase or steepest descent then the result is a “primitive” semiclassical approximation.20 (If the integrand is Gaussian, then the dynamics is linear and stationary phase is exact). This would typically be a less accurate approach, although the advantages of a coherent state (i.e., Gaussian) primitive representation may be retained at the cost of complex “roots” or stationary phase points. We discuss this next. 3.1. Coherent State Semiclassical Representation. The CD method implemented the first of two approaches that we mentioned for the application of the VVG propagator to wave packets, that is, a uniform and numerical approach. The stationary phase, fully primitive semiclassical approach is discussed next, leading to the coherent state semiclassical representation. Suppose we label the phase space location γt ) (pt,qt), and the Gaussian wave packet (coherent state) ψt(q) in eq 1 as ψt(q) ) 〈q|γt〉, where we have suppressed the spread parameter Rt ) at + ibt in the notation. We then can investigate the VVG propagation of |γ0〉 followed by its projection onto |γf〉:
〈γf|e-iHt/p|γ0〉 ≈
∫∫〈γf|q〉GVVG(q,q′,t)〈q′|γ0〉 dq dq′
(8)
VOL. 39, NO. 2, 2006 / ACCOUNTS OF CHEMICAL RESEARCH 129
Guided Gaussian Wave Packets Heller
with GVVG(q,q′,t) given by eq 4, and the integrals are done by stationary phase. Carrying out the stationary phase gives the simultaneous equations for the root solutions q and q′:
-R0(q - q0) - ip0 + ipt(q,q′,t) ) 0 -Rf(q - qf) + ipf -ipt(q,q′,t) ) 0
(9)
where
pt(q,q′,t) )
∂S(q,q′,t) ∂q
p0(q,q′,t) ) -
∂S(q,q′,t) ∂q
(10)
These equations, discussed first by Klauder8 and later by Weissman,7 have an interesting interpretation, as discussed in the first application of this approach to nontrivial potential problems.21,22 We recast eqs 9 in the form
j 0 + ip j 0 ) R0q0 + ip0 R0q Rfq j t - ip j t ) Rfqt - ipt
(11)
j 0 is the where q j t is the stationary phase value of q, and q stationary phase value of q′. Equations 11 have the interpretation of the permitted complexification of the Gaussian parameters (pt,qt) and (p0,q0), which leave the respective Gaussians untouched, except for an overall factor. That is, we may take q0 to an arbitrary real or complex value q j 0, and if p0 is adjusted to p j 0 according to eq 11, the Gaussian ψ0(q) stays the same. This freedom permits a root search, which may be viewed as finding an (complex) initial classical guiding trajectory with a j 0), becoming (p j t,q j t) under real time evolution. “center” (p j 0,q j t) is equivalent The complex root has been found if (p j t,q in the sense of eq 11 to the real values (pf,qf) of the final j 0) that satisfy, coherent state. The set of parameters (p j 0,q for example, the first eq 11 is called the equivalent initial complex manifold. The root search is then to find trajectories connecting the initial and final equivalent manifolds. This is subject to the usual difficulties of root searches and indeed of Stokes phenomena, as discussed in refs 21 and 22. Note that the essence of this approach is to tailor a propagator to the given initial and final states; there are many propagators in a certain sense. We call this approach the coherent state semiclassical approach, CSSA. It is easy to show that the initial conditions can become complex even in the case of linear dynamics, as in the harmonic oscillator potential. However in this case one gets the same (exact) result whether or not the complex guiding center is used. 3.2. Nearly Real Method. The notion that the important contributions to the primitive coherent state semiclassical amplitudes are classical trajectories that are nearly real was developed in ref 23. This works remarkably well and considerably simplifies implementation of the primitive coherent state approach. The approach was further refined, tested, and given a formal basis in terms of similarity transformed dynamics in ref 24. The procedure is to 130 ACCOUNTS OF CHEMICAL RESEARCH / VOL. 39, NO. 2, 2006
choose several initial conditions close to the classical center of the initial coherent state and propagate them forward in time using the classical Hamiltonian; then for each time that one of the above trajectories comes near the final point, perform a search to find the nearest solution to the boundary conditions, eq 11. Next, for each solution, one takes small steps forward and backward in time to trace out a branch of solutions. For each new time, the solution at the previous time is used as a guess for the nearly real solution. Finally, one adds all the contributions from the different branches. 3.3. Frozen Gaussians and Herman-Kluk. The frozen Gaussian approximation25 (FGA) was introduced as an attempt to simplify the numerics of the TGA and at the same time take nonlinearities into account. The numerical simplification was achieved by throwing out the van Vleck determinant and individual wave packet spreading, which costs some effort to compute, especially in many dimensions. This was patched up by expanding the initial Gaussian in terms of many other Gaussians, each of a different initial position and momentum. These are then propagated according to their different initial centers by the FGA, and the resulting sum of all the frozen Gaussians is taken to be the wave function. Nonlinearities are thus included, as different guiding trajectories follow the correct dynamics. The method is somewhat ad hoc and surely not accurate for long times. Nonetheless, it is reasonably accurate at short times and has the advantage that it is easy to implement for many body systems. It has been used, for example, to study electronic transitions in solution.26,27 Herman and Kluk28 (HK) found a new way to exploit frozen Gaussians, albeit at the expense of reintroducing elements of the van Vleck determinant. There has been some controversy surrounding the derivation, validity, and use16,29-31 of the HK, but it remains a very popular choice. It is essentially an “outside the box” semiclassical method, approaching the VVG as p f 0. This is an example of a very promising general plan: why not invent methods that agree with VVG in the limit (as they must, because VVG is correct in the limit p f 0) but are something different for finite p? Kay has investigated some of these issues and others surrounding IVR methods in an influential theoretical and numerical study.31,32 The interested reader is referred to these papers, and also refs 16, 29, and 30 for more details. Manolopoulos developed a hybrid of CD and HK in an interesting study of Franck-Condon spectra.33 Very sophisticated, chemically relevant, and remarkable “on the fly” wave packet approaches to multisurface nonadiabatic dynamics, coordinated with electronic structure calculations, have been made by Martinez and co-workers.1 3.4. Off-Center Guiding. 3.4.1. Background. We now are ready to discuss the rather different assumptions behind off-center guiding of Gaussian wave packets. The first thing to mention is that all the guided wave packet approaches mentioned so far have propagated Gaussians according to their center in phase space. As mentioned, even the CSSA adheres to this by complexifying the center
Guided Gaussian Wave Packets Heller
q)δ(p′ - p) is propagated, considering (p′,q′) as an initial condition, becoming δ(qt(p′,q′) - q)δ(pt(p′,q′) - p). Then the initial distribution Fi(p,q) ≈ exp[-R(q - qi)2 1/R(p - pi)2] becomes
Ft(p,q) ≈
∫∫ dp′ dq′ exp[-R(q′ - qi)2 1 (p′ - pi)2 δ(qt(p′,q′) - q)δ(pt(p′,q′) - p) (13) R
]
and the overlap with the final Gaussian density becomes
FIGURE 5. Off-center TGA approximation shown in gray, optimized for the overlap with the final Gaussian shown as a circle on the upper right. The trajectory at the center of the “+” was used to guide the Gaussian; it traveled to the “+” shown in the gray propagated state. It is a real trajectory. without changing the Gaussian. The complex center refers back to the original, real center through eq 11. Like the CSSA, off-center guiding tailors propagators to the amplitude being computed, but unlike the CSSA or NR-CSSA it uses purely real trajectories not related to the center of the wave packet around which to linearize the dynamics. The off-center approach has already achieved success under some of the most extreme circumstances. We discuss the applications the method has seen and the prospects for future work. 3.4.2. Off-Center Method. The basic idea motivating off-center guiding can be seen in Figure 5.2,3 The centerguided TGA linearization of the whole wave packet is a very poor approximation in the tails of the phase space distribution when nonlinear evolution has become important within the initial zone of size h surrounding the guiding trajectory. However the dynamics is still locally linear over domains smaller than h. Taking advantage of this permits a much better approximation to overlaps, such as that shown in Figure 5 between the distorted propagated state and the Gaussian used for final projection (black circle upper right). If the off-center trajectory marked with “+” is used as a guide of the initial wave packet (black circle to the left), it becomes the gray Gaussian state shown. This is actually a good approximation to the full nonlinear propagation of the Gaussian in the zone that matters, near the final Gaussian. The question immediately arises as to how to choose the off-center trajectory. We give here a heuristic derivation giving the best choice of guiding trajectories, which are taken to be purely real. We already know that a stationary phase attack on the problem starting with the VVG propagator results in complex trajectories and the CSSA approach, so this is not the approach we want here. Instead, we first construct the classical Gaussian phase space density corresponding to the initial and final Gaussian wave packets, for example,
1 Fi(p,q) ≈ exp -R(q - qi)2 - (p - pi)2 R
[
]
(12)
with a similar expression for Ff(p,q). (We assume for simplicity that the R parameter is real and the same for both initial and final state). Each initial trajectory δ(q′ -
O(t) ≈
∫∫ dp′ dq′ exp[-R(q′ - qi)2 - R1(p′ - pi)2] 1 exp -R(qt(p′,q′) - qf)2 - (pt(p′,q′) - pf)2 (14) R
[
]
We next evaluate this integral by steepest descent, which here amounts to finding the stationary points of the exponent and expanding quadratically around them, followed by doing the resulting Gaussian integral exactly. Let (pt+,qt+) be such a stationary point. It corresponds a local maximum in the product of the two densities, and is the choice for the end point of the guiding trajectory. Following it back in time to t ) 0, the point (p0+,q0+) is the guiding center around which the wave packet is propagated to get an approximation to that local contribution to the overlap. The main point of this derivation is that using the trajectory (p0+,q0+) as a guiding center for the OC-TGA gives a quantum amplitude with a magnitude that is the square root of the classical overlap as determined by steepest descent. This is easily seen by noting that quantum overlaps and Wigner phase space overlaps are identical for Gaussians and that steepest descent is equivalent to a quadratic expansion of the exponent. Since the classical procedure seems intuitively optimal for computing the overlap by steepest descent and the guiding center initial condition (p0+,q0+) corresponds to it in the usual sense of the square root of the classical probability, we have found an optimal trajectory for the OC-TGA. The approach is quite different than CSSA or the “nearly real” approach, since these refer back to the center of the initial Gaussian, whereas in the off-center guiding, we have abandoned the original center. However offcenter guiding shares the property that the propagator is optimized for the final state in question. 3.4.3. Numerical Applications. Off-center guiding was first used in connection with a completely chaotic stadium billiard problem.2,3 (It was further developed by Tomsovic and co-workers, who applied it to the highly singular and nonlinear Coulomb problem.4,5) The autocorrelation function of the wave packet launched in the stadium billiard, shown in Figure 6, was computed. The structure of phase space (Figure 7) is typical of chaotic systems, with a socalled “homoclinic tangle” surrounding the periodic orbit, which is at the center of the chosen wave packet. Homoclinic orbits are trajectories that start near the periodic orbit (and in fact asymptotically approach it going backward in time), then leave its vicinity, only to return and VOL. 39, NO. 2, 2006 / ACCOUNTS OF CHEMICAL RESEARCH 131
Guided Gaussian Wave Packets Heller
FIGURE 6. Schematic of the initial Gaussian wave packet |φ〉, which when propagated in the stadium billiard shown gives the autocorrelation function 〈φ|φ(t)〉 of Figure 9. The wave packet is initiated in the middle of the billiard with momentum in the direction of the arrow.
FIGURE 9. Autocorrelation function 〈φ|φ(t)〉 for the stadium wave packet shown in Figure 6: quantum result, solid line; semiclassical, dashed line.
FIGURE 10. Example of nonlinear dynamics and choice of guiding reference trajectories from ref 4. The ellipse is the region that is mapped into the gray area under the dynamics; the inner lobe is on its third return from the Coulomb singularity; the outer lobe is on its second return. The white dotted line within the black ellipse is mapped into the dotted lines within the gray zone. The two intersections are chosen as off-center guiding trajectories. FIGURE 7. Nonlinear homoclinic tangle developing for a generic chaotic system in the vicinity of a periodic orbit. The autocorrelation function can be represented as a sum of discrete overlaps with each of the approximately horizontal manifolds cutting through the initial Gaussian shown as a grayscale density; each of those contributions is in turn computed by the overlap with a thin Gaussian (an example is shown only at the bottom for clarity) obtained by the off-center propagation of the initial wave packet. For example, the contribution marked with an asterisk is computed from an off-center trajectory marked with a “+”.
FIGURE 8. A homoclinic orbit associated with the horizontal bounce unstable periodic orbit and the amplitude it contributes to the total autocorrelation function shown in Figure 9 are shown. From ref 3. asymptotically approach it again. That central periodic orbit is special only in that it never leaves; it’s contribution is added to the homoclinic ones. Each homoclinic orbit may be justified as a guiding center, along the lines just given. A single homoclinic orbit and the amplitude it contributes are seen in Figure 8. The exact quantum and semiclassical correlation functions are shown in Figure 9. The total time elapsed was enough to require about 30 000 homoclinic contributions, at the end of the run 132 ACCOUNTS OF CHEMICAL RESEARCH / VOL. 39, NO. 2, 2006
near t ) 6, which corresponds to about four periods of the central periodic orbit. There are several remarkable aspects of this calculation. Perhaps the most surprising is that the semiclassical result, no matter how it is computed, is so accurate or even close to the right result for such long times in a strongly chaotic system. After many years of controversy about whether the correspondence principle was applicable to chaotic systems, spurred, for example, by questions about the convergence of the Gutzwiller trace formula,34 in the energy domain, this direct construction of the quantum time-dependent amplitudes from the chaotic dynamics represented a new watershed. It was also possible to construct a nearly resolved quantum energy spectrum from the Fourier transform of the quantum autocorrelation function.2,3 Tomsovic and co-workers applied offcenter propagation to integrable nonlinear dynamics in the Coulomb potential.4,5 The method of choosing offcenter reference trajectories was simplified, as shown in Figure 10. Oddly, this system is a more severe test of the approach than the chaotic ones, since the stretching is not so strong as to make the overlap regions nearly linear. Nonetheless, the results are extremely good, as shown in Figure 11. 3.5. Comparison of Approaches. Suppose we attempted the stadium chaos calculation with the HK method. There would have to be five or so frozen Gaussians in each homoclinic branch near the initial state. If this was done by starting with enough trajectories initially, one would require at least a million frozen
Guided Gaussian Wave Packets Heller
FIGURE 11. Autocorrelation function |c(t)| for the wave packet shown in Figure 10 for the first 22 Kepler periods of motion. The exact quantum result is dashed. From ref 4.
FIGURE 12. The initial Gaussian (gray circle) and one of many homoclinic recurrences; the trajectory at q ) , p ) 0 returns to the vicinity at q ) 0, p ) -µ, bringing with it a stretched piece of the initial wave packet shown as a dashed line. The q-axis is taken to be the unstable manifold. Gaussians at each time, perhaps 10-50 times this altogether, in the regime where we required 30 000 guided Gaussians. Of course, one could begin to take the structure of phase space into account and drastically reduce this number. However, one of the advantages of the HK method would be lost, that is, it’s “set it and forget it” mode of running a single large set of trajectories and not worrying too much about where they go. Figure 13 illustrates some of the aspects of the HK and OC-TGA approaches. The CSSA, NR-CSSA, and OC-TGA approaches all tailor their propagators to the initial and final state pair. They do not provide a sum over wave packets, which is collectively the whole wave function. However, this difference is not as big as it might seem at first. If one wants to know the wave function at position b q, then only those frozen Gaussians with centers near b q contribute. One must make sure there are enough frozen Gaussians present to overlap each other in phase space, not just in coordinate space. It is reasonable to expect that the NR-CSSA would do just as well as the OC-TGA in the strong chaos limit. In fact, it is possible to show that the “nearly real” contributions are very nearly real and almost equivalent to the OCTGA, in the following way. Consider a homoclinic recurrence branch as shown in Figure 12. Without loss of
FIGURE 13. (A) The FGA and HK approaches expand the initial Gaussian (grey circle) in terms of others (open circles) typically as shown. The distorted gray area is the time evolution of the gray circle. The frozen Gaussians remain circular, their centers attracting along the stable axis of motion and expanding along the unstable axis. (B) The OC-TGA approximation typically uses the unstable axis and expands along that; expansion in the perpendicular direction is redundant. The pluses mark the off-center guiding trajectories. The TGA guided wave packets extend along the unstable axis; there are fewer of them to begin with and the growth in the number needed is much slower due to the extrapolation along the unstable direction. Two of the guided Gaussians are shown, one of which is the TGA center-guided Gaussian. Each Gaussian is used only when it is optimal in a given region; they are not added together. generality, we take the momentum axis as the stable axis of motion, along which trajectories are approaching each other exponentially, and the position axis as the unstable axis. The initial Gaussian we take at q ) p ) 0. The homoclinic trajectory shown arises from the point (q ) , 0); it guides a wave packet with real off-center trajectories to a recurrence at (q ) 0, -µ) along the stable axis. Expanding
qf ) 0 +
∂q ∂q (q - ) + (p - 0) + ... ∂q0 0 ∂p0 0
pf ) -µ +
∂p ∂p (q - ) + (p - 0) + ... ∂q0 0 ∂p0 0
(15)
where (q0,p0) are the initial position and momentum of the guiding trajectory. The alignment of the manifolds and their stabilities dictates however,
∂p ∂q ≈ 0; ≈0 ∂p0 ∂p0
(16)
Now let us inquire what the CSSA complex valued solution might be, corresponding to this branch of the j 0 ) iR, that dynamics. Using eq 11, we set q j 0 ) , then p is, the initial momentum is pure imaginary. However, eq VOL. 39, NO. 2, 2006 / ACCOUNTS OF CHEMICAL RESEARCH 133
Guided Gaussian Wave Packets Heller
15 together with eq 16 shows that this small momentum change has virtually no effect on pf or qf. The real OCguided trajectory is then necessarily almost the complex solution, which is thus nearly real.
4. Conclusion In this Account, we have emphasized the off-center guiding approach to semiclassical Gaussian wave packet dynamics (off-center thawed Gaussian approximation or OC-TGA). To develop the method and compare, we have derived and motivated several other strategies for semiclassical Gaussian propagation. The OC-TGA probably holds the record for being able to handle strong chaos and many (thousands of) branches to the classically contributing dynamics. Along with the nearly real approach to the CSSA, it deserves consideration as a method of choice for handling semiclassical approximations to nonlinear dynamical systems for long times.
References (1) See, for example, Coe, J. D.; Martinez, T. J. Competitive Decay at Two- and Three-State Conical Intersections in Excited-State Intramolecular Proton Transfer. J. Am. Chem. Soc. 2005, 127, 4560-4561 and references therein. (2) Tomsovic, S.; Heller, E. J. Semiclassical Dynamics in the Stadium: Unexpected Long Time Accuracy. Phys. Rev. Lett. 1991, 67, 664-667. (3) Tomsovic, S.; Heller, E. J. The Long-Time Semclassical Dynamics of Chaos: The Stadium Billiard. Phys. Rev. 1993, E47, 282-299. (4) Suarez Barnes, I. M.; Nauenberg, M.; Nockleby, M.; Tomsovic, S. Semiclassical theory of quantum propagation: The Coulomb potential. Phys. Rev. Lett. 1993, 71, 1961-1964. (5) Suarez Barnes, I. M.; Nauenberg, M.; Nockleby, M.; Tomsovic, S. Classical orbits and semiclassical wave packet propagation in the Coulomb potential. J. Phys. A 1994 , 27, 3299-3321. (6) Heller, E. J. Time-dependent Approach to Semiclassical Dynamics. J. Chem. Phys. 1975, 62, 1544-1555. (7) Weissman, Y. Semiclassical approximation in the coherent states representation. J. Chem. Phys. 1982, 76, 4067-4079. (8) Klauder, J. R. Path integrals and stationary-phase approximations. Phys. Rev. D 1979, 19, 2346-2356. (9) Heller, E. J. Spectra and Eigenstates Associated with Periodic Orbits. In Advances in Classical Trajectory Methods; Hase, W.; Ed.; JAI Press: Greenwich, CT, 1991; Vol. 1, pp 165-213. (10) Heller, E. J. Semiclassical Wave Packets. In The Physics and Chemistry of Wave Packets; Yeazel, J.; Uzer, T. eds, Wiley: New York, 1999. (11) Heller, E. J. Wave Packet Dynamics and Quantum Chaology. Lectures in the 1989 NATO Les Houches Summer School on Chaos and Quantum Physics; Giannoni, M.-J., Voros, A., ZinnJustin, J., Eds.; Elsevier Science Publishers B. V.: Amsterdam, 1991; p 547.
134 ACCOUNTS OF CHEMICAL RESEARCH / VOL. 39, NO. 2, 2006
(12) Heller, E. J. The Semiclassical Way to Molecular Spectroscopy. Acc. Chem. Res. 1981 14, 368-375. (13) van Vleck, J. H. The Correspondence Principle in the Statistical Interpretation of Quantum Mechanics. Proc. Natl. Acad. Sci. U.S.A. 1928, 14, 178-188. (14) Heller, E. J. Cellular Dynamics: A New Semiclassical Approach to Time Dependent Quantum Mechanics. J. Chem. Phys. 1991, 94, 2723-2729. (15) Miller, W. H. Comment on: Semiclassical time evolution without root searches. J. Chem. Phys. 1991, 95, 9428-9430. (16) Miller, W. H. The semiclassical initial value representation: a potentially practical way for adding quantum effects to classical molecular dynamics simulations. 2001, 19, 2942-2955. (17) Berry, M. V.; Balazs, N. L.; Tabor, M.; Voros, A. Quantum Maps. Ann. Phys. (N.Y.) 1979, 122, 26-63. (18) Sepœlveda, M. A.; Heller, E. J. Semiclassical calculation and analysis of dynamical systems with mixed phase space. J. Chem. Phys. 1994, 101, 8004-8015. (19) Sepœlveda, M. A.; Heller, E. J. Semiclassical analysis of hierarchical spectra. J. Chem. Phys. 1994, 101, 8016-8027. (20) Miller, W. H. Classical S Matrix: Numerical Application to Inelastic Collisions. J. Chem. Phys. 1970, 53, 3578-3587. (21) Huber, D.; Heller, E. J. Generalized Gaussian wave packet dynamics. J. Chem. Phys. 1987, 87, 5302-5311. (22) Huber, D.; Heller, E. J.; Littlejohn, R. G. Generalized Gaussian wave packet dynamics, Schro¨ dinger equation, and stationary phase approximation. J. Chem. Phys. 1988, 89, 2003-2014. (23) van Voorhis, T.; Heller, E. J. Nearly real trajectories in complex semiclassical dynamics. Phys. Rev. A 2002, 66, No. 050501. (24) van Voorhis, T.; Heller, E. J. Similarity transformed semiclassical dynamics. J. Chem. Phys. 2003, 119, 12126-12153. (25) Heller, E. J. Frozen Gaussians: A Very Simple Semiclassical Approximation. J. Chem. Phys. 1981, 75, 2923-2931. (26) Turi, L.; Rossky, P. J. Critical evaluation of approximate quantum decoherence rates for an electronic transition in methanol solution. J. Chem. Phys. 2004, 120, 3688-3699. (27) Neria, E.; Nitzan, A. Semiclassical evaluation of nonadiabatic rates in condensed phases. J. Chem. Phys. 1993, 99, 1109-1123. (28) Herman, M. F.; Kluk, E. A semiclassical justification for the use of nonspreading wave packets in dynamic calculation. Chem. Phys. 1984, 91, 27-34. (29) Baranger, M.; de Aguiar, M. A. M.; Keck, F.; Korsch, H. J.; Schellhaass, B. Semiclassical approximations in phase space with coherent states. J. Phys. A 2001, 34, 7227-7286. (30) Grossmann, F.; Herman, M. F. Comment on Semiclassical approximations in phase space with coherent states J. Phys. 2002, A35, 9489-9492. (31) Kay, K. G. Integral expressions for the semiclassical time dependent propagator. J. Chem. Phys. 1994, 100, 4377-4392. (32) Kay, K. G. Numerical study of semiclassical initial value methods for dynamics. J. Chem. Phys. 1994, 100, 4332-4445. (33) Walton, A. R.; Manolopoulos, D. E. A new semiclassical initial value method for Franck-Condon spectra. Mol. Phys. 1996, 87, 961-978. (34) Gutzwiller, M. C. Chaos in Classical and Quantum Mechanics; Springer-Verlag: New York, 1991.
AR040196Y
