Simulating two-dimensional correlation spectroscopies with third-order infrared and fifth-order infrared–Raman processes of liquid water

The delineation of spectroscopic line shapes of liquid water into vibrational modes and their interactions has been the subject of numerous experimental and theoretical studies.^1–4 The uniqueness of water arises from its high-frequency intramolecular modes that promote bond forming and bond breaking through complex hydrogen bonding (HB)^5–8 and its low-frequency intermolecular modes that realize irreversible nuclear motion.^9–12 The interplay between these modes plays key roles in chemical reaction processes.^13,14

Experimentally, infrared (IR)^15,16 and third-order off-resonant Raman spectroscopies^17–19 have been utilized for the analysis of vibrational motions of liquid water. The observables of these methods are defined in terms of the dipole $μ̂$ and polarizability $Π̂$ as $RIR(1)(t1)=i⟨[μ̂(t1),μ̂]⟩/ℏ$ and $RRaman(3)(t1)=i⟨[Π̂(t1),Π̂]⟩/ℏ$⁠, respectively, as functions of a single time duration t₁, which is referred to as one-dimensional (1D) spectroscopy.²⁰ The IR process has been utilized not only for investigations of the intramolecular OH stretching motion (∼3700 cm⁻¹) and the intramolecular HOH bending motion (∼1600 cm⁻¹) but also for those of the HB intermolecular translation (vibrational) motion (33–100 cm⁻¹) and inter-HB molecular vibrational motion (33–400 cm⁻¹), which are also regarded as THz processes. The Raman process, which exploits the large polarizability of water molecules, has been used to probe both intramolecular and intermolecular modes. This approach is particularly useful for investigations of intermolecular modes because of its superior time resolution and intensity in such a low-frequency region. Although 1D spectroscopic approaches are useful for the classification of vibrational line shapes of liquids into intermolecular and intramolecular modes, it is unclear whether the peaks that we investigate are mutually coupled and whether the widths of the peaks that we measure are of an inhomogeneous origin or a homogeneous origin as these peaks are usually broad and overlap.

Two-dimensional (2D) experiments, in which an additional time correlation is imprinted on the system response, are capable of determining anharmonic coupling among vibrational modes and homogeneous linewidths that are needed to understand the degree of coupling to the bath (i.e., the surrounding molecules with respect to the excited transition).^20–24 Recent developments of experimental techniques have significantly augmented the study of both intermolecular and intramolecular modes in the 0–4000 cm⁻¹ frequency range through 2D IR,^25–29 2D THz–Raman spectroscopies,^30–33 and 2D THz–IR–visible (TIV) spectroscopy,^34,35 which is equivalent to 2D IR–IR–Raman spectroscopy.^36,37 Although 2D IR and 2D THz–Raman experiments provide information about intramolecular and intermolecular degrees of freedom, respectively, neither approach can directly measure the anharmonic coupling between intermolecular and intramolecular modes. In contrast, 2D TIV, which utilizes light sources to excite both intermolecular and intramolecular modes, allows us to investigate their mode–mode coupling mechanism. Note that the analysis of both 2D IR–Raman and 2D THz–Raman spectroscopies defined in terms of three-body response functions, such as $RIIR(3)(t2,t1)=−⟨[[Π̂(t2+t1),μ̂(t1)],μ̂(0)]⟩/ℏ2$⁠, is far more complex than that of third-order 2D IR spectroscopy defined in terms of four-body correlation functions, such as $R2DIR(3)(t3,t2,t1)=−iμ̂(t3+t2+t1),μ̂(t2+t1),μ̂(t1),μ̂(0)/ℏ3$⁠, because contributions from the nonlinear polarizability and dipole moments involved in laser interactions give rise to additional peaks, and nondiagonal spectral peaks are not necessary to represent mode–mode coupling peaks, as in the case of 2D IR spectra.^38,39

In this study, we extend our analysis of 2D THz–Raman^40–43 and 2D IR–Raman spectroscopies^36,37 by including an additional time duration for laser excitation and calculate the fifth-order 2D IR–IR–Raman–Raman spectrum defined as $RIIRR(5)(t3,t2,t1)=−iΠ̂(t3+t2+t1),Π̂(t2+t1),μ̂(t1),μ̂(0)/ℏ3$⁠, in addition to the third-order 2D IR spectrum defined as $R2DIR(3)(t3,t2,t1)$⁠. Through simulation, we attempt to clarify the roles of anharmonic mode–mode coupling and vibrational dephasing, as well as thermal relaxation, among the intermolecular and intramolecular modes of balk water.

It should be noted that simulating the 2D spectrum of liquid water for a mixture of intramolecular and intermolecular modes is challenging, especially on account of the peak-splitting in third-order 2D IR spectra caused by the transitions between the |0⟩ → |1⟩ → |2⟩ and |0⟩ → |1⟩ → |0⟩ states, where |n⟩ is the nth vibrational energy state.^44–46 High-frequency intramolecular motions must be addressed within the framework of quantum dynamics theory, and their dynamics involve complex hydrogen bond networks for 2D spectroscopic measurements. These motions have only been studied from classical MD simulations^46,47 or classical MD-based quantum simulations, such as using stochastic theory⁴⁸ and wave function-based theory.^49,50 Therefore, peak profiles of experimentally obtained 2D IR spectra—particularly for mode–mode coupling peaks^25–28—have not been accurately predicted from a theoretical approach.

As a practical approach, the Brownian oscillator (BO) model has been employed for the analysis of both linear^51–53 and nonlinear spectra.^20–22 In this approach, vibrational modes representing spectroscopic properties of interest are described as functions of molecular coordinates, whereas environmental molecular motions are described using heat baths that exert thermal fluctuations and dissipation on the vibrational modes. Although analyses based on MD simulations have shown that vibrational relaxation and dephasing are important mechanisms for characterizing molecular motions,^54,55 the inclusion of nonlinear and non-Markovian interactions with the heat baths,^56–58 as well as the anharmonic mode–mode interactions,^22,59 is also significant. Therefore, for vibrational motions of liquid water, we introduced the multimode nonlinear BO model, which consists of interacting anharmonic oscillators coupled to their respective heat baths with linear–linear (LL) and square–linear (SL) system–bath (SB) interactions to describe vibrational relaxation and dephasing.^36,37,43 To guarantee the equilibrium state of the system at finite temperature, the model must be solved accurately such that the dephasing (fluctuations) and relaxation (dissipation) satisfy the quantum version of the fluctuation–dissipation theorem.^60–62 

The model parameters were then chosen to fit classical MD simulation for 2D IR–Raman spectra conducted using the polarizable water model for inter- and intramolecular vibrational spectroscopy (POLI2VS);⁶³ we then modified the anharmonicity of intramolecular modes to fit experimentally obtained 1D Raman and IR spectra,³⁷ which agree with quantum MD results obtained using the POLI2VS model⁶⁴ and with the ab initio MD results that have recently become available.^65–67 To simulate 2D correlation spectra that detect the relaxation time of vibrational dephasing, it is necessary to deal with non-perturbative and non-Markovian SB interactions, especially for intramolecular modes, in which vibrational motion and the environment are quantum-mechanically entangled.^44,45 Although various hierarchical equation of motion (HEOM) approaches have been developed for investigating nonlinear spectroscopies,^60–62 here, we employ discretized hierarchical equations of motion in mixed Liouville–Wigner space (DHEOM–MLWS) specifically developed for the water BO model.³⁷ 

This paper is organized as follows. In Sec. II, we explain our multimode LL+SL BO model and the parameter set used for simulations. In Sec. III, we give definitions of nonlinear optical observables to be computed. In Sec. IV, we present the results for third-order 2D correlation IR and fifth-order 2D correlation IR–IR–Raman–Raman spectra based on two-mode calculations. Section V is devoted to our concluding remarks.

The details of the model and the DHEOM–MLWS for the quantum case have been explained in Refs. 36 and 37. Here, we present an overview of the model because it is helpful for the analysis of calculated 2D spectra.

The multimode nonlinear BO model for liquid water consists of four primary oscillator modes representing (1) intramolecular OH stretching (“stretching”), (2) intramolecular HOH bending (“bending”), (3) HB–intermolecular librational (“librational”), and (4) HB–intermolecular translational (“translational”) motions.^36,37 Note that OH stretching modes consist of symmetry and anti-symmetry modes, but here, we treat the stretching modes as a single mode because we construct the model based on the broadened linear and nonlinear spectra as benchmarks, as in previous studies. Although it is possible to model the two modes separately by MD simulation with the machine learning approach,⁶⁸ solving such a model is computationally demanding, so here, we treat them as one mode; thus, they are described by dimensionless vibrational coordinates $q⃗=(…,qs,…)$⁠, with s = 1, …, 4 indexing the four vibrational modes. The total Hamiltonian is then expressed as

where

is the Hamiltonian for the sth mode with mass m_s, coordinate q_s, momentum p_s, fundamental frequency ω_s, and anharmonicity of potential $gs3$⁠; $gs2s′$ and $gss′2$ are the anharmonic coupling between modes s and s′;

is the bath Hamiltonian for the sth mode with the momentum, coordinate, mass, and frequency of the j_sth bath oscillator given by $pjs$⁠, $xjs$⁠, $mjs$ and $ωjs$⁠, respectively;

is the SB interaction, which consists of linear–linear (LL) and square–linear (SL) SB interactions; and $Vs(qs)≡VLL(s)qs+VSL(s)qs2/2$⁠, with coupling strengths of $VLL(s)$⁠, $VSL(s)$⁠, and $αjs$⁠.^{43–45,56–58} The last term of the bath Hamiltonian is a counter-term that maintains the translational symmetry of the system.⁶¹ While the LL interaction contributes mainly to energy relaxation, the SL system–bath interaction leads to vibrational dephasing in the slow modulation case due to the frequency fluctuation of system vibrations.^43–45 The sum of the bath coordinates, $Xs≡∑jsαjsxjs$⁠, acts as a collective coordinate that modulates mode s.^60–62 We introduce the spectral distribution function (SDF) $Js(ω)≡∑jsαjs2ℏδ(ω−ωjs)/2mjsωjs$ to characterize the environmental fluctuation. We assume that J_s(ω) has an Ohmic form with a Lorentzian cutoff,

Here, ζ_s is the system–bath coupling strength, and γ_s represents the width of the spectral distribution for mode s that relates to the vibrational dephasing time, defined as τ_s = 1/γ_s.

The dipole operator and polarizability for the sth mode are defined as

and

where μ_s and μ_ss′ are the linear and nonlinear elements of the dipole moment and Π_s and Π_ss′ are those of the polarizability, respectively. The vibrational modes interact through the mechanical anharmonic coupling (MAHC) described by $gs2s′$ and $gss′2$ and the electric anharmonic coupling (EAHC) described by μ_ss′ and Π_ss′. Although the effects of MAHC and EAHC in 2D spectra have been discussed analytically,^{38,39,69–72} no quantum theory exists for complex systems involving complex system–bath interactions, such as those in this study.

In the present case, the difference between fundamental frequencies of intermolecular and intramolecular modes is very large, and we can excite each mode separately by choosing the light sources: a high-frequency intramolecular mode is solely excited by an IR pulse, whereas a low-frequency intermolecular mode is solely excited by a Raman or THz pulse. Note that although the contribution from EAHC arises from nonlinear dipole elements or nonlinear polarizability elements are not important for a spectroscopic approach based on the third-order response functions for s = 1 and 2 and they can be ignored,^38,39 we used the same parameter values presented in Ref. 37 as given in Tables I and II because the computational cost in our approach is the same. Besides these, we used the same parameter values presented in Ref. 37. The element μ₂₃ is important even in 1D IR because this is the cause of the combination band at the frequency ω = 2100 cm⁻¹.

We chose a set of parameters for the model based on the MD results of 2D spectra and experimentally obtained 1D spectra, taking advantage of the fact that the 2D profiles differ significantly for different dynamic processes.^36,37

We consider 2D IR and IR–Raman spectroscopies whose pulse sequences are depicted in Fig. 1. These spectra can be calculated from the three-body nonlinear response function defined as^61,62

where the operators $Â$⁠, $B̂$⁠, $Ĉ$⁠, and $D̂$ are either the dipole moment $μ̂$ or the polarizability $Π̂$⁠. Here, we have employed the hyperoperator ^× defined as $Â×ρ̂̂≡[Â,ρ̂]$⁠, $G(t)$ is the Green’s function of the total Hamiltonian without a laser interaction, and $ρ̂eq$ is the equilibrium state. Using this expression, the third-order IR signal is calculated from $Â∼D̂=μ̂$ and the fifth-order IR–Raman signal is calculated from $Â=B̂=Π̂$ and $Ĉ=D̂=μ̂$⁠, respectively.^57–59 Note that, in our model calculation, the role of IR and Raman excitation for the intramolecular modes and that of THz and Raman excitation for the intermolecular modes are similar because nonlinear elements of each mode, such as μ_ss and Π_ss, do not play a major role in determining the profiles of the 2D spectrum, which is defined as Eq. (8). This is because nonlinear spectroscopy based on a three-body response function, such as $RIIR(3)(t2,t1)=−⟨[[Π̂(t2+t1),μ̂(t1)],μ̂(0)]⟩/ℏ2$⁠, is described as an odd number of optical operators, whereas that based on a four-body response function, such as $RIIRR(5)(t3,t2,t1)=−iΠ̂(t3+t2+t1),Π̂(t2+t1),μ̂(t1),μ̂(0)/ℏ3$⁠, is described as an even number of optical operators. When the anharmonicity of the mode described by $gs3qs3$ is weak, the thermal statistical average included in the calculation of the response function is performed as a Gaussian integral of q_s. Therefore, $RIIR(3)(t2,t1)$ must be of even order function of q_s so that it does not disappear in the Gaussian integral, and a contribution of Π_ss or μ_ss is required. In contrast, $RIIRR(5)(t3,t2,t1)$ is fourth order in q_s even without Π_ss or Π_ss, and if they contribute, it should be of the second order in Π_ss and/or μ_ss, such as Π_ssμ_ss.²¹ This indicates that the MAHC contribution, which is often referred to as the non-Condon effect,⁷³ in $RIIR(3)(t2,t1)$ and $RIIRR(5)(t3,t2,t1)$ measurements is minor. Thus, although it would be difficult to conduct an actual measurement, the profile of a seventh-order 2D Raman spectrum (i.e., $Â∼D̂=Π̂$⁠)^57,58,74 should be similar to that of third-order 2D IR.

The calculated signals are plotted as the 2D correlation spectrum that is obtained by adding the two terms corresponding to the rephasing and nonrephasing parts of R(t₃, t₂, t₁) in equal weights. A commonly used definition of a 2D correlation spectrum is expressed as^75–77 

where R_R and R_NR are the rephasing and nonrephasing responses that are experimentally detected separately using phase-matching conditions. Theoretically, although we can separate R_R and R_NR by choosing the specific Liouville paths in the energy state model, the process is not easy when we calculate response functions in the coordinate space. Therefore, here, we eliminate the undesired rephasing contribution by utilizing the Fourier transformation of t₂ for I_C(ω₃, t₂, ω₁) with R_R(t₃, t₂, t₁) = R_NR(t₃, t₂, t₁) = R(t₃, t₂, t₁) to remove the oscillating contribution in the t₂ period with frequency 2ω_s, where ω_s is the frequency of the targeting mode s.^78,79 Note that we can also separate the plain response function in terms of different phase-matching conditions using a phase-matching transformation,⁸⁰ but the Fourier-based method is simpler when only a 2D correlation spectrum is calculated.

After removing the dispersive contribution of the response function, we obtain the pure absorptive line shape. We observe a distorted line shape, called the phase-twisted line, which exhibits the non-Markovian and nonperturbative dynamics of the environment.^{22–24,75–77} To investigate not only anharmonic mode–mode coupling but also vibrational dephasing and thermal relaxation effects at finite temperature, we must employ nonperturbative and non-Markovian open quantum dynamics theory, which enables us to treat the quantum entanglement between the system and bath.⁶² Therefore, we employ discretized hierarchical equations of motion in mixed Liouville–Wigner space (DHEOM–MLWS): Lagrange–Hermite mesh discretization is employed in the Liouville space of intramolecular modes, and Lagrange–Hermite mesh discretization and Hermite discretization are employed in the Wigner space of intermolecular modes.³⁷ 

Although the DHEOM–MLWS for the water BO model was created to treat the above four modes simultaneously, the computational time required for 2D vibrational spectra is about 1000 times longer than that for linear spectra, even with a fixed t₂. Therefore, here, we calculate 2D spectra for mode–mode coupling, which can be investigated for two modes, as in combinations of stretching–bending (1–2) and bending–librational (2–3) for 2D correlation IR spectra and stretching–librational (1–3) and stretching–translational (1–4) for 2D correlation IR–Raman spectra. To calculate the spectrum, we integrate the DHEOM–MWLS using the same numerical setup as in Ref. 37 with the parameter set presented in Tables I and II. Thus, the calculated 1D IR and 1D Raman spectra (not shown) are identical to those presented in Ref. 37.

Note that the dipole moment and polarizability have similar functional forms in this model calculation. Therefore, if the same modes are excited, 2D Raman–Raman–IR–IR or 2D IR–Raman–IR–Raman results exhibit a profile that is similar to the present 2D IR–IR–Raman–Raman results.

Figure 2 illustrates 2D correlation spectra calculated for stretching–bending (1–2) modes for which experimental results are available.^26–28 Note that previous theoretical investigations for water 2D spectroscopy have been limited to stretching modes only for quantum calculations based on classical MD^48–50 and individual stretching and bending modes for classical MD.^46,47 The present investigations provide the first quantum simulation of 2D spectra for intramolecular–intramolecular and intramolecular–intermolecular mode–mode couplings. To conduct simulations, we set $μ̂=μ̂1+μ̂2$⁠. Each peak near (ω₁, ω₃) = (3400 cm⁻¹, 3400 cm⁻¹) in the upper panel represents the stretching mode; the red positive peaks and the blue negative peaks correspond to the transitions between the |0⟩_s → _s|1⟩_s → |0⟩_s and |0⟩_s → |1⟩_s → |2⟩_s states for s = 1, respectively, where |n⟩_s is the nth vibrational energy state of the sth mode.

In the lower panel of each panel of Fig. 2, we can see cross peaks representing the stretching → bending transition (e.g., |0⟩₁|0⟩₂ → |1⟩₁|0⟩₂ ∼ |0⟩₁|1⟩₂ → |0⟩₁|0⟩₂).²² The coupling peaks that appear at t₂ = 0 are due to coherent transfer between the two modes and quickly disappear at t₂ = 50 fs when coherence is lost. The coupling peaks that appear after t₂ ≥ 100 fs are due to population transfer, and their intensity increases with t₂.

The direction of the orange nodal lines represents the extent of noise correlation (non-Markovian effects) between the vibrational coherences of periods t₁ and t₃. The direction parallel to the ω₁ axis is the uncorrelated case, whereas the direction parallel to the ω₁ = ω₃ line is the fully correlated case.^75–77 The width of the peak parallel to the ω₁ = ω₃ direction corresponds to inhomogeneous broadening, whereas that perpendicular to ω₁ = ω₃ is homogeneous broadening.^22,81

The spectral width of the stretching peaks in the upper panel of Fig. 2 is broadened due to the dephasing (spectral diffusion regime) that arises from the non-perturbative and non-Markovian effects of the SL interaction.⁴⁴ When t₂ becomes larger than the correlation time of the noise τ₁ = 1/γ₁ ≈ 160 fs, the correlation begins to weaken, and at t₂ = 500 fs, the nodal line becomes horizontal. The SL interaction used in this model causes the same effect as in stochastic theory, where τ₁ determines the dephasing time.⁶¹ We found that our result is consistent with those of previous experiments^25–28,82 and theoretical calculations.^46–50 

Upon increasing t₂, it was observed experimentally that the blue 0–1–2 peak extended in the small ω₃ direction for t₂ ≤ 250 fs, and a peak, referred to as the hot ground state, appeared above the red 0–1–0 peak for t₂ ≥ 250 fs.^25–28 In contrast, our results show that the blue peak extends only slightly in the negative ω₃ for the t₂ ≥ 500 fs direction due to the population relaxation arising from the SL interaction,⁴⁴ while no hot ground state peak is observed. Note that the inclusion of LL interactions in the stretching^44,45 and/or the presence of non-radiative decay states⁶⁸ would enhance the elongation. To predict such phenomena with a current handy model calculation, one must include various hypothetical factors that are considered important in determining the vibrational process of water, with reference to the experimentally obtained 2D spectral profile. We have computed 2D spectra for the case of including quartic Darling–Dennison (DD) coupling described as $g̃1222q12g22$^22,72 and for the case of stronger cubic (Fermi) coupling $(g̃ss′2=0.5)$⁠. While inclusion of the DD coupling does not change the 2D profile much (not shown), we found that, when the Fermi coupling is strong, the 2D profile is relatively close to the experimental results (see  Appendix A). However, some other information, such as quantum MD results, is needed to explore such a possibility in greater detail.

The cross peaks representing the stretching → bending transition in the lower panel of Fig. 2 are uncorrelated due to the process involved in the bending mode having a very short noise correlation time, 1/γ₂ ≈ 40 fs, and the positive and negative peaks are elongated in the ω₁ direction.^26–28 In  Appendix B, we depict the 2D spectrum of the bending mode. It can be seen that the time scale of the spectral diffusion of the cross-peak is determined by the mode in which it is faster.

Next, we present the calculated result for the bending → librational (2 → 3) transition peak in Fig. 3. We chose the first two IR pulses for the excitation of the bending mode and the last two IR (or THz) pulses for the detection of the librational modes. In the 2D correlation spectrum based on the third-order (four body) response function, we can interpret the 2D spectral profile as investigating the excitation of each vibrational mode in the ω₁ direction in the IR region and then detecting the dynamics of that excited motion after t₂ as the absorption spectrum in the ω₃ direction in the THz region. In  Appendix B, we depict the 2D spectrum of the bending mode for the same conditions as in Fig. 3.

In the present case, the first two IR pulses excite the bending mode (ω₁ = 1600 cm⁻¹) and the combination band (around ω₁ = 2100 cm⁻¹) through μ₂ and μ₂₃, respectively. Because no population transfer for 2 → 3 occurs at time t₂ = 0, we cannot observe a prominent peak in the ω₃ direction along the line ω₁ = 1600 cm⁻¹. The positive and negative peaks of the combination band peak around ω₁ = 2100 cm⁻¹ appear because the IR excitations described by μ₂₃ cause transitions of the librational vibrational states, |0⟩₂|0⟩₃ → |1⟩₂|1⟩₃ → |1⟩₂|0⟩₃ and |0⟩₂|0⟩₃ → |1⟩₂|1⟩₃ → |1⟩₂|2⟩₃, which appear as the positive upper and negative lower peaks, respectively.

At t₂ = 20 fs, the elongated peak along line ω₁ = 1600 cm⁻¹ appears as the result of the population transfer from 2 → 3. The elongation of the peak along this line reflects the fact that the librational mode is spread out. Because of this, the positive and negative combination band peaks are also elongated in the ω₃ direction. The nodal line of the combination band peaks in the high-frequency region is twisted, indicating that the contribution from the librational mode, which seems to dominate the high-frequency part of the combination band, can be correlated with the IR excitation of the librational mode. Because these combination band peaks appear due to the coherent process |0⟩₂ → |1⟩₂, they decay quickly and almost disappear at t₂ = 80 fs, whereas the 2 → 3 transition peak along line ω₁ = 1600 cm⁻¹ grows because this process arises from the relaxation. The fast relaxation found between these two modes is also consistent with other classical MD analyses.⁵⁵ 

Next, we display our calculated results for the mode–mode coupling of the high-frequency stretching mode to low-frequency intramolecular modes. Here, we chose IR pulses for the excitation of the stretching mode and Raman pulses for the detection of intermolecular modes because we believe that such an experimental setup is feasible with the current technology. However, because the difference in the mathematical treatment of Raman and IR (or THz) excitation is minor in our model calculation, the spectral profile would be similar regardless of the laser configuration, for example, for a configuration where Raman pulses excite the stretching mode and THz pulses detect intermolecular modes.

We now present the simulation result of fifth-order 2D IR–IR–Raman–Raman spectra for the stretching–translational coupling case. Note that this process was investigated theoretically and experimentally through the use of 2D THz–IR–visible (or 2D IR–IR–Raman) measurement.^34–37 The cross peaks predicted from the 2D THz–IR–visible spectra arise from both MAHC and EAHC contributions of the stretching–translational modes, but in the 2D IR–IR–Raman–Raman spectra, the cross peaks always arise from MAHC, whereas those from EHAC appear at different positions or do not appear, as can be seen from analytical expressions of response functions.^{38,39,69–71}

Moreover, although 2D THz–IR–visible measurement is designed to detect the coherence between the two different modes through mode–mode interactions, the present 2D IR–IR–Raman–Raman spectroscopy is designed to detect the population transfer among modes under dephasing and relaxation through the use of period t₂. Thus, the information obtained from the two approaches is different and complimentary.

The results of 2D IR–IR–Raman–Raman (or 2D correlation IR–IR–Raman–Raman) spectra for different t₂ values are represented in Fig. 4. As in the bending–librational case in Fig. 3, the first two IR pulses applied at time duration t₁ excite not only the fast stretching mode but also the slow translational mode due to the presence of the large μ₁₄ element. As a result, the phase of the stretching motion (with large ω₁) is twisted by the slow translational motion. Thus, the intensities of the stretching–translational peaks vary as a function of t₂. Because the vibrational motion of the translational mode decays quickly due to the large LL interaction, the phase twist disappears around t₂ = 160 fs, and a single positive peak is observed along line ω₁ = 3400 cm⁻¹ at t₂ = 300 fs. This peak is elongated in the ω₃ direction, reflecting the fact that the translational mode extends into the high-frequency region.

Finally, 2D correlation IR–IR–Raman–Raman spectra for the stretching–librational case are plotted in Fig. 5. As shown in Table II, the intensity of the μ₁₃ element is comparable to μ₁. Thus, the first two infrared pulses excite not only the stretching mode but also the librational mode, and as in the stretching–translational case, the spectral profile changes in time rapidly for small t₂ values due to the phase twist. After librational oscillation decays around t₂ = 100 fs, we observe a single positive peak at the stretching–librational coupling position, which extends in the ω₃ direction due to the large broadening of the librational peak.

We calculated third-order 2D correlation IR and fifth-order 2D correlation IR–IR–Raman–Raman spectra for stretching–bending (1–2), bending–librational (2–3), stretching–librational (1–3), and stretching–translational (1–4) modes of liquid water using the DHEOM–MLWS approach based on the multimode LL BO model. Because the peaks that arise from EAHC play a minor role in four-body response function-based spectroscopies, we can easily discuss the role of mode–mode couplings that arise from MAHC processes. As the peaks arising from MAHC in this order of spectroscopy are associated with the population transfer process, they appear after the population relaxation occurs, whereas the peaks arising from EAHC are associated with the coherent transfer and appear at t₂ = 0 before decaying rapidly.

To reproduce a 2D spectral profile, the model must capture dynamical properties of vibrational motions accurately. Thus, the scientific scenario underlying the model can be validated by calculating 2D spectra and comparing them to those obtained from experiments or accurate simulations. Taking advantage of the low computational cost and simplicity of the model, we can easily examine the mechanism to determine the 2D profile of spectrum, as briefly demonstrated in  Appendix A. Although computationally expensive, within the framework of the HEOM formalism, it is possible to describe OH–stretching modes separately as symmetric and antisymmetric modes⁶⁸ or to introduce another optically inactive mode that leads the excitation of the stretching mode to the vibrational ground state.⁸³ However, the choice of models and model parameters should be justified by comparing results obtained by advanced experimental and simulation techniques.

With the limitation of CPU power, here, we considered only the combination of the two-mode cases. To investigate a pathway of energy relaxation from the OH–stretching mode to low-frequency intramolecular modes, we should consider at least three modes, such as the stretching → bending → librational (1 → 2 → 3) modes or the stretching → bending → translational (1 → 2 → 4) modes, including the anharmonicity of the three modes described as $g̃ss′s″qsqs′qs″$⁠. We leave such investigations to future studies.

The authors are thankful to Professor Shinji Saito and Professor Keisuke Tominaga for helpful discussions. Y.T. was supported by JSPS KAKENHI (Grant No. B21H01884). H.T. was supported by JST SPRING (Grant No. JPMJSP2110).

The authors have no conflicts to disclose.

Hideaki Takahashi: Data curation (lead); Investigation (equal); Software (lead); Validation (equal); Writing – original draft (equal). Yoshitaka Tanimura: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Project administration (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (lead).

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

In the investigation of the stretching mode, it was observed experimentally that as t₂ increased, the negative (blue) 0–1–2 peak extended in the small ω₃ direction, and a peak, referred to as the hot ground state, appeared above the positive (red) 0–1–0 peak.^25–28 To predict such phenomena with the current approach, one must include various hypothetical factors that are considered important in determining the vibrational process of water, with reference to the 2D spectral profile. Here, we examine the effects of the cubic coupling strength. Thus, we made $g̃122$ 2.5 times larger than the case in Fig. 2, leaving the other parameter values unchanged. Figure 6 illustrates the calculated results. Although less pronounced than the experimentally observed spectrum, the blue 0–1–2 peak extends toward small ω₃ direction from t₂ = 0, which indicates that this peak is instantaneously generated by the coherent transfer rather than a population transfer between the two modes. A small negative peak appears for t₂ ≥ 100 fs above the red 0–1–0 peak although its peak position is much higher than the experimentally observed hot grand state peak. The appearance of this peak indicates that this peak arises from a population transfer rather than coherent transfer as the stretching → bending peak for t₂ ≥ 100 fs in the lower panels. The vibrational coherence, represented by the angle of the orange nodal line, is slightly weaker than in the case in Fig. 2.

In Fig. 7, we present third-order 2D IR spectra of the bending mode calculated for the bending → librational (2 → 3) transition case under the same conditions as in Fig. 3. This spectrum has been measured experimentally.²⁹ In this case, we observe horizontally elongated positive (red) and negative(blue) peaks along ω₃ = 1650 cm⁻¹. The nodal line is always parallel to the ω₁ axis because the noise correlation time is very short (1/γ₂ ≈ 40 fs). Note that we have also computed 2D spectra for the case of a single bending mode (not shown) and we found that the inclusion of the bending → librational (2 → 3) coupling does not change the 2D profile much. This indicates that the effect of population and coherent transfer for the bending → librational transition is minor. As a result, the peak profile changes only slightly as t₂ increases.