Coupling between intra- and intermolecular motions in liquid water revealed by two-dimensional terahertz-infrared-visible spectroscopy

The interaction between intramolecular and intermolecular degrees of freedom in liquid water underlies fundamental chemical and physical phenomena such as energy dissipation and proton transfer. Yet, it has been challenging to elucidate the coupling between these different types of modes. Here, we report on the direct observation and quantification of the coupling between intermolecular and intramolecular coordinates using two-dimensional, ultra-broadband, terahertz-infrared-visible (2D TIRV) spectroscopy and molecular dynamics calculations. Our study reveals strong coupling of the O-H stretch vibration, independent of the degree of delocalization of this high-frequency mode, to low-frequency intermolecular motions over a wide frequency range from 50 to 250 cm−1, corresponding to both the intermolecular hydrogen bond bending (≈ 60 cm−1) and stretching (≈ 180 cm−1) modes. Our results provide mechanistic insights into the coupling of the O-H stretch vibration to collective, delocalized intermolecular modes.


Supplementary Note 1: Intensity spectrum of the THz laser pulse at the position of the sample
To characterize intensity spectrum of the THz laser pulse at the spatial overlap with the IR and VIS beams we employ terahertz-field-induced second harmonic generation (TFISH) in nitrogen gas. To this end, we block the IR beam and measure the spectrum of the UV light at the frequencies of 2 VIS ± THz (≈400 nm) produced by four wave mixing of the VIS and THz pulses in nitrogen. Intensity of such signal is THz , where (3) is the third-order nonlinear optical susceptibility of nitrogen, VIS and THz are intensities of the VIS and IR beams, respectively. Indeed, Supplementary Fig. 6a shows linear dependence of the signal intensity on the intensity of the THz pulse which confirms the nature of the signal. Magenta line in Supplementary Fig. 6b shows the intensity spectrum of the TFISH signal. Comparison of the TFISH spectrum with the second harmonic of the VIS pulse (black line in Supplementary Fig. 6b), which is generated at the mirrors (M1 and M2 in Supplementary Fig.4), demonstrates that the former is produced by 2 VIS − THz wave mixing. Thus, because we use narrowband VIS pulse the magenta trace in Supplementary Fig. 6b amounts to the THz spectrum shifted by 2 VIS . The THz spectrum derived from the TFISH measurement ( Fig. 2b) is in a good agreement with the spectrum measured previously by the air based coherent (ABC) detection ( Supplementary Fig. 6c,d) 1 . We note that the absolute-value 2D TIRV spectra are not sensitive to the variation of the phase across the THz spectrum (see Supplementary Note 2). We also stress that for the ABC detection the THz beam is overlapped in space with a sampling beam in the ABC detector. The spatial overlap of the THz and sampling beams in the ABC detector can be different from the overlap of the THz beam with the IR and VIS beams in the 2D TIRV measurement. This difference can be the reason, at least partially, for the difference in the spectra of the THz pulse measured by different methods. Figure 6. Intensity spectrum of the THz pulse. a, Intensity of the TFISH signal for different intensities of the THz laser pulse. The magenta trace shows experimental data and the black line shows a linear fit. b, Spectrum of the TFISH signal (magenta) and of the second harmonic of the VIS pulse (black). c, THz pulse measured using the air based coherent detection. d, Intensity spectrum (magenta) and phase (blue) of the THz pulse obtained by Fourier transform of the signal in c.

Supplementary Note 2: Perturbation formalism for the 2D TIRV spectroscopy
In this section, we derive the relationship between the signal measured by the 2D TIRV spectroscopy and the third-order nonlinear response function (3) of a sample. Electric field (3) of the 2D TIRV signal emitted by a sample after interactions with the THz, IR and VIS pulses is proportional to the third-order nonlinear polarization generated by these interactions. Thus, by employing the frequency domain formalism (Eq. (5.32) in Ref. 2 ) we obtain for time delay of the THz pulse: where frequency s = 1 + 2 + 3 ; THz ( 1 ), IR ( 2 ) and VIS ( 3 ) are electric fields of the THz, IR and VIS pulses at the frequencies 1 , 2 and 3 , respectively. Dispersion of the signal by the grating of the spectrometer is given by Fourier transform of the signal field over time : where is the Dirac delta function. Because after the grating positive and negative frequency components of light travel in the same direction the signal measured by the square-law detector (CCD camera) for frequency Ω 3 ≥ 0 in the presence of the local oscillator (LO) and integrated over time period is given by: where LO (Ω 3 ) is electric field of the LO at frequency Ω 3 . By squaring the term under integral and neglecting highly-oscillating terms we obtain for the heterodyne-detected signal (interference between the signal field and the LO): Heterodyne-detected signal depends parametrically on the time delay of the THz pulse, which is reflected by the time-domain 2D TIRV data in the experiment. In order to obtain 2D spectra we perform Fourier transform of this signal over : Using Supplementary Eq. (2) we obtain for (3) (Ω 3 , Ω 1 ): In the experiment, we use narrowband VIS pulse for which we approximate the spectrum by delta function at frequency Ω VIS > 0: By substituting Supplementary Eq. (7) into Supplementary Eq. (6) and integrating we obtain: Because in our experiment Ω 1 ≪ Ω 3 + Ω VIS the frequency Ω 3 + Ω VIS − Ω 1 is beyond the bandwidth of the IR pulse. Thus, IR (Ω 3 + Ω VIS − Ω 1 ) = 0 and the second term in Supplementary Eq. (8) vanishes. Similar, for (3) (−Ω 3 , Ω 1 ) we obtain: We introduce frequency Ω 2 = Ω 3 − Ω VIS to write Supplementary Eq. (10) in the alternative form: (Ω 2 + Ω VIS , Ω 1 ) ∝ (3) (Ω 2 + Ω VIS , Ω 2 , Ω 1 ) THz (Ω 1 ) IR (Ω 2 − Ω 1 ) LO (−Ω 2 − Ω VIS ) + (3) (−Ω 2 − Ω VIS , −Ω 2 , Ω 1 ) THz (Ω 1 ) IR (−Ω 2 − Ω 1 ) LO (Ω 2 + Ω VIS ). (11) Thus, in the experiment, the 2D TIRV spectrum is given by the product of the spectrum of the third-order nonlinear response function with the spectra of the laser pulses. The presence of coupled vibrational resonances at frequencies Ω 1 (THz frequency range) and Ω 2 (IR frequency range) results in the peaks in the spectrum of the (3) response function and, therefore, in the 2D TIRV spectrum.
Supplementary Equation (11) shows that 2D TIRV spectra in our experiment are produced by sum of the (3) spectra in two quadrants, (Ω 2 , Ω 1 ) and (−Ω 2 , Ω 1 ). This is due to the utilization of the spectrometer to measure the second frequency axis and is analogous to the summation of the rephasing and non-rephasing spectra in the 2D IR spectroscopy when a spectrometer and the pump-probe geometry are used.

Supplementary Note 3: Change of the 2D TIRV signal intensity with sample composition
The gain of the signal field sig produced by FWM of the THz, IR and VIS beams over the distance in the sample is given by 3 : where (3) is the third-order nonlinear optical susceptibility of the medium; THz ( ), IR ( ) and VIS ( ) are electric fields of the THz, IR and VIS pulses at the position inside the sample; ∆ = THz + IR + VIS − sig is the mismatch between the wavevector of the signal and those of the incident waves. Upon propagation into the sample the laser beams are depleted by absorption as well as conversion into the signal field. We neglect the latter because FWM in our experiment is much weaker than the linear absorption.
where , OH and OD and are the concentrations of all molecules, O-H and O-D oscillators, respectively. By integrating the differential Supplementary Eq. (13) we obtain signal emitted by the sample of thickness : Where we have introduced the total extinction coefficient = IR OH OH + IR OD OD + THz . For 1 mm thick samples that we use ≫ 1 for all isotope dilutions and thus the amplitude of the signal:  Fig. 3g is 0.95/1.07=0.89 of that in Fig. 3f. We use this coefficient to subtract the two spectra in order to eliminate the D 2 O signal in the 5% H/D sample and obtain 2D TIRV spectrum for the isolated O-H stretch oscillator ( Fig. 4a and Supplementary Fig. 3a).

Molecular dynamics (MD) simulation Protocols.
MD simulations were performed to calculate the 2D TIRV signals of water. We employed the POLI2VS water potential 4 . We used the cubic unit cells with the length of 12.426 Å and periodic boundary conditions in all directions. The non-electrostatic intermolecular interaction and the electrostatic interactions including quadrupole terms were smoothly cut off into zero from 6.0 Å to 6.1 Å. The other electrostatic interactions were calculated by the standard Ewald sum. The equations of motion were integrated by the r-RESPA algorithm 5 . A 0.5 fs time step was used for integrating the equations of motion for the intermolecular interactions with the 2nd order symplectic integrator, while a 0.25 fs time step was used for the intramolecular bonding interaction with the 6th order symplectic integrator 6 . In the simulation of the neat H 2 O, the system consisted of 64 water molecules. In the simulation for the mixture of 5% H 2 O and 95 % D 2 O, the system consisted of 58 D 2 O and 6 HOD molecules. We randomly generated 20 initial structures for both systems. The constant temperature (NVT) MD simulations were performed for 100 ps to equilibrate each system at 300 K, followed by the production runs. The Nose-Hoover chain algorithm was employed in the constant temperature (NVT) MD runs 7 . After the equilibration, the system temperatures were around 300 K in the constant energy (NVE) MD simulations. We employed the equilibrium and non-equilibrium hybrid algorithm to calculate 2D TIRV signals 8,9 . The IR laser pulses with the strengths of 0.1 V/Å were applied to the system for 0.5 fs. Note that the increase in temperature observed after the laser interactions was ~ 0.5 K. We collected total 10 6 nonequilibrium trajectories from the independent 20 equilibrium trajectories for each system. We performed the NEMD and EMD simulations for 250 fs and collected the polarizabilities and dipole moments every 1 fs. The response functions were calculated in the time periods of To take into account the limited bandwidth of the laser pulses we convolute the calculated response function with the THz and IR laser pulses, which is given by: