J . Phys. Chem. 1988, 92, 3225-3235 transitions in the absorption envelope; Le., eq 10. Calculating these transition energies quantum mechanically is extremely difficult due to the large number of states which must be handled. However, given the current understanding of classical and semiclassical mechanics, it may be fundamentally impossible to calculate the transition energies from a classical trajectory study. In general, the determination of semiclassical transition energies for polyatomic systems (i.e., more than three dimensions) has required quasi-periodic classical motion. For highly excited overtone states the classical intramolecular motion is usually chaotic, which is the case for all of the benzene CH(D) overtone states. From a time-dependent point of view, recurrences in the probability of occupying the initially prepared state gives rise to
3225
the fine structure in the overtone absorption spectrum. Though rudiments of these recurrences may be present in the short-time trajectory P(n,t), the chaotic classical motion destroys the longer time recurrences, which occur quantum mechanically. It is these latter recurrences which are needed to evaluate fine details in the absorption spectrum. Thus, as applied in this work, the classical trajectory method appears limited to the evaluation of low-resolution absorption spectra.
Acknowledgment. This research was supported by the National Science Foundation. We thank Professors G. S. Ezra and Y . T. Lee for very helpful discussions. Registry No. C6H6, 71-43-2; D,,7782-39-0.
On the Relationship between the Classical, Semiclassical, and Quantum Dynamics of a Morse Oscillator Jeffrey R. Reimerst Department of Theoretical Chemistry F l l , University of Sydney, Sydney, New South Wales 2006, Australia
and Eric J. Heller* Physics and Chemistry Departments BG- 10, University of Washington, Seattle, Washington 981 95 (Received: July 31, 1987)
The relationship between the classical, semiclassical, and quantum dynamics of the Morse oscillator is investigated, and the almost classical, frozen Gaussian approximation is shown to give an accurate solution of a differential equation similar to the Schriidinger equation that shares its eigenfunctions with the Schriidinger equation. This result is used to justify a posteriori the success of the De Leon-Heller spectral quantization method for determining eigenfunctions and eigenvalues from frozen Gaussian dynamics, and techniques are described for the optimization of such calculations. A close relationship is shown to exist between the various dynamical schemes, and physical problems that reveal or cloak this relationship are described.
I. Introduction A large part of the study of chemical physics is concerned with problems that are so complex that accurate quantum mechanical calculations are impractical. Our hope is that many systems might be treated by semiclassical means. The paradigm for the semiclassical correspondence is the harmonic oscillator. Here, the coherent stateI4 undergoes the same dynamics5-' using classical mechanics (Le., the Liouville equation) as using quantum mechanics (Le,, the Schrodinger equation), as exemplified by Ehrenfest's theorem.8 Much e f f ~ r t ~ *has ~ -been ' ~ extended to find analogous coherent states for anharmonic problems, but these states rarely have all of the desired properties. Problems for which such states have been found include the particle in a boxI7 and the two-dimensional rigid-rotor" problems. Another school of thought prefers to use the coherent states, or Gaussian wave packets. Properties of anharmonic quantum systems may be deduced from the classical or semiclassical dynamics of these wave packets. Widely differing views are held as to the merits of this approach, and it is the purpose of this paper to resolve some of these differences. We take a semiclassical approach, using the tools of Gaussian wave packet dynamicsIH1 to relate the classical to the quantum dynamics of the Morse oscillator, and show how it arises that some properties are poorly reproduced while other properties are closely related. This paper asserts that there does exist a strong correspondence and that this correspondence can be exploited to develop computational algorithms: it provides a n ~ t h e r ' ~ - ' *a*posteriori ~* justification of the De Leon-Heller spectral quantization method.31
The main argument against the use of Gaussian wave packets refers to their behavior on one-dimensional anharmonic potential (1) Glauber, R. J. Phys. Rev. Lett. 1963, 10, 84. (2) Glauber, R. J. Phys. Rev. 1963, 131, 2766. (3) Glauber, R. J. Phys. Rev. 1963, 130, 2529. (4) Schrodinger, E. Natunvissenschaften 1926, 14, 664. (5) Saxon, D. S. In Elementary Quantum Mechanics; Holden Day: San Francisco, 1964; p 158. (6) Nieto, M. M.; Simmons, L. M., Jr. Phys. Reu. A 1979, 19, 438. (7) Nieto, M. M.; Simmons, L. M.. Jr. Phvs. Reu. Lett. 1978. 41. 207. (8) Merzbacher, E. In Quantum Mechanics; Wiley: New York, 1970; p 41. (9) Nieto, M. M.; Simmons, L. M., Jr. Phys. Reu. D 1979, 20, 1321. (10) Nieto, M. M.; Simmons, L. M., Jr. Phys. Reu. D 1979, 20, 1332. (1 1) Nieto, M. M.; Simmons, L. M., Jr. Phys. Reu. D 1979, 20, 1342. (12) Nieto, M. M. Phys. Reu. D 1980, 22, 391. (13) Nieto, M. M.; Gutschick, V. P. Phys. Rev. D 1981, 23, 922. (14) Gutschick, V. P.; Nieto, M. M.; Simmons, L. M., Jr. Phys. Lett. A 1980, 76A, 15. (15) Perelomov, A. M. Commun. Math. Phys. 1972, 26, 222. (16) Barut, A. 0.;Girardello, L. Commun. Math. Phys. 1971, 21, 21. (17) Reimers, J. R.; Heller, E. J. J . Phys. A 1986, 19, 2559. (18) Reimers, J. R.; Heller, E. J. J . Chem. Phys. 1985, 83, 511. (19) Heller, E. J. J . Chem. Phys. 1975, 62, 1544. (20) Heller, E. J. Chem. Phys. Lett. 1975, 34, 321. (21) Heller, E. J. J . Chem. Phys. 1976, 64, 63. (22) Heller, E. J. J . Chem. Phys. 1976, 65, 4979. (23) Heller, E. J. J . Chem. Phys. 1977, 67, 3339. (24) Heller, E. J. J . Chem. Phys. 1978, 68, 3891. (25) Davis, M. J.; Heller, E. J. J . Chem. Phys. 1979, 71, 3383. (26) Heller, E. J.; Stechel, E. B.; Davis, M. J. J . Chem. Phys. 1980, 73, 4720. (27) Heller, E. J. J . Chem. Phys. 1981, 75, 2923. (28) Davis, M. J.; Heller, E. J. J . Chem. Phys. 1981, 75, 3916.
'Queen Elizabeth I1 Fellow.
0022-3654/88/2092-3225$01.50/0 , ~~I
I
1
~
0 1988 American Chemical Society
3226
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Reimers and Heller
surfaces. On harmonic surfaces wave packets move periodically phenomena such as gas-phase scattering,22neutron ~ c a t t e r i n g , ~ ~ and atom surface s ~ a t t e r i n g and , ~ ~ may even be applied at high in time, their centers moving in accordance with Newton’s laws: temperatures. It is amazing, however, that good results are also initial Gaussian form is maintained but they may oscillate in width obtained (using the FGA) for calculated eigenfunctions and ei(breathe) Here, two parallels between classical and quantum g e n v a l u e ~ ’ ~ . ’when ~ , ~ typically ~ , ~ ~ , ~the ~ semiclassical wave packet mechanics are drawn: first, the evolution of the phase space dynamics from which the eigenfunctions are derived (by Fourier probability distribution function is identical when propagated using either the time-dependent Schrodinger equation or the classical transform, (2.7)) poorly represents the exact wave packet dynamics. Accurate eigenfunctions are even calculated in the chaotic Liouville equation; second, the evolution of the quantum amplitude regime.s2,s3 How is it then that a technique apparently based function is given as a function of just four parameters, for a upon short-time assumptions can describe so many different one-dimensional Gaussian, all of which are calculable by using phenomena including some long-time effects? pure classical mechanics. The second parallel is the most imHerman and helped to resolve this paradox by expanding portant one as it allows the solution to any quantum problem to the initial wave packet in terms of a swarm of weighted frozen be expressed in terms of classical quantities; knowing just the time Gaussians. They showed that the errors introduced by keeping dependence of the probability distribution is not sufficient to do the wave packets frozen subsequently canceled so that a quite naive this, but it is a start.43 Fortunately, classical behavior is mainapproach gives good answers. Unfortunately, their results do not tained for short times on anharmonic surfaces, but at longer times account for the success of studies made using a single, unweighted the wave packets distort, interfere, bifurcate, and may even reform wave packet. at very long t i m e ~ . ~ ~ -Thus, ~ O in order to determine exactly the In some recent studies of the two-dimensional rigid r ~ t o r ’ ~ , ~ ~ quantum time dependence of the coherent states, a many-paand the particle in a box,17 it was demonstrated that the thawed rameter model is needed whose parameters do not necessarily have Gaussian approximation gives exact eigenvalues, eigenfunctions, classical analogues. As a result, the second parallel above appears and spectra. This is a surprising result in view of the fact that lost if the time scale for wave packet breakup is short compared a localized Gaussian wave packet spreads and becomes nonwith the time scale of physical interest (typically 0.1-10 periods). Gaussian rather quickly, but this problem is rectified by use of The correspondence between the time evolution of the probability the method of imagesss to enforce the boundary conditions. These distributions can last for much longer but even here studies also showed that the frozen Gaussian approximation gives quantum interference effects can cause discrepancies to appear at sort timess0 as well as bring about long-time r e c u ~ r e n c e s . ~ ~ - ~exact ~ spectral intensities and eigenfunctions, thoughjt incorrectly gives equally spaced eigenvalues. A method for determining For many real problems, it is often the case that the wave packet eigenvalues and eigenfunctions is the De Leon-Heller spectral breaks up within one p e r i ~ d , ~displaying ’,~~ interference effects, quantization ~ c h e m e , which ~’ determines the eigenfunctions I$n) and then it is tempting (though, as this paper shows, incorrect) semiclassically using frozen Gaussians and then deduces the eito conclude that the classical quantum correspondence is comwhere H i s the Hamiltonian. genvalues E, using E, = ($mlfl$n), pletely lost. Here, we show that the correspondence is merely Clearly, this method produces the exact eigenvalues for the obscured, not totally lost. two-dimensional rigid-rotor and particle in a box problems as for Semiclassical Gaussians are based upon classical guidance of these problems I$,) is exact and so E , must also be exact. This quantum wave packets, extending the analogy seen above for the is an amazing result which, as we shall see, has its analogue for harmonic oscillator. Here, two methods are in common use: the more general potentials as well. It is the basis of the success of frozen Gaussian approximation (FGA), in which the wave packets the De Leon-Heller method and justifies the use of frozen move but do not spread,27and the more general thawed Gaussian Gaussians to solve these problems. approximation (TGA), in which spreading is a l l ~ w e d . ’One ~ could, In section I1 of this paper, we construct the analogy between on the basis of the preceding discussion, expect these approxithe particle in a box and rigid-rotor problems and the problem mations to be useful only in restricted circumstances. Indeed, if of a continuously differentiable potential. Some test calculations, the “essential physics” of the process occurs in short times, then performed upon the Morse oscillator, are described in section 111, the TGA may be used directly, and this approach has worked well and the results are given in section IV. Section V concerns itself for molecular s p e c t r o s c ~ p yand ~ ~ for * ~short-lived ~ ~ ~ ~ ~ collision ~~ with optimizing calculations via Gaussian wave packet dynamics. I
11. The Analogy (29) Lee, S.-Y.;Heller, E. J. J . Chem. Phys. 1982, 76, 3035. (30) Tannor, D. J.; Heller, E. J. J . Chem. Phys. 1982, 77, 202. (31) De Leon, N.; Heller, E. J. J. Chem. Phys. 1983, 78, 4005. (32) Heller, E. J. Faraday Discuss. Chem. SOC.1983, 75, 141. (33) Blanco, M.; Heller, E. J . J. Chem. Phys. 1983, 78, 2504. (34) Reimers, J. R.; Wilson, K. R.; Heller, E. J. J . Chem. Phys. 1983, 79, 4749. (35) Drolshagen, G.; Heller, E. J. Chem. Phys. Lerr. 1984, 104, 129. (36) Heller, E. J.; Reimers, J. R.; Drolshagen, G. Phys. Rev. A , in press. (37) Littlejohn, R. G. Phys. Rep. 1986, 138, 193. (38) Littlejohn, R. G. Phys. Rev.Lett. 1986, 56, 2000. (39) Jackson, B.; Metiu, H. J. Chem. Phys. 1985,82, 5707. (40) Jackson, B.; Metiu, H . J. Chem. Phys. 1985, 83, 1952. (41) Jackson, B.; Metiu, H. J. Chem. Phys. 1986, 84, 3535. (42) Herman, M. F.; Kluk, E. Chem. Phys. 1984, 91, 27. (43) Brown, R. C.; Heller, E. J. J. Chem. Phys. 1981, 75, 186. (44) Alterman, E. B.; Tahk, C. T.; Wilson, D. J. J . Chem. Phys. 1966,44, 451. (45) Baetzold, R. C.; Tahk, C. T.; Wilson, D. J. J . Chem. Phys. 1966, 45, 4209. (46) Endres, P. F.; Wilson, D. J. J. Chem. Phys. 1967, 46, 425. (47) Stechel, E. B.; Schwartz, R. N. Chem. Phys. Lerr. 1981, 83, 350. (48) Brickmann, J.; Russegger, P. J. Chem. Phys. 1981, 75, 5744. (49) Brickmann, J. J. Chem. Phys. 1983, 78, 1884. (50) Thirumalai, D.; Bruskin, E. J.; Berne, B. J. J . Chem. Phys. 1985, 83, 230.
For the two-dimensional rigid rotor18 and the particle in a box,I7 the appropriate Gaussian-based coherent states I$) may be expressed in terms of the appropriate eigenstates Ix,) as
so that the exact Schriidinger dynamics of the coherent states I$s(O)) = I$) becomes
I$,(
t ) ) = Ccne-iEsnr/hlXn ) n
(2.2)
where the coefficients c, are time independent and Es, are the appropriate quadratically spaced energy levels. Frozen Gaussian dynamics, is not an exact solution of the Schradinger wave equation. However, a different wave e q ~ a t i o n ’ ~ can ~ ~ *always *~* be found whose solutions k ( t ) ) = Cc,e-lELnr’hlxn) n
(2.3)
(51) Reimers, J. R.; Wilson, K. R.; Heller, E. J.; Langhoff, S. R. J . Chem.
~~.
Phvs. ,-1985. - ~82. , 5064 . ~
~~~
(52) Heller, E. J. Phys. Rec;. Lert. 1984, 53, 1515. (53) Davis, M. J.; Heller, E. J . J. Chem. Phys. 1984, 80, 5036. (54) Reimers, J. R.; Heller, E. J. J. Chem. Phys. 1985, 83, 516. (55) Thomson, W. Br. Assoc. Adv. Sci. Rep. 1947, 17, 6
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3227
Dynamics of a Morse Oscillator share the same spatial eigenstates as the Schrodinger equation, and whose energy levels ELnare linearly spaced. The linear spacing forces the dynamics of this equation to be periodic, and indeed, this dynamics differs from the frozen Gaussian dynamics by a t most a phase factor y:
E F , = EL, = 7, I + d t ) ) = e-’Y‘/hl+L(t))
(2.4)
Clearly, the Schrodinger dynamics (+s(t) ), the enforced periodic dynamics I#L(t)),and the frozen Gaussian dynamics I+F(t)) are closely related as they all are based upon the same eigenstates. These eigenstates dephase differently as time evolves, and thus the observed dynamics, as seen via projections such as (qI+F(t)) and (ql+s(t)), appears unrelated and the link between the semiclassical and quantum dynamics is obscured. However, one can extract identical eigenstates from the results of either the Schriidinger dynamics (by inverting (2.2)) lxn) a
JmeiELnl/hl+s(r)) -m
(2.5)
the new enforced periodic dynamics (by inverting (2.3))
~ x , )a JmeiEL*f/hI+L(t))
(2.6)
-m
or the frozen Gaussian dynamics Ixn)
0:
Jme’EFnt/hl+F(t)) -_
(2.7)
and then, if desired, obtain the eigenvalues of the Schrodinger equation using
ESn = (XnlHIXn)
(2.8)
For more general one-dimensional problems, it is always possible to find a wave equation whose solutions can be written in the form of (2.3), where, as before, c, are constants, E,, are equally spaced, and Ix,) are the appropriate eigenfunctions. Again, I+L(t)) undergoes enforced periodic dynamics, with the period being 2 s / a where CY is the energy level spacing, but (2.4) in general will not hold for constant c, so that no analytical relationship exists between the frozen Gaussian dynamics and the enforced periodic dynamics. Even so, both dynamics are periodic, and if the spacing a is chosen so that the two periods are identical, then one expectss6that l+F(t)) = I+L(t)),and that the eigenfunctions evaluated from (2.6) and the energy levels evaluated from (2.8) are excellent approximations to the exact quantum solutions. In practice, it is not necessary to construct the wave equation that has as its solution the enforced periodic dynamics of (2.3). This construction may be performed, however, by expanding the eigenvalues ES, of the Schriidinger equation in a Taylor series of order m :
any two parameters at will: here, we select CY such that the enforced periodic dynamics has the same period as the frozen Gaussian dynamics, i.e. CY
= 27rh/T
(2.13)
where T i s the classical period; and as y serves only to offset the energy spectrum, its value does not effect the eigenfunctions calculated using (2.6). We choose to select y such that the strongest line in the spectrum has the same energy in both cases, Le., for this line, EL, = ES,. L i t t l e j ~ h n ’has ~ pointed out that there exists an alternate analogy to the one developed here, which for one-dimensional problems leads to the same results as above. It is, however, more easily generalizable to multidimensional problems. H e suggests the use of the good actions instead of the Hamiltonian H . If, for some continuous functionf, H L = f(H), then HL will have the same eigenfunctions as H itself, but the eigenvalues will differ so that different Schriidinger dynamics is produced. Choosing fcorrectly can cause the eigenvalues of HL to be linearly spaced, and one would expect frozen Gaussian dynamics to give an excellent description of this Schrodinger dynamics. The basic approximation in semiclassical mechanics is that the eigenvalues of the action operator I are linearly spaced, and so one should choose f such that Z = HL =f(H). In this paper, we do not take this approach but instead link the frozen Gaussian dynamics to the dynamics of a different wave equation and use the original Hamiltonian H . 111. Test Calculations for the Morse Oscillator
The Morse oscillator provides a good test problem as the q u a n t ~ m and~ classical ~ ~ ~ ~dynamics can both be solved analytically; the potential is applicable to most diatomic molecule^;^^^^^ and much is already known about the quantum dynamics of wave packets on this potential surface.6,46,48,49,61We choose the coefficients in the Morse potential so as to closely approximate62 the excited B electronic state of 1,63964and write the Hamiltonian as
H’ =
‘/2p’(f)2
+ D,’[1 - e-y(prd]2
(3.1)
where De’ = 4362.70 cm-’, p’ = 1.7523 .kl,re = 3.018 68 %.,and the reduced mass of I2 is p’ = 63.4523 amu. The calculations are more conveniently posed in terms of reduced variables, and we introduce reduced positions, momenta, and energies defined by q = ( r - re)(p’wgl/h)1/2,
P = p ’ ( p ’ ~ ~ ’ / h ) -=’ /(2D,)-1/2 ~ (3.2)
m
Es, = C Pknk
(2.9)
k=O
(for the Morse oscillator, the particle in a box, and the rigid rotor, m = 2; for the harmonic oscillator, m = l ) , and by expressing the wave equation in the form (2.10) The m
+ 1 coefficients ak are determined by writing EL,, = an + y m
(2.1 1)
k=O
m
+ y ) l k = k=O Cfiknk
(2.12)
+
for all n. Equating powers of n gives rise to m 1 equations relating the m 3 unknowns C Y , y, and ak to the coefficients Pk. The problem is underdetermined and thus we are free to choose
+
(56) Weissman, Y.; Jortner, J. Phys. Left. A 1979, 70A, 177.
(3.3)
H = H’/(hW,’), De= D,’/(hwQ’)
(3.4)
and
respectively, where wgl = (2D,’p’z/p’)1/2
and substituting this equation along with (2.3) and (2.9) into (2.10) yielding
Eak[-i(cYn
p = f(wglh/M’)-’/2
(3.5)
is the angular frequency of vibration for infinitesimal displacements from equilibrium, wgl/(2sc) = 119.31 cm-’. In reduced units, (57) Morse, P. M. Phys. Rev. 1929, 70, 222. (58) ter Haar, D. Phys. Rev. 1946, 34. 57. (59) Herzberg, G. In Molecular Spectra and Structure. I . Spectra of Diatomic Molecules, 2nd ed.; Van Nostrand: Princeton, 1965; p 101. (60) Lippincott, E. R.; Schroeder, R. J . Chem. Phys. 1955, 23, 113 1. (61) Walker, R. B.; Preston, R. K. J . Chem. Phys. 1977, 67, 2017. (62) Reimers, J. R.; Heller, E. J., manuscript in preparation. (63) Barrow, R. F.; Yee, K. K. J . Chem. Soc., Faraday Trans. 2 1973,69, 684. (64) Luc, P. J . Mol. Spectrosc. 1980, 80, 41.-
3228 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Reimers and Heller t-
.....,.......j II 1 J
E204 ................................... .........................
t i
I:
n
ii
t- T/2
j
................................................. De
t= 3T/4
l
t= T
~
A
I
t-0.
1
i
t-
.
4T
1 l
-
!
t- T/B
t- T/4
I
d
'-10
i\
II I
.......
I
t- T / 4
T/8
10
30 4
50
73
Figure 1. The cutoff Morse potential.
-
= h = wo = 1, and De = 36.5649. The Schrdinger equation has 73 bound-state eigenvalue^,^^*^^ and these are quadratically spaced with p
E,, = De - (2De - )/, - n),/(4De)
Schradinger dynamics.
(3.6)
Expressions for the associated eigenstates are available.s7@ In addition to these bound states, there exists a continuum of unbound eigenstates, and while these functions can be expressed in terms of Whitaker f ~ n c t i o n s , 6they ~ ~ are ~ ~ difficult to normalize and use. It is possible to expand the continuum states in terms of some (harmonic oscillator) basis set46 or to add extra terms to the potential48so as to produce discrete, normalizable states. We choose to select only those Whittaker functions that are stationary at some large distance qc = 74.6457 (rc = 8 A), truncating this set at n = 207, EZO7= 61.70. Figure 1 shows the potential function used, with these extremum points marked. Trajectories are run on this potential at two different classical energies, EL = 9.1412 and EH = 26.6905, and these energies and the associated classically allowed regions are also shown in Figure 1. The periods of trajectories at these energies are TL = 7.2553 and TH = 12.0910, respectively. At each energy, trajectories are run starting either on the inside wall qL, the equilibrium separation, qc = 0, or the outside wall, qo. For the trajectories at the lower energy, the second derivatives of the potential vanish a t q = qo, so significant anharmonicity is involved even for these trajectories. In all cases, the initial state is chosen to be a Gaussian wave packet,
( s / + ( O ) = exPNat(q - qJZ + P,(q - 4,) + r,l/hJ
Figure 2. Time evolution of the position space probability distribution for the HOW initial conditions: (- -) thawed Gaussian dynamics; (-) upper two rows are enforced periodic dynamics and lower two rows are
(3.7)
where (q,,p,) is the location of the center of the wave packet at time t , a,is the width parameter, and y r is a normalization and phase factor. a. is another arbitrary parameter, and we use three different values in these calculations corresponding to "wide", "harmonic", and "narrow" Gaussians. The "harmonic" Gaussians have a. = i/2, so that this function is the appropriate harmonic-oscillator coherent state for a harmonic oscillator fitted at the potential minimum. For the high-energy trajectories, the wide wave packet has a width that is matched to the period of these trajectories, TH. Note that in order to calculate the B X absorption spectrum of I, at 0 IC,one would run a trajectory at E = EH and q = qI,but the initial wave packet would have the width a. = 0.897ito match the width of the ground vibrational state of the X electronic state.
-
IV. Results Of the Dynamics' Figures and show the time position space probability density Obtained by using the various These are (a) exact quantum propagation (b) thawed Gaussian dynamics, and (c) enforced pe-
(65) Skodje, R. T . Chem. Phys. Lett. 1984, 109, 221. ( 6 6 ) Lennard-Jones, J. E.; Strachan, C. Proc. R. SOC.London, A 1935,
..
150, 442.
( 6 7 ) Strachan, C. Proc. R. SOC.London, A 1935, 150, 456. (68) Slater, L. J. In Handbook of Mathematical Functions; Abramowitz, M.; Stegun, I., Eds.; Dover Publications: New York, 1965; p 503.
i ;' n
0
0
20
0
0
20
28
0
20
i
4
Figure 3. Time evolution of the position space probability distribution for the HIW initial conditions: ( - - - ) thawed Gaussian dynamics; (-) upper two rows are enforced periodic dynamics and lower two rows are
Schradinger dynamics. riodic dynamics. The bottom two rows of each figure compares the Schrijdinger density, 1(q1+s(t))12(solid line), to the thawed Gaussian density, I( q1+T(t))12(broken line). Each frame shows the probability density at different times after the propagation started, and this time is expressed as a fraction of the period of the classical trajectory. Similarly, the top two rows of each figure compares the enforced periodic density, I(q ( $ L ( f ) ) (solid line), to the thawed Gaussian density. These two figures differ in that we varied the initial starting position and the initial width of the Gaussian. Of the nine calculations performed a t high energy ( E = EH), the calculation in which the propagation algorithms produce the most similar results is the trajectory HOW (High energy, Outside starting position, Wide initial Gaussian), and these results are shown in Figure 2. Here, the only significant differences are observed at times of half-integral periods, when the classical trajectory is at the inner turning point. In this region, the anharmonicity of the potential is large, and this generates interference effects in the Schradinger and enforced periodic dynamics which take the wave function well out of the domain of Gaussian functions. Some similaritv remains. however. and all functions are dominated by a sharpcentral peak. The'thawed Gaussian dynamics, however, overestimates the sharpness, and the peak cut off in Figure 2 actually goes far off scale. After integral numbers of periods, the thawed Gaussian to almost its original shape, and the enforced periodic dynamics to exactly its original shape, but the Schrijdinger dynamics continues to deteriorate with the wave packet breaking up significantly within four periods. The least similar results are found for the trajectory HIW (High energy, Inside starting position, Wide initial Gaussian), and the results are shown in Figure 3. Here, the thawed Gaussian tracks the Schrodinger dynamics quite well for the first half of a period,
The Journal of Physical Chemistry, Vol. 92, No. 11, I988 3229
1
2
3
1
2
3
1
2
3
t / T
Figure 6. Time evolution of the standard deviation in position for the low-energy trajectories: (-- -) thawed Gaussian dynamics; (-) enforced
Figure 4. Time evolution of the average position for the high-energy trajectories: (- - -) thawed Gaussian dynamics; (-) enforced periodic dynamics; (--) Schriidinger dynamics.
1
2
3
1
2 t/'
periodic dynamics; (- -) Schrodinger dynamics.
3
Figure 5. Time evolution of the standard deviation in position for the
high-energy trajectories: (- - -) thawed Gaussian dynamics; (-) enforced periodic dynamics; (--) Schrodinger dynamics.
and both wave functions become quite broad. Soon, the tail of these wave functions approaches the steep potential-wall part of the Morse function with the result that the Schrodinger dynamics banks the wave function up near the wall. The thawed Gaussian spreads in accordance with the local curvature of the potential a t the center of the wave packet. It does not recognize the wall at this time and continues to spread far into the classically forbidden region. This problem may be avoided by using a propagation algorithm69 for thawed Gaussian using the average potential, P = ($Iy$)and its derivatives instead of Vin the trajectory calculations, but this introduces new problems with a Schrixlinger wave function poorly represented within the important bowl region of the potential. Thus, a significant improvement is not obtained. Again, after one period the thawed Gaussian almost reforms to its original shape, but Schriidinger propagation has damaged the wave packet too severely and it continues to striate. Initially, the enforced periodic dynamics quickly breaks up the wave packet, with the high-energy components leaving the main body of the wave packet. That the high-energy components are more prominent here than in the Schriidinger dynamics is due to the energy level linearization process: high quantum states are pushed to much higher energy. The main body of the wave function does move coherently, however, and the entire wave packet reforms every period, so the enforced periodic dynamics remains much easier to model than the Schrodinger dynamics. For all of the high-energy trajectories, the average position
4 = (Ic.Iql$)
(4.1)
and its standard deviation
(69) van Leuven, P. In Dynamics of Wave Packets in Molecular and Nuclear Physics; van Leuven, P., Ed.; Springer-Verlag: Berlin, 1986; p 174.
.
are plotted as functions of time in Figures 4 and 5 , respectively. Because the continuum eigenstates are forced to be stationary at the boundary q = qB,some artificial reflection of the high-energy components of the Schrodinger dynamics occurs. 4 and oq are particularly sensitive to this, and thus we damp the contributions from the continuum eigenstates toward the classical values as the amplitude moves to large distance. The average positions are always well behaved: the thawed Gaussian and enforced periodic dynamics both give periodic results. These results differ most at half-periods, with the differences being slight for the HOW trajectory and largest for trajectories that have significant contributions from non-bound eigenstates. N o periodicity occurs in the Schrodinger dynamics as the everbroadening wave function continually occupies more and more of the potential well and the average position inside the well tends to a constant. Further, the dissociative states contribute a term which, after the bowl region is initially cleared, grows linearly with time. It is possible for the Schrodinger dynamics to show pse~dorecurrences~~ at very long times, but this usually does not happen on time scales of physical interest and is made less likely as the dissociative character of the initial wave packet increases. Looking at the standard deviations, Figure 5 , amplifies the effects seen in Figure 4. Now, only the enforced periodic dynamics is periodic, but, even with this restriction, large fluctuations of the width within one period may occur. The width produced from the Schrodinger dynamics approaches the width of the bowl of the potential if the wave packet is comprised largely of bound eigenstates, and otherwise increases linearly at long times. Unfortunately, the results from the thawed Gaussian approximation are erratic, showing, in the worst cases, wave packets changing width by X200 within TH/20units of time. Such behavior stems from the conditional stability (in a Floquet sense) of the classical equations of motion for a nearby trajectoryzz and gives rise to the linear increase in the amplitude of the oscillations in the width.50 From Figures 2-5 it is clear that thawed Gaussians describe the Schrodinger dynamics quite accurately for times of at least half of a vibrational period, but very poorly for longer times. This behavioral change over time is problem dependent, becoming shorter as the anharmonicity increases. In Figure 6 are given the standard deviations uq for the nine trajectories run at low energy E = EL. Here, the effective anharmonicity is considerably reduced, and the thawed Gaussian approximation can be quite accurate for times of up to four vibrational periods. Many problems often require trajectories of much lower energy than E,, probing regions of the potential surface for which the effective anharmonicity is quite small, and the thawed Gaussian approximation is then expected to be accurate for much longer times. The thawed Gaussian approximation will always be useful for problems, such as the calculation of absorption band contours, which require only short-time information, and for these problems it is useful for its speed, accuracy, and physical intuition.34 Frozen Gaussian dynamics is not as accurate as thawed Gaussian dynamics for short times, but it is much better behaved
3230
TABLE I traj HOW HOH HON HCW HCH HCN HIW HIH HIN LOW LOH LON LCW LCH LCN LIW LIH LIN
Reimers and Heller
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
E
40
26.6905 26.6905 26.6905 26.6905 26.6905 26.6905 26.6905 26.6905 26.6905 9.1412 9.1412 9.1412 9.1412 9.1412 9.1412 9.1412 9.1412 9.1412
16.4764 16.4764 16.4764 0
0 0 -5.281 1 -5.2811 -5.281 1 5.9277 5.9277 5.9277 0 0 0 -3.4674 -3.4674 -3.4674
PO
0 0 0 7.3062 7.3062 7.3062 0 0 0 0
0 0 4.2758 4.2758 4.2758 0 0 0
soli 0.2598 0.5 2.5115 0.2598 0.5 2.5115 0.2598 0.5 2.5115 0.26 0.5 1.5 0.26 0.5 1.5 0.26 0.5 1.5
V'o 1.06 1.06 1.06 0 0 0 -13.55 -13.55 -13.55 2.14 2.14 2.14 0 0 0 -6.41 -6.41 -6.41
V"o
(dao/dr)o
EL
Es
-0.10 -0.10 -0.10 1 1 1 5.02 5.02 5.02 0 0
0.19 0.55 12.67 -0.37 0 12.11 -2.38 -2.01 10.10 0.14 0.5 4.5 -0.37 0 4 -1.37 -1 3
26.8 26.9 28.0 27.7 28.5 34.7 35.9 34.1 29.3 9.3 9.4 9.9 9.8 9.7 10.2 10.9 10.3 10.2
26.8 26.9 27.9 27.3 27.2 27.7 27.9 28.0 28.2 9.3 9.4 9.9 9.8 9.6 9.9 10.5 10.1 10.1
0 1 1 1 3 3 3
at longer times. The average position, q, is identical for both propagation algorithms, and for the frozen Gaussian approximation, uqis a constant. After every period, the enforced periodic dynamics and the frozen Gaussian dynamics must produce identical results, and in between times, the differences between the two are typically relatively small. Note that the frozen Gaussian dynamics, when compared directly to the Schrodinger dynamics, appears poorly correlated. It is only when one remembers that the Schriidinger and enforced periodic dynamics are intimately related that one sees that a close relationship does in fact exist between the SchrMinger and frozen Gaussian dynamics. Another propagation algorithm that is often used to propagate wave functions is pure classical mechanics, in the form of the Liouville equation.46 We performed such calculations for these wave packets, but the results are not given here since much about this aspect of the relationship between classical and Schriidinger mechanics is already known.48-49 For these problems, both propagation algorithms give quite similar results for times of up to four periods. The agreement is best for wave packets that contain significant dissociative character, but note that if any quantum pseudorecurrences do occur, then the classical calcubtion can not reproduce them. This propagation algorithm is not investigated further as it is purely classical and the objects involved do not lend themselves to quantum theories as do Gaussian wave packets. B. Overlaps. While the projection of the dynamics onto position space is an important physical quantity, many physical properties such as spectra and scattering cross sections are more directly related to projections (or overlaps) of the wave functions onto other states. Thus,it is important to study such functions before making general conclusions about the relationship between classical and quantum mechanic^.^' One of the simplest overlap functions is the autocorrelation function (4.3)
whose Fourier transform a(w) = JieiW'C(t) dt
(4.4)
gives the eigenvalues of the Hamiltonian (see next section). This overlap function is also closely related to overlap functions that appear in the theory of one-photon s p e c t r o s ~ o p y . It ~ ~has been calculated for the 18 trajectories described earlier and is given in Figure 7 for the H O W and H I W trajectories. Here, the top panels compare the results from the enforced periodic and frozen Gaussian dynamics, while the bottom panels compare the results from the Schriidinger and thawed Gaussian dynamics. In the last section, the H O W trajectory was seen to give the most similar dynamics under the various propagation algorithms. (70)Heller, E. J . J . Chem. Phys. 1978, 68,2066.
GS GL
ETS
0.95 0.98 0.47 0.45 0.23 0.27 0.31 0.14 0.62 0.94 0.96 0.99 0.88 0.83 0.72 0.60 0.70 0.84
0.02 0.02 1.48 2.62 5.39 6.02 3.98 5.87 0.55 0.02 0.02 0.08 0.04 0.06 0.15 0.35 0.33 0.07
0.99 1.00 0.89 0.99 0.97 0.85 0.81 0.88 0.97 0.99 0.98 0.98 1.00 1.00 0.96 0.92 0.97 1.00
UTS
0.02 0.04 0.10 0.15 0.13 0.10 0.11 0.10 0.15 0.01 0.01 0.01 0.07 0.09 0.11 0.11 0.12 0.12
~ F L
0.014 0.040 0.080 0.002 0.008 0.014 0.009 0.006 0.002 0.002 0.014 0.076 -0.003 0.000 0.012 0.019 0.007 0.006
nmax
35 35 36 36 30 21 22 27 36 9 10 10 11 8 5 6 8 10
1.1 0.8 1.8 3.8 5.2 11.2 12.2 9.6 4.5 2.1 1.6 1.4 2.3 3.0 5.3 7.0 4.7 2.7
Figure 7. Time autocorrelation functions obtained by propagating the wave packets with the HOW or HIW initial conditions using the SchrGdinger, enforced periodic, thawed Gaussian, and frozen Gaussian
propagators. The Schrodinger autocorrelation function initially decays quite quickly i s the wave packet gains momentum, recurring at times near classical periods. Because the wave function slowly spreads in time, the recurrences in the overlap become less exact and last for longer times. Figure 2 shows that the thawed Gaussian wave packet mimics the SchrGdinger wave function quite poorly at half-integral periods, but quite well at integral periods. Here, the overlap at half-integral periods is very small in both cases, hiding the errors in the thawed Gaussian propagation. It is only after integral numbers of periods that the overlap is significant, and as this is where both algorithms give similar wave functions, the thawed Gaussians perform quite well for appreciable times. The same argument applies also to the comparison between the enforced periodic and frozen Gaussian overlap functions: the overlaps with the initial wave packet are significant only after integral numbers of periods, and it is at the time of maximum overlap that both propagation algorithms must give exactly the same results. Thus there is excellent agreement between the two overlap functions. Also, in the last section, the H I W trajectory gives the least similar dynamics under the various propagation algorithms. Here, this is reflected in the comparison of the Schriidinger and thawed Gaussian propagation algorithms. The Schriidinger wave packet breaks up quite quickly (see Figure 3), and after only half of a classical period it has significant amplitude in the region of the initial wave packet. Its autocorrelation function thus reflects the details of the quantum interference effects that take place near the potential wall. The thawed Gaussian is incapable of mimicking such effects, and its autocorrelation function is restricted to low-amplitude, low-width recurrences at times near the classical period. No such dramatic differences arise in the comparison of the results of the enforced periodic and frozen Gaussian wave functions, however. As is seen for the HOW trajectory, both wave
Dynamics of a Morse Oscillator
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3231 x
HIW
L
Ill
I
n 0
u 40
20
10
30
10
x)
50
60
30
W
Figure 8. Energy spectra obtained by Fourier transform of the autocorrelation functions shown in Figure 7.
functions are periodic and their maximum values must coincide; thus the agreement between the two propagation algorithms is, again, excellent. A summary of the results obtained from all 18 trajectories is given in Table I. Included here is the average unphased difference between the Schrijdinger and thawed Gaussian autocorrelation functions
and the average unphased difference between the enforced periodic and frozen Gaussian autocorrelation functions
No surprises are found in these results: the thawed Gaussians perform reasonably only when they are started from the outer side of the Morse potential bowl; and the frozen Gaussians perform poorly only when the dynamics leaves the initial region very slowly. C. Eigenvalues. The energy levels that generate the dynamics can be extracted from the dynamics by Fourier transforming the autocorrelation function, (4.4). Using (2.2) and (2.3), this transform is evaluated analytically to give US(@)
=
CICnI’ n
UL(@)
=
CICnI’ n
6 ( -~ E s n / h )
(4.7)
-ELn/h)
(4.8)
and
respectively. Both spectra consist of 6 functions located at the energy levels that generated the trajectory. As the frozen Gaussian dynamics is periodic in time with period T, its Fourier transform u ~ ( w must ) be linearly spaced with spacing 27r/T, and as the enforced periodic dynamics is constructed to include exactly this energy level spacing, the two spectra must be very similar. The energy levels differ by at most a constant, and to an excellent degree of approximation, this constant is attributable to the failure of the frozen Gaussian dynamics to deduce the effective zero-point energy, r / T . It is possible to correct for this effect, and OF(@)
= UF(W- n / T
+ Im ao)
(4.9)
has the same energy levels as uL(w),but this is at the price of displacing the band contour so that the intensity distribution is less accurate. Note that, in multidimensional problems, it is not always convenient to find the vibrational periods, and thus this correction is not so easy to apply. While the application of this correction increases our conceptual understanding, in the next section it is shown that one does not need to apply this correction in order to obtain good eigenvalues or eigenfunctions from frozen Gaussian dynamics. The results of numerically Fourier transforming the autocorrelation functions for the H O W and H I W trajectories are given in Figure 8. Here, the autocorrelation functions are truncated
at t = 20r, and a Gaussian spectral window is used in the transform to force a resolution of 0.1 between adjacent lines in the calculated spectra. For the HOW trajectory, the four calculated spectra are nearly indistinguishable. The thawed Gaussian and Schrodinger eigenvalues differ on average by only 0.01, well within the spectral resolution. However, it is clear from magnified spectra that the nonlinearity seen in the thawed Gaussian eigenvalues is significantly less than that seen in the Schrijdinger eigenvalues. All four band contours are very similar, a consequence of the very similar initial decays seen in Figure 7, and the fact that, for this trajectory, Im a. = 7r/T so that up(@) = u i w . Quite different results are obtained for the H I W trajectory. Here, the Schriidinger eigenvalues decrease in spacing, becoming unresolvable from each other at an energy just below the dissociation energy, De = 36.5649. Some structure reappears in this spectrum at very high energy because the continuum eigenfunctions are represented as a discrete (particle in a box type) set of eigenstates. The thawed Gaussian spectrum resembles this spectrum only in that the two band contours are nearly identical. Very few lines can be resolved in the thawed Gaussian spectrum, a fact attributable to the low area of the recurrences of the overlap seen in Figure 7. While the agreement for the H I W trajectory between the spectra from the enforced periodic and frozen Gaussian dynamics is not perfect, it is a dramatic improvement over the agreement seen between the Schrodinger and thawed Gaussian dynamics. Here, lines of the correct shape and location are produced, but it now is the band contour that is in error. Figure 8 shows how the linearization process increases the average energy of the enforced periodic wave packet and thus why, in Figure 3, the high-energy fragment of the wavefunction disappeared so quickly. Quantitatively, Table I gives values for the average energies
Es =
CICnI’ESn n
(4.10a)
EL =
CICnI’ELn
(4.10b)
and n
for the Schrodinger and enforced periodic dynamics, respectively, and these averages differ significantly for this trajectory. The skewness of Figure 8 means that the bulk of the dynamics takes place at low energy, and thus the bulk of the wave function leaves the initial region slower than classical mechanics indicates, and it oscillates with a smaller amplitude. In fact, the center of the main peak seen in Figure 3 moves like a classical particle a t E = 19, slaved to move with a longer period. The initial sluggishness of the motion leads to CL(I)decaying slower than CF(t), as seen in Figure 7, and, to turn this argument around full circle, to the sharper main band seen in Figure 8. A summary of the results obtained from all 18 trajectories is given in Table I. Two indices are used to gauge the accuracy of the thawed Gaussian results: ETS
=
CICnIZIESn n
- ETnl
(4.1 la)
is used to gauge the error in the calculated eigenvalues, where ET^ are the energies at the center of the peaks in u ~ ( w ) and , UTS
=
xcn*CTn
(4.1 1b)
n
for the error in the peak amplitudes, where CT: are the areas of the peaks in nT(w), and 1 = CnCTn2. The frozen Gaussian energies are in error by the constant a / T - Im a,,, and the index UFS
=
xcn*CFn n
(4.12)
is used to gauge the error in the peak amplitudes, where CF,Zare the areas of the peaks in C T ~ ( W ) ,and 1 = No surprises are seen in these results: always, the frozen Gaussians give good results whereas the thawed Gaussians give both good and poor results.
3232 The Journal of Physical Chemistry, Vol. 92, No. 1 I, 1988
Reimers and Heller
TABLE I1 traj HOW HOW HOW HOH HON HCW HCH HCN HIW HIW HIW HIW HIW HIH HIN LOW LOH LON LCW LCH LCH LCH LCH LCN LIW LIH LIN
Cn/Chx
n 33 35 37 3s 3s 3s 35 35 25 30 35 40 45 35 35
( XnlXFn)
0.91 0.92 0.89 0.92 0.87 0.93 0.93 0.88
0.74 1.oo 0.66 1.oo 1.oo 0.99 0.99 1.oo 0.93 0.98 1.oo 0.98 0.92 1.oo 1.oo 1 .oo 0.99 1.oo 0.96 0.62 1.oo 0.98 0.60 1.oo 1.oo 0.98 0.96
10 10 10 10 6 9
10 13 10
10 10 10
0.01 0.79 0.92 0.59 -0.19 0.94 0.98 0.97 0.98 0.97 0.98 0.94 0.98 0.98 0.90 0.97 0.97 0.98 0.97
Esn
EFn
-
EDHn
0.01
25.83 26.88 27.89 26.88 26.88 26.88 26.88 26.88 21.05 24.14 26.88 29.19
0.00 0.03 0.24 2.26 0.00 0.25 2.26 0.64 0.15 0.00 0.19 0.72 0.20 2.25 -0.17 0.07 1.07
31.36 26.88 26.88 9.75 9.75 9.75 9.75 6.21 8.88 9.75 12.25 9.7s 9.7s 9.75 9.75
0.63 0.07
0.51 -0.03 -0.03 0.07 -0.03 0.18 -0.01 -0.03 0.08 0.07
-0.03 -0.03 0.07 ' I \ '
"
7
13 17 11 12 12 9 10 48 21 13 29 67 8 11 8 6 6 9 8 3 6 15 6 9 7 6
0.44 0.02 0.08 0.35 -0.27 0.02 0.00 0.06
0.14 0.07 0.07 0.16 1.07 -0.17 0.07 1.07 '
14
0.06
-0.17
I
i-"
- ESn
0.00 0.00 -0.02
I
'
"
'
1
'
I
i
ill ili
A
-10
0
10
t
20
v I
8
,
10
# T20I
1
3
-10
20
10
10
20
9
Figure 9. Position space representationof the eigenfunctions calculated from the HIW initial conditions: (-- -) calculated using either the thawed or frozen Gaussian approximations; (-) exact eigenfunctions.
Figure 10, Position space representation of the eigenfunctions calculated from the HOW initial conditions: (- - -) calculated using either the thawed or frozen Gaussian approximations; (-) exact eigenfunctions.
D . Eigenfunctions. From the dynamics, the eigen functions of the Hamiltonian can be retrieved by projecting either the Schriidinger dynamics,
(4.13) or the enforced periodic dynamics, (4.14) a t the appropriate energy Es, or ELnof the peaks in as(w) and uL(w), respectively. Projecting the thawed and frozen Gaussian dynamics in a similar fashion3' leads to approximate expressions for the eigenstates
! -5
C
,
,
,
,I ,
5
1,
,
,
.
,
3
,
,
,
,
,
5
! , , , j :c
9
Figure 11. Position space representation of the eigenfunctions calculated from the LCH initial conditions: (-- -) calculated using either the thawed or frozen Gaussian approximations; (-) exact eigenfunctions.
and (4.16) respectively, where NT and NF are normalizing constants. It is important to note that the energies which appear in these two equations are precisely the values obtained from uF(w)and u ~ ( w ) , not the exact energy levels. Thus, errors made in calculating EF,
and ET^ do not directly contribute to these calculations, and corrections such as (4.9) must not be applied. Figures 9, 10, and 11 show eigenfunctions obtained form the HOW, HIW, and LCH trajectories, respectively. These figures show ( q l x n ) ,(qIXFn),and (qlXTn)for the eigenfunctions corresponding to the most intense lines in uF(w),and also ( q l X n ) and ( &Fn) for some other intense lines. Also, Table I1 gives details
Dynamics of a Morse Oscillator of the strength of each line n relative to the strongest line; the overlap ( xnlxFn),the Schrijdinger energy level Esn,and the average error in the position of the nodes relative to the average internodal spacing
where qruand are the positions of the ith nodes in (qlXn)and (qIXFn),respectively. For the H O W and LCH trajectories, accurate eigenfunctions are calculated from the frozen Gaussian dynamics for all states with strong spectral lines. The errors are least for the strongest lines and increase as the line strength decreases. This trend occurs because the strongest lines correspond to states with action similar to the action of the classical trajectory, and these semiclassical techniques effectively expand the energy as a linear function of action about the classical action.31 Always, the calculated eigenfunctions underestimate the amplitude in the region of the potential wall, a consequence of the wave packet penetrating this barrier rather than banking up at the wall. Also, the calculations overestimate the amplitude in the region of the outer turning point, a consequence of the classical dynamics spending too much of its time in this region. Both of these effects are known, and if necessary, corrections can be applied.28 The calculated nodal positions are, in general, just a fraction below the exact nodal positions: this effect is presumably a consequence of the misrepresentation of the outer regions of the eigenfunctions. The situation is somewhat different for the H I W trajectory as now the energy spectrum has a large bandwidth. States like n = 35 which lie close in energy to EH remain well represented, but 0 n)= 25 and n other states which have strong lines in ~ ~ (like = 45 are poorly represented, and it seems that n = 30 and n = 40 are the most distant states for which good eigenfunctions can be obtained. The problem arises because the distant eigenstates have very different spatial extents to the central eigenstate which is approximately the allowed region of the classical trajectory. This is demonstrated in Figures 9-1 1 by the arrows which indicate the classical turning points of the trajectories. In Figure 10, the eigenfunction with n = 45 has its outside lobes far past the classical turning point, and one could not expect good answers from the frozen Gaussian dynamics. In contrast, the state with n = 35 is considerable narrower than the range of the trajectory, and so there must be a complete cancellation of the frozen Gaussian density in the region where the frozen Gaussian spends most of its time, and again, one would not expeccreasonable results. That a state has an intense spectral line is thus a necessary condition for the frozen Gaussian dynamics to give accurate eigenfunctions, but it is not a sufficient condition and also the change in the quantum number must not be too great. For many practical problems, including most problems involving light atoms such as carbon and hydrogen, the energy spread is quite small and thus frozen Gaussians should give good results for all eigenstates with strong spectral lines. In Table I1 is given a summary of all of the eigenfunctions shown in Figures 9-1 1, and for the eigenfunctions with n = 35 obtained from all trajectories at high energy as well as for all of the eigenfunctions with n = 10 obtained from all of the low energy trajectories. For the central lines, the overlaps ( xnlxFn)are all quite good, being the poorest for the trajectories starting with narrow Gaussians. Generally speaking, the wide Gaussians give the best results at large q, the narrow Gaussians marginally give the best results at low q, and the Gaussians that start with the width of the harmonic oscillator fitted to the potential minimum give the best results in the intermediate regions. Analysis of the nodal errors verifies these conclusions. When eigenfunctions are extracted from the thawed Gaussian dynamics, only poor functions result, as seen in Figures 9-1 1. Always, the thawed Gaussian represents the region near the potential wall very badly, but the eigenfunctions do improve at large q. Because the thawed Gaussian approximation gets worse at long times, the results from thawed Gaussian calculations usually decrease in quality as the truncation time t , (see (4.17)) increases,
rn
The Journal of Physical Chemistry, Vol. 92, No. 11. 1988 3233 converging to very poor results. Even for the trajectories at low energy, the frozen Gaussian results are far superior to the thawed Gaussian results. E . The De Leon-Heller Spectral Quantization Method. This is a technique3I designed to calculate accurate eigenvalues and eigenfunctions for the Schrdinger equation from frozen Gaussian , dynamics. Autocorrelation functions CF(t),spectra u ~ ( w ) and eigenstates IxFn)are calculated by using the techniques previously described, and then accurate eigenvalues are obtained by using (4.18)
The basic idea is that, while the peaks in ~ ~ (are0linearly ) spaced and a poor approximation to Es,, IxFn)are an excellent approximation to Ix,), and thus EDHnare expected to be quite accurate. Table I1 gives the error in the eigenvalues calculated by using this equation, EDHn- Esm for all of the eigenfunctions calculated in the last section, as well as the error in the eigenvalues as calculated from u ~ ( w )EFn , - Esn. In general, a dramatic improvement in the calculated eigenvalues is observed when the De Leon-Heller method is applied. The errors in EDHn are all less than 0.1 (l%), except for the states from the HIW trajectory whose eigenfunctions are poorly represented (see Table I1 and Figure 10). The success of the De Leon-Heller method, demonstrated both here and e l ~ e w h e r e , ’ indicates ~ J ~ ~ ~ ~that a strong correlation between classical and quantum mechanics exists even for problems with anharmonic potentials.
V. Optimizing Calculations The frozen Gaussian approximation is exact when the potential function is where ai = Im ao,as, for this case, the frozen Gaussian is in fact the appropriate harmonic oscillator coherent state. If one is concerned with finding a good local propagator (that is, a propagator that, given a wave packet at time t , does the best possible job of constructing the wave packet at time t + dt), then this approximation will be good if the average potential difference (5.2) and it higher moments are small. At this, the thawed Gaussian approximation does much better than the frozen Gaussian approximation, but as the previous section indicates, such locally optimized propagators do not lead to the best long-time dynamics. To minimize the moments of the potential difference, one makes the wave packet as narrow as is possible in position space. A narrow position space distribution implies, however, a wide momentum space distribution, and so after a short period of time, the distribution becomes wide in position space and the moments of the potential difference become large. Thus, there exists some optimum initial wave packet width, and local considerations such as (5.2) cannot be used to find it. In section IV, the most accurate results always come from wave packets with narrow energy distributions, us(w): an index of this spread is the standard deviation of the energy US
= ([CIcn12Esn21- [ C I C ~ I ~ E S ~ I ~ I * ’ * (5.3) n
n
and this is given in Table I for the 18 trajectories run. When the energy spread is narrow, the enforced periodic wave function sticks together, its center moving along with the center of the frozen Gaussian. Then, the autocorrelation function decays quickly, hiding the mid-period differences between the enforced periodic and frozen Gaussian dynamics. Thus the calculated spectra are all rather similar, and all of the eigenstates calculated for states with strong spectral lines are accurately represented. In contrast, wave packets with a large energy spread have high-energy components which quickly leave the main body of the wave function, and the skewness of the energy distribution means that the bulk of the enforced periodic wave function moves with smaller amplitude and the guiding classical trajectory. This skewness is
3234 The Journal of Physical Chemistry, Vol. 92, No. 1 1 , 1988 indicated in Table I by the energy expectation values Es and E , as well as by nmax,the quantum number of the strongest line in u S ( w ) . When the energy distribution is narrow, E, and ELapproach E ai, and nmaxapproaches the central line, n = 35 for the high-energy trajectories or n = 10 for the low-energy trajectories. Because of the slackness of the main body of the enforced periodic wave function, its autocorrelation function decays slower, and its spectra are highly non-Gaussian in shape. Only eigenfunctions that lie close in quanta to the eigenfunction corresponding to the most intense spectral line can be determined accurately. Because the band contour, and thus us, is a function of only the short-time dynamics, the thawed Gaussian approximation gives an excellent approximation to us. If, for short times, the thawed Gaussian equations could be solved analytically, then formulas could be derived for us and then inverted to find the initial conditions that minimize us. Alas, this is a cumbersome task, but it is tractable for the case of frozen Gaussians. Let us expand the potential to third order about the initial position
Reimers and Heller
+
+ V,,’”(q - q0)3/6
V = Vo + Vo’(q - qo) + V,,”(q - q 0 ) 2 / 2
(5.4) and solve Hamilton’s equations of motion to this order giving
p , = po - Vir - poV/t2/2
+ (V{V,,”
- po2V,,’”)t3/6
(5.5)
and 41 = qo
+ pot - V,lt2/2 - V
