Notes on simulating two-dimensional Raman and terahertz-Raman signals with a full molecular dynamics simulation approach

Recent developments in two-dimensional (2D) THz-Raman and 2D Raman spectroscopies have created the possibility for quantitatively investigating the role of many dynamic and structural aspects of the molecular system. We explain the significant points for properly simulating 2D vibrational spectroscopic studies of intermolecular modes using the full molecular dynamics approach, in particular, regarding the system size, the treatment of the thermostat, and inclusion of an Ewald summation for the induced polarizability. Moreover, using the simulation results for water employing various polarization functions, we elucidate the roles of permanent and induced optical properties in determining the 2D profiles of the signal.


I. INTRODUCTION
Intermolecular vibrations of molecular liquids and biological material in the frequency range 0-700 cm À1 play an essential role in many chemical and biological processes, because they promote reactive dynamics via interactions through intramolecular modes and because they are active at room temperature. While one-dimensional (1D) IR, THz, and Raman spectroscopic approaches are versatile tools for investigating collective intermolecular motion, 1-6 using such approaches alone, it is not clear whether the vibrational modes that they investigate are mutually coupled, and whether the width of the vibrational mode peaks that they measure are from an inhomogeneous origin or a homogeneous origin, because these contributions are usually broadened and overlap in the 1D spectra. 7 To elucidate these points, in 1993, the fifthorder two-dimensional (2D) Raman spectroscopy approach was proposed to distinguish the contributions from inhomogeneous broadenings using three sets of Raman pulses. 8 This initiated the development of 2D infrared (IR) spectroscopy for intramolecular vibrational modes. [9][10][11] While the possibility for applying 2D Raman spectroscopy to the study of anharmonicity, [12][13][14][15][16] mode-mode coupling mechanisms, [17][18][19][20][21] and dephasing processes [22][23][24][25][26][27][28][29][30] of intermolecular modes has been explored theoretically, due to the unforeseen cascading effect of light emissions, [31][32][33] experimental signals have been obtained only for CS 2 , 34-37 benzene, 38 and formamide liquids. 39 However, recently developed single-beam spectrally controlled 2D Raman spectroscopy method overcomes the difficulty of cascading effects 40 and creates a new possibility for measuring intra-molecular interactions of liquids by means of 2D Raman spectroscopy.
While the 2D Raman signals of liquid water has not yet been observed, 2D THz-Raman (or 2D Raman-THz) spectroscopy was introduced in 2012 to measure the water signal. In this method, the cascading effects are suppressed using two THz pulses and one set of Raman a) Electronic mail: h.ito@kuchem.kyoto-u.ac.jp b) Electronic mail: ju8879@kuchem.kyoto-u.ac.jp c) Electronic mail: tanimura@kuchem.kyoto-u.ac.jp 2329-7778/2015/2(5)/054102/21 V C Author(s) 2015 2, 054102-1 pulses. [41][42][43][44][45][46][47] Because 2D THz-Raman utilizes the dipole moment, in addition to the polarizability, the applicability of this spectroscopy is different from that of 2D Raman spectroscopy. 45 Moreover, the information obtained from 2D Raman and 2D THz-Raman spectroscopies can be used in a complementary manner to investigate the fundamental nature of intermolecular interactions. 47 Theoretically, the linear absorption (1D IR or 1D THz spectroscopy) signal and the thirdorder Raman spectroscopy (1D Raman spectroscopy) signal are obtained from the linear response functions of the optical observables. There, the main contribution to a signal arises from harmonic vibrational motion as a function of molecular polarizability or the dipole moment. Contrastingly, in 2D Raman and 2D THz-Raman spectroscopic approaches, the threebody nonlinear response function of the molecular polarizability and/or the dipole moment is measured to monitor molecular motion. Because the complex profiles of such 2D signals depend on many dynamic and structural aspects of the molecular system, full molecular dynamics (MD) simulations for the nonlinear response function play important roles in the design of 2D spectroscopy experiments and the analysis of their results, particularly in regard to intermolecular vibrations. In the 2D case, the anharmonicity of potentials and the nonlinearities of the polarizability contribute significantly to the signals. 30,45,47 For these reasons, a reliable potential model and an accurate polarizability function are necessary to accurately predict quantitative features of 2D signals. Moreover, the theoretical details of the MD simulations, in particular, regarding the system size, the treatment of the thermostat, and the proper inclusion of an Ewald summation for the induced polarizability, must be carefully examined, because the 2D profile of a signal is also very sensitive to numerical errors. The purpose of this paper is to elucidate the significant points regarding the calculation of 2D Raman and THz-Raman signals using full MD approach while minimizing computational costs.
This paper is organized as follows. In Sec. II, we explain the methodology for simulating 2D Raman and 2D THz-Raman signals with full MD simulations. In Sec. III, we investigate the size dependence of simulations, the artifacts of the thermostat, sensitivity of the signals to the choice of the force field, and the importance of employing an Ewald sum of the induced polarizability function in calculating 2D Raman signals. In Sec. IV, we explain the role of the force field and the sensitivity of the signals to the choice of the polarizability functions, particularly in calculations of 2D THz-Raman signals for liquid water. Section V is devoted to concluding remarks.

II. FULL MD APPROACH FOR 2D RAMAN AND 2D THz-RAMAN SPECTROSCOPY
The optical observable in 2D spectroscopy is expressed in terms of the three-body response function as 8 whereÂ;B, andĈ can be the dipole moment, l, or the polarizability, P, of the molecules expressed as functions of the molecular position coordinates, q. The classical mechanical expression for the response functions can be obtained by replacing the commutator and operators with the Poisson bracket and c-number observables as The response function in the classical limit is then expressed as A. Equilibrium MD approach Three full MD approaches have been developed to this time for evaluating the above response function. The first approach is based on equilibrium MD simulations. We evaluate the outer Poison bracket in terms of the time derivative of the observable, employing the relation, fCðÀt 1 Þ; e ÀbH 0 ðp;qÞ g PB ¼ Àb _ CðÀt 1 Þe ÀbH 0 ðp;qÞ , with the molecular Hamiltonian H 0 ðp; qÞ, as follows: Here, b is the inverse temperature divided by the Boltzmann constant, and _ CðÀt 1 Þ is the time derivative of the observable CðÀt 1 Þ, defined as _ CðÀt 1 Þ ¼ dCðtÞ=dtj t¼Àt 1 . Then, we obtain the response function in terms of the stability matrix as We can calculate this response function using the time and coordinate derivatives of lðtÞ or PðtÞ evaluated from the molecular trajectories pðtÞ and qðtÞ that are obtained from the equilibrium MD simulations. [48][49][50] While the equilibrium MD approach is convenient for analyzing the 2D profile of a signal as the contribution of the anharmonicity and the nonlinear polarizability, 49,50 this requires a great deal of CPU time and memory due to the computational intensiveness of the treatment of the equation of motion used to calculate the Jacobian (or stability block) matrix, @qðt 2 Þ=@pð0Þ. Moreover, because of the stability matrix, the convergence of the signal becomes very slow, particularly for a large molecular system. 51,52 B. Non-equilibrium finite field MD approach The second approach that we consider, the non-equilibrium MD (NEMD) approach, was developed to evaluate double Poisson brackets from non-equilibrium trajectories in the case that there exist multiple external perturbations. 53,54 In this approach, the outer Poison bracket in Eq. (3) is evaluated from the relation fDðt 2 Þ; CðÀt 1 Þg PB ¼ ÀðD þCðÀt 1 Þ ðt 2 Þ À D ÀCðÀt 1 Þ ðt 2 ÞÞ=2FdðtÞ, where D 6CðÀt 1 Þ ðt 2 Þ is the expectation value corresponding to Dðt 2 Þ fAðt 2 Þ; Bð0Þg PB calculated from the trajectories subjected to weak perturbations 6ðÀFdðtÞCðÀt 1 ÞÞ, with the external field 6F acting on CðÀt 1 Þ. Then, the Poisson bracket Dðt 2 Þ is also evaluated from Dðt 2 Þ ¼ ÀðA þBð0Þ ðt 2 Þ À A ÀBð0Þ ðt 2 ÞÞ=2FdðtÞ, where A 6Bð0Þ ðt 2 Þ is defined in the same manner as D 6CðÀt 1 Þ ðt 2 Þ. This approach does not involve a convergence problem of the stability matrix. However, it is too computationally intensive, because we have to repeat the calculations four times in order to account for the configurations of the external perturbations.

C. Equilibrium-non-equilibrium hybrid MD approach
The third approach that we consider, the equilibrium-non-equilibrium hybrid approach, was developed to take advantage of the merits of both the equilibrium and non-equilibrium approaches. 55 In the hybrid approach, the outer Poison bracket is evaluated in terms of the time derivative of the observable, while the inter-Poison bracket is evaluated using the NEMD approach. Then, the response function is expressed as In this approach, we first obtain the time derivative of the dipole moment _ l eq ðÀt 1 Þ or polarizability _ P eq ðÀt 1 Þ from the equilibrium trajectories at time t ¼ Àt 1 . Next, we evaluate the dipole moments l 6lð0Þ ðt 2 Þ and l 6Pð0Þ ðt 2 Þ or the polarizabilities P 6lð0Þ ðt 2 Þ and P 6Pð0Þ ðt 2 Þ at time t ¼ t 2 from the non-equilibrium trajectories, which are generated by a perturbation at time t ¼ 0, either 6ðÀlð0ÞE 1 dðtÞÞ or 6ðÀE 1 Pð0ÞE 2 dðtÞ=2Þ, resulting from the external electric field of the jth pulse E j acting on the dipole moment lð0Þ or the polarizability Pð0Þ.

III. SIMULATING 2D RAMAN SIGNALS
Single-beam spectrally controlled 2D Raman spectroscopy, which was developed recently, allows us to obtain clear 2D spectra of molecular liquids in the THz region by suppressing cascading effects. 40 This motivated us to carry out MD simulations of 2D Raman spectroscopy, in addition to 2D THz-Raman spectroscopy. Here, we discuss the important points in carrying out full MD simulations of 2D Raman spectroscopy, specifically with regard to the conditions of the simulations and the verification of force fields and polarizability functions.
Unless otherwise noted, the MD simulations were carried out as follows. The 2D Raman signals were obtained using the hybrid MD approach described in Sec. II C. Such a signal is given by _ P eq ðÀt 1 ÞðP þPð0Þ ðt 2 Þ À P ÀPð0Þ ðt 2 ÞÞ E ; where Dt is the time step used in integrating the equations of motion and P eq ðÀt 1 Þ is the polarizability with the equilibrium trajectories at time Àt 1 . For comparison, we also calculated the 1D Raman spectra using the expression I Raman ðxÞ ¼ = In the cases of formamide and carbon disulfide (CS 2 ), we introduced a harmonic quantum correction factor and evaluated the signal from xI Raman ðxÞ. 4 The MD simulations were carried out with 108 molecules for the CS 2 liquid and 64 molecules for the other liquids. The simulations were carried out in a cubic box with periodic boundary conditions. The interaction potentials were cut off smoothly at a distance equal to half the box length using a switching function, and the long-range Coulomb interactions were calculated with the Ewald summation. 56 The intramolecular geometries were kept rigid throughout the simulations, using the RATTLE algorithm. The equations of motion were integrated using the velocity-Verlet algorithm with Dt ¼ 2.5 fs for water and Dt ¼ 5.0 fs for the other molecular liquids. The system volume and total energy were fixed (in accordance with a microcanonical simulation) after the completion of the isothermal simulations that were carried out in equilibration with a Nose-Hoover chain thermostat. 57 The volume was chosen to reproduce experimental densities: 0.997 g/cm 3 for water, 1.120 g/ cm 3 for formamide, 0.815 g/cm 3 for formaldehyde, and 1.270 g/cm 3 for CS 2 . The temperature was set to 300 K for water and formamide, 255 K for formaldehyde, and 270 K for the CS 2 . To estimate the polarizability of each liquid, we employed the dipole-induced-dipole (DID) model with an Ewald summation using the permanent molecular polarizability, which was determined from the Huiszoon polarizability for water, 58 the atomic polarizability for formamide and formaldehyde, 59 and experimental data for CS 2 . 60 In the hybrid MD approach, the NEMD part of the calculation was carried out with the Raman laser fields E ¼ 4.0 V/Å for water, E ¼ 1.0 V/Å for formamide and formaldehyde, and E ¼ 2.0 V/Å for CS 2 . These signals were obtained by averaging over 10 6 initial configurations.
A. Size dependence of the simulations Because our objective is to describe fast intermolecular modes, which arise from shortrange intermolecular interactions, it is not necessary to carry out large-scale simulations with many molecules. In general, 64 molecules are sufficient to obtain a reliable signal. To illustrate this point, we employed 2D Raman signals for water calculated with various system sizes. In these computations, the interactions between the molecules were modeled using the TIP4P/2005 potential, 61 and the full-order (see Appendix A 1) DID models were employed to evaluate the polarizability using an Ewald summation (see Sec. III D). 56 The computational results are presented in Fig. 1 with (a) 32, (b) 64, (c) 108, and (d) 216 water molecules. No clear size dependence is observed for water molecules. Although there is a possibility that the size dependence varies among the types of molecules, the present results suggest that a system size of 32 molecules is sufficient to elucidate the qualitative properties of 2D signals, at least for liquid water.

B. Microcanonical (NVE) and canonical (NVT) simulations
In MD simulations, a thermostat is often employed to maintain the system temperature (canonical simulation). 57 Because 2D spectroscopy is very sensitive to molecular motion, however, the thermostat may alter the 2D profiles of the signal. To demonstrate this point, we calculated and compared these 2D Raman signals for CS 2 liquid by using the Nose-Hoover chain thermostat, i.e., a canonical ensemble (NVT), and compared these results with those obtained from microcanonical (NVE) simulations. These simulations were carried out with the Lennard-Jones (LJ) model 62 using the polarizability described by the full-order DID model (see Appendix A 1). To evaluate the induced polarizability, we employed the Ewald summation. Figures 2(a)-2(c) illustrate the effects of the thermostat in the 1D and 2D Raman spectra. The broadened peaks near 50 cm À1 in Fig. 2(a) reflect the presence of intermolecular vibrations. While the 1D Raman signals in Fig. 2(a) are similar in both cases, there is a difference along the t 2 axis in the thermostatted case for the 2D Raman signals in Figs. 2(b) and (c). The elongation of the negative signal along the t 2 direction arises from the anharmonicity of the potential. 45,47,49 This indicates that the Nose-Hoover chain thermostat acts as an undesirable source of anharmonicity for the molecular dynamics. Because the Nose-Hoover chain thermostat was formulated to study thermal equilibrium states, the dynamics obtained using this thermostat are not necessarily accurate. The sensitivity of the 2D measurements may reveal such inaccuracy. Thus, it is most prudent to simulate the 2D spectrum using the NVE ensemble without a thermostat.

C. Choosing the force field
As was shown in a formamide case, the profiles of 2D Raman signals are very sensitive to the choice of the force field. 39 Here, we demonstrate this point using two types of force fields for CS 2 while keeping the molecular polarizability function fixed. The first model consists of Lennard-Jones interactions only (LJ model), in order to simulate intermolecular vibrational spectra, 62 whereas the second model also includes Coulomb interactions (LJ þ Coulomb model), in order to account for the structural information obtained from neutron and X-ray scattering experiments. 63 While almost all full MD simulations for 2D Raman spectroscopy of CS 2 have been carried out using the first model, 49,53-55 here we examine the validity of results obtained using a more reliable potential that includes the Coulomb interactions. In both cases, full-order DID models were employed to evaluate the polarizability using the Ewald summation.
While the LJ and LJ þ Coulomb results are similar in the 1D case considered in Fig. 2(a), we observe a difference in the 2D Raman case considered in Figs. 2(b) and 2(d). Because the anharmonicity of the potential is assumed to be large in the LJ þ Coulomb case, the negative signal along the t 2 axis does not decay even at 1 ps, while the elongation vanishes within 500 fs in the LJ case. It should be noted that the LJ þ Coulomb potential is more reliable than the LJ potential, because the LJ þ Coulomb potential reproduces the structure of the CS 2 liquid more accurately. 63 The present results indicate that an accurate force field is necessary in order to obtain the correct profile of a 2D Raman signal even in the case of a nonpolar molecular system, due to the sensitivity of the nonlinear response function. For this reason, the accuracy of the potential should be evaluated using the experimental results.
D. Ewald summation of the induced polarizability: Long-range effects By nature, the 2D Raman profile is very sensitive to the functional form of the polarizability. For this reason, the signal profile is significantly affected by the cutoff of the polarizability function that is introduced to reduce the computational intensiveness of the simulation. To demonstrate this point, we calculated the 2D Raman signals for water, formamide, and formaldehyde described by the TIP4P/2005 potential, 61 the modified T potential, 64 and the 4-site model, 65 respectively, with and without the Ewald summation for the induced polarizability. The DID model was employed to evaluate the total polarizability in all of the simulations. In the induced polarizability, we include the dipole-dipole interaction term defined by Eq. (A4), which includes the contributions from the infinite periodic images while taking into account the changes in the distances between all of the molecular pairs. Without using the Ewald summation, however, the long-range effect is ignored. Figure 3 displays the computational results of the 2D Raman signals for (i) water, (ii) formamide, and (iii) formaldehyde. The 2D Raman signals displayed in Fig. 3(a) were obtained using the Ewald summation, while the 2D Raman signals displayed in Fig. 3(b) were calculated by cutting off the dipole-dipole interaction at half of the box length without using the Ewald summation. All of the simulations whose results are presented here were carried out with 64 molecules. Note that we examined the size dependence of the results by considering systems with 64, 108, and 216 molecules without using the Ewald summation. We found that there were only minor differences among the signals obtained from these small systems. 45,55 While the 2D Raman signals were calculated from the same trajectories, they are significantly different for formamide and formaldehyde. By contrast, the 2D Raman signals for water exhibit only a slight difference near t 1 ¼ t 2 ¼ 0 fs, with the negative signals being weaker when computed using the Ewald summation. The cause of these differences is analyzed in Appendix B.
The results discussed above indicate that, although the effects depend on the type of molecule, the long-range contributions to the polarizability must be computed from the dipole-dipole interactions using the Ewald summation in the 2D Raman case.

IV. SIMULATING 2D THZ-RAMAN SIGNALS
Using the hybrid approach discussed in Sec. II C, the response functions for the 2D Raman-THz-THz (RTT), THz-Raman-THz (TRT), and THz-THz-Raman (TTR) spectroscopic approaches are calculated as respectively. Here, in addition to the 1D Raman spectrum defined by Eq. (8), we also evaluated 1D THz spectrum I THz ðxÞ expressed as I THz ðxÞ ¼ x= where the pre-factor is a harmonic quantum correction factor that must be applied to the classical calculations. 66 It should be noted that the types of information that we can obtain from the 2D Raman signal and each of the three THz-Raman signals are different, due to the role of the nonlinear polarizability. We can separate the inhomogeneous and anharmonic contributions to the signal more clearly in the cases of 2D THz-Raman spectroscopy than in the case of 2D Raman spectroscopy, because the inhomogeneous contribution arises from the TRT pulse configuration, while the anharmonic contribution arises from the RTT configuration. 45,47 The conditions of the 2D Raman simulation discussed in Sec. III, regarding the system size, usage of a thermostat, and the sensitivity to the force field, can be applied to the present 2D THz-Raman case as well. However, the 2D THz-Raman signals are insensitive to the use of the Ewald summation for the induced polarizability, because the 2D THz-Raman response functions account for only one Raman process that is sensitive to the induced polarizability

054102-7
Ito, Jo, and Tanimura Struct. Dyn. 2, 054102 (2015) contribution. 47 In this section, we explain the significant points involved in the simulation of 2D THz-Raman signals specifically for the case of water, focusing on the polarizability and potential model, because presently, water is of primary interest in experimental investigations employing 2D THz-Raman spectroscopy. 44 Unless otherwise noted, the MD simulations considered here were carried out for water under the same conditions as in the 2D Raman case. The NEMD part of the calculation involved in the TRT response was carried out with a Raman pulse intensity of E ¼ 4.0 V/Å , whereas those involved in the RTT and TTR responses were carried out with a THz pulse intensity of E ¼ 0.1 V/Å . The 2D signals were obtained by averaging over 10 6 initial configurations.

A. Choosing the permanent polarizability
Here, we first demonstrate the sensitivity of the 2D THz-Raman signals to the choice of the permanent polarizability. We compared the 2D profiles obtained using the Huiszoon permanent polarizability 58 and the coupled-cluster single-and double-excitation (CCSD) permanent polarizability (see Appendix A 4), while the TIP4P/2005 potential 61 was used for the force field in both cases. The total dipole moment and polarizability were calculated using the DID model with the Ewald summation.
The parameter values for the Huiszoon polarizability were set to a xx ¼ 1.626 Å 3 , a yy ¼ 1.495 Å 3 , and a zz ¼ 1.286 Å 3 , while those for the CCSD polarizability were set to a xx ¼ 1.442 Å 3 , a yy ¼ 1.375 Å 3 , and a zz ¼ 1.321 Å 3 . Here, the X axis is defined as that connecting the hydrogen atoms, the Y axis lies along the bisector of the H-O-H angle, and the Z axis is perpendicular to the XY plane.
The computational results for the 2D RTT, 2D TRT, and 2D TTR signals with the Huiszoon permanent polarizability and the CCSD permanent polarizability are displayed in modes, which are activated by the permanent and induced optical properties, respectively. We found that the anisotropy of the permanent polarizability in the Huiszoon case is larger than that in the CCSD case. For this reason, in the Huiszoon case, the signal intensity from the permanent optical properties was stronger, whereas in the CCSD case, the signal intensity from the induced optical properties was stronger.

B. Choosing the polarizability functions
Next, we elucidate the role of the polarizability functions in 2D THz-Raman signals. While the inter-molecular interaction potentials were modeled by the TIP4P/2005 potential, 61 we employed the full-order DID, atomic site dipole-induced-dipole (ASDID), and charge-flow dipole-induced-dipole (CFDID) polarizability function models to calculate the total dipole moment and polarizability with the Ewald summation. The definitions of these polarizability functions are presented in Appendix A. We display the computational results for 2D THz- Raman signals in Fig. 5. In the 2D RTT results, the first negative peak, close to t 1 ¼ t 2 ¼ 0, which arises from the librational motion, is prominent in the DID case, while it cannot be observed in the ASDID and CFDID cases. The 2D profiles of the RTT signals are similar, whereas those of the 2D TRT and TTR signals for these three cases differ significantly. This is due to the fact that the 2D TRT and TTR signals are more sensitive to the polarizability function, because the main contribution to the signal comes from the nonlinear element of the polarizability function. 45,47 More specifically, the sensitivities of the 2D signals are due to the induced polarizabilities, as illustrated in Appendix D. conclude that the difference between the 2D THz-Raman signals in the two cases arises from the dynamical aspects of the potential models. Because the largest contribution to the signal is from the nonlinear polarizability rather than the anharmonic motion of the molecules in the 2D TRT and TTR cases, 45,47 these 2D profiles are similar. By contrast, because the main contribution to the 2D RTT signals is from the anharmonicity motion of the molecules, these signals differ significantly.

V. CONCLUSION
In this paper, we elucidated the important points involved in full MD simulations of 2D Raman and THz-Raman spectroscopic approaches. We found that the 2D signals obtained in such approaches are sensitive to the nonlinearity of the polarizability and the anharmonicity of the force fields. For the purpose of obtaining accurate 2D profiles of the signals, the proper choice of the polarizability functions is more important than that of the force field, particularly in the 2D Raman, 2D TRT, and 2D TTR cases. Although we used a small system size, we found that even in this case, in order to obtain accurate signals, we had to employ the Ewald summation for both the force field and the induced polarizability function. Although a rigorous force field model is necessary for obtaining quantitatively accurate results, the use of such a model does not guarantee that we will obtain the correct 2D profile unless we also use the right polarizability function. 69 It should be noted that conventional MD methods employ both a force field and polarizability functions in a rather empirical manner. For this reason, the verification of the calculated results is not easy. Although computationally intensive, the ab initio MD approach may be more useful for calculating 2D spectra, because it utilizes not only a force field but also a dipole moment and a polarizability based on the electronic state of the molecules. 69,70 Even in this approach, however, many assumptions and approximations are involved. Thus, verification must be carried out by comparing the 2D profiles obtained with simulations and experiments.
In this paper, we restricted our investigation to 2D Raman and 2D THz-Raman spectroscopy, but most of the points discussed in this paper also apply to the full MD simulation of 2D IR spectroscopy. 71,72 For the purpose of studying high-frequency intra-molecular modes, in addition to inter-molecular modes, which is a primary application of 2D IR spectroscopy, a quantum treatment of the atomic motion within molecules is also important. 73,74 In any case, simulating multidimensional vibrational spectroscopy of molecular liquids is a stringent test to verify the accuracy of MD simulations.

ACKNOWLEDGMENTS
The authors would like to thank Dr. Taisuke Hasegawa for providing the source code for the POLI2VS potential model. Financial support provided by a Grant-in-Aid for Scientific Research (No. A26248005) from the Japan Society for the Promotion of Science is acknowledged.

APPENDIX A: POLARIZABILTY FUNCTIONS
The optical response of molecular liquid systems is calculated in terms of the total dipole moment and total polarizability defined by lðtÞ ¼ l P ðtÞ þ l I ðtÞ (A1) and PðtÞ ¼ P P ðtÞ þ P I ðtÞ; (A2) where we separate the dipole and polarizability into permanent and induced parts as l P ðtÞ P i l P i ðtÞ and l I ðtÞ P i l I i ðtÞ, and P P ðtÞ P i P P i ðtÞ and P I ðtÞ P i P I i ðtÞ. Here, l P i ðtÞ; l I i ðtÞ; P P i ðtÞ, and P I i ðtÞ are the permanent and induced dipole and polarizability of the ith molecule, respectively.
In principle, the molecular polarizability function should be part of a force field that contributes to the dynamics of the molecular motion. In practice, however, the polarizability is usually set independently of the force field, because a nonpolarizable potential model, for example, TIP4P/2005, is optimized with respect to thermodynamic properties without a polarizable force field, while the optical response of Raman excitation is characterized by the polarizability. In the present work, we set the force field and the polarizability function independently to demonstrate the importance of the polarizability function in 2D spectroscopy. As the polarizability function, we employed the DID, ASDID, and CFDID polarizability function models. These are described below.

Dipole-induced-dipole model
The full-order DID polarizability model includes the contributions to the polarizability from other molecules. In this model, interactions between molecules are defined at the centers of individual molecules. 59 The polarizability of the ith molecule in the laboratory frame is expressed in terms of the permanent and induced polarizability as where a i is identical to the permanent molecular polarizability of the ith molecule P P i and T ij is the dipole-dipole interaction tensor Here, r ij is the vector from the center of mass of molecule i to the center of mass of molecule j, with r ij jr ij j, and 1 and are the unit matrix and the tensor product, respectively. In this DID model, the permanent molecular polarizability is obtained from P P ¼ P i a i , where a i ¼ R i P P mol R T i and R i is the rotation matrix of molecule i from the molecular frame to the laboratory frame. We also estimate the dipole moment in terms of the permanent moment, l P ¼ P i l P i , and the induced dipole moment, Here, E P i ¼ P j6 ¼i P l q P j l r ij l =r 3 ij l is the electrostatic field at molecule i created by all of the other molecules, and r ij l is the vector between the center of mass of molecule i and that of atom l in molecule j.

Atomic site dipole-induced-dipole model
In the ASDID model, the polarizabilities are defined at the centers of atoms. The molecular polarizability, P i ðtÞ ¼ P k P i k , is then obtained from the atomic polarizabilities of the individual atoms k in molecule i, expressed as where a i k is the permanent atomic polarizability of the kth atom in the ith water molecule in isolation in the laboratory frame. The ASDID model includes a screening function for the dipoledipole tensor, T i k j l , between the kth atom of the ith water molecule and the lth atom of the jth water molecule. The dipole-dipole tensor is defined as where k intra n ðr i k i l Þ and k inter n ðr i k j l Þ are the damping functions for intramolecular interactions 75 and intermolecular interactions, 76 respectively. These are expressed as k intra 3 ðr i k i l Þ ¼ 1 À expðÀs intra kl r 4 i k i l Þ (A8) and k inter 3 ðr i k j l Þ ¼ 1 À expðÀs inter kl r 3 i k j l Þ: In both cases, the higher-order damping functions are evaluated from the relation k 5 ðuÞ ¼ k 3 ðuÞ À ðu=nÞðdk 3 ðuÞ=duÞ. Here, s intra kl ¼ a intra =ða k a l Þ 4=6 and s inter kl ¼ a inter =ða k a l Þ 3=6 are the screening parameters defined in terms of the damping parameters, a intra and a inter , and the isotropic atomic polarizabilities of atoms k and l, a k and a l . 75,76 It should be noted that the electrostatic interactions attenuate at short distances due to the overlap between electron densities in realistic situations. This behavior is accounted for by the damping functions.
In the ASDID model, the permanent and induce dipole moments in the MD simulation box are obtained as l P ¼ P i l P i and l I ¼ P i P k l I i k , where l I i k is the induced atomic dipole moment of atom k, expressed as Here, E P i k ¼ P j6 ¼i P l q P j l r i k j l k inter 3 ðr i k j l Þ=r 3 i k j l is the electrostatic field at atom k of molecule i generated by all of the other molecules.

Charge-flow dipole-induced-dipole model
The CFDID model was developed as a part of the polarizable water model for inter-and intra-molecular vibrational spectroscopy (POLI2VS). 67 In the flexible-POLI2VS model, the permanent molecular polarizability, a mol k , and the permanent charge-flow polarizability, a CFP , depend on the intra-molecular and inter-molecular structure of water molecules. In this paper, we employed the CFDID model with the TIP4P/2005 and rigid-POLI2VS potentials using a fixed intra-molecular geometry. We evaluated the induced charge, q I i k , at atom k due to the influence of the surrounding water molecules. In this model, the ith molecular polarizability, , is obtained from the kth atomic dipole polarizability and the charge-flow polarizability defined by and respectively, where T i k j l is defined by Eq. (A7) and Dr i H k ¼ r i H k À r i O is the 3 Â 1 column vector connecting the kth hydrogen atom and the oxygen atom in the ith water molecule. The charge- flow dipole polarizability (a 1 Â 3 row vector) of the kth hydrogen atom is expressed as In the CFDID model, the charge-dipole interaction, T DC i k j l , and the charge-charge interaction, T CC i k j l , between atom k of water molecule i and atom l of water molecule j are utilized to calculate the total polarizability, expressed as From the MD simulations, we also evaluated the permanent and induce dipole moments, is the induced atomic dipole moment and is the induced charge on the kth atom in the i th water molecule. Here, E P i k is the electrostatic field at the position of the kth atom, and the variable DV i H k ¼ V i O À V i H k is obtained from the local electric potentials at the positions of the oxygen atom, V i O , and the kth hydrogen atom, V i H k . The potential and electrostatic field at the position of atom k of water molecule i are expressed as V i k ¼ P j6 ¼i P l q P j l =r i k j l and E P i k ¼ P j6 ¼i P l q P j l r i k j l k inter 3 ðr i k j l Þ=r 3 i k j l for the TIP4P/2005 potential. However, these variables are evaluated using the permanent multipoles expressed as in Ref. 67 for the rigid-POLI2VS potential. The induced charge, q I i O , of the oxygen atom in the ith water molecule is obtained as q

Parameterizations
To evaluate the permanent dipole moment and the polarizability of a water molecule in the 2D spectroscopy studies discussed in IV A and IV B, we set the parameters of the DID, ASDID, and CFDID models on the basis of ab initio calculations using CCSD method with the aug-cc-pVQZ basis set for a monomer water in the gas phase.
The permanent molecular dipole moment of the ith water molecule, l P i ¼ P k q P i k r i k =c, was calculated from the permanent charges q P i k of the kth atom in the ith water molecule with the empirical parameter c ¼ 1.236. The value of the permanent dipole moment was then fitted to 1.865 Debye for the TIP4P/2005 potential model. The permanent molecular polarizability in the molecular frame, P P mol , was also set at the same level of the ab initio calculation (in units of Å 3 ) as where the X axis is defined as that connecting the hydrogen atoms, the Y axis lies along the bisector of the H-O-H angle, and the Z axis is perpendicular to the XY plane. Here, the positive X direction is from the H 2 atom to the H 1 atom, and the positive Y direction is opposite to the direction of the dipole moment of the water molecule. The ab initio calculations were performed using the Gaussian09 program package 77 employing the same geometry as in the case of the rigid TIP4P/2005 potential model. In the ASDID polarizability function model, the permanent molecular polarizability, P P i ¼ P k P P i k , is calculated from the intra-molecular atomic polarizability of the kth atom in the ith water molecule, where T i k i l is defined in Eq. (A7). We set the permanent atomic polarizability of the kth atom, Then, the permanent atomic polarizability, a i k , in the laboratory frame is obtained as a i k ¼ R i a mol k R T i with the rotation matrix. In the CFDID model, the permanent molecular polarizability of water molecule i is expressed as and P IMCFP i are the intra-molecular atomic dipole polarizability of atom k and the intra-molecular charge-flow polarizability of molecule i, defined by and respectively. These polarizabilities were obtained from the GDMA program package 78 We also set the permanent charge-flow polarizability, a CFP , as a CFP H 1 H 1 ¼ 0:530Å and a CFP H 1 H 2 ¼ À0:055Å. The dipole moments and the polarizabilities in the ASDID and CFDID cases were evaluated using a O ¼ 0.862 Å 3 and a H ¼ 0.514 Å 3 with the damping parameters a intra ¼ 0.410 and a inter ¼ 0.430. Using these parameters, we calculated the permanent molecular polarizabilities of the ASDID and CFDID models as in Eq. (A18).

Numerical calculations
When calculating either the dipole moments or the polarizabilites, we have to solve the selfconsistent equations of the DID, ASDID, or CFDID model. We can express these equations in matrix form as,  ; (A27) respectively. To solve these linear equations, we used the conjugate gradient (CG) method.

APPENDIX B: EWALD SUMMATION OF INDUCED POLARIZABILITY
To elucidate the importance of the long-range electrostatic interactions, we separated the polarizability function into the permanent and induced parts as in Eq. (A2): PðtÞ ¼ P P ðtÞ þ P I ðtÞ. Then, the 1D and 2D Raman signals were analyzed in terms of permanent, induced, and cross contributions. In the 1D Raman case defined by Eq. (8), for example, the signal was decomposed into h _ P P ð0ÞP P ðtÞi; h _ P I ð0ÞP I ðtÞi, and h _ P P ð0ÞP I ðtÞi þ h _ P I ð0ÞP P ðtÞi, respectively. Accordingly, the 2D Raman signal defined by Eq. (7) can be decomposed into 8 terms. The permanent contribution to 2D Raman signal is given by hffP P ðt 2 Þ; P P ð0Þg PB ; P P ðÀt 1 Þg PB i, whereas the induced contribution is given by hffP I ðt 2 Þ; P I ð0Þg PB ; P I ðÀt 1 Þg PB i. The cross contribution is estimated from the sum of the other 6 terms. Figure 7 displays 1D Raman spectra and 2D Raman signals in terms of the permanent, induced, and cross polarizability contributions in cases with and without the Ewald summation for formamide. Because we used identical trajectories, the computational results for the permanent contribution are identical in the cases with and without the Ewald summation. While the total contributions for the 1D Raman spectra are similar in the two cases, the cross term contributions, which have opposite signs as shown in Figs. 7(a-i) and (a-ii), are different. In the 2D Raman case, due to the long-range effects that are taken into account by the Ewald summation, the induced and cross contribution are dominant in case (b), whereas the permanent contribution is dominant in case (c). Figure 8 displays the 2D Raman signals in terms of the permanent, induced, and cross polarizability contributions in the cases with and without the Ewald summation for water. Here, the induced and cross contributions are dominant and are non-negligible, even in the case without the Ewald summation. Thus, a slightly negative signal near t 1 ¼ t 2 ¼ 0 fs, which reflects the permanent contribution, appears in the total 2D Raman signal in case (b). By contrast, the intensities of the first positive peaks of the induced and cross contributions are stronger in the case with the Ewald summation. For this reason, in the total 2D Raman signal, the negative signal fades in case (a).
The differences between the case with and without the Ewald summation indicate that the long-range effects of polarizability must be accounted for when calculating 1D and 2D Raman signals.

APPENDIX C: SENSITIVITY TO THE PERMANENT POLARIZABILITY
In addition to 1D Raman, we analyzed 1D THz spectra and 2D THz-Raman signals in terms of the permanent, induced, and cross contributions of the dipole and polarizability, respectively, as illustrated for the Raman case in Appendix B.
In Fig. 9, 1D THz and 1D Raman spectra are presented. Normal mode analysis of spectra obtained from molecular dynamics simulations indicates that the peaks close to 60,200 librational motion, respectively. [79][80][81] The broadened spectral peak near 600 cm À1 in the 1D THz spectra corresponds to the librational mode, and the dominant contribution of this peak arises from the permanent dipole contribution. In the 1D THz spectra calculated with the DID model, the peak near 200 cm À1 in the THz spectra represents the translational vibration modes of the hydrogen bond that arise from the induced dipole contribution. In the 1D Raman case, two distinct hydrogen-bond vibrational peaks corresponding to the hindered translational motion of O-O-O bending and O-H Á Á Á O stretching motion are observed near 60 and 200 cm À1 , respectively. Librational motion is only detected as the peak near 500 cm À1 , while its peak position is lower than in 1D THz case, due to the difference between the distributions of the IR active modes and Raman active modes. Because the Raman process is sensitive to the choice of the polarizability model, we observe a clear difference between the spectra in the two 1D Raman cases. Contrastingly, the 1D THz spectra are insensitive to the choice of the polarizability model. Figure 10 displays the computational results for the 2D RTT, 2D TRT, and 2D TTR signals obtained using Huiszoon polarizability and CCSD polarizability for water. Each signal is decomposed into the permanent, induced, and cross-term contributions of the dipole and polarizability. The induced and permanent contributions to the 2D RTT signal are, for example, expressed as hffl I ðt 2 Þ; l I ð0Þg PB ; P I ðÀt 1 Þg PB i and hffl P ðt 2 Þ; l P ð0Þg PB ; P P ðÀt 1 Þg PB i, respectively. The permanent and induced contributions arise from the libration and translation motion, while the cross contribution arises from the libration-translation coupling. 50,82 The permanent contribution to the 2D THz-Raman signal exhibits a characteristic profile that can be analyzed using a nonlinear Brownian model discussed in Refs. 45 and 47. Then, the induced part, for example, the 2D RTT signal, can be attributed to the contribution from the nonlinearity of the induced dipole moment and the anharmonic motion of the molecular vibration. It is important to note that, while the profile of each contribution presented in Fig. 10 is similar in the Huiszoon and CCDS polarizability cases, the total RTT, TRT, and TTR profiles presented in Fig. 4 differ significantly in these two cases. This is because the relative intensities of the permanent, induced, and cross contributions have different dependences on the polarizability function in the two cases. To obtain the correct 2D profile, we must describe both the dynamics and the observables of the molecules very accurately.

APPENDIX D: SENSITIVITY TO THE POLARIZABILITY FUNCTIONS
To elucidate the cause of the differences among the 2D profiles in the DID, ASDID, and CFDID cases depicted in Fig. 5, we evaluated the 2D RTT, TRT and TTR signals by separating the two induced dipole moments and the permanent polarizability (ID) part and by separating the two permanent dipole moments and the induced polarizability (IP) part. In the RTT case, for example, the ID and IP parts are defined by hffl I ðt 2 Þ; l I ð0Þg PB ; P P ðÀt 1 Þg PB i and hffl P ðt 2 Þ; l P ð0Þg PB ; P I ðÀt 1 Þg PB i. The sensitivities of the signals to the induced dipole moment and the induced polarizability were then analyzed by considering the ID and IP parts of the signal. It should be noted that because we fixed the permanent polarizabilities to have the values presented in Eq. (A18) for all three of the polarizability models, the differences among the 2D profiles arise only from differences in the induced dipole moment and the induced polarizability. The similarity of the 2D THz-Raman signals obtained with the three models presented in Figs. 11(a-i)-11(a-iii) indicates that these signals are not sensitive to differences in the induced dipole moment. In the 2D RTT case depicted in Fig. 11(b), because the signal profiles are determined by the linear elements of the dipole moments and the polarizabilities, while the nonlinear elements of the polarizability play only a minor role, the difference in the induced polarizability does not strongly influence the RTT signals. Contrastingly, because the signal profiles of the 2D TRT and TTR signals are determined by the nonlinear polarizabilitiy, 45,47 there are clear differences among the 2D TRT and TTR signals.
In the case of the DID model, the librational photon echo peak that appears in the 2D TRT signal is slightly blurred due to the anharmonicity of the molecular motion, 42,43,45,47 while in the CFDID case, there is a clearer negative photon echo peak along the t 1 ¼ t 2 direction arising from the IP contribution.