Energy relaxation of a dissipative quantum oscillator Pradeep Kumar and Eli Pollak Citation: The Journal of Chemical Physics 141, 234509 (2014); doi: 10.1063/1.4903809 View online: http://dx.doi.org/10.1063/1.4903809 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/23?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Femtosecond transient infrared and stimulated Raman spectroscopy shed light on the relaxation mechanisms of photo-excited peridinin J. Chem. Phys. 142, 212409 (2015); 10.1063/1.4915072 Probing non-adiabatic conical intersections using absorption, spontaneous Raman, and femtosecond stimulated Raman spectroscopy J. Chem. Phys. 139, 234101 (2013); 10.1063/1.4843395 Direct determination of state-to-state rotational energy transfer rate constants via a Raman-Raman double resonance technique: ortho-acetylene in v 2 = 1 at 155 K J. Chem. Phys. 132, 154303 (2010); 10.1063/1.3374031 Efficient stimulated Raman pumping for quantum state resolved surface reactivity measurements Rev. Sci. Instrum. 77, 054103 (2006); 10.1063/1.2200876 Role of electronic‐vibrational interaction in the short‐time coherent and long‐time dissipative dynamics of energy relaxation AIP Conf. Proc. 643, 110 (2002); 10.1063/1.1523790
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
THE JOURNAL OF CHEMICAL PHYSICS 141, 234509 (2014)
Energy relaxation of a dissipative quantum oscillator Pradeep Kumar and Eli Pollak Chemical Physics Department, Weizmann Institute of Science, 76100 Rehovot, Israel
(Received 29 September 2014; accepted 28 November 2014; published online 19 December 2014) The dissipative harmonic oscillator is studied as a model for vibrational relaxation in a liquid environment. Continuum limit expressions are derived for the time-dependent average energy, average width of the population, and the vibrational population itself. The effect of the magnitude of the solute-solvent interaction, expressed in terms of a friction coefficient, solvent temperature, and initial energy of the oscillator on the relaxation has been studied. These results shed light on the recent femtosecond stimulated Raman scattering probe of the 1570 cm−1 −C=C− stretching mode of trans-Stilbene in the first (S1 ) excited electronic state. When the oscillator is initially cold with respect to the bath temperature, its average energy and width increase in time. When it is initially hot, the average energy and width decrease with time in qualitative agreement with the experimental observations. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4903809] I. INTRODUCTION
A vibrationally excited molecule will eventually relax to an equilibrium state in the presence of a thermal bath. The evolution of this relaxation depends upon many factors such as solute-solvent coupling, intramolecular-vibrationalrelaxation, temperature of the bath and the initial state of the molecule.1–10 Many experiments have been carried out to probe the time-dependent vibrational relaxation of the Stilbene molecule.3–9, 11–17 The measured data on a given vibrational band are characteristically summarized in terms of the time evolution of the instantaneous maximum of the band and its full width at half maximum (FWHM). The early experimental studies were carried out by combining ultrafast spectroscopy techniques in a supersonic jet.13–15 Hamaguchi and co-workers4–9 used time-dependent Raman spectroscopy from the S1 to the S2 electronic states, to measure the time-dependent linewidth and peak position of the Raman spectra of the 1570 cm−1 −C=C− stretching mode of trans-Stilbene. They observed a linear correlation between the time-dependent evolution of the peak position and the linewidth. They found a narrowing of the linewidth in time until reaching its equilibrium value. The experimental studies of Hamaguchi and co-workers were limited to the picosecond domain. The advent of femtosecond-stimulated-Raman-spectroscopy (FSRS) enabled Ernsting and co-workers to extend the time-dependent study of the −C=C− stretching mode of trans-Stilbene to the femtosecond domain.3, 16, 17 In a recent study, Kovalenko et al.3 found that for a short time scale (0-3 ps) the time evolution of the line shape changes dramatically as a function of the initial excitation energy. When the initial energy is high, in the sense that it may excite a few quanta of the mode, the linewidth first broadens and only later narrows, as observed by Hamaguchi and co-workers. When the initial excitation populates only the ground state then the linewidth only broadens with time until reaching its equilibrium value. There have been a few attempts to understand these experimental observations theoretically.5, 10, 18–21 Hamaguchi 0021-9606/2014/141(23)/234509/11/$30.00
proposed a “dynamic polarization model” in which the time evolution is dominated by the stochastic polarization of the mode, induced by the solvent.18 Kovalenko et al.3 suggested that the change in the dynamics, from broadening of the FWHM at low excitation energies to narrowing at high excitation energies serves as an indication of the initial temperature of the excited mode. A highly excited mode couples with many solvent degrees of freedom so that its spectrum is broad. As it cools down, the number of coupled modes decreases and the width narrows. Conversely, when the mode is cold, it couples only to a few bath modes, subsequent heating by the bath then broadens the peak. This explanation, if correct should be of general validity and should not depend too much on the nature of the coupling of the mode to its surrounding. It is thus of interest to obtain deeper insight into the relaxation properties of an oscillator coupled to a heat bath. Such a study could then be used to quantify the notion of initially hot or cold modes and their evolution in time. In this paper, we use the simple model of a quantum harmonic oscillator coupled bilinearly to a dissipative heat bath to gain insight into the vibrational relaxation. In principle, this model is soluble analytically.22 Haake and Reibold derived a master equation for the dynamics of a dissipative harmonic oscillator23 which was rederived a few times by other groups.24, 25 There were also attempts to analyze the conditions for the existence of an exact Liouville operator for the dissipative harmonic oscillator.26 Ford and O’Connell derived an analytic solution for the master equation.27 Pollard and Friesner studied the relaxation of a dissipative harmonic oscillator using their generalized Redfield equation.28 Most of these studies assume factorized initial conditions for the oscillator and the bath. More recently, the quantum dynamics of the dissipative harmonic oscillator were studied in the presence of partial as well as unfactorized initial conditions.29, 30 However, to the best of our knowledge, none of the analytically exact solutions were used to study the time evolution of the vibrational distribution of the oscillator as a function of its initial state.
141, 234509-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-2
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
The purpose of the present paper is two fold. We show that indeed, the generic behavior of the dissipative oscillator conforms qualitatively with the experimental observations of Kovalenko et al.3 At the same time we show how these properties may be readily obtained in the continuum limit, using the normal mode representation of the dissipative harmonic oscillator.30–33 These results are organized as follows. The theoretical framework is discussed in Sec. II. Numerical results exemplifying the various relaxation routes are presented in Sec. III for an oscillator whose frequency mimics the 1570 cm−1 −C=C− stretching mode of trans-Stilbene in the first (S1 ) excited electronic state. We find that at room temperature, when the initial energy of the oscillator is substantially higher than the equilibrium thermal energy, the variance of the energy decreases as the system relaxes towards equilibrium. When the initial state is cold relative to the bath, the variance increases with time, in agreement with the experimental observations. We end in Sec. IV with a discussion of the limitations of the present harmonic model and methods to further extend it, making it more realistic, for example by taking into account the broadening of the intensity of individual lines through fluctuations and the anharmonicity of the vibrational mode.
or in more formal terms, the initial state remains uncoupled and the assumption of factorized initial conditions is valid. Thus, the system’s initial state is taken to be a Gaussian 1/4 0 q|(q0 , p0 ) = π¯ ip (q − q0 )2 + 0 (q − q0 ) × exp − 0 2¯ ¯ (2.7) whose frequency will be taken as the ground electronic state vibrational frequency 0 . The bath is thermally distributed so that the initial (factorized) density is ρˆ0 = ||
exp(−β Hˆ B ) . T r[exp(−β Hˆ B )]
(2.8)
To obtain the time dependent average energy and variance we let the initial density evolve up to time t under the action of the full Hamiltonian. The time evolved reduced density is then by definition, i ˆ i ˆ Ht . (2.9) ρˆt = T rB exp − H t ρˆ0 exp ¯ ¯
II. THEORY
We consider a photoexcitation process and use a dissipative Hamiltonian to describe the excited electronic state ⎡ ⎛ 2 ⎤⎞ N cj 1 2 2 2 2 2 ⎣pˆ j + ωj xˆj − Hˆ = ⎝pˆ q + 1 qˆ + qˆ ⎦⎠ , 2 ωj2 j =1 (2.1) where 1 is the frequency of the harmonic mode of the system in the excited electronic state. All momentum and position operators are mass weighted. The bath Hamiltonian is defined as 1 2 pˆ + ωj2 xˆj2 Hˆ B = 2 j =1 j
1 2 pˆ yj + λ2j yˆj2 . Hˆ = Hˆ NM = 2 j =0 N
(2.3)
|y = Uˆ |q, x.
(2.4)
This implies that the system and bath coordinates may be expressed in terms of the normal modes as (henceforth we suppressed the hat notation for the operators)
The bath is characterized by a spectral density N 2 π cj [δ(ωj − ω) − δ(ωj + ω)] J (ω) = 2 j =0 ωj
(2.10)
(2.2)
and the system harmonic oscillator Hamiltonian is
Hˆ 0 = Hˆ S + Hˆ B .
The dissipative Hamiltonian (Eq. (2.1)) is a quadratic form and so may be put into normal mode form by diagonalization,
Here pˆ y is the mass weighted momentum operator conjuj gate to the mass weighted normal mode coordinate yˆj of the jth normal mode and the λj ’s are the normal mode frequencies. The summation in the normal mode form is over N + 1 modes, since it includes also the system degree of freedom. The normal mode transformation is effected using the orthogonal transformation operator Uˆ such that23, 34, 35
N
1 Hˆ S = pˆ q2 + 21 qˆ 2 . 2 The uncoupled Hamiltonian is written as
A. The normal mode transformation
q=
(2.5)
which is related to the time dependent friction function γ (t), ∞ J (ω) = ω dt cos(ωt)γ (t). (2.6) 0
In the experiments the thermally equilibrated molecule in its ground electronic state is initially photo-excited to the electronically excited S1 state. In this fast excitation, it is reasonable to assume that initially the bath remains unperturbed,
N
(2.11)
uj 0 yj ,
(2.12)
uj k yj , k = 1, . . . , N,
(2.13)
j =0
xk =
N j =0
and conversely yj = uj 0 q +
N
uj k xk .
(2.14)
k=1
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-3
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
Since the normal modes are harmonic, their quantum time evolution is identical to the classical, py j yj (t) = yj cos(λj t) + sin(λj t), (2.15) λj py (t) = −λj yj sin(λj t) + py cos(λj t). j
j
(2.16)
Following the notation of Ref. 30 and using Eqs. (2.12), (2.13), (2.15) and (2.16) one finds that the time evolution of the system coordinate and momentum is given by N
qt = H˙0 (t)q + H0 (t)p +
H˙ l (t)xl +
l=1
N
pt = H¨0 (t)q + H˙0 (t)p +
H¨ l (t)xl +
l=1
N
N
uj 0 uj l
H˙ l (t)pl , (2.18)
l=1
sin(λj t) λj
j =0
N 2 π uj 0 [δ(λj − λ) − δ(λj + λ)] 2 j =0 λj
(2.20)
which is related to the spectral density of Eqs. (2.5) and (2.6) by33 J (λ) I (λ) = . 2 2 2 1 − λ − λI mγˆ (iλ) + J 2 (λ)
(2.21)
B. Evolution in time
where the functions Hl (t) are Hl (t) =
I (λ) =
Hl (t)pl ,
l=1
(2.17) N
To recast the normal mode representation in the continuum limit one defines a spectral density of normal modes32, 33
, l = 0, 1, . . . , N.
(2.19)
The time dependent average vibrational energy is by definition E(t) = T rS [Hˆ S ρˆt ] i i ˆ H t Hˆ S exp − Hˆ t ρˆ0 . (2.22) = T r exp ¯ ¯ Using the Wigner formulation36 this may be rewritten in phase space as
∞ ∞ N dpj dxj pt2 + 21 qt2 1 tanh(¯βωj /2) dpdq π ¯ −∞ π¯ 2 −∞ j =1 ⎛ ⎞ N 2 2 tanh(¯βω /2) (p − p (q − q ) ) j 0 0 − · exp ⎝− 0 − pj2 + ωj2 qj2 ⎠ , ¯ ¯0 ¯ω j j =1
E(t) =
where due to the harmonic nature of the model, the time evolution of the coordinates and momenta is given in Eqs. (2.17) and (2.18). The averaging over the bath coordinates is accomplished by introducing two delta function and their Fourier representations, ∞ dzδ(z − qt ), (2.24) 1= 1=
−∞ ∞
−∞
dpz δ(pz − pt )
˙ 2 + 2 G(t)2 G(t) 1 + 21 Q1 (t) + Q2 (t) (2.26) 2 such that the asymptotic value of the average energy is E(t) =
¯1 (Gc + Gs ). 4
(2.27)
Here, we used the notation Gc = 1
N
u2j 0
j =0
λj tanh(¯βλj /2)
,
(2.28)
N λj u2j 0 1 , 1 j =0 tanh(¯βλj /2)
G(t) = H˙0 (t)q0 + H0 (t)p0 , Q1 (t) =
(2.25)
into Eq. (2.23). With some further manipulation, one finds that the time dependent average energy of the oscillator is
E(t → ∞) =
Gs =
Q2 (t) =
¯ ˙ 2 H0 (t) + 21 H0 (t)2 40
N ¯ H˙ l (t)2 + ωl2 Hl (t)2 + , 4ωl tanh(¯βωl /2) l=1 ¯ ¨ 2 H0 (t) + 21 H˙0 (t)2 40
N ¯ H¨ l (t)2 + ωl2 H˙ l (t)2 , + 4ωl tanh(¯βωl /2) l=1
(2.23)
(2.29) (2.30)
(2.31)
(2.32)
1 ˙ D1 (t) = Q1 (t)Q2 (t) − Q (t)2 . (2.33) 4 1 All of these results can be expressed in terms of the spectral density of the normal modes and thus in the continuum limit. The derivation and the resulting expressions are provided in the Appendix. The continuum limit expressions obtained in
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-4
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
the Appendix for the functions Q1 (t) and Q2 (t) have not been derived previously. The time dependent width of the vibrational energy distribution of the system oscillator is obtained following similar manipulations, E(t)2 = E 2 (t) − E(t)2 = 241 Q1 (t)[Q1 (t) + G(t)2 ] ˙ 2] + 2Q2 (t)[Q2 (t) + G(t) ˙ ¯2 21 Q1 (t) 2 ˙ ˙ + G(t)G(t) − . + 21 Q1 (t) 2 4 (2.34)
1 S(ω, t) = 2π
∞ −∞
dt1
At equilibrium ¯2 21 E (t → ∞) − E(t → ∞) = 4 2
2
G2c + G2s −1 . 2 (2.35)
The functions appearing in this expression and their continuum limit form are also detailed and derived in the Appendix. The energy population spectrum of the density at time t is defined as Hˆ S (2.36) S (ω, t) = T rS δ ω − ρˆt . ¯ As shown in the Appendix, this may be expressed as
exp(iωt1 − (t)) 1/2 ,
2 2 cos(1 t1 /2) D1 (t) Ezz Ep p − Ezp z
(2.37)
z z
where the exponent is (t) =
tan(1 t1 /2) 2 4D1 (t) Ezz Ep p − Ezp z z z ˙ 2 (t) + 2 G2 (t) ˙ 2 (t) − Q ˙ 1 (t)G(t)G(t)] ˙ i G 4 tan(1 t1 /2)[Q2 (t)G2 (t) + Q1 (t)G 1 · − ¯1 ¯2
(2.38)
and 2 Ezz Ep p − Ezp = z z z
4 tan2 (1 t1 /2) i tan(1 t1 /2) 1 + Q2 (t) + 21 Q1 (t) − . 4D1 (t) ¯1 D1 (t) ¯2
III. NUMERICAL APPLICATION
The numerical values to be used will mimic as much as possible the dynamics of the −C=C− stretching mode of trans-Stilbene in the S1 excited electronic state. The mode frequency of the excited state is 1570 cm−1 corresponding to a time scale for the vibrational motion of roughly 20 fs. The characteristic time scale for vibrational motions of a molecule immersed in a liquid are of the same order of magnitude as the time scales of the vibrational motions of the liquid. It is therefore necessary to use a friction function with appropriate memory. This is also necessary, since, as is well known,37 purely Ohmic friction leads to divergence in the average kinetic energy. An additional time scale to be considered is the resolution time in the FSRS experiment which is of the order of 100 fs.3 With this in mind, we consider an exponential memory friction function defined as γ (t) =
γ exp (−t/τ ) , τ
(3.1)
where γ is the friction strength parameter and τ is the decay (memory) time. Its Laplace transform is γˆ (s) =
γ . sτ + 1
(3.2)
The spectral density of the bath is then (see Eq. (2.6)) γω J (ω) = . 1 + ω2 τ 2
(2.39)
(3.3)
We have chosen the experimentally measured frequencies3, 38 of the −C=C− mode which are 0 = 1596 cm−1 in the ground electronic state and 1 = 1570 cm−1 in the S1 state. The memory time τ (10 fs) and the friction coefficient γ (40 cm−1 ) were adjusted so as to obtain relaxation time scales that are in qualitative agreement with the experimentally measured ones (a few picoseconds). This implies that in reduced units 1 τ = 2.93 and γ /1 = 0.0254. At room temperature (T = 300 K) ¯1 β = 7.43. The range of initial energies in excess of the (0-0) transition frequency between the ground vibrational states of the S0 and S1 electronic states probed in the experiment was between 0–7000 cm−1 . Accordingly, the average energy of the initial Gaussian wavepacket (Eq. (2.7)) was adjusted by varying the central momentum p0 . The spatial center of the initial wavepacket was chosen as q0 = 0, in qualitative accordance with a vertical transition from the ground state. There is one significant difference between the present study and the experiment. Here, we are following the evolution in time of the whole vibrational population of the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-5
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
10000
4000
-1
(cm )
T=294K
v=5 v=4
(b) 3000
-1
6000
v=3 v=2
(a)
(cm )
8000
v=5 v=4
v=1 4000
v=3 v=2
T=294K 2000
v=1
1000
2000 0 0
2
4
6
8
10
12
14
0
2
time/ps
4
6
8
10
12
14
time/ps
FIG. 1. The relaxation of the average energy (E(t), panel (a) and the width ( E(t), panel (b) as a function of the initial average energy of the wavepacket (expressed in terms of the eigenstate energies of the system oscillator) at a bath temperature of 294 K.
oscillator, the FSRS experiment measures the instantaneous evolution in time of a single mode of the oscillator. This includes in it the direct measurement of the fluctuations induced by the bath (as further amplified in Sec. IV, see below). This difference implies, that in the present study, the standard deviation is large (a few thousand cm−1 ) while the experimental width of the measured FSRS spectra is of the order of 20–50 cm−1 . In the theory presented here the linewidth of any state is by definition zero. The effect of the initial energy on the relaxation path is illustrated in Fig. 1. The different “ν” values correspond to the different initial average energies of the wavepacket. The notation ν = n (n = 0 − 5) implies that the initial momentum of the wavepacket is chosen such that the initial energy of the wavepacket is equal to the energy of the nth eigenstate of the harmonic oscillator of the excited state. As can be seen from Figure 1, the average energy as well as the width E(t) decrease monotonically with time for all initial states whose initial average energy is larger than the thermal equilibrium limit. A more detailed picture of the relaxation is obtained by considering the time evolution of the vibrational population in the excited state. In Fig. 2 we plot this evolution for an initial energy corresponding to ν = 4. As may be seen in the figure, initially the population is Gaussian distributed among the states of the oscillator. In time the population of the excited states disappears until ultimately, at the bath temperature considered (T = 294 K), virtually all the population ends up in the ground state. It is also of interest to see the evolution of the energy when initially the mode is cold relative to its equilibrium limit. The effect is illustrated in Fig. 3 where the initial state is at the energy corresponding to ν = 1. The evolution of the average energy and the standard deviation is shown for three temperatures. At T = 300 and 600 K, the initial energy is higher than the equilibrium value and, as before, the energy and the standard deviation decay with time until reaching equilibrium. For T = 4000 K, the opposite occurs. Initially, the oscillator is cold and then it heats up due to the bath. The average energy and the width increase with time. This is qualitatively the same behavior found in the FSRS experiments. It is also of interest to study the dependence of the relaxation on the frictional strength. The time dependence
of the average energy and the width are plotted for strong (320 cm−1 ) moderate (160 cm−1 , 80 cm−1 ) and weak (40 cm−1 ) friction in Figure 4. As evident from the figure, the qualitative behavior of E(t) as well as E(t) is the same for all friction strengths. Increasing the magnitude of the friction leads, as expected, to a faster relaxation time. This result is qualitatively similar to the change in relaxation times observed in the FSRS experiments by changing the viscosity of the solvent.5 IV. DISCUSSION
Continuum limit expressions have been derived for the energy and population relaxation of a quantum damped harmonic oscillator using the normal mode representation. Application to a model for the −C=C− stretching mode of trans-Stilbene shows that when the initial preparation of the oscillator is such that its energy and width are larger than the thermal equilibrium limit, then the two decrease with time, until reaching the equilibrium values. Conversely, when the initial preparation is that of a cold oscillator relative to the thermal limit, the average energy and width increase with time until reaching the thermal limit. These results are in qualitative agreement with recent FSRS experiments, in which it was noted that at low photo-excitation energies, the width of the vibrational mode increases with time, while when the photoexcitation energy is high, the width decreases. In this sense, the present model supports the interpretation of the experimental time dependence of the width as a measure of the initial temperature of the nascent mode population, following photo-excitation. There is though one important qualitative difference between the model studied in this paper and the FSRS experiments. In the model, we follow the evolution of the vibrational population of all states of a given mode. The width of the distribution for a given mode is zero. The evolution is of the magnitude of the population in each state. In the experiment, one follows a single vibrational state which has a width and it is this width which is followed in time. Thus the next step is to consider a more refined model, which would also include a broadening of the spectral lines. In the present model we followed the evolution of the population defined by the system Hamiltonian as in Eq. (2.3). This Hamiltonian does not have a stochastic component to it, therefore the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-6
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
0.16
t=0 ps T=294 K
v=4
0.24 S(ω)
0.12 S(ω)
0.32
v=4
0.08 0.04
t=1 ps T=294 K
0.16 0.08
0
0 0
5
10
15 3
20
25
30
0
5
10
-1
3
ω(10 cm )
25
30
(b)
0.5
0.8
v=4
0.4
v=4
0.6 t=2 ps T=294 K
0.3
S(ω)
S(ω)
20 -1
ω(10 cm )
(a)
0.2
t=3 ps T=294 K
0.4 0.2
0.1 0
0 0
5
10
15 3
20
25
30
0
5
10
-1
15 3
20
25
30
-1
ω(10 cm )
ω(10 cm )
(c)
(d)
1
v=4
0.8 S(ω)
15
t=15 ps T=294 K
0.6 0.4 0.2 0 0
5
10
15 3
20
25
30
-1
ω(10 cm )
(e) FIG. 2. The time dependent vibrational population of the dissipative harmonic oscillator. Note the shift of the population towards the lower energy and the reduction of the variance of the distribution as equilibrium is reached. The width of each peak is a result of the finite in time numerical Fourier transform needed to compute the distributions (see Eq. (2.37)).
6000
5000 T=300K T=600K
(a)
4000
T=4000K
3000 2000
T=4000K
-1
4000
T=300K T=600K
(b) (cm )
(cm-1)
5000
v=1
3000 v=1
2000 1000
1000 0
0 0
2
4
6 8 10 time/ps
12
14
0
2
4
6 8 10 time/ps
12
14
FIG. 3. The evolution of the average energy (E(t), panel (a) and the width ( E(t), panel (b) are shown for three different bath temperatures and an initial energy corresponding to ν = 1. Note the increase of the energy and the width when initially the oscillator is cold relative to the bath temperature.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
P. Kumar and E. Pollak
8000
-1
γ=160 cm γ=320 cm-1 T=294K, v=4
4000
γ=40 cm-1 γ=80 cm-1
(b) 3000
-1
γ=160 cm γ=320 cm-1
-1
(cm-1)
4000
γ=40 cm-1 γ=80 cm-1
(a) 6000
J. Chem. Phys. 141, 234509 (2014)
(cm )
234509-7
2000
2000 T=294K, v=4 1000
0 0
2
4
6 8 10 time/ps
12
14
0
2
4
6 8 10 time/ps
12
14
FIG. 4. The relaxation of the average energy (E(t), panel (a) and the width ( E(t), panel (b) as a function of time for different value of the damping constant. In all cases, the bath temperature is fixed at T = 294 K. Initially, the oscillator is in the ν = 4 state and so is hot relative to the bath.
vibrational population consists of lines. One may include though a stochastic element in the definition of the system Hamiltonian by defining it, for example, as ⎞ ⎛ N N 2 c 1 1 j 2⎠ 2 ˆ ˆ q (4.1) Hˆ S = pˆ q2 + ⎝ + − cj xˆj q. 1 2 2 j =1 ωj2 j =1 One would then replace this Hamiltonian in the definitions of the time dependent evolution (see, e.g., Eq. (2.22)) and then one can follow the evolution of the distribution for a single state of the mode. Continuum limit analytic expressions may also be derived in this case although the algebra does become more involved. Finally, the FSRS experiments also measure a shift in time of the maximum of the distribution for a single mode. This shift has been interpreted in terms of the anharmonicity of the mode3, 4 which of course is not included in the present harmonic model. It should be possible to introduce anharmonicity and treat it with a perturbation theory or using numerically exact methods such as the multi-layer-multiconfiguration-time-dependent-hatree methodology with a dissipative bath.39, 40 These issues are topics for future work. ACKNOWLEDGMENTS
This work was supported by grants from the Israeli Science Foundation, the German Israel Foundation for basic research and the Minerva Foundation, Munich. P.K. acknowledges a fellowship of the Israel Council for Higher Education for Excellent Postdoctoral Research Fellows from India and China. APPENDIX: ANALYTIC DERIVATIONS
and original bath frequencies as31 uj l =
λ2j
cl u , l = 1, . . . , N − ωl2 j 0
(A1)
therefore N u2j 0 sin(λj t) c , l = 1, . . . , N. Hl (t) = λ λ2j − ωl2 l j =0 j
(A2)
To bring this expression to continuum limit form we note that the normal mode eigenvalue equation is31 λ2j = 21 1 +
N l=1
cl2 ωl2 ωl2 − λ2j
−1 (A3)
so that N l=1
c2 c2 l 2 l 2 = 21 − λ2j + 2 ω ωl − λj l l=1 N
(A4)
and it is this form which can then be used to remove any difficulties which might seemingly come from the denominator of Eq. (A2). We also note the relationship31 N l=1
cl2 1 2 = 2 − 1. 2 2 uj 0 ωl − λj
(A5)
The spectral density of normal modes and its continuum limit have been given in Eqs. (2.20) and (2.21). This then implies, for example, that (See Eqs. (2.28) and (2.29)): 21 ∞ I (λ) , (A6) dλ Gc = π 0 tanh(¯βλ/2)
1. Continuum limit
To express the results in the continuum limit, it is necessary to derive continuum limit expressions for the functions Q1 (t) and Q2 (t) defined in Eqs. (2.31) and (2.32). These are the central new results of this Appendix. This in turn implies obtaining continuum limit expressions for quantities such as Hl (t) (see Eq. (2.19)). For this purpose we note that the ujl off diagonal matrix elements of the normal mode transformation matrix can be expressed in terms of the coupling coefficients
Gs =
2 π 1
∞
dλ 0
λ2 I (λ) . tanh(¯βλ/2)
(A7)
For Ohmic friction, for large λ, I(λ) ∼ λ−3 so that in this case Gs would diverge, but introducing an exponential cutoff into the spectral density J(λ) eliminates it. In addition to the functions Hl (t) defined in Eq. (2.19), we use the following notation and using the normal mode spectral
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-8
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
and
density, note the continuum limit forms: H¯ (t) =
N
u2j 0 cos(λj t)
j =0
λ2j H0 (t) =
2 π
2 H˙0 (t) = π
2 = π
0
∞
2 H¨0 (t) = − π
I (λ) cos(λt) , (A8) dλ λ
∞
dλλ2 I (λ) sin(λt).
(A11)
0
We also note the expansion of the coth function41
∞
dλI (λ) sin(λt),
coth
(A9)
0
∞ 2 1 π ¯βωl 4ω = + l 2 2 2π ¯βωl ¯β k=1 ν k + ωl2
(A12)
and the definition of the Matsubara frequency
∞
dλλI (λ) cos(λt),
(A10)
ν=
0
2π . ¯β
(A13)
It is then a matter of some extensive algebra to show that (for l = 0),
N u2j 0 H˙ l (t)2 + ωl2 Hl (t)2 = ¯βλ j ωl tanh(¯βωl /2) j =0 λ tanh j
−
2
N u2j 0 sin(λj t) 2 2 ¯2 1 H (t) − H02 (t) − 2H0 (t) ¯βλ j ¯β j =0 tanh 2
⎛ ⎞2 ∞ N 2 u cos(λ t)
4 j j0 ⎠ 2 + ν 2 k 2 + νk γˆ (νk) ⎝ − 2 k 2 + λ2 ¯β k=1 1 ν j j =0 ⎛ ⎞2 ∞ N 2 4 ⎝ uj 0 λj sin(λj t) ⎠ + ¯β k=1 j =0 ν 2 k 2 + λ2j ⎛ ⎞2 ∞ N 2 u sin(λ t)
4 j j0 ⎠ . ν 2 k 2 21 +νk γˆ (νk) ⎝ + 2 k 2 +λ2 ¯β k=1 ν λ j j =0 j
(A14)
Here, we used the fact that the discretized form of the friction function is32, 33 γ (t) =
N cl2 cos(ωl t) ωl2 l=1
(A15)
so that its Laplace transform is
∞
γˆ (s) =
dt exp (−st) γ (t) =
0
N cl2 s 2 2 ωl s + ωl2 l=1
(A16)
or equivalently N l=1
c2 cl2 l = − s γˆ (s) . 2 2 s 2 + ωl ω l l=1 N
(A17)
With these preliminaries one finds that the continuum form of Q1 (t) (Eq. (2.31)) is ∞ ∞ 2 2 4 ¯ 2 4 dλI (λ)λ cos(λt) + 1 2 dλI (λ) sin(λt) Q1 (t) = 40 π 2 π 0 0 ⎡ ⎛ ⎞⎤ ∞ ∞ ∞ ¯ 2 I (λ) λI (λ) sin(λt) ⎠⎦ 8 − 2⎣ + dλ dλI (λ) sin(λt) ⎝ dλ ¯βλ 4 π 0 π 0 0 tanh tanh ¯βλ 2
2
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-9
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
∞ 2 ∞ 8 I (λ) cos(λt) 2 2 − 2 dλ − dλI (λ) sin(λt) 1 π ¯β λ 0 0 −
∞ ∞ 16 2 I (λ)λ cos(λt) 2 2 2 (νk) + ν k + νk γ ˆ dλ π 2 ¯β k=1 1 (ν 2 k 2 + λ2 ) 0
∞ ∞ 16 2 2 2 I (λ) sin(λt) 2 + 2 ν k 1 + νk γˆ (νk) dλ 2 2 π ¯β k=1 (ν k + λ2 ) 0 2 ∞ ∞ 16 I (λ)λ2 sin(λt) + 2 dλ 2 2 π ¯β k=1 0 (ν k + λ2 )
(A18)
and its time derivative follows immediately. Following the same steps one then finds that N ¨ N u2j 0 λj Hl (t)2 + ωl2 H˙ l (t)2 = ¯βλ j ωl tanh(¯βωl /2) l=1 j =0 tanh 2
⎡ −
⎛
N
2 ⎢ 2⎝ ⎣1 ¯β j =0
−2
N
u2j 0
⎞2 sin(λj t) λj
u2j 0 cos(λj t)
j =0
⎛
⎠ −⎝
N
⎞2 ⎤ ⎥ u2j 0 cos(λj t)⎠ ⎦
j =0
N λj u2j 0 cos(λj t) ¯βλ j j =0 tanh 2
⎛ ⎞2 ∞ N 2 u λ sin(λ t) 4 2 j j j 0 ⎠ − + ν 2 k 2 + νk γˆ (νk) ⎝ 2 k 2 + λ2 ¯β k=1 1 ν j j =0 ⎛ ⎞2 ∞ N 2 2 u λ cos(λ t) 4 j j0 j ⎝ ⎠ + ¯β k=1 j =0 ν 2 k 2 + λ2j ⎛ ⎞2 ∞ N u2j 0 cos(λj t) 4 2 2 2 ⎠ ν k 1 + νk γˆ (νk) ⎝ + ¯β k=1 ν 2 k 2 + λ2j j =0
(A19)
so that the continuum form of Q2 (t) (Eq. (2.32)) is ∞ ∞ 2 2 421 ¯ 4 2 dλI (λ)λ sin(λt) + 2 dλI (λ)λ cos(λt) Q2 (t) = − 2 40 π π 0 0 ⎡ ⎞⎤ ⎛ ∞ ∞ 2 ¯ 2 ∞ I (λ) cos(λt) I (λ)λ2 λ 8 ⎣ − 2 ⎠⎦ dλ dλI (λ)λ cos(λt) ⎝ dλ 4 π 0 π 0 0 tanh ¯βλ tanh ¯βλ −
8 π 2 ¯β
2
2
∞
21
dλI (λ) sin(λt) 0
−
2
∞
2
dλI (λ)λ cos(λt) 0
2 ∞ ∞ I (λ)λ2 sin(λt) 16 2 2 2 + ν k + νk γˆ (νk) dλ 2 2 − 2 π ¯β k=1 1 (ν k + λ2 ) 0 ∞ ∞ 16 2 2 2 I (λ)λ cos(λt) 2 + 2 ν k 1 + νk γˆ (νk) dλ 2 2 π ¯β k=1 (ν k + λ2 ) 0 2 ∞ ∞ 16 I (λ)λ3 cos(λt) + 2 dλ π ¯β k=1 0 (ν 2 k 2 + λ2 )
(A20)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-10
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
2. The quantum spectrum
To obtain the necessary expression for the quantum spectrum (Eqs. (2.36) and (2.37)), we use the Wigner formulation. The Wigner representation for the thermal distribution of the jth bath oscillator is tanh(¯βωj /2) 2 1 1 exp − pj + ωj2 qj2 (A21) exp(−βHj )W = 2π ¯ cosh(¯βωj /2) ¯ωj and T r[exp(−βHj )] =
1 . 2 sinh(¯βωj /2)
(A22)
The Wigner representation for the initial wavefunction is ||W
1 0 (q − q0 )2 (p − p0 )2 = . exp − − π¯ ¯ ¯0
(A23)
The Wigner representation for the system propagator is exp(−iHS t/¯)W The initial spectrum is then S(ω, t = 0) = 2 1 = 2π
1 i tan(1 t/2) 2 1 2 2 exp − = p + 1 q . 2π ¯ cos(1 t/2) ¯1
Hˆ dpdqδ ω − S ¯
W
∞
(p − p0 )2 0 (q − q0 )2 − exp − ¯ ¯0
1 dt1 exp(iωt1 ) cos(1 t1 /2) −∞
i tan(1 t1 /2) × exp − ¯
(A24)
0 1 1 + i0 tan(1 t1 /2)
1 0 + i1 tan(1 t1 /2)
0 p02 + 21 q02 + i1 tan(1 t/2) p02 + 20 q02 . [1 + i0 tan(1 t1 /2)][(0 + i1 tan(1 t1 /2))]
(A25)
The equilibrium spectrum is
Hˆ S ˆ exp(−β H ) T r[exp(−β Hˆ )] Sβ (ω) = T r δ ω − ¯ W ∞ exp(iωt1 ) 1 dt = 1/2 , 1 t1 1 t1 t t 2π −∞ 1 cos 2 + i sin 2 Gc cos 21 1 + i sin 21 1 Gs
(A26)
where Gc and Gs are defined in Eqs. (2.28) and (2.29). For the harmonic oscillator system the quantum time evolution of the Wigner distribution is the same as the classical so that ∞ i Hˆ i Hˆ 1 1 1 HS exp − = dt exp(iωt1 ) exp δ ω− ¯ ¯ ¯ 2π −∞ 1 2π ¯ cos(1 t1 /2) W i tan(1 t1 /2) 2 2 2 · exp − pt + 1 qt (A27) ¯1 and the expression for the spectrum is ∞ ∞ ∞ N dpj dxj 1 tanh(¯βωj /2) S (ω, t) = 2 dt1 dpdq 2π ¯ −∞ π¯ −∞ −∞ j =1 i tan(1 t1 /2) 2 exp(iωt1 ) 2 2 exp − (pt + 1 qt ) · cos(1 t1 /2) ¯1 ⎞ ⎛ N 2 2 tanh(¯βω /2) (p − p (q − q ) ) j 0 0 − pj2 + ωj2 qj2 ⎠ , − · exp ⎝− 0 ¯ ¯0 ¯ω j j =1
(A28)
where the time dependence is given explicitly in Eqs. (2.17) and (2.18). This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-11
P. Kumar and E. Pollak
To carry out the integrations in Eq. (A28) we introduce the Fourier representation of the two delta functions as in Eqs. (2.24) and (2.25). One then remains with Gaussian integrals which are readily carried out and the result is given in Eq. (2.37). 1 K.
A. Briney, L. Herman, D. S. Boucher, A. D. Dunkelberger, and F. F. Crim, J. Phys. Chem. A 114, 9788 (2010). 2 J. D. Leonard, Jr., and T. L. Gustafson, J. Phys. Chem. A 105, 1724 (2001). 3 S. Kovalenko, A. Dobryakov, E. Pollak, and N. Ernsting, J. Chem. Phys. 139, 011101 (2013). 4 K. Iwata and H. O. Hamaguchi, Chem. Phys. Lett. 196, 462 (1992). 5 K. Iwata, R. Ozawa, and H. O. Hamaguchi, J. Phys. Chem. A 106, 3614 (2002). 6 K. A. Iwata, K. B. Yoshida, Y. B. Takada, and H. O. Hamaguchi, Chem. Lett. 36, 504 (2007). 7 K. Iwata and H. O. Hamaguchi, J. Mol. Liq. 65–66, 417 (1995). 8 V. Deckert, K. Iwata, and H. O. Hamaguchi, J. Photochem. Photobiol. A 102, 35 (1996). 9 K. Iwata and H. O. Hamaguchi, J. Phys. Chem. A 101, 632 (1997). 10 C. Warmuth, F. Milota, H. F. Kauffmann, H. Wadi, and E. Pollak, J. Chem. Phys. 112, 3938 (2000). 11 D. Gegiou, K. A. Muszkat, and E. Fischer, J. Am. Chem. Soc. 90, 12 (1968). 12 R. M. Hochstrasser, Pure Appl. Chem. 52, 2683 (1980). 13 S. H. Courtney and G. R. Fleming, J. Chem. Phys. 83, 215 (1985). 14 M. W. Balk and G. R. Fleming, J. Phys. Chem. 90, 3975 (1986). 15 S. H. Courtney, M. W. Balk, L. A. Philips, S. P. Webb, D. Yang, D. H. Levy, and G. R. Fleming, J. Chem. Phys. 89, 6697 (1988). 16 S. A. Kovalenko, R. Schanz, H. Hennig, and N. P. Ernsting, J. Chem. Phys. 115, 3256 (2001).
J. Chem. Phys. 141, 234509 (2014) 17 S. A. Kovalenko, A. L. Dobryakov, I. Ioffe, and N. P. Ernsting, Chem. Phys.
Lett. 493, 255 (2010). O. Hamaguchi, Mol. Phys. 89, 463 (1996). 19 H. Wadi and E. Pollak, J. Chem. Phys. 110, 11890 (1999). 20 G. Gershinsky and E. Pollak, J. Chem. Phys. 107, 812 (1997). 21 J. Tatchen and E. Pollak, J. Chem. Phys. 128, 164303 (2008). 22 H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 168, 115 (1988). 23 F. Haake and R. Reibold, Phys. Rev. A 32, 2462 (1985). 24 B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D 45, 2843 (1992). 25 J. J. Halliwell and T. Yu, Phys. Rev. D 53, 2012 (1996). 26 R. Karrlein and H. Grabert, Phys. Rev. E 55, 153 (1997). 27 G. W. Ford and R. F. O’Connell, Phys. Rev. D 64, 105020 (2001). 28 W. T. Pollard and R. A. Friesner, J. Chem. Phys. 100, 5054 (1994). 29 A. Pereverzev, Phys. Rev. E 68, 026111 (2003). 30 E. Pollak, J. S. Shao, and D. H. Zhang, Phys. Rev. E 77, 021107 (2008). 31 A. M. Levine, M. Shapiro, and E. Pollak, J. Chem. Phys. 88, 1959 (1988). 32 I. Rips and E. Pollak, Phys. Rev. A 41, 5366 (1990). 33 E. Pollak, in Dynamics of Molecules and Chemical Reactions, edited by R. E. Wyatt, and J. Z. H. Zhang (Dekker, New York, 1996), p. 617. 34 E. Pollak, J. Chem. Phys. 85, 865 (1986). 35 E. Pollak, Phys. Rev. A 33, 4244 (1986). 36 E. Wigner, Phys. Rev. 40, 749 (1932). 37 U. Weiss, Quantum Dissipative Systems, 3rd ed. (World Scientific, 2008). 38 A. L. Dobryakov, I. Ioffe, A. A. Granovsky, N. P. Ernsting, and S. A. Kovalenko, J. Chem. Phys. 137, 244505 (2012). 39 H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003). 40 H. Wang and M. Thoss, in Multidimensional Quantum Dynamics: MCTDH Theory and Applications, edited by H.-D. Meyer, F. Gatti, and G. A. Worth (Wiley-VCH, 2009), p. 131. 41 I. S. Gradshteyn, and I. M. Ryzhik, in Table of Integrals, Series and Products, edited by A. Jeffrey (Academic Press, New York, 1965). 18 H.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
THE JOURNAL OF CHEMICAL PHYSICS 141, 234509 (2014)
Energy relaxation of a dissipative quantum oscillator Pradeep Kumar and Eli Pollak Chemical Physics Department, Weizmann Institute of Science, 76100 Rehovot, Israel
(Received 29 September 2014; accepted 28 November 2014; published online 19 December 2014) The dissipative harmonic oscillator is studied as a model for vibrational relaxation in a liquid environment. Continuum limit expressions are derived for the time-dependent average energy, average width of the population, and the vibrational population itself. The effect of the magnitude of the solute-solvent interaction, expressed in terms of a friction coefficient, solvent temperature, and initial energy of the oscillator on the relaxation has been studied. These results shed light on the recent femtosecond stimulated Raman scattering probe of the 1570 cm−1 −C=C− stretching mode of trans-Stilbene in the first (S1 ) excited electronic state. When the oscillator is initially cold with respect to the bath temperature, its average energy and width increase in time. When it is initially hot, the average energy and width decrease with time in qualitative agreement with the experimental observations. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4903809] I. INTRODUCTION
A vibrationally excited molecule will eventually relax to an equilibrium state in the presence of a thermal bath. The evolution of this relaxation depends upon many factors such as solute-solvent coupling, intramolecular-vibrationalrelaxation, temperature of the bath and the initial state of the molecule.1–10 Many experiments have been carried out to probe the time-dependent vibrational relaxation of the Stilbene molecule.3–9, 11–17 The measured data on a given vibrational band are characteristically summarized in terms of the time evolution of the instantaneous maximum of the band and its full width at half maximum (FWHM). The early experimental studies were carried out by combining ultrafast spectroscopy techniques in a supersonic jet.13–15 Hamaguchi and co-workers4–9 used time-dependent Raman spectroscopy from the S1 to the S2 electronic states, to measure the time-dependent linewidth and peak position of the Raman spectra of the 1570 cm−1 −C=C− stretching mode of trans-Stilbene. They observed a linear correlation between the time-dependent evolution of the peak position and the linewidth. They found a narrowing of the linewidth in time until reaching its equilibrium value. The experimental studies of Hamaguchi and co-workers were limited to the picosecond domain. The advent of femtosecond-stimulated-Raman-spectroscopy (FSRS) enabled Ernsting and co-workers to extend the time-dependent study of the −C=C− stretching mode of trans-Stilbene to the femtosecond domain.3, 16, 17 In a recent study, Kovalenko et al.3 found that for a short time scale (0-3 ps) the time evolution of the line shape changes dramatically as a function of the initial excitation energy. When the initial energy is high, in the sense that it may excite a few quanta of the mode, the linewidth first broadens and only later narrows, as observed by Hamaguchi and co-workers. When the initial excitation populates only the ground state then the linewidth only broadens with time until reaching its equilibrium value. There have been a few attempts to understand these experimental observations theoretically.5, 10, 18–21 Hamaguchi 0021-9606/2014/141(23)/234509/11/$30.00
proposed a “dynamic polarization model” in which the time evolution is dominated by the stochastic polarization of the mode, induced by the solvent.18 Kovalenko et al.3 suggested that the change in the dynamics, from broadening of the FWHM at low excitation energies to narrowing at high excitation energies serves as an indication of the initial temperature of the excited mode. A highly excited mode couples with many solvent degrees of freedom so that its spectrum is broad. As it cools down, the number of coupled modes decreases and the width narrows. Conversely, when the mode is cold, it couples only to a few bath modes, subsequent heating by the bath then broadens the peak. This explanation, if correct should be of general validity and should not depend too much on the nature of the coupling of the mode to its surrounding. It is thus of interest to obtain deeper insight into the relaxation properties of an oscillator coupled to a heat bath. Such a study could then be used to quantify the notion of initially hot or cold modes and their evolution in time. In this paper, we use the simple model of a quantum harmonic oscillator coupled bilinearly to a dissipative heat bath to gain insight into the vibrational relaxation. In principle, this model is soluble analytically.22 Haake and Reibold derived a master equation for the dynamics of a dissipative harmonic oscillator23 which was rederived a few times by other groups.24, 25 There were also attempts to analyze the conditions for the existence of an exact Liouville operator for the dissipative harmonic oscillator.26 Ford and O’Connell derived an analytic solution for the master equation.27 Pollard and Friesner studied the relaxation of a dissipative harmonic oscillator using their generalized Redfield equation.28 Most of these studies assume factorized initial conditions for the oscillator and the bath. More recently, the quantum dynamics of the dissipative harmonic oscillator were studied in the presence of partial as well as unfactorized initial conditions.29, 30 However, to the best of our knowledge, none of the analytically exact solutions were used to study the time evolution of the vibrational distribution of the oscillator as a function of its initial state.
141, 234509-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-2
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
The purpose of the present paper is two fold. We show that indeed, the generic behavior of the dissipative oscillator conforms qualitatively with the experimental observations of Kovalenko et al.3 At the same time we show how these properties may be readily obtained in the continuum limit, using the normal mode representation of the dissipative harmonic oscillator.30–33 These results are organized as follows. The theoretical framework is discussed in Sec. II. Numerical results exemplifying the various relaxation routes are presented in Sec. III for an oscillator whose frequency mimics the 1570 cm−1 −C=C− stretching mode of trans-Stilbene in the first (S1 ) excited electronic state. We find that at room temperature, when the initial energy of the oscillator is substantially higher than the equilibrium thermal energy, the variance of the energy decreases as the system relaxes towards equilibrium. When the initial state is cold relative to the bath, the variance increases with time, in agreement with the experimental observations. We end in Sec. IV with a discussion of the limitations of the present harmonic model and methods to further extend it, making it more realistic, for example by taking into account the broadening of the intensity of individual lines through fluctuations and the anharmonicity of the vibrational mode.
or in more formal terms, the initial state remains uncoupled and the assumption of factorized initial conditions is valid. Thus, the system’s initial state is taken to be a Gaussian 1/4 0 q|(q0 , p0 ) = π¯ ip (q − q0 )2 + 0 (q − q0 ) × exp − 0 2¯ ¯ (2.7) whose frequency will be taken as the ground electronic state vibrational frequency 0 . The bath is thermally distributed so that the initial (factorized) density is ρˆ0 = ||
exp(−β Hˆ B ) . T r[exp(−β Hˆ B )]
(2.8)
To obtain the time dependent average energy and variance we let the initial density evolve up to time t under the action of the full Hamiltonian. The time evolved reduced density is then by definition, i ˆ i ˆ Ht . (2.9) ρˆt = T rB exp − H t ρˆ0 exp ¯ ¯
II. THEORY
We consider a photoexcitation process and use a dissipative Hamiltonian to describe the excited electronic state ⎡ ⎛ 2 ⎤⎞ N cj 1 2 2 2 2 2 ⎣pˆ j + ωj xˆj − Hˆ = ⎝pˆ q + 1 qˆ + qˆ ⎦⎠ , 2 ωj2 j =1 (2.1) where 1 is the frequency of the harmonic mode of the system in the excited electronic state. All momentum and position operators are mass weighted. The bath Hamiltonian is defined as 1 2 pˆ + ωj2 xˆj2 Hˆ B = 2 j =1 j
1 2 pˆ yj + λ2j yˆj2 . Hˆ = Hˆ NM = 2 j =0 N
(2.3)
|y = Uˆ |q, x.
(2.4)
This implies that the system and bath coordinates may be expressed in terms of the normal modes as (henceforth we suppressed the hat notation for the operators)
The bath is characterized by a spectral density N 2 π cj [δ(ωj − ω) − δ(ωj + ω)] J (ω) = 2 j =0 ωj
(2.10)
(2.2)
and the system harmonic oscillator Hamiltonian is
Hˆ 0 = Hˆ S + Hˆ B .
The dissipative Hamiltonian (Eq. (2.1)) is a quadratic form and so may be put into normal mode form by diagonalization,
Here pˆ y is the mass weighted momentum operator conjuj gate to the mass weighted normal mode coordinate yˆj of the jth normal mode and the λj ’s are the normal mode frequencies. The summation in the normal mode form is over N + 1 modes, since it includes also the system degree of freedom. The normal mode transformation is effected using the orthogonal transformation operator Uˆ such that23, 34, 35
N
1 Hˆ S = pˆ q2 + 21 qˆ 2 . 2 The uncoupled Hamiltonian is written as
A. The normal mode transformation
q=
(2.5)
which is related to the time dependent friction function γ (t), ∞ J (ω) = ω dt cos(ωt)γ (t). (2.6) 0
In the experiments the thermally equilibrated molecule in its ground electronic state is initially photo-excited to the electronically excited S1 state. In this fast excitation, it is reasonable to assume that initially the bath remains unperturbed,
N
(2.11)
uj 0 yj ,
(2.12)
uj k yj , k = 1, . . . , N,
(2.13)
j =0
xk =
N j =0
and conversely yj = uj 0 q +
N
uj k xk .
(2.14)
k=1
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-3
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
Since the normal modes are harmonic, their quantum time evolution is identical to the classical, py j yj (t) = yj cos(λj t) + sin(λj t), (2.15) λj py (t) = −λj yj sin(λj t) + py cos(λj t). j
j
(2.16)
Following the notation of Ref. 30 and using Eqs. (2.12), (2.13), (2.15) and (2.16) one finds that the time evolution of the system coordinate and momentum is given by N
qt = H˙0 (t)q + H0 (t)p +
H˙ l (t)xl +
l=1
N
pt = H¨0 (t)q + H˙0 (t)p +
H¨ l (t)xl +
l=1
N
N
uj 0 uj l
H˙ l (t)pl , (2.18)
l=1
sin(λj t) λj
j =0
N 2 π uj 0 [δ(λj − λ) − δ(λj + λ)] 2 j =0 λj
(2.20)
which is related to the spectral density of Eqs. (2.5) and (2.6) by33 J (λ) I (λ) = . 2 2 2 1 − λ − λI mγˆ (iλ) + J 2 (λ)
(2.21)
B. Evolution in time
where the functions Hl (t) are Hl (t) =
I (λ) =
Hl (t)pl ,
l=1
(2.17) N
To recast the normal mode representation in the continuum limit one defines a spectral density of normal modes32, 33
, l = 0, 1, . . . , N.
(2.19)
The time dependent average vibrational energy is by definition E(t) = T rS [Hˆ S ρˆt ] i i ˆ H t Hˆ S exp − Hˆ t ρˆ0 . (2.22) = T r exp ¯ ¯ Using the Wigner formulation36 this may be rewritten in phase space as
∞ ∞ N dpj dxj pt2 + 21 qt2 1 tanh(¯βωj /2) dpdq π ¯ −∞ π¯ 2 −∞ j =1 ⎛ ⎞ N 2 2 tanh(¯βω /2) (p − p (q − q ) ) j 0 0 − · exp ⎝− 0 − pj2 + ωj2 qj2 ⎠ , ¯ ¯0 ¯ω j j =1
E(t) =
where due to the harmonic nature of the model, the time evolution of the coordinates and momenta is given in Eqs. (2.17) and (2.18). The averaging over the bath coordinates is accomplished by introducing two delta function and their Fourier representations, ∞ dzδ(z − qt ), (2.24) 1= 1=
−∞ ∞
−∞
dpz δ(pz − pt )
˙ 2 + 2 G(t)2 G(t) 1 + 21 Q1 (t) + Q2 (t) (2.26) 2 such that the asymptotic value of the average energy is E(t) =
¯1 (Gc + Gs ). 4
(2.27)
Here, we used the notation Gc = 1
N
u2j 0
j =0
λj tanh(¯βλj /2)
,
(2.28)
N λj u2j 0 1 , 1 j =0 tanh(¯βλj /2)
G(t) = H˙0 (t)q0 + H0 (t)p0 , Q1 (t) =
(2.25)
into Eq. (2.23). With some further manipulation, one finds that the time dependent average energy of the oscillator is
E(t → ∞) =
Gs =
Q2 (t) =
¯ ˙ 2 H0 (t) + 21 H0 (t)2 40
N ¯ H˙ l (t)2 + ωl2 Hl (t)2 + , 4ωl tanh(¯βωl /2) l=1 ¯ ¨ 2 H0 (t) + 21 H˙0 (t)2 40
N ¯ H¨ l (t)2 + ωl2 H˙ l (t)2 , + 4ωl tanh(¯βωl /2) l=1
(2.23)
(2.29) (2.30)
(2.31)
(2.32)
1 ˙ D1 (t) = Q1 (t)Q2 (t) − Q (t)2 . (2.33) 4 1 All of these results can be expressed in terms of the spectral density of the normal modes and thus in the continuum limit. The derivation and the resulting expressions are provided in the Appendix. The continuum limit expressions obtained in
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-4
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
the Appendix for the functions Q1 (t) and Q2 (t) have not been derived previously. The time dependent width of the vibrational energy distribution of the system oscillator is obtained following similar manipulations, E(t)2 = E 2 (t) − E(t)2 = 241 Q1 (t)[Q1 (t) + G(t)2 ] ˙ 2] + 2Q2 (t)[Q2 (t) + G(t) ˙ ¯2 21 Q1 (t) 2 ˙ ˙ + G(t)G(t) − . + 21 Q1 (t) 2 4 (2.34)
1 S(ω, t) = 2π
∞ −∞
dt1
At equilibrium ¯2 21 E (t → ∞) − E(t → ∞) = 4 2
2
G2c + G2s −1 . 2 (2.35)
The functions appearing in this expression and their continuum limit form are also detailed and derived in the Appendix. The energy population spectrum of the density at time t is defined as Hˆ S (2.36) S (ω, t) = T rS δ ω − ρˆt . ¯ As shown in the Appendix, this may be expressed as
exp(iωt1 − (t)) 1/2 ,
2 2 cos(1 t1 /2) D1 (t) Ezz Ep p − Ezp z
(2.37)
z z
where the exponent is (t) =
tan(1 t1 /2) 2 4D1 (t) Ezz Ep p − Ezp z z z ˙ 2 (t) + 2 G2 (t) ˙ 2 (t) − Q ˙ 1 (t)G(t)G(t)] ˙ i G 4 tan(1 t1 /2)[Q2 (t)G2 (t) + Q1 (t)G 1 · − ¯1 ¯2
(2.38)
and 2 Ezz Ep p − Ezp = z z z
4 tan2 (1 t1 /2) i tan(1 t1 /2) 1 + Q2 (t) + 21 Q1 (t) − . 4D1 (t) ¯1 D1 (t) ¯2
III. NUMERICAL APPLICATION
The numerical values to be used will mimic as much as possible the dynamics of the −C=C− stretching mode of trans-Stilbene in the S1 excited electronic state. The mode frequency of the excited state is 1570 cm−1 corresponding to a time scale for the vibrational motion of roughly 20 fs. The characteristic time scale for vibrational motions of a molecule immersed in a liquid are of the same order of magnitude as the time scales of the vibrational motions of the liquid. It is therefore necessary to use a friction function with appropriate memory. This is also necessary, since, as is well known,37 purely Ohmic friction leads to divergence in the average kinetic energy. An additional time scale to be considered is the resolution time in the FSRS experiment which is of the order of 100 fs.3 With this in mind, we consider an exponential memory friction function defined as γ (t) =
γ exp (−t/τ ) , τ
(3.1)
where γ is the friction strength parameter and τ is the decay (memory) time. Its Laplace transform is γˆ (s) =
γ . sτ + 1
(3.2)
The spectral density of the bath is then (see Eq. (2.6)) γω J (ω) = . 1 + ω2 τ 2
(2.39)
(3.3)
We have chosen the experimentally measured frequencies3, 38 of the −C=C− mode which are 0 = 1596 cm−1 in the ground electronic state and 1 = 1570 cm−1 in the S1 state. The memory time τ (10 fs) and the friction coefficient γ (40 cm−1 ) were adjusted so as to obtain relaxation time scales that are in qualitative agreement with the experimentally measured ones (a few picoseconds). This implies that in reduced units 1 τ = 2.93 and γ /1 = 0.0254. At room temperature (T = 300 K) ¯1 β = 7.43. The range of initial energies in excess of the (0-0) transition frequency between the ground vibrational states of the S0 and S1 electronic states probed in the experiment was between 0–7000 cm−1 . Accordingly, the average energy of the initial Gaussian wavepacket (Eq. (2.7)) was adjusted by varying the central momentum p0 . The spatial center of the initial wavepacket was chosen as q0 = 0, in qualitative accordance with a vertical transition from the ground state. There is one significant difference between the present study and the experiment. Here, we are following the evolution in time of the whole vibrational population of the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-5
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
10000
4000
-1
(cm )
T=294K
v=5 v=4
(b) 3000
-1
6000
v=3 v=2
(a)
(cm )
8000
v=5 v=4
v=1 4000
v=3 v=2
T=294K 2000
v=1
1000
2000 0 0
2
4
6
8
10
12
14
0
2
time/ps
4
6
8
10
12
14
time/ps
FIG. 1. The relaxation of the average energy (E(t), panel (a) and the width ( E(t), panel (b) as a function of the initial average energy of the wavepacket (expressed in terms of the eigenstate energies of the system oscillator) at a bath temperature of 294 K.
oscillator, the FSRS experiment measures the instantaneous evolution in time of a single mode of the oscillator. This includes in it the direct measurement of the fluctuations induced by the bath (as further amplified in Sec. IV, see below). This difference implies, that in the present study, the standard deviation is large (a few thousand cm−1 ) while the experimental width of the measured FSRS spectra is of the order of 20–50 cm−1 . In the theory presented here the linewidth of any state is by definition zero. The effect of the initial energy on the relaxation path is illustrated in Fig. 1. The different “ν” values correspond to the different initial average energies of the wavepacket. The notation ν = n (n = 0 − 5) implies that the initial momentum of the wavepacket is chosen such that the initial energy of the wavepacket is equal to the energy of the nth eigenstate of the harmonic oscillator of the excited state. As can be seen from Figure 1, the average energy as well as the width E(t) decrease monotonically with time for all initial states whose initial average energy is larger than the thermal equilibrium limit. A more detailed picture of the relaxation is obtained by considering the time evolution of the vibrational population in the excited state. In Fig. 2 we plot this evolution for an initial energy corresponding to ν = 4. As may be seen in the figure, initially the population is Gaussian distributed among the states of the oscillator. In time the population of the excited states disappears until ultimately, at the bath temperature considered (T = 294 K), virtually all the population ends up in the ground state. It is also of interest to see the evolution of the energy when initially the mode is cold relative to its equilibrium limit. The effect is illustrated in Fig. 3 where the initial state is at the energy corresponding to ν = 1. The evolution of the average energy and the standard deviation is shown for three temperatures. At T = 300 and 600 K, the initial energy is higher than the equilibrium value and, as before, the energy and the standard deviation decay with time until reaching equilibrium. For T = 4000 K, the opposite occurs. Initially, the oscillator is cold and then it heats up due to the bath. The average energy and the width increase with time. This is qualitatively the same behavior found in the FSRS experiments. It is also of interest to study the dependence of the relaxation on the frictional strength. The time dependence
of the average energy and the width are plotted for strong (320 cm−1 ) moderate (160 cm−1 , 80 cm−1 ) and weak (40 cm−1 ) friction in Figure 4. As evident from the figure, the qualitative behavior of E(t) as well as E(t) is the same for all friction strengths. Increasing the magnitude of the friction leads, as expected, to a faster relaxation time. This result is qualitatively similar to the change in relaxation times observed in the FSRS experiments by changing the viscosity of the solvent.5 IV. DISCUSSION
Continuum limit expressions have been derived for the energy and population relaxation of a quantum damped harmonic oscillator using the normal mode representation. Application to a model for the −C=C− stretching mode of trans-Stilbene shows that when the initial preparation of the oscillator is such that its energy and width are larger than the thermal equilibrium limit, then the two decrease with time, until reaching the equilibrium values. Conversely, when the initial preparation is that of a cold oscillator relative to the thermal limit, the average energy and width increase with time until reaching the thermal limit. These results are in qualitative agreement with recent FSRS experiments, in which it was noted that at low photo-excitation energies, the width of the vibrational mode increases with time, while when the photoexcitation energy is high, the width decreases. In this sense, the present model supports the interpretation of the experimental time dependence of the width as a measure of the initial temperature of the nascent mode population, following photo-excitation. There is though one important qualitative difference between the model studied in this paper and the FSRS experiments. In the model, we follow the evolution of the vibrational population of all states of a given mode. The width of the distribution for a given mode is zero. The evolution is of the magnitude of the population in each state. In the experiment, one follows a single vibrational state which has a width and it is this width which is followed in time. Thus the next step is to consider a more refined model, which would also include a broadening of the spectral lines. In the present model we followed the evolution of the population defined by the system Hamiltonian as in Eq. (2.3). This Hamiltonian does not have a stochastic component to it, therefore the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-6
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
0.16
t=0 ps T=294 K
v=4
0.24 S(ω)
0.12 S(ω)
0.32
v=4
0.08 0.04
t=1 ps T=294 K
0.16 0.08
0
0 0
5
10
15 3
20
25
30
0
5
10
-1
3
ω(10 cm )
25
30
(b)
0.5
0.8
v=4
0.4
v=4
0.6 t=2 ps T=294 K
0.3
S(ω)
S(ω)
20 -1
ω(10 cm )
(a)
0.2
t=3 ps T=294 K
0.4 0.2
0.1 0
0 0
5
10
15 3
20
25
30
0
5
10
-1
15 3
20
25
30
-1
ω(10 cm )
ω(10 cm )
(c)
(d)
1
v=4
0.8 S(ω)
15
t=15 ps T=294 K
0.6 0.4 0.2 0 0
5
10
15 3
20
25
30
-1
ω(10 cm )
(e) FIG. 2. The time dependent vibrational population of the dissipative harmonic oscillator. Note the shift of the population towards the lower energy and the reduction of the variance of the distribution as equilibrium is reached. The width of each peak is a result of the finite in time numerical Fourier transform needed to compute the distributions (see Eq. (2.37)).
6000
5000 T=300K T=600K
(a)
4000
T=4000K
3000 2000
T=4000K
-1
4000
T=300K T=600K
(b) (cm )
(cm-1)
5000
v=1
3000 v=1
2000 1000
1000 0
0 0
2
4
6 8 10 time/ps
12
14
0
2
4
6 8 10 time/ps
12
14
FIG. 3. The evolution of the average energy (E(t), panel (a) and the width ( E(t), panel (b) are shown for three different bath temperatures and an initial energy corresponding to ν = 1. Note the increase of the energy and the width when initially the oscillator is cold relative to the bath temperature.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
P. Kumar and E. Pollak
8000
-1
γ=160 cm γ=320 cm-1 T=294K, v=4
4000
γ=40 cm-1 γ=80 cm-1
(b) 3000
-1
γ=160 cm γ=320 cm-1
-1
(cm-1)
4000
γ=40 cm-1 γ=80 cm-1
(a) 6000
J. Chem. Phys. 141, 234509 (2014)
(cm )
234509-7
2000
2000 T=294K, v=4 1000
0 0
2
4
6 8 10 time/ps
12
14
0
2
4
6 8 10 time/ps
12
14
FIG. 4. The relaxation of the average energy (E(t), panel (a) and the width ( E(t), panel (b) as a function of time for different value of the damping constant. In all cases, the bath temperature is fixed at T = 294 K. Initially, the oscillator is in the ν = 4 state and so is hot relative to the bath.
vibrational population consists of lines. One may include though a stochastic element in the definition of the system Hamiltonian by defining it, for example, as ⎞ ⎛ N N 2 c 1 1 j 2⎠ 2 ˆ ˆ q (4.1) Hˆ S = pˆ q2 + ⎝ + − cj xˆj q. 1 2 2 j =1 ωj2 j =1 One would then replace this Hamiltonian in the definitions of the time dependent evolution (see, e.g., Eq. (2.22)) and then one can follow the evolution of the distribution for a single state of the mode. Continuum limit analytic expressions may also be derived in this case although the algebra does become more involved. Finally, the FSRS experiments also measure a shift in time of the maximum of the distribution for a single mode. This shift has been interpreted in terms of the anharmonicity of the mode3, 4 which of course is not included in the present harmonic model. It should be possible to introduce anharmonicity and treat it with a perturbation theory or using numerically exact methods such as the multi-layer-multiconfiguration-time-dependent-hatree methodology with a dissipative bath.39, 40 These issues are topics for future work. ACKNOWLEDGMENTS
This work was supported by grants from the Israeli Science Foundation, the German Israel Foundation for basic research and the Minerva Foundation, Munich. P.K. acknowledges a fellowship of the Israel Council for Higher Education for Excellent Postdoctoral Research Fellows from India and China. APPENDIX: ANALYTIC DERIVATIONS
and original bath frequencies as31 uj l =
λ2j
cl u , l = 1, . . . , N − ωl2 j 0
(A1)
therefore N u2j 0 sin(λj t) c , l = 1, . . . , N. Hl (t) = λ λ2j − ωl2 l j =0 j
(A2)
To bring this expression to continuum limit form we note that the normal mode eigenvalue equation is31 λ2j = 21 1 +
N l=1
cl2 ωl2 ωl2 − λ2j
−1 (A3)
so that N l=1
c2 c2 l 2 l 2 = 21 − λ2j + 2 ω ωl − λj l l=1 N
(A4)
and it is this form which can then be used to remove any difficulties which might seemingly come from the denominator of Eq. (A2). We also note the relationship31 N l=1
cl2 1 2 = 2 − 1. 2 2 uj 0 ωl − λj
(A5)
The spectral density of normal modes and its continuum limit have been given in Eqs. (2.20) and (2.21). This then implies, for example, that (See Eqs. (2.28) and (2.29)): 21 ∞ I (λ) , (A6) dλ Gc = π 0 tanh(¯βλ/2)
1. Continuum limit
To express the results in the continuum limit, it is necessary to derive continuum limit expressions for the functions Q1 (t) and Q2 (t) defined in Eqs. (2.31) and (2.32). These are the central new results of this Appendix. This in turn implies obtaining continuum limit expressions for quantities such as Hl (t) (see Eq. (2.19)). For this purpose we note that the ujl off diagonal matrix elements of the normal mode transformation matrix can be expressed in terms of the coupling coefficients
Gs =
2 π 1
∞
dλ 0
λ2 I (λ) . tanh(¯βλ/2)
(A7)
For Ohmic friction, for large λ, I(λ) ∼ λ−3 so that in this case Gs would diverge, but introducing an exponential cutoff into the spectral density J(λ) eliminates it. In addition to the functions Hl (t) defined in Eq. (2.19), we use the following notation and using the normal mode spectral
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-8
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
and
density, note the continuum limit forms: H¯ (t) =
N
u2j 0 cos(λj t)
j =0
λ2j H0 (t) =
2 π
2 H˙0 (t) = π
2 = π
0
∞
2 H¨0 (t) = − π
I (λ) cos(λt) , (A8) dλ λ
∞
dλλ2 I (λ) sin(λt).
(A11)
0
We also note the expansion of the coth function41
∞
dλI (λ) sin(λt),
coth
(A9)
0
∞ 2 1 π ¯βωl 4ω = + l 2 2 2π ¯βωl ¯β k=1 ν k + ωl2
(A12)
and the definition of the Matsubara frequency
∞
dλλI (λ) cos(λt),
(A10)
ν=
0
2π . ¯β
(A13)
It is then a matter of some extensive algebra to show that (for l = 0),
N u2j 0 H˙ l (t)2 + ωl2 Hl (t)2 = ¯βλ j ωl tanh(¯βωl /2) j =0 λ tanh j
−
2
N u2j 0 sin(λj t) 2 2 ¯2 1 H (t) − H02 (t) − 2H0 (t) ¯βλ j ¯β j =0 tanh 2
⎛ ⎞2 ∞ N 2 u cos(λ t)
4 j j0 ⎠ 2 + ν 2 k 2 + νk γˆ (νk) ⎝ − 2 k 2 + λ2 ¯β k=1 1 ν j j =0 ⎛ ⎞2 ∞ N 2 4 ⎝ uj 0 λj sin(λj t) ⎠ + ¯β k=1 j =0 ν 2 k 2 + λ2j ⎛ ⎞2 ∞ N 2 u sin(λ t)
4 j j0 ⎠ . ν 2 k 2 21 +νk γˆ (νk) ⎝ + 2 k 2 +λ2 ¯β k=1 ν λ j j =0 j
(A14)
Here, we used the fact that the discretized form of the friction function is32, 33 γ (t) =
N cl2 cos(ωl t) ωl2 l=1
(A15)
so that its Laplace transform is
∞
γˆ (s) =
dt exp (−st) γ (t) =
0
N cl2 s 2 2 ωl s + ωl2 l=1
(A16)
or equivalently N l=1
c2 cl2 l = − s γˆ (s) . 2 2 s 2 + ωl ω l l=1 N
(A17)
With these preliminaries one finds that the continuum form of Q1 (t) (Eq. (2.31)) is ∞ ∞ 2 2 4 ¯ 2 4 dλI (λ)λ cos(λt) + 1 2 dλI (λ) sin(λt) Q1 (t) = 40 π 2 π 0 0 ⎡ ⎛ ⎞⎤ ∞ ∞ ∞ ¯ 2 I (λ) λI (λ) sin(λt) ⎠⎦ 8 − 2⎣ + dλ dλI (λ) sin(λt) ⎝ dλ ¯βλ 4 π 0 π 0 0 tanh tanh ¯βλ 2
2
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-9
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
∞ 2 ∞ 8 I (λ) cos(λt) 2 2 − 2 dλ − dλI (λ) sin(λt) 1 π ¯β λ 0 0 −
∞ ∞ 16 2 I (λ)λ cos(λt) 2 2 2 (νk) + ν k + νk γ ˆ dλ π 2 ¯β k=1 1 (ν 2 k 2 + λ2 ) 0
∞ ∞ 16 2 2 2 I (λ) sin(λt) 2 + 2 ν k 1 + νk γˆ (νk) dλ 2 2 π ¯β k=1 (ν k + λ2 ) 0 2 ∞ ∞ 16 I (λ)λ2 sin(λt) + 2 dλ 2 2 π ¯β k=1 0 (ν k + λ2 )
(A18)
and its time derivative follows immediately. Following the same steps one then finds that N ¨ N u2j 0 λj Hl (t)2 + ωl2 H˙ l (t)2 = ¯βλ j ωl tanh(¯βωl /2) l=1 j =0 tanh 2
⎡ −
⎛
N
2 ⎢ 2⎝ ⎣1 ¯β j =0
−2
N
u2j 0
⎞2 sin(λj t) λj
u2j 0 cos(λj t)
j =0
⎛
⎠ −⎝
N
⎞2 ⎤ ⎥ u2j 0 cos(λj t)⎠ ⎦
j =0
N λj u2j 0 cos(λj t) ¯βλ j j =0 tanh 2
⎛ ⎞2 ∞ N 2 u λ sin(λ t) 4 2 j j j 0 ⎠ − + ν 2 k 2 + νk γˆ (νk) ⎝ 2 k 2 + λ2 ¯β k=1 1 ν j j =0 ⎛ ⎞2 ∞ N 2 2 u λ cos(λ t) 4 j j0 j ⎝ ⎠ + ¯β k=1 j =0 ν 2 k 2 + λ2j ⎛ ⎞2 ∞ N u2j 0 cos(λj t) 4 2 2 2 ⎠ ν k 1 + νk γˆ (νk) ⎝ + ¯β k=1 ν 2 k 2 + λ2j j =0
(A19)
so that the continuum form of Q2 (t) (Eq. (2.32)) is ∞ ∞ 2 2 421 ¯ 4 2 dλI (λ)λ sin(λt) + 2 dλI (λ)λ cos(λt) Q2 (t) = − 2 40 π π 0 0 ⎡ ⎞⎤ ⎛ ∞ ∞ 2 ¯ 2 ∞ I (λ) cos(λt) I (λ)λ2 λ 8 ⎣ − 2 ⎠⎦ dλ dλI (λ)λ cos(λt) ⎝ dλ 4 π 0 π 0 0 tanh ¯βλ tanh ¯βλ −
8 π 2 ¯β
2
2
∞
21
dλI (λ) sin(λt) 0
−
2
∞
2
dλI (λ)λ cos(λt) 0
2 ∞ ∞ I (λ)λ2 sin(λt) 16 2 2 2 + ν k + νk γˆ (νk) dλ 2 2 − 2 π ¯β k=1 1 (ν k + λ2 ) 0 ∞ ∞ 16 2 2 2 I (λ)λ cos(λt) 2 + 2 ν k 1 + νk γˆ (νk) dλ 2 2 π ¯β k=1 (ν k + λ2 ) 0 2 ∞ ∞ 16 I (λ)λ3 cos(λt) + 2 dλ π ¯β k=1 0 (ν 2 k 2 + λ2 )
(A20)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-10
P. Kumar and E. Pollak
J. Chem. Phys. 141, 234509 (2014)
2. The quantum spectrum
To obtain the necessary expression for the quantum spectrum (Eqs. (2.36) and (2.37)), we use the Wigner formulation. The Wigner representation for the thermal distribution of the jth bath oscillator is tanh(¯βωj /2) 2 1 1 exp − pj + ωj2 qj2 (A21) exp(−βHj )W = 2π ¯ cosh(¯βωj /2) ¯ωj and T r[exp(−βHj )] =
1 . 2 sinh(¯βωj /2)
(A22)
The Wigner representation for the initial wavefunction is ||W
1 0 (q − q0 )2 (p − p0 )2 = . exp − − π¯ ¯ ¯0
(A23)
The Wigner representation for the system propagator is exp(−iHS t/¯)W The initial spectrum is then S(ω, t = 0) = 2 1 = 2π
1 i tan(1 t/2) 2 1 2 2 exp − = p + 1 q . 2π ¯ cos(1 t/2) ¯1
Hˆ dpdqδ ω − S ¯
W
∞
(p − p0 )2 0 (q − q0 )2 − exp − ¯ ¯0
1 dt1 exp(iωt1 ) cos(1 t1 /2) −∞
i tan(1 t1 /2) × exp − ¯
(A24)
0 1 1 + i0 tan(1 t1 /2)
1 0 + i1 tan(1 t1 /2)
0 p02 + 21 q02 + i1 tan(1 t/2) p02 + 20 q02 . [1 + i0 tan(1 t1 /2)][(0 + i1 tan(1 t1 /2))]
(A25)
The equilibrium spectrum is
Hˆ S ˆ exp(−β H ) T r[exp(−β Hˆ )] Sβ (ω) = T r δ ω − ¯ W ∞ exp(iωt1 ) 1 dt = 1/2 , 1 t1 1 t1 t t 2π −∞ 1 cos 2 + i sin 2 Gc cos 21 1 + i sin 21 1 Gs
(A26)
where Gc and Gs are defined in Eqs. (2.28) and (2.29). For the harmonic oscillator system the quantum time evolution of the Wigner distribution is the same as the classical so that ∞ i Hˆ i Hˆ 1 1 1 HS exp − = dt exp(iωt1 ) exp δ ω− ¯ ¯ ¯ 2π −∞ 1 2π ¯ cos(1 t1 /2) W i tan(1 t1 /2) 2 2 2 · exp − pt + 1 qt (A27) ¯1 and the expression for the spectrum is ∞ ∞ ∞ N dpj dxj 1 tanh(¯βωj /2) S (ω, t) = 2 dt1 dpdq 2π ¯ −∞ π¯ −∞ −∞ j =1 i tan(1 t1 /2) 2 exp(iωt1 ) 2 2 exp − (pt + 1 qt ) · cos(1 t1 /2) ¯1 ⎞ ⎛ N 2 2 tanh(¯βω /2) (p − p (q − q ) ) j 0 0 − pj2 + ωj2 qj2 ⎠ , − · exp ⎝− 0 ¯ ¯0 ¯ω j j =1
(A28)
where the time dependence is given explicitly in Eqs. (2.17) and (2.18). This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
234509-11
P. Kumar and E. Pollak
To carry out the integrations in Eq. (A28) we introduce the Fourier representation of the two delta functions as in Eqs. (2.24) and (2.25). One then remains with Gaussian integrals which are readily carried out and the result is given in Eq. (2.37). 1 K.
A. Briney, L. Herman, D. S. Boucher, A. D. Dunkelberger, and F. F. Crim, J. Phys. Chem. A 114, 9788 (2010). 2 J. D. Leonard, Jr., and T. L. Gustafson, J. Phys. Chem. A 105, 1724 (2001). 3 S. Kovalenko, A. Dobryakov, E. Pollak, and N. Ernsting, J. Chem. Phys. 139, 011101 (2013). 4 K. Iwata and H. O. Hamaguchi, Chem. Phys. Lett. 196, 462 (1992). 5 K. Iwata, R. Ozawa, and H. O. Hamaguchi, J. Phys. Chem. A 106, 3614 (2002). 6 K. A. Iwata, K. B. Yoshida, Y. B. Takada, and H. O. Hamaguchi, Chem. Lett. 36, 504 (2007). 7 K. Iwata and H. O. Hamaguchi, J. Mol. Liq. 65–66, 417 (1995). 8 V. Deckert, K. Iwata, and H. O. Hamaguchi, J. Photochem. Photobiol. A 102, 35 (1996). 9 K. Iwata and H. O. Hamaguchi, J. Phys. Chem. A 101, 632 (1997). 10 C. Warmuth, F. Milota, H. F. Kauffmann, H. Wadi, and E. Pollak, J. Chem. Phys. 112, 3938 (2000). 11 D. Gegiou, K. A. Muszkat, and E. Fischer, J. Am. Chem. Soc. 90, 12 (1968). 12 R. M. Hochstrasser, Pure Appl. Chem. 52, 2683 (1980). 13 S. H. Courtney and G. R. Fleming, J. Chem. Phys. 83, 215 (1985). 14 M. W. Balk and G. R. Fleming, J. Phys. Chem. 90, 3975 (1986). 15 S. H. Courtney, M. W. Balk, L. A. Philips, S. P. Webb, D. Yang, D. H. Levy, and G. R. Fleming, J. Chem. Phys. 89, 6697 (1988). 16 S. A. Kovalenko, R. Schanz, H. Hennig, and N. P. Ernsting, J. Chem. Phys. 115, 3256 (2001).
J. Chem. Phys. 141, 234509 (2014) 17 S. A. Kovalenko, A. L. Dobryakov, I. Ioffe, and N. P. Ernsting, Chem. Phys.
Lett. 493, 255 (2010). O. Hamaguchi, Mol. Phys. 89, 463 (1996). 19 H. Wadi and E. Pollak, J. Chem. Phys. 110, 11890 (1999). 20 G. Gershinsky and E. Pollak, J. Chem. Phys. 107, 812 (1997). 21 J. Tatchen and E. Pollak, J. Chem. Phys. 128, 164303 (2008). 22 H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 168, 115 (1988). 23 F. Haake and R. Reibold, Phys. Rev. A 32, 2462 (1985). 24 B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D 45, 2843 (1992). 25 J. J. Halliwell and T. Yu, Phys. Rev. D 53, 2012 (1996). 26 R. Karrlein and H. Grabert, Phys. Rev. E 55, 153 (1997). 27 G. W. Ford and R. F. O’Connell, Phys. Rev. D 64, 105020 (2001). 28 W. T. Pollard and R. A. Friesner, J. Chem. Phys. 100, 5054 (1994). 29 A. Pereverzev, Phys. Rev. E 68, 026111 (2003). 30 E. Pollak, J. S. Shao, and D. H. Zhang, Phys. Rev. E 77, 021107 (2008). 31 A. M. Levine, M. Shapiro, and E. Pollak, J. Chem. Phys. 88, 1959 (1988). 32 I. Rips and E. Pollak, Phys. Rev. A 41, 5366 (1990). 33 E. Pollak, in Dynamics of Molecules and Chemical Reactions, edited by R. E. Wyatt, and J. Z. H. Zhang (Dekker, New York, 1996), p. 617. 34 E. Pollak, J. Chem. Phys. 85, 865 (1986). 35 E. Pollak, Phys. Rev. A 33, 4244 (1986). 36 E. Wigner, Phys. Rev. 40, 749 (1932). 37 U. Weiss, Quantum Dissipative Systems, 3rd ed. (World Scientific, 2008). 38 A. L. Dobryakov, I. Ioffe, A. A. Granovsky, N. P. Ernsting, and S. A. Kovalenko, J. Chem. Phys. 137, 244505 (2012). 39 H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003). 40 H. Wang and M. Thoss, in Multidimensional Quantum Dynamics: MCTDH Theory and Applications, edited by H.-D. Meyer, F. Gatti, and G. A. Worth (Wiley-VCH, 2009), p. 131. 41 I. S. Gradshteyn, and I. M. Ryzhik, in Table of Integrals, Series and Products, edited by A. Jeffrey (Academic Press, New York, 1965). 18 H.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.188.128.74 On: Wed, 03 Jun 2015 12:24:27
