Condensed-Phase Relaxation of Multilevel Quantum Systems. I. An

18758

J. Phys. Chem. B 2006, 110, 18758-18763

Condensed-Phase Relaxation of Multilevel Quantum Systems. I. An Exactly Solvable Model† Simone Peter,‡ Deborah G. Evans,*,§ and Rob D. Coalson‡ Department of Chemistry, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15260, and Department of Chemistry, UniVersity of New Mexico, Albuquerque, New Mexico 87131 ReceiVed: February 24, 2006; In Final Form: May 13, 2006

An analytically solvable model of multilevel condensed-phase quantum dynamics relevant to vibrational relaxation and electron transfer is presented. Exact solutions are derived for the reduced system density matrix dynamics of a degenerate N-level quantum system characterized by nearest-neighbor hopping and off-diagonal coupling (which is linear in the bath coordinates) to a harmonic oscillator bath. We demonstrate that for N > 2 the long-time steady-state system site occupation probabilities are not the same for all sites; that is, they are distributed in a non-Boltzmann manner, which depends on the initial conditions and the number of levels in the system. Although the system-bath Hamiltonian considered here is restricted in form, the availability of an exact solution enables us to study the model in all regions of an extensive parameter space.

1. Introduction Models of a small quantum system coupled to a condensed phase bath have played an important role in understanding a host of problems in chemical physics, ranging from spectral line broadening of a solute1 in the solution phase to exciton hopping rates in solid-state materials and energy transfer in molecular assemblies.2-6 Because of the wide variety of practical applications, the general theory of quantum relaxation in a system coupled to a bath has received considerable attention over the last three decades. A series of seminal papers by Silbey, Harris, and many others7-15 demonstrated the utility of approximate methods for solving such many-body quantum problems. Applications of perturbative approaches first introduced by Kubo16 to problems in solid-state physics have been applied to a range of problems in physical chemistry and molecular spectroscopy. The development of numerous approximate methods, notably the Redfield approach,17-21 Pauli master equations,22-27 the Nakajima-Zwanzig16,28,29 projection operator formalism, and new path integral algorithms,30,31 have allowed theorists to compare and test the efficacy of these methods for accurately describing condensed-phase dynamical processes. In a recent paper,32 we focused attention on a class of nontrivial N-level Hamiltonians coupled to a dissipative environment (or “bath”) for which the time evolution of the composite quantum many-body system can be calculated exactly. In particular, we showed that the reduced system density matrix of an N-level condensed-phase system can be obtained as a linear combination of N2 time-dependent functions, which can be calculated exactly for certain types of baths and systembath couplings. This technique has been used recently by Cook, Evans, and Coalson33 to study non-Condon effects in condensedphase electron transfer and constitutes a generalization of simpler Hamiltonians discussed by Creechly and Dahnovsky34 and Gayen et al.35 in a different context. Although the range of exactly solvable Hamiltonians is limited with respect to the form of the system-bath coupling and the nature of the bath modes, the advantages of having an †

Part of the special issue “Robert J. Silbey Festschrift”. * Corresponding author. E-mail: [email protected]. ‡ University of Pittsburgh. § University of New Mexico.

exactly solvable class of system-bath Hamiltonians for a condensed-phase problem are obvious: such exact solutions provide a testing ground for various approximate solutions to the same problem. In this way, the accuracy of approximate condensedphase quantum dynamics algorithms can be assessed and the exact solutions can be used to study time evolution under particular system-bath Hamiltonians in all regions of parameter space, beyond the weak coupling limit most frequently tackled by approximate methods. In this paper, we derive a set of analytical formulas to describe quantum relaxation of a degenerate N-level system in the presence of a particular type of dissipative bath at finite temperature. Given an initially prepared system state, the evolution of the reduced density matrix elements of the system is calculated exactly and explicitly for all system-bath coupling strengths. There are several potentially important applications of this work to model physical situations where a general N-level system Hamiltonian is coupled to a high-dimensional bath. For example, the proposed models can be used to understand bathmodulated coupling in physical situations including, for example, vibrational relaxation of an N-level system in a solvent or other dissipative medium,36,37 electron transfer rates in multisite “donor-bridge-acceptor” systems where the bridge vibrations are strongly coupled to the electronic Hamiltonian,27,38-40 and in chromophore assemblies where the interchromophore distances are modulated by low-frequency phonons in the surrounding medium.5,41 The outline of the present paper is as follows: In Section 2, the model system-bath Hamiltonian is introduced, and the general method for solving the associated Schro¨dinger equation is reviewed. In Section 3, exact formulas are derived for the dynamics of a degenerate N-level system as a function of initialstate preparation and the explicit form of the system-bath coupling. These formulas are then used to study the quantum dynamics of both an isolated system and a system coupled to a dissipative bath. Finally, Section 4 contains further discussion and conclusions. 2. Model Hamiltonian and Physical Significance Consider a chain of N sites, placed adjacent to each other at equidistant intervals. Label these sites according to their position

10.1021/jp061198y CCC: $33.50 © 2006 American Chemical Society Published on Web 08/09/2006

Relaxation of Multilevel Quantum Systems. I.

J. Phys. Chem. B, Vol. 110, No. 38, 2006 18759

along the chain: 1,2,...N. All sites correspond to the same energy, which can be taken without loss of generality to be 0. Nearest-neighbor sites are coupled by a hopping matrix element, R. For example, for a three-site chain (N ) 3): H ˆ s ) R Rˆ s, where

( )

0 1 0 Rˆ s ) 1 0 1 0 1 0

(1)

In general, we want to compute the quantum dynamics that takes place when this N-site system is coupled to a condensed phase environment. The total system-bath Hamiltonian, H ˆ , can then be written as

ˆb H ˆ )H ˆ s + Rˆ s Rˆ b + H

(2)

where H ˆ b is the isolated bath Hamiltonian and Rˆ b is the bath operator factor in the system-bath coupling term. Equation 2 appears to be rather general in structure, but in fact we will consider here only nearest-neighbor hopping operators of the type indicated in eq 1 for Rˆ s and will take H ˆ s to be proportional to Rˆ s. (The general class of exactly solvable models that can be treated by the analysis of ref 32 requires that [H ˆ s, Rˆ s] ) 0.) Furthermore, we will choose H ˆ b to be a collection of uncoupled unit mass harmonic oscillators and Rˆ b to be a linear combination of bath coordinates; that is

H ˆb )

1 2

∑j (

pˆ 2j

+

ω2j

qˆ 2j )

and Rˆ b )

∑j cj qˆ j

where each oscillator, characterized by frequency ωj, is coupled to the system coordinate via cj as prescribed by a specified spectral density function

J(ω) )

π 2

c2j

∑j ω δ(ω - ωj) j

where δ is the Dirac delta function. Again, these features of the bath Hamiltonian and the bath factor in the system-bath coupling operator are required for analytical solubility of the quantum dynamics of the multidimensional model. Note that because the system factor in the system-bath coupling term, Rˆ s, consists of nearest-neighbor hopping matrix elements, the physical effect of the system-bath coupling under consideration is to add fluctuations to the hopping matrix elements. The details of the fluctuations depend on the details of H ˆ b and Rˆ b, and the role of the bath in the reduced system density matrix is entirely captured by the bath Heisenberg autocorrelation function 〈Rˆ b(t) Rˆ b〉,42 which for the class of Hamiltonians under consideration here is determined by the spectral density function defined above. As noted above, the key feature that provides the exact solution is the condition that H ˆ s and Rˆ s commute.32 In this paper, we consider initial density operators of the form Γˆ (0) ) Fˆ β X Fˆ (0), where Γˆ (0) is the complete system-bath density matrix and Fˆ (0) (unsubscripted for notational convenience) is the initial system density operator; furthermore, Fˆ β is the thermal density operator for the bath at inverse temperature β ) 1/kBT; that is, Fˆ β ) exp(-βH ˆ b)/Zb with Zb ) trb{exp(-βH ˆ b)} (Note: kB is the Boltzmann constant, T is the absolute temperature, and trb indicates a trace over all bath degrees of freedom.) In this paper, we seek to extract elements of the reduced system density operator Fˆ (t) ) trb[Γˆ (t)], where Γˆ (t) is the time evolution of the complete system-bath density operator under the evolu-

tion of the full Hamiltonian; that is, Γˆ (t) ) exp(-iH ˆ t) Γˆ (0) exp(iH ˆ t) (we set p ) 1 throughout). Although the reduced system density operator can in principle be expressed in any basis of N mutually orthogonal vectors that span the N-level system Hilbert space, for the remainder of this paper we will use the notation Fˆ (t) to designate the N × N reduced system density matrix in the basis of site-localized states |J〉, that is, [Fˆ (t)]JK ) Fˆ JK(t) ) trb(〈J| Γˆ (t)|K〉). It is convenient to transform to an interaction picture representation with respect to the system Hamiltonian in order to carry out the evaluation of the desired time propagation. Let Fˆ I(t) be the resultant interaction-picture reduced system density matrix in the site basis, defined by Fˆ I(t) ) exp(iH ˆ st) Fˆ (t) exp(-iH ˆ st), with H ˆ s represented in the site basis (as in eq 1). Further, it proves useful to transform this interaction picture density matrix to the basis of eigenstates of Rˆ s. Specifically, if we denote the matrix of unit-normed eigenvectors of Rˆ s as U ˆ (with one real-valued eigenvector of Rˆ s in each column of U ˆ ), then F˜ I(t), the interaction picture reduced system density matrix expressed in the basis of Rˆ s eigenstates, is given by F˜ I(t) ) U ˆ TFˆ I(t)U ˆ . As shown in ref 32, the time evolution of the elements of F˜ I(t) is explicitly prescribed as

F˜ µνI(t) ) Aµ,ν(t) F˜ µνI(0)

(3)

with t

Aµ,ν(t) ) e

-(λµ-λν)2

t′

∫ ∫ dt′

0

t

dt′′κ(t′-t′′)+i(λ2µ-λ2ν)

0

t′

∫ ∫dt′′ξ(t′-t′′) dt′

0

(4)

0

where κ(t) - iξ(t) ) 〈Rˆ b(t)Rˆ b〉, and λµ is the eigenvalue corresponding to the µth eigenvector of Rˆ s. Given F˜ Iµν(t), it is straightforward to obtain the reduced system density matrix in the Schro¨dinger picture and site basis by reversing the transformations noted above. Again, we stress that this constitutes an exact solution of the time-evolution of the reduced system density matrix for a harmonic oscillator bath and a bath coupling operator Rˆ b which is linear in the oscillator coordinates. 3. Quantum Evolution in Multisite Systems The basic question of interest in this paper is the following: if we prepare the system with the population localized at site I (1 < I < N) and the bath in thermal equilibrium, then how does the reduced system density matrix evolve in time? In this Section, we will focus on deriving exact formulas for the reduced system density matrix in the site basis, FJK(t). To do this, we need to evaluate the eigenvalues and eigenvectors of Rˆ s. For the model directly under consideration in this paper (cf. eq 1), Rˆ s is a Jacobi matrix. Its eigenvalues are given by43

λk ) -2 cos(kθ), θ ≡

π ; k ) 1, 2, ..., N N+1

(5)

The corresponding matrix of eigenvalues (with eigenvector k in column k) is given by

Uj,k )

xN +2 1 sin(jkθ);

j, k ) 1, 2, ..., N

(6)

Note that this U ˆ has the convenient (if somewhat unusual) property that U ˆ )U ˆT ) U ˆ -1. For concreteness, we focus on the case that the system is prepared in pure state I (i.e., localized at site I) so that FJK(0) ) δI,JδI,K, with δ being the Kronecker delta symbol. The exact reduced system density matrix at time t is then given for any

18760 J. Phys. Chem. B, Vol. 110, No. 38, 2006

Peter et al.

prescribed oscillator bath (defined by its spectral density) as

FLM(t) )

( ) ∑∑ 2 N

2

N+1

N

e-iR(λj-λr)t Aj,r(t) sin(Ljθ)

j)1 r)1

sin(Mrθ) sin(Ijθ) sin(Irθ) (7) where Aj,r (t) is the function prescribed in eq 4, with λk taken from eq 5. This formula includes the case of no dissipation as one obvious limit (Aj,r (t) ) 1), which we will examine in Section 3.1, before moving on to study the dissipative case in Section 3.2. 3.1. Evolution of the Isolated N-Site System. In the absence of coupling to a dissipative environment, the site occupation probabilities PK(t) ≡ FKK(t) are given (cf. eq 7) by

PK(t) )

( )

2 N

2

N+1

N

∑ ∑ cos[R(λj - λr)t] sin(Kjθ) j)1 r)1 sin(Krθ) sin(Ijθ) sin(Irθ) (8)

Figure 1. Time evolution of probability distribution PJ(t) ) FJJ(t) of a nondissipative N-level system for initial conditions FJK(0) ) δI,JδI,K and thus PJ(0) ) δI,J. Here N ) 49, I ) 25, and R ) 1.

From these, one can compute the time-dependent mean square site displacement: N

〈[K - I]2(t)〉 )

∑ (K - I)2 PK(t)

(9)

K)1

The case of motion according to H ˆ s ) R Rˆ s in eq 1 is intimately connected to that of free-particle wave-packet motion in 1D, because of the form of the lattice discretized second derivative operator, that is, ψ′′j = (ψj+1 - 2ψj + ψj-1)/2, where j labels the lattice position of the x coordinate (xj ) j,  being the lattice spacing). The discretized free-particle Hamiltonian is thus equivalent (to within a physically irrelevant overall phase factor) to H ˆ s in eqs 1 and 2. One feature that slightly blurs the analogy between the sitehopping dynamics under study and free-particle motion of a 1D wave packet is our choice of initial conditions, which correspond approximately to δ(x - xI), with δ here indicating the Dirac delta function. If we were working with a smoothedout initial Gaussian wave packet (spanning several sites), then we could appeal directly to the well-known expression for the time evolution of an initially Gaussian wave packet under the free-particle Hamiltonian,44 which remains Gaussian (implying a Gaussian probability distribution) for all times. In the limit that the initial wave packet collapses to an infinitely narrow delta function, the relevant time evolution is that of the freeparticle propagator, which for t > 0 extends infinitely in space and is characterized by rapid oscillation as one moves away from the location of the t ) 0 delta function blip. Because we are working with a coarse-grained lattice, it might be expected that the behavior of our lattice-hopping Hamiltonian, H ˆ s, would be intermediate between these extremes, and indeed, that is what is found numerically. As illustrated in Figure 1, the probability distribution is confined at early times to an envelope that spreads out around the initial point of localization; within this envelope, the probability distribution oscillates rapidly. Eventually, the envelope expands until it “hits the walls” at the end of the lattice, and then probability is reflected back into the interior of the lattice. As is generally true of discrete N-level quantum systems, the site probabilities never reach a long-time steady state, but rather fluctuate indefinitely. To illustrate some of these properties, it is useful to plot the mean square site displacement as a function of time. This is done in Figure 2. Two features are particularly worthy of note. First, the mean squared fluctuation grows quadratically in time

Figure 2. Time evolution of mean square deviation from initially prepared state I (i.e., FJK(0) ) δI,JδI,K) for an isolated N-level system. Shown are results for N ) 39, I ) 20 and N ) 49, I ) 25, as well as the short-time quadratic approximation: 〈[K - I]2(t)〉 = 2R2t2.

according to the formula 〈[K - I]2(t)〉 = 2R2t2 until the wavepacket edges reach the wall, after which the mean square displacement associated with the probability distribution begins to oscillate. Second, as the number of sites, N, is increased, the early-time behavior is the same as that for smaller N, but the “break time” is larger, corresponding to the longer time that it takes for the wave packet to reach the ends of the lattice. To deduce the curvature coefficient in the quadratic time dependence of the early time dynamics of the mean square fluctuation, it is instructive to write

cos[R(λj - λr)t] ) 1 + (cos[R(λj - λr)t] - 1) in eq 8. Then, making the short time approximation

cos[R(λj - λr)t] - 1 = -[R(λj - λr)t]2/2 the following approximation to the probability distribution is obtained

PK(t) = δI,K + (Rt)2[δK,I+1 - 2δK,I + δK,I-1]

(10)

which implies the quadratic formula for the mean square displacement given above. Clearly, eq 10 is truly a short-time approximation, which becomes inaccurate for times after which more than just the first nearest neighbors adjacent to site I are populated (for N ) 49, the case illustrated in Figure 1, this corresponds to Rt ≈ 0.5). Curiously, however, the mean squared displacement follows the quadratic dependence derived here for

Relaxation of Multilevel Quantum Systems. I.

J. Phys. Chem. B, Vol. 110, No. 38, 2006 18761

Figure 3. Time evolution of site population distribution PJ(t) ) FJJ(t) for a dissipative N-level system, given initial system preparation on site I (i.e., FJK(0) ) δI,JδI,K). Here N ) 49, I ) 25 with R ) 0, Γ ) Ω ) 1.

short time motions up to much longer times, namely, until the system wave packet reaches the edges of the lattice (for N ) 49, this corresponds to Rt ≈ 13). 3.2. Evolution of an N-Site System Coupled to a Dissipative Bath. Returning to eq 7, when the N-site system is coupled to the dissipative bath as specified in the discussion following eq 2, the details of the quantum dynamics depend on the Heisenberg correlation function 〈Rˆ b(t) Rˆ b〉. Here we invoke the following real-valued correlation function

〈Rˆ b(t) Rˆ b〉 ) Ω2e-Γ |t|

(11)

with parameters Ω, Γ > 0, which can be regarded as the high T limit associated with an appropriate oscillator spectral density (obtained essentially by Fourier transformation of eq 11; see ref 32 for full details). Clearly, the value of R determines the balance between isolated N-site dynamics and relaxation dynamics induced by coupling to the harmonic oscillator bath. Because we wish here to study effects of coupling to dissipative bath, we shall focus primarily on the case that R ) 0 (the case of nonzero R will be considered briefly at the end of this subsection). Adopting the form of the Heisenberg correlation function given in eq 11 and setting R ) 0, the site occupation probabilities at time t are given by

PK(t) )

( ) ∑∑ 2

N+1

2 N

N

e-(λj-λr) F(t) sin(Kjθ) sin(Krθ) 2

Figure 4. Time evolution of mean square deviation from initial state I (i.e., FJK(0) ) δI,JδI,K) of a dissipative N-level system for N ) 45, I ) 23 and N ) 39, I ) 20, with R ) 0, Γ ) Ω ) 1. Also shown is the short-time approximant 〈[K - I]2(t)〉 = 4F(t).

Figure 3 illustrates the time evolution of the probability distribution characterizing the system for the case that N ) 49 and I ) 25 (corresponding to initial preparation in the center of the lattice). The probability distribution spreads smoothly as time progresses, until the wings of the distribution reach the edges of the lattice. At this point in time, the system begins the slow process of “settling in” to the asymptotic steady-state distribution, which is characterized almost [though not quite (vide infra)] by equal population of all sites. To further characterize this dynamics, we examine the mean square fluctuation of site occupation probabilities. The time evolution of this quantity for N ) 45, I ) 23 and N ) 39, I ) 20 is shown in Figure 4. Note that it has quadratic time dependence at very short times (t < 2) and then exhibits a linear dependence at intermediate times, until it levels off to a longtime limiting value. This limit value is given approximately by N2/12, as discussed in detail below. As is clear from the figure, the presaturation time evolution of the mean square displacement is well-approximated by the following functional form:

〈[K - I ]2(t)〉 = 4F(t)

(13)

At early times this corresponds to the quadratic time dependence 2Ω2t2, whereas at long times the rhs of eq 13 becomes linear in time, namely, 4Ω2t/Γ. The mean square fluctuation behavior just described can be anticipated by carrying out a short-time approximation to the site probability distribution analogous to the procedure considered in the absence of dissipative coupling (but with nonzero R) above. Here we refer to eq 12, and write

e-(λj-λr) F(t) ) 1 + (e-(λj-λr) F(t) - 1) 2

j)1 r)1

sin(Ijθ) sin(Irθ) (12a) with

2

Then, making the short time approximation

e-(λj-λr) F(t) -1 = -(λj - λr)2F(t) 2

Ω2 F(t) ≡ 2 (e-Γt + Γt - 1) Γ

(12b)

For such a system, the dynamics is purely overdamped (with no oscillatory characteristics), leading ultimately to a welldefined long-time steady state. Full details of the time evolution depend on the choice of Ω and Γ, but because F(t) increases monotonically from 0 to ∞, then for any {Γ,Ω} the system traces through the same sequence of reduced system density matrix states (including probability distributions). (Obviously, increasing Ω for fixed Γ accelerates progress toward the long-time steady state, etc.) So, for simplicity, we will present numerical results for the case where Ω ) Γ ) 1.

the following approximation to the probability distribution is obtained

PK(t) = δI,K + 2F(t)[δK,I+1 - 2δK,I + δK,I-1]

(14)

which implies 〈[K - I]2(t)〉 ) 4F(t), as noted above. As in the nondissipative case, we caution that eq 14 holds only at very short times (up to ca. t ) 0.5 for the case of N ≈ 50 considered in Figures 3 and 4). However, the time-dependence of the mean squared site fluctuation prescribed by eq 13 appears to hold for considerably longer times, in particular up to the time where the site occupation probability distribution spreads to the edge

18762 J. Phys. Chem. B, Vol. 110, No. 38, 2006

Peter et al.

Figure 5. Time evolution of mean squared deviation from initial state I (i.e., FJK(0) ) δI,JδI,K) for N ) 17, I ) 9 in a dissipative N-level system (Γ ) Ω ) 1), for several values of R: (i) R ) 0 (-[-), (ii) R ) 1 (-9-), and (iii) R ) 2 (-2-).

of the site lattice (after which the mean squared fluctuation begins to saturate to a long-time limiting value). The details of the long-time-limit site probabilities merit further discussion. If there is any dissipation at all (i.e., 〈Rˆ b(t)Rˆ b〉 * 0), then Aj,r(t) f δj,r as t f ∞ so that the exact long-time site populations are given by

PK(∞) )

( )∑ 2 N

2

N+1

sin2(Kjθ) sin2(Ijθ)

(15)

j)1

As is apparent from Figure 3, the probability to be at any site becomes nearly identical for large N, that is, Pk(∞) = 1/N. In this approximation, taking for concreteness the initial condition I ) (N + 1)/2 (preparation at the center of the lattice), then

〈(

K-

)〉

N+1 2

2

)

1

(N-1)/2

∑

N J)-(N-1)/2

J2 ) (N2 - 1)/12 = N2/12 (16)

as noted above. In fact, for finite N, eq 15 implies small corrections to this result. For I ) (N + 1)/2, it turns out that all site populations are equally populated, except for site I, which has twice the occupation probability as any of the others. (The perturbation from uniformity is slightly different for I * (N + 1)/2.) However, as N f ∞ these tiny “glitches” in the probability distribution have a negligible effect on its global characteristics: in particular, eq 16 is well-satisfied. Finally, we briefly consider the case that R * 0 so that

e-(λj-λr) F(t) f e-(λj-λr) F(t) cos(R[λj - λr]t) 2

2

in eq 12. This mixes an oscillatory component into the dynamics of the site probability distribution (reminiscent of the nondissipative limit case). However, as is clear from the discussion above, the long-time limit of the distribution is independent of the value of R. To illustrate this behavior, we show in Figure 5 the time evolution of 〈[K - I ]2(t)〉 for the case N ) 17, I ) 9, Γ ) Ω ) 1, and several values of R. 4. Discussion and Conclusions In this paper, we have analyzed an exactly solvable model of a degenerate N-level quantum system coupled to a dissipative bath. This model is a many-level generalization of exactly solvable two-level condensed-phase system models studied

recently,32,34 and for N ) 2, the analytical formulas derived here reduce appropriately. We have shown that for N > 2 a nonuniform (thus non-Boltzmann) site population is obtained at long time, with the precise populations depending on the details of initial system density matrix (and thus, among other factors, on the number of levels in the system). Although the exactly soluble system-bath Hamiltonians considered here are restricted to certain forms, we are able to study the quantum dynamics of systems described by such Hamiltonians in a wide range of parameter space for all times. In a companion paper (Peter, S.; Evans, D. G.; Coalson, R. D. J. Phys. Chem. B. 2006, 110, 18764.), we explore the generalization of this model to a nondegenerate N-level system coupled in a similar fashion to a harmonic oscillator bath, which is a multidimensional quantum dynamics problem that does not admit analytical solutions, but where numerically exact path integral calculations can be used to compare the reduced system density matrix dynamics with various perturbative approximations to it. The analytical solutions presented here provide impetus for examining these more complicated models. Acknowledgment. With best wishes to Prof. Bob Silbey on occasion of his 60th birthday and in recognition of his seminal contributions to the field of condensed-phase quantum dynamics. D.G.E. is a Cottrell Scholar of the Research Corporation and a Camille Dreyfus Teacher-Scholar. D.G.E.’s work is partially supported by a grant from the Petroleum Research Fund (PRF no. 43207-AC10) administered by the ACS. R.D.C.’s work was partially supported by a grant from the National Science Foundation. References and Notes (1) Oxtoby, D. W. Annu. ReV. Phys. Chem. 1981, 32, 77-101. (2) Jean, J. M. J. Chem. Phys. 1994, 101, 10464-10473. (3) Jean, J. M.; Fleming, G. R. J. Chem. Phys. 1995, 103, 2092-2101. (4) Jean, J. M. J. Phys. Chem. A 1998, 102, 7549-7557. (5) Cina, J. A.; Fleming, G. R. J. Phys. Chem. A 2004, 108, 1119611208. (6) Pereverzev, A.; Bittner, E. R. J. Chem. Phys. 2005, 123, 244903. (7) Harris, R. A.; Silbey, R. J. Chem. Phys. 1983, 78, 7330. (8) Parris, P. E.; Silbey, R. J. Chem. Phys 1985, 83, 5619. (9) Rackovsky, S.; Silbey, R. Mol. Phys. 1973, 25, 61. (10) Reichman, D. R.; Silbey, R. J. J. Chem. Phys. 1996, 104, 1506. (11) Silbey, R. Annu. ReV. Phys. Chem. 1976, 27, 203. (12) Silbey, R.; Harris, R. A. J. Chem. Phys. 1984, 80, 2615. (13) Silbey, R.; Harris, R. A. J. Phys. Chem. 1989, 93, 7062-7071. (14) Jang, S.; Cao, J. S.; Silbey, R. J. J. Chem. Phys. 2002, 116, 27052717. (15) Jang, S.; Newton, M. D. J. Chem. Phys. 2005, 122, 024501. (16) Kubo, R.; Toyozawa, Y. Prog. Theor. Phys. 1955, 13, 160-182. (17) Redfield, A. AdV. Magn. Reson. 1965, 1, 1. (18) Jean, J. M.; Friesner, R. A.; Fleming, G. R. J. Chem. Phys. 1992, 96, 5827-5842. (19) Pollard, W. T.; Felts, A. K.; Friesner, R. A. AdV. Chem. Phys., Vol XCIII 1996, 93, 77-134. (20) Walsh, A. M.; Coalson, R. D. Chem. Phys. Lett. 1992, 198, 293299. (21) Gaspard, P.; Nagaoka, M. J. Chem. Phys. 1999, 111, 5668-5675. (22) Sparpaglione, M.; Mukamel, S. J. Chem. Phys. 1988, 88, 32633280. (23) Golosov, A. A.; Friesner, R. A.; Pechukas, P. J. Chem. Phys. 1999, 110, 138-146. (24) Laird, B. B.; Budimir, J.; Skinner, J. L. J. Chem. Phys. 1991, 94, 4391-4404. (25) Davies, E. B. Commun. Math. Phys. 1974, 9, 91-110. (26) Dumcke, R.; Spohn, H. Z. Phys. B: Condens. Matter 1979, 34, 419-422. (27) Egorova, D.; Kuhl, A.; Domcke, W. Chem. Phys. 2001, 268, 105120. (28) Zwanzig, R. Physica 1964, 30, 1109. (29) Kohen, D.; Marston, C. C.; Tannor, D. J. J. Chem. Phys. 1997, 107, 5236-5253. (30) Wilkie, J. J. Chem. Phys. 2001, 114, 7736-7745.

Relaxation of Multilevel Quantum Systems. I. (31) Yan, Y. J.; Shuang, F.; Xu, R. X.; Cheng, J. X.; Li, X. Q.; Yang, C.; Zhang, H. Y. J. Chem. Phys. 2000, 113, 2068-2078. (32) Coalson, R. D.; Evans, D. G. Chem. Phys. 2004, 296, 117. (33) Cook, W. R.; Evans, D. G.; Coalson, R. D. Chem. Phys. Lett. 2006, 420, 362-366. (34) Creechley, J.; Dahnovsky, Y. Chem. Phys. 2004, 296, 171. (35) Gayen, T.; Mcdowell, K.; Burns, A. J. Chem. Phys. 2000, 112, 4310-4320. (36) Nitzan, A.; Jortner, J. Chem. Phys. Lett. 1972, 15, 350-&. (37) Nitzan, A.; Jortner, J. J. Chem. Phys. 1973, 58, 2412-2434. (38) Skourtis, S.; Archontis, G.; Xie, Q. J. Chem. Phys. 2001, 115, 9444.

J. Phys. Chem. B, Vol. 110, No. 38, 2006 18763 (39) Garg, A.; Onuchic, J. N.; Ambegaokar, V. J. Chem. Phys. 1985, 83, 4491-4503. (40) Scherer, T.; Vanstokkum, I.; Brouwer, A.; Verhoeven, J. J. Phys. Chem. 1994, 98, 10539. (41) Jang, S.; Jung, Y. J.; Silbey, R. J. Chem. Phys. 2002, 275, 319332. (42) Schatz, G. C.; Ratner, M. A. Quantum Mechanics in Chemistry; Dover Publications: Mineola, NY, 2002. (43) Coulson, C. A. Proc. R. Soc. (London) 1938, A164, 383. (44) Heller, E. J. J. Chem. Phys. 1981, 75, 2923.