Computational Study of the One and Two Dimensional Infrared

Sep 18, 2008 - Gabriel Hanna* and Eitan Geva*. Department of Chemistry and FOCUS center, UniVersity of Michigan, Ann Arbor, Michigan 48109-1055...
0 downloads 0 Views 4MB Size
J. Phys. Chem. B 2008, 112, 12991–13004

12991

Computational Study of the One and Two Dimensional Infrared Spectra of a Vibrational Mode Strongly Coupled to Its Environment: Beyond the Cumulant and Condon Approximations Gabriel Hanna* and Eitan Geva* Department of Chemistry and FOCUS center, UniVersity of Michigan, Ann Arbor, Michigan 48109-1055 ReceiVed: May 9, 2008; ReVised Manuscript ReceiVed: July 31, 2008

The effect of the commonly employed Condon and second-order cumulant approximations on one- and twodimensional infrared spectra is examined in the case of a vibrational mode which is strongly coupled to its environment. The analysis is performed within the context of the hydrogen stretch of a moderately strong hydrogen-bonded complex dissolved in a dipolar liquid. The IR spectra are calculated using an adiabatic mixed quantum-classical approach that treats the hydrogen quantum-mechanically, while the remaining degrees of freedom are treated classically. While the cumulant and Condon treatments are seen to produce extremely broad and rather structureless spectra, the non-Condon spectra are found to consist of several relatively narrow bands that can be traced back to subsets of bath configurations with large transition dipole moments. Thus, although the cumulant and Condon approximations can capture some general qualitative spectral trends and are able to reproduce some highly averaged quantities such as the photon-echo peak shift, they fail to reproduce many important features of the spectra. We show that the great sensitivity of the transition dipole moment to the bath configuration provides new means for decongesting the spectra, probing statistically unfavorable bath configurations, and obtaining unique information regarding the dynamics of individual subsets of bath configurations and of the rates of transitions between them. I. Introduction The Condon and second-order cumulant (2OC) approximations are often invoked in modeling infrared (IR) spectra of vibrational modes embedded in condensed phase hosts.1 Within the Condon approximation, one assumes that the transition dipole moments of the vibrational modes are independent of the configuration of the host. This assumption can be justified in the limit of weak coupling between the vibrational mode and its environment since in this case the vibrational wave functions are insensitive to the host configuration. However, this is not the case for strongly coupled systems where the vibrational wave functions are expected to be influenced by the host configuration. The Condon approximation is often accompanied by the 2OC approximation, which is based on the additional assumption that when performing the cumulant expansion of the optical response functions (ORFs), cumulants of third- and higher-order can be neglected. Within this Gaussian approximation, one can express the ORFs in terms of the average transition frequencies and two-time frequency-frequency correlation functions. This approximation happens to be exact when the adiabatic potentials that describe the dynamics of the environment at different vibrational states are given by shifted harmonic surfaces. It is also believed to be justified in systems where the coupling between the vibrational mode and its environment arises from a large number of independent interactions. However, the validity of this approximation becomes highly questionable in the presence of strong coupling between the vibrational mode and its environment, since anharmonic effects and correlated dynamics are expected to be prevalent. In this paper, we examine the spectral signatures of breaking down the Condon and 2OC approximations in the case of a * Corresponding authors. E-mail: [email protected]; [email protected].

vibrational mode strongly coupled to its environment. We do so within the context of the Azzouz-Borgis (AB) model of a moderately strong hydrogen-bonded (H-bonded) phenol-amine complex in liquid methyl-chloride,2 which, as we will show below, is highly suitable for this purpose. H-bonded complexes play an important role in chemistry and biology.3-6 The high sensitivity of the vibrational frequency of H-bonds to the structure and dynamics of their environment has established IR spectroscopy as a powerful tool in H-bond science.4,5 More recently, it has become possible to probe vibrational dynamics via two-dimensional IR (2DIR) techniques in a variety of different systems,7-45 including H-bonded systems.46-60 These experimental advances were accompanied by substantial theoretical and computational studies aimed at understanding the signature of the underlying molecular structure and dynamics on the corresponding IR spectra.1,51,61-111 However, relatively little attention has been given to the unique features of IR spectra in the presence of strong coupling between the IR-active vibrational mode and its environment and, in particular, to the validity of the Condon and 2OC approximations and the implications of avoiding them. The main goal of this paper is to provide a systematic examination of these questions in the context of an atomistically detailed model of a H-bonded complex in liquid solution. The remainder of this paper is organized as follows. The model of an H-bonded complex in a dipolar liquid solution is described in section 2. The theoretical framework for calculating 1D and 2D IR spectra and related observables via different approximate schemes is outlined in section 3. The simulation techniques are described in section 4. The results are presented and discussed in section 5. Finally, the main conclusions are summarized in section 6.

10.1021/jp804120v CCC: $40.75  2008 American Chemical Society Published on Web 09/18/2008

12992 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Hanna and Geva

ˆ (qˆ, pˆ, Q, P) ) KB(P) + KˆP(pˆ) + Vˆ(Q, qˆ) H

II. Model In this section, we provide a short outline of the AB model as used in the present study (a more detailed description is available in refs 2 and 112-114). Within this model, it is assumed that the hydrogen moves along a one-dimensional (1D) axis that connects the donor and acceptor. The proton donor, A, and acceptor, B, are modeled as single particles and parametrized to represent phenol and trimethylamine, respectively. The functional form of the intramolecular potential surface as a function of the A-B and A-H distances is adopted from ref 115. Polarizability effects are accounted for by assuming that the charges on A and B are explicitly dependent on the position of the proton. Increasing the equilibrium A-B distance modifies the potential along the proton coordinate from a single well at short A-B distances to a double well at large A-B distances. The double-well form of the potential gives rise to tautomeric equilibrium between a covalent form and an ionic form of the complex: AH-B h A--H+B. In the covalent form, the proton is attached to the phenol (A) and the complex has a small (2.5 D) dipole moment due to the polarity of the A-H bond. In the ionic form, the proton is attached to trimethylamine (B), and the resulting charge separation gives rise to a significantly larger (10.5 D) dipole moment. In this study, we set the equilibrium A-B distance to 2.7 Å, which corresponds to a moderately strong H-bond (i.e., a double-well profile with a relatively low barrier along the proton coordinate). The solvent consists of a liquid of methyl-chloride molecules, which are modeled as rigid dipolar diatomic molecules. The complex-solvent and solventsolvent interactions are modeled in terms of site-site LennardJones (LJ) and Coulomb interactions, and the corresponding force field parameters are adopted from ref 116. Although the covalent form, AH-B, is favored in the absence of the solvent, the ionic form, A-H+B, becomes more stable when the complex is dissolved in a dipolar solvent. Most previous studies based on the AB model have focused on calculating the rates for proton transfer2,112,114,116-123 and H-stretch vibrational energy relaxation.113,124-126 A calculation of the absorption IR spectrum of the H-stretch for a symmetrical H-bonded complex model of the AB type, in the slowmodulation limit, was reported in ref 127. A related study by the same group focused on the absorption IR spectrum of the hydroxyl stretch of H-bonded methanol dimer and trimer solvated in liquid CCl4.128 However, it should be noted that the solvent-solute interactions are weaker in this system and that the authors of ref 128 only considered the absorption IR spectrum within the framework of the Condon approximation. More recently, Skinner and co-workers have modeled 1D and 2D IR spectra of the H-stretch in water and its isotopomers via a variety of methods.42,55,61-64,67-72 The study in ref 64, in particular, provided a clear demonstration of the limitations of the Condon and 2OC approximations in modeling 1D and 2D IR spectra of the hydrogen stretch in water. This study also provided the impetus for examining those effects in greater detail and in other model systems. III. Theory A. Adiabatic Mixed Quantum-Classical Treatment. A fully quantum-mechanical treatment of the AB model is not feasible. However, one may take advantage of the light mass of the proton relative to the other particles in order to employ a mixed quantum-classical treatment, where the hydrogen is treated quantum-mechanically while the solvent and the A and B groups are treated in a classical-like manner. The mixed quantum-classical Hamiltonian of the overall system is given by

(1)

Here, as in the rest of the paper, vectors are represented by boldface symbols (e.g., A) while operators are capped (e.g., Aˆ). P ≡ (P1, P2, . . .) and Q ≡ (Q1, Q2, . . .) are the coordinates and momenta of the solvent and A and B groups (collectively referred to as the bath degrees of freedom (DOF) from now on), qˆ and pˆ are the coordinate and momentum operators of the proton, KB(P) is the kinetic energy of the bath DOF, KˆP(pˆ) is the kinetic energy operator of the proton, and Vˆ(Q, qˆ) is the overall potential energy operator. The above-mentioned mass scale separation also implies that one may assume that the vibrational energy levels of the proton follow the bath DOF adiabatically. It is therefore reasonable to identify the vibrational energy levels and wave functions of the proton stretch with the eigenvalues and eigenfunctions of ˆ P(qˆ, pˆ, Q) ≡ KˆP(pˆ) + Vˆ(Q, qˆ): the protonic Hamiltonian, H

ˆ P|j;Q 〉 ) Ej(Q)|j;Q〉 H

(2)

Here, {| j; Q〉} are the adiabatic vibrational states and {Ej(Q)} are the corresponding vibrational energy levels (both explicitly dependent on the coordinates of the bath DOF). B. One-Dimensional IR Spectroscopy. Within the standard perturbative treatment,1 the 1D IR spectrum is given by

I(ω) ) Re

∫0∞ dt J(t)eiωt

(3)

where J(t) is the linear ORF. In this paper, we assume that J(t) is given by the following commonly used expression for the classical limit of J(t):1,64,129-131

J(t) )

∫ dQ0 ∫ dP0

exp[-βH0(Q0, P0)] V01(t)V10(0) × Z0 exp[-i

≡ 〈V01(t)V10(0) exp[-i

∫0t dτ ω10(τ)]

∫0t dτ ω10(τ)]〉

(4)

Here, H0(Q0, P0) ) KB(P0) + E0(Q0) and Z0 ) ∫dQ0 ∫dP0 exp[-βH0(Q0, P0)] are the classical bath Hamiltonian and partition function that correspond to the ground vibrational state; Vij(t) ) 〈i; Qt|qˆ|j; Qt〉 is the transition dipole moment; pωij(t) ) Ei(Qt)-Ej(Qt) is the energy gap between vibrational states i and j; Qt represents the trajectory of the bath DOF obtained via classical propagation on the ground state potential surface, E0(Q), with the initial conditions (Q0, P0). It should be noted that J(t) actually corresponds to the diagonal element of the linear optical response tensor along the body-fixed proton displacement coordinate. Since the rotational relaxation in this model system is an order of magnitude slower than the decay time of J(t), one can obtain the linear ORF in the laboratory frame from that in the body-fixed frame simply by multiplying it by 1/3.132 One limit of interest corresponds to the case where ω10(t) and V10(t) change on a time scale which is slow compared to the time scale on which J(t) decays to zero. The linear ORF in this slow modulation (or inhomogeneous) limit is given by

Jinh(t) )

∫ dQ0 ∫ dP0

exp[-βH0(Q0, P0)] |V10(Q0)|2 × Z0 exp[-iω10(Q0)t]〉

≡ 〈|V10(Q0)|2 exp[-iω10(Q0)t]〉 The corresponding inhomogeneous line shape is given by

(5)

1D and 2D IR Spectra of a Vibrational Mode

J. Phys. Chem. B, Vol. 112, No. 41, 2008 12993

Iinh(ω) ) π〈|V10(Q0)|2δ(ω - ω10(Q0))〉

(6)

It should be noted that Iinh(ω) is proportional to the distribution of ω10 weighted by the modulus square of the transition dipole moment, |V10(Q0)|2. It is important to point out that the sensitivity of the linear ORF in eq 4 to the bath configuration arises from its dependence on both vibrational energy levels (through the transition frequencies) and wave functions (through the transition dipole moments). The Condon approximation amounts to neglecting the fluctuations in the transition dipole moments due to the dependence of the wave functions on the bath configuration. It corresponds to assuming that the transition dipole moments are constant and can therefore be substituted by their average value:

Vij(t) ≈ 〈Vij 〉

(7)

The resulting approximation for the linear ORF is given by

∫ dQ0∫ dP0

JC(t) ) |〈V01 〉 |2

exp[-βH0(Q0, P0)] × Z0 exp[-i

≡ |〈V10 〉 |2〈exp[-i

Rr(t3, t2, t1) ) 2〈V01(0)V10(t1)V10(t1 + t2)V01(t1 + t2 + t3) × exp[i

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω10(τ)]〉 1

∫0t dτω10(τ)]〉

(8)

The 2OC approximation corresponds to truncating the cumulant expansion of 〈exp[-i∫t0dτ ω10(Qτ)]〉 at second-order to yield the following expression for the linear ORF:1

J2OC(t) ) |〈V01 〉 |2 exp[-i〈ω10 〉 t - g10,10(t)]

exp[i

3

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω21(τ)]〉 (12) 1

1

1

2

3

2

respectively. It should be noted that the first terms on the righthand side (RHS) of eqs 11 and 12 correspond to contributions from Liouville pathways1 that involve only the fundamental transition between the ground and first excited state, while the second terms correspond to contributions from Liouville pathways that also involve the overtone transition between the first and second excited states. It should also be noted that these two terms have opposite signs. Finally, the 2DIR spectrum is obtained by combining the double Fourier-Laplace transforms of Rnr(t3, t2, t1) and Rr(t3, t2, t1) with respect to t1 and t3 as follows18,64,133

∫0∞ dt1∫0∞ dt3 {ei(ω t +ω t )Rnr(t3, t2, t1) + 11

33

ei(-ω1t1+ω3t3)Rr(t3, t2, t1)} (13) At t2 longer than the correlation time of the bath, the rephasing and nonrephasing signals become uncorrelated, so that

Rnr(t3, t2, t1) f 2〈V10(0)V01(t1) exp[-i 〈V10(0)V01(t3) exp[-i exp[-i

(9)

where

∫0t dτ ω10(τ)]〉 × 1

∫0t dτ ω10(τ)]〉 - 〈V10(0)V01(t1) × 3

∫0t dτ ω10(τ)]〉〈V21(0)V12(t3) × t exp[-i∫0 dτ ω21(τ)] 〉 (14) 1

3

∫0t dτ (t - τ)[〈ωij(τ)ωkl(0) 〉 -〈ωij 〉 〈ωkl〉] t ≡ ∫0 dτ (t - τ)Cij,kl(τ)

gij,kl(t) )

and

(10)

Here, Cij,kl(τ) is an equilibrium two-time frequency-frequency correlation function. When ij ) kl, it reduces to the autocorrelation function of the transition frequency ωij. C. Two-Dimensional IR Spectroscopy. The measurement of 2DIR spectra 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 third-order 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. The commonly used expressions for the nonrephasing and rephasing signals in the classical limit (assuming weak and impulsive pulses and groundstate dynamics) are given by64

Rnr(t3, t2, t1) ) 2〈V10(0)V01(t1)V10(t1 + t2)V01(t1 + t2 + t3) × exp[-i

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω10(τ)]〉 1

1

1

2

3

Rr(t3, t2, t1) f 2〈V01(0)V10(t1) exp[i

exp[-i

1

2

3

3

(15) respectively. Thus, the 2DIR spectrum at long t2 consists of a positive diagonal feature which is proportional to I(ω1) × I(ω3) and a negative off-diagonal feature (if ω10(τ) and ω21(τ) are significantly different) which is proportional to I(ω1) × I′(ω3), where I(ω) is as in eq 3 and

I′(ω) ) Re

∫0∞ dt eiωt〈V21(t)V12(0) exp[-i∫0t dτ ω21(τ)]〉 (16)

Within the Condon approximation (see eq 7), the nonrephasing and rephasing signals are given by

RCnr(t3, t2, t1) ) 2|〈V10 〉 |4〈exp[-i i

3

1

and

1

2

3

2

〈 exp[-i

2

∫0t dτ ω10(τ) -

∫t t+t+t +t dτ ω10(τ)] 〉 -|〈V10 〉 |2|〈V21 〉 |2 × 1

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω21(τ)]〉 (11) 1

1

1

2

1

∫0t dτ ω10(τ)]〉 ×

∫0t dτ ω10(τ)]〉 - 〈V01(0)V10(t1) × t t exp[i∫0 dτ ω10(τ)]〉〈V21(0)V12(t3) exp[-i∫0 dτ ω21(τ)]〉

〈V10(0)V01(t3) exp[-i

〈V10(0)V01(t1)V21(t1 + t2)V12(t1 + t2 + t3) ×

and

2

2

〈V01(0)V10(t1)V21(t1 + t2)V12(t1 + t2 + t3) ×

I(ω1, t2, ω3) ≡ Re

∫0t dτω10(τ)]

1

1

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω21(τ)]〉 (17) 1

1

1

2

2

3

12994 J. Phys. Chem. B, Vol. 112, No. 41, 2008

RCr (t3, t2, t1) ) 2|〈V10 〉 |4〈exp[i i

Hanna and Geva

∫0t dτ ω10(τ) -

IV. Simulation Techniques

1

∫t t+t+t +t dτ ω10(τ)] 〉 -|〈V10 〉 |2|〈V21 〉 |2 × 1

1

2

3

2

〈exp[i

∫0t dτ ω10(τ) - i∫t t+t+t +t dτ ω21(τ)]〉 (18) 1

1

1

2

3

2

respectively. Within the 2OC approximation, the nonrephasing and rephasing signals are further simplified to give 4 R2OC nr (t3, t2, t1) ) 2|〈V10 〉 | exp[-i〈ω10 〉 (t1 + t3) -

g10,10(t1) - g10,10(t2) + g10,10(t3)+ g10,10(t1 + t2) + g10,10(t2 + t3) - g10,10(t1 + t2 + t3)] - |〈V10 〉 |2|〈V21 〉 |2 × exp[-i〈ω10 〉 t1 - i〈ω21 〉 t3 - g10,10(t1) - g10,21(t2) g21,21(t3)+ g10,21(t1 + t2) + g10,21(t2 + t3) g10,21(t1 + t2 + t3)] (19) and

R2OC (t3, t2, t1) ) 2|〈V10 〉 |4 exp[i〈ω10 〉 (t1 - t3) - g10,10(t1) + r g10,10(t2) - g10,10(t3) - g10,10(t1 + t2) - g10,10(t2 + t3) + g10,10(t1 + t2 + t3)] - |〈V10 〉 |2|〈V21 〉 |2 exp[i〈ω10 〉 t1 i〈ω21 〉 t3 - g10,10(t1) + g10,21(t2) - g21,21(t3) - g10,21(t1 + t2) g10,21(t2 + t3) + g10,21(t1 + t2 + t3)] (20) respectively. Another quantity of interest is the three-pulse photon-echo peak shift (3PEPS). It can be defined in terms of the homodynedetected and integrated (over t3) contribution to the rephasing signal that arises from Liouville pathways that only involve the fundamental transition:

I10 r (t2, t1) )

∫0∞ dt3 |R10r (t3, t2, t1)|2

(21)

The peak shift, t/1(t2), is then defined as the value of t1 at which Ir10 is maximized for a given value of t2. Using this definition, it can be shown that (within the Condon and 2OC approximations, while ignoring bath dynamics that takes place during t1 and t3 and assuming that t2 is long in comparison to the time scale on which C10,10(t) decays):134

C10,10(t2) ∼ t/1(t2)

(22)

Thus, the 3PEPS signal, t/1(t2), provides a convenient way for determining the long-time tail of the correlation function C10,10(t), provided that the above assumptions are satisfied. An alternative way for determining C10,10(t) is based on the short-time slope of the integrated and homodyne-detected threepulse photon-echo (S3PE):134

The simulation techniques used in this paper are similar to those described in refs 112 and 113 and will therefore be only briefly outlined below. Simulations were performed with a single AHB complex and 255 methyl-chloride molecules in a cubic simulation box with periodic boundary conditions at a temperature and density of 250 K and 0.012 Å-3, respectively. The bond length of the methyl chloride molecules was constrained at 1.781 Å using the SHAKE135 and RATTLE136 algorithms. 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 adiabatic surfaces were computed “on-the-fly” by diagonalizing the protonic Hamiltonian for different instantaneous bath configurations. To this end, we first constructed the protonic Hamiltonian matrix in terms of a basis that consisted of 12 harmonic oscillator basis functions with the first set of 6 basis functions centered at the ionic well of the gas-phase intramolecular potential surface and the other set of 6 basis functions centered at the covalent well. The matrix was then diagonalized by using the LAPACK DSYGV routine. The bath dynamics was simulated by solving the corresponding classical equations of motion on the ground-state potential surface via the velocity Verlet algorithm with a time step of 1 fs. Initial bath configurations were generated via a constant temperature MD simulation, using a Nose´-Hoover thermostat.142 The results reported below were obtained by averaging over trajectories of length 1-2 ns. Some of the results were decomposed into their covalent and ionic contributions by averaging over many segments of these trajectories where the system is solely in the covalent or ionic configuration, respectively. The criterion which determines whether the system is in the covalent or ionic configuration is described in the following section. Finally, the 1D and 2D Fourier transforms required for computing the 1D and 2D IR spectra, respectively, were carried out numerically via the FFT routine. The 1D spectrum was generated based on a 2048-point time grid with a time interval of 1 fs. The 2D spectrum was generated based on a 2048 × 2048 (t1, t3) time grid with a time interval of 1 fs. V. Results and Discussion A. 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 the first to second excited states of the H-stretch (see section III). The corresponding vibrational energy levels and wave functions are explicitly dependent on the bath coordinates, Q. A convenient collective solvent coordinate is provided by the solvent polarization, which is defined as the difference between the solvent electrical potentials at the minima of the covalent and ionic wells:143,144

∆E(Q) ) ∂ C10,10(t2) ∼ I10 |(t , t )| ∂t1 r 2 1 t1)0

(23)

The relationship in eq 23 is more accurate than that in eq 22 since it does not rely on assuming that t2 is long in comparison to the time scale on which C10,10(t) decays. It can therefore provide information on the behavior of C10,10(t) at short times. However, eq 23 still relies on the Condon and 2OC approximations and on ignoring bath dynamics within the time intervals t1 and t3.

(

∑ zi,ae |Qa1- s| - |Qa 1- s′| i,a

i

i

)

(24)

Here, zae (e ) 1.602 × 10-19 C) is the charge on atom a in solvent molecule i, and s and s′ correspond to the center of mass of the complex and a point displaced by 0.56Å relative to it toward the B group, respectively.112,114 Figure 1 shows the ground, first-excited, and second-excited free energy surfaces (FESs) as a function of the solvent polarization. The three FESs are clearly very different, which is a manifestation of the strong coupling between the vibrational

1D and 2D IR Spectra of a Vibrational Mode

J. Phys. Chem. B, Vol. 112, No. 41, 2008 12995

Figure 1. Ground (0, blue), first-excited (1, red), and second-excited (2, green) free energy surfaces as a function of the solvent polarization.

Figure 2. Equilibrium distributions of the transition frequencies ω10 and ω21.

mode and the bath. The ground-state FES has a double-well form while the excited-state FESs have a single-well form. In the ground state, the solvent is seen to favor configurations that correspond to either low polarizations or high polarizations, with a low likelihood for configurations with intermediate values of the polarization. This behavior is a signature of the abovementioned tautomeric equilibrium between the covalent form of the complex, which favors low-polarization configurations (∆E < 0.0138 eC/Å), and the ionic form of the complex, which favors high-polarization configurations (∆E > 0.0138 eC/Å). It should be noted that significant subpopulations of both forms coexist at equilibrium (the equilibrium composition is 65% ionic and 35% covalent). It has also been shown that proton transfer, which converts the complex between the covalent and ionic forms, is driven by the very same fluctuations in the solvent polarization (for this model, the covalent-to-ionic and ionic-tocovalent free energy barrier heights are given by 2.09kBT and 2.88kBT, respectively114).112,145 The fluctuations in the vibrational energy levels also modulate the vibrational transition frequencies. In Figure 2, we show the distributions of the fundamental and overtone transition frequencies (ω10 and ω21, respectively). The distributions are seen to be pronouncedly asymmetrical (with peaks at ∼2500 and ∼500 cm-1, respectively), very broad (over 1000 cm-1), and very different from each other, all of which are manifestations of the strong coupling between the H-stretch and the bath and signal the breakdown of the Condon and 2OC approximations. B. One-Dimensional Spectra. The 1D IR spectra of the H-stretch as obtained within the non-Condon, Condon, and 2OC treatments are shown in Figure 3. Also shown in this figure are the spectra in the limit of inhomogeneous broadening and the individual contributions from the ionic and covalent subpopulations, defined based on whether the solvent polarization is larger or smaller than its transition state value, respectively (see Figure 1). We start by considering the non-Condon 1D spectrum (see lower panel in Figure 3), which is seen to consist of three major

Figure 3. 1D IR spectra (black) as obtained within the 2OC, Condon, and non-Condon treatments (normalized so that the area under the curve is equal to 1). Also shown are the corresponding spectra in the limit of inhomogeneous broadening (blue) and the relative contributions from the ionic (green) and covalent (red) tautomers.

Figure 4. Transition dipole moments V10 (blue) and V21 (red) averaged over all the configurations that correspond to a given value of the transition frequency.

bands centered at ∼200, ∼2250, and ∼2500 cm-1. The two, partially overlapping, high-frequency bands can be assigned to the covalent (∼2250 cm-1) and ionic (∼2500 cm-1) subpopulations. It should be noted that these bands cannot be resolved in the inhomogeneous spectrum, which implies that the line shape in the high frequency region is at least somewhat affected by motional narrowing. As we will show below, the low frequency band results from a strong enhancement of the transition dipole moment in the low frequency region. Importantly, the bath configurations that give rise to those low transition frequencies 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 the bath configurations at the transition state. We next turn to the linear absorption spectrum as obtained within the framework of the Condon approximation (see middle panel in Figure 3). The main difference between the Condon and non-Condon absorption spectra is the appearance of a very broad shoulder in the frequency range ∼500-2000 cm-1 in the Condon line shape, which is missing from the non-Condon line shape. This range of transition frequencies corresponds to bath

12996 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Figure 5. (upper panels) Potential energy as a function of the proton displacement, q, and the ground (blue), first-excited (red), and secondexcited (green) energy levels of the proton at representative bath configurations that correspond to different transition frequencies. The bath configurations in panels A (∆E ) 0 eC/Å) and E (∆E ) 0.027 eC/Å) are in the vicinity of the covalent and ionic wells, respectively; that in panel C (∆E ) 0.0138 eC/Å) is in the transition state region; that in panel B (∆E ) 0.009 eC/Å) is intermediate between the covalent well and transition state; and that in panel D (∆E ) 0.018 eC/Å) is intermediate between the transition state and ionic well. (lower panels) Vibrational wave functions that correspond to the ground (blue), firstexcited (red), and second-excited (green) vibrational energy levels.

Figure 6. Autocorrelation function of the transition frequency between the ground and first-excited states, ω10 (black). Also shown are the same correlation functions when the A-B distance is held fixed (blue) and the contributions from the ionic (green) and covalent (red) subpopulations.

Figure 7. Autocorrelation function of the transition dipole moment V10 (black). Also shown are the same correlation functions when the A-B distance is held fixed (blue) and the contributions from the ionic (green) and covalent (red) subpopulations.

configurations that are intermediate between the stable ionic and covalent tautomers and the unstable transition state. The lack of this shoulder in the non-Condon spectrum signals the breakdown of the Condon approximation and can be traced back to the strong dependence of the fundamental transition dipole moment V10 on the bath configuration. Indeed, the transition dipole moment averaged over all the configurations that j 10(ω10), correspond to a given value of the transition frequency, V shows a wide minimum within the frequency range 500-2000 cm-1, with an increase by factors of ∼3 and ∼10 at the

Hanna and Geva frequency ranges ∼2000-2500 and ∼0-500 cm-1, respectively (see Figure 4). This behavior gives rise to a solVent-induced selection rule, which renders transitions unfaVorable within the intermediate frequency range ∼500-1500 cm-1. The three-band structure of the non-Condon absorption line shape shown in Figure 3 and, in particular, the existence of the low frequency transition state band are a direct consequence of this selection rule. Further insight into the origin of the non-Condon effects in this system can be obtained by examining the protonic potential energy and wave functions as a function of the proton displacement, q, at representative bath configurations that correspond to different values of ω10 (see upper and lower panels of Figure 5, respectively). To this end, we first consider bath configurations which correspond to the vicinity of the covalent and ionic minima (see panels A and E in Figure 5, respectively). In those cases, the potential energy as a function of the proton displacement has the form of a highly asymmetric double well. The two minima in those double-well profiles correspond to the proton displacements for the covalent and ionic species, respectively. As expected, the asymmetry is tilted strongly in favor of the covalent tautomer when the solvent polarization is small (panel A), or in favor of the ionic tautomer when the solvent polarization is large (panel E). Because of the pronounced asymmetry, the ground and first-excited energy levels are seen to lie below the minimum of the less stable species. As a result, the corresponding wave functions are localized in the well of the more stable species and therefore occupy the same region in space. This in turn leads to a sizable transition dipole, V10 ) 〈1; Q|q|0; Q〉. It should be noted that the bath configurations that correspond to those values of the solvent polarization also give rise to large transition frequencies. The two high frequency bands in the non-Condon spectrum should therefore be attributed to the following two factors: (1) The underlying bath configurations are rather favorable, since they correspond to the vicinity of the covalent and ionic minima. (2) The transition dipole moments in this frequency range are rather large. Next, consider a bath configuration in the vicinity of the transition state (see panel C in Figure 5). In this case, the potential energy as a function of the proton displacement has the form of a symmetrical double-well. The ground and firstexcited energy levels are almost degenerate and lie below the barrier top. Since the potential is symmetrical, the corresponding wave functions are delocalized over both wells and therefore once again occupy the same region in space. This results in sizable values of the transition dipole moment V10, which are in fact even larger than those in the high frequency region (see Figure 4). The existence of a low frequency band in the nonCondon spectrum, despite the fact that the corresponding bath configurations are unfavorable, can be traced back to this enhancement in the transition dipole moment V10 at transition state bath configurations. Finally, consider the bath configurations that give rise to intermediate transition frequencies (see panels B and D in Figure 5). In this case, the potential energy as a function of the proton displacement has the form of a moderately asymmetric doublewell. As a result, the ground and first-excited energy leVels lie below and aboVe the minimum that corresponds to the less stable species. The corresponding wave functions are therefore localized in different wells and therefore occupy different regions in space. The resulting transition dipole moments are therefore significantly smaller in comparison to the high and low frequency regions. Thus, the absence of the broad shoulder

1D and 2D IR Spectra of a Vibrational Mode

J. Phys. Chem. B, Vol. 112, No. 41, 2008 12997

Figure 8. 2DIR relaxation spectra of the H-stretch as obtained within the 2OC treatment, for the indicated values of the waiting time, t2. (upper panels) Overall spectra. (middle panels) Positive contributions from Liouville pathways that only involve the fundamental transition. (lower panels) Negative contributions from Liouville pathways that also involve the overtone transition.

Figure 9. Same as Figure 8 within the Condon treatment.

at the intermediate frequency region in the non-Condon spectrum results from the combinations of the following two factors: (1) These configurations are less favorable in comparison to the ones corresponding to the ionic and covalent wells. (2) They give rise to significantly smaller values of the transition dipole moment V10 in comparison to the solvent configurations in the transition state region. We next consider the absorption spectrum as obtained within the 2OC treatment (see upper panel in Figure 3). The 2OC

spectrum is seen to be inconsistent in almost every respect when compared to the linear spectra obtained by avoiding the 2OC approximation. More specifically, the 2OC line shape is inhomogeneously broadened and consists of a single Gaussian band whose peak location and width deviate by hundreds of wave-numbers in comparison to the non-2OC lineshapes. Also, although the separate contributions from the covalent and ionic subpopulations can be differentiated in terms of their peak locations (i.e., the covalent peak occurs at a slightly lower

12998 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Hanna and Geva

Figure 10. Same as Figure 8 within the non-Condon treatment.

Figure 11. Scatter plot of the correlation between the fundamental transition frequency, ω10, and the overtone transition frequency, ω21.

frequency than the ionic peak), they cannot be distinguished in the overall spectrum. These observations provide strong evidence for the non-Gaussian nature of the molecular dynamics in this strongly coupled anharmonic and highly correlated system. Despite the failure of the 2OC approximation in reproducing the absorption spectrum, the frequency-frequency correlation function still provides useful diagnostics into the nature of the underlying molecular dynamics. The autocorrelation of the transition frequency between the ground and first excited states, ω10, is shown in Figure 6. It is seen to exhibit damped oscillations with a ∼80 fs period, which is consistent with the period of the low frequency A-B stretch. Indeed, the oscillations disappear upon constraining the A-B distance to its equilibrium value and can therefore be attributed to a modulation of the H-stretch frequency via intramolecular coupling to the A-B stretch. Further insight can be obtained by considering the individual contributions to the frequency-frequency correlation function from the ionic and covalent tautomers. While the covalent contribution exhibits damped oscillatory behavior, those oscillations are hardly noticeable in the case of the ionic contribution. This difference in behavior can be attributed to the fact that the bath configurations favored by the ionic

tautomer are more rigid and therefore inhibit the vibrational dynamics of the A-B stretch. Finally, the autocorrelation function of the transition dipole moment V10 is shown in Figure 7. It exhibits oscillations and decays on a time scale of 400 fs, which is similar to that of the transition frequency autocorrelation function. Again, the oscillations disappear upon constraining the A-B distance to its equilibrium value and the individual contributions from the ionic and covalent subpopulations are seen to behave similarly to their counterparts in the case of the frequency autocorrelation function. However, the fact that there are noticeable oscillations in the ionic contribution (although less pronounced than those in the covalent one) implies that the transition dipole moment is more sensitive to fluctuations in the bath DOF than the frequency. These observations suggest that the transition frequencies and dipole moments are affected by the same bath DOF and that a self-consistent description must account for the fluctuations in both. C. Two-Dimensional Spectra. The 2DIR spectra of the H-stretch as obtained within the 2OC, Condon, and non-Condon treatments are shown in Figures 8-10, respectively. In those figures, we show the overall spectra (upper panels) as well as the individual contributions from Liouville pathways that only involve the fundamental transition (middle panels) and Liouville pathways that also involve the overtone transition (lower panels), at different values of the waiting time, t2, within the range 0-250 fs. We start by considering the 2D spectra obtained within the 2OC treatment (see Figure 8). The overall spectra are seen to consist of a positive diagonal feature and a negative off-diagonal feature. These features originate from the first and second terms, respectively, on the RHS of eqs 19 and 20. The fact that |〈V21〉| is significantly larger than |〈V10〉| (see Figures 4 and 5) makes the negative off-diagonal feature more prominent than the positive diagonal feature. At t2 ) 0, the diagonal feature is elongated along the diagonal axis, while the off-diagonal feature is elongated along an axis with a negative slope that crosses

1D and 2D IR Spectra of a Vibrational Mode

J. Phys. Chem. B, Vol. 112, No. 41, 2008 12999

Figure 12. Closeups on the positive contribution to the 2DIR relaxation spectra from Liouville pathways that only involve the fundamental transition in the specified frequency windows, as obtained via the non-Condon treatment at the specified values of the waiting time, t2.

the diagonal axis. The fact that the diagonal and off-diagonal features lie along axes with opposite slopes is consistent with the fact that the values of ω10 within the time period (0, t1) are correlated with the values of ω10 and anticorrelated with the values of ω21 within the subsequent time period (t1, t1 + t3). The observation that ω10 and ω21 are anticorrelated reflects the fact that the bath configurations that give rise to large values of ω10 also give rise to small values of ω21, and vice versa (see Figure 11). It should also be noted that although ω21 is smaller than ω10 at bath configurations in the vicinity of the ionic and covalent wells, the trend reverses as one approaches the transition state, where ω21 is actually larger than ω10. As a result, the diagonal and off-diagonal features cross within the frequency range ω1, ω3∼1000-1500 cm-1. The fact that |〈V21〉| is significantly larger than |〈V10〉| leads to the cancelation of the positive diagonal feature by the significantly larger negative offdiagonal feature in the overlap region. The behavior of the 2OC 2D spectrum at longer t2 is consistent with the loss of correlation with increasing t2 between the values of ω10 within the time period (0, t1) and the values of ω10 and ω21 within the time period (t1 + t2, t1 + t2 + t3). The loss of correlation is manifested by broadening of the diagonal and off-diagonal features along the ω3 axis and asymptotically approaching the uncorrelated 2D spectrum (see eqs 14 and 15). It should be noted that the negative off-diagonal feature remains dominant throughout and leads to cancelation of the positive diagonal feature in the overlap region. The 2D spectra obtained within the Condon treatment, Figure 9, are quantitatively different from those obtained within the 2OC approximation. The differences are consistent with the corresponding differences between the 2OC and Condon linear

spectra (see Figure 3). The 2OC predictions of the location and spectral widths differ by hundreds of wavenumbers from those obtained via the Condon approximation. Furthermore, the 2OC treatment is unable to capture such basic properties of the spectrum as its inherently asymmetrical shape and bimodal structure. Despite this, the Condon 2D spectra are seen to follow general trends which are qualitatively similar to those observed in the 2OC spectra. Namely, they consist of a positive diagonal feature and a negative off-diagonal feature that are elongated along axes with opposite slopes at t2 ) 0 and broaden as t2 increases. Here too, the negative off-diagonal feature is significantly more prominent due to the fact that |〈V21〉| is significantly larger than |〈V10〉| and therefore leads to the cancelation of the diagonal feature in the overlap region. The non-Condon 2D spectra are shown in Figure 10 at different values of t2 within the range 0-250 fs. In addition, we show in Figures 12 and 13 closeups of the positive and negative spectral features of the non-Condon spectra in certain frequency windows and over a range of t2 values that extends up to 2 ps. The non-Condon 2D spectra are seen to be distinctly different from those obtained within the Condon approximation, thereby demonstrating the importance of accounting for nonCondon effects in this system. The positive and negative features no longer overlap effectively due to the small value of the transition dipole moments in the overlap region. Thus, the lack of signal in the intermediate frequency region should be attributed to non-Condon effects, rather than to those two features canceling each other out. The greatly enhanced transition dipole moments in the transition state region are seen to give rise to a low frequency positive diagonal feature at t2 ) 0 (see Figure 12) that is centered

13000 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Hanna and Geva

Figure 13. Closeups on the negative contribution to the 2DIR relaxation spectra from Liouville pathways that also involve the overtone transition in the specified frequency windows, as obtained via the non-Condon treatment at the specified values of the waiting time, t2.

at (ω1,ω3) ∼ (200, 200)cm-1 and to a negative off-diagonal feature (see Figure 13) that is centered at (ω1, ω3) ∼ (200, 1500)cm-1. Those features are seen to diminish with increasing t2 (from 0-100 fs) due to spontaneous solvation processes into more favorable bath configurations that correspond to intermediate transition frequencies and vanishingly small values of the transition dipole moment. At even longer t2 (from 750 fs-2 ps), one observes the emergence of two new positive offdiagonal peaks that are centered at (ω1, ω3) ∼ (200, 2250) and ∼ (200, 2500)cm-1. The emergence of those peaks signals the realization of bath configurations associated with stable covalent and ionic tautomers, respectively, which can be associated with high transition frequencies and sizable transition dipole moments. Moreover, one observes the emergence of a negative off-diagonal peak at (ω1, ω3) ∼ (200, 500)cm-1 that can also be traced back to the emergence of stable ionic and covalent tautomers, although in this case one cannot resolve the individual contributions of those two species. The value of t2 at which these off-diagonal positive and negative peaks emerge provides a direct and rather unique probe of the rate of solvation from the transition state region to the stable product regions. Furthermore, the fact that the ionic and covalent products give rise to two spectrally distinct positive peaks implies that one can distinguish between the solvation pathways that lead to the different tautomers. Indeed, the solvation times obtained in this way are consistent with those reported in ref 113 (∼1 ps) for the same model system. Furthermore, the observation that the covalent peak emerges earlier than the ionic one is consistent with the previously reported observation that solvation into the covalent well occurs more rapidly.113 Finally, it should be noted that the emergence of the covalent and ionic off-diagonal

positive peaks at long t2 is accompanied by the emergence of two additional off-diagonal positive peaks at (ω1, ω3) ∼ (2250, 200) and ∼ (2500, 200)cm-1. These peaks are associated with the subset of trajectories that start out with solvent configurations that correspond to the stable covalent and ionic tautomers, respectively, and end up in the transition state region. Unfortunately, these peaks are masked in the overall 2D spectrum by the negative term which is dominant in this spectral range. The 2D spectra can also provide a powerful probe of the rate of the proton transfer reaction between the ionic and covalent tautomers. More specifically, the emergence of cross peaks in the high frequency positive diagonal feature at values of t2 between 5 and 7 ps can be attributed to trajectories that started ionic and ended up covalent, or vice versa (see Figure 14). It should be noted that those cross peaks cannot be resolved within the 2OC treatment and that the timescale on which those peaks emerge in the Condon and non-Condon spectra is consistent with the adiabatic rate constant of 0.158 ps-1 recently reported for this proton transfer reaction.114 Finally, we show in Figure 15 the autocorrelation function of the transition frequency ω10 as predicted via the 3PEPS and S3PE methods when applied to the simulated rephasing signals calculated via the Condon and non-Condon treatments. In the case of the Condon treatment (upper panel), one observes that both the 3PEPS and S3PE methods reproduce the overall time scale on which the correlation function decays rather well. Also, as expected, the agreement between the correlation function and that obtained via the 3PEPS method is rather poor at short t2 but improves with increasing t2. In contrast, the agreement between the correlation function obtained via the S3PE method and the actual one is seen to be rather good even at short t2.

1D and 2D IR Spectra of a Vibrational Mode

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13001

Figure 14. High frequency positive diagonal feature in the 2DIR spectra at long t2. The emergence of off-diagonal peaks at t2 ) 5-7 ps is the signature of proton transfer.

Figure 15. Autocorrelation function of ω10, C10,10(t) (black), and as obtained via the 3PEPS (red) and S3PE (green) methods within the Condon (upper panel) and non-Condon (lower panel) treatments. All functions were normalized as in ref 134.

Both methods are, however, unable to reproduce the oscillations in the correlation function. The reason for this can be traced back to the fact that the derivation of the relationships between the peak shift and the initial slope to the correlation function as given in eqs 22 and 23 relies on neglecting the dynamics on times scales comparable to t1 and t3.134 Since the period of these oscillations (∼80 fs) is within the range of t1 and t3 in this case, they are missed by both the 3PEPS and S3PE methods, which can only reproduce the overall envelope of the correlation function. In the case of the non-Condon treatment (lower panel), one observes that both the 3PEPS and S3PE methods do not reproduce the overall time scale on which the correlation function decays as well as in the Condon treatment. Also, the

agreement between the correlation function and those obtained via the 3PEPS and S3PE methods at short t2 is poor and improves with increasing t2 only in the case of the 3PEPS result. In contrast to the Condon results, the non-Condon results do exhibit oscillations, but they differ slightly in phase and amplitude from those in the correlation function. This difference is consistent with the slight deviations observed between the dynamical timescales of the transition frequencies and dipole moments (see Figures 6, 7, and 15). The fact that the 3PEPS and S3PE methods can reproduce the frequency correlation function rather well in the Condon treatment is somewhat surprising in light of the fact that the overall 2DIR spectrum is reproduced rather poorly within the 2OC approximation. However, it should be noted that the 3PEPS and S3PE methods are based on measuring the homodynedetected and t3-integrated rephasing signal, rather than the heterodyne-detected time-resolved signal underlying the measurement of the 2DIR spectrum. It therefore appears that the extensive averaging involved in those methods can therefore make them rather insensitive to the non-Gaussian attributes of the underlying dynamics. On the other hand, it is not surprising that within the non-Condon treatment the 3PEPS and S3PE methods do not reproduce the frequency correlation function as well as in the Condon treatment, since the non-Condon effects in this model are quite pronounced. Furthermore, the fact that oscillations appear in the non-Condon results suggests that, despite the extensive averaging involved in constructing the signal, the signal remains somewhat sensitive to the fluctuations in the transition dipole moments. This in turn suggests that the transition dipole moments, through their dependence on the wave functions, can sometimes probe elements of the bath dynamics that the transition frequencies are insensitive to. VI. Concluding Remarks In this paper, we have explored the rather unique signatures of non-Gaussian and non-Condon effects on the 1DIR and 2DIR spectra of a vibrational mode strongly coupled to its environment. Although the study was focused on a specific model system, namely the AB model of a moderately strong hydrogenbonded complex solvated in an aprotic dipolar liquid, we believe that the physical mechanisms underlying those effects and their impact on the spectra are not unique to this model system and

13002 J. Phys. Chem. B, Vol. 112, No. 41, 2008 will prove to be relevant in other strongly coupled systems. More specifically, one generally expects the vibrational energy levels and wave functions to be very sensitive to the bath configurations in strongly coupled systems, thereby giving rise to a wide range of transition frequencies and transition dipole moments. Thus, even though we are dealing with a single mode, one can expect spectra that consist of multiple bands scattered over a wide range of frequencies, which represent subsets of bath configurations that are either highly probable and/or associated with large transition dipole moments. For instance, this concept may prove useful in the context of the spectroscopy of transition states in systems which exhibit strong nonadiabatic coupling in the vicinity of a barrier-top. This situation should be contrasted with that in the case of a vibrational mode which is weakly coupled to a bath, where each mode gives rise to a single band whose width is determined by a relatively narrow distribution of transition frequencies and a limited range of spectral diffusion. Our main conclusion is that while the 2OC and Condon approximations may be able to capture some general qualitative trends and reproduce some highly averaged quantities such as the peak shift, they do rather poorly when it comes to reproducing some of the most interesting and informative features of the spectra. More specifically, whereas the 2OC and Condon spectra are seen to be incredibly broad and rather structureless, the corresponding non-Condon spectra are seen to consist of several relatively narrow bands that can be associated with subsets of bath configurations that give rise to large transition dipole moments. In some respects, the diminished signal from those bath configurations that give rise to vanishingly small transition dipole moments represents a disadvantage since it implies that one cannot probe those bath configurations directly. At the same time, the very same sensitivity of the transition dipole moment to the bath configuration represents an advantage since it leads to less congested spectra and makes it possible to probe statistically unfavorable bath configurations. In fact, due to its dependence on the wave functions, the transition dipole moment may sometimes be a more sensitive probe of bath dynamics than the transition frequency. For example, in the case of the AB model, the sensitivity of the transition dipole moment translates into spectra that consist of three distinct bands that represent the stable ionic and covalent tautomers and the transition state between them. Thus, one can obtain unique information regarding the rate of spectral diffusion within each subset of bath configurations and the rates of spontaneous solvation and chemical exchange processes that take the system between them. The work reported in this paper can be extended in several directions. First, a more accurate treatment of the long-range electrostatic interactions is desirable. In this context, it should be noted that a straightforward application of Ewald summations within the framework of a mixed-quantum classical simulation is complicated by the fact that one needs to deal with periodic replicas of the proton in calculating the proton-proton and proton-solvent Coulombic forces. Thus, implementing techniques that circumvent this problem, such as those proposed in refs 137-141 would be worthwhile. Another limitation of the treatment is associated with the underlying assumption that the bath dynamics occurs on the ground-state adiabatic potential energy surface. As noted by several authors, this choice is not unique.1,129-131,134,146-150 The main advantage of this popular choice is that the ORFs can be calculated from equilibrium MD simulations on the ground-state potential surface, while other choices require costly nonequilibrium simulations that start at

Hanna and Geva equilibrium on the ground-state potential surface but proceed, at least part of the time, on the excited-state surfaces. However, the effect of the choice of potential energy surface on the spectra is expected to be pronounced in the case of strong system-bath coupling, since the potential energy surfaces can be very sensitive to the vibrational state. This is particularly true for the dependence of the third-order ORFs on t2, since it is typically much longer than t1 and t3. More specifically, the rephasing and nonrephasing signals are known to consist of contributions from several Liouville pathways which differ with respect to the vibrational state during time t2 (i.e., the ground-state or the first excited state). Thus, the 2DIR spectra reported in this paper represent only one part of the actual 2DIR spectra which arises from Liouville pathways that put the system in the ground vibrational state during t2. Accounting for the contributions from Liouville pathways that involve nonequilibrium excited-state dynamics to the 2DIR spectra in a system of the level of complexity of the AB model represents a rather challenging problem and will therefore be addressed in a separate and forthcoming paper. Acknowledgment. This work was supported by the National Science Foundation, through grants PHY-0114336 and CHE0809506, and by the Petroleum Research Fund. References and Notes (1) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford: New York, 1995. (2) Azzouz, H.; Borgis, D. J. Chem. Phys. 1993, 98, 7361. (3) Schuster, P.; Zundel, G.; Sandorfy, C. The hydrogen bond: Recent deVelopment in theory and experiment (Vols, I-III); North Holland: Amsterdam, 1976. (4) Nibbering, E. T. J.; Elsasser, T. Chem. ReV. 2004, 104, 1887. (5) Giese, K.; Petkovic, M.; Naundorf, H.; Ku¨hn, O. Phys. Rep. 2006, 430, 211. (6) Pimentel, G. C.; McClellan, A. L. The Hydrogen bond; Freeman and Co.: San Francisco, CA, 1960. (7) Jonas, D. M. Annu. ReV. Phys. Chem. 2003, 54, 425. (8) Hamm, P.; Lim, M.; DeGrado, W. F.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. USA 1999, 96, 2036. (9) Hamm, P.; Lim, M.; DeGrado, W. F.; Hochstrasser, R. M. J. Chem. Phys. 2000, 112, 1907. (10) Asplund, M. C.; Zanni, M. T.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. USA 2000, 97, 8219. (11) Zanni, M. T.; Hochstrasser, R. M. Curr. Opin. Struct. Biol. 2001, 11, 516. (12) Zanni, M. T.; Ge, N.-H.; Kim, Y. S.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. USA 2001, 98, 11265. (13) Zanni, M. T.; Asplund, M. C.; Hochstrasser, R. M. J. Chem. Phys. 2001, 114, 4579. (14) Gnanakaran, S.; Hochstrasser, R. M. J. Am. Chem. Soc. 2001, 123, 12886. (15) Ge, N.-H.; Zanni, M. T.; Hochstrasser, R. M. J. Phys. Chem. B 2001, 11, 516. (16) Rubtsov, I. V.; Hochstrasser, R. M. J. Phys. Chem. B 2002, 106, 9165. (17) Rubtsov, I. V.; Wang, J.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. USA 2003, 100, 5601. (18) Khalil, M.; Demirodoven, N.; Tokmakoff, A. J. Phys. Chem. A 2003, 107, 5258. (19) Khalil, M.; Demirodoven, N.; Tokmakoff, A. J. Chem. Phys. 2004, 121, 362. (20) Krummel, A. T.; Mukherjee, P.; Zanni, M. T. J. Phys. Chem. B 2003, 107, 9165. (21) Mukherjee, P.; Krummel, A. T.; Fulmer, E. C.; Kass, I.; Arkin, I. T.; Zanni, M. T. J. Chem. Phys. 2004, 120, 10215. (22) Fulmer, E. C.; Ding, F.; Zanni, M. T. J. Chem. Phys. 2005, 122, 034302. (23) Krummel, A. T.; Zanni, M. T. J. Phys. Chem. B 2006, 110, 913991. (24) Mukherjee, P.; Kass, I.; Arkin, I. T.; Zanni, M. T. J. Phys. Chem. B 2006, 110, 24740. (25) Bredenbeck, J.; Helbing, J.; Behrendt, R.; Renner, C.; Moroder, L.; Wachtveitl, J.; Hamm, P. J. Phys. Chem. B 2003, 107, 8654. (26) Bredenbeck, J.; Hamm, P. J. Chem. Phys. 2003, 119, 1578.

1D and 2D IR Spectra of a Vibrational Mode (27) Bredenbeck, J.; J, H.; Hamm, P. J. Am. Chem. Soc. 2004, 126, 990. (28) Woutersen, S.; Hamm, P. J. Phys. Chem. B 2000, 104, 11316. (29) Woutersendagger, S.; Mu, Y.; Stock, G.; Hamm, P. Proc. Natl. Acad. Sci. USA 2001, 98, 11254. (30) Woutersen, S.; Hamm, P. J. Chem. Phys. 2001, 115, 7737. (31) Woutersen, S.; Hamm, P. J. Chem. Phys. 2001, 114, 2727. (32) Merchant, K. A.; Thompson, D. E.; Fayer, M. D. Phys. ReV. Lett. 2001, 86, 3899. (33) Thompson, D. E.; Merchant, K. A.; Fayer, M. D. Chem. Phys. Lett. 2001, 340, 267. (34) Thompson, D. E.; Merchant, K. A.; Fayer, M. D. J. Chem. Phys. 2001, 115, 317. (35) Hochstrasser, R. M. Proc. Natl. Acad. Sci. USA 2007, 104, 14189. (36) Shim, S.-H.; Strasfeld, D. B.; Ling, Y. L.; Zanni, M. T. Proc. Natl. Acad. Sci. USA 2007, 104, 14197. (37) Kwak, K.; Park, S.; Fayer, M. D. Proc. Natl. Acad. Sci. USA 2007, 104, 14221. (38) Kurochkin, D. V.; Naraharisetty, S. R. G.; Rubtsov, I. V. Proc. Natl. Acad. Sci. USA 2007, 104, 14209. (39) Roberts, S. T.; Loparo, J. J.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 084502. (40) Smith, A. W.; Tokmakoff, A. J. Chem. Phys. 2007, 126, 045109. (41) Nee, M. J.; McCanne, R.; Kubarych, K. J.; Joffre, M. Opt. Lett. 2007, 127, 054504. (42) Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Proc. Natl. Acad. Sci. USA 2007, 104, 14215. (43) Zheng, J.; Kwak, K.; Asbury, J.; Chen, X.; Piletic, I. R.; Fayer, M. D. Science 2005, 309, 1338. (44) Chung, H. S.; Ganim, Z.; Jones, K. C.; Tokmakoff, A. Proc. Natl. Acad. Sci. USA 2007, 104, 14237. (45) Bredenbeck, J.; Helbing, J.; Nienhaus, K.; Ulrich Nienhaus, G.; Hamm, P. Proc. Natl. Acad. Sci. USA 2007, 104, 14243. (46) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194521. (47) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194522. (48) Gundogdu, K.; Bandaria, J.; Nydegger, M.; Rock, W.; Cheatum, C. M. J. Chem. Phys. 2007, 127, 044501. (49) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698. (50) Heyne, K.; Huse, N.; Nibbering, E. T. J.; Elsaesser, T. Chem. Phys. Lett. 2003, 382, 19. (51) Heyne, K.; Huse, N.; Dreyer, J.; Nibbering, E. T. J.; Elsaesser, T.; Mukamel, S. J. Chem. Phys. 2004, 121, 902. (52) 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. (53) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Goun, A.; Fayer, M. D. Phys. ReV. Lett. 2003, 91, 23742. (54) Asbury, J. B.; Steinel, T.; Fayer, M. D. J. Phys. Chem. B 2004, 108, 6544. (55) Steinel, T.; Asbury, J. B.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. Chem. Phys. Lett. 2004, 386, 295. (56) Kolano, C.; Helbing, J.; Kozinski, M.; Sander, W.; Hamm, P. Nature 2006, 444, 469. (57) Woutersen, S.; Mu, Y.; Stock, G.; Hamm, P. Chem. Phys. 2004, 266, 137. (58) 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. (59) Park, S.; Fayer, M. D. Proc. Natl. Acad. Sci. USA 2007, 104, 16731. (60) Cho, M. Chem. ReV. 2008, 108, 1331. (61) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 5827. (62) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 264. (63) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 3840. (64) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2005, 123, 044513. (65) Diraison, M.; Guissani, Y.; Licknam, J. C.; Bratos, S. Chem. Phys. Lett. 1996, 258, 348. (66) Rey, R.; Moller, K. B.; Hynes, J. T. J. Phys. Chem. A 2002, 106, 11993. (67) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 8847. (68) Piryatinski, A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 9664. (69) Piryatinski, A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 9672. (70) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 1623. (71) Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2004, 120, 8107. (72) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2004, 121, 887. (73) Li, S.; Schmidt, J. R.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2006, 124, 204110.

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13003 (74) Scheurer, C.; Piryatinski, A.; Mukamel, S. J. Am. Chem. Soc. 2001, 123, 3114. (75) Scheurer, C.; Mukamel, S. J. Chem. Phys. 2002, 116, 6803. (76) Venkatramani, R.; Mukamel, S. J. Chem. Phys. 2002, 117, 11089. (77) Moran, A. M.; Park, S.-M.; Mukamel, S. J. Chem. Phys. 2003, 118, 9971. (78) Hayashi, T.; Mukamel, S. J. Phys. Chem. A 2003, 107, 9113. (79) Zhuang, W.; Abramavicius, D.; Hayashi, T.; Mukamel, S. J. Phys. Chem. B 2006, 110, 3362. (80) Zhuang, W.; Abramavicius, D.; Voronine, D. V.; Mukamel, S. Proc. Natl. Acad. Sci. USA 2007, 104, 14233. (81) anda, F.; Mukamel, S. J. Chem. Phys. 2006, 124, 124103. (82) Hayashi, T.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 9747. (83) Hayashi, T.; la Cour Jansen, T.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 64. (84) la Cour Jansen, T.; Hayashi, T.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2005, 123, 114504. (85) Hayashi, T.; Mukamel, S. J. Chem. Phys. 2006, 125, 194510. (86) Williams, R. B.; Loring, R. F. J. Chem. Phys. 2000, 113, 10651. (87) Williams, R. B.; Loring, R. F. J. Chem. Phys. 2000, 113, 1932. (88) Williams, R. B.; Loring, R. F. Chem. Phys. 2001, 266, 167. (89) Akiyama, R.; Loring, R. F. J. Chem. Phys. 2002, 116, 4655. (90) Akiyama, R.; Loring, R. F. J. Phys. Chem. A 2003, 107, 8024. (91) Noid, W. G.; Loring, R. F. J. Chem. Phys. 2004, 121, 7057. (92) Noid, W. G.; Ezra, G. S.; Loring, R. F. J. Chem. Phys. 2004, 120, 1491. (93) Noid, W. G.; Loring, R. F. J. Chem. Phys. 2005, 122, 174507. (94) Cho, M. J. Chem. Phys. 2001, 115, 4424. (95) Kwac, K.; Cho, M. J. Chem. Phys. 2003, 119, 2256. (96) Kwac, K.; Lee, H.; Cho, M. J. Chem. Phys. 2004, 120, 1477. (97) Sung, J.; Silbey, R. J. J. Chem. Phys. 2001, 115, 9266. (98) Sung, J.; Silbey, R. J. J. Chem. Phys. 2003, 118, 2443. (99) Kryvohuz, M.; Cao, J. Phys. ReV. Lett. 2005, 95, 180405. (100) Kryvohuz, M.; Cao, J. Phys. ReV. Lett. 2006, 96, 030403. (101) Egorova, D.; Gelin, M. F.; Domcke, W. J. Chem. Phys. 2007, 126, 074314. (102) Steffen, T.; Fourkas, J. T.; Duppen, K. J. Chem. Phys. 1996, 105, 7364. (103) Piryatinski, A.; Tretiak, S.; Chernyak, V.; Mukamel, S. J. Raman Spectrosc. 2000, 31, 125. (104) Mukamel, S.; Abramavicius, D. Chem. ReV. 2004, 104, 2073. (105) Zhang, W. M.; Chernyak, V.; Mukamel, S. J. Chem. Phys. 1999, 110, 5011. (106) Dreyer, J. J. Quantum Chem. 2005, 104, 782. (107) DeVane, R.; Space, B.; Perry, A.; Neipert, C.; Ridley, C.; Keyes, T. J. Chem. Phys. 2004, 121, 3688. (108) Lazonder, K.; Pshenichnikov, M. S.; Wiersma, D. A. Opt. Lett. 2006, 31, 3354. (109) Hamm, P. J. Chem. Phys. 2006, 124, 124506. (110) Ishizaki, A.; Tanimura, Y. J. Chem. Phys. 2006, 125, 084501. (111) Kato, T.; Tanimura, Y. J. Chem. Phys. 2004, 120, 260. (112) Hanna, G.; Kapral, R. J. Chem. Phys. 2005, 122, 244505. (113) Hanna, G.; Geva, E. J. Phys. Chem. B 2008, 112, 4048. (114) Hanna, G.; Kapral, R. J. Chem. Phys. 2008, 128, 164520. (115) Lippincott, E. R.; Schroeder, R. J. Chem. Phys. 1955, 23, 1099. (116) Hammes-Schiffer, S.; Tully, J. C. J. Chem. Phys. 1994, 101, 4657. (117) Staib, A.; Borgis, D.; Hynes, J. T. J. Chem. Phys. 1995, 102, 2487. (118) Antoniou, D.; Schwartz, S. D. J. Chem. Phys. 1999, 110, 7359. (119) Antoniou, D.; Schwartz, S. D. J. Chem. Phys. 1999, 110, 465. (120) McRae, R. P.; Schenter, G. K.; Garrett, B. C.; Svetlicic, Z.; Truhlar, D. G. J. Chem. Phys. 2001, 115, 8460. (121) Kim, S. Y.; Hammes-Schiffer, S. J. Chem. Phys. 2003, 119 (8), 4389. (122) Yamamoto, T.; Miller, W. H. J. Chem. Phys. 2005, 122, 044106. (123) Hanna, G.; Kapral, R. Acc. Chem. Res. 2006, 39, 21. (124) Bruehl, M.; Hynes, J. T. Chem. Phys. 1993, 175, 205. (125) Kim, H. J.; Staib, A.; Hynes, J. T. In Femtochemistry and femtobiology: Ultrafast reaction dynamics at atomic scale resolution; Sundsto¨rm, V., Ed.; Imperial College Press: United Kingdom, 1997; p 510. (126) Hammes-Schiffer, S.; Tully, J. C. J. Phys. Chem. 1995, 99, 5793. (127) Borgis, D.; Tarjus, G.; Azzouz, H. J. Chem. Phys. 1992, 97, 1390. (128) Staib, A.; Borgis, D. Chem. Phys. Lett. 1997, 271, 232. (129) Saven, J. G.; Skinner, J. L. J. Chem. Phys. 1993, 99 (6), 4391. (130) Ka, B. J.; Geva, E. J. Chem. Phys. 2006, 125, 214501. (131) Stephens, M. D.; Saven, J. G.; Skinner, J. L. J. Chem. Phys. 1997, 106, 2129. (132) Golonzka, O.; Tokmakoff, A. J. Chem. Phys. 2001, 115, 297. (133) Khalil, M.; Demirodoven, N.; Tokmakoff, A. Phys. ReV. Lett. 2003, 90, 47401. (134) Everitt, K. F.; Geva, E.; Skinner, J. L. J. Chem. Phys. 2001, 114, 1326.

13004 J. Phys. Chem. B, Vol. 112, No. 41, 2008 (135) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (136) Andersen, H. C. J. Comp. Phys. 1983, 52, 24. (137) Staib, A.; Borgis, D. J. Chem. Phys. 1995, 103, 2642. (138) Yang, C.-Y.; Wong, K. F.; Skaf, M.; Rossky, P. J. J. Chem. Phys. 2001, 114, 3598. (139) Turi, L.; Borgis, D. J. Chem. Phys. 2002, 117, 6186. (140) Nicolas, C.; Boutin, A.; Le´vy, B.; Borgis, D. J. Chem. Phys. 2003, 118, 9689. (141) Borgis, D.; Rossky, P. J.; Turi, L. J. Chem. Phys. 2006, 125, 064501. (142) Evans, D. J.; Holian, B. L. J. Chem. Phys. 1985, 83, 4069.

Hanna and Geva (143) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265– 322. (144) Warshel, A. J. Phys. Chem. 1982, 86, 2218–2224. (145) Kiefer, P. M.; Hynes, J. T. Solid State Ionics 2004, 168, 219– 224. (146) Mukamel, S. J. Chem. Phys. 1982, 77, 173. (147) Shemetulskis, N. E.; Loring, R. F. J. Chem. Phys. 1992, 97, 1217. (148) Bursulaya, B. D.; Kim, H. J. J. Phys. Chem. 1996, 100, 16451. (149) Fried, L. E.; Bernstein, N. B.; Mukamel, S. Phys. ReV. Lett. 1992, 68, 1842. (150) Egorov, S. A.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1998, 108, 1407.

JP804120V