Iterative path integral calculation of quantum ... - Semantic Scholar

Chemical Physics 268 (2001) 1±10

www.elsevier.nl/locate/chemphys

Iterative path integral calculation of quantum correlation functions for dissipative systems Jiushu Shao *, Nancy Makri Departments of Chemistry and Physics, University of Illinois, 601 S. Goodwin Avenue, Urbana, IL 61801, USA Received 13 December 2000

Abstract An iterative path integral procedure is developed for calculating equilibrium two-time correlation functions of quantum dissipative systems. Its numerical feasibility is demonstrated by studying the position±position correlation function of two-level systems at various temperatures. It is also shown that canonical partition functions are also amenable to an iterative treatment. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction Understanding the quantum dynamics of complex systems is essential to elucidate many processes such as electron or proton transfer, vibrational relaxation, and electron transport which are ubiquitous in biology, chemistry and physics. In the so-called dissipative systems, the interesting physical quantities that can directly be measured are observables pertaining to a system consisting of only a few degrees of freedom. The surrounding medium or bath, whose coupling to the system results in a behavior generally referred to as dissipation, is not probed directly. Thus, the relevant dynamics is characterized by the reduced density matrix, which is de®ned as the trace of the density matrix for the total system with respect to the degrees of the freedom of the bath. In FeynmanÕs path integral treatment of the dynamics [1,2], all the characteristics of the bath that are

*

Corresponding author.

relevant to the evolution of the systemÕs reduced density matrix are captured in the in¯uence functional [3]. A successful paradigm of dissipative dynamics is the Caldeira±Leggett model [4] in which the bath is mimicked by an in®nite number of harmonic oscillators. Using the path integral approach, Caldeira and Leggett showed that the collective e€ects of the bath on the system of interest can be described in terms of a spectral density. When the system±bath interaction is linear in the bath coordinates, the in¯uence functional assumes a Gaussian form that involves integrals of the spectral density and introduces nonlocal e€ects to the path integral expression for the reduced density matrix of the system. This nonlocality de®es a straightforward use of iterative algorithms that is routinely used for following the quantum dynamics of non-dissipative systems with few degrees of freedom. Several approximate theories have been developed for calculating the dynamics in harmonic dissipative environments with di€erent systems and in diverse parameter regimes. Among the most common treatments are the traditional resolvent

0301-0104/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 1 ) 0 0 2 8 6 - 5

2

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

operator method [5] and Red®eld theory [6]. The non-interacting blip approximation [7], a powerful tool developed in the context of the dissipative two-level system, revealed a host of interesting behaviors ranging from exponential decay to damped oscillatory evolution and even complete localization. Numerically exact evaluation of the path integral entails performing a summation over all paths of the probed system. The number of paths grows exponentially with the number of time steps, and thus direct summation of the path integral [8,9] is practical only at short times. Monte Carlo sampling is hampered by the sign problem and therefore useful only for short-time calculations under strong damping conditions [10±15]. Accurate long-time path integral calculations have become possible with the aid of an iterative procedure developed earlier in our group [16±19]. This is based on the observation that the nonlocality of the e€ective action spans a ®nite time domain, just as in the case of memory in the classical generalized Langevin equation. The algorithm developed in the context of the reduced density matrix relies on the assumption of ``factorized initial conditions'' [7], i.e., the system and environment are uncorrelated at the start of the propagation. Within this framework, it is possible to obtain the exact long-time dynamics of few-level systems interacting with harmonic environments of moderate memory span [20]. When sluggish solvents are involved, the medium-induced memory can become long, necessitating the use of a Monte Carlo selection process to identify path segments of appreciable weight [21] or the use of memory equation algorithm [22]. Finally, it is worth noting that an iterative decomposition of the path integral remains well-posed even in cases of anharmonic environments [23]. Spectroscopy provides the convenient experimental tool for studying quantum dynamics in dissipative media. In general, spectroscopic observables can be cast in terms of Fourier integrals of various time correlation function. In the case of linear spectroscopy, the relevant correlation functions in the canonical ensemble take the form GAB …t† ˆ

1  Tr e Z

bH

AeiHt=h Be

iHt= h



;

…1:1†

where H is the Hamiltonian of the total system, 1 b ˆ …kB T † the reciprocal temperature T scaled by the Boltzmann constant kB , Z the partition function at the given temperature, and A, B the operators of the observed system. Time correlation functions are similar in structure to the reduced density matrix. However, correlation functions measure the dynamical ¯uctuations of the system close to equilibrium, implying that the initial density must involve the Boltzmann operator for the fully coupled Hamiltonian. Since the matrix elements of the latter are not available in closed form, one must evaluate them via a separate imaginary time path integral calculation [24]. Another possibility is to focus on the symmetrized correlation function CAB …t† ˆ

1  iH s =h Tr Ae Be Z

iH s= h



;

…1:2†

where the complex time s ˆ t ihb=2 is introduced. The two correlation functions contain the same dynamical information because their Foue AB …x† ˆ rier transforms satisfy the relation G hbx=2 e  e C AB …x†. Notice that the symmetrized correlation function CAB …t† is purely real. It is also useful to note that ¯ux correlation functions [25], which provide direct information about reaction rate constants, are expressed in the symmetrized form of Eq. (1.2). We envision the calculation of reaction rates in dissipative environments as one of the principal applications of the methodology presented in this paper. Section 2 reviews the quasi-adiabatic propagator path integral representation of the symmetrized correlation function and investigates the behavior of the in¯uence functional. We show that the nonlocal interactions of the latter decay along all directions in the complex plane if the environment has a broad spectrum. This decay is a consequence of dynamical ¯uctuations exploited in earlier work as well as thermal ¯uctuations that are unique to expressions involving propagators in complex time. Based on these features we present an iterative algorithm for evaluating correlation functions in harmonic dissipative environments. Section 3 illustrates the methodology with numerical examples on a dissipative two-level system

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

Hamiltonian, and Section 4 summarizes with a brief discussion and some concluding remarks. 2. Theory and methodology 2.1. Symmetrized correlation function and path integral representation Throughout this article we focus on the dynamics of quantum dissipative systems described by the Hamiltonian 8 !2 9 = X< pj2 1 c f …s† j H ˆ H0 …s; ps † ‡ ‡ mj x2j xj ; : 2mj 2 ; mj x2j j …2:1† where s denotes the coordinate(s) of the observable system and xj are the coordinates harmonic phonon-like degrees of freedom that constitute the dissipative bath. Note that the interaction of the system with the bath is linear in xj and depends on the system position through an arbitrary function f …s†. The functional form of f …s† has implications on the spectrum [26]. Since A and B are operators that depend only on the coordinate and/or the momentum of the system, the symmetrized correlation function CAB …t† can be written as Z Z Z ds ds‡ CAB …t† ˆ Z 1 ds‡ 0 0 N Z

 ‡   sN B sN K…s  dsN s‡ 0 A s0 0 ; sN †; …2:2† where

iH s =h  ‡  s s e K…s 0 ; sN † ˆ Trbath s0 e N N



iH s= h ‡ s0

…2:3† is the trace of a product of complex time evolution operators over the bath and the notation s N is introduced in anticipation of the discretized path integral representation for the time evolution operators in Eq. (2.3). Note that the kernel K is the same as that in Ref. [9]. Now the propagators are represented in terms of discretized path integrals in which the complex time s is divided into N slices.

3

Earlier work has shown that the quasi-adiabatic propagator splitting [27], in which the short time propagators are split symmetrically into a reference part pertaining to the one-dimensional Hamiltonian H0 along the adiabatic path xj ˆ cj f …s†= mj x2j and a part comprising displaced bath oscillators, is advantageous as it is suciently accurate with moderately large time steps. The quasi-adiabatic propagator path integral representation of the kernel takes the form [28] Z Z N Y 

    K…s ; s † ˆ ds


ds sk 1 eiH0 Ds =h sk 0 N 1 N 1

 s‡ k e



kˆ1

iH0 Ds= h ‡ sk 1



  F …s 0 ; s1 ; . . . ; sN †;

…2:4† where Ds ˆ s=N and F is the quasi-adiabatic propagator version of Feynman and VernonÕs in¯uence functional [3]. The latter factorizes by virtue of the separable form of the bath Hamiltonian. It can be shown that the in¯uence functional in Eq. (2.4) is given by the expression
2 ! Z Y c2 f s…s0 † iX 0 0 j F ˆ Zj exp ds h j C 2mj x2j j Z Z X 1  exp ds0 ds00 4hmj xj C C j   cos xj …sa sb ‡ ihb=2†  sinh …hbxj =2† !  c2j f …s…s0 ††f …s…s00 †† ;

…2:5†

where Zj0 is the canonical partition function of the jth harmonic oscillator, and sa , sb are the earlier and later of the times s0 , s00 along the time contour C displayed in Fig. 1. The total in¯uence functional thus reads Z

X c2j f s…s0 † ds 2mj x2j C j j  Z Z 1  exp ds0 ds00 a…sa h C C 

sb †f s…s0 † f s…s00 † ;

Y F ˆ Zj0 exp

i h


2 !

0

…2:6†

4

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

Fig. 1. The contour of integration in the complex time plane.

where a…s† is the bath response to the feedback of the system, given by the expression   X c2j cos xj …s ‡ i hb=2† : …2:7† a…s† ˆ 2mj xj sinh …hbxj =2† j It is well known that all characteristics of the bath pertaining to the dynamics of the observable system are speci®ed by the spectral density function J …x† ˆ

p X c2j d…x 2 j mj xj

xj †:

…2:8†

2.2. Response function and memory in complex time The spectral density characterizing dissipative environments is a broad function. Under these

conditions, and when the canonical density operator is replaced by a ``factorized'' initial condition (in which case the integration contour in Fig. 1 runs along the real time axis) the response function decays irreversibly to zero, implying that the nonlocality of the in¯uence functional spans a ®nite range. This observation provides the basis of an iterative calculation of the dynamics as detailed in Refs. [17,19]. It is easy to see that this behavior persists in the present case, where the integration contour is complex. In fact, there is now an additional cause for decay of the response function: thermal ¯uctuations, namely the imaginary part of the argument in Eq. (2.7), lead to an exponential decay of a…s† away from the vertices of the triangle de®ned by the integration contour. Since hb 6 Im s 6 0, the response function attains its largest values when the imaginary part of the complex time t is near one of the two limits, 0 and hb. This behavior is illustrated very clearly in Figs. 2 and 3. As shown in Fig. 2, even when the bath consists of a single harmonic oscillator and when the real time equals zero, the response function decays along the imaginary axis. This ``imaginary time decoherence'' is entirely due to thermal ¯uctuations, and the response function is entirely nonlocal along the real time axis. If many degrees of freedom are added to the bath, leading to a broad spectral density, the decay of the response function becomes faster and now extends to all directions away from vertices in the complex time plane. These e€ects are illustrated in Fig. 3, which shows

Fig. 2. The response function a…s† in the case where the bath consists of a single oscillator of frequency x such that hbx ˆ 10: (a) real part; (b) imaginary part.

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

5

Fig. 3. The response function a…s† for an ohmic bath with hbxc ˆ 7:5: (a) real part; (b) imaginary part.

the response function in the case of an Ohmic bath corresponding to the spectral density J …x† ˆ pnxe

x=xc

:

…2:9†

It is seen that the nonlocality of the response function and therefore the in¯uence functional is ®nite along any direction in the complex plane. These features allow evaluation of the correlation function via an iterative procedure which is described below. 2.3. Iterative propagation of the dynamical factor With Eq. (2.7) the discretized in¯uence functional in Eq. (2.6) becomes  N Y 1X 0 F ˆ Zj exp a‡‡ f ‡2 ‡ ak fk 2 h kˆ0 k k j ‡ ‡ a‡ k fk fk

ˆ a‡ k




‡ bkk0 fk fk0 ‡

N X k  1X b‡‡0 f ‡ f ‡0 h kˆ0 k0 ˆ0 kk k k 

b‡ kk 0

fk‡ fk0

‡

bkk‡0 fk

fk‡0



 sin …x…sk‡1 ‡ sk

;

where fk ˆ f …s k † and the coecients are given by the expressions 

…2:11†

sk‡1

sk

ihb†=2†; …2:12†



b‡‡ kk 0 ˆ …bkk 0 † Z 2 1 J …x† sin …x…sk‡1 sk †=2† ˆ dx 2 p 1 x sinh …hbx=2†
 sin x…sk0 ‡1 sk0 †=2  cosx…sk‡1 ‡ sk sk0 ‡1 sk0 ‡ ihb†=2†; …2:13† and 

‡ b‡ kk 0 ˆ …bkk 0 † Z 1 2 J …x† sin …x…sk‡1 sk †=2† ˆ dx 2 p 1 x sinh …hbx=2†

 sin …x…sk0

)

…2:10†

a‡‡ ˆ …ak † k Z 1 1 J …x† sin …x…sk‡1 sk †=2† ˆ dx 2 p 1 x sinh … hbx=2†  sin …x…sk‡1 sk ‡ i hb†=2†;

Z 2 1 J …x† sin …x…sk‡1 sk †=2† dx 2 p 1 x sinh …hbx=2†
   sin x…sk sk‡1 †=2



cos…x…sk0 ‡1

sk0 ‡1 †=2† ‡ sk0

sk‡1

sk

ihb†=2†; …2:14†

where the spectral density for negative frequencies is de®ned as J … x† ˆ J …x†. Except for the ®rst and last time points, the nonlocal terms in Eq. (2.10) within the forward or the backward branch of the time contour depend on the ``distance'' …k k 0 † Ds of the integration variables. Such a property no longer holds for the interactions of points in the forward segment with coordinates in the backward branch of the time contour. The

6

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

decay of the response function implies that the coecients of the nonlocal terms become practically zero when jk k 0 j exceeds a certain value Dkmax , and thus the second summation in the exponent of Eq. (2.10) can be truncated without loss of accuracy. The ``memory length'' Dkmax jDsj is treated as a convergence parameter in the calculation. To implement the iteration procedure in the calculation of the real time dynamics of symmetrized correlation function (4), we recast the in¯uence functional (12) as a product of functions corresponding to one-point and two-point interactions, namely, F ˆ

Y j



Zj0

N Y

N Y1

kˆ0

kˆ0

N Y Dkmax kˆ0

F0 …s k †

 F1 …s k ; sk‡1 †




 FDkmax …s k ; sk‡Dkmax †;

…2:15†

where the one-point and two-point functions are  
1 ‡‡ ‡2  2 ‡ ‡ a f ‡ ak fk ‡ ak fk fk F0 …sk † ˆ exp h k k  …2:16† and   FDk …s k ; sk‡Dk † ˆ exp

  ˆ K0 …s k ; sk‡1 †F0 …sk †

minfDkY max ;N kg mˆ1

 Fm …s k ; sk‡m †;

…2:19† where

‡  K0 …s k ; sk‡1 † ˆ sk‡1 e

 iH Ds =h s s e 0

iH0 Ds= h ‡ sk

k

 k‡1

…2:20† is the forward±backward propagator of the bare system. The operation  R…s k‡1 ;


; sk‡Dkmax ; …k ‡ 1†Dt† Z   ˆ ds k K…sk ; . . . ; sk‡Dkmax †   R…s k ; . . . ; sk

1‡Dkmax ; k Dt†

…2:21†

propagates the array R forward by one time step. It is straightforward to check that the symmetrized correlation function at the time t ˆ N Dt is given by CAB …t† ˆ Z 

1

Z

‡ 
ds N sN B sN R sN ; 0; . . . ; 0; N Dt : …2:22†

1 h

‡ ‡ b‡‡ k‡Dk;k fk‡Dk fk

2.4. Iterative evaluation of the quantum partition function

‡ bk‡Dk;k fk‡Dk fk ‡ ‡ b‡ k‡Dk;k fk‡Dk fk

‡

  K…s k ; sk‡1 ; . . . ; sk‡Dkmax ; kDt†

‡ bk‡Dk;k fk‡Dk fk‡



:

…2:17†

With the above de®nitions, the calculation of the symmetrized correlation function proceeds as described below. We de®ne a Dkmax -dimensional ar ray R…s with the initial k ; . . . ; sk 1‡Dkmax ; kDt† condition

‡   R…s …2:18† 0 ; . . . ; sDkmax 1 ; 0† ˆ s0 A s0 ; whose elements are system path segments that span the memory length Dkmax Dt in real time. We also de®ne a …Dkmax ‡ 1†-dimensional array K, which is de®ned as

Apart from the dynamical part considered above, calculation of the symmetrized correlation function also requires knowledge of the quantum partition function for the full Hamiltonian. Canonical partition function can be routinely evaluated by means of imaginary time path integral methods which employ Monte Carlo sampling to perform the multidimensional integration. One of the goals of the present paper is to point out that the partition function itself is also amenable to iterative evaluation. This is a consequence of the decay properties that the bath response function exhibits along the imaginary time axis, as described above and illustrated in Figs. 2 and 3. To this end, notice that the partition function can be obtained from the dynamical factor in the

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

symmetrized correlation function by setting t ˆ 0 and A ˆ B ˆ 1. Thus, its evaluation follows precisely the procedure developed in the previous subsection for real time equal to zero, in which case the contour is a straight line along the imaginary axis up to the point ihb. However, according to the analysis presented earlier, the response function (and thus the nonlocality of the in¯uence functional) may, depending on the temperature and the nature of the bath spectrum, decay faster along other directions in the complex time plane. The change of memory length with time is displayed in Fig. 4 where the number of evolution steps in the complex plane is ®xed for all times. It is seen that the real part of the response function drops o€ more rapidly at longer times. Provided one does not cross any singularities, one is allowed to deform the contour of integration in any desired fashion. If the dynamical ¯uctuations of the bath are more e€ective than thermal ¯uctuations in quenching the memory, it is advantageous to modify the integration contour to include a real time component. Thus, one may

7

Fig. 4. The real part of the response function a…t; n† for a given number of complex time slices, N ˆ 80, at di€erent times with the parameters of Fig. 3.

calculate the quantum partition function via an iterative procedure along time contours of the type displayed in Fig. 1. The results of such a calculation are shown in Fig. 5, which displays the partition function for a symmetric TLS coupled to an ohmic bath with n ˆ 0:1 and hbX ˆ 0:2 (see Section 3) as a function of Dkmax with a time step equal to 0:1=X. Clearly, faster convergence for a

Fig. 5. Comparison of the results of calculating partition function at di€erent times. The solid and hollow circles show the results of calculation with t ˆ 0 and t ˆ 2X 1 , respectively.

8

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

given complex time step is achieved when the contour of integration goes o€ the imaginary axis. The possibility to calculate the partition function iteratively may come as a surprise, since correlations along the imaginary time axis between the times 0 and ihb are indeed very important. The important feature that may have escaped attention in earlier works is the decay of these correlations halfway between these two time points, toward the imaginary time i hb=2. Thus, by pairing coordinates corresponding to time points located symmetrically with respect to i hb=2 and proceeding inward toward the midpoint we have been able to devise an iterative algorithm. In practice, we have observed that the highest rate of convergence is achieved when we evaluate the partition function at each time point along the same complex time contour involved in the dynamical part and with the same complex time step and memory length. This simple rule is a consequence of error cancellation.

3. Application to dissipative two-level systems In this section we apply the iterative procedure developed in this paper to calculate the position correlation function of symmetric and asymmetric two-level systems coupled to a generic dissipative bath of Ohmic spectral density. The Hamiltonian

is given by Eq. (2.1) where the system coordinate is replaced by the 2  2 Pauli spin matrix rz , i.e., ( ) X pj2 1 2 2 ‡ mj xj xj rz cj xj ; H ˆ hXrx ‡  herz ‡ 2mj 2 j …3:1†

where 2hX is the tunneling splitting and the bath spectral density is given by Eq. (2.9) with xc =X ˆ 7:5. Below we report path integral results for the symmetrized correlation function
C…t† ˆ Z 1 Tr rz eiHt=h bH =2 rz e iHt=h bH =2 : …3:2† Fig. 6 displays C…t† and the original correlation function G…t† at the high temperature b 1 ˆ 5hX for the symmetric case e ˆ 0 and n ˆ 0:1. Fig. 7 shows analogous results at a lower temperature, b 1 ˆ hX. In both cases, the absolute time step, namely, jt ihb=2j=N is chosen as 0:1X 1 . The convergence of the results is excellent with Dkmax ˆ 6. Because of the simple relationship e e G…x† ˆ ehbx=2 C…x† between the two correlation functions, G…t† is feasibly obtained from C…t† through Fourier transforms. As one expects, at high temperatures the symmetrized correlation function C…t† is similar to the real part of the original one G…t†, which is indeed the case for b 1 ˆ 5hX as shown in Fig. 6. Note also that at this temperature quantum coherence, which is characterized by the imaginary part of

Fig. 6. Correlation functions at b 1 ˆ 5hX. (a) Symmetrized correlation function. Solid line and markers correspond to Dkmax ˆ 6 and 9, respectively. (b) Original correlation function. Solid and dashed lines show real and imaginary parts calculated from the data for C…t† with Dkmax ˆ 6.

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

(a)

9

(b)

Fig. 7. Correlation functions at b 1 ˆ hX. (a) Symmetrized correlation function. Solid line and markers correspond to Dkmax ˆ 6 and 9, respectively. (b) Original correlation function. Solid and dashed lines show real and imaginary parts calculated from the data for C…t† with Dkmax ˆ 6.

Fig. 8. Correlation functions for an asymmetric TLS with e ˆ X at b 1 ˆ hX. (a) Symmetrized correlation function. Solid line and markers correspond to Dkmax ˆ 9 and 10, respectively. (b) Original correlation function. Solid and dashed lines show real and imaginary parts calculated from the data for C…t† with Dkmax ˆ 9.

G…t†, is still observable because the dissipation is weak. As the temperature decreases, the real part of G…t† di€ers drastically from C…t† and quantum coherence becomes more prominent. All these effects are clearly exhibited in Fig. 7. The dynamical behavior of the asymmetric system (e ˆ X) is displayed in Fig. 8 under the same conditions. One can see that while the real part of G…t† assumes a non-vanishing stationary value, its imaginary part goes to zero after several oscillations. For the asymmetric system convergence is achieved with a slightly larger value Dkmax .

4. Discussion Just as in situations where a macroscopic bath is initially at thermal equilibrium, the length of nonlocal interactions in the in¯uence functional is ®nite also when the initial condition corresponds to the Boltzmann operator for the coupled system± bath Hamiltonian. This convenient behavior arises because the in¯uence functional for an equilibrium correlation function is determined by the response function of the bath in complex time. Under conditions resulting in dissipative dynamics for the

10

J. Shao, N. Makri / Chemical Physics 268 (2001) 1±10

system, the response function of the environment decays irreversibly in real time. This feature alone would be sucient to imply decay of the nonlocal interactions along any complex time contour. However, the thermal ¯uctuations of the environment represent an additional source of memory quenching along the imaginary time axis, such that the response function exhibits nonlocal behavior of a ®nite span even at zero time and even under non-dissipative conditions, i.e., if the bath consists of a single oscillator. These characteristics imply that the partition function itself is amenable to an iterative algorithm. Together, the decay of in¯uence functional correlations due to thermal ¯uctuations and also due to dynamical decoherence in dissipative baths allow evaluation of the correlation function by an iterative procedure. By reordering the evolution contour of the numerical path integral we have recently developed a more ecient algorithm for the calculation of multi-time equilibrium correlation functions. Acknowledgements This work has been supported by the National Science Foundation under award no. NSF CHE9877050. References [1] R.P. Feynman, Rev. Mod. Phys. 20 (1948) 367±387. [2] R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965.

[3] R.P. Feynman, J.F.L. Vernon, Ann. Phys. 24 (1963) 118± 173. [4] A.O. Caldeira, A.J. Leggett, Physica A 121 (1983) 587±616. [5] S. Dattagupta, Relaxation Phenomena in Condensed Matter Physics, Academic Press, New York, 1987. [6] A.G. Red®eld, Adv. Mag. Reson. 1 (1965) 1. [7] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fisher, A. Garg, M. Zwerger, Rev. Mod. Phys. 59 (1987) 1±85. [8] R.D. Coalson, J. Chem. Phys. 86 (1987) 995±1009. [9] M. Topaler, N. Makri, J. Chem. Phys. 101 (1994) 7500± 7519. [10] N. Makri, W.H. Miller, J. Chem. Phys. 89 (1988) 2170± 2177. [11] J.D. Doll, D.L. Freeman, Adv. Chem. Phys. 73 (1988) 289± 304. [12] B.A. Mason, K. Hess, Phys. Rev. B 39 (1989) 5051±5069. [13] C.H. Mak, D. Chandler, Phys. Rev. A 44 (1991) 2352± 2369. [14] C.H. Mak, Phys. Rev. Lett. 68 (1992) 899±902. [15] N. Makri, Comp. Phys. Commun. 63 (1991) 389±414. [16] D.E. Makarov, N. Makri, Chem. Phys. Lett. 221 (1994) 482±491. [17] N. Makri, D.E. Makarov, J. Chem. Phys. 102 (1995) 4600± 4610. [18] N. Makri, D.E. Makarov, J. Chem. Phys. 102 (1995) 4611± 4618. [19] E. Sim, N. Makri, Comp. Phys. Commun. 99 (1997) 335± 354. [20] N. Makri, J. Phys. Chem. 102 (1998) 4414±4427. [21] E. Sim, N. Makri, Chem. Phys. Lett. 249 (1996) 224±230. [22] A.A. Golosov, R.A. Friesner, P. Pechukas, J. Chem. Phys. 110 (1999) 138±146. [23] N. Makri, J. Chem. Phys. 111 (1999) 6164±6167. [24] U. Weiss, Quantum Dissipative Systems, World Scienti®c, Singapore, 1993. [25] W.H. Miller, S.D. Schwartz, J.W. Tromp, J. Chem. Phys. 79 (1983) 4889±4898. [26] K. Okumura, Y. Tanimura, Phys. Rev. E 56 (1997) 2747. [27] N. Makri, Chem. Phys. Lett. 193 (1992) 435±444. [28] M. Topaler, N. Makri, Chem. Phys. Lett. 210 (1993) 285± 293.