Nonstationary time series analysis of intramolecular energy transfer

Volume 214. number 2

CHEMICAL PHYSICS LETTERS

29 October I993

Nonstationary time series analysis of intramolecular energy transfer Lisa M. Finney,

Angela Borrmann

and Craig C. Martens

Departmentof Chemistry,Universityof California,Irvine, CA92717-2025.USA Received

I July 1993

The dynamics of intramolecular energy transfer are investigated using an application of nonstationary time series methodology to the analysis of classical trajectory simulations. The approach is based on calculation of a joint time-frequency Wiper-Ville distribution from the evolving particle velocities. The resulting quantity has the interpretation of a time- and frequency-dependent effective temperature; the evolution of the distribution provides a concise view of energy flow between vibrational modes and highlights the importance of nonlinear frequency resonances. The approach is applied to overtone relaxation in a model fouratom chain and the unimolecular decomposition of hydrogen peroxide solvated in an Ar cluster.

1. Introduction The dynamics of intramolecular vibrational energy redistribution (IVR) is a subject of great current interest and research activity. Recent experimental studies have provided detailed information about the states of highly excited molecules (for reviews, see ref. [ 1 ] ), while theoretical investigations have elucidated the dynamical mechanisms of energy transfer in polyatomic systems (for a recent review, see ref. [ 21). Both experiment and theory have highlighted the important role played by low-order nonlinear resonant interactions in enhancing energy transfer in excited polyatomic molecules. Due to anharmonicity, the frequencies of zeroth-order modes depend on their energies, resulting in dynamically evolving resonance conditions as energy flows within the system. Energy transfer is enhanced or inhibited as the frequency ratios of the anharmonic modes “tune” into or out of resonance. A great deal of insight into these processes has been gained by studying classical models of molecular systems from the perspective of nonlinear dynamics (see, for example, ref. [ 3 ] ), which emphasizes the importance of frequency commensurabilities. The corresponding quantum spectral and dynamical consequences of classical nonlinear mechanics have important fundamental relevance 0009-2614/93/%

to understanding the semiclassical limit (see, for example, ref. [4] ). In this Letter, we describe an approach to the theoretical analysis of intramolecular energy transfer in polyatomic molecules. The method we present is based on computing joint time-frequency distributions from the evolving particle velocities generated by classical trajectory simulation. The result is an evolutionary spectrum, which has the interpretation of a time- and frequency-dependent effective temperature. By visualizing the time dependence of the spectrum, a concise representation of energy transfer between “modes”, defined by their intrinsic frequencies, is obtained. Unlike other approaches to the qualitative analysis of nonlinear systems, such as the Poincare surface of ref. [ 3 1, this method can be generalized to systems with many degrees of freedom. A brief report of the method and its application to the many-body dynamics of energy relaxation in van der Waals clusters has been given previously [ 5 1. Here, we apply the approach to the analysis of few-mode molecular vibrations. A number of researchers have previously used spectral analysis methods to characterize IVR dynamics [ 6- 10 1. The rest of this Letter is organized as follows: section 2 describes the methodology employed to analyze classical trajectory data. In section 3, the approach is illustrated by treating energy transfer in a

06.00 0 1993 Elsevier Science Publishers B.V. All ri&ts reserved.

159

Volume2 14,number 2

model four-atom system and the IVR dynamics accompanying the overtone-induced unimolecular decomposition of hydrogen peroxide associated with an Ar cluster. A discussion is given in section 4.

The connection between 8 and more conventional statistical mechanical quantities used to describe equilibrium systems [ 13,141 can be made by suitable averaging. An integration over frequency yields the instantaneous kinetic energy of the system, 1 O” @((,w) do= j$, fmj$(t)=K(t) YG I -m

2. Method As dictated by the uncertainty principle, it is impossible to have complete simultaneous time and frequency resolution of a signal. However, it is possible to define joint time-frequency distributions which balance the unavoidable uncertainty in a useful way. One example is the Wigner-Ville distribution [ 111 (see, for example, ref. [ 121). This quantity is well-known in quantum mechanics [ 11,131, where, with the replacement of time and frequency by position and momentum, a phase space-like distribution known as the Wigner function is obtained. Wigner-Ville distributions have also been employed in a number of other applications in the field of signal processing [ 121. Our method for nonstationary spectral analysis of molecular dynamics simulation data is based on the construction of a joint time-frequency distribution function from the evolving dynamical variables. A Wigner-Ville distribution is calculated from the Cartesian velocities of the particles in the system. We define the function e( t, w) by

where m Vj(t’tT)‘V,(t+tT)

ej(t9m)=tm, I

-uJ

(2)

In eq. (2), mj and v, are the mass and velocity of atom i, respectively. The sum in eq. (1) is over the number of particles N. @(t, w) is thus the sum of single particle contributions. The units of 8 are [energy] X [time], or equivalently, energy per unit frequency, and the quantity is interpreted as a dynarnically evolving spectral decomposition of the instantaneous kinetic energy of the system - a timeand frequency-dependent effective “temperature”. 160

,

(3)

where the integral representation of the &function, 1 Oz eiU7dw=S( r) , zn 5 -Q)

(4)

has been used. An integration over time yields a spectral density, which is defined as the Fourier transform of the mass-averaged velocity autocotrelation function,

xeiwrdr ,

(5)

where the angular brackets ( > denote time averaging. Finally, integrating over both frequency and time yields the time-averaged kinetic energy, which is proportional to the temperature T of the system, T m @(t,m)dtdw=t

$ mj(r$) j-1

=(K)

xeiWrdr.

29 October 1993

CHEMICALPHYSICS LETTERS

=;NkT,

(6)

where k is Boltzmann’s constant. Eqs. ( 3 )- ( 6 ) establish the connection between the distribution @(t, w) and statistical mechanical or thermodynamic properties of equilibrium systems by integrating over one or both of the independent variables (t, w ). For treating energy transfer and relaxation in nonequilibrium systems, however, the integrations over time or frequency need not be performed, allowing the time evolution of the partitioning of energy as a function of frequency to be monitored. As we shall illustrate in the next section, this provides a concise representation of the energy relaxation pathways and dynamical mechanisms underlying IVR in polyatomic molecules. In numerical applications, the Wigner function itself has two undesirable properties: e(t, w) is not

Volume2 14, number2

CHEMICALPHYSICSLETTERS

positive definite, complicating its interpretation as a temperature, and can experience rapid oscillations as a function of time or frequency. To avoid these complications, we calculate a local Gaussian smoothing of the Wigner distribution. This function is always positive, free of oscillations, and can be computed numerically in terms of Gaussian windowed Fourier transforms. The resulting function @n( t, co) is a Husimi distribution [ IS], or equivalently, a coherent state representation [ 161, of the time series. It is given by

-m

-a3

xexp[ - (s-t)2/a2]

exp[ -02(a-o)2]

,

Table 1 Parameters for four-atom chain Hamiltonian, in dimensionless units ml m2 m3 m, DL2 4.3 D 3.4

81.1 Bz.3 83.4 Xl,2 XT.3 x5.4

0.7275 2.9100 2.9100 1.4550 14.110 7.901 14.11 0.1384 0.1364 0.1384 18.24 27.44 18.24

(7)

where Q is a width parameter determining the relative smoothing windows in frequency and time. The results presented below are obtained using the smoothed Husimi distribution.

3. Results 3. I. Intramolecular energy redistribution in a model four-atom chain We first illustrate the method by considering the dynamics of intramolecular energy transfer in a model system, consisting of a chain of four atoms moving in one dimension and bound by pairwise Morse potentials. The Hamiltonian is given in Cartesian coordinates by

X{l-exP[-B,j+l(X,+I-X,-~,+1)1}*.

29 October1993

(8)

This system has three vibrational degrees of freedom, consisting of the three anharmonic local bond modes. Although there is no potential coupling in this model, the bond modes can exchange energy due to kinetic coupling. The particle masses and potential parameters, and the resulting harmonic normal mode frequencies, are given in table 1; the normal modes of the system correspond closely to the local bond modes and are identified as such in the table. This

system is chosen as a model problem, and all quantities are in dimensionless units. In fig. 1, we show data for a trajectory of the system, integrated using a sixth-order hybrid Gear routine [ 171, with initial conditions corresponding to a large excitation of mode 1, with the other degrees of freedom initially cold. The energies of modes 1, 2, and 3 are displayed in figs. la-lc, respectively. The initial excitation is trapped in mode 1 until time t= tl, when energy is suddenly transferred to the adjacent bond, mode 2. Somewhat later, at time t2, energy escapes from mode 2 and flows into the final bond, mode 3. This particular trajectory thus exhibits a sequential flow of energy down the bonds of the chain. In fig. 2, we show the results of applying time-frequency analysis to the trajectory illustrated in the previous figure. The figure gives a family of displaced frequency spectra, corresponding to constant time slices of the distribution function &( t, co), with the vertical displacement of the base line indicating the time of the spectrum. The absolute amplitude of the distribution is arbitrary; the relative amplitudes at different times are accurately represented. Each frequency spectrum shows three major peaks, corresponding to the three vibrational frequencies of the system. The identity of each major peak is indicated at the top of the figure by labels w,(t) Cj= 1, 2, 3). The peak positions and amplitudes evolve with time with the dynamics of the trajectory. At early times (near the bottom of the figure), the time-frequency distribution is characterized by one dominant peak around w= 0.84. This corresponds to 161

Volume214.number2

CHEMICAL PHYSICSLETTERS

29 October1993

5000

4000

.g

3000

b o

2000

F.,,,,,,,,,,,,,,,,.,,..,.i

0

1000

2000

3000

4000

5000 1000

time

0

5

5

b)

t2

0

0.2

0.4

0.6

0.6

1

12

Frequency Fig. 2. Time-frequency distribution &.,(t, w) for the trajectory shown in fa. 1. See text for discussion.

0’ ,“l,,,,“,,“.,‘,‘,‘.” 0

1000

2000

3000

4000

5000

4000

5000

time

4 ti 2

3

w2 1 0 0

1000

2000

3000

time Fig. 1.Localmodeenergiesversustime for a singletrajectoryof the modelfour-atomchainsystem.(a) Mode1. (b) Mode2. (c) Mode3.

the large amplitude initial excitation of mode 1. The peak frequency is shifted by anharmonicity from the harmonic value of 07 = 0.977. Much smaller peaks are visible at approximately 0=0.34 and w=O.74. These are close to the harmonic frequencies &=0.366 and wi=O.774 of modes 2 and 3, respectively (see table 1). The spectrum remains quite stationary until t= t,. At that point, which is marked on fig. 2, the amplitude of the mode 1 peak drops 162

sharply, with an accompanying increase in its frequency. This is indicative of sudden energy transfer out of mode 1; the increase in frequency is due to the anharmonicity of the mode. Simultaneously, the amplitude of the peak corresponding to mode 2 increases, indicating energy flow into the center bond of the chain. The value of o2 (t) decreases, again due to the anharmonicity of the mode. During these processes, the amplitude corresponding to mode 3 remains small. At time tz, a second sudden energy transfer event occurs. Here, the energy exchange is out of mode 2 and into mode 3. This is reflected in fig. 2 by a drop in the mode 2 amplitude together with an increase in the frequency and a growth of the amplitude of the mode 3 peak. The qualitative view of mode-mode energy transfer provided by time-frequency analysis is consistent with the local mode energies versus time shown in fig. 1. In addition, the joint time-frequency distribution of fig. 2 gives a view of the underlying mechanism of energy flow. The first energy transfer event occurring at time t, is due to a 3 : 1 resonance between the anharmonic modes 1 and 2. This manifests itself in the Husimi distribution by a near coincidence of the third mechanical overtone of the receiving mode 2 (labeled R on fig. 2) and the fundamental of the excited bond 1. The resonance condition is satisfied during the energy transfer, and as energy flows into mode 2, the overtone Q= 3~~ ( t )

Volume214,number2

CHEMICALPHYSICSLETTERS

emerges from under the fundamental w, (1). A 3: 1 commensurability is also responsible to the mode 2+mode 3 energy flow occurring at time tZ, which corresponds to a proximity of a and the fundamental us(t) at the point of energy transfer. 3.2. Overtone-inducedunimoleculardecomposition of H20-Ar13 As a second illustration of the method, we consider the intramolecular dynamics and unimolecular decomposition of H202 following overtone excitation. The example shown corresponds to the molecule solvated in an Arl 3 cluster, which we have studied previously as a model of the effect of microsolvation on chemical reaction dynamics [ 18 1. The properties of this system have been described in detail in ref. [ 18 1, and we refer the reader there for a discussion of the potential energy surface employed and the method used to generate ensembles of initial conditions. Here, we apply our method of nonstationary time series analysis to gain a complementary perspective of the IVR dynamics. In fig. 3, the dynamics of a typical dissociative trajectory from an ensemble with initial OH excitation corresponding to the v=9 state are presented. The remaining intramolecular and cluster degrees of freedom are characterized by an initial temperature of T= 10 K. In the figure, we show the time dependence of the internal vibrational degrees of freedom, calculated as described in ref. [ 18 1. Figs. 3a-3f give the energies of the initially excited OH stretch, the 00 stretch, the initially unexcited OH stretch, HO0 bend 1, OOH bend 2, and the torsion, respectively. The duration of the trajectory integration is 12 ps. The initial OH stretch energy is over 25000 cm-‘. This excitation remains trapped for approximately 3.5 ps, at which time energy transfers into the bend modes and the 00 stretch. Vibrational redistribution then continues until unimolecular dissociation occurs. The distant OH stretch and the torsional mode do not participate strongly in the IVR dynamics. The decomposition of Hz02 occurs at around t = 8.5 ps, and is reflected in fig. 3 by a sharp increase in the 00 energy to above the dissociation threshold D,= 17340 cm-l. Energy flow between the OH bonds ceases as energetically isolated fragments are created. In fig. 4, the time-frequency distribution &( t, CO)

29 October1993

calculated from the velocity time series for this trajectory is shown. The major features of the distribution have been labeled in the figure. At early times, the energy of the system is localized in the excited OH stretch. The predominant features in the spectrum for the first 3 ps are the large fundamental frequency peak at around 1700 cm-’ and a series of mechanical overtones at integer multiples of the fundamental. The frequency shift from the harmonic limit of ~~3880 cm-’ and the pronounced overtones are due to the very strong anharmonicity of the OH oscillator at v= 9. (The harmonic frequencies of this model of hydrogen peroxide [ 181 are given in table 2). The spectrum undergoes a marked change at an elapsed time of approximately 3.5 ps, with the appearance of a number of new features; these peaks have also been identified on fig. 4. As indicated by the mode energies shown in fig. 3, a sudden transfer of excitation from the OH stretch into the bend modes and the 00 stretch occurs at this time. This is reflected in the appearance of peaks in the timefrequency distribution corresponding to the bend modes, which are located at around W= 1300 cm-‘, and 00 stretch, which fluctuates between approximately 500 and 800 cm-‘. This is accompanied by a simultaneous reduction in the OH stretch amplitude and a shift of the frequency upward to around 3000 cm-‘. The deviations of the frequencies from the corresponding normal frequencies (see table 2) are due to anharmonicity. At t = 8.5 ps, unimolecular decomposition occurs. This event is reflected in the time-frequency distribution by an increase in the zero frequency component, corresponding to translational (i.e. zero frequency) motion accompanying OH and Ar atom dissociation, and liquid-like motion of the Ar cluster remnant [ 191, together with a disappearance of the 00 stretch and bend features in the spectrum. In addition, energy transfer involving the OH stretches ceases, leaving two isolated and relatively cold hydroxyl radicals with spectral features around 3500 cm-‘, which are modulated by the rotational motion of the fragments.

163

Volume 2 14,number 2 30000.0

.' "

-

2500D.o

I 7

%Y

20000.0

29 October1993

CHEMICALPHYSICSLETTERS I

I.

I "

1

'. ' a>

30000.0

'+

_ ':

25000.0

ti 20000.0 v

~15000.0

:

!f E 10000.0

:

3

15000.0

r; 10000.0 w

W 5000.0

:

5000.0

,

r

o.ot.“““l.““““““.” 0

2

0.0

4

6

6

10

0

12

2

4

6

6

10

12

10

12

10

12

time (psec)

time (psec)

_ v

25000.0

E

20000.0

z G

15000.0

5 c; 10000.0 W 5000.0 0.0 0

2

4

6

6

10

0

12

2

time (psec)

25000.0

E Y

20000.0

6

1000.0

h t? 15000.0 2

6

time (psec)

30000.0

_

4

10000.0

-

600.0

?g

600.0

i

400.0

2

W

W 5000.0 0.0

200.0

00 0

2

4

6

'8

10

12

time (psec)

0

2

4

6

6

time (psec)

Fig. 3. Time dependence of mode energies in internal coordinates for a single trajectory of the H20Z-AI13system. (a) Initially excited OH stretch. (b) 00 stretch. (c) Initially unexcited OH stretch. (d) Bend 1. (e) Bend 2. (f ) Torsion.

4. Discussion In this Letter, we have presented a computational tool for analyzing classical trajectory simulations of intramolecular dynamics, based on an application of nonstationary time series analysis. A joint time-frequency distribution is defined by computing Gaussian-smoothed Wigner-Ville distributions of the 164

evolving particle velocities, providing a dynamic spectral decomposition of the instantaneous kinetic energy of the system. The result is a time- and fre quency-dependent effective temperature, and the evolution of the distribution gives a concise visual representation of energy flow within the molecule. The results of this analysis were compared with conventional model energy data, and the view provided

CHEMICALPHYSICSLETTERS

Volume2 14,number2

29 October1993

Foundation and the Of&x of Naval Research. We also acknowledge the UC1 Ofi& of Academic Computing for an allocation of computer time.

References [ I ] S. Leone, Ann. Rev.Phys.Chem.35 (1984) 109: F.F.Crim,Ann. Rev.Phys.Chem.35 (1984) 657; C.E.Hamilton,J.L.Kinscyand R.W.Field,Ann.Rev.Phys.

D

1000

2000

3000

4000

5000

6000

,000

Frequency (cm.‘) Fig. 4. Time-frequencydistribution0&t, o) for the trajectory of H,02-Arl,! shown in tig. 3. Seetextfor discussion. Table2 Hydrogen peroxide normal mode frequencies Mode

0 (cm-‘)

OH symm. str.

OH asym.str. OOHsymm.bend

3883 3882 1474

OOH asym. bend 00 str. torsion

1402 937 406

by the time-frequency distribution in terms of timedependent peak amplitudes and positions is consistent with the evolving mode energies in the two examples considered. An advantage of the timefrequency method is that a set of zeroth-order modes need not be defined in order to monitor energy transfer. From the Cartesian velocity coordinates, a unique spectral decomposition of the kinetic energy is performed, in effect determining the energy partitioning over modes characterized by their instantaneous frequencies. This approach is advantageous in situations where a good zeroth-order representation is not known a priori.

Acknowledgement This work was supported by the National Science

Chem. 37 (1986) 493; M. Quack,Ann. Rev. Phys. Chem. 41 (1990) 839. [2] T. Uzer and W.H. Miller, Phys. Rept. 199 (1991) 73. [3] V.I.Arnold,Mathematicalmethodsof classicalmechanics

(Springer,Berlin,1978); A.J. Lichtenbergand M.A. Lieberman, Regular and chaotic dynamics, 2nd Ed. (Springer, Berlin, 1992); S.N. Raaband, Chaotic dynamics of nonlinear systems (Wiley,New York, 1990). [4] A.M. Ozorio de Almeida, Hamiltonian systems: chaos and quantization (Cambridge Univ. Press, Cambridge, 1988); M.C. Gutzwiller,Chaos in classicaland quantum mechanica (Springer, Berlin, 1990); MS. Child, Semiclassical mechanics with molecular applications (Clarendon Press, Oxford, 1991). [ 51C.C. Martens, Phys. Rev. A 45 (1992) 6914. [61 J.D. McDonaldand R.A.Marcus,J. Chem. Phys. 65 ( 1976) 2180. [7] C.C. Martens, M.J. Davis and G.S. Ezra, Chem. Phys. Letters 142 (1987) 519. [ 81 R. Roy, B.G. Sumpter, GA. Pfeffer, SK. Gray and D.W. Noid, Phys. Rept. 205 ( 1991) 11I I (9lT.D. Sewell, D.L. Thompson and R.D. Levine, J. Phys. Chem. 96 (1992) 8006; X.Y. Chang, D.L. Thompson and L.M. RatT,Chem. Phys. Letters 206 (1993) 137. [ lo] D. Kleinhesselink and M. Wolfsberg,Chem. Phys. Letters 205 (1993) 461. [ 111E.P. Wigner, Phys. Rev. 40 (1932) 749. [ 121M.B.Priestley, Spectral analysisand time series (Academic Press, New York, 198I ); J.M.Comb, A. Grossmannand Ph. Tchamitchi, e&., Wavelets(Springer, Berlin, 1989). [I31 L.E. Reich& A modem course in statistical physics (University of TexasPress, Austin, 1980). [ 141D.A.McQuarrie, Statistical mechanics (Harper and Row, New York, 1976). [ 151K. Husimi, Proc. Phys. Math. Sot. Japan 22 (1940) 264. [ 161J.R.IUauderand B. Skagerstam,eds., coherentstates (World Scientific,Singapore, 1985). [ I7 I C.W. Gear, SIAMJ. Num. Anal. 2B ( 1964) 69. [ 181L.M. Finney and CC. Martens, J. Phys. Chem. 96 (1992) 10626. [ 191R.S. Berry, T.L. Beckand H.L. Davis, Advan. Chem. Phys. 70 (1988) 75.

165