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.

## I. INTRODUCTION

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} + *iγ*_{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 spectroscopy^{1–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.

**m**generally depend on the vibrational coordinate

*q*, they can be Taylor-expanded in the form

**m**(

*t*) =

**m**

_{0}+ (

*∂*

**m**/

*∂q*)

_{0}

*q*+ $\cdots \u2009$ 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

*∂*

**μ**/

*∂q*)

_{0}and (

*∂*

**α**/

*∂q*)

_{0}are the transition dipole and polarizability, respectively, and the bar over, for example, $\u2202\mu /\u2202q02\xaf$ 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 Raman^{19} 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.

## II. FUNDAMENTALS

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.

### A. Third-order optical response function

**E**(

**r**,

*t*) is the superposition of the electric fields of the three incident (IR, THz, visible, or X-ray) pulses denoted as

**E**

_{1},

**E**

_{2}, and

**E**

_{3}. The total Hamiltonian of the system is the sum of the system Hamiltonian

*H*

_{0}in the absence of radiation and

*H*

_{int}(

*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., $\u2212\alpha ^:E(r,t)E(r,t)$ for Raman spectroscopy.

*ρ*(

*t*) of the system as follows:

*A*(

*t*) through the expectation value $Tr\xc2\rho (t)$, where Tr denotes the trace of a matrix and $\xc2$ 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}≡ (

*E*

_{a}−

*E*

_{b})/

*ℏ*determined by the energy difference of the two states.

*H*

_{int}(

*t*) as the perturbation to the reference Hamiltonian

*H*

_{0}, 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

*H*

_{int}(

*t*) and is given by

*G*

_{0}(

*t*) = exp(−

*iL*

_{0}

*t*/

*ℏ*) is the time evolution operator in the absence of radiation and the Liouville operators are defined as

*L*

_{a}

*A*= [

*H*

_{a},

*A*] for

*a*= 0 or int. According to Eq. (4), the system initially defined by

*ρ*(

*t*

_{0}) evolves freely without perturbation for

*τ*

_{1}−

*t*

_{0}, as given by

*G*

_{0}(

*τ*

_{1}−

*t*

_{0}), and then interacts with the radiation at time

*τ*

_{1}, as given by

*L*

_{int}(

*τ*

_{1}). This sequence is repeated

*n*times until the final field-matter interaction at

*τ*

_{n}, as given by

*L*

_{int}(

*τ*

_{n}). Finally, the system evolves freely until the observation time

*t*according to

*G*

_{0}(

*t*−

*τ*

_{n}). The multiple integrals over

*τ*

_{1}, …,

*τ*

_{n}account for all possible interaction times under the time ordering condition

*t*

_{0}≤

*τ*

_{1}≤ … ≤

*τ*

_{n}≤

*t*.

*ρ*(

*t*) in Eq. (4) gives rise to the corresponding polarization $P(n)(r,t)=Tr\mu ^\rho (n)(t)$ in the system as follows:

*n*th order optical response function given by

*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)}(

*t*

_{3}, …,

*t*

_{1}), the latter being the fourth rank tensor. Note that the time variables

*t*

_{1}, …,

*t*

_{n−1}in Eqs. (5) and (6) are time intervals between consecutive field-matter interactions related to

*τ*

_{1}, …,

*τ*

_{n}in Eq. (4) as

*t*

_{m}=

*τ*

_{m+1}−

*τ*

_{m}(1 ≤

*m*≤

*n*− 1), while

*t*

_{n}=

*t*−

*τ*

_{n}is the time elapsed since the last field-matter interaction. Therefore,

*t*

_{1}, …,

*t*

_{n}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.

**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 as

^{30,31}

*n*(

*ω*) is the refractive index of the medium and $Ps(n)(t)$ is the polarization component propagating with wave vector

**k**

_{s}and frequency

*ω*

_{s}that are one of the combinations ±

**k**

_{1}±

**k**

_{2}± ⋯ ±

**k**

_{n}and ±

*ω*

_{1}±

*ω*

_{2}± ⋯ ±

*ω*

_{n}, respectively. Note that Eq. (7) gives the approximate signal field arising from a single Fourier component of the

*n*th order polarization expanded as

^{8,31}

**k**

_{s}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.

### B. Response function of a three-level system

2D vibrational spectroscopy usually induces transitions up to the second vibrational excited state. Therefore, a three-level system with eigenstates $g$, $e$, and $f\u2009\u2009$ 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}.

*ℏω*

_{m}is the energy of state

*m*in the absence of bath,

*V*

_{m}(

**q**) is the chromophore-bath interaction energy of the state

*m*that depends on the bath degrees of freedom

**q**,

*H*

_{B}(

**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 terms

^{31,34}

**R**

_{i}(

*t*

_{3},

*t*

_{2},

*t*

_{1}) are given by

^{8,35}

*g*. Here,

**μ**

_{ab}is the transition dipole between states

*a*and

*b*obtained by the repeated insertion of the closure relation $\u2211mmm=1$ in Eq. (6), which is assumed to be independent of the bath coordinate (Condon approximation), $\u210f\omega \xafab=\u210f\omega a\u2212\omega b+Va(q)\u2212Vb(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)\u2212Vm(q)B$. Therefore, the response function is composed of multiple quantum transition pathways represented by individual

**R**

_{i}, 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 (

*F*

_{1−4}). Throughout this paper, third order double-quantum or zero-quantum spectroscopies are not taken into consideration.

*U*

_{ab}(

**q**) =

*U*

_{a}(

**q**) −

*U*

_{b}(

**q**). Alternatively, we can invoke the second-order cumulant expansion approximation, which becomes exact when the fluctuation of the energy gap $\u210f\omega \xafab$ obeys the Gaussian statistics, to obtain

^{8}

**R**

_{3}and

**R**

_{4}in Eq. (11), that represent coherence evolution during

*t*

_{2}, 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.,

*U*

_{fg}(

*t*) ≅ 2

*U*

_{eg}(

*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

*t*

_{1}and

*t*

_{3}have the opposite or the same signs, respectively.

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.

## III. CLASSICAL MECHANICAL APPROACHES TO NONLINEAR RESPONSE FUNCTIONS

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.

### A. Equilibrium molecular dynamics approach

^{38}

^{,}

*X*and

*Y*are physical variables and the Poisson bracket is defined as follows:

*Y*is the equilibrium distribution function of the system

*ρ*

_{eq}=

*e*

^{−βH}/

*Z*, where Z is the canonical partition function, Eq. (15) becomes ${X,\rho eq}PB=\u2212\beta X\u0307\rho eq$. Here,

*β*is the reciprocal temperature of the system and $X\u0307$ 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

*B*

^{7}by using Eq. (14) in the quantum linear response function $R(1)(t)=i\u210f\theta (t)TrA(t)[B(0),\rho eq]$, which corresponds to

*n*= 1 in Eq. (6),

^{39}

*α*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.

^{19,38}which is the second-order response function of the polarizability, is given as

^{39–45}

^{39–41,46}

^{47–49}

*μ*denotes the dipole moment of the system.

^{43,50}(of which the Newtonian time evolution is an example) and the chain rule of the form $\u2202\mu (t\u2032)/\u2202p\alpha (t)=\u2211\beta \u2202\mu (t\u2032)/\u2202q\beta (t\u2032)\u2202q\beta (t\u2032)/\u2202p\alpha (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 as

^{38,39}

**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}

### B. Nonequilibrium molecular dynamics approach

#### 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 H_{2}O. 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.

*B*

_{−}and a dimensionless parameter

*ε*, we have the following identity:

^{56}

*R*

^{(1)}(

*t*) (for

*t*> 0) of a physical variable

*A*to a perturbation

*B*, we have

^{56–58}

*A*(

*t*) on the perturbed trajectory given by

*H*

_{0}−

*εBδ*(

*t*), whereas the second term is the expectation value of

*A*on the trajectory expressed by the unperturbed Hamiltonian

*H*

_{0}.

^{46}

*et al.*have applied the nonequilibrium finite-field method to the 1D and 2D Raman spectra of liquids.

^{57,58}

^{48}

*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:

#### 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.

^{60}In the hybrid method, the second-order response function for the 2D Raman spectra given by Eqs. (18) and (24) can be written as

^{61}

*t*

_{1}+

*t*

_{2}) on the nonequilibrium trajectory generated by the external fields applied at

*t*

_{1}.

^{61}

^{47}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

*μ*(

*t*

_{2}+

*t*

_{3})

_{E(t2)}and $\mu \u0307(\u2212t1)E(0)$ denote the dipole moment at

*t*

_{2}+

*t*

_{3}on the perturbed trajectory with the Hamiltonian

*H*

_{0}−

*εμEδ*(

*t*−

*t*

_{2}) and the time derivative of dipole moment at −

*t*

_{1}on the backward perturbed trajectory

*H*

_{0}−

*εμEδ*(

*t*), respectively, while $\mu \u0307(\u2212t1)\mu \u0307(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

*t*

_{2}values because the third-order response functions at different

*t*

_{2}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 **k**_{I} = −**k**_{1} + **k**_{2} + **k**_{3} and **k**_{II} = +**k**_{1} − **k**_{2} + **k**_{3}. However, the calculated third-order response function consists of all the wave vector components generated by three incident electric fields, i.e., **k**_{III} = +**k**_{1} + **k**_{2} − **k**_{3} and **k**_{IV} = +**k**_{1} + **k**_{2} + **k**_{3} as well as **k**_{I} and **k**_{II}. Therefore, the components with **k**_{III} and **k**_{IV} 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, **k**_{III} and **k**_{IV} components have been eliminated by exploiting the fact that the three-dimensional Fourier transformed spectra of **k**_{III} and **k**_{IV} components show higher frequency oscillation in *ω*_{2} than those of **k**_{I} and **k**_{II}.^{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.

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}

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 Saito^{65,66} and Jeon and Cho.^{67,68}

Recently, the nonequilibrium molecular dynamics method has also been applied to the intramolecular vibrational modes of H_{2}O 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 H_{2}O in water.

### C. Limitations

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 *t*_{3}.

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}

## IV. NUMERICAL INTEGRATION OF THE SCHRÖDINGER EQUATION

^{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

*B*

_{m}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, $\mu \u2192n(t)$, transition polarizabilities, $\alpha \xaf\xafn(t)$, and anharmonicities, Δ

_{n}(

*t*). The different local vibrations are mixed by their mutual coupling,

*J*

_{nm}(

*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\u2192(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 ices^{94,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}

### A. Hamiltonian parameterization, frequency maps

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 potential^{104–108} or electric field^{102,109–125} and sometimes the electric field gradient^{126–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.

### B. Mappings for interactions and the related techniques

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 reorientation^{153–156} and vibrational energy transfer^{72,153,157} can be monitored by measuring 2D IR spectra using different polarizations^{158,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 geometries^{160–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 generation^{17,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 effect^{147,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}

### C. Limitations—Bandwidth, high-frequency modes

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 dynamics^{152,178–182} and chemical exchange^{101,155,179,183–192} arising when the bath coordinates are not harmonic. Commonly used methods invoking the second-order cumulant approximation^{34,193} or methods assuming a harmonic bath spectral density^{194–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.

## V. PERSPECTIVES AND A FEW CONCLUDING REMARKS

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 experimentally^{202} and simulated^{203,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 spectroscopy^{205,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 models^{21,209} or quantum mechanical potentials. In this regard, recently proposed *ab initio* theories of vibrational solvatochromism^{135–137} and direct QM/MM or full *ab initio* MD simulations of vibrational spectra^{210–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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

*Two-Dimensional Optical Spectroscopy*

*Advances in Multi-Photon Processes and Spectroscopy*

*Principles of Nonlinear Optical Spectroscopy*

*Regular and Chaotic Dynamics*