Two-Dimensional Vibrational Spectroscopy of a Dissipative System with the Optimized Mean-Trajectory Approximation.

The optimized mean-trajectory (OMT) approximation is a semiclassical method for computing vibrational response functions from action-quantized classical trajectories connected by discrete transitions representing radiation-matter interactions. Here we apply this method to an anharmonic chromophore coupled to a harmonic bath. A forward-backward trajectory implementation of the OMT method is described that addresses the numerical challenges of applying the OMT to large systems with disparate frequency scales. The OMT is shown to well reproduce line shapes and waiting time dynamics in the pure dephasing limit of weak coupling to an off-resonant bath. The OMT is also shown to describe a case where energy transfer is the predominant source of line broadening.

Introduction

Multidimensional infrared spectroscopy probes nuclear motions in condensed phase and biomolecular systems through the dynamics associated with vibrational couplings.[1−4] Maximal interpretation of multidimensional vibrational spectra in terms of molecular structure and dynamical fluctuations relies on simulating observables with atomistic models. Fully quantum dynamical calculations of nonlinear response functions are generally not practical for large anharmonic systems, while purely classical response function calculations[5−12] can deviate qualitatively from quantum mechanical results at long times. Several general strategies[13,14] have been successfully used to calculate the response functions that underlie multidimensional vibrational spectra. One of these is to divide a large system into a quantum subsystem coupled to surroundings that are treated classically or semiclassically.[15−27] Another approach is to use quantum calculations on isolated subsystems to develop mappings between vibrational frequencies and collective coordinates such as values of the local electrostatic potential and its gradients. These mappings are then applied to molecular dynamics simulations of the full system to generate a time-varying vibrational frequency, from which the response may be calculated.[20,28−34] A third approach for delocalized chromophores such as amide modes[35] is to apply a tight-binding exciton Hamiltonian with structure dependent couplings to an ensemble of structures generated in a molecular dynamics simulation.[17,36−38] Here we follow a distinct strategy of using a semiclassical approximation[39−49] to quantum vibrational dynamics that in principle can be applied uniformly to all degrees of freedom. This approach has the goal of approximating quantum vibrational response functions from classical trajectories on a specified potential surface.

The optimized mean trajectory (OMT) approximation[50,51] is based on the approximate identification of semiclassical paths with pairs of quantum-mechanical double-sided Feynman diagrams,[1,3,52,53] which depict nonlinear optical processes from a density operator perspective. These semiclassical paths are composed of classical trajectories linked by discontinuities representing the effect of radiation–matter interactions.[54−57] The OMT approach relies on the identification of classical-mechanical action and angle variables;[58] the trajectories are calculated at quantized values of action and the discontinuities are transitions in action at fixed angle values. Such action quantization also features in previous semiclassical approximations to vibrational response functions.[59−61] An advantage of the OMT method is the absence of classical stability matrices that occur in other semiclassical propagation methods,[62−65] while a limitation is the requirement that action and angle variables exist and the necessity of performing canonical transformations between these variables and Cartesian coordinates and momenta.

We have applied the OMT method to compute the third-order vibrational response of thermal distributions of anharmonic oscillators[50] and of coupled pairs[51] of anharmonic oscillators. For these readily assessed cases, the OMT third-order response functions accurately reproduce the time dependences associated with each of the three time variables, describing the dynamics of populations and one- and two-quantum coherences. For our calculations on pairs of strongly coupled, near-resonant anharmonic oscillators, we developed an implementation of the OMT method that employed a small number of precalculated classical trajectories[51] that are efficiently reused to compute dynamics during each of the three time periods. We will refer to this implementation as the “fixed-trajectory” approach. This implementation of the OMT approximation becomes impractical for a large number of degrees of freedom particularly when there is a disparity in frequency scales among different motions.

Here we use the OMT to calculate two-dimensional infrared (2DIR) spectra for a high-frequency anharmonic chromophore coupled to a harmonic dissipative medium.[66−69] Because this model contains both of the features just mentioned, a large number of degrees of freedom and separation of frequency scales between chromophore and medium, the implementation of the OMT used previously[51] is no longer appropriate. We present here a forward–backward implementation[43,49,70] of the OMT, which is better suited to treat systems characterized by these features.

In section II, we review the form of the third-order response function and specify the model. The OMT approximation and its fixed-trajectory implementation are summarized in section III. In section IV, the forward–backward implementation of the OMT, suitable for larger systems with disparities in frequency scales, is presented. Numerical calculations are reported in section V, and are compared to results from the fluctuating frequency approximation[3,53,71−73] and from previous quantum calculations.[69] Conclusions are drawn in section VI.

RESPONSE FUNCTION AND MODEL

The quantum mechanical third order vibrational response function for the signal with wavevector ks = αk1 + βk2 + γk3 iswith α, β, γ, and δ either + or – and the restriction that two signs are + and two −. The rotating-wave approximation has been applied.[53] A single chromophore mode labeled a interacts with the field through its dipole operator, μ̂. We take the dipole operator to be proportional to the coordinate q̂ with the constant of proportionality suppressed. The operators μ̂ ±  can then be expressed in terms of the boson creation and annihilation operators, b̂† and b̂, as μ̂+ = b̂†(ℏ/2mω)1/2 and μ̂– = (μ̂+)†. The equilibrium density operator is given by ρ̂ and (t)Â = e–iÂei, with Ĥ the Hamiltonian without the electric field. The purely absorptive (correlation) spectrum,[3,69]Rabs(ω3,ω1;t2), is obtained by combining the rephasing, RI = R++–, and nonrephasing, RII = R+–+ responseswith R̃(3)(ω3,ω1;t2) the one-sided Fourier transform of R(3)(t3,t2,t1) with respect to t1 and t3.

The chromophore is coupled to a bath of harmonic oscillators[66−68]Here Ĥ is the chromophore Hamiltonian, Nb is the number of bath modes, m and ω are the mass and frequency of bath mode j, and c quantifies the coupling between mode j and the chromophore. The dependence of the coupling on the chromophore coordinate is determined by V(q̂), which has dimensions of length. We will consider this function up to second order in chromophore coordinate[69]with vLL and vSL quantifying linear–linear (LL) and square-linear (SL) coupling, respectively. This Hamiltonian with LL coupling has been used to quantify the effects of a dissipative bath on quantum corrections to the classical vibrational linear response function.[61]

The bath coupling coefficients and frequencies can be obtained from a specified spectral density, J(ω) ≡ ∑c2/(2mω)δ(ω – ω). We take J(ω) in the continuum limit to be Ohmic with a Lorentzian cutoff[69]where γ controls the width of J(ω) and η is the classical friction coefficient for vLL = 1. This continuous distribution can be approximated by a discrete bath of oscillators with linearly spaced frequencies, ω = jΩ/Nb, and maximum frequency Ω by taking the coupling constants c to be[74]

To better approximate a continuous bath, an additional zero frequency mode was included in eq 3.[74] This mode contributes a frequency shift to the chromophore. The coupling coefficient of this zero frequency mode and of mode Nb with frequency Ω is given by half of the right side of eq 6. We describe the chromophore–bath coupling in terms of the dimensionless parameters νLL = vLL(η/mω)1/2 and νSL = vSL(ηℏ/m2ω2)1/2. Quantum calculations of purely absorptive spectra for the model of eqs 3–5 with a continuous bath were performed by Ishizaki and Tanimura[69] using a quantum Fokker–Planck equation approach.[75,76] Our semiclassical results are compared to these quantum calculations in section V.

OPTIMIZED MEAN-TRAJECTORY APPROXIMATION

The OMT approximation for the nonlinear vibrational response of a single oscillator[50] and its extension to a collection of coupled oscillators[51] have been previously described. Here the OMT will be reviewed, focusing on computing the rephasing and nonrephasing responses for a chromophore weakly coupled to a bath. The OMT results from translating pairs of quantum mechanical double-sided Feynman diagrams into semiclassical OMT paths. Double-sided Feynman diagrams portray the perturbative time evolution of the density operator throughout a spectroscopic experiment.[1,3,52,53,73] The identification of double-sided Feynman diagrams with OMT paths relies on evaluating the quantum diagrams in a harmonic approximation to the energy eigenstate basis. It further rests on the analogy between energy eigenstates in quantum mechanics and the invariant tori defined by the classical-mechanical action–angle variables,[58] in cases where the latter can be defined. In analogy to free propagation in the energy basis interrupted by the radiation-induced transitions between eigenstates in the quantum diagrams, OMT paths consist of classical trajectories propagated at quantized action values connected by jumps in action at constant angle representing radiation–matter interactions. For a single degree of freedom, the evolution of the density operator |n⟩⟨n| in a double-sided Feynman diagram is represented in an OMT path as a classical trajectory with action ℏ(n + n + 1)/2, the mean action associated with the harmonic energies of the bra and ket aspects of the density operator. Interactions with the electric field are similarly treated within this harmonic approximation as jumps in action of magnitude ℏ/2 at constant angle. Although these quantization rules are derived within a harmonic approximation, they are applied to the good action variables of the system, thereby including anharmonic effects.

Within the quasiperiodic regime[77,78] where action–angle variables can be defined, the OMT for one oscillator can be extended to multiple degrees of freedom.[51] Anharmonic terms in the Hamiltonian are initially ignored, and the harmonic approximations of the OMT for a single degree of freedom are applied to each normal mode. The quantization rules are then applied to the action variables of the full coupled Hamiltonian. This approach relies on each action variable being unambiguously associated with one normal mode. This will hold if the anharmonic couplings between normal modes in the full Hamiltonian are sufficiently small that the transformation to good action–angle variables is well approximated by perturbation theory in anharmonicity.[51,77] For both single and multimode models, OMT paths can be attributed to specific physical processes using the connection between these semiclassical diagrams and double-sided Feynman diagrams.

For f coupled oscillators the OMT approximation for the quantum mechanical response in eq 1 is[51]Here κ is the overlap of normal mode r with the local chromophore mode a, given by the expansion q = ∑ κx , with {x} the normal mode coordinates. In general, contributions to the system response result from either one normal mode interacting with the electric field, as in the first sum in eq 7 or two normal modes, as in the second sum. The relative contribution of each term in eq 7 is determined by the factors κ4 or κ2κ2. For the couplings investigated here, a single normal mode will have significant overlap with the chromophore. We will refer to this mode as the “system” normal mode labeled 1 and the others as “bath” normal modes. Because κ2 ≫ κ2 for r ≠ 1 the OMT calculation can be greatly simplifiedThe superscript has been omitted from ρ because there is only one contributing term. Within this approximation only action jumps in the system normal mode are allowed. This corresponds to computing contributions to the two-dimensional spectrum from the diagonal peak near the fundamental frequency of the system normal mode and the peaks associated with overtone transitions in this mode.

For the response function components in eq 2, ρ is given by[51]

The index p in eq 9 sums over semiclassical paths, composed of classical trajectories linked by jumps in the system action at fixed angle values. There are four such paths, corresponding to the combinations of either increasing or decreasing the system’s action at the second and third interactions with the field. Whether the action is increased or decreased following a particular trajectory in the OMT path is controlled by the subscripts σ, which can be either + or −, in the terms Δσ( defined in eq 13. Superscripts s = 1, 2, 3 label the three classical trajectories of the OMT path. In eqs 9–14, z((t) refers to the phase space variables z for trajectory s of the OMT path at time t; if no time is specified, t = 0. Subscripts on action values, J, and coordinates and momenta, x and p, in eqs 10–14 indicate normal modes. As defined in eq 12, F is the classical distribution function renormalized to reflect the quantization rules imposed on the action values at equilibrium.

Although the expression for ρ+±∓ is written as three 2f-dimensional integrals, delta functions in angle in eq 9 and in action in eq 13 simplify the integrations over z(2) and z(3), so that only the integration over z(1) needs to be performed numerically. The integration over action variables, J(1), is further reduced to f sums by Γ(z(1)) in eq 10. Therefore, it is convenient to visualize the OMT calculation as a summation of OMT paths with different initial conditions z(1). For the rephasing and nonrephasing responses the four possible paths are shown in Figure 1.[50,51]

Figure 1

The four OMT paths contributing to the rephasing and nonrephasing signals are shown for system normal mode 1 having action J1(1) following the first interaction with the field. Time increases from left to right. Interactions with the field are treated as ℏ/2 jumps in the system action with all angle values held constant. Colored dots represent states on the OMT path used in the calculation of Q σ factors in eq 9. Green squares represent a set of z(3)(t3) values used in the calculation of Q–(z(3)(t3)).

Figure 1 shows the action of the system normal mode as a function of time; the action values of the bath normal modes are not shown and remain constant throughout the path. Dotted horizontal lines indicate half-odd-integer multiples of ℏ, corresponding to the system density matrix being diagonal in the eigenstate basis. Integer multiples of ℏ correspond to single-quantum coherences and are indicated by dashed horizontal lines. OMT paths begin after the first radiation–matter interaction. The first sum in eq 10 restricts the action of the system normal mode in trajectory 1, J1(1), to be an integer multiple of ℏ. The dots at the beginning and end of this segment represent points in the path where factors of Q σ in eq 9 are evaluated. The blue dot corresponds to z(1), and the red dot to z(1)(t1). At time t1 the system mode interacts with the field, either increasing or decreasing its action by ℏ/2. This transition is controlled by the factor Δσ(1) in eq 9. All angle values are unchanged during this jump in action as required by δ(ϕ(2) – ϕ(1)(t1)). The second trajectory is propagated for time t2 and another ℏ/2 jump in action is performed at constant angle according to the factors Δσ(2) and δ(ϕ(3) – ϕ(2)(t2)). The purple dot at the start of trajectory 3 represents z(3). The green squares represent z(3)(t3) values used to compute the response as a function of t3 for a particular t1 and t2 value.

The statistical weight of each path is determined by the renormalized classical equilibrium distribution function, F, in eq 12, evaluated at the two allowed states prior to the first interaction with the field. The actions of these states are obtained by replacing the system action, J1(1), in J(1) with J1(1) – ℏ/2 or J1(1) + ℏ/2 and the angle values are the same as after the first field interaction. This total weight is given by eq 11, where the terms have opposite sign due to the innermost commutator in eq 1. In eq 9 ϵ is an overall sign arising from the remaining commutators. For the two paths with action J1(1) during the third trajectory ϵ = −1, while ϵ = 1 for the paths finishing with system action J1(1) + ℏ or J1(1) – ℏ.

Three numerical challenges must be addressed in applying the OMT method. First, trajectories need to be propagated with their initial phase space points z( specified in action–angle variables. Second, jumps in action at constant angle must be performed. Third, a 2f-dimensional phase space integration over z(1) must be performed.

For two coupled near-resonant oscillators, we developed a fixed-trajectory implementation that addressed each of these challenges.[51] The minimum number of trajectories is required in this implementation because only one classical trajectory at each set of actions reached during any OMT path is computed. These trajectories are precalculated at the start of the computation and are reused to determine time-evolution for each appropriate segment of an OMT path. Initial conditions for these trajectories are determined using a perturbative approximation to the canonical transformation between Cartesian coordinates and momenta and action–angle variables.[77] Computing this transformation to second order in cubic anharmonicity and first order in quartic anharmonicity was sufficient to reproduce the quantum mechanical response for a chromophore with relatively high anharmonicity.[51] To carry out action jumps, approximate constant-angle mappings between the fixed classical trajectories are also calculated at the start of the computation. To determine a constant-angle jump between an initial and final trajectory, a target state is computed by harmonically scaling the initial coordinate and momentum of the normal mode interacting with the field by the square root of the ratio of that mode’s final and initial actions. With the unscaled phase space points of the other normal modes, this defines a 2f-dimensional target phase space point that would represent a jump in action at constant-angle variables in the absence of anharmonicity. The closest point in phase space on the final trajectory to this target state is used as the end point of the approximate constant-angle jump. For relatively high-frequency oscillators, few quantum states are thermally accessible so that the sum over initial actions requires a small number of terms. The angle average is performed in the fixed-trajectory implementation as a time average by varying the starting point in trajectory 1.

Although this method is well suited to treat a few high-frequency oscillators,[51] there are difficulties when it is applied to larger systems with disparate frequency scales. The simplicity of the initial action sum is lost when the OMT is applied to low frequency oscillators because they thermally access a large range of action values, making it impractical to compute the response function from one fixed set of trajectories. A second challenge caused by a low frequency bath is that very long trajectories must be computed to sample all combinations of bath normal mode angles. The constant-angle jumps are also more challenging for larger systems because a 2f-dimensional minimization must be performed to determine constant-angle mappings between the trajectories. These challenges motivated the development of an alternative OMT implementation described in the following section.

FORWARD–BACKWARD OMT IMPLEMENTATION

A forward–backward implementation of the OMT was developed to address the challenges of applying the OMT to large systems including low-frequency oscillations. This implementation has numerical advantages similar to the doorway- and window-function factorization[79,80] described by Hasegawa and Tanimura for the computation of nonlinear spectra with nonequilibrium molecular dynamics simulations.[70] Our method is illustrated in Figure 2 for the OMT path in Figure 1 in which both interactions with the field increase the system action. In this implementation we begin OMT paths at z(1)(t1), represented by the red dot, instead of z(1). Replacing z(1) by z(1)(t1) in the factors Γ and ΔF leaves eq 9 unchanged because both points lie on the same constant-action trajectory and in the good action-angle variables of the coupled Hamiltonian, F has no angle dependence. The contribution from this path to ρ+±∓ in eq 9 can be written in the formIn eq 15, the integrand is factored into two expressions, grouped in square brackets. The first factor is independent of t2 and t3, and the second is independent of t1. Therefore, the integrand can be computed as the outer product of these two terms.

Figure 2

Forward–backward implementation of an OMT path for a fixed value of t2 is shown. The path is initialized at z(1)(t1), represented by a red dot, with known ϕ(1)(t1) and J(1). The t1 trajectory is propagated backward from this point to obtain a set of z(1) values indicated by blue squares. A trajectory for t2 is propagated forward in time from the phase space point z(2), with known initial angles ϕ(2) = ϕ(1)(t1). For each t2 value the point after an ℏ/2 jump in action z(3), shown as a purple dot, must be determined. From this initial condition a t3 trajectory is propagated, giving a set of z(3)(t3) values, shown as green squares.

In this implementation, the phase space points z(1) are obtained by propagating the initial point z(1)(t1) backward. Values of z(1) are represented by blue squares in Figure 2. When values of z(1)(t1) are sampled directly, ϕ(1)(t1) is known. Therefore, computing z(2) requires no approximations in addition to those used to compute z(1)(t1). This is indicated by the arrowless dashed line connecting trajectories 1 and 2 in this path. The phase space point z(2) is then propagated forward in time to z(2)(t2). For each t2 value a jump in system action at constant-angle values is approximately performed to obtain initial conditions z(3), which are propagated forward in time to obtain a set of z(3)(t3) values.

To compute the contribution to the system response from each OMT path in the forward–backward implementation a single t1 and a single t2 trajectory are propagated. Then, for each t2 value, a distinct t3 trajectory is propagated. Therefore, the total number of trajectories scales as the number of t2 values, n. The number of computations required to perform constant-angle jumps also scales with this number. By comparison, for the implementation suggested by Figure 1, where z(1) serves as the initial point instead of z(1)(t1), a single t1 trajectory is propagated for each path. For each t1 value a distinct t2 trajectory is propagated and for each combination of t1 and t2 values a distinct t3 trajectory is propagated. In this implementation, the total number of propagated trajectories scales as n(1 + n), as does the number of calculations needed to perform constant-angle jumps. Typically, 2DIR spectra are plotted as Fourier transforms with respect to t1 and t3 for a small number of t2 times, so that n ≪ n. The forward–backward implementation is therefore significantly more efficient than the implementation of the OMT suggested by Figure 1, although it does require more trajectories than the fixed-trajectory implementation described in section III.

For the fixed-trajectory implementation, initial Cartesian coordinates and momenta were determined in terms of action and angle variables using low-order perturbation theory applied to the normal modes.[51,77] For the systems investigated here, this canonical transformation was approximated to zeroth order in the anharmonic couplings involving bath normal modes and to first order in cubic anharmonicity for the system normal mode. These results were used to transform the action–angle variables into Cartesian coordinates and momenta, and the full anharmonic Hamiltonian was used to propagate all trajectories.

The forward–backward implementation avoids difficulties associated with applying the fixed-trajectory implementation to systems with low-frequency oscillations through its treatment of action jumps. Because of the relatively small number of jumps in the forward–backward implementation, these transitions can be treated within the same perturbative framework used to compute initial phase space conditions. For the results shown in section V the perturbative coordinate and momentum expressions were numerically inverted to determine ϕ(2)(t2) = ϕ(3), so that initial conditions z(3) could be computed.

The final challenge to be addressed is performing the phase space integration. Low-frequency bath oscillators thermally sample a large number of action values so that it is impractical to directly sum all combinations of these actions. Therefore, initial conditions for z(1)(t1) were sampled from the equilibrium distribution function F in eq 12 using a Metropolis criterion. The difference in eq 11 was approximated as ΔF(z(1)(t1)) ≈ F(z(1)(t1))|, which is valid for βℏω1 ≳ 1. For sampling, the dependence of the full Hamiltonian on action variables was approximated as a sum of uncoupled contributions with the system mode treated perturbatively to second order in cubic anharmonicity and the bath modes treated harmonically. The forward–backward implementation was used to compute all OMT results in section V.

Results

We treat the Hamiltonian in eq 3, with the chromophore mode taken to be a Morse oscillator defined by dimensionless parameters, βD = 391 and βℏω = 7.75. Here β ≡ 1/kBT, D is the Morse oscillator well depth, and ω is its harmonic frequency. At 300 K for the fundamental frequency ω10 = ω – Δanh = 1600 cm–1 this corresponds to anharmonicity Δanh = ℏω2/2D = 16 cm–1. These parameters are appropriate to describe the amide I band[14,35] and are used in ref (69). This set of chromophore parameters is used in all figures.

We compare our results to the quantum calculations of Ishizaki and Tanimura[69] and to the widely applied fluctuating frequency approximation,[71,72] where the effect of the bath interaction is approximated by Gaussian fluctuations in the chromophore frequency. This approximation is expected to be valid for weak coupling to an off-resonant bath where dissipation is not significant.[69] The rephasing and nonrephasing system responses for this approximation[3,53,73] are given in eqs 2.11a and 2.11b of ref (69). The response functions are determined by the frequency autocorrelation function C(t) ≡ ⟨δω10(t) δω10(0)⟩, which in the present model is given in terms of the classical friction kernel η(t) for vLL = 1 asThis is the discrete bath analog of the results in eqs 2.15 and 3.18b of ref (69) for a continuum bath, with coordinate matrix elements evaluated to lowest order in anharmonicity.

In Figure 3 purely absorptive spectra are shown for bath parameters in the pure dephasing regime in which the fluctuating frequency approximation is expected to accurately reproduce the response.[69] The width of the spectral density in eq 5 is taken to be γ = 4.95 × 10–3ω , with Nb = 20 and maximum bath frequency Ω = 0.02ω. The chromophore–bath coupling is bilinear with νLL = 1.41 and νSL = 0. These parameters correspond to a finite bath version of the calculation in Figure 5(a-i) of ref (69). Purely absorptive spectra for the fluctuating frequency approximation are shown in column (a) and for the OMT approximation using 5000 initial conditions in column (b). Rows (i)–(iii) show spectra for ωt2 = 0, 150, and 1200, respectively. For typical amide I vibrational frequencies[14,35] the waiting times in rows (ii) and (iii) correspond to approximately 0.5 and 4 ps.

Figure 3

Rabs(ω3,ω1;t2) is calculated for an anharmonic oscillator interacting bilinearly with a medium. Column (a) shows results of the fluctuating frequency approximation, and OMT results are shown in (b). Spectra are shown for ωt2 = 0 in row (i), ωt2 = 150 in row (ii), and for ωt2 = 1200 in row (iii). Results for each approximation are normalized to the maximum absolute value at t2 = 0. Six contours equally spaced between −1 and 0 and between 0 and +1 are shown, with negative contours in blue and positive in red.

The fluctuating frequency approximation spectrum at t2 = 0 in Figure 3(ai) shows diagonal elongation, indicating inhomogeneous broadening that is characteristic of waiting times short relative to relaxation time scales. The OMT spectrum at t2 = 0 in (bi) agrees qualitatively with both the fluctuating frequency result in (ai) and the corresponding quantum calculation in ref (69). As the waiting time is increased in rows (ii) and (iii), changes to the relative contributions of the rephasing and nonrephasing responses cause the peaks to become more symmetric, indicating homogeneous broadening. Both the degree of broadening and the crossover between diagonally broadened and symmetric peaks are accurately reproduced by the OMT results, as seen from comparison of the two columns in Figure 3. The t2 dependence was not considered in ref (69), so no comparison to quantum mechanical results is made. In the pure dephasing limit, all broadening from LL coupling is a result of the system anharmonicity, as is evident in eq 16 where the LL term is proportional to D–1/2. Accurately reproducing the broadening present in the fluctuating frequency approximation demonstrates that important anharmonic effects are incorporated into the OMT description by propagating trajectories with the full Hamiltonian. Overall, the fluctuating frequency and OMT approximation results are similar for this set of parameters, except for a bath-induced shift in the centers of both peaks in the OMT approximation which is also present in the t2 = 0 quantum results of ref (69).

In Figure 4, the calculations in Figure 3 are repeated for chromophore–bath coupling that is quadratic in the chromophore coordinate, νLL = 0, νSL = 0.704. These correspond to the coupling parameters used in Figure 5(a-ii) of ref (69). The OMT results in column (b) of Figure 4 were computed using 35 000 initial conditions, although qualitative features were apparent with a few thousand initial conditions. The waiting time dynamics of the spectra in Figure 4 are qualitatively similar to the dynamics in the LL coupling case. At t2 = 0 peaks show inhomogeneous broadening and become homogeneously broadened as the waiting time increases. There is greater line broadening for the SL coupling in Figure 4 than for the LL coupling shown in Figure 3, in agreement with the t2 = 0 calculations in ref (69). The OMT approximation reproduces the line shapes of the fluctuating frequency approximation, including the decay in the peak amplitude as a function of t2 as well as the degree of dephasing relative to the LL case. Unlike LL coupling terms, anharmonic SL coupling terms do not enter in determining normal modes. The SL coupling terms are only incorporated in the OMT approximation through their presence in the full Hamiltonian used to propagate trajectories. The results in Figure 4 again demonstrate the capacity of the OMT to reproduce anharmonic effects, even when the action–angle variables are crudely approximated.

Figure 4

Purely absorptive spectra are calculated for the parameters in Figure 3 but with quadratic coupling, νLL = 0, νSL = 0.704. Fluctuating frequency approximation results are shown in column (a), and OMT approximation results in column (b). The t2 values are the same as in Figure 3.

We have further investigated the waiting time dynamics predicted by the OMT approximation. Figure 5 shows the waiting time dynamics of the absolute value of the rephasing response in (a) and of the nonrephasing response in (b) for the same parameters as in Figure 4, with all signals normalized to their t2 = 0 value. Fluctuating frequency results are shown as dashed lines and OMT results are shown as solid lines. Three values of t = t1 = t3 are shown for each signal, ωt = 30 (blue), 90 (red), and 180 (purple). For a fundamental transition frequency of 1600 cm–1 these times correspond to approximately 100, 300, and 600 fs, respectively.

Figure 5

|RI/II(3)(t,t2,t)/RI/II(3)(t,0,t)| is calculated for the same parameters as in Figure 4. The rephasing response is shown in panel (a), and the nonrephasing response in panel (b). Three t = t1 = t3 values are shown in each panel: ωt = 30 (blue), 90 (red), and 180 (purple). Fluctuating frequency approximation results are shown with dashed lines, and OMT results with solid lines.

The fluctuating frequency approximation rephasing responses in Figure 5(a) show an overall decay with waiting time, with greater relative decay for larger ωt values. The rephasing and nonrephasing signals in (a) and (b) are reciprocals in the fluctuating frequency approximation, so that the nonrephasing results show corresponding increases with waiting time. The OMT results share these features. Relative to the fluctuating frequency approximation, the OMT results consistently show less decay in the rephasing signal and smaller increases in the nonrephasing signal, so that the OMT results in (a) and (b) are also approximately reciprocals. In (b) qualitatively similar small amplitude oscillations, caused by the finite bath, are apparent in both results, especially for ωt = 180 shown in purple. This is a rigorous test of the OMT because all t2 dynamics of the response functions in Figure 5 are the result of chromophore–bath couplings. Results in ref (69) were computed at t2 = 0, so no comparison to quantum calculations is made.

Figures 3–5 show results in the pure dephasing regime where the fluctuating frequency approximation is expected to accurately reproduce the response function.[69] This approximation will not describe a case in which energy transfer between the system and bath is significant. Figure 6 shows 2DIR spectra for such a model, where the width of the bath spectral density has been increased relative to that of Figures 3–5 and the chromophore–bath coupling is bilinear, facilitating single quantum excitation transfer between the chromophore and bath. The bath parameters are γ = 9.90 × 10–2ωa, Nb = 125, and Ω = 1.4654ω , with the maximum bath frequency chosen to avoid resonances with the chromophore mode. The coupling parameters are νLL = 0.222 and νSL = 0. These are the same coupling strengths used in row (i) of Figure 3 in ref (69) but here the width of the spectral distribution is a factor of 5 smaller to reduce the number of oscillators in the finite bath. Fluctuating frequency results are shown in column (a) of Figure 6, and OMT results computed from 5000 initial conditions are shown in column (b). Rows show results for the same waiting times as in the corresponding rows of Figures 3 and 4. Time-domain results were not fully decayed so, to reduce artifacts[12] caused by taking the discrete Fourier transform of aperiodic data, the response functions used to compute Figure 6 were multiplied by the product of one-sided cosine-squared window functions for the t1 and t3 time variables. Applying this window function to the time-domain results used to compute Figures 3 and 4 did not result in significant additional broadening.

Figure 6

Rabs(ω3,ω1;t2) for a thermal ensemble of Morse oscillators with bilinear coupling to a harmonic bath is shown as a function of ωt2. Fluctuating frequency approximation results are shown in panel (a), and OMT results are shown in panel (b). Spectra at ωt2 = 0 are shown in row (i), at ωt2 = 150 in row (ii), and at ωt2 = 1200 in row (iii). All spectra are normalized to the maximum absolute value at t2 = 0. Six contours equally spaced between −1 and 0 and between 0 and +1 are shown, with negative contours in blue and positive in red.

In column (a) of Figure 6 the fluctuating frequency approximation results show no significant broadening or waiting time dynamics at the figure resolution. This indicates minimal pure dephasing for this set of parameters. In contrast, the OMT results in column (b) show significant line broadening. At t2 = 0 the OMT result in Figure 6(bi) shows inhomogeneous broadening, while both peaks at finite t2 are homogeneously broadened. Fluctuating frequency and OMT spectra calculated with the same parameters as in Figure 6 but with quadratic coupling νLL = 0 and νSL = 0.222, as in Figure 3(ii) in ref (69), show minimal broadening for all waiting times. Quadratic coupling facilitates exchange of two system quanta and one bath quantum and so is unlikely to produce line broadening for this spectral density. We do not make a direct comparison to Figure 3 of ref (69) because the results here are for a smaller γ and, with relatively high Ω/Nb, the finite bath does not well represent the continuum. However, the spectra in the energy transfer regime shown in Figure 3 of ref (69) display some of the same features as the results here. First, the fluctuating frequency result for LL coupling has no significant broadening. Second, the quantum LL result shows no apparent inhomogeneous broadening at t2 = 0, while our results show diminished inhomogeneous broadening at finite t2 compared to the pure dephasing results in Figures 3 and 4. Finally, the minimal broadening with SL coupling predicted by the OMT is consistent with Figure 3(a-ii) in ref (69). This qualitative agreement indicates that the OMT treatment is applicable in regimes in which energy transfer processes are a significant source of line broadening.

Conclusions

The OMT approximation to nonlinear vibrational response functions results from the assignment of semiclassical paths[50,51] to sums of pairs of double-sided Feynman diagrams,[3,52,53,73] which represent additive contributions to exact quantum response functions. The OMT procedure allows nonlinear vibrational response functions to be calculated from classical trajectories linked by transitions representing radiation–matter interactions. Implementing the method requires approximating the canonical transformation between Cartesian coordinates and momenta and good action and angle variables, with the assumption that the latter exist. This transformation is necessary to select initial conditions that correspond to quantized action values and also to execute transitions between trajectories that correspond to discrete changes in action at constant angle. We have shown here and previously[50,51] that this canonical transformation can be performed with relatively crude approximations such as low-order perturbation theory in anharmonicity. The calculated response functions, however, are nonperturbative in anharmonicity, because all trajectories are propagated numerically using the full classical Hamiltonian. For the model studied here, an anharmonic oscillator coupled to a harmonic bath, the necessary canonical transformations were carried out as follows. A normal mode transformation was applied to the full anharmonic Hamiltonian. The canonical transformation to action and angle variables was then approximated to zeroth order in all anharmonic couplings involving bath normal modes and to first order in cubic anharmonicity in the system normal mode. Yet because correct numerical trajectories were used, the calculations in section V correctly reproduce anharmonic effects of chromophore-bath coupling in the regimes of pure dephasing (Figures 3–5) and energy transfer (Figure 6).

An advantage of the OMT approach is that different peaks in a 2DIR spectrum are associated with different semiclassical paths, just as they are associated with different double-sided Feynman diagrams. Using this association, contributions from different physical processes can be identified and calculated separately. For example, if the system mode is restricted to the ground state at thermal equilibrium, the overtone peaks in section V result solely from the semiclassical diagram where the system action increases with the second and third field interactions, as in Figure 2. The diagonal peaks would then result from the sum of the other two allowed semiclassical diagrams in Figure 1, where the system action is the same in the first and third trajectories. This association of semiclassical paths with spectral peaks emphasizes the physical significance that may be assigned to these paths.

A naive implementation of the OMT as suggested in Figure 1 is numerically challenging, with the required number of trajectory propagations and action transitions for each OMT path scaling as n(1 + n). Here n is the number of response function values computed during the t time interval. This previously motivated the development of a highly efficient fixed-trajectory implementation,[51] which is unsuited for large systems with disparate frequency scales. We present here a forward–backward implementation that is well suited for such systems. In the forward–backward implementation the number of trajectories and action transitions required to calculate the system response contribution from an OMT path scales as n. Because calculating 2DIR spectra requires a Fourier transform in t1 for a small number of waiting times, n ≫ n, so that the forward–backward implementation is significantly more efficient than the naive implementation. Here the forward–backward implementation is shown to be both accurate and practical for computing 2DIR spectra.

We have applied the forward–backward implementation of the OMT approximation to compute purely absorptive spectra for an anharmonic vibrational chromophore interacting with a harmonic bath with couplings that are linear in bath coordinates and either linear or quadratic in the chromophore coordinate. 2DIR spectra at t2 = 0 for this model with a continuum bath have been computed with a quantum Fokker–Planck equation approach by Ishizaki and Tanimura,[69] providing a basis for assessment of the OMT approach. In the pure dephasing limit of a weakly coupled off-resonant bath, our results agree both with the fluctuating frequency approximation, established to work well in this limit and with the quantum calculations of Ishizaki and Tanimura[69] for both forms of chromophore–bath interactions. Though this weak-coupling limit is theoretically simple, it poses a significant numerical challenge for the OMT. First, there is a large disparity in frequency scales, one of the challenges that motivated the development of the forward–backward implementation. Second, achieving correct time dependences of the response function requires extensive averaging over initial conditions and adequate treatment of transitions in system action at constant-angle variables. We have also computed spectra beyond this pure dephasing regime, in which energy transfer between system and bath influences the time and frequency dependences of the spectra. In this regime, our results agree qualitatively with the quantum calculations of ref (69). Although the interpretation of these spectra is more complex than in the pure dephasing limit, these OMT calculations are less numerically demanding, because the bath modes sample a smaller range of action values. We have used the OMT method to calculate waiting time dependences of 2DIR spectra. The line-shape dynamics in the pure dephasing regime are qualitatively similar to those obtained from the fluctuating frequency approximation. The t2 dependence in the energy transfer regime is plausible, but to our knowledge quantum calculations are not available for comparison.

Our results suggest the level of computational effort that will be required to apply the OMT approach to more general models, by identifying simplifying approximations that are likely to be valid. First, for purposes of generating initial classical states and for executing transitions, solvent motions weakly coupled to vibrational chromophores may be treated as independent. Second, depending on the spectral region of interest in the 2DIR spectrum, quantized action jumps need only be performed for a limited number of degrees of freedom. Last, for motions requiring these transitions, relatively crude approximations, such as low-order perturbation theory in anharmonicity are likely to suffice. Quantization of relatively low-frequency modes is unimportant, so that with the first simplifying approximation, low-frequency modes, such as solvent modes, may not need to be expressed in action–angle variables. Initial conditions could be sampled directly in Cartesian coordinates and momenta with constant angles maintained by simply fixing these variables during transitions. Thus, our findings have positive practical implications for the application of the OMT approach to anharmonic systems of a size that can be treated with conventional molecular dynamics simulations on the necessary time scales.