Two-dimensional infrared-Raman spectroscopy as a probe of water's tetrahedrality

Two-dimensional spectroscopic techniques combining terahertz (THz), infrared (IR), and visible pulses offer a wealth of information about coupling among vibrational modes in molecular liquids, thus providing a promising probe of their local structure. However, the capabilities of these spectroscopies are still largely unexplored due to experimental limitations and inherently weak nonlinear signals. Here, through a combination of equilibrium-nonequilibrium molecular dynamics (MD) and a tailored spectrum decomposition scheme, we identify a relationship between the tetrahedral order of liquid water and its two-dimensional IR-IR-Raman (IIR) spectrum. The structure-spectrum relationship can explain the temperature dependence of the spectral features corresponding to the anharmonic coupling between low-frequency intermolecular and high-frequency intramolecular vibrational modes of water. In light of these results, we propose new experiments and discuss the implications for the study of tetrahedrality of liquid water.

visible pulses, such as the THz-THz-Raman, [31][32][33][34] THz-Raman-THz, Raman-THz-THz, [35][36][37][38] and THz-IR-Raman (also called THz-IR-visible 39,40 or TIRV). Their development was enabled by the recent advances in the generation of strong THz pulses that are needed to induce a nonlinear light-matter interaction. 41 The THz-THz-Raman and TIRV methods are related to other two-dimensional IR-Raman techniques, namely the two-dimensional IR doubly vibrationally enhanced 42 and IR-IR-visible sum-frequency generation spectroscopies. 43 For example, TIRV experiments have revealed unambiguous spectral signatures of coupling between the intramolecular O-H stretch and intermolecular hydrogen bond bending and stretching modes. 39 Theoretical simulations by Ito and Tanimura 44 predicted such spectral features and assigned them to both mechanical and electrical anharmonic coupling between the said vibrational modes. Similarly, THz-THz-Raman spectroscopy recently revealed signatures of anharmonic coupling between phonons of ionic solid LiNbO 3 . 45 Finally, Raman-THz-THz and THz-Raman-THz spectroscopies have been used to study the inhomogeneity of liquid water and aqueous solutions, 46,47 as well as the coupling among intermolecular and intramolecular modes in liquid and solid bromoform. 48,49 Even so, due to a limited availability of efficient THz emitter materials, not all frequencies have been covered by the reported techniques. Specifically, most two-dimensional hybrid THz-Raman spectroscopies of liquid water targeted hydrogen-bond bending and stretching modes, i.e., frequencies up to 400 cm −1 , leaving the water librational dynamics largely unexplored.
Here, we aim to provide new insights into the capabilities of two-dimensional hybrid IR-Raman vibrational spectroscopies. To this end, we study the temperature dependence of the two-dimensional IR-IR-Raman (IIR) spectrum, which is given by the double Fourier (or sine) transform of an appropriate third-order, two-time response function. Since the computational model involves all vibrational modes of the system, the response function covers a broad range of frequencies, which, in practice, can be mapped out only through separate experiments. To date, only the TIRV frequency region has been experimentally measured, although its temperature dependence was not studied. Second, following Ref. 26, we then establish a structure-spectrum relationship by separating the spectral contributions from molecules exhibiting low or high tetrahedral coordination. Third, to justify the molecular dynamics (MD) results at low temperature, we analyze whether nuclear quantum effects are discernible in the two-dimensional IIR spectrum.

Theoretical model
We simulated the IIR response function 33,44 using the equilibrium-nonequilibrium MD approach, in which the quantum-mechanical trace is replaced by a classical average 50-53 whereρ is the thermal density operator,μ = µ(q) is the dipole moment operator, and Π = Π(q) is the polarizability. q t denotes the position of a classical trajectory initiated at (q 0 , p 0 ), whereas q ±,t corresponds to the initial conditions (q ±,0 = q 0 , p ±,0 = p 0 ± εµ (q 0 )/2) after an instantaneous interaction with the electric field. ε is a free parameter in the calculations and corresponds to the magnitude of the external electric field integrated over the short interaction time. The two-dimensional spectra are computed through a double sine To model water, we used a fast and accurate qTIP4P/F force field, 55 which has been well studied for spectroscopic simulations. 25 To allow for nonlinear dependence of the dipole moments and polarizabilities on nuclear coordinates, we employed the truncated dipoleinduced-dipole (DID) model. Following Hamm, 37 each water molecule was amended with a permanent anisotropic polarizability, which was used for the evaluation of the induced contributions to the dipoles and polarizabilities. In contrast to the rigid water simulations of Ref. 37, here the permanent polarizability was an explicit function of intramolecular degrees of freedom, according to Ref. 56. We note that more advanced models for the potential energy, dipoles, and polarizabilities of liquid water exist and have been used in the simulation of one-and two-dimensional spectroscopies. 19,22,24,26,37,57 Nevertheless, our approach is sufficiently accurate for the analysis that is presented in this work. In Fig. 1 Most of the results presented below rely on the validity of the classical MD approach, which neglects quantum-mechanical properties of atomic nuclei. However, due to the presence of light hydrogen atoms, their effect might not be negligible. 55 Although such nuclear quantum effects on one-dimensional IR and Raman spectra have been well studied using approximate but reliable classical-like methods, no tools similar to the equilibriumnonequilibrium MD have been available to study nuclear quantum effects on the twodimensional IR-Raman spectra. [58][59][60] Recently, we have developed a new RPMD approach, 61 which reproduces correctly the quantum-mechanical thermal distribution and preserves it during approximate classical dynamics. In this way, the equilibrium-nonequilibrium RPMD method can simulate, at least approximately, the nuclear quantum effects on the twodimensional IIR spectra. Briefly, the RPMD method 62 replaces the original quantummechanical problem by an extended classical system consisting of N replicas (beads) of the original system connected by harmonic springs. In the limit of N → ∞, the extended classical system reproduces exact quantum-mechanical thermal distribution, while if N = 1, RPMD reduces to classical MD. In our simulations, we used N = 32, which is sufficiently large for liquid water in the studied temperature range. 25 Since RPMD is known to suffer from the spurious resonance issue, where unphysical peaks due to artificial harmonic springs appear in the spectra, we employed its thermostatted version (TRPMD). 63,64 IR absorption and anisotropic Raman spectra simulated with TRPMD are presented in Fig. 1.
Two-dimensional IIR spectrum of liquid water  which we will refer to as the TIRV region, are a product of interplay between mechanical and electrical anharmonicity. Namely, the mechanical anharmonicity leads to a shape comprising a positive and a negative lobe above and below the central ω 2 frequency, whereas the electrical anharmonicity produces an approximately symmetric feature. In the lowfrequency region of interest, we observe a strong peak at about (600 cm −1 , 1650 cm −1 ) due to anharmonic coupling between librations (400 − 1000 cm −1 ) and intramolecular bending mode. 65 The same coupling mechanism is responsible for the combination transition at around 2150 cm −1 , also known as the "association" band, 9,20 in the one-dimensional spectra.
Interestingly, this combination band is not captured in our model (see Fig. 1), even though the corresponding two-dimensional spectral feature clearly appears in the IIR spectrum.
Temperature dependence of the IIR spectrum We now turn to the temperature dependence of the two IIR spectral regions (Fig. 3). The TIRV part of the spectrum experiences several changes as the temperature is increased from 280 K to 360 K. At higher temperatures the peaks are blue-shifted along ω 2 , which agrees with the observed trends in the experimental Raman spectra and simulated vibrational density of states. 26 In addition, the shape of the spectral peaks changes drastically. The most prominent change in the spectrum is the disappearance of the negative peak around (250 cm −1 , 3600 cm −1 ). According to the interpretation based on model simulations, this indicates that the electrical anharmonicity contribution to the spectrum, marked by strong symmetric lineshape, increases compared to that of the mechanical anharmonicity. Unlike the frequency shift, this feature cannot be accessed through steady-state spectroscopy. Finally, we note that a peak around (400 cm −1 , 3700 cm −1 ) appears at elevated temperature, while the other peak at about (700 cm −1 , 3400 cm −1 ) almost disappears. Interestingly, both of these coupling terms have been discussed as possible origins of the libration-stretch combination band appearing in the anisotropic Raman spectrum of liquid water at 4100 cm −1 . 26,66 The same combination band has also been shown to play an important role in second-order vibrational sum-frequency spectroscopy of interfacial water. 15 The low-frequency region also exhibits strong temperature dependence. For example, the results (see also difference spectra in the Supplementary Fig. 5) clearly indicate a red shift of the libration-bend peak along ω 1 with increasing temperature. This is a straightforward consequence of the equivalent frequency shift found in the linear IR absorption spectrum (see Supplementary Fig. 6), which also agrees with the experimentally observed trend. 17 Furthermore, the two-dimensional IIR spectrum at 280 K contains peaks close to diagonal, around ω 1 = ω 2 = 250 cm −1 , due to mechanical anharmonicity of the hydrogen-bond stretching modes. These seem to progressively disappear at elevated temperatures. Importantly, this change in intensity could not be simply predicted from the one-dimensional spectra, which show little change in the intensity of the 250 cm −1 band. We note, however, that the features in this congested spectral region could be affected by the short times available from our simulations and by artificial broadening. Longer simulations are possible with rigid water models, which have been used to study explicitly the intermolecular modes and corresponding two-dimensional THz-Raman signals in the time domain. 37,51 Spectral signatures of water's tetrahedrality in the IIR spectrum To understand the temperature-dependent spectral features, we follow the work of Morawietz et al., 26 where the temperature effects were analyzed in relation to the local structuring around individual water molecules. This can be quantified by the local tetrahedral order parameter Q, 67-69 which measures how the arrangement of the four neighboring water molecules deviates from the ideal tetrahedral arrangement around the central one. By convention, Q = 1 for a perfect tetrahedral arrangement, while Q = 0 for an ideal gas. The distribution of Q at thermal equilibrium exhibits a strong temperature dependence and a clear isosbestic point at Q ≈ 0.67 (see Fig. 4a), attributed to the existence of two components whose populations change with temperature. 26 To directly correlate the local tetrahedral order parameter with the spectral features observed in the two-dimensional IIR spectra, we decomposed the spectrum into contributions of individual molecules (see Supplementary Section VII). Then, we simulated the spectra (see Fig. 4b) originating predominantly from the molecules with either high (Q > 0.72) or low (Q < 0.62) local tetrahedral order parameter. The two IIR spectra align remarkably well with the observed temperature-dependent changes. Specifically, the high-order TIRV spectrum (Fig. 4b, top right) exhibits a clear mechanical anharmonic coupling feature with a positive and a negative lobe at ω 1 ≈ 250 cm −1 , whereas the corresponding low-order spectrum contains a strong symmetric peak, indicating electrical anharmonic coupling between the hydrogen-bond modes and O-H stretch. In addition, the libration-stretch peak at (400 cm −1 , 3700 cm −1 ) appears almost exclusively in molecules with a low tetrahedral order.
Similarly, the libration-bend peak in the low-frequency part the spectrum (bottom panels of Fig. 4b) appears at lower ω 1 frequencies for low Q. Overall, the temperature dependence of the IIR spectra can be almost exclusively assigned to the changes in the distribution of the local tetrahedral order parameter and the effect of the local structure on the spectral features. Therefore, TIRV spectroscopy expanded into the water libration frequency range, i.e., with ω 1 covering up to ≈ 800 cm −1 could provide an alternative probe of the local molecular structure in liquid water and aqueous solutions. Importantly, the spectral changes in IIR are more drastic, and therefore more sensitive to local ordering, than the frequency shifts observed in conventional, one-dimensional IR and Raman spectroscopies. Nuclear quantum effects on the two-dimensional IIR spectrum Finally, we report the first TRPMD simulation of the two-dimensional IIR spectrum of liquid water and compare it with the MD simulation in Fig. 5. Some nuclear quantum effects can be easily explained in terms of the differences between one-dimensional spectra (see Supplementary Fig. 7). Namely, the red-shifted O-H stretch band in the TRPMD steady-state spectra reflects in a red shift along ω 2 of the spectral features in the TIRV spectrum (Fig. 5, top). Furthermore, the libration-stretch peak shifts from ω 1 ≈ 700 cm −1 in the MD spectrum to ω 1 ≈ 800 cm −1 in TRPMD, which aligns with the differences between MD and TRPMD IR absorption spectra. In the low-frequency part of the spectrum, the libration-bend peak is analogously shifted along ω 1 . However, some features cannot be understood from comparison with the one-dimensional spectra, such as the absence of the strong negative feature at (ω 1 ≈ 750 cm −1 , ω 2 ≈ 50cm −1 ) in TRPMD or the change in the spectral lineshapes in the TIRV region at ω 1 > 400 cm −1 . Nevertheless, these differences are not sufficiently large to alter the above conclusions based on classical MD simulations.

DISCUSSION
To conclude, we have presented MD and TRPMD simulations of the two-dimensional IIR spectra of liquid water. By analyzing the temperature dependence of the spectrum, we have identified features that report on the degree of local tetrahedral ordering in the first coordination shell. We envision that this structure-spectrum relationship could also be used to probe the structure of complex aqueous systems, such as confined water and aqueous solutions. Further computational work is needed to confirm that IIR spectra of such systems can be broadly mapped to the local tetrahedral order parameter. Our work demonstrates that the temperature dependence of such two-dimensional spectra can be used in combination with computational modeling to gain insight into the structure of complex liquids.
Furthermore, the simulations presented here imply that the contributions of mechanical and electrical anharmonic coupling change as a function of local tetrahedrality, which cannot be studied with conventional, one-dimensional spectroscopy. Specifically, the electrical anharmonicity dominates at higher temperatures, where the tetrahedral order parameter is low, whereas the mechanical coupling, due to the anharmonic terms in the potential energy surface, are pronounced at lower temperature and higher tetrahedrality. The spectral region in which we identified the main signature of tetrahedrality has been studied experimentally. 39 However, the experiments could only measure an absolute value Fourier transform spectrum (see Supplementary Section III and Supplementary Fig. 3), which is quite different from the real sine transform discussed here and in Ref. 44. Hence, further experimental work is needed to successfully measure the features that, according to our prediction, correlate with the local tetrahedral order parameter.
Finally, we observe that the 400 cm −1 < ω 1 < 1000 cm −1 region of the IIR spectrum, corresponding to the librations of water molecules, carries rich information about liquid water and should be further explored experimentally. In particular, such studies could provide additional information on the combination bands appearing in IR absorption and Raman scattering spectra that have challenged physical chemists for decades.

Dipole moment and polarizability models
The induced dipole moments and polarizabilities were modeled as 37 where E ij = a∈{O, H1, H2} Q ja r ja,iO /r 3 ja,iO is the electric field produced by molecule j on the oxygen atom of molecule i, T ij = T( r jO,iO ), T( r) = (r 2 − 3 r ⊗ r)/r 5 is the dipole-dipole interaction tensor, r ja,iO = r ja − r iO , r ja denotes the position vector of atom a ∈ {O, H1, H2} of molecule j, Q ja is its partial qTIP4P/F charge, r = r denotes the norm of vector r, and ⊗ is the outer product of two vectors. Permanent dipole moment of molecule i was µ i = a∈{O, H1, H2} Q ia r ia /c, where c = 1.3 is a scaling factor that reduces the dipole moment of water to its gas-phase value. 37,52 The permanent polarizability of each molecule i, Π i ≡ Π i ( r iO , r iH1 , r iH2 ), was computed according to Ref. 56 in the molecular reference frame and rotated into the laboratory frame. Then, the total dipole moment and polarizability are given by Finally, we note that the sums over i in all expressions above are taken over the molecules in a unit cell, while the sums over j extend to infinity beyond the central unit cell. Due to the slow convergence of Coulomb interactions in the direct space, we used Ewald summation for the electric field E ij and tensor T ij . 22,52 Computational details Classical thermal average of equation (2)  Finally, starting from positions at which the dipole derivatives were evaluated, "nonequilibrium" trajectories, i.e., with modified initial momenta p ±,0 (see main text), were propagated for 1000 steps (250 fs). This generated 2.56 × 10 5 samples for equation (2) of the main text, which is comparable to the numbers used in Refs. 39,51,52, and the response function was computed for 250 fs in both t 1 and t 2 . The statistical error was estimated using bootstrapping and is shown in Supplementary Fig. 1. We set ε = 2, which converts into the electric field magnitude of E 0 = ε/∆t ≈ 10 V/Å. Weaker electric field was also tested (ε = 0. Linear absorption and anisotropic Raman spectra were computed from the equilibrium trajectories according to: where µ is the dipole moment vector of the cell,

DATA AVAILABILITY
The data that support the findings of this study, together with the input files and scripts for plotting the data, are openly available in Zenodo repository with the digital object identifier (DOI): 10.5281/zenodo.7265859.

CODE AVAILABILITY
Modified i-PI code used to run MD and TRPMD simulations is available at https://github.com/tbegusic/i-pi.git; encorr code used for processing the outputs is available at https://github.com/tbegusic/encorr.git.

ACKNOWLEDGMENTS
The

III. COMPARISON WITH THE TIRV EXPERIMENT
To compare with the TIRV experiment and simulations of Ref. 39, we also evaluated where E THz (ω) and E IR (ω) are the electric fields obtained as square roots of the intensity spectra presented in Fig. 2b,c of Ref.

IV. TWO-DIMENSIONAL MODEL SYSTEM
Following Ito and Tanimura, 44 we simulated the spectra of two different two-dimensional models with harmonic Hamiltonian ω 1 = 0.5, ω 2 = 2, and µ 1 = q 1 and µ 2 = q 2 , where µ 1 and µ 2 are the dipole moments for the interaction with the first and second light pulses. For the model with mechanical anharmonicity, we set the Hamiltonian to where α i = 2.5 × 10 −5 ω 4 i and λ = 0.1, and we set Π = q 2 . The weak quartic terms are added to ensure that the potential is bound. For the model with electrical anharmonicity, the Hamiltonian was quadratic, H = H 0 , while Π = −0.05q 1 q 2 . The spectra were simulated with exact quantum mechanics on a grid in position representation, similar to the benchmark simulations of Ref. 61. The response function was simulated for maximum t 1 and t 2 times of 50 with a time step of 0.25. Exponential damping exp(−(t 1 + t 2 )/τ ) with τ = 1 was applied to the response function before evaluating the discrete sine transform.

VII. SPECTRUM DECOMPOSITION
Assuming sufficiently small ε, we can expand Π(q ±,t ) = Π(q t ) ± ε 2 dΠ(q t ) dp 0 · dµ(q 0 ) dq 0 (5) to first order in ε, where the sum goes over N mol molecules, and q i and p i denote positions and momenta of atoms in molecule i. Then, it can be shown that where and q (i) ±,t is the position vector (all 3N atom coordinates) of a trajectory with initial momentum p (i) i.e., a trajectory with the electric field applied only on molecule i. In our simulations, we did not decompose the spectrum into individual molecules, but into two groups of molecules exhibiting low and high tetrahedral order parameter, which was computed at time zero.
Spectra evaluated in this way were divided by the fractions of molecules of a given order in the thermal distribution, which were x(Q < 0.62) =