Two-dimensional (2D) vibrational spectroscopy has emerged as one of the most important experimental techniques useful to study the molecular structure and dynamics in condensed phases. Theory and computation have also played essential and integral roles in its development through the nonlinear optical response theory and computational methods such as molecular dynamics (MD) simulations and electronic structure calculations. In this article, we present the fundamental theory of coherent 2D vibrational spectroscopy and describe computational approaches to simulate the 2D vibrational spectra. The classical approximation to the quantum mechanical nonlinear response function is invoked from the outset. It is shown that the third-order response function can be evaluated in that classical limit by using equilibrium or non-equilibrium MD simulation trajectories. Another simulation method is based on the assumptions that the molecular vibrations can still be described quantum mechanically and that the relevant molecular response functions are evaluated by the numerical integration of the Schrödinger equation. A few application examples are presented to help the researchers in this and related areas to understand the fundamental principles and to use these methods for their studies with 2D vibrational spectroscopic techniques. In summary, this exposition provides an overview of current theoretical efforts to understand the 2D vibrational spectra and an outlook for future developments.

Resonance is a phenomenon ubiquitous in nature. Molecular spectroscopy is an experimental tool that utilizes one of the resonance phenomena that involves an interaction of oscillating charged particles in a given molecule with an external electromagnetic field whose frequency (ωwave) is close to that (ωmolecule) of a vibrational or electronic oscillation determined by the associated interaction potential between constituent charged particles in a given polyatomic molecule in condensed phases. Such resonant field-matter interaction is manifested in the complex quantum transition amplitude that is approximately proportional to (ωwaveωmolecule + molecule)−1, where γmolecule determines the linewidth of such a transition band that reflects complicated molecular interactions with fluctuating internal and bath degrees of freedom. As the wave frequency ωwave approaches the molecular oscillation frequency ωmolecule, the corresponding quantum transition amplitude increases, which has been referred to as a quantum resonant field-matter interaction process.

Vibrational spectroscopy, one of the molecular spectroscopic techniques, involves quantum transitions of vibrational degrees of freedom. The most widely used infrared (IR) spectroscopy measures the vibrational transition amplitude of normal modes via their resonant interaction with infrared (e.g., near IR, mid-IR, far-IR, and THz) radiation. The vibrational Raman spectroscopy is a resonance-enhanced scattering of the incident electronically non-resonant field when the beat frequency between the incident radiation field and the inelastic Raman scattering field is close to molecular vibrational frequency. Although IR-vis sum-frequency-generation (SFG) spectroscopy for molecular systems on the surface or at the interface with broken centrosymmetry is a nonlinear optical (three-wave-mixing) spectroscopy, it still involves a single vibrational transition whose amplitude is essentially determined by the quantum correlation between the transition dipole and the transition polarizability of a given vibrational mode.

Coherent multidimensional vibrational spectroscopy1–5 is an extension of these linear spectroscopic methods and involves multiple vibrational transitions that are temporally separated. To observe and analyze the correlation between the distinctive vibrational transitions, one usually applies a series of coherent laser pulses with specific phase relations to induce such transitions. The multiple transitions are detected and presented as a function of corresponding frequency variables, making the spectrum multidimensional. In particular, coherent two-dimensional (2D) vibrational spectroscopy employs three femtosecond laser pulses in the infrared (IR) or visible frequency range to induce the third-order polarization in molecular systems. The signal due to the polarization is measured and presented in two frequency dimensions conjugates to the time intervals between the first and second pulses (τ) and between the third pulse and the detection (t), respectively. A series of 2D spectra at different time intervals between the second and the third pulse (T) reveals the molecular dynamics (MD) taking place during that time interval T that is called the waiting or population time.

2D vibrational spectroscopy provides ultrafast time resolution along the waiting time (T) axis as well as high spectral resolution of individual peaks along the diagonal in the 2D frequency space and their couplings in the form of off-diagonal cross peaks. Therefore, the method conveys rich information on molecular systems such as homogeneous (anti-diagonal) and inhomogeneous (diagonal) spectral broadening effects, vibrational anharmonicity, spectral diffusion,6 and intermode coupling strength and its temporal variation. Over the past two decades, 2D vibrational spectroscopy has been extensively used to study the structure and dynamics of small peptides, proteins, DNA, and lipid bilayers, energy transfer dynamics, hydrogen-bonding (H-bonding) structure and dynamics of liquid water, and its isotopologues, and configurational and H-bonding dynamics of biomolecules.

The amplitudes of IR and Raman transitions critically depend on the transition dipole and polarizability, respectively. Since the molecular moments m generally depend on the vibrational coordinate q, they can be Taylor-expanded in the form m(t) = m0 + (m/∂q)0q + with respect to the equilibrium point denoted by subscript 0. Then, when the rotational motion of each vibrational chromophore is assumed to be very slow compared to vibrational dephasing rates, according to the time correlation function formalism of molecular spectroscopy,7 the corresponding line shapes can be approximately written as
(1)
where (μ/∂q)0 and (α/∂q)0 are the transition dipole and polarizability, respectively, and the bar over, for example, μ/q02¯ denotes the orientational average. Due to the first-order coordinate dependences of transition dipole and transition polarizability, linear spectra (IR, Raman, IR-vis SFG), which are in principle related to the dipole-dipole, polarizability-polarizability, dipole-polarizability correlation functions, are determined by the vibrational coordinate-coordinate correlation function q(t)q(0).8 

For coupled multi-oscillator systems, vibrational couplings through space via intermolecular interactions or through bond via anharmonicities in the multi-dimensional potential energy surface are crucial in understanding vibrational dynamics.9,10 However, since the linear vibrational spectrum is mainly determined by the harmonic properties, e.g., vibrational frequency and transition dipole moment, of the oscillators, it is difficult to quantitatively extract such weak features, e.g., vibrational coupling constants and potential anharmonic coefficients of coupled oscillators in condensed phases, from those linear spectra.8 On the other hand, the coherent multi-dimensional vibrational spectroscopic methods have been found to be of exceptional use because they enable one to explore these mode-mode couplings and provide crucial information on the structural dynamics through the time-dependent changes of the spatial proximity and relative orientation of chromophores. Note that the nonlinear vibrational response functions determining amplitudes and lineshapes of coherent multidimensional vibrational spectra vanish for perfect harmonic oscillators.9,11–14 One of the most popular techniques is the 2D IR spectroscopy,1 which is a four-wave-mixing method capable of providing information on the molecular nonlinear response function that is given by multi-time correlation functions of transition dipole moments. Again, using the linear expansion form of electric dipole moment with respect to vibrational coordinates and employing perturbation theory treating electric and mechanical anharmonicities as weak perturbations to harmonic oscillators, it was possible to rewrite the nonlinear response functions in terms of vibrational Green functions and such perturbation terms.11,12,15,16 This method was known as a linked diagram theory for the multidimensional vibrational response function. Although such theoretical approaches were popular in the early time of theory on multidimensional vibrational spectroscopy, the intrinsic difficulties in describing the anharmonicity-induced frequency shift of the excited state absorption contribution to the 2D IR spectrum, for example, and in accurately calculating vibrational (both mechanical and electrical) anharmonicities for oscillators in solutions, prohibited its wide use later.

Although 2D IR spectroscopy employing three incident IR laser pulses is one of the most widely used forms of coherent multidimensional vibrational spectroscopy, different varieties of the method have also been developed and applied for specific purposes. Examples include surface specific 2D SFG spectroscopy,17,18 2D Raman19 and terahertz spectroscopy,20–24 2D IR-IR-visible spectroscopy,3,9,13,25–29 and so on. In this article, because the basic theories for these variants are similar to one another, we mainly focus on the theory and computation of 2D IR spectroscopy and reference is made to these variants when appropriate.

By necessity, the exposition in this paper cannot be exhaustive. However, it is hoped to serve as an up-to-date introduction and a basis for elaborate research especially with the aid of references provided. In Sec. II, we briefly introduce the third-order response function formalism, which has the central importance in the description of 2D vibrational spectroscopy. In Sec. III, classical mechanical approaches to calculate 2D vibrational spectra are presented based on the classical limit of the quantum vibrational third-order response function. In Sec. IV, we present a quantum mechanical (QM) method to simulate the spectra by numerical integration of the Schrödinger equation (NISE). Finally, we conclude with perspective on future development.

In general, spectroscopic measurements are conducted in two steps: the preparation step where the molecular system is excited by incident radiations and the detection step where the generated signal is measured and presented. In 2D vibrational spectroscopy, the system is irradiated with three coherent laser pulses, usually in the infrared (IR) frequency region, in the preparation step and then the generated signal field is detected and presented in two frequency dimensions, representing two distinct coherence oscillations separated by a waiting (population) time. It is a kind of four-wave-mixing spectroscopy due to the fact that the signal field arises from three preceding field-matter interactions that are each linear in the field. In each of the four field-matter interaction events, a quantum transition takes place between vibrational states of the system, on either the ket or the bra side of the system density operator (see below). Depending on the configuration of the optical laser pulses such as the frequency, the direction of propagation (wave vector), and polarization, as well as the detection methods, different quantum transition pathways can be selectively generated and measured.8 

The spectroscopic observables are determined by the third-order optical response function of the molecular system that includes all distinctive quantum transition sequences consistent with the incident radiations and the detection method. The response function naturally emerges from the quantum mechanical time-dependent perturbation theory treatment of the molecular system in the presence of the three perturbative light-matter interactions of the preparation step. In this section, we sketch the theoretical procedure and present key results relevant to 2D vibrational spectroscopy.

In 2D IR spectroscopy, the molecular system interacts with the incident electric field and, in the electric dipole approximation, the interaction Hamiltonian can be written as
(2)
where μ^ is the electric dipole operator and E(r, t) is the superposition of the electric fields of the three incident (IR, THz, visible, or X-ray) pulses denoted as E1, E2, and E3. The total Hamiltonian of the system is the sum of the system Hamiltonian H0 in the absence of radiation and Hint(t). In the cases of the other 2D vibrational spectroscopy utilizing visible, UV, etc., one needs to consider the corresponding effective field-matter interaction Hamiltonian, e.g., α^:E(r,t)E(r,t) for Raman spectroscopy.
The system evolves in time according to the quantum Liouville equation for the density operator ρ(t) of the system as follows:
(3)
The solution of this equation provides quantitative information about any physical observable of the system A(t) through the expectation value TrÂρ(t), where Tr denotes the trace of a matrix and  is the operator for observable A. A diagonal element ρaa of the density matrix in a basis set a represents the probability that the system is in state a or the population of the system in state a. The off-diagonal element ρab of the density matrix, which is related to coherence or super-position state evolution of two states a and b, gives rise to the temporal oscillation of the aforementioned probability with a frequency ωωab ≡ (EaEb)/ determined by the energy difference of the two states.
Treating Hint(t) as the perturbation to the reference Hamiltonian H0, Eq. (3) can be solved by applying the time-dependent perturbation theory. The solution is expressed as a power series expansion of ρ(t), the zeroth-order term of which is the equilibrium density operator for the unperturbed system ρ(0)(t) = ρeq. Each of the higher-order terms ρ(n)(t) contains n factors of Hint(t) and is given by
(4)
where G0(t) = exp(−iL0t/) is the time evolution operator in the absence of radiation and the Liouville operators are defined as LaA = [Ha, A] for a = 0 or int. According to Eq. (4), the system initially defined by ρ(t0) evolves freely without perturbation for τ1t0, as given by G0(τ1t0), and then interacts with the radiation at time τ1, as given by Lint(τ1). This sequence is repeated n times until the final field-matter interaction at τn, as given by Lint(τn). Finally, the system evolves freely until the observation time t according to G0(tτn). The multiple integrals over τ1, …, τn account for all possible interaction times under the time ordering condition t0τ1 ≤ … ≤ τnt.
Each term of the power series expansion of ρ(t) in Eq. (4) gives rise to the corresponding polarization P(n)(r,t)=Trμ^ρ(n)(t) in the system as follows:
(5)
in terms of the nth order optical response function given by
(6)
where μ(t)=exp(iH0t/)μ^exp(iH0t/) is the dipole operator in the interaction picture and the angular bracket denotes the trace of a matrix. We obtain the linear response by setting n = 1 in Eqs. (5) and (6). The signal of 2D IR spectroscopy is determined by the third order polarization P(3)(r, t) and the third-order response function R(3)(t3, …, t1), the latter being the fourth rank tensor. Note that the time variables t1, …, tn−1 in Eqs. (5) and (6) are time intervals between consecutive field-matter interactions related to τ1, …, τn in Eq. (4) as tm = τm+1τm (1 ≤ mn − 1), while tn = tτn is the time elapsed since the last field-matter interaction. Therefore, t1, …, tn are all positive, and the response function must vanish if any of its time arguments are negative in accordance with the causality principle, as imposed by the Heaviside step functions θ(t) in Eq. (6). In addition, R(n) is a real function because P(n)(r, t) and E(r, t) in Eq. (5) are both real quantities, although individual terms comprising R(n) are complex in general and represent different quantum transition pathways.
The signal electric field Es(n)(r,t) detected in nonlinear spectroscopy is obtained by solving Maxwell’s equation taking the nonlinear polarization P(n)(r, t) as the source. After making simplifying assumptions that (i) the signal field is only weakly absorbed by the medium, (ii) the envelopes of polarization and signal fields vary slowly in time compared to the optical period, (iii) the signal field envelope spatially varies slowly compared to its wave length, (iv) the frequency dispersion of the medium refractive index is weak, the approximate solution can be obtained as30,31
(7)
Here, n(ω) is the refractive index of the medium and Ps(n)(t) is the polarization component propagating with wave vector ks and frequency ωs that are one of the combinations ±k1 ± k2 ± ⋯ ±kn and ±ω1 ± ω2 ± ⋯ ±ωn, respectively. Note that Eq. (7) gives the approximate signal field arising from a single Fourier component of the nth order polarization expanded as8,31
(8)
By changing the location of the detector appropriately, individual components of the polarization with different ks can be selectively measured. Note that the assumption (ii) above could become invalid for the far-IR and THz spectroscopy. See Ref. 32 and references therein for more general approaches to the nonlinear signal field calculation with wider applicability.

2D vibrational spectroscopy usually induces transitions up to the second vibrational excited state. Therefore, a three-level system with eigenstates g, e, and f   is a useful model for the response function relevant to 2D vibrational spectroscopy. Because the third-order response function vanishes for a harmonic oscillator, the model system must represent an anharmonic oscillator where the fundamental transition frequency ωeg is slightly larger than ωfe.

The evaluation of a realistic response function critically depends on the accurate description of the system-bath interaction that is responsible for important spectroscopic phenomena such as dephasing, relaxation, reorientation, spectral diffusion, and population and coherence transfers. Methods to incorporate the effect of environment as well as the multimode vibrational coupling are discussed in Secs. III and IV. In this section, to highlight the structure of the response function, we consider a simple model where a single three-level chromophore interacts with the environment according to the following Hamiltonian:
(9)
Here, ℏωm is the energy of state m in the absence of bath, Vm(q) is the chromophore-bath interaction energy of the state m that depends on the bath degrees of freedom q, HB(q) is the energy of the bath, and the basis states m (m = g, e, f) are chosen as eigenstates of an isolated chromophore. Note that the off-diagonal elements of the chromophore-bath interaction are assumed negligible for the sake of simplicity.33 Using this Hamiltonian, the three nested commutators in the response function in Eq. (6) can be expanded as the sum of eight terms31,34
(10)
where the components Ri(t3, t2, t1) are given by8,35
(11)
assuming that the system is initially in the ground state g. Here, μab is the transition dipole between states a and b obtained by the repeated insertion of the closure relation mmm=1 in Eq. (6), which is assumed to be independent of the bath coordinate (Condon approximation), ω¯ab=ωaωb+Va(q)Vb(q)B is the energy gap averaged over bath degrees of freedom, and Fngabc(t3,t2,t1) is the line shape function expressed in terms of time-ordered exponentials of the fluctuations in the system-bath interactions, Um(q)=Vm(q)Vm(q)B. Therefore, the response function is composed of multiple quantum transition pathways represented by individual Ri, each of which is the product of three factors determining the transition strength (products of transition moments), the transition frequency (coherence oscillation), and the line shape (F1−4). Throughout this paper, third order double-quantum or zero-quantum spectroscopies are not taken into consideration.
To facilitate computation of Fngabc(t3,t2,t1), we can approximately replace the time-ordered exponential operators with normal exponential functions containing difference potential energies Uab(q) = Ua(q) − Ub(q). Alternatively, we can invoke the second-order cumulant expansion approximation, which becomes exact when the fluctuation of the energy gap ω¯ab obeys the Gaussian statistics, to obtain8 
(12)
where the auxiliary functions are given by
(13)
The second terms of R3 and R4 in Eq. (11), that represent coherence evolution during t2, are not included in Eq. (12) and the energy fluctuation between states g and f is assumed to be twice that between g and e, i.e., Ufg(t) ≅ 2Ueg(t). A few approximate expressions for the line-broadening function g(t) can be found in other review articles and books. The six response functions in Eq. (12) can be illustrated using double-sided Feynman diagrams (Fig. 1), which highlight their physical interpretation in terms of ground state bleach, stimulated emission, and excited state absorption contributions. These are further classified as rephasing (top row in Fig. 1) and non-rephasing (bottom row in Fig. 1) versions depending on whether the coherence oscillations during t1 and t3 have the opposite or the same signs, respectively.
FIG. 1.

Double-sided Feynman diagrams illustrating the six contributions of Eq. (12). Time is running from bottom to the top with the vertical lines illustrating the ket and the bra wave functions. Interactions with applied laser fields are illustrated with full arrows, while emitted light is illustrated with dashed arrows. The numbers (zero, one, and two) indicate the number of vibrations excited. The response functions R2A and R4 are identified as ground state bleach as the system is in the ground state during the waiting time t2, while R3 and R1A are stimulated emission contributions, and R1B and R2B are excited state absorption contributions.

FIG. 1.

Double-sided Feynman diagrams illustrating the six contributions of Eq. (12). Time is running from bottom to the top with the vertical lines illustrating the ket and the bra wave functions. Interactions with applied laser fields are illustrated with full arrows, while emitted light is illustrated with dashed arrows. The numbers (zero, one, and two) indicate the number of vibrations excited. The response functions R2A and R4 are identified as ground state bleach as the system is in the ground state during the waiting time t2, while R3 and R1A are stimulated emission contributions, and R1B and R2B are excited state absorption contributions.

Close modal

This general formulation for three-level systems is very powerful and can be used to understand the effect of dynamics on the two-dimensional lineshapes,36 which, for example, allows the distinction between homogeneous and inhomogeneous line-broadenings and defines rules for extracting the correlation functions describing the bath dynamics in the limit of the Gaussian dynamics, for example, using the nodal or center line slope.6,37 This is illustrated in Fig. 2.

FIG. 2.

Illustration of the 2D IR spectra with slow (left) and fast (right) dynamics. The corresponding linear spectra, which have little sensitivity to dynamics, are shown on the top. The diagonal peaks shown with red contours are bleaching signals (R1A, R2A, R3, and R4), while the peaks with blue contours, which are shifted by the anharmonicity below the diagonal, arise from the absorption of the excited state (R1B and R2B). When the bath dynamics scramble the vibrational frequencies faster than the waiting time, all correlation between excitation and detection frequencies is lost and round peaks are obtained, while elongated peaks arise from residual memory.

FIG. 2.

Illustration of the 2D IR spectra with slow (left) and fast (right) dynamics. The corresponding linear spectra, which have little sensitivity to dynamics, are shown on the top. The diagonal peaks shown with red contours are bleaching signals (R1A, R2A, R3, and R4), while the peaks with blue contours, which are shifted by the anharmonicity below the diagonal, arise from the absorption of the excited state (R1B and R2B). When the bath dynamics scramble the vibrational frequencies faster than the waiting time, all correlation between excitation and detection frequencies is lost and round peaks are obtained, while elongated peaks arise from residual memory.

Close modal

In Sec. II, we derived the quantum mechanical linear and nonlinear response functions. Although these equations are formally exact, fully quantum mechanical simulations are still impractical for systems with many degrees of freedom even with state-of-the-art supercomputers. Thus, applicable and efficient methods are required for the calculation of nonlinear response functions for such systems. In this section, we first summarize the classical mechanical approaches to the calculation of the linear and nonlinear response functions, which can be expressed as (multi-)time correlation functions on equilibrium and nonequilibrium trajectories, and then present some results of the 2D IR and pump-probe spectra obtained from the approaches.

The classical mechanical response functions can be derived by using the relationship between the commutator and the Poisson bracket,38,
(14)
Here, X and Y are physical variables and the Poisson bracket is defined as follows:
(15)
When Y is the equilibrium distribution function of the system ρeq = eβH/Z, where Z is the canonical partition function, Eq. (15) becomes {X,ρeq}PB=βẊρeq. Here, β is the reciprocal temperature of the system and Ẋ is the time derivative of X. As a result, we obtain the following well-known general expression for the classical linear response function of a physical quantity A to a perturbation B7 by using Eq. (14) in the quantum linear response function R(1)(t)=iθ(t)TrA(t)[B(0),ρeq], which corresponds to n = 1 in Eq. (6),
(16)
Here, the angular bracket indicates ensemble average. For example, the classical expression of the Raman response function is given as39 
(17)
where α is the total polarizability of the system. Equation (17) shows that the response function for the optical Kerr effect is expressed as the time derivative of time correlation of the polarizability along the equilibrium classical trajectory.
By using Eq. (14), we can also derive any classical nonlinear response function. For example, the fifth-order nonlinear Raman spectroscopy,19,38 which is the second-order response function of the polarizability, is given as39–45 
(18)
These equations have been used for the analyses of the fifth-order nonlinear Raman spectra of liquids.39–41,46
Similarly, the classical third-order response function for 2D IR spectroscopy is expressed as47–49 
(19)
where μ denotes the dipole moment of the system.
It should be noted here that the Poisson brackets of physical variables at different times, e.g., {μ(t),μ(t)}PB, are required for the calculation of nonlinear response functions based on the equilibrium molecular dynamics approach. The explicit expression of {μ(t),μ(t)}PB is
(20)
where the dipole moment is assumed to be a function of particle positions only and we utilized the invariance of Poisson brackets under canonical transformation43,50 (of which the Newtonian time evolution is an example) and the chain rule of the form μ(t)/pα(t)=βμ(t)/qβ(t)qβ(t)/pα(t). The partial derivative matrix, q(t′)/p(t), represents the variation of the position at t′ induced by a small change of momentum at t. This matrix is a sub-matrix of the so-called stability matrix, J(t) = γ(t)/γ(0), where γ(t) = (q(t), p(t)). The time-evolution of the stability matrix is expressed as38,39
(21)
with the initial condition of J(0) = 1. The stability matrix represents the transformation of the phase space along the trajectory and plays an important role in nonlinear mechanics and semiclassical theory.51–53 The nonlinear response functions are highly sensitive to the trajectory, and thus, they can provide more detailed information on the dynamics than the linear response function.39 

Jeon and Cho have investigated the 2D IR spectra of intramolecular vibrations in water.49,54 In their studies, the quantum mechanical/molecular mechanical (QM/MM) approach has been employed for the accurate description of the intramolecular vibrations of solute molecules: deuterated N-methylacetamide (d-NMA) and HOD molecules have been described with the semiempirical PM3 method and scc-DFTB potential, respectively. The flexible TIP3P and SPC/Fw models have been used for solvent water molecules, respectively. Figure 3 shows the 2D IR spectra of the OD stretch of the HOD molecule in a water cluster. In these calculations,54 Jeon and Cho have improved the computational efficiency by exploiting the time reversibility of trajectory and have successfully calculated the 2D IR spectra. The positive and negative peaks are found at (∼2650 cm−1, ∼2700 cm−1) and (∼2600 cm−1, ∼2550 cm−1), respectively. The positive peak corresponds to the stimulated emission and the ground state bleaching of the amide I band, whereas the negative peak is related to the excited state absorption. Note that various general features of the experimental 2D IR spectra are reproduced in the calculated 2D spectra based on classical mechanics: the vertical splitting of positive and negative peaks due to the vibrational anharmonicity, the diagonal elongation of the signal reflecting the inhomogeneity of the solvation environment, the change in lineshape with waiting time due to the spectral diffusion, and the decrease in the tilt angle of the nodal line. The classical 2D spectra also show that the main decay component of the nodal line with a ∼1.6 ps time scale obtained from the calculated 2D spectra is comparable to experimental results.55 

FIG. 3.

Absorptive 2D IR spectra of the OD stretch of the hydrated HOD molecule at different waiting times T. All figures are drawn on the same amplitude scale. Diminishing tilt angles of the nodal line and the signal amplitudes can be noticed with increasing waiting time T, indicating the loss of frequency memory and the vibrational population relaxation, respectively. Reprinted with permission from J. Jeon and M. Cho, J. Phys. Chem. B 118, 8148 (2014). Copyright 2014 American Chemical Society.

FIG. 3.

Absorptive 2D IR spectra of the OD stretch of the hydrated HOD molecule at different waiting times T. All figures are drawn on the same amplitude scale. Diminishing tilt angles of the nodal line and the signal amplitudes can be noticed with increasing waiting time T, indicating the loss of frequency memory and the vibrational population relaxation, respectively. Reprinted with permission from J. Jeon and M. Cho, J. Phys. Chem. B 118, 8148 (2014). Copyright 2014 American Chemical Society.

Close modal

1. Nonequilibrium finite-field method

In Subsection III A, we described the prescription for the linear and nonlinear response functions based on the equilibrium molecular dynamics approach. Since the stability matrix is required in the approach, high computational costs and numerical instability would make it difficult to calculate nonlinear response functions, when many degrees of freedom have to be considered explicitly, for example, the low frequency intermolecular motions and the OH stretches in liquid H2O. Thus, other approaches circumventing the use of the stability matrix would be required for such systems. In this subsection, we introduce the nonequilibrium molecular dynamics approach in which the response functions are evaluated by considering (sequential) external fields, as in real experiments.

By introducing a Liouville operator B and a dimensionless parameter ε, we have the following identity:56 
(22)
By applying Eq. (22) to the linear response function R(1)(t) (for t > 0) of a physical variable A to a perturbation B, we have56–58 
(23)
Here, the first term is the expectation value of A(t) on the perturbed trajectory given by H0εBδ(t), whereas the second term is the expectation value of A on the trajectory expressed by the unperturbed Hamiltonian H0.
In the nonequilibrium finite-field method, the second-order response function of 2D Raman spectroscopy, corresponding to Eq. (18) in the equilibrium molecular dynamics approach, is46 
(24)
Jansen et al. have applied the nonequilibrium finite-field method to the 1D and 2D Raman spectra of liquids.57,58
Based on the nonequilibrium finite-field method, the third-order response function for 2D IR spectroscopy is expressed as48 
(25)
Later, Jansen et al., have developed the efficient method to subtract higher-order responses from calculated signals by combining the nonequilibrium molecular dynamics simulations with the positive and negative electric fields.59 By using this method, the third-order nonlinear response function can be recast in the following form:
(26)
where the bar on the subscript denotes the application of the external field with the opposite direction.

2. Hybrid approach combining equilibrium and nonequilibrium finite-field methods

The nonequilibrium finite-field method does not need the stability matrix for the calculation of nonlinear responses. However, the approach still demands high computational costs; i.e., extra three and seven nonequilibrium trajectories as well as the equilibrium trajectory are required for the second- and third-order nonlinear response functions, respectively.

In order to reduce the number of trajectories, Hasegawa and Tanimura have developed an efficient computational method for higher-order response functions, i.e., the hybrid method combining the equilibrium molecular dynamics and nonequilibrium finite-field methods.60 In the hybrid method, the second-order response function for the 2D Raman spectra given by Eqs. (18) and (24) can be written as61 
(27)
In Eq. (27), the inverse force method is also used, as indicated by the overbar in the last subscript. Equation (27) shows that, in the hybrid method, the second-order response function can be calculated from the time correlation function of Π̇(0) on the equilibrium trajectory and Π(t1 + t2) on the nonequilibrium trajectory generated by the external fields applied at t1.
The third-order response function of the 2D IR spectra can be calculated with the following four terms:61 
(28)
where the inverse forces are denoted with overbars.
Later, Hasegawa and Tanimura47 have proposed another efficient method to calculate the third-order response function by exploiting backward nonequilibrium trajectories. In this method, the third-order response function is re-expressed as
(29)
where μ(t2 + t3)E(t2) and μ̇(t1)E(0) denote the dipole moment at t2 + t3 on the perturbed trajectory with the Hamiltonian H0εμEδ(tt2) and the time derivative of dipole moment at −t1 on the backward perturbed trajectory H0εμEδ(t), respectively, while μ̇(t1)μ̇(0) is evaluated from the equilibrium trajectory. In practical computations, the second term in Eq. (29) is required due to the breakdown of equipartition.61 The backward-forward sampling method given as Eq. (29) is particularly efficient for the sampling of the response function at non-zero t2 values because the third-order response functions at different t2 values can be efficiently calculated by choosing the pair of initial configurations along the forward and backward trajectories.

It should be noted that the experimental third-order nonlinear spectra are obtained under different phase matching conditions [see Eq. (8)]. For example, the 2D IR spectra consist of the responses emitted along the wave vector directions kI = −k1 + k2 + k3 and kII = +k1k2 + k3. However, the calculated third-order response function consists of all the wave vector components generated by three incident electric fields, i.e., kIII = +k1 + k2k3 and kIV = +k1 + k2 + k3 as well as kI and kII. Therefore, the components with kIII and kIV have to be eliminated from the calculated third-order responses for the 2D IR spectrum. The responses arising from the individual phase matching conditions can easily be obtained by using the nonequilibrium molecular dynamics approach in which external electric fields are explicitly considered. Thus, various third-order nonlinear spectroscopies, such as the three-pulse photon echo peak shift and the pump-probe spectrum, can also be calculated.61 In the earlier studies, kIII and kIV components have been eliminated by exploiting the fact that the three-dimensional Fourier transformed spectra of kIII and kIV components show higher frequency oscillation in ω2 than those of kI and kII.47–49,54

Hasegawa and Tanimura have applied the backward-forward sampling method to the 2D IR spectra of liquid HF.47 As shown in Fig. 4, the positive and negative peaks of intermolecular libration motion are located at (550 cm−1, 550 cm−1) and (550 cm−1, 400 cm−1), respectively. The frequency correlation of the libration motion of liquid HF was found to be lost with a time scale of ∼200 fs. It was also found that the width of the anti-diagonal line of the peak is narrow compared with that of liquid water, indicating the suppressed homogeneous broadening in liquid HF.

FIG. 4.

Correlation spectra of liquid HF at t2 = (a) 0 fs, (b) 20 fs, (c) 100 fs, (d) 200 fs, (e) 300 fs, and (f) 400 fs. The diagonal and off-diagonal peaks in each spectrum arise from fundamental and anharmonic oscillations, respectively, of the hydrogen bond network in liquid HF. Reprinted with permission from T. Hasegawa and Y. Tanimura, J. Chem. Phys. 128, 064511 (2008). Copyright 2008 AIP Publishing LLC.

FIG. 4.

Correlation spectra of liquid HF at t2 = (a) 0 fs, (b) 20 fs, (c) 100 fs, (d) 200 fs, (e) 300 fs, and (f) 400 fs. The diagonal and off-diagonal peaks in each spectrum arise from fundamental and anharmonic oscillations, respectively, of the hydrogen bond network in liquid HF. Reprinted with permission from T. Hasegawa and Y. Tanimura, J. Chem. Phys. 128, 064511 (2008). Copyright 2008 AIP Publishing LLC.

Close modal

Yagasaki and Saito investigated the fluctuation and relaxation of the intermolecular motions in liquid water by using the hybrid method with the SPC/E model for water molecules.48  Figure 5 shows the 2D IR spectra of the intermolecular motions of liquid water at several waiting times. The intermolecular translational and librational motions are found below and above 300 cm−1, respectively. As shown in Fig. 5, the tilt angle of the nodal line between the positive and negative peaks of the librational motion decays with a time constant of ∼110 fs. Yagasaki and Saito found that the frequency fluctuation of librational motion becomes three times slower when the translational motions of individual molecules are frozen, indicating that the ultrafast frequency fluctuation of the librational motion arises from the coupling between the translational and librational motions. It was found that the three-pulse stimulated photon-echo signal of the librational motion in water decays with time scales of ∼20 and ∼100 fs, in which the slower time scale is similar to that obtained from the change of the nodal line tilt angle.48 

FIG. 5.

2D IR correlation spectra of intermolecular motion of liquid water at several waiting times. The tilt angle is defined as the angle between the ν1 axis and the nodal line of the libration peak as shown in (a) and is found to decay with a time scale of 115 fs. Reprinted with permission from T. Yagasaki and S. Saito, J. Chem. Phys. 128, 154521 (2008). Copyright 2008 AIP Publishing LLC.

FIG. 5.

2D IR correlation spectra of intermolecular motion of liquid water at several waiting times. The tilt angle is defined as the angle between the ν1 axis and the nodal line of the libration peak as shown in (a) and is found to decay with a time scale of 115 fs. Reprinted with permission from T. Yagasaki and S. Saito, J. Chem. Phys. 128, 154521 (2008). Copyright 2008 AIP Publishing LLC.

Close modal

The energy relaxation process of the librational modes has been examined by the pump-probe signals calculated from the third-order response functions.61 It was found that the absorption changes at 700 and 800 cm−1 are induced by the initial decrease of absorption due to the ground state bleach and stimulated emission followed by the fast energy relaxation of the librational motion to low frequency modes within ∼60 fs and the subsequent slow thermalization that corresponds to its energy relaxation to hot ground states occurs in 500 fs. These time scales of the two relaxation processes are consistent with the experimental results.62–64 

As mentioned above, once third-order response functions with individual phase-matching conditions are obtained, any third-order nonlinear spectrum can be calculated in principle. However, despite these developments of efficient computational approaches, the calculation of the third-order response function is still expensive and computationally demanding. In this regard, it should be noted that other efficient computational methods to estimate mode-to-mode vibrational energy relaxation rates have been developed by Yagasaki and Saito65,66 and Jeon and Cho.67,68

Recently, the nonequilibrium molecular dynamics method has also been applied to the intramolecular vibrational modes of H2O in water.69 In this study, the 2D IR spectra of the OH stretch and the HOH bending modes were calculated with the TTM3-F interaction potential, which is the ab initio-based charge transferable, flexible, and polarizable Thole-type model to describe the inter- and intra-molecular interactions in liquid water. The waiting time-dependent changes of the OH stretch spectra indicate the heterogeneous dynamics of the local H-bonding network which has also been found in an experimental result.70 The spectrum loses diagonal correlation at waiting times larger than 200 fs, which is consistent with the other quantum mechanical calculation results.71,72 The 2D IR spectra of the HOH bend are found to have a smaller nodal line tilt angle at T = 0 fs compared to the OH stretch spectra and become parallel to the ω1 axis within ∼400 fs.

The pump-probe spectra of the OH stretch and the HOH bend in water have also been calculated with the TTM3-F interaction potential.69 The positive and negative peaks of the pump-probe spectra of the OH stretch at (ωpump, ωprobe) = (3600 cm−1, 3225 cm−1) and (3600 cm−1, 3600 cm−1), corresponding to the 1 → 2 and 0 → 1 transitions, show the initial decay with a time constant of 240 fs followed by the slow relaxation to the hot-ground state. The calculated time scale of the initial decay is in agreement with experimental results of 200-270 fs.73–75 The pump-probe spectra of the HOH bend at (1650 cm−1, 1475 cm−1) and (1650 cm−1, 1650 cm−1) decay with a time constant of 250 fs, which is in good agreement with the experimental results of 170–260 fs.62,64 Furthermore, the frequency-resolved kinetic energy analysis revealed the detailed relaxation processes after the excitations of the OH stretch and HOH bend of H2O in water.

In this section, the classical simulation methods for the nonlinear response functions were presented with their applications to the calculations of the 2D IR and pump-probe spectra of liquids. Before closing this section, we summarize the limitations of the classical treatment of nonlinear spectral simulations.

The validity of the classical 2D IR spectra has been investigated by several groups.76–81 It should be noted that the classical nonlinear response functions are not stable, for example, for integrable systems and systems without dissipation.77,78 Sakurai and Tanimura examined the quantum effects on the simulated IR and 2D IR spectra of a Morse oscillator in a harmonic bath by solving the quantum and classical hierarchical equations of motion (HEOM).80 They found that the classical 2D IR spectra are a good approximation of the quantum 2D IR spectra when the system (vibration) is largely modulated by the bath, i.e., via a strong system-bath coupling or a fast bath modulation even in a weak system-bath coupling. On the other hand, there are significant differences between the quantum and classical mechanically calculated spectra when the modulation due to the system-bath interactions is weak or slow. Very recently, Reppert and Brumer have also shown that the classical 2D IR spectra of a Morse oscillator, which mimics the amide I vibration, can reproduce the qualitative features of the quantum 2D IR spectra very well.81 As shown in these studies, the validity of the classical mechanically calculated spectra depends on the system-bath coupling. Notably, in their study, it was shown that the anharmonicity of the classical signal is determined almost entirely by the frequency resolution afforded by the simulated scan time t3.

Classical nonlinear spectral simulations can reproduce spectra calculated with quantum methods qualitatively well when the system vibration is largely modulated by the bath, such as those examples discussed in this section. It should be noted, however, that any vibrational spectra are described as transitions from one vibrational level to another in quantum mechanics, whereas those are obtained from the fluctuation of dipole moment in classical mechanics. In quantum mechanics, wave functions or density matrices are determined by non-local information; i.e., the kinetic energy operator is expressed as the second derivative of coordinates in the coordinate representation. The anharmonic shift in quantum 2D IR spectra is thus determined by information about the whole potential energy profile. On the other hand, no discrete vibrational levels are considered in any classical approaches. In classical mechanics, a trajectory is determined by local information on coordinates and momenta. Consequently, the anharmonic shift corresponding to the frequency difference between the positive and negative peaks in the classical 2D IR spectra arises from the difference between the curvatures of trajectories perturbed by one and two electric fields. As a result, in many examples including those shown in this section, the anharmonic shifts in the classical 2D IR spectra tend to be smaller than those in the experimental or quantum mechanically calculated 2D IR spectra. In addition to the potential profile, the transition dipole moment is also important quantity for quantitative simulations of the IR and 2D IR spectra. In quantum mechanics, a spectral intensity is related to the transition dipole moments between discrete vibrational levels. On the other hand, the classical spectra are obtained by the dipole moment induced by vibrational and conformational changes. Accurate model potential and dipole moment are therefore extremely important to accurately simulate the 2D IR spectra of complicated molecular systems.

A fully quantum mechanical simulation of nuclear degrees of freedom including all the effects of surrounding thermal bath remains unpromising, despite the technological advancement of computers. Therefore, semiclassical approaches, such as (linearized) semiclassical-initial value representation,82 centroid molecular dynamics,83 and ring-polymer molecular dynamics,84 have been developed to accurately calculate the static and dynamic properties of molecules in condensed phases; for example, the IR spectrum of liquid water.85 A formalism for two-time correlation functions has also been developed.86 It is hoped that these semiclassical approaches will be able to provide a novel framework to calculate the 2D IR spectra incorporating the nuclear quantum effects by developing the formalism of nonlinear three-time correlation functions.87 

The response functions described in Sec. II can be calculated with quantum-classical simulations. This essentially involves solving the time-dependent Schrödinger equation. In the first application of the method to systems involving more than one vibration, it was referred to as the Numerical Integration of the Schrödinger Equation (NISE) approach.88–90 To evaluate the response functions, one must identify the essential coordinates directly coupled with light and treat those as weakly anharmonic oscillators represented by coupled three-level systems. All other coordinates are treated classically and only included through their time-dependent modulation of the parameters in the quantum system. In reality, the full vibrational quantum Hamiltonian will contain anharmonicities of all orders. For efficient calculations, one typically employs an effective model Hamiltonian, where all anharmonicities are collected in one effective quartic anharmonicity term sufficient to describe the energy fluctuations of the three energy levels crucial for the 2D IR spectroscopy. This approximation leads to an effective expression without anharmonicity terms that can cause relaxation between the three vibrational levels. The quantum system is, thus, typically described by the time-dependent effective model Hamiltonian
(30)
Here, Bn and Bm are the Bosonic creation and annihilation operators for the N vibrations considered quantum mechanically, which are numbered with indices n and m, respectively. The individual local vibrations are characterized by their frequencies, ωn(t), transition dipoles, μn(t), transition polarizabilities, α¯¯n(t), and anharmonicities, Δn(t). The different local vibrations are mixed by their mutual coupling, Jnm(t). This is a generalization of Eq. (9) applicable to coupled multi-chromophore systems. The time-dependence arises from the couplings of the quantum oscillators with the classical bath coordinates. The last two terms account for the interaction with the applied electric field(s), E(t), and are included through the time-dependent perturbation theory that results in the response function formulation. Determining the fluctuating quantities in the Hamiltonian depends on the system under consideration and will be discussed in Sec. IV A. If we assume that these are known, the response functions can be calculated through the calculation of the time-evolution operators. To determine these, the solution of the time-dependent Schrödinger equation is needed. However, if time is divided into short enough intervals (Δt) that the Hamiltonian can be considered constant the trivial solution of the time-evolution during one time interval is given by the solution of the time-dependent Schrödinger equation with a time-independent Hamiltonian and the time-evolution operators for longer time intervals are obtained by successive application of time-evolution operators for neighboring time intervals.

In practical calculations, the Hamiltonian in Eq. (30) only includes a change in the number of vibrations in the terms involving the coupling with the external electric field(s) and thermalization is not included. Therefore, the remaining part of the Hamiltonian is block diagonal. This allows treating the ground state, single excited states, and double excited states, separately. As there is only one ground state and the energy of this is set to zero, the time-evolution operator is the unit operator. When N vibrational degrees of freedom are treated, there are N single excited states and N(N + 1)/2 double excited states. The time-evolution in each excitation manifold can be evaluated independently in the corresponding harmonic basis, which, in principle, requires the evaluation of matrix exponentials of the corresponding dimensionalities, which is typically done by diagonalizing the Hamiltonian, multiplying the eigenvalues with −iΔt/, taking the exponent and transforming back to the original basis. However, in practice, the time-evolution matrix for the double excited states can be made more efficient as the time required for the direct evaluation of the time-evolution operator scales as N,6 which is a costly procedure to repeat even for moderately sized systems. Among the solutions to this problem are the nonlinear exciton propagation scheme,91 which utilizes an anharmonic perturbation of the harmonic solution, and an approach based on the Trotter algorithm making use of the sparse nature of the Hamiltonian of the doubly excited states.92,93 The latter allowed the calculation of the 2D IR spectra of systems up to 864 coupled vibrations for water ices94,95 as well as applications to the amide I band of full proteins.96–98 Numerous applications of the NISE method for 2D vibrational spectroscopy have been implemented when only one vibration is involved.99–102 A few different implementations applicable to coupled systems have been presented.88–93,103

The terms in the Hamiltonian in Eq. (30) have different origins, and their bath dependence can be parameterized in different ways. In the following, the different terms will be considered one by one. First, the vibrational frequency of a given mode depends on the local environment. The well-known Stark shift is a simple way to describe such frequency shift and fluctuation. Local electric fields modulate the vibrational frequency, and this type of effect can be parameterized by ab initio calculations of the vibrational frequency in different point charge environments of inhomogeneous electric fields. In this way, the solvent dependence of the vibrational frequencies of a number of common vibrations has been parameterized expressing the vibrational frequency in an expansion of the electric potential104–108 or electric field102,109–125 and sometimes the electric field gradient126–130 on specific points inside a molecule. The parameterization of the solvatochromic vibrational frequency shift parameters can either be obtained from ab initio calculations or empirical fitting to experimental data. These mappings then allow obtaining a frequency trajectory from a molecular dynamics trajectory. As typical frequencies fluctuate on a sub-picosecond time scale, the trajectories need to be sampled and saved every ∼20 fs and typical lengths are on the order of nanoseconds.

In reality, the solvatochromic vibrational frequency shift is not just an electrostatic effect. Charge transfer, dispersion forces, Pauli repulsion, polarization, and multipole effects may further play a role.121,125,131–137 Empirical models may compensate for this as long as other effects are correlated with electrostatic potential/field fluctuations. Ab initio calculations based on clusters, where the electrons of the solvent molecules are treated explicitly, may also be able to account for such effects in an efficient way. Detailed analysis of the amide I vibrations has shown that the various other effects tend to cancel out.132,136 For other systems, this may not be the case.

The anharmonicity also depends on the solvent, and this dependence has been included in mappings.102,126,128 Still, for many systems such as the amide I vibrations in proteins, the variation in the anharmonicity is small relative to the variation in the frequencies and the solvent dependence can be neglected using a constant anharmonicity value for each type of vibrational mode.

The couplings between different local modes, for long-range interactions, are well approximated by the transition-dipole coupling model.138–143 Transition charge coupling models have been developed for intermediate range couplings.130,144,145 These models account for some multipole contributions to the electrostatic interaction. For vibrational modes directly connected through covalent bonds, ab initio based coupling models have been developed. This has both been done for the OH stretch vibrations in the water molecule,146 for neighboring amide I modes in a peptide chain,117,118,130,145,147,148 and amide I and II modes in the same and neighboring peptide units of proteins.129,149 Ab initio modeling of the vibrational couplings in small peptides and DNA has also been reported based on the polarizable continuum model (PCM) description of solvent.150,151 For the peptide units, this coupling depends on the local configuration, as characterized by the Ramachandran angles. The developed mappings were thus made by varying these angles systematically for small di- or tri-peptides and then calculate the vibrational coupling through ab initio calculations. Just as for the vibrational frequencies, the couplings can be extracted from molecular dynamics trajectories using these mappings.

Figure 6 illustrates a comparison between the spectra simulated for the OH stretch of water molecules diluted in acetonitrile and experimental data. The evolution of cross-peaks between the high-frequency asymmetric OH-stretch and the low-frequency symmetric OH-stretch is well described by the theory. Furthermore, the difference in the dynamics of the two peaks, which was the focus of the study,152 is revealed by the center line slope analysis. From the simulations, this could be interpreted in terms of non-Gaussian dynamics as the bath dynamics for strong hydrogen bond environments is faster than for weak hydrogen bond environments. As strong hydrogen bonded OH-stretches absorb at lower frequencies, the fast dynamics contributes more to the low frequency symmetric vibration peak. While the anharmonicity of OH-stretch modes is around 200 cm−1, a sequential transition is observed with a much smaller apparent anharmonicity. This phenomenon arises as the large anharmonicity largely localizes the double excited states.

FIG. 6.

2D IR spectra for parallel laser polarization of the OH-stretch of water dissolved in acetonitrile at different waiting times t23 reproduced with permission from Jansen et al. J. Phys. Chem. A 113, 6260 (2009), Copyright 2009 American Chemical Society. (a) and (b) present experimental and simulated data, respectively. Dotted lines show the frequencies of symmetric and asymmetric stretching modes of the H2O molecule. The 2D spectra are normalized to the maximal amplitude. Red contours illustrate bleach signals, while the cyan contours illustrate excited state absorption. The contour lines are drawn at 10% steps of the maximal amplitude in each individual plot. Solid black illustrated the max lines.

FIG. 6.

2D IR spectra for parallel laser polarization of the OH-stretch of water dissolved in acetonitrile at different waiting times t23 reproduced with permission from Jansen et al. J. Phys. Chem. A 113, 6260 (2009), Copyright 2009 American Chemical Society. (a) and (b) present experimental and simulated data, respectively. Dotted lines show the frequencies of symmetric and asymmetric stretching modes of the H2O molecule. The 2D spectra are normalized to the maximal amplitude. Red contours illustrate bleach signals, while the cyan contours illustrate excited state absorption. The contour lines are drawn at 10% steps of the maximal amplitude in each individual plot. Solid black illustrated the max lines.

Close modal

A linear absorption spectrum can be calculated from the one-time (two-point) transition-dipole response function, while the 2D IR spectra are governed by a combination of six three-time (four-point) transition-dipole response functions related to rephasing and non-rephasing versions of the ground state bleach, stimulated emission, and excited state absorption pathways (Sec. II). The transition-dipoles of different states are reflected in the relative intensity of the peaks. Ultrafast reorientation of the transition-dipoles resulting from molecular reorientation153–156 and vibrational energy transfer72,153,157 can be monitored by measuring 2D IR spectra using different polarizations158,159 of the applied laser fields. Cross peak intensities in the 2D IR spectra depend on the relative orientation of the transition dipoles of the two involved resonances. This allows for the determination of molecular geometries160–162 and is powerful in the peak assignment. For many vibrations, the transition-dipole magnitude is fairly independent of the solvent environment. A clear exception of this is the OH-stretch vibration, which in water shows a five-fold change across the spectrum,163 with the very weak dipole strength of the free OH-stretch vibration and very strong absorption of a strongly H-bonded OH species. Therefore, it is crucial to include the variation of the transition-dipole in the modeling. This has been done very similar to how the solvent effects are accounted for through mapping for the vibrational frequencies.102,128,130,163

Raman spectra,102,146,164,165 sum-frequency generation,166–169 and 2D sum-frequency generation17,18,168,170–173 spectroscopy signals can be calculated using similar response functions as for linear absorption and 2D IR spectroscopies. The only major change is that one needs to replace transition-dipoles in the expressions with transition-polarizabilities. The possible polarization setups available increase as the transition-polarizabilities represent two interactions with external visible laser fields, while the transition-dipoles represent one interaction with an infrared laser field. This opens a possibility of directly probing molecular orientational distribution. Most importantly in the sum-frequency generation type experiments, which are only sensitive to non-centrosymmetric signal contributions,167,168,174,175 the orientation of the vibrations relative to the sample surface can be determined. Just as the transition-dipole, the transition-polarizability depends on the environment. This has been studied for the OH-stretch of water using an electric field mapping. Still, the solvent effect both for the OH-stretch and other systems is relatively small.102,146 The difference between the transition-dipole and the transition-polarizability solvent effects results in a non-coincidence effect147,176,177 making the infrared and Raman type spectra exhibit significantly different lineshapes. In particular, Raman and SFG spectra in the OH-stretch region are much more sensitive to free OH-species than infrared spectra.146 

The quantum-classical methods discussed here have a number of important advantages and limitations. The methods are powerful in providing a quantum description of a given system and allow for treating essentially arbitrary bath dynamics, as provided by the combination of mappings and molecular dynamics simulations. For example, the methods allow describing non-Gaussian dynamics152,178–182 and chemical exchange101,155,179,183–192 arising when the bath coordinates are not harmonic. Commonly used methods invoking the second-order cumulant approximation34,193 or methods assuming a harmonic bath spectral density194–197 cannot account for such effects. The methods can use input from molecular dynamics simulations and are thus able to predict spectra starting from first principles. Sometimes it is not desirable to perform a full molecular dynamics simulation, and then stochastic models can be employed instead.183,198

The Hamiltonian in Eq. (30) does not allow for relaxations between the different excitation manifolds. This means that vibrational energy relaxation is intrinsically neglected. Furthermore, the methods account for the effect that the bath exerts on the system, while the feedback of the system to the bath when the system is on an excited state is not included.140,199–201 This results in a lack of correct thermalization in the quantum system, which leads to an equal population of all quantum states after equilibration.201 In practice, this means that artifacts can be expected at low temperatures, and when broad spectral ranges are considered, where the thermal energy is comparable to the width of the considered spectral region. In systems like water and ice, the OH-stretch band is broader than the thermal energy at room temperature, but still fairly good approximations of the spectra are found for such systems as long as the waiting time is not too long.72 For longer waiting times, persistent ground state bleach signals are observed in experiments resulting from the thermalization effect. A number of attempts to overcome this issue have been reported;140–143 however, including feedback to the bath and thermalization effects that, generally, require simultaneous propagation of the bath modes and the system Hamiltonian resulting in the need for extensive computational time.

Another limitation of the method is that the modes treated quantum mechanically should be well defined. OH-stretch vibrations are an example of such well-defined modes as the OH-bond under normal circumstances does not break at ambient temperatures. On the contrary, hydrogen bonds form and break all the time and such low-frequency volatile bond vibrations are not suited for this quantum-classical treatment as their nature is dynamically changed on the short time-scale of the spectroscopy. If a vibration of interest is strongly coupled with many other modes, this also poses a problem as all the strongly coupled modes would need to be included in the system Hamiltonian for an accurate description requiring extensive parameterization and long simulation times. In such case, the methods discussed in Sec. III of this paper may be worth considering.

While 2D IR spectra contain valuable information on couplings, anharmonicity, and frequency correlation, more information can be obtained through higher-order methods generating multidimensional spectra. Three-dimensional infrared spectroscopy has been explored experimentally202 and simulated203,204 with a method based on that described in Sec. IV. Still the interpretation of such spectra remains very challenging. Infrared methods can also be combined with electronic spectroscopy205,206 to reveal the vibrations on the electronic state or to utilize the longer lifetime of the electronic state to probe reaction dynamics. With fluorescence detection,207 2D IR spectroscopy and microscopy would also be possible.208 

In the NISE method described in Sec. IV, an essential step is the modeling and parameterization of the multi-chromophore Hamiltonian from electronic structure calculations on their static structures. By contrast, the classical mechanical approach in Sec. III relies on the accurate description of the molecular vibrational property in the context of classical dynamics. This often requires the use of sophisticated molecular mechanics potential models21,209 or quantum mechanical potentials. In this regard, recently proposed ab initio theories of vibrational solvatochromism135–137 and direct QM/MM or full ab initio MD simulations of vibrational spectra210–213 can be fruitfully adapted for either approach to improve their efficiency or accuracy.

2D vibrational spectra can be calculated using the exact hierarchical equations of motion (HEOM) approach.80,214–215 It should be, however, noted that this HEOM approach is based on the assumption of a spectral density description of the bath, which means that the bath is modelled as a collection of independent harmonic oscillators, i.e., the Caldeira-Leggett quantum dissipative bath model.216 Under this assumption, the quantum correlation between the system and the bath can be fully accounted for. However, this method scales rather unfavorably, limiting its application to relatively small systems even though a number of approximations have been developed to improve efficiency.217 Another limitation is that the method requires the parameterization of the spectral density.

In summary, we have presented a brief introduction to the fundamental theory and computational methods used to numerically simulate linear and nonlinear vibrational response functions and the corresponding spectra. One of the widely used methods is to approximate the associated quantum mechanical nonlinear response function with invoking classical approximations so that the linear and nonlinear vibrational response functions can be evaluated by using classical MD simulation trajectories of either equilibrium or non-equilibrium systems. The second method of choice is to solve the time-dependent Schrödinger equation of coupled oscillators where their frequencies, coupling constants, and anharmonicities that are fluctuating in time due to the system-bath interactions can be obtained from independent computational methods. Here, a few representative cases were presented and discussed. Over the past decade, a variety of 2D vibrational spectroscopic techniques that utilize femtosecond IR, THz, and visible pulses have been developed and used to study large-scale delocalized modes in condensed phases and ultrafast reaction dynamics during chemical and biological reactions. Thus, we anticipate that the computational and theoretical methods for accurate calculations of coherent multidimensional vibrational spectroscopic signals of coupled multi-chromophore systems could be incorporated into various MD simulation and quantum chemistry calculation program packages to help the experimentalists in this and related areas to interpret their experimental results and to further understand underlying principles, mechanisms, and functions of materials and molecules in condensed phases.

This work was supported by IBS-R023-D1 (M.C.) and JSPS KAKENHI under Grant No. JP16H02254 (S.S.).

The authors declare no competing financial interests.

1.
P.
Hamm
,
M.
Lim
, and
R. M.
Hochstrasser
,
J. Phys. Chem. B
102
,
6123
(
1998
).
2.
M.
Khalil
and
A.
Tokmakoff
,
Chem. Phys.
266
,
213
(
2001
).
3.
J. C.
Wright
,
Int. Rev. Phys. Chem.
21
,
185
(
2002
).
4.
M.
Khalil
,
N.
Demirdöven
, and
A.
Tokmakoff
,
J. Phys. Chem. A
107
,
5258
(
2003
).
6.
S. T.
Roberts
,
J. J.
Loparo
, and
A.
Tokmakoff
,
J. Chem. Phys.
125
,
084502
(
2006
).
7.
D. A.
McQuarrie
,
Statistical Mechanics
(
Harper & Row
,
New York
,
1976
),
8.
M.
Cho
,
Two-Dimensional Optical Spectroscopy
(
CRC Press
,
Boca Raton
,
2009
).
9.
M.
Cho
, in
Advances in Multi-Photon Processes and Spectroscopy
, edited by
S. H.
Lin
,
A. A.
Villaeys
, and
Y.
Fujimura
(
World Scientific Publishing Co.
,
Singapore
,
1999
), p.
229
.
11.
K.
Okumura
and
Y.
Tanimura
,
J. Chem. Phys.
106
,
1687
(
1997
).
12.
K.
Okumura
and
Y.
Tanimura
,
J. Chem. Phys.
107
,
2267
(
1997
).
13.
K.
Park
and
M.
Cho
,
J. Chem. Phys.
109
,
10559
(
1998
).
14.
J.
Sung
and
M.
Cho
,
J. Chem. Phys.
113
,
7072
(
2000
).
15.
K.
Park
and
M.
Cho
,
J. Chem. Phys.
112
,
5021
(
2000
).
16.
K.
Park
and
M.
Cho
,
J. Chem. Phys.
112
,
10496
(
2000
).
17.
J.
Bredenbeck
,
A.
Ghosh
,
M.
Smits
, and
M.
Bonn
,
J. Am. Chem. Soc.
130
,
2152
(
2008
).
18.
W.
Xiong
,
J. E.
Laaser
,
R. D.
Mehlenbacher
, and
M. T.
Zanni
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
20902
(
2011
).
19.
Y.
Tanimura
and
S.
Mukamel
,
J. Chem. Phys.
99
,
9496
(
1993
).
20.
W.
Kuehn
,
K.
Reimann
,
M.
Woerner
,
T.
Elsaesser
,
R.
Hey
, and
U.
Schade
,
Phys. Rev. Lett.
107
,
067401
(
2011
).
21.
J.
Savolainen
,
S.
Ahmed
, and
P.
Hamm
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
20402
(
2013
).
22.
A.
Shalit
,
S.
Ahmed
,
J.
Savolainen
, and
P.
Hamm
,
Nat. Chem.
9
,
273
(
2016
).
23.
H.
Ito
,
T.
Hasegawa
, and
Y.
Tanimura
,
J. Chem. Phys.
141
,
124503
(
2014
).
24.
H.
Ito
,
J.-Y.
Jo
, and
Y.
Tanimura
,
Struct. Dyn.
2
,
054102
(
2015
).
25.
K.
Park
,
M.
Cho
,
S.
Hahn
, and
D.
Kim
,
J. Chem. Phys.
111
,
4131
(
1999
).
27.
W.
Zhao
and
J. C.
Wright
,
Phys. Rev. Lett.
84
,
1411
(
2000
).
28.
P. M.
Donaldson
,
R.
Guo
,
F.
Fournier
,
E. M.
Gardner
,
L. M. C.
Barter
,
C. J.
Barnett
,
I. R.
Gould
,
D. R.
Klug
,
D. J.
Palmer
, and
K. R.
Willison
,
J. Chem. Phys.
127
,
114513
(
2007
).
29.
F.
Fournier
,
R.
Guo
,
E. M.
Gardner
,
P. M.
Donaldson
,
C.
Loeffeld
,
I. R.
Gould
,
K. R.
Willison
, and
D. R.
Klug
,
Acc. Chem. Res.
42
,
1322
(
2009
).
30.
Y. R.
Shen
,
The Principle of Nonlinear Optics
(
Wiley
,
New York
,
1984
).
31.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
Oxford
,
1995
).
32.
N.
Belabas
and
D. M.
Jonas
,
J. Opt. Soc. Am. B
22
,
655
(
2005
).
33.
J.
Sung
,
R. J.
Silbey
, and
M.
Cho
,
J. Chem. Phys.
115
,
1422
(
2001
).
35.
36.
K.
Kwac
and
M.
Cho
,
J. Phys. Chem. A
107
,
5903
(
2003
).
37.
K.
Lazonder
,
M. S.
Pshenichnikov
, and
D. A.
Wiersma
,
Opt. Lett.
31
,
3354
(
2006
).
38.
S.
Mukamel
,
V.
Khidekel
, and
V.
Chernyak
,
Phys. Rev. E
53
,
R1
(
1996
).
39.
S.
Saito
and
I.
Ohmine
,
J. Chem. Phys.
108
,
240
(
1998
).
40.
A.
Ma
and
R. M.
Stratt
,
Phys. Rev. Lett.
85
,
1004
(
2000
).
41.
S.
Saito
and
I.
Ohmine
,
Phys. Rev. Lett.
88
,
207401
(
2002
).
42.
J.
Cao
,
S.
Yang
, and
J.
Wu
,
J. Chem. Phys.
116
,
3760
(
2002
).
43.
A.
Ma
and
R. M.
Stratt
,
J. Chem. Phys.
116
,
4962
(
2002
).
44.
R. A.
Denny
and
D. R.
Reichman
,
J. Chem. Phys.
116
,
1987
(
2002
).
45.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
J. Chem. Phys.
119
,
6073
(
2003
).
46.
S.
Saito
and
I.
Ohmine
,
J. Chem. Phys.
125
,
084506
(
2006
).
47.
T.
Hasegawa
and
Y.
Tanimura
,
J. Chem. Phys.
128
,
064511
(
2008
).
48.
T.
Yagasaki
and
S.
Saito
,
J. Chem. Phys.
128
,
154521
(
2008
).
49.
J.
Jeon
and
M.
Cho
,
New J. Phys.
12
,
065001
(
2010
).
50.
H.
Goldstein
,
Classical Mechanics
, 2nd ed. (
Addison-Wesley
,
Reading
,
1980
).
51.
A.
Lichtenberg
and
M.
Lieberman
,
Regular and Chaotic Dynamics
(
Springer-Verlag
,
New York
,
1992
).
52.
W. H.
Miller
,
J. Chem. Phys.
53
,
1949
(
1970
).
53.
H.
Nakamura
,
J. Phys. Chem. A
110
,
10929
(
2006
).
54.
J.
Jeon
and
M.
Cho
,
J. Phys. Chem. B
118
,
8148
(
2014
).
55.
J. B.
Asbury
,
T.
Steinel
,
K.
Kwak
,
S. A.
Corcelli
,
C. P.
Lawrence
,
J. L.
Skinner
, and
M. D.
Fayer
,
J. Chem. Phys.
121
,
12431
(
2004
).
56.
S.
Mukamel
and
J. B.
Maddox
,
J. Chem. Phys.
121
,
36
(
2004
).
57.
T. L. C.
Jansen
,
J. G.
Snijders
, and
K.
Duppen
,
J. Chem. Phys.
113
,
307
(
2000
).
58.
T. L. C.
Jansen
,
J. G.
Snijders
, and
K.
Duppen
,
J. Chem. Phys.
114
,
10910
(
2001
).
59.
T. l. C.
Jansen
,
K.
Duppen
, and
J. G.
Snijders
,
Phys. Rev. B
67
,
134206
(
2003
).
60.
T.
Hasegawa
and
Y.
Tanimura
,
J. Chem. Phys.
125
,
074512
(
2006
).
61.
T.
Yagasaki
and
S.
Saito
,
Acc. Chem. Res.
42
,
1250
(
2009
).
62.
S.
Ashihara
,
N.
Huse
,
A.
Espagna
,
E. T. J.
Nibbering
, and
T.
Elsaesser
,
J. Phys. Chem. A
111
,
743
(
2007
).
63.
L.
Chieffo
,
J.
Shattuck
,
J. J.
Amsden
,
S.
Erramilli
, and
L. D.
Ziegler
,
Chem. Phys.
341
,
71
(
2007
).
64.
J.
Lindner
,
D.
Cringus
,
M. S.
Pshenichnikov
, and
P.
Vöhringer
,
Chem. Phys.
341
,
326
(
2007
).
65.
T.
Yagasaki
and
S.
Saito
,
J. Chem. Phys.
134
,
184503
(
2011
).
66.
T.
Yagasaki
and
S.
Saito
,
J. Chem. Phys.
135
,
244511
(
2011
).
67.
J.
Jeon
and
M.
Cho
,
J. Chem. Phys.
135
,
214504
(
2011
).
68.
J.
Jeon
,
J. H.
Lim
,
S.
Kim
,
H.
Kim
, and
M.
Cho
,
J. Phys. Chem. A
119
,
5356
(
2015
).
69.
S.
Imoto
,
S. S.
Xantheas
, and
S.
Saito
,
J. Chem. Phys.
139
,
044503
(
2013
).
70.
J. D.
Eaves
,
J. J.
Loparo
,
C. J.
Fecko
,
S. T.
Roberts
,
A.
Tokmakoff
, and
P. L.
Geissler
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
13019
(
2005
).
71.
A.
Paarmann
,
T.
Hayashi
,
S.
Mukamel
, and
R. J. D.
Miller
,
J. Chem. Phys.
128
,
191103
(
2008
).
72.
T. L. C.
Jansen
,
B. M.
Auer
,
M.
Yang
, and
J. L.
Skinner
,
J. Chem. Phys.
132
,
224503
(
2010
).
73.
A. J.
Lock
and
H. J.
Bakker
,
J. Chem. Phys.
117
,
1708
(
2002
).
74.
D.
Kraemer
,
M. L.
Cowan
,
A.
Paarmann
,
N.
Huse
,
E. T.
Nibbering
,
T.
Elsaesser
, and
R. J.
Miller
,
Proc. Natl. Acad. Sci. U. S. A.
105
,
437
(
2008
).
75.
K.
Ramasesha
,
L. D.
Marco
,
A.
Mandal
, and
A.
Tokmakoff
,
Nat. Chem.
5
,
935
(
2013
).
76.
C.
Dellago
and
S.
Mukamel
,
Phys. Rev. E
67
,
035205
(
2003
).
77.
W. G.
Noid
,
G. S.
Ezra
, and
R. F.
Loring
,
J. Phys. Chem. B
108
,
6536
(
2004
).
78.
M.
Kryvohuz
and
J.
Cao
,
Phys. Rev. Lett.
96
,
030403
(
2006
).
79.
S. V.
Malinin
and
V. Y.
Chernyak
,
Phys. Rev. E
77
,
056201
(
2008
).
80.
A.
Sakurai
and
Y.
Tanimura
,
J. Phys. Chem. A
115
,
4009
(
2011
).
81.
M.
Reppert
and
P.
Brumer
,
J. Chem. Phys.
148
,
064101
(
2018
).
82.
W. H.
Miller
,
J. Phys. Chem. A
105
,
2942
(
2001
).
83.
G. A.
Voth
, “
Path-integral centroid methods in quantum statistical mechanics and dynamics
,” in
New Methods in Computational Quantum Mechanics
,
Advances in Chemical Physics
Vol. 93, edited by
I.
Prigogine
and
S. A.
Rice
(
Wiley, 1996
), pp.
135-218
.
84.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
85.
J.
Liu
,
W. H.
Miller
,
G. S.
Fanourgakis
,
S. S.
Xantheas
,
S.
Imoto
, and
S.
Saito
,
J. Chem. Phys.
135
,
244503
(
2011
).
86.
D. R.
Reichman
,
P.-N.
Roy
,
S.
Jang
, and
G. A.
Voth
,
J. Chem. Phys.
113
,
919
(
2000
).
87.
D. R.
Moberg
,
M.
Alemi
, and
R. F.
Loring
,
J. Chem. Phys.
143
,
084101
(
2015
).
88.
T. l. C.
Jansen
and
J.
Knoester
,
J. Phys. Chem. B
110
,
22910
(
2006
).
90.
M.
Kobus
,
R. D.
Gorbunov
,
P. H.
Nguyen
, and
G.
Stock
,
Chem. Phys.
347
,
208
(
2008
).
91.
C.
Falvo
,
B.
Palmieri
, and
S.
Mukamel
,
J. Chem. Phys.
130
,
184501
(
2009
).
92.
A.
Paarmann
,
T.
Hayashi
,
S.
Mukamel
, and
R. J. D.
Miller
,
J. Chem. Phys.
130
,
204110
(
2009
).
93.
C.
Liang
and
T. L. C.
Jansen
,
J. Chem. Theory Comput.
8
,
1706
(
2012
).
94.
L.
Shi
,
J. L.
Skinner
, and
T. L. C.
Jansen
,
Phys. Chem. Chem. Phys.
18
,
3772
(
2016
).
95.
H.
Tran
,
A. V.
Cunha
,
J. J.
Shephard
,
A.
Shalit
,
P.
Hamm
,
T. L. C.
Jansen
, and
C. G.
Salzmann
,
J. Chem. Phys.
147
,
144501
(
2017
).
96.
C.
Liang
,
M.
Louhivuori
,
S. J.
Marrink
,
T. L. C.
Jansen
, and
J.
Knoester
,
J. Phys. Chem. Lett.
4
,
448
(
2013
).
97.
A. S.
Bondarenko
and
T. L. C.
Jansen
,
J. Chem. Phys.
142
,
212437
(
2015
).
98.
A. V.
Cunha
,
A. S.
Bondarenko
, and
T. L. C.
Jansen
,
J. Chem. Theory Comput.
12
,
3982
(
2016
).
99.
K.
Kwac
and
M.
Cho
,
J. Chem. Phys.
119
,
2256
(
2003
).
100.
T. l. C.
Jansen
,
W.
Zhuang
, and
S.
Mukamel
,
J. Chem. Phys.
121
,
10577
(
2004
).
101.
M. F.
DeCamp
,
L.
DeFlores
,
J. M.
McCracken
,
A.
Tokmakoff
,
K.
Kwac
, and
M.
Cho
,
J. Phys. Chem. B
109
,
11016
(
2005
).
102.
B.
Auer
,
R.
Kumar
,
J. R.
Schmidt
, and
J. L.
Skinner
,
Proc. Natl. Acad. Sci. U. S. A.
104
,
14215
(
2007
).
103.
C. R.
Baiz
,
K. J.
Kubarych
,
E.
Geva
, and
E. L.
Sibert
,
J. Phys. Chem. A
115
,
5354
(
2011
).
104.
S.
Ham
and
M.
Cho
,
J. Chem. Phys.
118
,
6915
(
2003
).
105.
K.
Kwac
and
M.
Cho
,
J. Chem. Phys.
119
,
2247
(
2003
).
106.
J.-H.
Choi
,
S.
Ham
, and
M.
Cho
,
J. Phys. Chem. B
107
,
9132
(
2003
).
107.
J.-H.
Choi
,
K.-I.
Oh
,
H.
Lee
,
C.
Lee
, and
M.
Cho
,
J. Chem. Phys.
128
,
134506
(
2008
).
108.
J.-H.
Choi
and
M.
Cho
,
J. Chem. Phys.
134
,
154513
(
2011
).
109.
S. A.
Corcelli
,
C. P.
Lawrence
, and
J. L.
Skinner
,
J. Chem. Phys.
120
,
8107
(
2004
).
110.
Y. S.
Lin
,
J. M.
Shorb
,
P.
Mukherjee
,
M. T.
Zanni
, and
J. L.
Skinner
,
J. Phys. Chem. B
113
,
592
(
2009
).
111.
L.
Wang
,
C. T.
Middleton
,
M. T.
Zanni
, and
J. L.
Skinner
,
J. Phys. Chem. B
115
,
3713
(
2011
).
112.
M.
Reppert
and
A.
Tokmakoff
,
J. Chem. Phys.
138
,
134116
(
2013
).
113.
M.
Reppert
and
A.
Tokmakoff
,
J. Chem. Phys.
143
,
061102
(
2015
).
114.
T. M.
Watson
and
J. D.
Hirst
,
J. Phys. Chem. A
107
,
6843
(
2003
).
115.
T. M.
Watson
and
J. D.
Hirst
,
Mol. Phys.
103
,
1531
(
2005
).
116.
R. D.
Gorbunov
,
D. S.
Kosov
, and
G.
Stock
,
J. Chem. Phys.
122
,
224904
(
2005
).
117.
R. D.
Gorbunov
and
G.
Stock
,
Chem. Phys. Lett.
437
,
272
(
2007
).
118.
H.
Maekawa
and
N.-H.
Ge
,
J. Phys. Chem. B
114
,
1434
(
2010
).
119.
K. E.
Furse
and
S. A.
Corcelli
,
J. Chem. Theory Comput.
5
,
1959
(
2009
).
120.
C. S.
Miller
and
S. A.
Corcelli
,
J. Phys. Chem. B
114
,
8565
(
2010
).
121.
S.
Li
,
J. R.
Schmidt
,
S. A.
Corcelli
,
C. P.
Lawrence
, and
J. L.
Skinner
,
J. Chem. Phys.
124
,
204110
(
2006
).
122.
N. M.
Levinson
,
E. E.
Bolte
,
C. S.
Miller
,
S. A.
Corcelli
, and
S. G.
Boxer
,
J. Am. Chem. Soc.
133
,
13236
(
2011
).
123.
C.
Liang
,
K.
Kwak
, and
M.
Cho
,
J. Phys. Chem. Lett.
8
,
5779
(
2017
).
124.
K.
Cai
,
C.
Han
, and
J.
Wang
,
Phys. Chem. Chem. Phys.
11
,
9149
(
2009
).
125.
E.
Małolepsza
and
J. E.
Straub
,
J. Phys. Chem. B
118
,
7848
(
2014
).
126.
T.
Hayashi
,
T. l. C.
Jansen
,
W.
Zhuang
, and
S.
Mukamel
,
J. Phys. Chem. A
109
,
64
(
2005
).
127.
T.
Hayashi
,
W.
Zhuang
, and
S.
Mukamel
,
J. Phys. Chem. A
109
,
9747
(
2005
).
128.
T. I. C.
Jansen
and
J.
Knoester
,
J. Chem. Phys.
124
,
044502
(
2006
).
129.
R.
Bloem
,
A. G.
Dijkstra
,
T. l. C.
Jansen
, and
J.
Knoester
,
J. Chem. Phys.
129
,
055101
(
2008
).
130.
S.
Roy
,
J.
Lessing
,
G.
Meisl
,
Z.
Ganim
,
A.
Tokmakoff
,
J.
Knoester
, and
T. l. C.
Jansen
,
J. Chem. Phys.
135
,
234507
(
2011
).
131.
T. L. C.
Jansen
,
J. Phys. Chem. B
118
,
8162
(
2014
).
133.
H.
Torii
,
J. Phys. Chem. B
114
,
13403
(
2010
).
134.
135.
B.
Błasiak
,
H.
Lee
, and
M.
Cho
,
J. Chem. Phys.
139
,
044111
(
2013
).
136.
B.
Błasiak
and
M.
Cho
,
J. Chem. Phys.
140
,
164107
(
2014
).
137.
B.
Błasiak
and
M.
Cho
,
J. Chem. Phys.
143
,
164111
(
2015
).
138.
S.
Krimm
and
J.
Bandekar
,
Adv. Protein Chem.
38
,
181
(
1986
).
139.
H.
Torii
and
M.
Tasumi
,
J. Chem. Phys.
96
,
3379
(
1992
).
140.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
129
,
124505
(
2008
).
141.
G.
Hanna
and
E.
Geva
,
J. Phys. Chem. B
113
,
9278
(
2009
).
142.
K.
Kwac
and
E.
Geva
,
J. Phys. Chem. B
116
,
2856
(
2012
).
143.
C. P.
van der Vegte
,
S.
Knop
,
P.
Vöhringer
,
J.
Knoester
, and
T. L. C.
Jansen
,
J. Phys. Chem. B
118
,
6256
(
2014
).
144.
P.
Hamm
and
S.
Woutersen
,
Bull. Chem. Soc. Jpn.
75
,
985
(
2002
).
145.
T. l. C.
Jansen
,
A. G.
Dijkstra
,
T. M.
Watson
,
J. D.
Hirst
, and
J.
Knoester
,
J. Chem. Phys.
125
,
044312
(
2006
).
146.
B. M.
Auer
and
J. L.
Skinner
,
J. Chem. Phys.
128
,
224511
(
2008
).
148.
S.
Ham
,
S.
Cha
,
J.-H.
Choi
, and
M.
Cho
,
J. Chem. Phys.
119
,
1451
(
2003
).
149.
A. G.
Dijkstra
,
T. L. C.
Jansen
, and
J.
Knoester
,
J. Phys. Chem. A
114
,
7315
(
2010
).
150.
A.
Biancardi
,
C.
Cappelli
,
B.
Mennucci
, and
R.
Cammi
,
J. Phys. Chem. B
114
,
4924
(
2010
).
151.
A.
Biancardi
,
R.
Cammi
,
C.
Cappelli
,
B.
Mennucci
, and
J.
Tomasi
,
Theor. Chem. Acc.
131
,
1157
(
2012
).
152.
T. l. C.
Jansen
,
D.
Cringus
, and
M. S.
Pshenichnikov
,
J. Phys. Chem. A
113
,
6260
(
2009
).
153.
K. J.
Gaffney
,
I. R.
Piletic
, and
M. D.
Fayer
,
J. Chem. Phys.
118
,
2270
(
2003
).
154.
A. A.
Bakulin
,
O.
Selig
,
H. J.
Bakker
,
Y. L. A.
Rezus
,
C.
Müller
,
T.
Glaser
,
R.
Lovrincic
,
Z.
Sun
,
Z.
Chen
,
A.
Walsh
,
J. M.
Frost
, and
T. L. C.
Jansen
,
J. Phys. Chem. Lett.
6
,
3663
(
2015
).
155.
M.
Ji
and
K. J.
Gaffney
,
J. Chem. Phys.
134
,
044516
(
2011
).
156.
Y. S.
Lin
,
P. A.
Pieniazek
,
M.
Yang
, and
J. L.
Skinner
,
J. Chem. Phys.
132
,
174505
(
2010
).
157.
T. l. C.
Jansen
and
J.
Knoester
,
Biophys. J.
94
,
1818
(
2008
).
158.
159.
A.
Tokmakoff
,
J. Chem. Phys.
105
,
1
(
1996
).
160.
O.
Golonzka
,
M.
Khalil
,
N.
Demirdöven
, and
A.
Tokmakoff
,
J. Chem. Phys.
115
,
10814
(
2001
).
161.
M. T.
Zanni
,
N.-H.
Ge
,
Y. S.
Kim
, and
R. M.
Hochstrasser
,
Proc. Natl. Acad. Sci. U. S. A.
98
,
11265
(
2001
).
162.
K.-K.
Lee
,
K.-H.
Park
,
S.
Park
,
S.-J.
Jeon
, and
M.
Cho
,
J. Phys. Chem. B
115
,
5456
(
2011
).
163.
J. R.
Schmidt
,
S. A.
Corcelli
, and
J. L.
Skinner
,
J. Chem. Phys.
123
,
044513
(
2005
).
165.
M.
Tsuboi
and
G. J.
Thomas
,
Appl. Spectrosc. Rev.
32
,
263
(
1997
).
166.
B. M.
Auer
and
J. L.
Skinner
,
J. Phys. Chem. B
113
,
4125
(
2009
).
167.
S. J.
Roeters
,
C. N.
van Dijk
,
A.
Torres-Knoop
,
E. H. G.
Backus
,
R. K.
Campen
,
M.
Bonn
, and
S.
Woutersen
,
J. Phys. Chem. A
117
,
6311
(
2013
).
168.
C.
Liang
and
T. l. C.
Jansen
,
J. Phys. Chem. B
117
,
6937
(
2013
).
169.
A.
Morita
,
J. Phys. Chem. B
110
,
3158
(
2006
).
170.
P. C.
Singh
,
S.
Nihonyanagi
,
S.
Yamaguchi
, and
T.
Tahara
,
J. Chem. Phys.
137
,
094706
(
2012
).
171.
J. E.
Laaser
and
M. T.
Zanni
,
J. Phys. Chem. A
117
,
5875
(
2013
).
172.
Y.
Ni
,
S. M.
Gruenbaum
, and
J. L.
Skinner
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
1992
(
2013
).
173.
S.
Roy
,
S. M.
Gruenbaum
, and
J. L.
Skinner
,
J. Chem. Phys.
141
,
22D505
(
2014
).
174.
J. E.
Laaser
,
D. R.
Skoff
,
J.-J.
Ho
,
Y.
Joo
,
A. L.
Serrano
,
J. D.
Steinkruger
,
P.
Gopalan
,
S. H.
Gellman
, and
M. T.
Zanni
,
J. Am. Chem. Soc.
136
,
956
(
2014
).
175.
E. H. G.
Backus
,
N.
Garcia-Araez
,
M.
Bonn
, and
H. J.
Bakker
,
J. Phys. Chem. C
116
,
23351
(
2012
).
177.
M.
Musso
,
H.
Torii
,
P.
Ottaviani
,
A.
Asenbaum
, and
M. G.
Giorgini
,
J. Phys. Chem. A
106
,
10152
(
2002
).
178.
S.
Roy
,
M. S.
Pshenichnikov
, and
T. L. C.
Jansen
,
J. Phys. Chem. B
115
,
5431
(
2011
).
179.
K.
Kwac
,
H.
Lee
, and
M.
Cho
,
J. Chem. Phys.
120
,
1477
(
2004
).
180.
M. H.
Farag
,
B. J.
Hoenders
,
J.
Knoester
, and
T. L. C.
Jansen
,
J. Chem. Phys.
146
,
234201
(
2017
).
181.
M.
Dinpajooh
and
D. V.
Matyushov
,
J. Phys. Chem. B
118
,
7925
(
2014
).
182.
A.
Anda
,
L.
De Vico
,
T.
Hansen
, and
D.
Abramavičius
,
J. Chem. Theory Comput.
12
,
5979
(
2016
).
183.
T. L. C.
Jansen
and
J.
Knoester
,
J. Chem. Phys.
127
,
234502
(
2007
).
184.
J. F.
Cahoon
,
K. R.
Sawyer
,
J. P.
Schlegel
, and
C. B.
Harris
,
Science
319
,
1820
(
2008
).
185.
F.
Šanda
and
S.
Mukamel
,
J. Chem. Phys.
125
,
014507
(
2006
).
186.
Y. S.
Kim
and
R. M.
Hochstrasser
,
J. Phys. Chem. B
110
,
8531
(
2006
).
187.
J.
Zheng
,
K.
Kwak
,
J. B.
Asbury
,
X.
Chen
,
I. R.
Piletic
, and
M. D.
Fayer
,
Science
309
,
1338
(
2005
).
188.
Y. S.
Kim
and
R. M.
Hochstrasser
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
11185
(
2005
).
189.
S.
Woutersen
,
Y.
Mu
,
G.
Stock
, and
P.
Hamm
,
Chem. Phys.
266
,
137
(
2001
).
190.
C.
Scheurer
and
T.
Steinel
,
ChemPhysChem
8
,
503
(
2007
).
191.
K.
Kwak
,
J.
Zheng
,
H.
Cang
, and
M. D.
Fayer
,
J. Phys. Chem. B
110
,
19998
(
2006
).
192.
O.
Kühn
and
Y.
Tanimura
,
J. Chem. Phys.
119
,
2155
(
2003
).
193.
R.
Kubo
,
J. Phys. Soc. Jpn.
17
,
1100
(
1962
).
194.
T.
Steffen
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
69
,
3115
(
2000
).
195.
Y.
Tanimura
and
T.
Steffen
,
J. Phys. Soc. Jpn.
69
,
4095
(
2000
).
196.
Y.
Suzuki
and
Y.
Tanimura
,
J. Chem. Phys.
119
,
1650
(
2003
).
197.
T.
Kato
and
Y.
Tanimura
,
J. Chem. Phys.
120
,
260
(
2003
).
198.
D.
Cringus
,
T. l. C.
Jansen
,
M. S.
Pshenichnikov
, and
D. A.
Wiersma
,
J. Chem. Phys.
127
,
084507
(
2007
).
199.
C. P.
van der Vegte
,
A. G.
Dijkstra
,
J.
Knoester
, and
T. L. C.
Jansen
,
J. Phys. Chem. A
117
,
5970
(
2013
).
200.
P. L.
McRobbie
,
G.
Hanna
,
Q.
Shi
, and
E.
Geva
,
Acc. Chem. Res.
42
,
1299
(
2009
).
201.
T. L. C.
Jansen
,
J. Phys. Chem. A
122
,
172
(
2018
).
202.
F.
Ding
and
M. T.
Zanni
,
Chem. Phys.
341
,
95
(
2007
).
203.
S.
Garrett-Roe
and
P.
Hamm
,
J. Chem. Phys.
130
,
164510
(
2009
).
204.
S.
Garrett-Roe
and
P.
Hamm
,
Acc. Chem. Res.
42
,
1412
(
2009
).
205.
L. J. G. W.
van Wilderen
,
A. T.
Messmer
, and
J.
Bredenbeck
,
Angew. Chem., Int. Ed.
53
,
2667
(
2014
).
206.
L. J. G. W.
van Wilderen
and
J.
Bredenbeck
,
Angew. Chem., Int. Ed.
54
,
11624
(
2015
).
207.
J. N.
Mastron
and
A.
Tokmakoff
,
J. Phys. Chem. A
120
,
9178
(
2016
).
208.
C. R.
Baiz
,
D.
Schach
, and
A.
Tokmakoff
,
Opt. Express
22
,
18724
(
2014
).
209.
T.
Hasegawa
and
Y.
Tanimura
,
J. Phys. Chem. B
115
,
5545
(
2011
).
210.
P. L.
Silvestrelli
,
M.
Bernasconi
, and
M.
Parrinello
,
Chem. Phys. Lett.
277
,
478
(
1997
).
211.
M. K.
Ghosh
,
J.
Lee
,
C. H.
Choi
, and
M.
Cho
,
J. Phys. Chem. A
116
,
8965
(
2012
).
212.
M.
Sulpizi
,
M.
Salanne
,
M.
Sprik
, and
M.-P.
Gaigeot
,
J. Phys. Chem. Lett.
4
,
83
(
2013
).
213.
M.
Thomas
,
M.
Brehm
,
R.
Fligg
,
P.
Vöhringer
, and
B.
Kirchner
,
Phys. Chem. Chem. Phys.
15
,
6608
(
2013
).
214.
A.
Ishizaki
and
Y.
Tanimura
,
J. Chem. Phys.
125
,
084501
(
2006
).
215.
A.
Ishizaki
and
Y.
Tanimura
,
J. Phys. Chem. A
111
,
9269
(
2007
).
216.
A. O.
Caldeira
and
A. J.
Leggett
,
Physica A
121
,
587
(
1983
).
217.
J.-J.
Ding
,
J.
Xu
,
J.
Hu
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
135
,
164107
(
2011
).