Discretized hierarchical equations of motion in mixed Liouville--Wigner space for two-dimensional vibrational spectroscopies of liquid water

A model of a bulk water system describing the vibrational motion of intramolecular and intermolecular modes is constructed, enabling analysis of its linear and nonlinear vibrational spectra, as well as the energy transfer processes between the vibrational modes. The model is described as a system of four interacting anharmonic oscillators nonlinearly coupled to their respective heat baths. To perform a rigorous numerical investigation of the non-Markovian and nonperturbative quantum dissipative dynamics of the model, we derive discretized hierarchical equations of motion in mixed Liouville-Wigner space (DHEOM-MLWS), with Lagrange-Hermite mesh discretization being employed in the Liouville space of the intramolecular modes and Lagrange-Hermite mesh discretization and Hermite discretization in the Wigner space of the intermolecular modes. One-dimensional infrared and Raman spectra and two-dimensional terahertz-infrared-visible and infrared-infrared-Raman spectra are computed as demonstrations of the quantum dissipative description provided by our model.


I. INTRODUCTION
Water in the condensed phase is the mother of chemistry, providing an environment that makes a variety of chemical and biological reaction processes possible. [1][2][3] One of the key challenges in investigating the physical chemistry aspects of water is the study of irreversible energy transfer processes arising from the highfrequency intramolecular modes that promote bond formation and bond breaking through complex hydrogen bonding, 4,5 and the low-frequency intermolecular modes that realize irreversible nuclear motion. [6][7][8] The interplay between these modes plays key roles in chemical reaction processes. 9 Experimentally, such phenomena have been investigated by infrared (IR) 10,11 and third-order off-resonant Raman spectroscopies. 12,13 Because water exhibits very great inhomogeneity, laser measurements based on linear responses, such as one-dimensional (1D) IR and Raman spectroscopies, yield broadened spectral peaks, as a consequence of which it is difficult to use these measurements to investigate the underlying mechanisms of complex molecular dynamics.
Thus, ultrafast two-dimensional (2D) spectroscopies have been proposed in which experiments are conducted by varying the time intervals in laser pulse trains applied to the molecules. [14][15][16][17] By plotting the nonlinear optical response as a function of the time intervals between the laser pulses, it is possible to obtain detailed information about the molecular dynamics in condensed phases. [18][19][20][21][22][23][24][25][26][27] However, analysis of these 2D signals is difficult, because they arise from complex motions described by nonlinear response functions exhibiting complicated spectral profiles that are very sensitive to the physical conditions a) Author to whom correspondence should be addressed: tanimura.yoshitaka.5w@kyoto-u.jp and the setup of the experimental system. Therefore, to enable this type of experiment to be conducted successfully, support from theoretical analysis based on a predictive model is crucial. 15,28 Molecular dynamics (MD) simulations have proved to be powerful means for such spectral analysis, but only a few approaches are able to cover the wide frequency range needed to analyze the role of energy relaxation between the intramolecular and intermolecular interactions of water molecules. [29][30][31][32] In particular, simulations of the intramolecular modes are difficult, because their motions have to be treated quantum mechanically. [33][34][35][36] Because the computational time required for 2D vibrational spectra is about 1000 times as long as that for linear spectra, full MD simulations of 2D vibrational spectroscopies have been conducted mainly for the intermolecular modes, where classical descriptions work reasonably well, [23][24][25][37][38][39][40][41][42] although full quantum MD simulations for liquid water have recently become possible. [43][44][45] As a practical approach, the Brownian oscillator (BO) model has been employed for the analysis of both linear [46][47][48] and nonlinear spectra. 15,49 In this approach, the vibrational modes representing the spectroscopic properties of interest are described as functions of molecular coordinates, while the environmental molecular motions are described using heat baths that exert thermal fluctuations and dissipation on the vibrational modes. While analysis based on 2D spectroscopy has shown that vibrational relaxation and dephasing are important mechanisms for characterizing molecular motions, 42,50 the inclusion of linear-linear (LL) and square-linear (SL) non-Markovian system-bath (SB) interactions is significant, 51 in addition to anharmonic mode-mode interactions (the LL+SL BO model). 49,52 Then, to obtain an accurate numerical solution of the model in a nonperturbative regime, various hierarchical equations of motion (HEOM) approaches have been employed. 53,54 The multimode LL+SL BO model was employed for the analysis of classical MD results of 2D THz-Raman 55 and 2D IR-Raman spectra 56 with the use of the classical hierarchical Fokker-Planck equations (CHFPE), which are the classical limit of the quantum hierarchical Fokker-Planck equations (QHFPE) in the Wigner space representation. 57 The 2D spectra obtained from the CHFPE accurately reproduce the 2D profile of the classical MD results. It should be noted that because the 2D spectral profiles are extremely sensitive to the essential features of the intermolecular and intramolecular motion, it is not possible to reproduce them satisfactorily without capturing the key features of the molecular vibrational motion that are necessary to reproduce the complex 2D profile from a simple model. Such a BO model analysis of the 2D IR-Raman spectrum predicted the presence of a cross peak representing the modemode interactions between intramolecular OH stretching motion and hydrogen-bonded (HB)-intermolecular vibrational (translational) motion, 56 which was later observed experimentally. 26,27 Although the classical MD results for 2D water spectra explain the dynamical properties of water reasonably well, the intramolecular motion of water inherently requires a quantum treatment, especially to account for the peak splitting of stretching modes in the 2D IR spectrum caused by transition between the 0-1-0 and 0-1-2 vibrational levels. 34 To aid in the construction of a reliable model, it would be helpful to have as a reference either experimental or theoretical results for 2D spectra covering the entire vibrational spectral region. However, at present, no such results are available, and to overcome this difficulty, we attempt here to extend the classically constructed model 55,56 to the quantum case and solve it accurately in the framework of open quantum dynamics theory. It is not certain that our extended model in its present form will be able to reproduce the experimental spectra if they are obtained, but modifications of the model parameters or model interactions to take account experimental and quantum simulation results of 2D spectrum should be straightforward. Once such a model has been successfully established, it could be used as a convenient tool for analyzing the energy relaxation process of water molecules, even in the quantum regime. Moreover, the model itself could be used as a quantum heat bath describing a complex water environment. Thus, the purpose of this study is to provide a basis for further analysis of complex water dynamics.
The remainder of this paper is organized as follows. In Sec. II, we present the model Hamiltonian and HEOM. We then introduce the discretized HEOM in a mixed Liouville-Wigner space in Sec. III. Although there are four major modes of water, two modes are sufficient for the calculation of a linear (1D) spectrum with a combination band peak and a 2D spectrum for detection of anharmonic mode-mode coupling. Thus, in Sec. IV, we present the calculated results for 1D IR and 1D Raman, 2D THz-IR-Visible and 2D IR-IR-Raman spectra based on two-mode calculations. Section V is devoted to concluding remarks. The computer codes for the DHEOM-MLWS used in the present calculations are provided as supplementary material.

II. MULTIMODE LL+SL BO MODEL AND HEOM
We consider a model that consists of four primary oscillator modes of liquid water representing (1) intramolecular OH stretching ("stretching"), (2) intramolecular HOH bending ("bending"), (3) hydrogen-bonded (HB)intermolecular librational ("librational"), and (4) HBintermolecular translational ("translational") motions. They are described by dimensionless vibrational coordinates q = (. . . , q s , . . . ) with s = 1, . . . , 4 indexing the four vibrational modes. Although the present model is constructed to simulate these four modes simultaneously, the calculation of the 2D spectrum takes about 1000 times longer than that of the 1D spectrum, and, in this study, we limit our analysis to various combinations of two-mode cases. The Hamiltonian of the sth mode is expressed as 55,56Ĥ where m s andp s are the mass and momentum for the sth modes with s = 1, . . . , 4. The potential of the sth mode and the interaction between the modes s and s are denoted by U s (q s ) and U ss (q s ,q s ), respectively. Each mode is independently coupled to the optically inactive vibrational modes, which are regarded as a bath system. This bath system is represented by an ensemble of harmonic oscillators. The total Hamiltonian with counterterms is then expressed as 49,52-56 where the momentum, coordinate, mass, frequency, and coupling strength of the j s th oscillator for the sth bath are given by p js , x js , m js , ω js , and α js , respectively. The bath dynamics can be characterized by the spectral distribution function (SDF), defined as J s (ω) ≡ js α 2 js δ(ω − ω js )/2m s ω js , and the inverse temperature β = 1/k B T , where k B is Boltzmann's constant and T is the thermodynamic temperature. The thermal properties of the bath are then characterized by the symmetrized correlation function (SCF) and the relaxation function (RF), expressed as C s (t) = ∞ 0 dω J s (ω) coth(β ω/2) cos(ωt) and Ψ(t) = 2 ∞ 0 dω J s (ω) cos(ωt)/ω, which are related by the fluctuation-dissipation theorem. 53,58 For the vibrational modes of water, the anharmonicity of the potential and the mode-mode interactions are weak. Thus, we assume the potential of the sth mode and the interaction between the sth and s th modes aŝ respectively, where ω s is the frequency of the s th mode, and g s 3 , g s 2 s , and g ss 2 represent the third-order anharmonicity. The dipole operator and polarizability are defined asμ respectively, where µ s and µ ss are the linear and nonlinear elements, respectively, of the dipole moment, and Π s and Π ss are those of the polarizability. Then, the vibrational modes interact through the mechanical anharmonic coupling (MAHC) described by g s 2 s and g ss 2 and the electric anharmonic coupling (EAHC) described by µ ss and Π ss . The system part of the SB interactions is denoted aŝ V s (q s ), which consists of linear-linear (LL) and squarelinear (SL) SB interactions as follows: with coupling strengths V (s) LL and V (s) SL . 33,34,[59][60][61] In this study, we assume the SDF in the Drude form as where ζ s is the SB coupling strength and γ s represents the inverse correlation time of the sth bath noise. We then have and Ψ s (t) = r (s) where ν s 0 ≡ γ s , and c (s) δ is a renormalization factor. 33,53,54 The vibrational modes of the water are then described by the reduced density matrix elements in Liouville space, ρ(t) = ρ({q s , q s }; t), with s = 1, . . . , 4. The HEOM in its Liouville-space expression for this system is then given by 53,62 Here, we have introduced the superoperator notation Ks ) consists of nonnegative integers, and e k is the unit vector of the kth element. The zeroth element ρ 0,..., 0 (t) = ρ({q s , q s }; t) corresponds to the original density element.

III. DISCRETIZED HEOM IN MIXED LIOUVILLE-WIGNER SPACE
The density matrix in the present model is a function of the eight-dimensional elements in coordinate space, and the above HEOM cannot be solved easily using currently available computer resources. To reduce computational costs, we discretize the HEOM, taking into account the characteristics of the vibrational modes of water.
Hereinafter, we distinguish the intramolecular modes and the intermolecular modes by s = 1 and 2, ands = 3 and 4, respectively, and express the reduced density matrix as ρ(t) = ρ({q s , q s }; {qs, q s }; t).

A. Lagrange-Hermite mesh discretization for intramolecular modes
In the simulation of high-frequency intramolecular modes ( ω s k B T with s = 1 and 2), the quantum nature of the system, which is commonly described in terms of energy eigenstates, plays an important role. 33,34,63 However, to accurately treat the highly excited vibrational states, which are only slightly populated, we have to deal with many low-temperature correction terms (LTCTs) in the HEOM formalism.
We then find that by using the coordinate-space representation, we do not have to treat the excited states with low population explicitly, and we can reduce the LTCT elements dramatically, while maintaining numerical accuracy. Hence, we describe the system using the Liouville operator in coordinate space and then employ the basis functions to discretize. All of the system operators are then expressed in matrix form as the elements of the basis set. Basis function approaches such as the discrete variable representation (DVR) method 64 and the Lagrange-mesh method (LMM) 65,66 have been used to discretize the wave function in coordinate space. While the DVR method has been applied to the HEOM for the investigation of the electron transfer problem of metallic surfaces, 67 we found that the Lagrange-Hermite mesh method (LHMM) is more efficient for the description of high-frequency vibrational modes.
The LHMM is a variational method utilizing Gaussian basis functions. The wave function is then described as and N s is the total number of the basis set for the sth mode, H Ns (x) is the N s th Hermite polynomial, and h Ns is the squared norm of H Ns (x), which is given by h Ns = π2 Ns N s !. We use the scaling factor defined by b s ≡ /m s ω s , where ω s is the characteristic frequency of the oscillator. The total number of mesh points in the LHMM is then N s × N s .
The reduced density matrix elements for intramolecular modes are now expressed in terms of the discretized Liouville space elements as where ρ {js,j s } (t) ≡ α js (t)α j s (t). The superoperators of the system in Liouville space for the sth mode,X(q s ), as functions of q s , are expressed aŝ where and Here, we define j l ≡ i/N s and j r ≡ 1 + (N s − 1) mod i. The superoperator for the commutator of the squared momentum operator p 2 s × can also be expressed as where Here, B. Lagrange-Hermite mesh discretization and Hermite discretization for intermolecular modes In general, the equations of motion for the density operator in Liouville space expressed as ρ({q s , q s }) are simple, since the Liouvillian is a local operator in coordinate space. Nevertheless, a phase-space (p-q space)-like description such as the Wigner representation has advantages, in particular for the investigation of intermolecular modes (s = 3 and 4), because such a distribution is a real function and is localized more or less centered at ps = 0 and qs = 0 for a system with weak anharmonicity, and this allows better usage of computer memory. Moreover, for low-frequency vibrational modes, the heat baths are considered to be at a high temperature, and a semiclassical treatment can be justified. Conversely, the density matrix element of {q s , q s } is a complex and nonlocalized function.
Thus, we introduce the Wigner distribution function (WDF), defined as fors = 3 and 4, while we fix the elements of s = 1 and 2. For the expression of system operators in Wigner space, we employ the star operator , which represents the Moyal product defined as where we have introduced the differentiation operations from the left and right, which are defined by In this expression, the commutator and anticommutator are then replaced byÔ ×ρ → O ρ−ρ O for any operator O in the Wigner representation.
For the LL+SL SB interaction, the relaxation operators in the Wigner representation for the HEOM have been evaluated and presented in Refs. 59-61, and 68. Note that, as can be seen from Eq. (22), higher-order terms can be omitted when the wavepackets in momentum space are nearly Gaussian or the anharmonicity of the potential is weak. Moreover, in the case of the lowfrequency intramolecular modes ( ωs ≤ k B T withs = 3 and 4), the quantum effects of the system and the bath are minor, and a semiclassical or even classical description of the system is reasonably accurate. Thus, here we omit higher-order Kramers-Moyal expansion terms and employ the classical Liouvillian expressed as The intramolecular and intermolecular anharmonic couplings are expressed aŝ L (s,s) In terms of the Lagrange-Hermite functions and Hermite functions, the WDF can be discretized as 69 Here, Ks and Ns are the total numbers of basis functions for momentum and coordinate spaces, respectively. The function H k (x) is the kth Hermite polynomial, and h N is the squared norm of H N (x), which is given by We use scaling factors defined as as ≡ 2ms/β and bs ≡ 1/ βmsω 2 for momentum and coordinate spaces, respectively.
Then, any function of qs for thesth mode Z (s) (qs) can be expressed as where Z (s) is,js = Z(q js )δ is,js , Hereinafter, we denote the tensor expression for the first derivative of a function Z(qs) with respect to qs as The first derivative ∂/∂qs can be expressed as where Here, we introduce the creation and annihilation operatorsbs andb † s , which act onρ {ks,js} to decrease and increase the number of ks as follows: With these operators, we can express the momentum operator as and The LHM discretization and Hermite discretization are efficient because the WDF is localized in q s and p s space for a system at high temperature (semiclassical) or in an overdamped condition, while the SB interaction operatorV (s) (q s ) is still in diagonal form as V (s) js,j s ≈ V (s) (q js )δ js,j s . Moreover, in general, the LMM is numerically stable in comparison with the finite-difference method. 65 Such features allow us to dramatically reduce the computational time required to integrate the HEOM, in particular as the size of a system increases. Note that for an unbounded system and a rotationally invariant system, a Lagrange-Fourier mesh method 70 and a discrete Wigner function method 71 are respectively more efficient.

C. DHEOM-MLWS
The elements of the discretized reduced density matrix are now expressed asρ {ksjs} {js,j s } (t) ≡ ρ {js,j s } (t)ρ {ks,js} (t). The reduced density operator is then expressed in tensor form asρ The discretized HEOM are then expressed as where and L (s,s) I Here, intra−inter denotes summation with respect to coupling between intramolecular and intermolecular mode and, to keep the notation simple, the unit matrix is not denoted, The expressions for other auxiliary operators are presented in Appendix A

IV. NUMERICAL DEMONSTRATIONS
We now report the results of our numerical computations of the DHEOM-MLWS. We employed the parameter values of the multimode LL+SL BO model chosen to reproduce 2D IR-Raman spectra obtained from classical MD simulations 56 with the use of the POLI2VS force fields, 29 which possess the essential capability of simulating both IR and Raman spectra. We then modified the anharmonicity and bath parameters of the intramolecular modes to fit an experimentally obtained IR spectrum 11 and Raman spectra 12,13 that are consistent with the 1D spectra obtained from quantum MD simulations with the POLI2VS force fields. 43 That is, the anharmonicity of the intramolecular modes were modified by factors of 10 and 5 for the OH stretching mode (s = 1) and HOH bending mode (s = 2), respectively, while the bath coupling strength and inverse noise correlation time (ζ s , γ s ) were modified by factors of (1.5, 2) for s = 1 and (0.7, 1.1) for s = 2 from the classical values presented in Ref. 56. Moreover, we enhanced the optical propertiesμ ss and Π ss for the 1-4 and 2-3 mode-mode couplings by factors of 2 and about 100, respectively, to reproduce their overtone peaks. The anharmonicity of the potential and the mode-mode coupling strength are listed in Tables I  and II. Here, we employ the normalized parameters to compare the effect of anharmonicity with respect to the potential for each mode and mode-mode coupling. The bath temperature was set to T = 300 K (β ω 0 ≈ 19.2), with fundamental frequency ω 0 = 4000 cm −1 , which was chosen as a frequency close to the OH stretching mode.
A. Linear response: 1D IR and 1D Raman spectra In Figs. 1(a) and 1(b), we present the calculated 1D IR and 1D Raman spectra, respectively, from CHFPE and DHEOM-MLWS using the same parameter values of the present model. We obtained these spectra by combining the results from a single-mode model with s = 1 and 2 ands = 3 and 4 and those from a two-mode model with the bending mode (s = 2) and stretching mode (s = 4), because the calculations for four modes are computationally expensive and because the effects of mode-mode coupling are important only for the combination band among s = 2 ands = 4.
In both the IR and Raman cases, the classical and quantum results agree for the low-frequency intermolecular modes. However, for the high-frequency intramolecular modes, the stretching and bending peaks in the classical case are blue-shifted owing to the quantum effects arise from the anharmonicity. Our DHEOM-MLWS results reproduce the weak bending-librational combination band at 2130 cm −1 accurately, while the classical MD results 29 and the classical BO results, 56 including the present CHFPE results, underestimate the peak position and peak intensity.
Note that the present LL+SL BO model was constructed based on classical 2D IR-Raman simulations, but the force field used in the simulation (POLI2VS) was developed for quantum MD simulations. Because of this, although some modifications of the parameter values for the intramolecular modes are necessary, such an LL+SL BO model can predict a reasonably accurate vibrational spectrum when we conduct quantum HEOM calculations, as is also the case with quantum MD calculations using the POLI2VS force fields. 43 This indicates the possibility of constructing a quantum BO model from 2D spectra obtained from first-principles classical MD simulations, in which the nuclear motion of the molecules is classical. 32 To elucidate how the effect of quantum dissipative dynamics is manifested in mode-mode coupling peaks, we next present numerical results for the 2D THz-IR-visible (2D TIV) spectrum 26,27 and 2D IR-IR-Raman (2D IIR) spectrum (the observable part of the 2D TIV spectrum is equivalent to part of the 2D IIR spectrum presented in Ref. 56). Note that in the 2D TIV and 2D IIR spectra expressed in terms of the three-body correlation functions of optical observables, the nondiagonal spectral peaks are not necessary to represent mode-mode coupling peaks as in the case of the third-order 2D IR spectrum expressed in terms of the four-body correlation functions of the dipole moment, 15-17 because the signal from the EAHC appears at a similar location to that from the MAHC. 25,74 Such 2D experiments have been conducted 26,27 on the basis of classical MD simulations, 56 but the results are not in good agreement with theoretical predictions, partly because of the classical description of the system. Hence, here we calculate and compare the 2D spectra for the classical and quantum cases using the same BO model, although the difference from the MD results may also be due to the limitation of the MD description. Note that although here we consider the 2D TIV and 2D IIR cases, 2D spectra computed from different pulse configurations such as 2D IR-Raman-IR and 2D Raman-IR-IR spectra exhibit similar profiles, because the difference between the IR and Raman spectra determined from the EAHC is minor in our calculations based on the BO model.
In Figs. 2(a) and 2(b), we compare 2D TIV results under the same conditions, calculated for the quantum (DHEOM-MLWS) and classical (CHFPE) cases, respectively. Characteristic features of anharmonicity and nonlinear polarizability on such 2D spectral profiles in a single-mode case and a two-mode case described by the BO model were elucidated in Ref. 56. From that analysis, the negative peak [at (ω 1 , ω 2 ) = (150 cm −1 , 3700 cm −1 ) in the quantum case] arises only from the MAHC between the stretching-translational modes, whereas the positive peak [around (ω 1 , ω 2 ) = (150 cm −1 , 3400 cm −1 ) in the quantum case] arises from contributions from MAHC and EAHC. This can be easily confirmed by comparing the same calculation without the nonlinear polarizability (Π s,s = 0). Figure 3 reveals a negative peak and a positive peak whose node lines are centered at the resonant frequency, whereas we observe only a positive peak at the resonant frequency in the pure EAHC case (see Ref. 56). As this fictitious model analysis has demonstrated, we can easily identify the key dynamics of a liquid water system that determine the 2D spectral profiles obtained from experiments and complex MD simulations. This is because, to reproduce a complex 2D spectral profile accurately from a simple model, the model must capture the dynamical properties of the system correctly. 49,53,54 Compared with previous classical MD and BO model calculations, in the present calculation we observe peaks near (ω 1 , ω 2 ) = (0 cm −1 , 3700 cm −1 ) and (ω 1 , ω 2 ) = (0 cm −1 , 3400 cm −1 ) in both the classical and quantum cases. Such peaks arise for a vibrational system strongly coupled to an Ohmic bath, as has been demonstrated from the analytical expressions for 1D and 2D spectra. 51,75 Although the BO model that we have employed here is similar to that used to analyze classical MD results, such low-frequency peaks could not be observed in the previous studies, because their MD simulation period in the t 1 direction was too short. 56 We next demonstrate the description of the stretching- bending (1-2) modes. In Fig. 4, we depict the 2D IIR spectrum for the 1-2 modes. In this figure, there are both positive and negative peaks in the stretching, bending and their cross peak positions, indicating that these arise from the MAHC. This is because the EAHC contribution in the 1-2 modes is small, as indicated in Table II.
Here we observe the bending peak around (ω 1 , ω 2 ) = (1600 cm −1 , 1600 cm −1 ), but this peak overlaps with the bending-librational EHAC peak, as shown below, and cannot be identified. In Fig. 5, we depict the 2D IIR spectrum calculated for the bending-librational (2-3) modes. As in the previ-ous 1-4 case, because we perform an accurate quantum mechanical evaluation of the bending mode, the 2-3 coupling peaks around (ω 1 , ω 2 ) = (600 cm −1 , 1600 cm −1 ) and (ω 1 , ω 2 ) = (1600 cm −1 , 600 cm −1 ) are blue-shifted in comparison with the classical results. 56 The peak intensities of the 2-3 coupling peaks are much larger than the previous result 56 because we enhance the MAHC to reproduce the 2-3 combination band in the 1D IR spectra. The peaks that appear positively and negatively around ω 2 = 1600cm −1 across ω 1 = 1600cm −1 are caused by the EAHC of the 2-3 modes. Because we set Π 2,3 much larger value than classical case, 56 the bending peak displayed in Fig. 4 is completely covered by these EAHC peaks.

V. CONCLUSION
We have developed a model to analyze 1D and 2D vibrational spectra for both intramolecular and intermolecular vibrational modes involving all of their mode-mode interactions, taking into account the effects of energy relaxation and vibrational dephasing. To compute 2D signals from the model system, it is important to adopt a quantum-mechanically consistent treatment of the system and bath, in particular for the intramolecular modes, because the quantum entanglement between system and bath plays an essential role. Thus, we adapted the HEOM formalism here, enabling us to perform rigorous numerical calculations of linear and nonlinear spectra. Because integrating the HEOM for a multimode system is computationally expensive, we developed the DHEOM-MLWS approach to maintain the accuracy of the numerical calculation.
The description of the multimode LL+SL BO model with the use of the DHEOM-MLWS was investigated by calculating 1D and 2D spectra. From calculations of linear and nonlinear spectra, we obtained accurate predictions of the positions of the frequency stretching and bending peaks, for which the classical results are redshifted.
The parameter values of our model were first chosen by solving the classical HEOM to fit the classical MD results for 1D and 2D spectra, 56 with the anharmonicity of the potentials being modified using the experimentally obtained 1D spectrum. We found that by using the POLI2VS force field in classical MD simulations, we could obtain a reasonable parameter value set for quantum HEOM simulation. This is because the POLI2VS force field was developed for quantum MD simulations, 43 and the complexity of molecular interactions, which is important in describing water spectra, is not directly related to the issue of quantum effects. This also indicates that, even using first-principles MD results, 32 in which the nuclear motion of the molecules is classical, we may construct a quantum SB model that includes complex anharmonic and bath interactions.
To reproduce a 2D spectral profile, the model must capture the dynamical properties of the vibrational mo-tions correctly. Taking advantage of the low computational cost and simplicity of the model, we can easily examine, for example, the effects of higher-order anharmonicity on 2D spectra. The ability of the quantum mechanical model to calculate 2D spectra provides the possibility of directly analyzing experimentally obtained spectra. Once the model has been fully established, we can use it to investigate energy and excitation transfer processes in liquid water. Moreover, we can employ the LL+SL BO model as a heat bath to study the spectra and energy relaxation of liquids containing ions.
Extensions of the present model, for example, to describe symmetric and antisymmetric OH stretching modes separately and to employ the SDF to include the effects of optically inactive modes, are also possible. 76 The present model with the HEOM approach provides a platform for analyzing novel experimental and simulation results. We leave such extensions to future studies, depending on progress in experimental and simulation techniques.

SUPPLEMENTARY MATERIAL
See the supplementary material for the computer codes for the DHEOM-MLWS used in the present calculations. and For intermolecular modess = 3 and 4, we have and In the classical limit → 0, Eqs. (A4)-(A7) reduce to Ξ (s) ≡ 0, and Λ (s) 2 = 0.
Appendix B: Truncated Padé spectral decomposition When quantum effects described by an SB model becomes important, we have to take into account many LTCTs involved in the HEOM formalism, which makes the integration of the HEOM computationally very expensive. Thus, to reduce the number of LTCTs, a Padé spectral decomposition (PSD) scheme has been developed. 77 We can further reduce the number of hierarchical elements by incorporating into the PSD the balanced truncation method (BTM), which was originally developed as a model order reduction (MOR). 78 Here, we adapt the algorithm developed in Ref. 79 for the LTCTs and demonstrate the efficiency of the truncated PSD (TPSD) method.
For a desired accuracy > 0, we consider the condition where C PSD (t) is the SCF described by PSD, and C TPSD (t) is the SCF described by PSD with BTM. In Fig. 6 we depict the time evolution of the groundstate population of the OH stretching mode (s = 1) described using the LL+SL BO model with parameters values as listed in Table I. We then employ the energyeigenstate representation and integrate the HEOM with PSD and with TPSD. As depicted in Fig. 6, the results calculated with TPSD converge faster than those without TPSD, while the number of hierarchical terms is fewer in the case without PSD. The improvement becomes significant for a system strongly coupled to a heat bath at low temperatures. This approach is particularly beneficial when we deal with multiple heat baths. Note that this method can also be combined with the NZ2 truncation method. 80 Appendix C: Linear and nonlinear spectra Because the HEOM formalism is able to take accurate account of the quantum entanglement between system and bath, it is possible to calculate linear nonlinear response functions. To compute a nonlinear spectrum in the HEOM approach, we express the response functions in terms of the time-propagation operator. For example, a 1D spectrum defined by first-order response functions is expressed in terms of the two-body correlation function as where we have employed the hyperoperator × defined aŝ A ×ρ ≡ [Â,ρ] for the Liouville space representation and A ×ρ ≡Â Ŵ −Ŵ Â for the Wigner space representation, G(t) is the Green's function of the system Hamiltonian without a laser interaction, andρ eq is the equilibrium state. For 1D IR or 1D Raman calculations, we chosê A =μ orÂ =Π given in Eq. (5) or (6). Accordingly, the 2D spectrum defined by the secondand third-order response functions is expressed in terms of the three-and four-body correlation functions of optical observables as follows: 53,81 R (2) (t 2 , t 1 ) = i The above equations represent the time evolution of the system under laser excitation. For example, Eq. (C2) can be interpreted as follows. The system is initially in the equilibrium stateρ eq and is then modified as a result of the first laser pulse via the dipole interaction bŷ C. It then propagates for time t 1 under G(t 1 ). The system is next excited through the second laser pulse bŷ B and propagates for time t 2 under G(t 2 ). Finally, the expectation value of the polarizability at t 1 + t 2 is generated through the laser pulses byÂ. 53,54 The 2D THz-IR-visible signal can be computed from R  From the second-and third-order response functions, the 2D TIV and 2D IR spectra, for example, are evaluated as TIV (t 2 , t 1 ) sin(ω 1 t 1 ) sin(ω 2 t 2 ) (C4) and I (3) IR (t 3 , t 2 , t 1 ) sin(ω 1 t 1 ) sin(ω 3 t 3 ). (C5) This 2D sine-Fourier representation is more intuitive than the real part of the 2D-Fourier representation, since it can extract only absorptive components.

Appendix D: Estimation of signal parameters via rotational invariance techniques
Calculating 2D spectra is computationally expensive, because we have to repeat the dynamics calculations for different t 1 and t 2 , or t 1 , t 2 , and t 3 . Estimation of signal parameters via rotational invariance techniques (ES-PRIT) leads to a dramatic reduction in the computational cost, because it allows us to find an optimal explored solution of a targeting signal as a linear combination of complex exponentials.
Thorough the use of ESPRIT, the second-order response function as a function of t 1 is, for example, expressed as where a i (t 2 ) and b i (t 2 ) are complex functions of t 2 that are chosen to optimize R(t 1 , t 2 ). Although the real part of b i (t 2 ) must be positive to avoid divergence in the t 1 direction, we can eliminate this limitation with the use of a Fourier-Laplace transform. The 2D Fourier transform of Eq. (D1) is then expressed as where F represents the discrete Fourier transform (DFT) on t 2 .
Using ESPRIT, we obtained the 2D IIR spectrum for the stretching-bending modes shown in Fig. 4. While conventional calculations require 192 sample points in the t 1 direction, we can reduce the necessary calculations four times (48 sampling points) with the use of ESPRIT, giving almost identical results.
Note that while Prony's method has been incorporated into the HEOM formalism to eliminate instabilities arising from the discrete-bath HEOM 82 and to apply the time-domain Prony fitting decomposition (t-PFD) scheme as an efficient description of SDF,, 83 the same method can also be used to evaluate a 2D spectral profile efficiently. However, we have found ESPRIT to be more convenient and stable, in particular when the signal is mixed with noise.