Two-dimensional Raman and hybrid terahertz-Raman spectroscopic techniques provide invaluable insight into molecular structures and dynamics of condensed-phase systems. However, corroborating experimental results with theory is difficult due to the high computational cost of incorporating quantum-mechanical effects in the simulations. Here, we present the equilibrium–nonequilibrium ring-polymer molecular dynamics (RPMD), a practical computational method that can account for nuclear quantum effects on the two-time response function of nonlinear optical spectroscopy. Unlike a recently developed approach based on the double Kubo transformed (DKT) correlation function, our method is exact in the classical limit, where it reduces to the established equilibrium-nonequilibrium classical molecular dynamics method. Using benchmark model calculations, we demonstrate the advantages of the equilibrium–nonequilibrium RPMD over classical and DKT-based approaches. Importantly, its derivation, which is based on the nonequilibrium RPMD, obviates the need for identifying an appropriate Kubo transformed correlation function and paves the way for applying real-time path-integral techniques to multidimensional spectroscopy.

Two-dimensional vibrational spectroscopy is a versatile technique to study microscopic interactions at their natural, femtosecond time scales.1 Recently, a series of hybrid spectroscopic experiments2–5 involving mid-infrared, far-infrared (or terahertz), and visible (Raman) pulses, have been developed to study electrical and mechanical anharmonicities,6–9 structural heterogeneities of liquids,10–12 and the couplings between intermolecular and intramolecular vibrational modes.13,14 However, the interpretation of such spectra is still an open question and requires adequate simulation methods.12,15–17

Several computational methods have been proposed for simulating two-dimensional off-resonant Raman and hybrid terhertz–Raman spectra. Model-based approaches aim at constructing simplified few-dimensional model systems that can be solved fully quantum mechanically in the presence of a harmonic oscillator bath.18,19 These approaches can disentangle contributions to spectra and relate them to different physical effects encoded in the model parameters. The parameters can be obtained by fitting to the existing experiment8 or molecular dynamics (MD) simulations.20 Alternatively, MD can be used to simulate spectra directly,6,21–31 as a way to cross-validate the simplified model and further test the molecular mechanics force field or the ab initio quantum chemistry method used for describing the forces between atoms, dipole moments, and polarizabilities.32–34 However, one improvement of MD is highly desired—MD simulations are based on classical atomic nuclei and neglect nuclear quantum effects, which can significantly modify the spectra. Indeed, recent two-dimensional Raman–THz–THz spectra of H2O and D2O revealed the effect of isotopic substitution beyond what would be expected from classical mechanics, pointing at nuclear quantum effects of the light hydrogen atoms.12 

Quantum simulations of the condensed phase have been enabled by various semiclassical methods based on the classical MD framework, including linearized semiclassical initial value representation (LSC-IVR),35,36 path-integral Liouville dynamics,37 centroid molecular dynamics (CMD),38 and ring-polymer molecular dynamics (RPMD).39 These methods have proven useful for the computation of one-time correlation functions40 related to various dynamical properties, including reaction rates,35,41–45 diffusion coefficients,46,47 and one-dimensional vibrational spectra.48–53 Some of them have also been applied to two-dimensional electronic54–58 and infrared spectroscopies,59–62 in which the quantum subsystem, consisting of, e.g., electronic or high-frequency vibrational degrees of freedom, can be well defined.63–65 In contrast, their application to two-dimensional Raman and hybrid terahertz-Raman spectroscopic techniques, in which all vibrational degrees of freedom should be treated on an equal footing, has been limited. Recently, the group of Batista66–68 proposed a set of methods that approximates the symmetric contribution to the so-called double Kubo transformed (DKT) correlation function, which is a two-time extension of the original one-time Kubo transformed correlation function.40 However, the relation between the symmetrized DKT correlation function and the spectroscopic response function is only approximate.

Here, we present an alternative RPMD approach to two-dimensional spectroscopy, which directly considers the two-time response function and does not rely on the additional approximation related to the DKT correlation function. To derive the proposed method, we start from the recently developed nonequilibrium RPMD69–73 and employ classical response theory. We then explore its validity and limitations both theoretically and numerically.

In two-dimensional Raman74,75 and hybrid terahertz-Raman spectroscopies,3–5 the signal is measured as a function of two delay times between three ultrashort light pulses, whose specific sequence determines the type of spectroscopy. To keep the discussion general, we consider the time-dependent Hamiltonian

Ĥtot=ĤÂF1(t)B̂F2(t),
(1)

comprised of the system’s field-free Hamiltonian Ĥ and two interactions with the first two ultrashort pulses, controlled, respectively, by coordinate-dependent operators  and B̂ and by time-dependent functions F1,2(t) representing the pulse shapes. For terahertz pulses, the interaction operators are the system’s dipole moments and functions F(t) are the electric fields of the light, whereas for the visible/near-infrared (Raman) pulses or sum-frequency terahertz excitation,9 the operators are the system’s polarizabilities and the time-dependent functions correspond to the squares of the electric fields. F1 is centered at t = −t1, F2 is centered at t = 0, and the signal is measured at t = t2, where t1 and t2 represent the time delays between the three light pulses. The recorded signal is proportional to the expectation value7 

Ĉ(t2)=Tr[Ĉρ̂(t2)]
(2)

of another coordinate-dependent operator Ĉ, which is again either the polarizability or dipole moment operator of the system. In Eq. (2), the system’s equilibrium density operator ρ̂=exp(βĤ)/Tr[exp(βĤ)] evolves under the time-dependent Hamiltonian (1), i.e.,

ρ̂(t)=Ûtot(t,)ρ̂Ûtot(t,),
(3)
Ûtot(tf,ti)=TexpititfdτĤtot(τ).
(4)

Here, β is the inverse temperature, and T is the time ordering operator. In practice, the experiments are designed to extract the components of the signal that are proportional to F1 and F2, which can be accounted for by evaluating the signal as8 

S(t2)=Ĉ(t2)++Ĉ(t2)+Ĉ(t2)++Ĉ(t2),
(5)

where ⟨·⟩jk represents the expectation value (2) of the system evolved under the Hamiltonian (1) with F1(t) → (j/2)F1(t) and F2(t) → (k/2)F2(t). For sufficiently weak external fields, we can further invoke the second-order time-dependent perturbation theory and recover the well-known result74 

S(t2)=0dτ20dτ1R(τ2,τ1)F2(t2τ2)F1(t2τ2τ1),
(6)

where

R(τ2,τ1)=12TrĈ(τ2+τ1)[B̂(τ1),[Â,ρ̂]]
(7)

is the corresponding response function, which depends only on the system’s properties. Because it can be easily convolved with different experimental pulses to recover the observed experimental signals (6), the response function is often the direct target of many theoretical studies. Here, we also intend to focus on the response function but take a slightly different route in developing the computational approach. Specifically, we note that for delta pulses, F1(t) = ɛ1δ(t + t1) and F2(t) = ɛ2δ(t), the signal measured after the second time delay t2 is proportional to the response function,

S(t2)=ε1ε2R(t2,t1),
(8)

where ɛ1,2 control the amplitude of the external fields. In the following, we derive an approximation to S(t2) under weak delta pulses and relate it to R(t2, t1) using Eq. (8).

To develop a tractable theory for condensed-phase simulations of spectra, we replace the exact quantum-mechanical expression (2) by its nonequilibrium RPMD69 approximation,

Ĉ(t)jkC(t)jkRP=dqdpCN(qjk,t)ρN(q,p),
(9)

where q and p are the positions and momenta of the extended system comprised of N replicas of the original D-dimensional system,

ρN(q,p)=eβNHN(q,p)dqdpeβNHN(q,p)
(10)

is the corresponding phase-space distribution,

HN(q,p)=i=1N12piTm1pi+12(βN)2(qiqi1)Tm(qiqi1)+V(qi)
(11)

is the system’s field-free ring-polymer Hamiltonian,40 βN = β/N,

CN(q)=1Ni=1NC(qi).
(12)

Here, m is the symmetric mass matrix of the system and q0 = qN. The classical time evolution of qjk,t and pjk,t is governed by the time-dependent ring-polymer Hamiltonian,

Hjk,N(q,p,t)=HN(q,p)+jVA,N(q)δ(t+t1)+kVB,N(q)δ(t),
(13)
VA,N(q)=ε12NAN(q),VB,N(q)=ε22NBN(q),
(14)

with j, k ∈ {+, −}. Nonequilibrium RPMD, in which the dynamics and initial distribution depend on different Hamiltonians [see Eqs. (9)(11) and (13)], has been rigorously justified as an approximation to the more general real-time nonequilibrium Matsubara dynamics.69 Although its original derivation involved only time-independent Hamiltonians, this assumption was not used explicitly at any step of the derivation, justifying the use of nonequilibrium RPMD even in the time-dependent setting. At this point, we recognize that the equations above already formulate a valid computational technique for evaluating the response function, which, in the limit of N = 1, is equivalent to the finite-field nonequilibrium MD method of Jansen, Snijders, and Duppen.21,22 In contrast to their classical approach, the nonequilibrium RPMD theory can include, at least approximately, the nuclear quantum effects on two-dimensional spectra when the path-integral continuum limit is reached. However, these nonequilibrium (RP)MD methods are impractical because a different set of trajectories is needed for each choice of the delay time t1. To address this problem, we employ the equilibrium–nonequilibrium approach, originally developed by Hasegawa and Tanimura26 in the context of classical MD simulations, in which one of the two pulses is treated perturbatively, and combine it with the RPMD to account for nuclear quantum effects.

We start by rewriting Eq. (9) as

C(t)jkRP=dqdpCN(q)Ujk,N(t,)ρN(q,p),
(15)

where

Ujk,N(tf,ti)=TexptitfdτLjk,N(τ)
(16)

governs the time evolution under the Hamiltonian Hjk,N,

Ljk,N(t)={,Hjk,N(t)}=LN+jδ(t+t1)LA,N+kδ(t)LB,N
(17)

is the corresponding Liouvillian, LN = {·, HN}, LA,N = {·, VA,N}, LB,N = {·, VB,N}, and {·, ·} denotes the Poisson bracket. Next, we expand76,77

Ujk,N(t,)Uk,N(t,)jtdτUk,N(t,τ)LA,Nδ(τ+t1)Uk,N(τ,)
(18)
=Uk,N(t,)jUk,N(t,t1)LA,NUN(t1,)
(19)

to the first order in ɛ1 to obtain

C(t)jkRPC(t)kRP,(0)+C(t)jkRP,(1),
(20)

where

C(t)kRP,(0)=dqdpCN(q)Uk,N(t,)ρN(q,p)
(21)

and

C(t)jkRP,(1)=jdqdpCN(q)Uk,N(t,t1)LA,N×UN(t1,)ρN(q,p).
(22)

In Eq. (18), we introduced a shorthand notation

Uk,N(tf,ti)=Texptitfdτ[LN+kδ(τ)LB,N]
(23)

for the evolution operator that involves only the second pulse. In going from Eq. (18) to Eq. (19), we assumed that t, t1 > 0 and used

Uk,N(tf,ti)=UN(tf,ti),0(ti,tf),Uk,N(tf,0)UN(0,ti)otherwise,
(24)

with

UN(tf,ti)=eLN(tfti).
(25)

Then, we derive

C(t)jkRP,(1)=jβε12dqdpCN(qk,t)ȦN(qt1)ρN(q,p)
(26)

from Eq. (22) by applying UN(−t1, −)ρN(q, p) = ρN(q, p),

LA,NρN(q,p)=VA,N(q)qTpρN(q,p)
(27)
=βNε12i=1NA(qi)qiTq̇iρN(q,p)
(28)
=βε12ȦN(q)ρN(q,p),
(29)

where ḟ=f/t denotes the time derivative, and Eq. (24) in combination with CN(qk,t) = CN(q)Uk,N(t, 0) and ȦN(qt1)=UN(0,t1)ȦN(q). qt1 is the position along a backward equilibrium trajectory. qk,t2 corresponds to a nonequilibrium trajectory evolved under the Hamiltonian that involves the interaction with the second pulse, which is equivalent to the field-free evolution with the initial conditions qk,0 = q and pk,0 = p + (2/2)N [∂BN(q)/∂q].

Combining Eq. (26) with Eqs. (5), (8), and (20) yields the final result

R(t2,t1)βε2dqdpCN(q+,t2)CN(q,t2)ȦN(qt1)ρN(q,p),
(30)

assuming ɛ2 is sufficiently small to guarantee the validity of the perturbative expression (8). It follows that the response function R(t2, t1) can be evaluated as an equilibrium ensemble average of an estimator based on two nonequilibrium trajectories (q±,t2) and one backward equilibrium trajectory (qt1). At high temperatures, where β → 0, or for N = 1, Eq. (30) reduces to the equilibrium–nonequilibrium classical MD approach.26 In contrast to the classical approach, the RPMD method captures nuclear quantum effects in accordance with the well-known advantages and limitations of the RPMD for one-time correlation functions.40 Furthermore, in the supplementary material, we show that the RPMD is exact for the harmonic oscillator, when two out of the three operators Â, B̂, and Ĉ are linear functions of position. Finally, the proposed equilibrium–nonequilibrium RPMD method always satisfies R(t2 = 0, t1) = 0 (because q+,0 = q−,0), which holds for the exact quantum response function (7), but is not guaranteed by some other approximate methods (as discussed below).

Before demonstrating these properties numerically, we briefly review the only other RPMD-based method proposed to date for simulating two-dimensional vibrational spectra. In 2018, Jung, Videla, and Batista66 showed that the response function (7) is related to the DKT correlation function,78 

K(t2,t1)=1β20βdλ10λ1dλ2×TrρÂ(iλ1)B̂(t1iλ2)Ĉ(t1+t2),
(31)

through the frequency-domain expression

R(ω2,ω1)=Q+(ω2,ω1)Ksym(ω2,ω1)+Q(ω2,ω1)Kasym(ω2,ω1),
(32)

where

Ksym(t2,t1)=2Re[K(t2,t1)],Kasym(t2,t1)=2iIm[K(t2,t1)],
(33)
f(ω2,ω1)=dt1dt2f(t2,t1)eiω1t1iω2t2,
(34)

and Q± are some known functions of ω1 and ω2.78 Although Eq. (32) is exact, an additional approximation is needed before K(t2, t1) can be evaluated using the RPMD. Specifically, the authors employed the symmetrized DKT approximation

R(ω2,ω1)Q+(ω2,ω1)Ksym(ω2,ω1),
(35)

neglecting the asymmetric contribution to the DKT correlation function, to make use of the RPMD approximation

Ksym(t2,t1)dqdpAN(q)BN(qt1)CN(qt1+t2)ρN(q,p).
(36)

Thus, the symmetrized DKT approximation, and, as a consequence, the RPMD DKT method, is not exact in the classical limit. In addition, the corresponding response functions need not be zero at t2 = 0,78 except in the high-temperature limit, in which case the RPMD DKT method reduces to the classical correlation function approach of DeVane, Ridley, Space, and Keyes.25,79–81

In this section, the above considerations about the classical, RPMD, and RPMD DKT methods for simulating two-dimensional vibrational spectra are demonstrated numerically on an anharmonic model system determined by the Hamiltonian

H1D(q,p)=12(p2+q2)+aq3+a2q4,
(37)

where a is a tunable parameter that controls the degree of anharmonicity. This system can be solved numerically exactly in a finite basis and was previously used with a = 0.1 to evaluate the accuracy of RPMD and nonequilibrium RPMD methods.39,69 Details of the exact quantum and approximate trajectory-based simulations can be found in the supplementary material.

Figure 1 shows that the equilibrium-nonequilibrium RPMD approach performs better than the RPMD DKT method at both low and high temperatures. Notably, the newly proposed RPMD approach exhibits the expected short-time accuracy, while most of the errors in the RPMD DKT method can be attributed to the neglect of the asymmetric part of the full DKT correlation function (compare columns 3 and 4 of Fig. 1).

FIG. 1.

Errors of approximate methods (columns 2–4) for evaluating the two-time response function (shown in the left-most column): New equilibrium–nonequilibrium RPMD method (labeled “RPMD”), RPMD DKT approximation,66 and the exact quantum symmetrized DKT correlation function (“Symmetrized DKT”).66 Results are shown at both high (β = 1, top) and low (β = 8, bottom) temperatures, with the anharmonicity parameter a = 0.1 [Eq. (37)], and operators Â=B̂=q̂ and Ĉ=q̂2/2.

FIG. 1.

Errors of approximate methods (columns 2–4) for evaluating the two-time response function (shown in the left-most column): New equilibrium–nonequilibrium RPMD method (labeled “RPMD”), RPMD DKT approximation,66 and the exact quantum symmetrized DKT correlation function (“Symmetrized DKT”).66 Results are shown at both high (β = 1, top) and low (β = 8, bottom) temperatures, with the anharmonicity parameter a = 0.1 [Eq. (37)], and operators Â=B̂=q̂ and Ĉ=q̂2/2.

Close modal

Figure 2 shows the time slices of R(t2, t1) along t2 (left panels) and t1 (right panels). At low temperature (Fig. 2, bottom), equilibrium–nonequilibrium RPMD is more accurate than the classical method due to the inclusion of nuclear quantum effects, while the two are similar at high temperatures (Fig. 2, top). Although the RPMD response function agrees with the exact result at short time (see Fig. 2, bottom left), its long-time oscillation frequency deviates from the exact because RPMD is not expected to recover quantum coherence effects. In this example, the classical simulation is more accurate than the RPMD DKT method, showing that the uncontrolled error of the symmetrized DKT approximation can outweigh the benefits of including nuclear quantum effects. In fact, similar findings were implied, even if not explicitly discussed, by the earlier reports on the symmetrized DKT approximation.66,78 In addition, both classical and RPMD equilibrium–nonequilibrium methods are exact (zero) at t2 = 0 (see the left panels of Fig. 2), which does not hold for the RPMD DKT approach [see also Figs. 6(e) and 8(e) of Ref. 78].

FIG. 2.

Cuts along t1 = 10 (left) and t2 = 10 (right) of the two-dimensional response function at β = 1 (top) and β = 8 (bottom) evaluated with the exact quantum approach, with the RPMD and classical equilibrium–nonequilibrium methods, and with the RPMD DKT approach. System’s parameters and operators were the same as in Fig. 1.

FIG. 2.

Cuts along t1 = 10 (left) and t2 = 10 (right) of the two-dimensional response function at β = 1 (top) and β = 8 (bottom) evaluated with the exact quantum approach, with the RPMD and classical equilibrium–nonequilibrium methods, and with the RPMD DKT approach. System’s parameters and operators were the same as in Fig. 1.

Close modal

To demonstrate numerically that the proposed approach inherits some well-known properties of RPMD, we choose to quantify the error of an approximate response function Rapprox as

Error=RapproxRexactmax(Rapprox,Rexact),
(38)

where ‖R2 = ∫dt1∫dt2|R (t2, t1)|2 and Rexact is the exact response function. Figure 3 shows that the equilibrium–nonequilibrium RPMD consistently outperforms the classical approach for systems with different degrees of anharmonicity (a), at different temperatures (b), and for different combinations of operators (c). As discussed earlier, the RPMD method is exact for the harmonic oscillator [a = 0, see Fig. 3(a)] and certain combinations of operators Â, B̂, and Ĉ. Interestingly, the classical approach is also exact in some of those limiting cases (see the supplementary material). Nevertheless, as soon as the system is even weakly anharmonic, RPMD is clearly superior to the classical approach. As expected, both RPMD and classical methods perform worse as the anharmonicity increases. Furthermore, the equilibrium–nonequilibrium RPMD method converges to the exact (classical) result in the high-temperature limit [β → 0, Fig. 3(b)]. Again, as the temperature decreases, the RPMD approach becomes less accurate. Finally, the proposed RPMD approach is more accurate for linear compared to nonlinear operators [Fig. 3(c)].

FIG. 3.

Errors based on Eq. (38) of the approximate two-time response functions for systems with different degrees of anharmonicity (a), at different temperatures (b), and for different combinations of linear and quadratic operators Â, B̂, and Ĉ (c). Apart from the parameters that are being varied or explicitly indicated in each plot, the default parameters are a = 0.1 [Eq. (37)], β = 8, Â=B̂=q̂, and Ĉ=q̂2/2. In panel (c), “1” denotes linear operator q̂ and “2” denotes quadratic operator q̂2/2. For example, “122” means Â=q̂ and B̂=Ĉ=q̂2/2. Operator combinations are ordered according to the error of the RPMD simulation at β = 8.

FIG. 3.

Errors based on Eq. (38) of the approximate two-time response functions for systems with different degrees of anharmonicity (a), at different temperatures (b), and for different combinations of linear and quadratic operators Â, B̂, and Ĉ (c). Apart from the parameters that are being varied or explicitly indicated in each plot, the default parameters are a = 0.1 [Eq. (37)], β = 8, Â=B̂=q̂, and Ĉ=q̂2/2. In panel (c), “1” denotes linear operator q̂ and “2” denotes quadratic operator q̂2/2. For example, “122” means Â=q̂ and B̂=Ĉ=q̂2/2. Operator combinations are ordered according to the error of the RPMD simulation at β = 8.

Close modal

We further consider a two-dimensional model Hamiltonian

H2D(q1,p1,q2,p2)=H1D(Ω1q1,p1)+H1D(Ω2q2,p2)+λq1q2
(39)

composed of two anharmonic oscillators [Eq. (37), a = 0.2] with different central frequencies (Ω1 = 0.5, Ω2 = 2) and a linear–linear coupling term proportional to λ = 0.1. The light–matter interaction operators are set to Â=q̂1 and B̂=Ĉ=q̂2, reflecting the terahertz-infrared-Raman pulse sequence,5,14 in which the first terahertz pulse interacts only with the low-frequency mode q1, while the infrared and off-resonant Raman interactions probe the high-frequency mode q2. Figure 4 shows that the two-dimensional spectrum, obtained as the double cosine transform of the two-time response function (see the supplementary material for further details), exhibits off-diagonal peaks at (Ω1, Ω2) and (Ω1, Ω2 ± Ω1), as predicted by the Feynman diagrams of Refs. 14 and 17. The signal vanishes in the absence of coupling, λ = 0, or anharmonicity, a = 0 (see Sec. III of the supplementary material). Equilibrium–nonequilibrium RPMD and classical MD spectra reproduce the shape of the exact spectrum. However, at low temperature (Fig. 4, bottom), the classical method largely overestimates the magnitudes of the off-diagonal peaks due to the neglect of the nuclear quantum effects. In our simplified case, this results mostly in an overall scaling factor but, in general, could affect the relative intensities of the spectral features.

FIG. 4.

Exact (left), equilibrium–nonequilibrium RPMD (center), and equilibrium–nonequilibrium classical MD (right), two-dimensional spectra of the two-dimensional model system described by the Hamiltonian (39) with Â=q̂1 and B̂=Ĉ=q̂2, simulated at high (β = 1, top) and low (β = 8, bottom) temperatures.

FIG. 4.

Exact (left), equilibrium–nonequilibrium RPMD (center), and equilibrium–nonequilibrium classical MD (right), two-dimensional spectra of the two-dimensional model system described by the Hamiltonian (39) with Â=q̂1 and B̂=Ĉ=q̂2, simulated at high (β = 1, top) and low (β = 8, bottom) temperatures.

Close modal

To conclude, we have introduced an RPMD-based method that captures nuclear quantum effects on three-pulse two-dimensional vibrational spectra at the cost of classical MD simulations. The proposed approach was shown to perform better than the classical and recently developed RPMD DKT methods. Although our discussion focused only on RPMD, the overall scheme could be extended to the more general Matsubara dynamics82–84 and to its efficient approximations, including thermostated RPMD,49 CMD,38 and quasi-CMD.51 It would be interesting to see how these approximations perform when applied to multidimensional spectroscopy because it is expected to provide a more stringent test for dynamical approximations compared to linear spectroscopy, for which similar studies exist.50 Because two-dimensional Raman and hybrid THz–Raman spectra vanish for harmonic potentials with linear dipole moment or polarizability operators,74 adequate computational approaches must capture the effect of electrical and mechanical anharmonicities, which poses a challenge for RPMD and the above-mentioned methods.77,85

To derive the equilibrium–nonequilibrium RPMD, we did not invoke the concept of the Kubo transformation, which is traditionally done for one-time correlation functions,40 but instead combined nonequilibrium RPMD with classical response theory. This new way of deriving real-time path-integral methods will be of special interest for simulating time-resolved spectroscopic experiments, but also in other contexts where multi-time correlation or response functions appear or if a derivation of the appropriate Kubo transformed correlation function is not obvious. The theory is general and enables new applications of established quantum dynamical methods.

See the supplementary material for the computational details, analytical expressions in the harmonic limit, and further discussion of the two-mode system in the absence of coupling or anharmonicity.

The authors thank Kenneth A. Jung, Roman Korol, and Jorge L. Rosa-Raíces for helpful discussions. T.B. acknowledges financial support from the Swiss National Science Foundation through the Early Postdoc Mobility Fellowship (Grant No. P2ELP2-199757). G.A.B. and T.F.M. gratefully acknowledge support from the National Science Foundation Chemical Structure, Dynamics and Mechanisms program (Grant No. CHE-1665467). The computations presented here were conducted in the Resnick High Performance Computing Center, a facility supported by the Resnick Sustainability Institute at the California Institute of Technology.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
P.
Hamm
and
M.
Zanni
,
Concepts and Methods of 2D Infrared Spectroscopy
(
Cambridge University Press
,
2011
).
2.
P.
Hamm
and
J.
Savolainen
,
J. Chem. Phys.
136
,
094516
(
2012
).
3.
J.
Savolainen
,
S.
Ahmed
, and
P.
Hamm
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
20402
(
2013
).
4.
I. A.
Finneran
,
R.
Welsch
,
M. A.
Allodi
,
T. F.
Miller
, and
G. A.
Blake
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
6857
(
2016
).
5.
M.
Grechko
,
T.
Hasegawa
,
F.
D’Angelo
,
H.
Ito
,
D.
Turchinovich
,
Y.
Nagata
, and
M.
Bonn
,
Nat. Commun.
9
,
885
(
2018
).
6.
P.
Hamm
,
J. Chem. Phys.
141
,
184201
(
2014
).
7.
I. A.
Finneran
,
R.
Welsch
,
M. A.
Allodi
,
T. F.
Miller
, and
G. A.
Blake
,
J. Phys. Chem. Lett.
8
,
4640
(
2017
).
8.
I. B.
Magdǎu
,
G. J.
Mead
,
G. A.
Blake
, and
T. F.
Miller
,
J. Phys. Chem. A
123
,
7278
(
2019
).
9.
G.
Mead
,
H.-W.
Lin
,
I.-B.
Magdău
,
T. F.
Miller
, and
G. A.
Blake
,
J. Phys. Chem. B
124
,
8904
(
2020
).
10.
A.
Shalit
,
S.
Ahmed
,
J.
Savolainen
, and
P.
Hamm
,
Nat. Chem.
9
,
273
(
2017
).
11.
P.
Hamm
and
A.
Shalit
,
J. Chem. Phys.
146
,
130901
(
2017
).
12.
A.
Berger
,
G.
Ciardi
,
D.
Sidler
,
P.
Hamm
, and
A.
Shalit
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
2458
(
2019
).
13.
G.
Ciardi
,
A.
Berger
,
P.
Hamm
, and
A.
Shalit
,
J. Phys. Chem. Lett.
10
,
4463
(
2019
).
14.
L.
Vietze
,
E. H. G.
Backus
,
M.
Bonn
, and
M.
Grechko
,
J. Chem. Phys.
154
,
174201
(
2021
).
15.
P.
Hamm
,
J. Chem. Phys.
151
,
054505
(
2019
).
16.
D.
Sidler
and
P.
Hamm
,
J. Chem. Phys.
150
,
044202
(
2019
).
17.
D.
Sidler
and
P.
Hamm
,
J. Chem. Phys.
153
,
044502
(
2020
).
18.
Y.
Tanimura
and
A.
Ishizaki
,
Acc. Chem. Res.
42
,
1270
(
2009
).
19.
T.
Ikeda
,
H.
Ito
, and
Y.
Tanimura
,
J. Chem. Phys.
142
,
212421
(
2015
).
20.
S.
Ueno
and
Y.
Tanimura
,
J. Chem. Theory Comput.
16
,
2099
(
2020
).
21.
T. l. C.
Jansen
,
J. G.
Snijders
, and
K.
Duppen
,
J. Chem. Phys.
113
,
307
(
2000
).
22.
T. l. C.
Jansen
,
J. G.
Snijders
, and
K.
Duppen
,
J. Chem. Phys.
114
,
10910
(
2001
).
23.
S.
Saito
and
I.
Ohmine
,
Phys. Rev. Lett.
88
,
207401
(
2002
).
24.
S.
Saito
and
I.
Ohmine
,
J. Chem. Phys.
119
,
9073
(
2003
).
25.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
J. Chem. Phys.
119
,
6073
(
2003
).
26.
T.
Hasegawa
and
Y.
Tanimura
,
J. Chem. Phys.
125
,
074512
(
2006
).
27.
T.
Hasegawa
and
Y.
Tanimura
,
J. Chem. Phys.
128
,
064511
(
2008
).
28.
T.
Yagasaki
and
S.
Saito
,
Acc. Chem. Res.
42
,
1250
(
2009
).
29.
H.
Ito
,
T.
Hasegawa
, and
Y.
Tanimura
,
J. Chem. Phys.
141
,
124503
(
2014
).
30.
X.
Sun
,
J. Chem. Phys.
151
,
194507
(
2019
).
31.
T. L. C.
Jansen
,
S.
Saito
,
J.
Jeon
, and
M.
Cho
,
J. Chem. Phys.
150
,
100901
(
2019
).
32.
T.
Hasegawa
and
Y.
Tanimura
,
J. Phys. Chem. B
115
,
5545
(
2011
).
33.
H.
Ito
,
T.
Hasegawa
, and
Y.
Tanimura
,
J. Phys. Chem. Lett.
7
,
4147
(
2016
).
34.
D.
Sidler
,
M.
Meuwly
, and
P.
Hamm
,
J. Chem. Phys.
148
,
244504
(
2018
).
35.
H.
Wang
,
X.
Sun
, and
W. H.
Miller
,
J. Chem. Phys.
108
,
9726
(
1998
).
36.
J.
Liu
and
W. H.
Miller
,
J. Chem. Phys.
125
,
224104
(
2006
).
37.
J.
Liu
,
J. Chem. Phys.
140
,
224107
(
2014
).
38.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
111
,
2371
(
1999
).
39.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
40.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
41.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
112
,
8747
(
2000
).
42.
N.
Boekelheide
,
R.
Salomón-Ferrer
, and
T. F.
Miller
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
16159
(
2011
).
43.
Y. V.
Suleimanov
,
F. J.
Aoiz
, and
H.
Guo
,
J. Phys. Chem. A
120
,
8488
(
2016
).
44.
X.
Tao
,
P.
Shushkov
, and
T. F.
Miller
,
J. Phys. Chem. A
123
,
3013
(
2019
).
45.
X.
Tao
,
P.
Shushkov
, and
T. F.
Miller
,
J. Chem. Phys.
152
,
124117
(
2020
).
46.
T. F.
Miller
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
184503
(
2005
).
47.
T. F.
Miller
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
154504
(
2005
).
48.
A.
Witt
,
S. D.
Ivanov
,
M.
Shiga
,
H.
Forbert
, and
D.
Marx
,
J. Chem. Phys.
130
,
194510
(
2009
).
49.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
140
,
234116
(
2014
).
50.
R. L.
Benson
,
G.
Trenins
, and
S. C.
Althorpe
,
Faraday Discuss.
221
,
350
(
2019
).
51.
G.
Trenins
,
M. J.
Willatt
, and
S. C.
Althorpe
,
J. Chem. Phys.
151
,
054109
(
2019
).
52.
R.
Korol
,
J. L.
Rosa-Raíces
,
N.
Bou-Rabee
, and
T. F.
Miller
,
J. Chem. Phys.
152
,
104102
(
2020
).
53.
J. L.
Rosa-Raíces
,
J.
Sun
,
N.
Bou-Rabee
, and
T. F.
Miller
,
J. Chem. Phys.
154
,
024106
(
2021
).
54.
R. F.
Loring
,
J. Chem. Phys.
146
,
144106
(
2017
).
55.
K.
Polley
and
R. F.
Loring
,
J. Phys. Chem. B
124
,
9913
(
2020
).
56.
J.
Provazza
,
F.
Segatta
,
M.
Garavelli
, and
D. F.
Coker
,
J. Chem. Theory Comput.
14
,
856
(
2018
).
57.
J.
Provazza
,
F.
Segatta
, and
D. F.
Coker
,
J. Chem. Theory Comput.
17
,
29
(
2020
).
58.
X.
Gao
and
E.
Geva
,
J. Chem. Theory Comput.
16
,
6491
(
2020
).
59.
S. M.
Gruenbaum
and
R. F.
Loring
,
J. Chem. Phys.
131
,
204504
(
2009
).
60.
M.
Gerace
and
R. F.
Loring
,
J. Chem. Phys.
138
,
124104
(
2013
).
61.
M.
Alemi
and
R. F.
Loring
,
J. Chem. Phys.
142
,
212417
(
2015
).
62.
K.
Kwac
and
E.
Geva
,
J. Phys. Chem. B
117
,
7737
(
2013
).
63.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
129
,
124505
(
2008
).
64.
P. L.
McRobbie
,
G.
Hanna
,
Q.
Shi
, and
E.
Geva
,
Acc. Chem. Res.
42
,
1299
(
2009
).
65.
P. L.
McRobbie
and
E.
Geva
,
J. Phys. Chem. A
113
,
10425
(
2009
).
66.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
148
,
244105
(
2018
).
67.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
151
,
034108
(
2019
).
68.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
153
,
124112
(
2020
).
69.
R.
Welsch
,
K.
Song
,
Q.
Shi
,
S. C.
Althorpe
, and
T. F.
Miller
,
J. Chem. Phys.
145
,
204118
(
2016
).
70.
A.
Marjollet
and
R.
Welsch
,
J. Chem. Phys.
152
,
194113
(
2020
).
71.
A.
Marjollet
and
R.
Welsch
,
Int. J. Quantum Chem.
121
,
e26447
(
2021
).
72.
H.
Jiang
,
X.
Tao
,
M.
Kammler
,
F.
Ding
,
A. M.
Wodtke
,
A.
Kandratsenka
,
T. F.
Miller
, and
O.
Bünermann
,
J. Phys. Chem. Lett.
12
,
1991
(
2021
).
73.
A.
Marjollet
,
L.
Inhester
, and
R.
Welsch
,
J. Chem. Phys.
156
,
044101
(
2022
).
74.
Y.
Tanimura
and
S.
Mukamel
,
J. Chem. Phys.
99
,
9496
(
1993
).
75.
S.
Palese
,
J. T.
Buontempo
,
L.
Schilling
,
W. T.
Lotshaw
,
Y.
Tanimura
,
S.
Mukamel
, and
R. J. D.
Miller
,
J. Phys. Chem.
98
,
12466
(
1994
).
76.
B. L.
Holian
and
D. J.
Evans
,
J. Chem. Phys.
83
,
3560
(
1985
).
77.
T.
Plé
,
S.
Huppert
,
F.
Finocchi
,
P.
Depondt
, and
S.
Bonella
,
J. Chem. Phys.
155
,
104108
(
2021
).
78.
Z.
Tong
,
P. E.
Videla
,
K. A.
Jung
,
V. S.
Batista
, and
X.
Sun
,
J. Chem. Phys.
153
,
034117
(
2020
).
79.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
Phys. Rev. E
70
,
050101(R)
(
2004
).
80.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
J. Chem. Phys.
123
,
194507
(
2005
).
81.
R.
DeVane
,
C.
Kasprzyk
,
B.
Space
, and
T.
Keyes
,
J. Phys. Chem. B
110
,
3773
(
2006
).
82.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
,
J. Chem. Phys.
142
,
134103
(
2015
).
83.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
,
J. Chem. Phys.
142
,
191101
(
2015
).
84.
85.
R. L.
Benson
and
S. C.
Althorpe
,
J. Chem. Phys.
155
,
104107
(
2021
).

Supplementary Material