Dissipative wave packet dynamics of the intramolecular hydrogen bond in o-phthalic acid monomethylester

Abstract

We investigate the infrared laser-driven ultrafast dynamics of the carboxy-deuterated title compound in the condensed phase using a system–bath approach. The three-dimensional relevant system comprises the OD stretching and bending motion as well as a low-frequency out-of-plane twisting vibration modulating the hydrogen-bond strength. The dominant relaxation and dephasing channels are identified giving rise to three terms which contribute to the system–bath interaction: (i) A linear system-solvent coupling leading to the relaxation of the low-frequency mode. Here, a classical molecular dynamics simulation is performed to obtain the spectral density and thus the relevant relaxation time scale for this mode. (ii) A coupling which is quadratic in the system coordinates and responsible for pure dephasing. (iii) A fourth-order coupling which involves the release of vibrational energy into the environment by means of the simultaneous excitation of two intramolecular bath vibrations and a solvent mode. Using a quantum master equation approach it is demonstrated that this model is in accord with the results of recent infrared pump-probe and four-wave mixing experiments. The dynamics is discussed in terms of a diabatic representation defined with respect to the high-frequency stretching and bending modes. The behaviour of the low-frequency mode can be characterized as a dissipative wave packet motion in the diabatic ground state which contains contributions from bath-induced coherence transfer out of the initially excited diabatic state.

---

1 Introduction

Recent advances in ultrafast infrared (IR) spectroscopy paved the road to a better understanding of vibrational molecular dynamics in the condensed phase. Special attention has been focussed on H-bonded systems with liquid water being the most prominent example.^1–8 The well-defined intramolecular H-bonds of the type A–H⋯B, on the other hand, are conceptually much simpler than the H-bonded network which characterizes liquid water. Nevertheless the interpretation of IR line shapes is notoriously difficult due to, for instance, the anharmonic coupling between the fast A–H vibration and the slow A⋯B motion which modulates the strength of the H-bond.^9,10 In addition the coupling between A–H stretching and bending vibrations can give rise to typical Fermi-resonance features such as the so-called Evans window, see refs. 11–13. The richness of intramolecular H-bond dynamics which in principle can be burried underneath a broad IR absorption band was demonstrated by the Elsaesser group^14–18 using the femtosecond IR spectroscopy. For carboxy-deuterated o-phthalic acid monomethylester (PMME-D, see Fig. 1)^14,15 which exhibits an intramolecular H-bond of medium strength, the decay of the excited state absorption of the OD stretching vibration in C₂Cl₄ could be characterized by a time scale of 400 fs. A second slower component of 20 ps was attributed to vibrational cooling of the anharmonically coupled modes in the OD vibrational ground state. As a clear indication of the effect of anharmonic coupling between the fast OD stretching and a low-frequency mode modulating the H-bond geometry, an oscillatory component was observed in the signal. It was attributed to damped wave packet motion with respect to the low-frequency mode with the related dephasing time being about 600 fs. For the fundamental OH stretching transition in PMME-H, on the other hand, photon echo studies showed an extremely short dephasing time of about 40 fs.¹⁷

Fig. 1 Most stable structure of PMME-D as obtained from a DFT/B3LYP (6-31+G(d,p)) geometry optimization. The arrows indicate the displacement vector for the low-frequency mode ν_low which corresponds to an out-of-plane torsional motion. The O–O distance of 2.56 Å and the single minimum anharmonic shape of the potential surface indicate a H-bond of moderate strength.

The theoretical investigation of intramolecular IR laser field induced H-bond dynamics requires a multidimensional treatment. For a molecule as large as PMME-D it is not clear a priori what is the relative importance of intramolecular and intermolecular processes for the observed relaxation and dephasing phenomena. Concerning the intramolecular dynamics we have shown that in principle a full-dimensional but approximate quantum simulation is possible if one uses the ab initio reaction surface Hamiltonian approach.^19,20 The Schrödinger equation for this Hamiltonian can then be solved using the time-dependent Hartree method either in the mean-field²¹ or the numerically exact multiconfiguration^22,23 limit. The latter of course requires to select a subset of important modes. For PMME-D this parameter-free approach was capable of predicting the character of the modes which are anharmonically coupled to the OD stretching vibration. Specifically, the mode responsible for the periodic low-frequency modulation of the H-bond dynamics, seen in the experiment, could be identified. However, the estimated decay time of the OD stretching excitation due to intramolecular vibrational energy randomization was about 20 ps, i.e. much longer than the experimental time of 400 fs. (The coincidence with the time scale of vibrational cooling is somewhat accidental as we will see below.) In other words, the interaction with the solvent must play an important role for the relaxation dynamics. In ref. 24 we have developed a system–bath model of PMME-D in the condensed phase. This approach combines an accurate description of the most relevant degrees of freedom, i.e. the OD stretch and bending modes as well as a low-frequency out-of-plane twisting mode, with a perturbational treatment of their interaction with the remaining intramolecular modes and the solvent. The key feature of the model was a fourth-order system–bath coupling which, by virtue of the interplay between two specific intramolecular and a solvent mode, could account for the ultrafast relaxation of the OD stretching vibration. This model has to be distinguished from previous considerations which were aimed mostly at the description of the dephasing of the A–H fundamental transition. Here two models have been discussed which differ in the way the solvent couples to the stretching mode. In the direct coupling case the related dipole moment interacts with the fluctuating field of the polar surroundings.²⁵ Alternatively, an indirect coupling may occur via the damping of the low-frequency mode which has been described either quantum mechanically^26,27 or stochastically.²⁸ On the other hand, Rey and Hynes have considered a linear system–bath coupling model for vibrational relaxation of HOD in D₂O leading, however, to a step-wise mechanism involving the HOD bending fundamental and first overtone transitions.^29,30 Recently, this mechanism was shown to be dominant even if higher-order coupling terms are present.³¹ In the case of PMME-D the fourth-order order coupling opens a direct relaxation channel for the OD stretching vibration. Although of higher order it is a resonant process which can be rate-determining provided that lower order processes are not effective. This has been shown, although in a different context, in refs. 32–34.

Besides the fourth-order coupling our system–bath model includes a linear and a quadratic (in the system coordinate) interaction. From the latter we have taken only the contribution it gives to pure dephasing. While the dynamics of the OD stretching mode and the dephasing of the low-frequency mode could be fitted to the experimental results, the vibrational relaxation of the low-frequency mode could not be modelled without ambiguity. In ref. 24 we have tacitly assumed that the relaxation time for the low-frequency mode will be on the order of the vibrational cooling time of about 20 ps. However, our high-dimensional simulations of the intramolecular dynamics²² suggest that the latter may take place on a similar time scale. In other words, to assume that the 20 ps vibrational cooling time involves the release of energy from the low-frequency mode into the solvent only, is not straightforward. Clearly, to complete the model there is a need for an independent determination of the population relaxation time for the low-frequency mode.

Classical molecular dynamics simulations have been shown to provide an effective way for obtaining relaxation rates for complex solute-solvent situations. The conceptual framework goes back to early work of Landau and Teller and has been reviewed, e.g., in refs. 35 and 36 (see also ref. 37). The link between solute and solvent dynamics is established via time correlation functions of the force which is exerted on the relevant system coordinate by the solvent interactions. In a number of publications such correlation functions have been used to determine population relaxation rates as such^{29,30,38,39,40,41} or for use within density matrix theory.^42,43 The validity of this approach was scrutinized, for instance, in refs. 44–47. In particular, for a bilinear coupling between an otherwise harmonic system and bath, Bader and Berne⁴⁴ showed that classical and quantum relaxation rates coincide while a mixed quantum-classical treatment fails for high-frequency solutes. The conditions for the validity of the linear response theory assumed in this context where investigated by Makri on the basis of an influence functional approach.⁴⁷ These considerations lead to the notion of an effective harmonic bath whose interaction with the relevant system, as manifested in the spectral density function, resembles the real anharmonic solvent environment.

In the present contribution we will extend our previous studies²⁴ and combine the density matrix formulation of the dissipative H-bond dynamics in PMME-D with a molecular dynamics simulation of the solvent spectral density for the low-frequency system mode. This way we will be able to fix the time scale for the relaxation of the latter which will turn out to be much faster than the 20 ps required for the overall vibrational cooling. We start with reviewing our model as well as the link to the classical molecular dynamics simulation of spectral densities in Section 2. The determination of the system–bath interaction parameters is outlined in Section 3.1. In Section 3.2–3.3 the laser-driven dynamics is discussed in terms of a diabatic state representation with respect to the stretching and bending modes. The paper is summarized in Section 4.

2 Theory

2.1 System–bath Hamiltonian

In the following we briefly review the model for the dissipative nuclear dynamics in the adiabatic electronic ground state of PMME-D as put forward in ref. 24. All calculations have been performed using the DFT/B3LYP level of theory with a Gaussian 6-31+G(d,p) basis set.⁴⁸ The most stable geometry is shown in Fig. 1. The potential energy surface for the O–D–O intramolecular H-bond is rather anharmonic, but has a single minimum only, i.e. the H-bond is of moderate strength. Previous investigations on the basis of a high-dimensional reaction surface Hamiltonian have led to the conclusion that the laser-driven dynamics takes place in the vicinity of the most stable configuration.^21–23 This suggests to describe the nuclear motions in terms of normal mode coordinates. Specifically we will include three normal mode vibrational coordinates, {Q_i} = {Q_s,Q_b,Q_low}, into the relevant system, i.e. those of the OD stretching (ν_s) and bending (ν_b) modes and a low-frequency mode (ν_low) which modulates the strength of the H-bond as shown in Fig. 1. The choice of Q_b and Q_low has been motivated by calculations of the forces acting on the normal modes given an elongated Q_s mode.²⁴ Additional support for the model comes from the double peak structure of the IR absorption spectrum in the OD stretch region (cf.Fig. 3 below) as well as from the experimental observation of low-frequency wave packet motion in ref. 15.

The potential energy surface for the relevant system is obtained by performing additional single point calculations according to

V(Q_s, Q_b, Q_low) ≡ V_sys({Q_i}) = V_s(Q_s) + V_sb(Q_s,Q_b) + V_sl(Q_s,Q_low) (1)

Here, the V_s(Q_s), V_sb(Q_s,Q_b), and V_sl(Q_s, Q_low) are obtained with all other coordinates held fixed at their equilibrium positions, e.g., V_s(Q_s) = V(Q_s, Q_b = 0, Q_low = 0). Notice that eqn. (1) is an approximation because instead of a full 3D grid only two 2D grids are spanned. (We used 18 points along Q_s and 7 to 11 points along Q_b and Q_low, depending on the shape of the potential. Subsequently the potential surface was fitted to a high-order polynomial.) We have checked that the direct coupling between modes ν_b and ν_low is negligible for moderate elongations along Q_s.

The interaction of the relevant 3D system with the external electric field E(t) has been included within the dipole approximation via the Hamiltonian

H_f(t) = −E(t)d(Q_s) (2)

The gradient of the dipole moment along the other two modes is more than an order of magnitude smaller than for ν_s, i.e. it can be safely neglected. In other words, in harmonic approximation only the stretching mode is IR active; the actual dipole moment function was fitted to a quadratic function on the Q_s grid.

The 3D system Hamiltonian H_sys = T_sys + V_sys has been diagonalized using an expansion of the total state vector in terms of zero-order states, |φ^(k)_i〉, i.e.

The zero-order states were defined with respect to the reference potentials V_s(Q_s), V_sb(Q_s = 0, Q_b), and V_sl(Q_s = 0, Q_low) using the Fourier grid Hamiltonian method.⁴⁹ To obtain convergence for the lower part of the spectrum up to energies corresponding to the ν_s overtone transition at about 4400 cm⁻¹ it was sufficient to use 5, 8, and 25 basis states for the expansion in the Q_s, Q_b, and Q_low direction, respectively.

The influence of the remaining intramolecular, {q_m}, and solvent, {Z_ξ}, modes on the dynamics of the relevant system is taken into account using the formalism of dissipative quantum dynamics.⁵⁰ Here, the system–bath interaction Hamiltonian is written as

with the system and bath operators K^(u) and Φ^(u), respectively. In principle H_SB can be obtained from a Taylor expansion of the global potential energy surface in terms of the coordinates {Q_i}, {q_m}, and {Z_ξ}. For the solvent such a treatment might imply, e.g., the introduction of instantaneous normal modes⁵¹ or the mapping onto an effective harmonic bath (see below).

A related question concerns the choice of terms to be kept from the Taylor expansion. In ref. 24 we have shown that the experimental observations can be modelled by retaining three terms, i.e.H_SB = H^(I)_SB + H^(II)_SB + H^(III)_SB, which open rather distinct phase and energy relaxation channels.

The bilinear system–solvent coupling

is mainly responsible for the relaxation of the low-frequency mode ν_low. The contribution quadratic in the system coordinates accounts for the ultrafast pure dephasing in all system modes. The rapid relaxation of the OD stretching vibration, on the other hand, can be modelled by a fourth-order coupling which involves both intramolecular and solvent modes:

In eqns. (5)–(7) we introduced the dimensionless coupling strengths g^(…)_i,…(ξ) as well as the frequency of the ξth solvent bath mode Ω_ξ. Further, we use dimensionless coordinates according to Q̃_i = Q_i/λ_i. Below we will take λ_i = (ħ/2ω_i)^1/2 with ω_i being the fundamental transition frequency of the uncoupled mode ν_i. An analogous factor is introduced for the harmonic bath modes.

As pointed out in ref. 24 this form of H_SB represents a minimum model for the envisaged relaxation dynamics. The actual choice has been guided by the “rule” that only lowest-order terms for a particular mechanism should be taken into account (see also refs. 32 and 33).

2.2 Equations of motion

In the following we will focus on the time evolution of the reduced statistical operator ρ = tr_Bρ_tot. Assuming a weak system–bath interaction and the validity of the Markov approximation for the system–bath dynamics the equation of motion for ρ reads⁵⁰where the Liouville superoperators for the coherent evolution have been introduced as L_sys• = [H_sys,•]/ħ and L_f(t)• = [H_f(t),•]ħ. The relaxation and dephasing due to the system–bath interaction is contained in the (Redfield) superoperator R. It can be expressed in the energy representation given by {|a〉} as⁵⁰with the damping matrixbeing defined in terms of the Fourier transform of the bath correlation function aswhere 〈•〉_B denotes the equilibrium expectation value of the bath.

Neglecting the mixing between the three different contribution to H_SB, eqns. (5)–(7) the damping matrix can be split into three terms which have been derived in ref. 24. For later reference we give the contribution of the bilinear coupling which reads:

with

C^(I)_ij(ω) = πω²(1 + n(ω))[J^(I)_ij(ω) − J^(I)_ij(−ω)] (13)

where we defined the spectral densityNote that the second line involves an approximation, i.e. the coupling strength has been expressed in terms of the single parameter g^(I)_i and the new spectral density j^(I)(ω) has been introduced as being equal for all system coordinates.⁵⁰

The quadratic coupling is treated in a similar fashion. However, we will only consider its contribution to pure dephasing, i.e. to those elements of the Redfield tensor which arise from Γ^(II)_ab,cd(ω) with the frequency argument being equal to zero. From the ω → 0 limit of the correlation function one arrives at the following parameterization:^52,53C^(II)_ij(ω → 0) = g^(II)_ig^(II)_jγ_pd, i.e. the magnitude of the pure dephasing term is scaled by the parameter γ_pd which has the dimension of an inverse time.

The correlation function for the fourth-order relaxation contains again a spectral density which is approximated as J^(III)_ij(ω) ≈ g^(III)_ig^(III)_jj^(III)(ω). This implies that g^(III)_i,mn(ξ) ≈ g^(III)_i(ξ), i.e. there is no difference in the coupling of the intramolecular modes. Following ref. 24 we will incorporate two intramolecular modes into eqn. (7) which are selected according to their coupling strength to the OD stretching vibration and the simultaneous requirement that, in combination with a low-frequency bath mode, they should provide a resonant transition for both the ν_s fundamental and the ν_b overtone excitation. The appropriate modes have frequencies of 620 cm⁻¹ and 1329 cm⁻¹.

2.3 Spectral densities

The effect of the solvent on the system dynamics can be entirely characterized by the spectral densities introduced in Section 2.2. This quantity is defined by the density of states of the solvent weighted by the respective coupling strength for the considered mode. Eqn. (14) assumes the validity of a harmonic description of the solvent. Although real solvents are harmonic only under some limiting conditions, e.g., on a short time scale,⁵¹ the actual influence of the solvent can often be mapped onto that of an effective harmonic bath. This can be shown, for instance, by considering a second-order cumulant expansion of the influence functional within the path integral formalism.⁴⁷ This so-called linear response treatment is possible for an arbitrary nonlinear solvent provided that quantum effects can be neglected and that the system–bath coupling is not solvent mode specific, i.e. there is a uniform distribution of coupling strengths over a large number of solvent modes. In the present context these conditions are most likely fulfilled for the low-frequency mode ν_low. The respective correlation function for the low-frequency mode in eqn. (13) then reads⁴⁷

C^(I)_low(ω) = π(1 + n(ω))[J^(eff)_low(ω) − J^(eff)_low(−ω)] (15)

with Here C^(cl)_low(t) is the classical correlation function of the force acting on the system coordinate Q_low.

On the other hand, it is common practice to assume the validity of an effective harmonic bath description and to use a model spectral density within a fully quantum mechanical formulation. For the case of a rather unspecific low-frequency solvent environment the Ohmic spectral density with a cut-off frequency ω^(…)_c is widely used,^50,54i.e.

ω²j^(…)(ω) = Θ(ω)ωexp(−ω/ω^(…)_c) (17)

Note that we included the factor ω² into the definition of the spectral density as it is usually done in literature. Below we will use this form in our simulations. For the case of the bilinear coupling to the low-frequency mode ν_low a comparison with an effective spectral density obtained from classical molecular dynamics simulations will be provided (see Section 3.1).

3 Numerical results

3.1 System–bath interaction

Our relaxation model is characterized by twelve parameters, i.e. the nine coupling strengths g^(…)_i, the scaling factor for pure dephasing γ_pd, and the two cut-off frequencies ω^(I)_c and ω^(III)_c. From the experiment^15,17 it is known that (i) the OD stretching vibration relaxes within 400 fs, (ii) the dephasing time for the OH vibration in the normal species PMME-H is about 40 fs (iii) the observed low-frequency wave packet dynamics is subject to a dephasing on a timescale of 600 fs, and (iv) the overall vibrational cooling takes place in about 20 ps. In addition one has the linear IR absorption spectrum which shows two distinct peaks in the OD stretching region.

In ref. 24 it has been demonstrated that these experimental observations can be reasonably reproduced by the present model. Particularly it should be emphasized that the available parameter space is considerably reduced by the fact that each term in H_SB has a distinct influence on the relaxation and dephasing rates for a given system mode. According to our model the IR spectrum in the OD stretch region is influenced by a Fermi resonance between the ν_s fundamental and the ν_b overtone transition (see Fig. 3 below). Its line width is mostly determined by the pure dephasing and the dephasing coming from the fourth-order coupling acting on the ν_s and ν_b modes. The relative importance of these dephasing channels can be scrutinized by the requirement for a 400 fs population relaxation time for the ν_s excitation. This in turn gives a dephasing time of about 100 fs for the ν_s 0 → 1 coherence.

The situation is a little more complicated for the dynamics of the low-frequency mode. In ref. 24 we have assumed that the equilibration of ν_low proceeds on a time scale longer than 10 ps, i.e. it is related to the vibrational cooling time of 20 ps observed in the experiment. The 600 fs dephasing time for this mode therefore is coming from the pure dephasing term only. However, in the light of our previous investigations²² one might argue that on a picosecond time scale intramolecular vibrational energy redistribution processes take place which are not captured by the present model for the system–bath interaction. In order to obtain a time scale for the population relaxation of the low-frequency mode due to the interaction with the solvent we have performed a classical molecular dynamics calculation of PMME-D in CCl₄ solution at room temperature.†

A single PMME-D molecule frozen at its equilibrium geometry was put into a periodically continued simulation box containing 724 CCl₄ molecules such that the solvent density at room temperature (1.5833 g ml⁻¹) is reproduced. All calculations have been performed using GROMACS 2.0.^55,56 The GROMACS force field is taken for the solute–solvent interactions and the intramolecular solvent dynamics. Electrostatic interactions are not considered. Also the relaxation due to Coriolis and centrifugal coupling is neglected. They have been shown to give a small contribution for HOD in D₂O,²⁹i.e. their effect should be negligible for the much heavier PMME-D. After a 6 ps equilibration of the solvent, a 1 ns production run was performed from which the components of the forces on the solute atoms, F_n(t), are obtained. Given the normal mode transformation matrix U the force on the ith mode can be obtained from where the M_n are the respective masses (see also ref. 29).

From these data the spectral density of the effective harmonic bath can be determined according to eqn. (16). The result is shown in Fig. 2 for the low-frequency region. Here J^(eff)_low(ω) is compared with an Ohmic spectral density (eqn. (17)) having a cut-off frequency of ω_c/2πc = 50 cm⁻¹. The spectral densities are rather similar, with the notable exception that J^(eff)_low has a small shoulder around 200 cm⁻¹ which is due to an intramolecular CCl₄ vibration.‡ As indicated by the arrow in Fig. 2 the relaxation of the ν_low mode does not involve a direct excitation of this solvent mode. Having J^(eff)_low one can calculate state-to-state relaxation rates from k_ab = 2Γ_ab,ba(ω_ab). For the transition between the two lowest states, |a = 2〉 → |b = 1〉, which are almost pure zero-order states of ν_low one obtains 1/k₂₁ = 1.67 ps, a value which is substantially shorter than the observed vibrational cooling time.

Fig. 2 Normalized spectral density for the interaction between the solvent and the low-frequency system mode ν_low; solid-molecular dynamics result (eqn. (16), smoothed for visual clarity), dashed-Ohmic form with ω^(I)_c/2πc = 50 cm⁻¹ (eqn. (17)). The arrow indicates the value of the transition frequency where the spectral density is “probed”.

It has already been emphasized that J^(eff)_low is close to the Ohmic spectral density in the range around 65 cm⁻¹ which is relevant for the ν_low relaxation. Therefore we will stay with the Ohmic form, eqn. (17) for all relaxation processes as in ref. 24. However, the parameters given in ref. 24 will be slightly modified to be conform with the present findings. Specifically, we chose ω^(I)_c/2πc = 50 cm⁻¹ and g^(I)_i = 0.042 which reproduces the |2〉 → |1〉 relaxation rate of the molecular dynamics simulation. In addition, the scaling parameter for the pure dephasing γ_pd is chosen as 0.5 s⁻¹ instead of 0.6 s⁻¹²⁴ in order to have a 500–600 fs dephasing time for the excited low-frequency oscillations. These parameter changes have only a small influence on the relaxation of the Q_s and Q_b modes, which is dominated by the (unchanged) fourth-order coupling.

3.2 Diabatic state dynamics

The time evolution of the relevant system is conveniently followed employing a diabatic state representation of the Hamiltonian H_sys which disentangles the information contained in the eigenstates. A suitable diabatic basis includes the two high-frequency vibrations ν_s and ν_b, i.e.

[T_s + T_b + V(Q_s,Q_b,Q_low = 0)]|ϕ_m〉 = E_m|ϕ_m〉 (18)

Hence the diabatic representation of H_sys readsThe diabatic potential energy curves for the motion of the slow coordinate

U_m(Q_low) = E_m + 〈ϕ_m|V_sys({Q_i}) − V(Q_s,Q_b,Q_low = 0)|ϕ_m〉 (20)

are shown in Fig. 3 for the four lowest states. The state couplings in eqn. (19) are defined as V_mn = 〈ϕ_m|V_sys({Q_i}) − V(Q_s,Q_b,Q_low = 0)|ϕ_n〉. This figure also nicely illustrates the analogy to the treatment of electron-vibrational problems in optical spectroscopy. Stressing this point, the diabatic state populations and the vibrational dynamics along Q_low should be reflected in the IR pump–probe measurements (see, also ref. 15).

Fig. 3 Diabatic potential curves (A) and diabatic state coupling (B) along the low-frequency system coordinate Q_low. In the left part of panel (A) the IR absorption spectrum, I(ω), in the OD stretch region is shown (for a detailed discussion see ref. 24). The parameters of the relaxation model are: g^(I)_i = 0.042, g^(II)_s = 0.005, g^(II)_b,low = 0.0001, γ_pd = 0.5 s⁻¹, g^(III)_i = 0.02, ω^(I)_c/2πc = 50 cm⁻¹, and ω^(III)_c/2πc = 250 cm⁻¹.

In ref. 24 it has been shown that the linear IR absorption spectrum in the OD-stretching region is dominated by a Fermi-resonance with the bending overtone transition. In terms of the present diabatic basis this is reflected in the fact that upon expansion in zero-order states as

one finds that |ϕ₃〉 ∝ −0.32|φ^(s)₁〉|φ^(b)₀〉 + 0.94|φ^(s)₀〉|ε^(b)₂〉 and |ϕ₄〉 ∝ 0.94|φ^(s)₁〉|φ^(b)₀〉 + 0.32|φ^(s)₀〉|φ^(b)₂〉. In addition the diabatic states |ϕ₃〉 and |ϕ₄〉 are coupled with a strength depending on the value of the low-frequency coordinate, see Fig. 3B. In the energy representation of H_sys (cf. Table 2 of ref. 24) this coupling is manifested in the various combination transitions contributing to the IR line shape of Fig. 3A.

The vibrational states of the low-frequency mode pertaining to a particular diabatic state |ϕ_m〉, i.e. the set {|χ_Mm〉} can be obtained from diagonalization of the Hamiltonian T_low + U_m(Q_low) which is conveniently done by expansion in the zero-order basis {|φ^(low)_k〉}, i.e.

The propagation of the reduced density matrix is performed in the energy representation of the total system Hamiltonian, {|a〉}. For an analysis in terms of the more intuitive diabatic states one needs the diabatic representation of ρ, i.e.ρ_mMm,nNn, which is readily obtained from ρ_ab using the transformation

where |m,M_m〉 = |ϕ_m〉|χ_Mm〉 and for the transformation matrix one has

The equations of motion for the reduced density matrix have been solved using the Runge–Kutta–Verner fifth and sixth-order method with adaptive step size control. Following ref. 24 we invoked the Bloch approximation, where the dynamics of populations and coherences is decoupled;⁵⁰ see Section 3.3. All 57 eigenstates up to a transition frequency of 3000 cm⁻¹ have been included in the propagation. In Fig. 4 we have shown the diabatic state populations, P_m(t) = ∑_Mmρ_mMm,mMm(t) after resonant excitation of the ν_s dominated absorption peak at 2288 cm⁻¹ with a laser field having a Gaussian envelope of width 130 fs and being centered at t = 0 fs (see Fig. 3A). For this excitation condition the bending type diabatic state |ϕ₃〉 is only marginally excited via those ν_low states which overlap with the potential curve U₄(Q_low). In the dissipation-free case, Fig. 4(A), there is a coherent flow of population between the stretching fundamental and the bending overtone states |ϕ₄〉 and |ϕ₃〉, respectively. The relatively long time required for maximum population transfer reflects the average state coupling as well as the excitation conditions. There is no noticeable population of the state |ϕ₂〉 (bending fundamental). Switching on the system–bath interaction changes the situation dramatically, see Fig. 4(B). The fourth-order coupling causes a rapid relaxation of the |ϕ₄〉 population back to |ϕ₁〉. Note that the time scale for the relaxation of the diabatic state population P₄(t) is about 340 fs, i.e. a little shorter than the overall relaxation time of the excited stretching and bending type eigenstates which is 400 fs.²⁴ In the diabatic picture the linear system–bath coupling leads to a population transfer between |ϕ₄〉 and the bending overtone state |ϕ₃〉 as seen in Fig. 4(B). The relaxation of the latter back to the ground state is slower by a factor of about two. It should be noticed that also the bending fundamental state |ϕ₂〉 is somewhat populated due to relaxation. It has, however, a small population already in the initial equilibrium state due to a mixing with the overall eigenstates.

Fig. 4 Diabatic state populations after excitation with a laser field having a Gaussian envelope of width 130 fs and being centered at t = 0 fs. The excitation is resonant with the main ν_s dominated absorption peak at 2288 cm⁻¹ and the field strength is taken as 0.001 E_h/(ea₀√amu). Panel (A) and (B) correspond to the coherent and dissipative evolution, respectively, using the parameters given in Fig. 3.

A more detailed look at the dynamics in the diabatic potentials can be obtained from the wave packets with respect to the low-frequency mode. They are obtained from

and are shown in Fig. 5. In panels (A) and (C) the results are given for the state |ϕ₄〉 in the coherent and dissipative case, respectively. The wave packet motion triggered by the laser pulse in the resonantly excited |ϕ₄〉 state is clearly visible in the coherent case (panel(A)). Switching on the dissipation leads to a rapid depopulation of this state and the damped oscillatory motion is hardly discernible. This is in accord with the result of ref. 24 as well as with the experimental findings,¹⁵i.e. the wave packet motion takes place mainly in the ground state |ϕ₁〉, see Fig. 6. The diabatic representation, however, enables one also to investigate the extent of wave packet motion in the bending overtone potential U₃(Q_low). In the coherent case (Fig. 5B) the |ϕ₃〉 wave packet is rather structured, reflecting the overlap integrals 〈χ_M4|V₄₃(Q_low)|χ_N3〉 which enters the diabatic state coupling, but also the resonance condition which facilitates efficient transfer to the states |3,M₃ = 4〉 and |3,M₃ = 5〉. In the dissipative case (Fig. 5D) essentially only these two states are substantially populated. However, this population starts to emerge already during the action of the laser pulse, i.e. in the diabatic state picture it is mostly due to interstate relaxation. Subsequently, the |ϕ₃〉 wave packet is depleted due to the fourth-order relaxation back to the ground state |ϕ₁〉.

Fig. 5 Diabatic state wave packet dynamics along the low-frequency mode according to eqn. (25) ((A) and (B) coherent limit, (C) and (D) dissipative limit). (Note that the grid along the time axis is wider in panel (C) for visual clarity.)

Fig. 6 Time evolution of the coherence density matrix element ρ_1,1,1,2 obtained using the Bloch (solid), secular (long dashed), and full (short dashed) Redfield theory.

3.3 Beyond the Bloch limit

The Redfield equations of motion are commonly used in three limits: full Redfield tensor, secular and Bloch approximation. The Bloch limit discussed so far provides a convenient way for simulating the laser-driven dissipative dynamics of PMME-D. As already indicated the Bloch approximation implies a decoupling between coherence and population dynamics, i.e. ∑_cdR_ab,cdρ_cd = (1 − δ_ab)R_ab,abρ_ab + δ_ab∑_cR_aa,ccρ_cc. It follows as an additional simplification of the Redfield equations in secular approximation.⁵⁰ In the latter it is assumed that the actual relaxation dynamics of the relevant system is much slower than the time scale given by |ω_ab − ω_cd|⁻¹. In secular approximation the elements of R_ab,cd of those pairs of transitions for which ω_ab − ω_cd ≠ 0 holds are neglected. Apparently, this condition is fulfilled for the Bloch limit. However, it can also hold, e.g., for systems having almost harmonic potentials. In such a case bath-induced coherence transfer, ρ_ab → ρ_cd, becomes possible (see, e.g., refs. 57 and 58. Beyond the secular approximation the conversion of populations into coherences, ρ_aa → ρ_cd, and vice versa becomes possible as well (see, e.g., ref. 59). In passing we note that an intuitive phase-space picture of the effect of the secular approximation has been given by Tannor and co-workers in ref. 60.

In order to study the dynamics beyond the Bloch approximation in the present case we have solved the full as well as the secular Redfield equations using the same parameters as before. Because the present system is sufficiently anharmonic such that there is no exact degeneracy between pairs of transitions we have used the criteria that |ω_ab − ω_cd|/2πc < 10 cm⁻¹ for the secular limit. We find that the overall relaxation dynamics of P₄(t) is somewhat accelerated as compared to the Bloch case. There is no noticeable difference between the secular and the full Redfield dynamics for this quantity. A more drastic influence, however, is observed for the dynamics of the coherences and here in particular for the diabatic ground state |ϕ₁〉. This is conveniently illustrated by the matrix elements ρ_1M1,1N1(t) which are responsible for nuclear wave packet motion along Q_low in state |ϕ₁〉. In Fig. 6 we compare ρ_1,1,1,2(t) for the full and approximate Redfield dynamics. Overall the magnitude of the coherence is much larger and the dephasing time is much longer when going beyond the Bloch limit. Apparently, this is due to the bath-induced coherence transfer, e.g., from the pair of diabatic states |4,2〉–|4,1〉 to |1,2〉–|1,1〉. Switching to the eigenstate picture the states having the largest overlap with these diabatic states obey |ω_ab − ω_cd|/2πc = 2 cm⁻¹, which corresponds to a time scale much longer than the relaxation time of |ϕ₄〉. In addition the respective matrix element of the Redfield tensor yields a time constant of 375 fs for the coherence transfer. In other words, this type of coherence transfer is likely to play a role in the present model of PMME-D. Comparing the secular approximation with the full Redfield result we notice that nonsecular contributions have an influence on the coherence dynamics as well. Here the ∼400 fs |ϕ₄〉 → |ϕ₁〉 relaxation time competes with the time scale of ∼500 fs set by the difference |ω_ab − ω_cd|/2πc = 65 cm⁻¹ for transitions, e.g., between the diabatic state pairs |4,1〉–|4,1〉 and |1,2〉–|1,1〉. However, the time scale for this conversion from a population to a coherence in terms of the respective Redfield tensor element is with about 1.4 ps much longer.

Finally, it should be emphasized that the experimental dephasing time for the low-frequency wave packet motion is about 600 fs. The results from this section demonstrate that bath-induced coherence transfer is an important part of the overall dynamics. Of course, in order to establish the relation between the full Redfield dynamics and the experiment would require to modify the parameter set entering our model. Because the population dynamics of the excited diabatic states is only slightly affected by non-Bloch terms, it is primarily the pure dephasing contribution for the low-frequency mode which would become more important in a revised parameter set.

4 Conclusions

The theoretical simulation of condensed phase intramolecular ultrafast H-bond dynamics demands for the combination of methods from quantum chemistry, molecular dynamics, and nonequilibrium statistical mechanics. This high level of description appears to be adequate to match the insight provided by current infrared nonlinear spectroscopies. For the specific example of carboxy-deuterated o-phthalic acid monomethylester we have shown that experimentally observed time scales and phenomena can not only be reproduced, but also related to microscopic processes involving nuclear wave packet motion accompanied by energy and phase relaxation. In the present contribution we have extended our previous model²⁴ to account more accurately for the spectral density describing the interaction between the low-frequency mode modulating the strength of the H-bond and the solvent environment. As a result of a classical molecular dynamics simulation the relaxation time for this mode was determined to be about 1.7 ps. Therefore, it can be concluded that the vibrational cooling which is observed on a 20 ps time scale most likely involves many intramolecular degrees of freedom, in accord with our previous gas phase simulations.²² Thus, the applicability of the present model is limited to ultrafast processes as they occur immediately after laser excitation, i.e. approximately during the first 2 ps.

The dynamics triggered by the laser field is conveniently analyzed in terms of a diabatic representation with respect to the fast stretching and bending modes. This picture facilitates an intuitive illustration of nuclear wave packet motion of the low-frequency coordinate. In the dissipation-free limit, the diabatic state coupling leads to a coherent flow of population between the stretching fundamental and bending overtone state on a time scale of 1–2 ps. Thus, the wave packet motion takes place on both potentials. Switching on the interaction with the bath leads to a rapid (∼400 fs) decay of the excited diabatic state populations and the wave packet motion occurs in the diabatic ground state. Here it dephases with a time constant of 500–600 fs. In this respect we have shown that, in principle, in the course of the excited state relaxation bath-induced coherence transfer and a conversion of populations into coherences can contribute to the ground state wave packet dynamics.

The present choice of expressing the relevant system potential in terms of normal mode coordinates relies on the fact that the H-bond in PMME-D is of moderate strength, i.e. there is no possibility for proton transfer. The theoretical description of proton transfer reactions, of course, requires the use of large amplitude coordinates at least for the moving proton. This can be achieved, for instance, by using a reaction surface description. Recently we have shown that an all-Cartesian formulation provides a convenient way for establishing a system–bath Hamiltonian and for identifying the contribution of specific intramolecular degrees of freedom to dephasing and relaxation rates.⁶¹ Combined with a linear response treatment of the interaction with the solvent this approach should allow for a realistic simulation of the IR spectroscopy of proton transfer systems.

Acknowledgements

This work has been supported by the Deutsche Forschungsgemeinschaft through the Sfb450 and by the Fonds der Chemischen Industrie (O.K.). The authors gratefully acknowledge stimulating discussions with J. Dreyer, T. Elsaesser, and E. Nibbering from the Max Born Institute Berlin.

References

1. H. Graener, G. Seifert and A. Laubereau, Phys. Rev. Lett., 1991, 66, 2092.
2. S. Woutersen, U. Emmerichs and H. J. Bakker, Nature, 1997, 278, 658.
3. S. Bratos and A. Laubereau, in Theoretical Treatments of Hydrogen Bonding, ed. D. Hadži, John Wiley, Chichester, 1997, p. 187.
4. G. M. Gale, G. Gallot, F. Hache, N. Lascoux, S. Bratos and J.-C. Leicknam, Phys. Rev. Lett., 1999, 82, 1068.
5. J. C. Deak, S. T. Rhea, L. K. Iwaki and D. D. Dlott, J. Phys. Chem. A, 2000, 104, 4866.
6. J. Stenger, D. Madsen, P. Hamm, E. T. J. Nibbering and T. Elsaesser, Phys. Rev. Lett., 2001, 87, 2 7401.
7. J. Stenger, D. Madsen, P. Hamm, E. T. J. Nibbering and T. Elsaesser, J. Phys. Chem. A, 2002, 106, 2341.
8. H. J. Bakker and H. K. Nienhuys, Science, 2002, 297, 587.
9. Y. Marechal and A. Witkowski, J. Chem. Phys., 1968, 48, 3697.
10. O. Henri-Rousseau and P. Blaise, Adv. Chem. Phys., 1998, 103, 1.
11. A. Witkowski and M. Wójcik, Chem. Phys., 1973, 1, 9.
12. S. Bratos, J. Chem. Phys., 1975, 63, 3499.
13. O. Henri-Rousseau, P. Blaise and D. Chamma, Adv. Chem. Phys., 2002, 121, 241.
14. J. Stenger, D. Madsen, J. Dreyer, E. T. J. Nibbering, P. Hamm and T. Elsaesser, in Ultrafast Phenomena XII, Springer Series in Chemical Physics, ed. T. Elsaesser, S. Mukamel, M. Murnane and N. Scherer, Springer, New York, 2000, p. 542.
15. J. Stenger, D. Madsen, J. Dreyer, E. T. J. Nibbering, P. Hamm and T. Elsaesser, J. Phys. Chem. A, 2001, 105, 2929.
16. D. Madsen, J. Stenger, J. Dreyer, E. T. J. Nibbering, P. Hamm and T. Elsaesser, Chem. Phys. Lett., 2001, 341, 56.
17. J. Stenger, D. Madsen, J. Dreyer, P. Hamm, E. T. J. Nibbering and T. Elsaesser, Chem. Phys. Lett., 2002, 354, 256.
18. D. Madsen, J. Stenger, J. Dreyer, P. Hamm, E. T. J. Nibbering and T. Elsaesser, Bull. Chem. Soc. Jpn., 2002, 75, 909.
19. B. A. Ruf and W. H. Miller, J. Chem. Soc., Faraday Trans. 2, 1988, 84, 1523.
20. H. Naundorf, J. A. Organero, A. Douhal and O. Kühn, J. Chem. Phys., 1999, 110, 11 286.
21. G. K. Paramonov, H. Naundorf and O. Kühn, Eur. J. Phys. D, 2001, 14, 205.
22. H. Naundorf, G. A. Worth, H.-D. Meyer and O. Kühn, J. Phys. Chem. A, 2002, 106, 719.
23. H. Naundorf and O. Kühn, in Femtochemistry and Femtobiology, eds. A. Douhal and J. Santamaria, World Scientific, Singapore, 2002, p. 438.
24. O. Kühn, J. Phys. Chem. A, 2002, 106, 7671.
25. N. Rösch and M. A. Ratner, J. Chem. Phys., 1974, 61, 3344.
26. B. Boulil, O. Henri-Rousseau and P. Blaise, Chem. Phys., 1988, 126, 263.
27. B. Boulil, J.-L. Dejardin, N. E. Ghandour and O. Henri-Rousseau, J. Mol. Struct. (THEOCHEM), 1994, 314, 83.
28. G. N. Robertson and J. Yarwood, Chem. Phys., 1978, 32, 267.
29. R. Rey and J. T. Hynes, J. Chem. Phys., 1996, 104, 2356.
30. J. T. Hynes and R. Rey, in Ultrafast Infrared and Raman Spectroscopy, ed. M. D. Fayer, Marcel Dekker, New York, 2001, p. 593.
31. C. P. Lawrence and J. L. Skinner, J. Chem. Phys., 2002, 117, 5827.
32. V. M. Kenkre, A. Tokmakoff and M. D. Fayer, J. Chem. Phys., 1994, 101, 10 618.
33. A. Tokmakoff, B. Sauter and M. D. Fayer, J. Chem. Phys., 1994, 100, 9035.
34. S. A. Egorov and J. L. Skinner, J. Chem. Phys., 2000, 112, 275.
35. D. W. Oxtoby, Adv. Chem. Phys., 1981, 47, 487.
36. J. C. Owrutsky, D. Raftery and R. M. Hochstrasser, Annu. Rev. Phys. Chem., 1994, 45, 519.
37. R. Zwanzig, J. Chem. Phys., 1961, 34, 1931.
38. R. M. Withnell, K. R. Wilson and J. T. Hynes, J. Phys. Chem., 1990, 94, 8625.
39. M. Bruehl and J. T. Hynes, Chem. Phys., 1993, 175, 205.
40. S. Gnanakaran and R. M. Hochstrasser, J. Chem. Phys., 1996, 105, 3486.
41. S. A. Egorov and J. L. Skinner, J. Chem. Phys., 1996, 105, 7047.
42. F. E. Figueirido and R. M. Levy, J. Chem. Phys., 1992, 97, 703.
43. H. Gai and G. A. Voth, J. Chem. Phys., 1993, 99, 740.
44. J. S. Bader and B. J. Berne, J. Chem. Phys., 1994, 100, 8359.
45. S. A. Egorov and B. J. Berne, J. Chem. Phys., 1997, 107, 6050.
46. S. A. Egorov, K. F. Everitt and J. L. Skinner, J. Phys. Chem. A, 1999, 103, 9494.
47. N. Makri, J. Phys. Chem. B, 1999, 103, 2823.
48. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. C. V. G. Zakrzewski, J. A. Montgomery, R. E. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P. Y. Ayala, Q. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P. M. W. Gill, B. G. Johnson, W. Chen, M. W. Wong, J. L. Andres, M. Head-Gordon, E. S. Replogle and J. A. Pople, Gaussian 98 (Revision A.7), Gaussian Inc., Pittsburgh, PA, 1998.
49. C. C. Marston and G. G. Balint-Kurti, J. Chem. Phys., 1989, 91, 3571.
50. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, Berlin, 2000.
51. G. Goodyear and R. M. Stratt, J. Chem. Phys., 1996, 105, 10 050.
52. D. Reichman, R. J. Silbey and A. Suarez, J. Chem. Phys., 1996, 105, 10 500.
53. Y. Yan, F. Shuang, R. Xu, J. Cheng, X.-Q. Li, C. Yang and H. Zhang, J. Chem. Phys., 2000, 113, 2068.
54. A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg and M. Zwerger, Rev. Mod. Phys., 1987, 59, 1.
55. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43.
56. D. van der Spoel, A. R. van Buuren, E. Apol, P. J. Meulenhoff, D. P. Tieleman, A. L. T. M. Sijbers, B. Hess, K. A. Feenstra, E. Lindahl, R. van Drunen and H. J. C. Berendsen, Gromacs User Manual version 2.0, Nijenborgh 4, 9747 AG Groningen, The Netherlands, http://md.chem.rug.nl./∼gmx/, 1999.
57. J. M. Jean and G. R. Fleming, J. Chem. Phys., 1995, 103, 2092.
58. D. Egorova, A. Kuhl and W. Domcke, Chem. Phys., 2001, 268, 105.
59. O. Kühn, V. Sundström and T. Pullerits, Chem. Phys., 2002, 275, 15.
60. D. Kohen, C. C. Marston and D. J. Tannor, J. Chem. Phys., 1997, 107, 5236.
61. R. Xu, Y.-J. Yan and O. Kühn, Eur. J. Phys. D, 2002, 19, 293.

---

Footnotes

† Note that the experiment reported in ref. 15 has been performed in C₂Cl₄ solution. Concerning the main results, however, there is no substantial difference compared to CCl₄ (E. T. J. Nibbering, personal communication).

‡ For the isolated CCl₄ one obtains this mode at 225 cm⁻¹ (MP2/6-31+G(d,p)).

---