9278
J. Phys. Chem. B 2009, 113, 9278–9288
Multidimensional Spectra via the Mixed Quantum-Classical Liouville Method: Signatures of Nonequilibrium Dynamics Gabriel Hanna* and Eitan Geva* Department of Chemistry, UniVersity of Michigan, Ann Arbor, Michigan 48109-1055 ReceiVed: March 27, 2009; ReVised Manuscript ReceiVed: May 15, 2009
Multidimensional optical spectra are often expressed in terms of optical response functions. These optical response functions consist of contributions from a number of Liouville pathways that differ with respect to the chromophore’s quantum state during the time intervals between light-matter interactions. The dynamics of the photoinactive degrees of freedom during those time intervals are dictated by potential energy surfaces that are explicitly dependent on the chromophore’s quantum state. One therefore expects the system to hop between potential surfaces in a manner dictated by the Liouville pathways and the spectra to reflect the dynamics during the resulting nonequilibrium process. However, the approach commonly used to model spectra of complex condensed-phase systems is based on the ad hoc assumption that the photoinactive degrees of freedom undergo equilibrium dynamics on the potential surface that corresponds to the chromophore’s ground state. In this paper, we formulate optical response in terms of mixed quantum-classical Liouville dynamics, which inherently accounts for the underlying nonequilibrium dynamics. It is shown that, when nonadiabatic transitions are neglected, the resulting formulation is equivalent to that obtained via the linearized semiclassical approximation. We demonstrate the feasibility and utility of the approach by using it to calculate the oneand two-dimensional infrared spectra of the hydrogen stretch of a moderately strong hydrogen-bonded complex dissolved in a dipolar liquid. The results are compared with previously reported spectra that were calculated within the framework of the standard equilibrium ground-state dynamical approach [Hanna and Geva J. Phys. Chem. B 2008, 112, 12991.], thereby shedding light on the spectral signatures of nonequilibrium dynamics in systems of this type. I. Introduction
ˆ jt/p] ) exp[-iH ˆ reft/p] exp+[-i exp[-iH
Linear and nonlinear optical spectroscopy have long been recognized as providing extremely powerful and sensitive probes of the structure and dynamics of molecular systems.1–6 It is also widely accepted that molecularly detailed and dynamically accurate models are necessary in order to take full advantage of these capabilities. However, accomplishing this objective requires overcoming a number of nontrivial theoretical and computational challenges, including the design of spectroscopically accurate force fields, the choice of a proper dynamical treatment, and the development of schemes for calculating spectral signals which are both feasible and self-consistent. Modeling of optical spectra often starts by expressing them in terms of optical response functions (ORFs).1 However, an exact calculation of the fully quantum-mechanical ORFs is not feasible in most cases of practical interest. Mixed quantumclassical methods, which are based on treating a small subset of degrees of freedom (DOF) quantum-mechanically while the rest are treated in a classical-like manner, represent attractive alternatives.1,7–9 However, as is well-known, taking the classical limit of the ORFs with respect to a subset of DOF in a direct manner can lead to expressions which are not unique.1,9–19 In this respect, the choice of a working expression has often been obtained in an ad hoc manner, motivated by computational convenience or speculation rather than by rigorous reasoning. More specifically, one can consider the following identity for the quantum mechanical propagator:1 * To whom correspondence should be addressed.
∫0t dt' Uˆj(t')/p] (1)
where exp+ corresponds to a positively time-ordered exponenˆ ref is an arbitrarily chosen reference Hamiltonian and tial, H
ˆ j(t) ) exp[iH ˆ reft/p](H ˆj - H ˆ ref) exp[-iH ˆ reft/p] U
(2)
ˆ ref ) Applying the identity in eq 1, with the specific choice of H ˆ g, to an ORF will yield an exact expression involving timeH ordered exponentials and time-evolved operators of the form ˆ gt/p]. A “classical limit” of this ˆ gt/p]Aˆ exp[-iH Aˆ(t) ) exp[iH expression is then taken by replacing Aˆ(t) by its classical analogue A[Qt], where Qt is a classical trajectory on the groundstate potential energy surface (PES), and the time-ordered exponentials by regular exponentials. The resulting working expression for the ORF reflects only ground-state equilibrium dynamics, which can be traced back to the above ad hoc choice of reference Hamiltonian. However, other equally justifiable choices will yield different working expressions of the “classical limit” that involve nonequilibrium dynamics on the reference potential of choice. In order to avoid making such ad hoc approximations and to arrive at more rigorous expressions, several approaches for calculating linear and nonlinear ORFs have been proposed based on semiclassical treatments, including the cellular dynamics method of Mukamel and co-workers,1,109 the mixed-state propagation method of Loring and co-workers,110,111 linearization
10.1021/jp902797z CCC: $40.75 2009 American Chemical Society Published on Web 06/16/2009
Mixed Quantum-Classical Liouville Method schemes,10,16,17,19 and methods based on the forward-backward initial value representation (FB-IVR).17,19 However, these methods were mostly applied to relatively simple model systems. At the same time, applications to complex condensed-phase systems containing many DOF have often been based on equilibrium classical molecular dynamics simulations of the photoinactive DOF on the PES that corresponds to the chromophore’s ground state. In this paper, we will refer to this practice as the standard approach. The main advantage of the standard approach is that ORFs can be efficiently calculated from equilibrium molecular dynamics (MD) simulations on the ground-state PES, as opposed to nonequilibrium MD simulations involving multiple PESs which are computationally more expensive. More specifically, the fact that the time origin can be set anywhere along the equilibrium trajectory enhances the computational efficiency by making it possible to average over a large number of, possibly overlapping, segments of this trajectory. The standard approach can be justified in the case of weak coupling between the chromophore and its environment since in this case the PES that dictates the dynamics of the photoinactive DOF is expected to be insensitive to the quantum state of the chromophore. However, this is not the case when the chromophore is coupled strongly to its environment, so that the PESs differ significantly from one quantum state of the chromophore to another. As a result, the photoinactive DOF tend to undergo significant rearrangements upon changing the chromophore’s quantum state, which is expected to give rise to measurable spectral signatures. In this paper, we formulate an alternative approach toward calculating ORFs, which is based on mixed quantum-classical Liouville (MQCL) dynamics.20–35 Importantly, MQCL dynamics provides a rigorous route for deriving unambiguous mixed quantum-classical expressions for the ORFs that reflect nonequilibrium dynamics on multiple PESs, rather than equilibrium ground-state dynamics. We show that the MQCL-based algorithm for calculating ORFs is equivalent to that obtained via the linearized semiclassical (LSC) approximation17,19 in the absence of nonadiabatic transitions. The applicability and usefulness of this approach is demonstrated in the case of the Azzouz-Borgis (AB) model36 for a moderately strong hydrogenbonded phenol-amine complex dissolved in liquid methyl chloride, previously studied by our group within the framework of the standard equilibrium ground-state dynamics methodology.37 In this context, it should be noted that two-dimensional infrared (2DIR) spectroscopy has been recently used to probe various types of hydrogen-bonded systems.6,38–51 These experiments were accompanied by substantial theoretical and computational studies aimed at understanding the signature of the underlying molecular structure and dynamics on the spectra.1,43,52–102 However, relatively little attention has been given to the signature of the above-mentioned nonequilibrium processes on these spectra. Thus, although the present study does not aim at reproducing specific experimental spectra, we hope that the general approach as well as its application to the AB model will shed light on previously unexplored aspects of IR spectroscopy in systems of this type. The goal of this paper is two-fold: (1) formulate the optical response of a driven system within the framework of MQCL dynamics and show its equivalence to the LSC approximation in the absence of nonadiabatic transitions; (2) use the MQCL scheme in order to gain insight into the spectral signature of nonequilibrium solvation processes on the multidimensional spectra of hydrogen-bonded systems. The remainder of this paper is arranged as follows. In section II, we consider the
J. Phys. Chem. B, Vol. 113, No. 27, 2009 9279 MQCL dynamics of a system driven by an explicitly timedependent classical radiation field. In section III, we present a formulation of optical response in terms of MQCL dynamics, with emphasis on one-dimensional (1D) and two-dimensional (2D) spectroscopy of a three-state system. In section IV, we demonstrate the applicability and usefulness of the MQCL formalism by using it to compute the 1D and 2D IR spectra of the hydrogen stretch of a hydrogen-bonded complex dissolved in a dipolar liquid. We conclude with a summary of the main results in section V. II. Mixed Quantum-Classical Liouville Dynamics of a System Driven by a Classical Radiation Field Consider a quantum system driven by a classical radiation field, whose Hamiltonian is explicitly time dependent and given by
pˆj2 + 2mj j)1 n
ˆ (t) ) H
∑
Pˆj2 ˆ ) + W(rˆ, t) + V(rˆ, R 2Mj J)1 N
∑
(3)
Here, {m1, ...,mn}, rˆ ) (rˆ1, ...,rˆn), and pˆ ) (pˆ1, ...,pˆn) are the masses, positions, and momenta of the photoactiVe DOF (i.e., ˆ ) (R ˆ 1, ...,R ˆ n), the chromophore), respectively; {M1, ...,MN}, R ˆ ˆ ˆ and P ) (P1, ...,PN) are the masses, positions, and momenta of ˆ ) is the photoinactiVe DOF (i.e., the bath), respectively; V(rˆ,R the field-free potential energy, and W(rˆ,t) is the field-matter interaction term. The latter is assumed to be explicitly given by
W(rˆ, t) ) -µ(rˆ) · E(t) cos(ωt - φ)
(4)
where E(t), ω, and φ are the envelope, leading frequency, and phase of the driving field, respectively, and µ(rˆ) is the chromophore’s dipole moment operator [throughout this paper, vectors are represented by boldface symbols (e.g., A) while operators are capped (e.g., Aˆ)]. The state of the overall system (i.e., photoactive and photoinactive DOF) can be described by a density operator, Fˆ (t), whose time evolution is dictated by the quantum Liouville equation:
d i ˆ Fˆ (t) ) - [H (t), Fˆ (t)] dt p
(5)
The state of the system can also be represented, without loss of generality, by the corresponding partial Wigner transform103 over the photoinactive DOF of Fˆ (t):
Fˆ W(R, P;t) )
( 2πp1 ) ∫ dZ e N
iP · Z/p
〈R - Z/2|Fˆ (t)|R + Z/2〉
(6) where Fˆ W(R,P;t) is an operator in the chromophore’s Hilbert space. The partial Wigner representation of an observable represented by the operator Aˆ is similarly defined by
AˆW(R, P) )
∫ dZ eiP · Z/p〈R - Z/2|Aˆ|R + Z/2〉
so that expectation values can be obtained as follows:
(7)
9280
J. Phys. Chem. B, Vol. 113, No. 27, 2009
〈Aˆ〉(t) )
∫ dR ∫ dP Tr[Fˆ W(R, P;t)AˆW(R, P)]
Hanna and Geva
(8)
where Tr denotes a trace over the photoactive DOF. The MQCL equation describes the dynamics of Fˆ W(R,P;t) to first order in (mj/MJ)1/2, and is therefore expected to be accurate when mj/MJ , 1:22,25–30
∂ i ˆ (R, P;t), Fˆ W(R, P;t)] + Fˆ (R, P;t) ) - [H ∂t W p W 1 ˆ ˆ W(R, P;t)}) ({HW(R, P;t), Fˆ W(R, P;t)} - {Fˆ W(R, P;t), H 2 (9)
∂ RR′ F (R, P;t) ) -iωRR′(R)FRR′ W (R, P;t) ∂t W f FRR(R) + FR′R′(R) f RR′ - U · ∇RFRR′ · ∇PFW (R, P;t) W (R, P;t) 2 i Rβ βR′ [WRβ(R;t)FβR′ W (R, P;t) - W (R;t)FW (R, P;t)] p β f 1 * * dR′β · U 1 + SR′β (R) · ∇P FRβ W (R, P;t) 2 β*R′ f 1 dRβ · U 1 + SRβ(R) · ∇P FβR′ W (R, P;t) 2 β*R (14)
∑
∑
∑
r where ∇R/P and ∇R/P correspond to taking the gradient with respect to R/P of the term to the right and to the left, respectively. In keeping with the mass separation between the photoactive and photoinactive DOF, we choose to cast the operator MQCL equation, eq 9, in terms of a representation that consists of chromophore states that follow the photoinactive DOF adiabatically. To this end, we first define the field-free adiabatic Hamiltonian: f
pˆj2 + V(rˆ, R) 2mj j)1
[
]
]
f f ˙ , FRR(R) ) -∇ Here, U ) R RER(R), dRβ(R) ≡ 〈R(R)|∇R|β(R)〉 ) -d*βR(R), ωRR′(R) ) [ER(R) - ER′(R)]/p, SRβ ) pωRβ(R)dRβ(R)/ dRβ(R) · U and
Here, { · · · } corresponds to the Poisson bracket which is explicitly given by
f r ∇ {AˆW(R, P;t), BˆW(R, P;t)} ) AˆW(R, P;t)(∇ P R f r ∇ ˆ ∇ R P)BW(R, P;t) (10)
[
WRβ(R;t) ) -µRβ(R) · E(t) cos(ωt - φ)
(15)
It should be noted that the dependence of µRβ(R) ≡ 〈R(R)|µ(rˆ)|β(R)〉 on R reflects the parametric dependence of the adiabatic functions on R, rather than any explicit dependence of the dipole moment operator, µ(rˆ), on R. The RHS of eq 14 contains adiabatic terms, which couple RR′ (R,P;t) to itself, as well as the phase space density FW nonadiabatic coupling terms, which couple FRR′ W (R,P;t) to either Rβ βR′ (R,P;t) (where β * R′) or FW (R,P;t) (where β * R). The FW nonadiabatic coupling terms are responsible for population relaxation processes which typically occur on time scales longer than that of the dynamics described by the adiabatic terms. Also, these nonadiabatic transitions can often be described in terms of rate constants (within Fermi’s golden rule) and are known to diminish the spectral signal. In the present study, we will neglect the nonadiabatic terms, thereby restricting ourselves to time scales that are shorter than the time scale required for them to affect the spectra. Under those conditions, eq 14 simplifies considerably, to give
n
ˆ ad ) H
∑
(11)
We denote the eigenstates and eigenvalues of the adiabatic Hamiltonian, which are parametrically dependent on R, by {|R(R)〉} and {ER(R)}, respectively, so that
ˆ ad |R(R)〉 ) ER(R)|R(R)〉 H
(12)
The partially Wigner transformed density matrix elements in the adiabatic representation are then given by
FRR′ ˆ W(R, P;t)|R′(R)〉 W (R, P;t) ) 〈R(R)|F
(13)
It should be noted that the dependence of FRR′ W (R,P;t) on R arises from the R dependence of the adiabatic states as well as the fact that Fˆ W(R,P;t) is explicitly R-dependent. The equation of motion for FRR′ W (R,P;t), which can be obtained from eq 9, is given by22,25–30
∂ RR′ F (R, P;t) ) -iωRR′(R)FRR′ W (R, P;t) ∂t W RR R′R′ b FRR′(R, P;t) - F (R) + F (R) · ∇ b FRR′(R, P;t) - U·∇ R W P W 2 i βR′ Rβ Rβ βR′ [W (R;t)FW (R, P;t) - W (R;t)FW (R, P;t)] p β (16)
∑
It should be noted that, unlike the field-free MQCL equation, eq 16 still allows for coupling between adiabatic states via the matrix elements of the field-matter interaction, {WRβ(R;t)}. However, unlike the above-mentioned irradiative nonadiabatic transitions, the timing of the field-induced transitions can be controlled by employing a carefully selected sequence of impulsive light pulses (see section III). III. Optical Response Functions from MQCL Dynamics A. Preliminary Considerations. In this section, we will introduce a general theoretical framework for calculating ORFs in the case of a system whose dynamics is described by eq 16. For the sake of concreteness, we will do so in the context of a three-state system. The three states will be denoted by |0(R)〉, |1(R)〉, and |2(R)〉, corresponding to the ground, first excited, and second excited adiabatic states, respectively. We will also
Mixed Quantum-Classical Liouville Method
J. Phys. Chem. B, Vol. 113, No. 27, 2009 9281
assume that the only allowed transitions are between |0(R)〉 and |1(R)〉 and between |1(R)〉 and |2(R)〉. Finally, we will make the rotating wave approximation (RWA), so that the field-matter interaction term can be put in the following form:
ˆ (t) ) - pχ01(R, t)e-iφeiωt |0(R)〉〈1(R)| W - pχ10(R, t)eiφe-iωt |1(R)〉〈0(R)| - pχ12(R, t)e-iφeiωt |1(R)〉〈2(R)| - pχ21(R, t)eiφe-iωt |2(R)〉〈1(R)|
(17)
where the Rabi frequencies are given by * χRβ(R, t) ) χβR (R, t) )
1 µRβ(R) · E(t) 2 p
(18)
We assume that the phase of the field is given by φ ) k · Rc, where k is the wave vector that dictates the directionality of the light pulse and Rc represents the location of the chromophore in the lab frame. Solving eq 16 with the field-matter interaction Hamiltonian described in eq 17 is facilitated by making the transformation to the rotating frame:
ˆ rott/p]Fˆ W(R, P;t) exp[-iH ˆ rott/p] F ˜W(R, P;t) ) exp[iH (19) ˆ rot ) pω|1(R)〉〈1(R)| + 2pω|2(R)〉〈2(R)|. This implies where H that
F ˜RR′ ˆ RR′ W (R, P;t) ) exp(i(R - R′)ωt)F W (R, P;t)
{
E0 t0 - τ0 /2 e t e t0 + τ0 /2 0 otherwise
∂ RR′ F˜ (R, P;t) ) -i∆RR′(R)F˜ RR′ W (R, P;t) ∂t W RR R′R′ b F˜ RR′(R, P;t) - F (R) + F (R) · ∇ b F˜ RR′(R, P;t) - U·∇ R W P W 2 i ˜ Rβ(R;t)F˜ βR′ ˜ βR′ [W ˜ Rβ W (R, P;t) - W (R;t)F W (R, P;t)] p β (21)
∑
where the detuning is given by
(24)
Since we have applied the RWA (i.e., assumed that the leading frequency of the pulse is in resonance with the 0 T 1 and 1 T 2 transitions), we can assume that ∆RR′(R) ≈ 0 when the pulse actually impacts the system. It should be noted that, in cases where these transitions give rise to multiple bands, several experiments, in which the leading frequency of each pulse is in resonance with one band at a time, may be necessary in order to obtain the spectrum over the entire spectral range. Setting ∆RR′(R) ) 0 and assuming that the pulse is impulsive (i.e., that no bath dynamics takes place during τ0) leads to the following equation that prescribes the evolution during the pulse:
∂ RR′ i ˜ (R, P;t) ) F ∂t W p
∑ [W˜Rβ(R;t)FβR′ W (R, P;t) β
˜ βR′(R;t)FRβ W W (R, P;t)] (25) which can be equivalently written in the following operator form:
d i ˜ ˜W(R, P;t) ) - [W ,F ˜W(R, P;t)] F dt p
(26)
˜ is time-independent and assuming that the light Noting that W pulses are weak, the impact of the pulse on the state of the system can be described by ˜
(20)
which upon substitution into eq 16 yields the following MQCL equation of motion in the rotating frame:
∆RR′ ) ωRR′ - (R - R′)ω
E(t) )
˜
F ˜W(R, P;t0+) ) e-iWτ0/pF ˜W(R, P;t0-)eiWτ0/p ˜τ0 /p]F ˜τ0 /p] (27) ≈ [1ˆ - iW ˜W(R, P;t0-)[1ˆ + iW where t0( ) t0 ( τ0/2. Finally, we note that the field-free dynamics in between impulsive pulses can be described by eq ˜ Rβ ) 0. 21 with W B. Linear Optical Response. The first scenario that we will consider corresponds to the measurement of the system’s linear response and involves an interaction with a single impulsive pulse at time t ) 0. We assume that prior to the interaction with the pulse, the chromophore is in its ground state and the photoinactive DOF are in the corresponding thermal equilibrium state: 00 F ˜W(R, P;0-) ) |0(R)〉F ˜eq,W (R, P)〈0(R)|
(28)
(22)
and the nonvanishing matrix elements of the field-matter ˜ (t), are given by interaction term in the rotating frame, W
00 Here, F˜ eq,W (R,P) is the Wigner transform of the ground-state equilibrium density operator:
ˆ
˜ 01(R;t) ) [W ˜ 10(R;t)]* ) -pχ01(R, t)e-iφ, W ˜ 21(R;t)]* ) -pχ12(R, t)e-iφ ˜ 12(R;t) ) [W W
Fˆ 00 eq
)
F ˜00 eq
)
e-H0/kBT ˆ
Tr[e-H0/kBT]
(29)
(23)
We next consider the interaction of the system with an impulsiVe light pulse with a square envelope of width which is centered at some time t0:
where kB is the Boltzmann constant, T is the absolute temperN ˆ ). In what follows, we ˆ 0 ) ∑J)1 (Pˆ2J )/(2MJ) + E0(R ature, and H 00 will also assume that F˜ eq,W(R,P) can be well approximated by its classical limit:
9282
J. Phys. Chem. B, Vol. 113, No. 27, 2009 00 F ˜eq,W (R, P) ≈
e-H0(R,P)/kBT
∫ dR ∫ dP e-H (R,P)/k T 0
Hanna and Geva
(30)
B
Using eq 27, the impact of the pulse is described by 00 F ˜W(R, P;0+) ) |0(R)〉F ˜eq,W (R, P)〈0(R)| 00 ˜eq,W (R, P)〈1(R)| - iχ01(R)τ0e-iφ |0(R)〉F 00 ˜eq,W (R, P)〈0(R)| + iχ10(R)τ0eiφ |1(R)〉F 00 ˜eq,W (R, P)〈1(R)| + χ10(R)χ01(R)τ20 |1(R)〉F
(31)
Thus, the pulse is seen to give rise to a coherence between the ground and first excited states which is manifested by a reshuffling of the phase space densities, so that the phase space density that corresponded initially to the 00 density matrix element is split between the 00, 01, 10, and 11 density matrix elements. It should be noted that each phase space point associated with the 10 or the 01 density matrix elements also acquires the weighting factors iχ10(R)τ0eiφ or -iχ01(R)τ0e-iφ, respectively, as a result of the interaction with the light pulse. The field-free dynamics following the pulse is described by ˜ Rβ ) 0. This equation dictates that the phase space eq 21 with W densities associated with the 00, 01, 10, and 11 density matrix elements are to be propagated classically and independently from one another on the corresponding PESs. In particular, according to eq 21, the phase space densities associated with the 10 and 01 density matrix elements are to be propagated on the aVerage 10 01 (Rt,Pt;t) and F˜ W (Rt,Pt;t) can PES, [E0(R) + E1(R)]/2. Thus, F˜ W 01 (R P ;0 ) and F ˜ (R P ;0 ), by propagating be obtained from F˜ 10 W 0, 0 + W 0, 0 + each phase space point on the average PES for a period of time t and associating with it a phase factor of the form exp[-i∫0t dτ ∆10(Rτ)] or exp[i∫0t dτ ∆10(Rτ)], respectively. The signal field is assumed to be proportional to the expectation value of the dipole moment operator at time t in the direction of the incoming field that is associated with the phase factor e-iφ ) e-ik · Rc. After a transformation back to the nonrotating frame, this yields the following MQCL expression for the linear ORF:
JMQCL(t) )
|
∫0t dτ ω10(Rτ)] 00-01
≡ 〈χ10(R0)χ01(Rt) exp[i
∫0t dτ ω10(Rτ)]〉00-01
∫0∞ dt e-iωtJMQCL(t)
t1
t1+t2+t3
t1
t1+t2+t3
× e(i ∫ 0 dτω10(Rτ)-i ∫ t1+t2 dτω10(Rτ)〉00-01-11-01 - 〈χ10(R0)χ01(Rt1)χ12(Rt1+t2)χ21(Rt1+t2+t3)
(32)
Here, the subscript 00-10 indicates that Rτ is obtained by classical propagation on the aVerage PES, [E0(R) + E1(R)]/2, with initial conditions {R0,P0}, which are sampled based on the classical equilibrium ground-state phase space density, 00 (R0,P0). Finally, the linear absorption spectrum is defined FW,eq by1
I(ω) ) [1 - e-βω]Re
Rr/nr(t3, t2, t1) ) 〈χ01(R0)χ10(Rt1)χ10(Rt1+t2)χ01(Rt1+t2+t3) × e(i ∫ 0 dτω10(Rτ)-i ∫ t1+t2 dτω10(Rτ)〉00-01-00-01 + 〈χ01(R0)χ10(Rt1)χ10(Rt1+t2)χ01(Rt1+t2+t3)
∫ dR0 ∫ dP0 F00eq(R0, P0)χ10(R0)χ01(Rt) × exp[i
the coherences is actually dictated by the MQCL equation, which should be contrasted with the ad hoc choice of the groundstate PES for performing the dynamics within the standard approach. The difference between the two approaches will obviously be negligible in the limit of inhomogeneous broadening, where the linear ORF decays to zero on time scales faster than those associated with motions of the photoinactive DOF, or when the ground and excited PESs are similar (e.g., in the case of weak coupling between the photoactive and photoinactive DOF). However, in systems where this is not the case, eq 32 may give rise to a pronounced spectral signature of the nonequilibrium dynamics on the average PES on the absorption spectrum. C. Third-Order Optical Response. The second scenario that we will consider corresponds to the measurement of 2D spectra, which involves using three subsequent laser pulses with wave vectors ka, kb, and kc, and time delays t1 (between pulses a and b) and t2 (between pulses b and c). These pulses create a thirdorder polarization in the sample, which gives rise to a signal field that is heterodyne detected at a time interval t3 after pulse c, in the background-free directions kr ) -ka + kb + kc and knr ) ka - kb + kc, corresponding to the rephasing and nonrephasing signals, respectively. We therefore consider a scenario where a system, whose initial state is given by eq 28, interacts with a sequence of three impulsive square pulses at times 0, t1 and t1 + t2, followed by the detection of the signal field at time t1 + t2 + t3. As for linear response, the field-free dynamics between pulses is assumed to be governed by eq 21 ˜ Rβ ) 0. The three pulses are also assumed to differ with with W respect to their phase, so that φa, φb, and φc correspond to the phases of the first, second, and third pulses, respectively. We will also focus our attention on the rephasing and nonrephasing signals which can be distinguished by the fact that they correspond to the different overall phases φr ) -φa + φb + φc and φnr ) φa - φb + φc, respectively. Following the same procedure as in the case of linear response leads to the following MQCL expressions for the rephasing and nonrephasing signals:
(33)
Although eq 32 resembles the standard mixed quantumclassical expression for the linear ORF, it is distinctly different from it by the fact that the dynamics takes place on the aVerage PES, rather than on the ground-state PES. Furthermore, the use of the average PES for the phase space densities associated with
t1
× e(i ∫ 0
+t2+t3 dτω10(Rτ)-i ∫ tt1+t dτω21(Rτ) 1
2
〉00-01-11-12
(34)
Here, ( corresponds to the rephasing (r f +) and nonrephasing (nr f -) signals, and the subscripts 00-01-00-01, 00-0111-01, and 00-01-11-12 correspond to the sequence of PESs used for propagating the phase space point during the various time intervals t1, t2, and t3. For example, 00-01-11-12 implies sampling from an equilibrium distribution on the ground state PES (00), followed by dynamics on the average PES [E0(R) + E1(R)]/2 during t1 (01), the first excited state PES E1(R) during t2 (11), and on the average PES [E1(R) + E2(R)]/2 during t3 (12). Finally, the two-dimensional spectrum is defined by:5,55 I(ω3, t2, ω1) ≡ Re
∫
∞
0
dt1
∫
∞
0
dt3{ei(ω1t1+ω3t3)Rnr(t3, t2, t1) + ei(-ω1t1+ω3t3)Rr(t3, t2, t1)}
(35)
Mixed Quantum-Classical Liouville Method
J. Phys. Chem. B, Vol. 113, No. 27, 2009 9283
Although eqs 34 resemble the standard mixed quantumclassical expressions for the rephasing and nonrephasing signals,37,55 they are distinctly different from them by the fact that the latter are calculated based on equilibrium classical trajectories that take place solely on the ground-state PES. It should also be noted that, as in the case of linear response, the sequence of PESs indicated for the different terms in eq 34 is actually dictated by the MQCL equation, rather than chosen in an ad hoc manner. Finally, as t2 f τeq (where τeq collectively refers to the time scales on which equilibrium is reached on the ground and first excited state PESs), one expects the dynamics during t1 to become uncorrelated with that during t3, so that
IV. Application to the Azzouz-Borgis Model of a Moderately Strong Hydrogen-Bonded Complex in a Dipolar Liquid
(i ∫ t01 dτω10(Rτ)
Rr/nr(t3, t2, t1) f 〈χ01(R0)χ10(Rt1)e × 〈χ10(R0)χ01(Rt3)e-i
+ 〈χ01(R0)χ10(Rt1)e(i × 〈χ10(R0)χ01(Rt3)e-i
〉00-01
+t2+t3 ∫tt1+t 1 2
∫ t01
- 〈χ10(R0)χ01(Rt1)e(i
∫ t01
× 〈χ12(R0)χ21(Rt3)e-i
〉00-01
dτω10(Rτ)
〉00-01
dτω10(Rτ)
+t2+t3 ∫ tt1+t 1 2
〉11-01
dτω10(Rτ)
〉00-01
dτω10(Rτ)
+t2+t3 ∫ tt1+t 1 2
dτω21(Rτ)
〉11-12
arithmetic mean of the ground- and excited-state surfaces. The latter approximation is based on replacing the product of quantum propagators by the corresponding forward-backward semiclassical propagators. In ref 19, Geva and co-workers extended their analysis in ref 17 by deriving expressions for the LSC and FB-IVR 3PE ORFs, within the Condon approximation and in the context of a two-state system. In the LSC expression, the photoinactive DOF hop between PESs in the same fashion as dictated by the MQCL approach. This implies that the MQCL approach for calculating ORFs is equivalent to the LSC approach in the absence of nonadiabatic transitions.
(36)
Thus, within the MQCL scheme, the asymptotic 2D spectrum A A at t2 ) τeq is proportional to the real part of I01 (ω1) × [I01 (ω3)+ E A A E A (ω3) - I12 (ω3)], where I01 (ω), I01 (ω) and I12 (ω) correspond I01 to the MQCL absorption spectrum for the 0 f 1 transition, the emission spectrum for the 1 f 0 transition and the absorption spectrum for the 1 f 2 transition, respectively (where we start at thermal equilibrium on the first excited state PES in the case of the 1 f 0 and 1 f 2 transitions). In contrast, the asymptotic 2D spectrum predicted by the standard approach is proportional A A A A A (ω1) × [2Ij01 (ω3) - jI12 (ω3)], where jI01 (ω) and jI12 (ω) to jI01 correspond to the standard absorption spectra for the 0 f 1 and the 1 f 2 transitions, respectively, and where we always remain at thermal equilibrium on the ground-state PES. We conclude this section by comparing our theoretical results to those previously obtained via semiclassical and mixed quantum-classical approaches. The first- and third-order ORFs derived above (eqs 32 and 34) apply to absorption and threepulse photon echo (3PE) spectroscopy in three-state systems and incorporate non-Condon effects. It should be noted that similar expressions exist in the literature, but they are not identical to those derived herein. Loring and co-workers10 derived a two-pulse photon echo (2PE) ORF (employing the Condon approximation), which dictates that the dynamics of the photoinactive DOF be propagated on the arithmetic mean of the ground- and excited-state surfaces. Kim and co-workers11 started with a mixed quantum-classical expression for the 2PE ORF involving propagation with a reference Hamiltonian and made two different ad hoc choices for it: the ground-state Hamiltonian and the arithmetic mean of the ground- and excitedstate Hamiltonians. Along the same lines, ref 1 suggests an ad hoc prescription for calculating the 3PE ORFs for a two-state system. In ref 17, Geva and co-workers present results for the LSC and FB-IVR approximations of the linear and 2PE ORFs, employing the Condon approximation. The former approximation is obtained by linearizing the forward-backward action in the path integral expressions for the ORFs, with respect to the difference between the forward and backward paths. Within this approach, the photoinactive DOF are also propagated on the
In this section, we will use the MQCL methodology outlined in the previous section in order to elucidate the signatures of nonequilibrium dynamics on the 1D and 2D IR spectra of the hydrogen (H) stretch of a moderately strong H-bonded complex dissolved in a dipolar liquid. To this end, we will compare the spectra obtained via the MQCL methodology to the recently reported spectra obtained within the framework of the standard equilibrium ground state approach.37 A. Model and Simulation Techniques. In this subsection, we provide a brief outline of the Azzouz-Borgis (AB) model of a moderately strong H-bonded complex, AHB, dissolved in a dipolar liquid (a more detailed description is available in refs 32 and 34–36). Within this model, it is assumed that the proton moves along a 1D axis connecting the donor and acceptor. The donor, A, and acceptor, B, are modeled as single particles and parametrized to represent phenol and trimethylamine, respectively. The intramolecular PES as a function of the A-B and A-H distance is the same as that used in refs 34, 37, and 104 with the equilibrium A-B distance set to 2.7 Å, which corresponds to a moderately strong H-bond. The charges on A and B are assumed to be explicitly dependent on the position of the proton. The intramolecular PES has a double-well profile as a function of the proton displacement, thereby giving rise to tautomeric equilibrium between covalent and ionic forms: AH - B h A- - H+B. The solvent is assumed to consist of a liquid of methylchloride molecules, which are modeled as rigid dipolar diatomic molecules. The complex-solvent and solvent-solvent interactions are modeled in terms of site-site Lennard-Jones (LJ) and Coulomb interactions, and the corresponding force field parameters are as in refs 34, 37, and 104. Although the covalent form, AH-B, is favored in the absence of the solvent, the ionic form, A-H+B, becomes more stable in the dipolar solvent. The photoactive quantum DOF is taken to be the highfrequency H-stretch, and we assume that the only allowed transitions are between the ground to the first excited and between the first excited to the second excited vibrational states of the H-stretch. The associated transition dipole moments are given by µRβ(R) ∝ 〈R(R)|rˆ|β(R)〉, where rˆ ) rˆuAB is the position operator corresponding to the H-stretch and uAB is the unit vector between A and B. It should also be noted that, to a good approximation, uAB remains parallel to E(t) in the body-fixed frame of the complex. The remaining DOF, which constitute the bath, are treated classically and are assumed to be photoinactive and to follow the H-stretch adiabatically. The vibrational energy levels and wave functions are obtained by diagonalizing the adiabatic protonic Hamiltonian on the fly using the LAPACK DSYGV routine.34,37,104
9284
J. Phys. Chem. B, Vol. 113, No. 27, 2009
Figure 1. Ground (00, blue), first excited (11, red), and second excited (22, green) free energy surfaces as a function of the solvent polarization. Also shown are the averaged free energy surfaces of the ground and first excited states (01, magenta) and the first excited and second excited states (12, orange).
Molecular dynamics simulations were performed via the velocity Verlet algorithm with a time step of 1 fs. The simulations were performed within a cubic simulation box with periodic boundary conditions, containing a single AHB complex and 255 methyl chloride molecules at a temperature and density of 250 K and 0.012 Å-3, respectively. The LJ potentials were spherically truncated at Rc ) 13.8 Å and shifted accordingly. The Coulomb potentials were smoothly truncated to zero at Rc ) 13.8 Å. The linear ORF was calculated based on eq 32 by averaging over 16 000 nonequilibrium trajectories of length t ) 250 fs, giving rise to a 250-point time grid with a time interval of 1 fs, which was then padded with zeros to generate a 2048-point time grid. The 1D Fourier transform required for computing the 1D IR spectrum, eq 33, was carried out numerically via the FFT routine. The third-order ORFs in eq 34 were obtained by averaging over 16 000 nonequilibrium trajectories for each point on a (t1,t3) time grid with 0 < t1 < 250 fs (time interval of 5 fs) and 0 < t3 < 250 fs (time interval of 1 fs) for various values of t2. The 50 × 250 (t1,t3) 2D time grid was padded with zeros to generate a 2048 × 2048 time grid. It should be noted that it is not necessary to run a separate nonequilibrium trajectory for each point on the (t1,t3) grid, which would have amounted to an overall number of 50 × 250 × 16 000 ) 2 × 108 trajectories. Instead, one can obtain all the necessary information for all t3 grid points, for a given value of t1, from a single trajectory, which reduces the overall number of trajectories to 50 × 16 000 ) 8 × 105. Finally, the 2D Fourier transform required for computing the 2D IR spectrum was carried out numerically via the FFT routine. B. Free Energy Surfaces and Transition Frequencies. The 1D and 2D IR spectra reported below directly probe the transitions between the ground and first excited states and between the first and second excited states of the H-stretch. The corresponding vibrational energy levels and wave functions follow the bath DOF adiabatically and are therefore explicitly dependent on the bath coordinates. It is convenient to consider this dependence in terms of a collective order parameter, such as the solvent polarization, ∆E, which is defined as the difference between the solvent electrical potentials at the minima of the covalent and ionic wells.32,34,35 Figure 1 shows the ground (00), first excited (11), and second excited (22) free energy surfaces (FESs), as well as the averages of the ground and first excited state FESs (01) and the first and second excited FESs (12), as a function of the solvent polarization. The fact that the FESs are very different suggests that the same is true for the corresponding multidimensional PESs and serves as an indication for the strong coupling between the H-stretch and the bath.
Hanna and Geva
Figure 2. Averaged fundamental transition frequency, ω10, along the Liouville pathways 00 f 10 f 00 f 10 (solid line) and 00 f 10 f 11 f 10 (dashed line). The results were generated for t1 ) 100 fs, t2 ) 300 fs, and t3 ) 100 fs.
The ground-state FES is seen to have a double-well form, which reflects coexistence between the ionic and covalent tautomers (equilibrium composition: 65% ionic and 35% covalent; proton transfer rate constant: 0.16 ps-1).35 In contrast, the first excited FES has the shape of a symmetric single well whose minimum is centered in the vicinity of the ground-state transition state between the ionic and covalent tautomers. It is important to note that the ionic/covalent tautomers and transition state correspond to stable and unstable solvent configurations, respectively, when the chromophore is in its ground state. In contrast, the ionic/covalent tautomers and transition state correspond to unstable and stable solvent configurations, respectively, when the chromophore is in its first excited state. Finally, we note that the second excited FES has the shape of an asymmetric single well, which is broader than that of the first excited FES. The rather pronounced differences between the various surfaces shown in Figure 1 suggest that any transition between them will be accompanied by similarly pronounced solvent reorganization. This, in turn, implies that one can expect significant signatures of nonequilibrium dynamics on the 1D and 2D IR spectra of the H-stretch, thereby making this an appropriate model system for exploring those effects. This is demonstrated in Figure 2, which shows the average fundamental transition frequency, ω10, along two of the three Liouville pathways which contribute to the 2D spectra, that is, 00 f 01 f 00 f 01 (solid line) and 00 f 01 f 11 f 01 (dashed line). The trajectories that correspond to both pathways are seen to be different, drifting upward and downward rather than fluctuating around a fixed equilibrium value. This provides a clear indication of the nonequilibrium nature of the underlying dynamics. The oscillations observed along the nonequilibrium spectral trajectories can be traced back to the coupling with the low-frequency intramolecular A-B stretch, suggesting that the equilibrium position of this vibrational mode also changes significantly with the chromophore’s quantum state. Finally, in the upper panel of Figure 3, we show the distributions of the fundamental transition frequency, ω10, that correspond to equilibrium on the ground-state PES (red) and first excited state PES (green). We also show the distribution of the overtone transition frequency, ω21, that corresponds to equilibrium on the first excited state PES (blue). The transition frequency distributions are seen to be rather broad and pronouncedly different in range and in shape, all of which are manifestations of the strong coupling between the H-stretch and the bath. In the lower panel of Figure 3, we show the dependence of the fundamental transition dipole moment, µ10, averaged over all the configurations that correspond to a given value of the transition frequency, on the fundamental transition frequency,
Mixed Quantum-Classical Liouville Method
Figure 3. Upper panel: The distributions of the fundamental transition frequency, ω10, corresponding to equilibrium on the ground-state PES (red) and first excited state PES (green), and the distribution of the overtone transition frequency, ω21, corresponding to equilibrium on the first excited state PES (blue). Lower panel: The average transition dipole moment, µ10, as a function of the fundamental transition frequency, ω10, at equilibrium on the ground-state PES (red) and first excited state PES (green), and the average overtone transition dipole moment, µ12, as a function of the corresponding transition frequency, ω21, at equilibrium on the first excited state PES (blue).
Figure 4. 1D IR spectra of the H-stretch calculated using the standard ground-state equilibrium dynamical approach (dashed line) and the nonequilibrium MQCL approach (solid line).
ω10, at equilibrium on the ground state PES (red) and first excited state PES (green). We also show the average overtone transition dipole moment as a function of the corresponding transition frequency, ω21, at equilibrium on the first excited state PES (blue). The rather dramatic enhancement of the transition dipole moments in the low-frequency region reflects the great sensitivity of the adiabatic state wave functions to the configuration of the bath in this region and is consistent with the recent observation of very pronounced non-Condon effects in this system.37 C. One-Dimensional Spectra. The 1D IR spectra of the H-stretch calculated using the standard ground-state equilibrium dynamical approach (dashed line) and the nonequilibrium MQCL approach (solid line) are shown in Figure 4. Both spectra are seen to consist of three major bands centered at ∼200, ∼2150-2250, and ∼2500 cm-1, which can be assigned to the transition state, covalent, and ionic configurations, respectively.37 It should be noted that these bands cannot be resolved in the fundamental transition frequency distribution (see upper panel of Figure 3). This implies that the line shape in the highfrequency region is at least somewhat affected by motional narrowing. Furthermore, the lack of absorption in the intermediate-frequency region and the emergence of a low-frequency band
J. Phys. Chem. B, Vol. 113, No. 27, 2009 9285 result from weakening of the fundamental transition dipole moment in the intermediate region and a strong enhancement of the fundamental transition dipole moment in the lowfrequency region (see lower panel of Figure 3). Importantly, the bath configurations that give rise to the low-frequency band correspond to the transition state of the proton transfer reaction that interconverts between the ionic and covalent tautomers. Thus, the low-frequency band provides a rather unique and direct probe of transition state bath configurations (albeit in an inconveniently located spectral range). Although the spectra obtained via the two approaches are similar, there are still significant differences between them. More specifically, the ionic and covalent bands in the MQCL-based spectrum are red-shifted and broader and the transition-state band is more pronounced in comparison to the spectrum generated via the standard approach. These differences represent the signatures of solvation on the averaged PES during t1 (see eq 32), and can be traced back to the fact that the average PES tends to destabilize the ionic and covalent bath configurations and stabilize the transition-state configurations, thereby redshifting and broadening the ionic and covalent bands, while sharpening and intensifying the transition-state band at their expense. D. Two-Dimensional Spectra. The 2D IR spectra of the H-stretch calculated at four different values of t2 using the standard ground-state equilibrium approach (top panels) and nonequilibrium MQCL approach (bottom panels) are shown in Figure 5. A detailed discussion of the different spectral features of the 2D spectra as obtained within the standard approach can be found in ref 37. We therefore focus our attention here on the differences between the 2D spectra predicted by the standard and MQCL approaches, which correspond to signatures of nonequilibrium processes accounted for within the MQCL approach, but not in the standard approach. At t2 ) 0 fs, both approaches predict 2D spectra which consist of three positive diagonal peaks that correspond to the ionic, covalent, and transition state bands in the absorption spectrum (see Figure 4), and negative off-diagonal peaks which represent contributions from the Liouville pathway associated with the overtone transition, namely 00 f 01 f 11 f 12 (see eq 34). The relative similarity of the spectra at t2 ) 0 can be traced back to the fact that the extent of solvation on the average 01 and 12 PESs during the coherence periods t1 and t3 is limited by the short dephasing time (∼150 fs). Nevertheless, the effect of nonequilibrium dynamics on the average PES during these time intervals is discernible even at t2 ) 0. More specifically, nonequilibrium dynamics on the average 01 PES intensifies the diagonal transition-state peak at the expense of the ionic and covalent peaks. The latter peaks also broaden and red-shift as a result of the fact that the 01 PES drives the system away from the ionic and covalent configurations and toward the transitionstate configurations. The negative feature at the lower right corner of the spectrum also broadens, red-shifting along the ω1 axis and blue-shifting along the ω3 axis. This is consistent with the fact that the fundamental frequency, ω10, decreases during solvation on the average 01 PES within the time interval t1, while the overtone frequency, ω21, increases during solvation on the average 12 PES within the time interval t3 (see Figure 1). Finally, the negative peak at (ω1,ω3) ∼ (200,1500), which is only discernible in the MQCL 2D spectrum, can be traced back to the fact that solvation on the 01 PES during t1 drives the system toward transition-state configurations, thereby giving rise to absorption at the overtone frequency of 1500 cm-1 during the time interval t3.
9286
J. Phys. Chem. B, Vol. 113, No. 27, 2009
Hanna and Geva
Figure 5. 2D IR spectra of the H-stretch as obtained via the standard ground-state equilibrium dynamical approach (upper panels) and the nonequilibrium MQCL approach (lower panels) for the following values of t2: 0, 125, and 250 fs and τeq.
The deviations between the spectra obtained via the two approaches become significantly larger with increasing t2. Within MQCL, starting at equilibrium on the ground-state PES, the first pulse is most likely to excite the system within the ionic, covalent, and transition-state bands. As discussed above, the state of the system changes relatively little during the short time it spends on the average 01 PES between the first and second pulse (t1). The second pulse then either returns the system to the ground-state PES or places it on the first excited state PES. Thus, if the system is returned to the ground state at t1, the ensuing nonequilibrium dynamics during t2 will drive the system away from transition-state configurations toward ionic and covalent configurations. In contrast, if the system is promoted to the first excited PES, the ensuing nonequilibrium dynamics during t2 drives the system away from the ionic and covalent configurations and toward the transition-state configurations. Since the intermediate spectral range between the high-frequency ionic/covalent bands and low-frequency transition-state band is significantly less photoactive, solvation actually weakens the signal at intermediate times (t2 ) 125 fs in Figure 5). However, a revival of the signal occurs when the system reaches either the spectral region of the ionic/covalent bands, in the case of ground state solvation, or the transition-state band, in the case of first excited state solvation (t2 ) 250 fs in Figure 5). At this point, the resulting MQCL-based 2D spectrum looks rather different in comparison to the standard one, as it starts taking on its asymptotic post-solvation shape (eq 36) and t2 ) τeq in Figure 5). More specifically, the emerging 2D spectrum consists of three spectral features: (1) a diagonal positive feature A A (ω1) × I01 (ω3) and is similar to the which corresponds to I01 A A j jI01 (ω1) × I01(ω3) spectral feature in the standard 2D spectrum; (2) an elongated positive spectral feature parallel to the ω1 axis and centered around ω3 ) 200 cm-1, which corresponds to A E (ω1) × I01 I01 (ω3)and replaces the spectral feature due to the A A j (ω3) term in the standard 2D spectrum; (3) second I01(ω1) × jI01 an elongated negative spectral feature parallel to the ω1 axis and centered around ω3 ) 1500 cm-1, which corresponds to A A A A (ω1) × I12 (ω3) and replaces the term -Ij01 (ω1) × jI12 (ω3) -I01 term in the standard 2D spectrum. It should be noted that excited-state solvation ends with the system in a region where nonadiabatic transitions are likely. The spectral effects of such
transitions are not accounted for in the present formalism, but will be addressed in a separate paper. V. Concluding Remarks In this paper, we have derived mixed quantum-classical expressions for the linear and third-order ORFs based on MQCL dynamics. On time scales where nonadiabatic transitions can be neglected, the expressions were found to be equivalent to those obtained via the LSC method. Within this approach, which is based on taking the classical limit in a rigorous manner, the temporal behavior of the ORFs reflects the nonequilibrium classical dynamics of the photoinactive DOF in response to a sequence of light-matter interactions. More specifically, the photoinactive DOF alternate between the adiabatic PESs and their averages, as prescribed by the Liouville pathways that are dictated by each ORF. This should be contrasted with the standard approach, which is based on taking the classical limit in an ad hoc manner, thereby giving rise to ORFs that reflect equilibrium dynamics that takes place solely on the groundstate PES. The feasibility and usefulness of the MQCL approach were demonstrated on the Azzouz-Borgis model of a moderately strong H-bonded complex dissolved in a dipolar liquid. The nonequilibrium nature of the underlying dynamics was shown to have a pronounced effect on both the 1D and 2D IR spectra of the H-stretch in this rather nontrivial model system. This is particularly true in the case of the 2D spectra at long values of the waiting time, t, where the MQCL-based 2D spectra were found to be qualitatively different from those predicted by the standard approach. The MQCL approach described in this paper will be particularly useful for complex systems, such as liquid solutions and biomolecular assemblies, in cases where there are reasons to expect that the excited-state PESs may be significantly different from the ground-state PES. For example, recent pump-probe and 2D IR spectra of methanol oligomers in liquid CCl4 have revealed that H-bonds break at a faster rate during the nonequilibrium solvation process that follows photoexcitation of the OH stretch.45,46,105–108 It was also found that strong H-bonds on the red side of the absorption band break faster
Mixed Quantum-Classical Liouville Method than weak H-bonds on the blue side.108 While the standard approach would not be able to account for such effects, the MQCL approach may. Work toward the highly desirable extension of the MQCL approach to include the effect of nonadiabatic transitions and its application to systems like methanol/CCl4 is currently underway in our group and will be reported in future publications. Acknowledgment. This work was supported by the National Science Foundation through grants CHE-0809506 and PHY0114336. References and Notes (1) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford: New York, 1995. (2) Fleming, G. R.; Cho, M. Annu. ReV. Phys. Chem. 1996, 47, 109. (3) Joo, T.; Jia, Y.; Yu, J.; Lang, M. J.; Fleming, G. R. J. Chem. Phys. 1996, 104, 6089. (4) Jonas, D. M. Annu. ReV. Phys. Chem. 2003, 54, 425. (5) Khalil, M.; Demirodoven, N.; Tokmakoff, A. J. Phys. Chem. A 2003, 107, 5258. (6) Cho, M. Chem. ReV. 2008, 108, 1331. (7) Skinner, J. L. J. Chem. Phys. 1982, 77, 3398. (8) Walsh, A. M.; Loring, R. F. Chem. Phys. Lett. 1991, 186, 77. (9) Mukamel, S. J. Chem. Phys. 1982, 77, 173. (10) Shemetulskis, N. E.; Loring, R. F. J. Chem. Phys. 1992, 97, 1217. (11) Bursulaya, B. D.; Kim, H. J. J. Phys. Chem. 1996, 100, 16451. (12) Fried, L. E.; Bernstein, N. B.; Mukamel, S. Phys. ReV. Lett. 1992, 68, 1842. (13) Saven, J. G.; Skinner, J. L. J. Chem. Phys. 1993, 99 (6), 4391. (14) Stephens, M. D.; Saven, J. G.; Skinner, J. L. J. Chem. Phys. 1997, 106, 2129. (15) Everitt, K. F.; Geva, E.; Skinner, J. L. J. Chem. Phys. 2001, 114, 1326. (16) Egorov, S. A.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1998, 108, 1407. (17) Shi, Q.; Geva, E. J. Chem. Phys. 2005, 122, 064506. (18) Ka, B. J.; Geva, E. J. Chem. Phys. 2006, 125, 214501. (19) Shi, Q.; Geva, E. J. Chem. Phys. 2008, 129, 124505. (20) Martens, C. C.; Fang, J. J. Chem. Phys. 1997, 106, 4918. (21) Donoso, A.; Martens, C. C. J. Chem. Phys. 1998, 102, 4291. (22) Donoso, A.; Martens, C. C. J. Chem. Phys. 2000, 112, 3980. (23) Donoso, A.; Kohen, D.; Martens, C. C. J. Chem. Phys. 2000, 112, 7345. (24) Schu¨tte, C. Konard-Zuse-Zentrum fu¨r informationstechnik Berlin, June 1999; pages Preprint SC 99-10. (25) Santer, M.; Manthe, U.; Stock, G. J. Chem. Phys. 2001, 114, 2001. (26) Kapral, R.; Ciccotti, G. J. Chem. Phys. 1999, 110, 8919. (27) Nielsen, S.; Kapral, R.; Ciccotti, G. J. Chem. Phys. 2000, 112, 6543. (28) Wan, C.; Schofield, J. J. Chem. Phys. 2000, 113, 7047. (29) Wan, C.; Schofield, J. J. Chem. Phys. 2000, 112, 4447. (30) Kapral, R. Annu. ReV. Phys. Chem. 2006, 57, 129. (31) Shi, Q.; Geva, E. J. Chem. Phys. 2004, 121, 3393. (32) Hanna, G.; Kapral, R. J. Chem. Phys. 2005, 122, 244505. (33) Hanna, G.; Kapral, R. Acc. Chem. Res. 2006, 39, 21. (34) Hanna, G.; Geva, E. J. Phys. Chem. B 2008, 112, 4048–4058. (35) Hanna, G.; Kapral, R. J. Chem. Phys. 2008, 128, 164520. (36) Azzouz, H.; Borgis, D. J. Chem. Phys. 1993, 98, 7361. (37) Hanna, G.; Geva, E. J. Phys. Chem. B 2008, 112, 12991–13004. (38) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194521. (39) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194522. (40) Gundogdu, K.; Bandaria, J.; Nydegger, M.; Rock, W.; Cheatum, C. M. J. Chem. Phys. 2007, 127, 044501. (41) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698. (42) Heyne, K.; Huse, N.; Nibbering, E. T. J.; Elsaesser, T. Chem. Phys. Lett. 2003, 382, 19. (43) Heyne, K.; Huse, N.; Dreyer, J.; Nibbering, E. T. J.; Elsaesser, T.; Mukamel, S. J. Chem. Phys. 2004, 121, 902. (44) Cowan, M. L.; Bruner, B. D.; Huse, N.; Dwyer, J. R.; Chugh, B.; Nibbering, E. T. J.; Elsaesser, T.; Miller, R. J. D. Nature 2005, 434, 199. (45) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Goun, A.; Fayer, M. D. Phys. ReV. Lett. 2003, 91, 23742. (46) Asbury, J. B.; Steinel, T.; Fayer, M. D. J. Phys. Chem. B 2004, 108, 6544–6554.
J. Phys. Chem. B, Vol. 113, No. 27, 2009 9287 (47) Steinel, T.; Asbury, J. B.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. Chem. Phys. Lett. 2004, 386, 295. (48) Kolano, C.; Helbing, J.; Kozinski, M.; Sander, W.; Hamm, P. Nature 2006, 444, 469. (49) Woutersen, S.; Mu, Y.; Stock, G.; Hamm, P. Chem. Phys. 2004, 266, 137. (50) Huse, N.; Bruner, B. D.; Cowan, M. L.; Dreyer, J.; Nibbering, E. T. J.; Miller, R. J. D.; Elsaesser, T. Phys. ReV. Lett. 2005, 95, 147402. (51) Park, S.; Fayer, M. D. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 16731. (52) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 5827. (53) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 264. (54) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 3840. (55) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2005, 123, 044513. (56) Diraison, M.; Guissani, Y.; Licknam, J. C.; Bratos, S. Chem. Phys. Lett. 1996, 258, 348. (57) Rey, R.; Moller, K. B.; Hynes, J. T. J. Phys. Chem. A 2002, 106, 11993. (58) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 8847. (59) Piryatinski, A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 9664. (60) Piryatinski, A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 9672. (61) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 1623. (62) Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2004, 120, 8107–8117. (63) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2004, 121, 8887–8896. (64) Li, S.; Schmidt, J. R.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2006, 124, 204110. (65) Scheurer, C.; Piryatinski, A.; Mukamel, S. J. Am. Chem. Soc. 2001, 123, 3114. (66) Scheurer, C.; Mukamel, S. J. Chem. Phys. 2002, 116, 6803. (67) Venkatramani, R.; Mukamel, S. J. Chem. Phys. 2002, 117, 11089. (68) Moran, A. M.; Park, S.-M.; Mukamel, S. J. Chem. Phys. 2003, 118, 9971. (69) Hayashi, T.; Mukamel, S. J. Phys. Chem. A 2003, 107, 9113. (70) Zhuang, W.; Abramavicius, D.; Hayashi, T.; Mukamel, S. J. Phys. Chem. B 2006, 110, 3362. (71) Zhuang, W.; Abramavicius, D.; Voronine, D. V.; Mukamel, S. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14233. (72) Sˇanda, F.; Mukamel, S. J. Chem. Phys. 2006, 124, 124103. (73) Hayashi, T.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 9747. (74) Hayashi, T.; la Cour Jansen, T.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 64. (75) la Cour Jansen, T.; Hayashi, T.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2005, 123, 114504. (76) Hayashi, T.; Mukamel, S. J. Chem. Phys. 2006, 125, 194510. (77) Williams, R. B.; Loring, R. F. J. Chem. Phys. 2000, 113, 10651. (78) Williams, R. B.; Loring, R. F. J. Chem. Phys. 2000, 113, 1932. (79) Williams, R. B.; Loring, R. F. Chem. Phys. 2001, 266, 167. (80) Akiyama, R.; Loring, R. F. J. Chem. Phys. 2002, 116, 4655. (81) Akiyama, R.; Loring, R. F. J. Phys. Chem. A 2003, 107, 8024. (82) Noid, W. G.; Loring, R. F. J. Chem. Phys. 2004, 121, 7057. (83) Noid, W. G.; Ezra, G. S.; Loring, R. F. J. Chem. Phys. 2004, 120, 1491. (84) Noid, W. G.; Loring, R. F. J. Chem. Phys. 2005, 122, 174507. (85) Cho, M. J. Chem. Phys. 2001, 115, 4424. (86) Kwac, K.; Cho, M. J. Chem. Phys. 2003, 119, 2256. (87) Kwac, K.; Lee, H.; Cho, M. J. Chem. Phys. 2004, 120, 1477. (88) Sung, J.; Silbey, R. J. J. Chem. Phys. 2001, 115, 9266. (89) Sung, J.; Silbey, R. J. J. Chem. Phys. 2003, 118, 2443. (90) Kryvohuz, M.; Cao, J. Phys. ReV. Lett. 2005, 95, 180405. (91) Kryvohuz, M.; Cao, J. Phys. ReV. Lett. 2006, 96, 030403. (92) Egorova, D.; Gelin, M. F.; Domcke, W. J. Chem. Phys. 2007, 126, 074314. (93) Steffen, T.; Fourkas, J. T.; Duppen, K. J. Chem. Phys. 1996, 105, 7364. (94) Piryatinski, A.; Tretiak, S.; Chernyak, V.; Mukamel, S. J. Raman Spectrosc. 2000, 31, 125. (95) Mukamel, S.; Abramavicius, D. Chem. ReV. 2004, 104, 2073. (96) Zhang, W. M.; Chernyak, V.; Mukamel, S. J. Chem. Phys. 1999, 110, 5011. (97) Dreyer, J. J. Quantum Chem. 2005, 104, 782. (98) DeVane, R.; Space, B.; Perry, A.; Neipert, C.; Ridley, C.; Keyes, T. J. Chem. Phys. 2004, 121, 3688. (99) Lazonder, K.; Pshenichnikov, M. S.; Wiersma, D. A. Opt. Lett. 2006, 31, 3354. (100) Hamm, P. J. Chem. Phys. 2006, 124, 124506. (101) Ishizaki, A.; Tanimura, Y. J. Chem. Phys. 2006, 125, 084501. (102) Kato, T.; Tanimura, Y. J. Chem. Phys. 2004, 120, 260.
9288
J. Phys. Chem. B, Vol. 113, No. 27, 2009
(103) Wigner, E. Phys. ReV. 1932, 40, 749. (104) Hanna, G.; Geva, E. J. Phys. Chem. B 2008, 112, 15793. (105) Graener, H.; Ye, T. Q.; Laubereau, A. J. Chem. Phys. 1989, 90, 3413–3416. (106) Levinger, N. E.; Davis, P. H.; Fayer, M. D. J. Chem. Phys. 2001, 115, 9352–9360. (107) Gaffney, K. J.; Davis, P. H.; Piletic, I. R.; Levinger, N. E.; Fayer, M. D. J. Phys. Chem. A 2002, 106, 12012–12023.
Hanna and Geva (108) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Fayer, M. D. J. Chem. Phys. 2003, 119, 12981–12997. (109) Sepulveda, A.; Mukamel, S. J. Chem. Phys. 1995, 102, 9327. (110) Spencer, F.; Loring, R. F. J. Chem. Phys. 1996, 105, 6596. (111) Pentidis, A.; Loring, R. F. Chem. Phys. Lett. 1998, 287, 217.
JP902797Z
