Anharmonic exciton dynamics and energy dissipation in liquid water from two-dimensional infrared spectroscopy

It is well understood that water’s remarkable attributes and dynamical complexity come from its rapidly fluctuating extended hydrogen-bonding network.¹ Although no current technique possesses the ability to directly measure molecular structure in liquids on femtosecond time scales, nonlinear time-resolved infrared (IR) spectroscopy has proven itself to be highly successful in the study of aqueous systems due to its sensitivity to hydrogen bonding.^2,3 The majority of these experiments, however, have been conducted on isotopically dilute water with comparatively few focusing on isotopically pure water.^4–6 In isotopically dilute systems, it is well understood that the frequency of the dilute oscillator is highly correlated with the electric field projected along the bond, which is almost entirely determined by its hydrogen-bond configuration.^7–10 Red-shifted absorption corresponds to strong hydrogen bonds, while blue-shifted absorption corresponds to weak ones. The strong correlation between frequency and structure in addition to sub-picosecond time resolution has made ultrafast IR spectroscopy immensely successful in the study of these systems.

One is naturally led to ask whether such simple correlation between frequency and structure exists for neat water. In isotopically pure water, the high density of resonant oscillators results in vibrations that are delocalized due to strong coupling between them. At the least, we may consider symmetric and anti-symmetric stretches of a single water molecule, but a more realistic picture contains excitonic vibrations delocalized over multiple molecules. Therefore, while it will certainly be true that, on average, strong hydrogen bonding leads to a red shift, the exact frequency of an excitonic vibration depends on the configuration of an entire cluster over which a vibration is delocalized.

The earliest femtosecond IR spectroscopy studies on H₂O employed transient absorption (TA) spectroscopy. By comparing polarization anisotropy measurement to a series of isotope dilutions, it was concluded that resonant energy transfer between O–H oscillators was the dominant effect in the vibrational dynamics.^11,12 The results were interpreted in terms of Förster resonance energy transfer and neglected the role of strong intermolecular coupling between oscillators. Follow-up studies probed the temperature dependence of the lifetime of the O–H stretch and found it to increase modestly with temperature from 260 fs at 298 K to 320 fs at 358 K.¹³ This was attributed to reduced overlap between the bend overtone and stretch fundamental, since it was found using IR-Raman spectroscopy that excitation of the stretch resulted in a significant response from the bending mode.¹⁴ In these early studies, spectral diffusion was found to occur on a 0.5–0.7 ps time scale.^13,15,16

With the first 2D IR measurements^17,18 and additional TA measurements,^19,20 it became increasingly clear that delocalization effects were particularly important in the vibrational spectroscopy of water and the concept of vibrational excitons in water was introduced.²¹ These studies were at odds with the previous ones employing the Förster mechanism of incoherent vibrational energy transfer as the main driver of vibrational dynamics. Furthermore, these interpretations were supported by theoretical considerations that showed that excitonic delocalization is indeed important^21–23 and that the Förster model is not valid in this limit.²⁴ 

With recent advances in IR pulse generation, 2D IR experiments were performed on H₂O using sub-50 fs pulses and a probe pulse with a 2000 cm⁻¹ bandwidth that spans the entire spectral range of the O–H stretching band in water as well as the HOH bending mode.²⁵ This study corroborated the picture of delocalized vibrational excitons in H₂O but emphasized the importance of mixing of the O–H stretch with bending and low-frequency modes (LFMs). However, the broader bandwidth illustrated two points that were at odds with previous literature. With an overall time resolution of ∼70 fs, it was shown that the time scale for spectral diffusion is 175 fs, which is much longer than the value of 50 fs obtained with a longer time resolution and narrower bandwidth.¹⁷ Additionally, these data showed direct and ultrafast relaxation to LFMs as opposed to the cascading stretch-to-bend mechanism proposed based on transient absorption data,¹⁹ which is at odds with previous TA studies, but is consistent with previous IR-Raman measurements.^14,16 Importantly, this study showed the clear influence of the temporal and spectral profiles of the pulses on the conclusions drawn by experiments.

In this work, we extend our previous results on broadband 2D IR spectroscopy of H₂O to capture a broader range of dynamic phenomena in water. We present high signal-to-noise 2D IR and TA measurements of H₂O’s vibrational dynamics, acquired as a function of polarization and temperature, and interpret these results in terms of the highly anharmonic nuclear potential energy surface (NPES) due to hydrogen bonding. We use 2D IR depolarization measurements to further strengthen the case for delocalized excitonic vibrations, discuss the effect of delocalization on the dynamic parameters that can be extracted from the 2D IR spectrum, and show that the vibrational dynamics are largely insensitive to temperature. Finally, we address the non-equilibrium nature of the so-called heating, or hot ground state (HGS), commonly observed in time-resolved IR spectroscopy of water. These results allow us to refine our picture of molecular dynamics in liquid water and provide stringent constraints on the simulation and interpretation of vibrational spectra of H₂O.

2D IR spectroscopy was used to characterize the time-dependent spectral correlations between water vibrations excited with an intense pump beam resonant with the O–H stretching mode of liquid H₂O at 3400 cm⁻¹ and a weak broadband IR probe pulse with detection frequencies between 1500 and 3700 cm⁻¹. Many of the details of the experiment have been described previously, which we will summarize briefly in addition to describing new methods.

The output of a commercial 1 kHz Ti:Sapphire regenerative amplifier (Coherent Legend USX; 5 mJ energy, λ = 800 nm, Δλ = 70 nm) is split into two beams with energies of 500 and 600 μJ and fed into independent grating compressors, where they are compressed to a final pulse length of less than 25 fs. The two outputs drive a homebuilt optical parametric amplifier (OPA) and broadband continuum source. Pump pulses centered at 3400 cm⁻¹ are generated in a two-stage BBO–KNbO₃ OPA, which has been described previously.²⁶ The pulses have an energy of 3 μJ, a center wavelength and bandwidth of 3400 cm⁻¹ and 380 cm⁻¹, respectively, and a pulse length of 45 fs. Broad bandwidth IR probe pulses are generated by using 800 nm light and its first two harmonics to ionize nitrogen gas.^27,28 The resulting plasma radiates broadband IR light, which is centered at 2000 cm⁻¹ with a bandwidth of ∼2000 cm⁻¹, is roughly 50 fs long, horizontally polarized, and has an energy below 10 nJ per pulse.

The output of the 3 μm OPA is fed into a Mach-Zehnder interferometer equipped with CaF₂ beam splitters and high-accuracy translation stages (Aerotech ANT-25L). The interferometer generates two pulse pairs, which are separated by an experimentally controlled time delay with a resolution of 0.3 fs. One of the pulse pairs is sent through a tunable zero-order waveplate (Alphalas) and wire-grid polarizer (Specac) before being sent to the sample. The other pulse is sent into a room-temperature mercury-cadmium-telluride (MCT) detector for timing and automatic phasing of 2D IR spectra.²⁹ 

The 2D IR experiments are performed in the pump–probe geometry,^30,31 with the two pump pulses collinear and the probe crossing them in the sample at a small angle. The pulses are focused to a 100 μm spot in the sample with a gold parabolic mirror (EFL = 10 cm). After passing through an analyzing polarizer, the transmitted probe is dispersed by a 50 groove/mm grating onto one stripe of a 2 × 64 pixel MCT array detector (Infrared Systems). A small portion of the probe from before the sample is dispersed onto the other stripe as a reference.

The 2D IR signal is measured in the frequency domain on the array (ω₃) by continuously sweeping³² the delay between the pump pulses (τ₁) for a fixed delay between the fixed arm and the probe pulse (τ₂). Each signal shot is divided in ω₃ by the reference to obtain a change in absorption (ΔOD) between irradiation with both pump pulses and only one pump pulse. In post-processing, for each value of ω₃, the DC background is subtracted. The free induction decay (FID) is first phase-corrected using the Mertz method³³ with the phase calculated from the single-channel detector in the Mach-Zehnder interferometer. The phase-corrected FID is apodized with a Hann window, zero-padded, and Fourier transformed to generate an ω₁ axis. Transient absorption spectra are collected on the same instrument with the movable arm of the interferometer blocked. In this case, the change in absorption in ω₃ is measured as a function of τ₂. The polarization for 2D IR and TA measurements is adjusted by rotating the polarization of the pump-pulse pair to be parallel (S_∥, ZZZZ), perpendicular (S_⊥, ZZYY), or at the magic angle (54.7°) with respect to the analyzing probe polarization. TA polarization anisotropy measurements are calculated from r(τ₂) = (S_∥ − S_⊥)/(S_∥ + 2S_⊥), and 2D depolarization spectra are obtained from ρ(ω₁, ω₃) = S_⊥/S_∥.

Samples were prepared by sandwiching 18.2 MΩ ⋅ cm resistivity water, purified by reverse osmosis, between two 1 mm thick CaF₂ windows. This resulted in an optical density of ∼0.6 in the O–H stretching region, corresponding to a path length of ∼1 μm between the windows. To avoid the nonresonant electronic response of the CaF₂ windows, which arises when the pump and probe are overlapped in time, data were analyzed only for τ₂ ≥ 100 fs.³⁴ The temperature of the sample-cell holder was regulated with a recirculating bath. The temperature of the sample cell was monitored by affixing a type-K thermocouple to the sample cell with a silver compound thermal paste.

2D IR spectra of H₂O, taken at various waiting times in the magic-angle polarization scheme, are shown in Fig. 1. In the early-time 2D IR spectrum, the O–H stretching mode shows a ground-state bleach (GSB) along the diagonal. As a function of ω₃, the width of the GSB reflects that of the linear spectrum; however, the width in ω₁ is limited by the excitation pulse bandwidth (∼380 cm⁻¹). The GSB at 100 fs is elongated along the diagonal, with a center-line slope (CLS) of ∼0.45, which reflects that the excitation has not had sufficient time to randomize its frequency by this time. We note that the center of the bleach of the O–H stretch lies beneath the diagonal, which is contrary to typical 2D IR spectra. This is due to the vibrational non-Condon effect in which strongly hydrogen-bonded (i.e., most red-shifted) oscillators have a larger transition dipole moment than weakly hydrogen bonded ones.^35,36 Shifted to lower frequencies in ω₃ lies the excited state absorption (ESA), which corresponds to the transition from the first excited state to higher lying excited states. Unlike the GSB, it is significantly broader in ω₃ than either the fundamental or overtone transition in the linear absorption spectrum (3405 and 6800 cm⁻¹, respectively), and it even continues past our lowest detection frequency.

In addition to the GSB and ESA in the O–H stretching region of the spectrum, a bleach of the HOH bending mode at (ω₁, ω₃) = (3400 cm⁻¹, 1675 cm⁻¹) is evident, consistent with previous two-color TA measurements^19,20 and mixed IR–Raman experiments.^14,37 The presence of a cross peak between the stretching mode and the bending mode at the earliest waiting times indicates that the two vibrations are coupled. Physically, this implies that the stretching mode has a certain amount of bend character to it, and the excitation of a local stretching mode would result in the coherent exchange of energy between the stretch and the bend. The bend cross peak shows only a weak slope, which means that there is little frequency correlation with the stretching mode. Unlike usual cross peaks, there is no clear ESA peak associated with the bleach, which would normally correspond to the transition to the stretch–bend combination band. We note that the combination band in the linear spectrum peaks at a frequency of ∼5200 cm⁻¹,³⁸ which suggests that the stretch-to-combination transition ought to occur at a frequency of 5200–3400 = 1800 cm⁻¹. This is unusual in that the combination band is normally shifted to lower frequencies than the sum of the fundamentals, and this represents an atypical example of a negative off-diagonal anharmonicity.³⁹ 

Figure 2 shows the subtle changes in the 2D IR spectrum associated with changing the relative polarization between the excitation and probe pulses at a waiting time of τ₂ = 100 fs. The peaks retain the same general features but change considerably in shape. The parallel spectrum appears more elongated along the diagonal compared to the perpendicular spectrum, whereas the latter shows more overall intensity off the diagonal. The difference in off-diagonal intensity can be quantified by the zero-time CLS, which is 0.55 for the parallel spectrum while it is 0.35 for the perpendicular one. In the parallel spectrum, features of the long-time spectrum are evident as early as 50 fs, which are indicated by the negative feature in the upper-left part of the bleach. However, the extra intensity off the diagonal in the perpendicular spectrum suppresses the early-time signatures of the long-time spectrum. In addition, the central peak of the perpendicular spectrum is somewhat more structured than the parallel spectrum, qualitatively reflecting the O–H stretching band observed in the linear spectrum.

The 2D IR spectra as a function of waiting time, shown in Fig. 1, provide insight into the molecular dynamics of the liquid. Figure 3(a) (blue squares) shows the normalized ESA signal obtained by integrating the 2D IR spectrum in the ω₁ = 3100–3700 cm⁻¹ and ω₃ = 2500–3000 cm⁻¹ regions. With increasing τ₂, the broad ESA decays on a τ_ESA = 260 fs time scale, corresponding to the relaxation of the system out of the initially excited state. From Fig. 1, it is evident that the ESA decays more quickly at low ω₁ than at high ω₁. For example, in the τ₂ = 500 fs surface, there is clearly more intensity on the blue side of the ESA while the red side has decayed away. Integrating the ESA for ω₁ < 3400 cm⁻¹ and ω₁ > 3400 cm⁻¹ yields decays characterized by exponential time constants with a value of $τESAR=220$ fs on the red side and $τESAB=310$ fs on the blue side.

As frequency memory is lost with increasing τ₂, the initially elongated GSB becomes symmetric about ω₁, which we characterize by the decay in CLS.^40,41 While there are many ways of quantifying the loss in frequency memory, for H₂O, the CLS is best since it is least susceptible to the growth of the long-time spectrum which grows in without a slope, though small interference effects will be present in the decay time scale obtained. The CLS for a given value of τ₂ is calculated by fitting the first moment of the bleach to a straight line. Fitting the CLS, shown in Fig. 3(b), as a function of τ₂ yields a decay time scale of τ_CLS = 165 fs, which implies that the environment of an excitation is randomized on this time scale. Calculating the CLS decay for the red side or the blue side alone does not yield significantly different values, suggesting that the measured CLS decay is due primarily to energy exchange and not hydrogen-bond reorganization.

After τ₂ ≈ 1 ps, the 2D IR spectrum is essentially unchanged with increasing waiting time and features in ω₁ and ω₃ are entirely uncorrelated (see Fig. 1). In previous literature,^15,34,42 this long-time spectrum has been attributed to the heating of the liquid upon rapid relaxation of the excitation into intermolecular modes, and as such, has been referred to as the “hot ground state” (HGS). In the HGS spectrum, we observe an overall blue shift of the stretching mode along ω₃, implying that hydrogen bonds have weakened compared to the initial spectrum. Though there is interference from the early-time GSB, integrating the blue-shifted absorption in the region ω₁ = 3100–3700 cm⁻¹ and ω₃ = 3600–3800 cm⁻¹ provides a time scale for the growth of the HGS of τ_HGS = 700 fs. The bend is also red-shifted and narrowed, which is consistent with weaker hydrogen bonds. These features are in accord with what is observed in a heated linear difference spectrum³⁸ and so has led to the natural interpretation that the sample is at a higher temperature than at the outset of the experiment. Indeed, the projection of the 2D IR spectrum onto the ω₃ axis qualitatively matches a linear thermal difference spectrum (see Fig. 4(b)). For consistency with the literature, we will use the term “HGS” throughout, although we emphasize that the τ₂ = 5 ps HGS spectrum is not an equilibrium system, as discussed below.

Keeping in mind that the 2D IR spectrum projected onto ω₃ axis is equivalent to the TA spectrum for the same waiting time, the magic-angle TA surface of H₂O, shown in Fig. 4(a), reveals the continuous dynamical evolution of the nonlinear spectrum. A slice through the τ₂ = 100 fs TA spectrum of H₂O (Fig. 4(b), blue curve) shows a GSB of the fundamental O–H stretch and its associated ESA, which is as broad as it is in the 2D IR spectrum. Sitting atop the ESA at 1675 cm⁻¹ is the cross peak to the bending mode. The TA spectrum continuously evolves to the HGS spectrum as time increases (Fig. 4(b), red curve). The HGS spectrum in Fig. 4(b) can be compared to a linear thermal difference spectrum (Fig. 4(b), black curve) for a temperature change of ΔT = 15 K. The linear thermal difference spectrum qualitatively matches the HGS spectrum over the entire mid-IR though there are some slight differences in amplitude and linewidth.

As in the 2D IR spectrum, we integrate the ESA over the 2500–3000 cm⁻¹ region to extract τ_ESA and over the 3600–3800 cm⁻¹ region to extract τ_HGS. The evolution of these features is shown in Fig. 3(a) (blue and red dots, respectively), and corresponding integration patches are shown in Fig. 4(b). Integration of the ESA yields a biexponential decay with time constants of 250 and 700 fs. We attribute the slow time scale to the growth of the HGS, while we associate the fast time scale with τ_ESA. Integration of the HGS feature yields a decay of 610 fs. These values are in good agreement with previous TA measurements^15,17,19 and deviate only slightly from the values obtained from 2D IR. Previously, these time scales were extracted by singular value decomposition of the TA spectrum²⁵ yielding τ_ESA = 275 fs and τ_HGS = 720 fs. The differences between the direct fitting of the experimental data and the fitting of the components of the singular value decomposition are minor and well within the experimental noise.

Polarization anisotropy measurements, r(τ₂), report on the orientational correlation function of the excited transition dipole moment. The τ₂ = 0 value of r reports on the angle between the transition dipole of the two vibrations being excited and probed, varying from r(0) = 0.4 for parallel transition dipoles to r(0) = − 0.2 for perpendicular transition dipoles.⁴³ The polarization anisotropy for the O–H stretch, as measured by integrating the TA signals over the 3300–3500 cm⁻¹ window, is shown in Fig. 5(a) (blue dots). An exponential fit (Fig. 5(a), black curve) yields an amplitude of r(0) = 0.39 and an orientational relaxation time of τ_r = 70 fs. The time scale is consistent with other measurements in the literature.¹⁷ This is essentially the fastest time scale that can be resolved by our instrument and so it should be taken as an upper bound. However, a time scale of 70 fs corresponds to a frequency of ∼500 cm⁻¹, which is consistent with the librational motion in H₂O. Polarization measurements revealed that the τ₂ = 100 fs value of the CLS decay varied from ∼0.40 to 0.25 from parallel to perpendicular probing, but τ_CLS did not change.

The cross-anisotropy measurement between the stretch and the bend, which reports on the correlation between the angle of the transition dipole of the stretch and the bend, is shown in Fig. 5(a) (red dots). Since the bend cross peak sits atop the broad ESA, it is necessary to correct for this broad background. Practically, it is found that the cross-peak anisotropy decay is insensitive to the correction, and a linear correction was used. Unlike the diagonal anisotropy measurement, the stretch-bend anisotropy begins at a small value and rapidly decays to zero in a seemingly oscillatory fashion, which we attribute to non-resonant response of the CaF₂ windows in this region. This implies that any correlation between the transition dipole of the stretch excitation and the bend is both weak and short lived.

In addition to the simple polarization anisotropy measurement, we may perform a 2D depolarization measurement (ρ = S_⊥/S_∥), which reports on orientational correlations in a frequency dependent manner.^44,45 In the absence of vibrational relaxation, a value of ρ = 1/3 indicates that excited and probed transition dipoles remain aligned, while ρ = 1 indicates that correlations have decayed. Figure 5(b) shows the 2D depolarization surface taken at a waiting time of τ₂ = 100 fs. Tracking the behavior along the diagonal, the center of the band (ω = 3400 cm⁻¹) shows ρ ≈ 0.8, while the edges show a depolarization that is significantly lower. The center-band depolarization of 0.8 is consistent with the polarization anisotropy measurement of ∼0.1 at τ₂ = 100 fs shown in Fig. 5(a). At the red edge of the spectrum ρ ≈ 0.33, indicating that orientational correlations have not decayed significantly for the most red-shifted oscillators. Oscillators on the blue side also show a lower depolarization ratio; however, it is larger than that of the oscillators on the red side, which demonstrates that the most red-shifted oscillators retain the strongest correlations.

Changing the temperature of the sample tunes the interactions between water molecules by changing the average hydrogen-bond strength. The temperature-dependent relaxation time scales of the system are therefore a direct reporter on how intermolecular interactions affect the dynamics of the system. Over 1300 2D IR spectra were measured as a function of waiting time for temperatures between 265 and 325 K and analyzed for four vibrational and orientational relaxation parameters (τ_ESA, τ_HGS, τ_CLS, τ_r) at each temperature. Representative 2D IR spectra are presented in the supplementary material and representative curves for the ESA, HGS, and CLS, extracted as described above, are shown in Figs. 6(a)–6(c), respectively. In addition, polarization anisotropy measurements from TA spectra were also measured. It is worth noting that with such a large data set taken continuously over the course of several days, the signal-to-noise ratio on any given spectrum is not particularly large, fluctuations in laser power caused drifting of the overall signal, and there were issues with pump scatter at high temperature. However, the effective averaging that is a result of extracting time scales results in reliable curves and time scales as evidenced by the data quality in Fig. 6. Figure 7 shows the dynamic observables, τ_ESA, τ_HGS, τ_CLS, and τ_r, as a function of temperature. Strikingly, no significant dependence of the lifetime of the O–H stretch, τ_ESA, is seen with temperature over the range studied (Fig. 7(a), blue circles); almost all of the points lie within a 14 fs standard deviation centered about a mean value of 260 fs. The growth of the HGS, τ_HGS, shows more significant differences upon increasing temperature (Fig. 7(a), red diamonds). Qualitatively, a weak dependence is apparent, with the HGS growing in more slowly at low temperatures and more quickly at temperatures up to ∼295 K at which point, the growth seems to slow down again. The HGS shows a mean of 740 fs across the entire range but shows a minimum value of about 600 fs.

The decay of the CLS (Fig. 7(b)) also shows no significant temperature dependence and with an average value of τ_CLS = 170 fs and a standard deviation of 16 fs. This is in contrast with previous studies, which found that spectral diffusion occurs more slowly at 274 K compared to 304 K.¹⁸ On the other hand, the polarization anisotropy decay shows a modest increase in τ_r as a function of temperature (Fig. 7(c)), increasing from ∼80 fs at our lowest temperature to ∼160 fs at our highest with a slope of ∼0.85 fs/K. While the anisotropy measurement is inherently quite noisy, this increase is outside the range of statistical error. This trend is at odds with intuition in which we would expect to lose orientational correlations more rapidly at higher temperature.

Heating evidently plays an important role in nonlinear spectroscopy, since ∼1 ps after excitation, the spectrum resembles that of a heated equilibrium liquid. While the correspondence between the thermal difference spectrum and the HGS spectrum, shown in Fig. 4(b), has been consistently made in the literature for the O–H stretching region, we emphasize here that the correspondence holds across the entire mid-IR spectrum. To correlate the spectral changes in nonlinear spectra with equilibrium temperature differences, a series of linear temperature difference spectra were measured between 265 and 285 K using the same instrument and sample as the nonlinear spectra. For each temperature, the transmission of the probe pulse (with pump pulses blocked) was measured. The magnitude of each difference spectrum increases with increasing temperature difference while the shape of each spectrum remains roughly the same.

To quantify the changes, we empirically found that the magnitude of the differential absorption change in the O–H stretching region δ(T_i, T_f) scales linearly with the temperature difference ΔT = T_f − T_i and is independent of initial temperature, as illustrated for various T_i in Fig. 8. A linear fit to the data results in a temperature difference that scales as ΔT = δ/6.3 K. For this plot, we did not correct δ for sample path length (ℓ), since an accurate determination is difficult. However, one can estimate the path length as ℓ = 1.1 μm using the optical density of 0.6 and the 25 °C O–H stretch molar extinction coefficient of 100.6 M⁻¹ cm⁻¹.^46,47 Then an approximate expression for the temperature change that accounts for path length is ΔT = δ/5.7ℓK μm.

Under the assumption that the energy has thermalized at long waiting times, this linear fit provides a means for quantifying the temperature increase in a nonlinear IR experiment. In typical experiments, we find that δ ranges between 45 and 100, depending on the pump characteristics (i.e., energy, center wavelength, bandwidth), corresponding to temperature changes ranging from 7 to 16 K. For example, for the thermal difference spectrum shown in the inset of Fig. 8, we can estimate a temperature rise of ∼16 K after excitation. These values are significantly larger than the values of ∼3 K which are typically reported in the nonlinear spectroscopy of dilute HOD;⁴² however, this can easily be rationalized by the fact that the excitation density is significantly higher in H₂O than in HOD. For comparison, we may also predict the temperature jump from energy conservation. Using the overlap of the water absorption with the spectral energy density of our pump pulses, the excitation volume in our sample and the heat capacity and density of water, we calculate a temperature rise of 15.0 K, which in good agreement with the linear calibration.

A comparison to other hydrogen-bonding liquids⁴⁸ shows that the decay of the CLS in water is a factor of five faster than other hydrogen-bonding liquids, such as methanol (τ_CLS = 1.08 ps, see supplementary material for details and a discussion of viscosity effects). In water, hydrogen bonds are formed in a geometry that is tetrahedral on average, allowing for the formations of extended networks. On the other hand, the ability to form an extended network in methanol is limited by the crowding due to the methyl groups. We must conclude that water’s rapid relaxation dynamics stem from its ability to form such networks of resonant oscillators, which indicates that vibrations are delocalized over multiple molecules. The evolution of these vibrational excitons gives rise to the rapid relaxation observed. That the dynamics we observe are largely dominated by excitonic effects is a testament to the strong intermolecular coupling resulting from the hydrogen-bonding interaction.

The vibrational eigenstates are primarily determined by the instantaneous coupling between local O–H stretching modes. If the intramolecular coupling was much larger than the intermolecular coupling, then vibrations would be mostly localized to a single molecule. However, we must consider that the eigenstates of the system are excitonic. In other words, strong intermolecular coupling, effected by the hydrogen-bonding interaction, directly results in the eigenstates being delocalized. Indeed, calculations have shown that in the local-mode basis both the inter and intramolecular coupling are comparable,²² having a value of about −30 cm⁻¹ with complex eigenstates spread over ∼400 cm⁻¹ of bandwidth; furthermore, even the vibrations in the symmetric/antisymmetric basis show a broad distribution of coupling constants.⁴⁹ The details of the strength of the coupling between local O–H vibrations is determined by the NPES through the local hydrogen-bond configuration, and its large anharmonicity is essential for strong intermolecular coupling.

The 2D depolarization measurement provides further insight into the delocalization of the O–H stretching vibrations through spectral resolution. If we associate low initial orientational correlation with vibrational delocalization, the depolarization surface provides a metric for the extent of delocalization as a function of frequency. In particular, stretching vibrations at the center of the band are the most delocalized, while vibrations at the edges of the O–H stretching band (in particular the red edge) are most localized. These results are consistent with previous theoretical predictions²² in which the most strongly hydrogen-bonded O–H oscillators were found to be the most localized. In this study, vibrations at the center of the band, 3400 cm⁻¹, were delocalized over up to six molecules, while those at 3150 cm⁻¹ were localized to a single molecule. Going forward, the 2D depolarization measurement is a promising avenue for experimentally characterizing the evolution of initially delocalized of excitonic vibrations in water in a time-dependent manner.

The anharmonic interaction between the excitonic stretching vibration and the HOH bending mode is evidenced by the stretch–bend cross peak in the 2D IR spectrum. Traditionally, these couplings are discussed in terms of an anharmonic expansion of the local mode stretch-bend interaction potential^50–52 with even-order terms coupling vibrations within an excitation manifold and odd-order terms coupling adjacent excitation manifolds. In particular, the Fermi-resonance interaction^53–55 that originates from the ∼2:1 frequency ratio between the stretch and bend overtone originates from the cubic and other odd-order terms and is commonly observed for hydride vibrations such as C–H,^56–58 N–H,^59–61 and O–H.^62,63 In liquid H₂O, such an expansion is not particularly useful since the highly anharmonic potential cannot be simplified to a small number of degrees of freedom without sacrificing crucial details of the system. The excitonic nature of vibrations in liquid water necessitates the consideration of many coordinates. Furthermore, in H₂O, we do not expect a “classic” two-peak signature of a Fermi resonance in which the bend overtone gains intensity and shifts from the stretch.⁶² Rather, the bend overtone lies entirely in the breadth of the stretch fundamental, and so the shifting most likely appears as broadening.

The presence of the bend cross peak indicates that the eigenstates at 3400 cm⁻¹ have a character which is an admixture of stretch and bend. This cross peak does not show significant frequency correlation with the O–H stretch, which is expected since the frequency of the bend is determined by the hydrogen-bonding configuration of both O–H oscillators, resulting in a weak correlation. Only a bleach is evident for the bend cross peak; however, upon excitation of the stretching mode, we expect to observe the stretch-to-combination transition, which would correspond to the ESA of the bend cross peak. The reason for this absence is unclear, but it is most likely rooted in the fact that ω_s ≈ 2ω_b in H₂O and that excitation of the O–H stretch necessarily results in either the direct or indirect excitation of the bend overtone. Since we excite the bend overtone along with the stretch fundamental, we expect to observe the ν = 2 → 3 ESA and the ν = 2 → 1 bleach of the bend. The strong anharmonicity of the potential also means that two-quantum transitions are not forbidden, and so we expect a bleach at the difference frequency between the stretch and bend, which arises from stimulated emission from the former to the latter. All of these features are expected to be at roughly the same frequency, broad, and highly overlapping. It is therefore unsurprising that there is no clear ESA associated with the bend cross-peak.

Perhaps the most striking feature of the early-time 2D IR spectrum of H₂O is the broad ESA. ESA peaks normally correspond to features in the linear overtone (near-IR) spectrum; however, while the near-IR spectrum of liquid H₂O has little optical density at ∼6000 cm⁻¹,^38,64 the ESA easily spans transitions covering this region. Furthermore, simulations of the 2D IR spectrum of liquid water capture the width of the GSB but cannot reproduce the width of the ESA.^21,65 If line-broadening in ω₃ is due to pure dephasing or vibrational relaxation, harmonic scaling arguments can be used to predict that the ESA will be a factor of 1 or 2^1/2 broader than the GSB, respectively.⁶⁶ For H₂O, the ESA is far broader still, which implies that states are accessed in the 2D IR experiment that are not in the linear experiment.

To understand this effect, we consider the NPES and its relation to linear and nonlinear spectroscopy. In linear spectroscopy, transitions originate from the ground state of the high-frequency mode, while in 2D IR spectroscopy, the final transition originates from the first excited state. The ground state experiences a more harmonic part of the NPES since displacements from the minimum are on average small at low energy. On the other hand, in the first excited state, the system samples a highly anharmonic part of the potential since displacements from the minimum are much larger. A transition originating in a more anharmonic state implies stronger coupling between modes (or equivalently, weakened selection rules) such that there are more transitions with nonzero oscillator strength available. At this level of excitation then, it stands to reason that the modes are mixed enough such that it is no longer appropriate to identify these features as arising solely from O–H stretching-type motion. The effects of this extreme anharmonicity have been discussed previously in the context of ice^67,68 and was used to explain the qualitative differences between H₂O and D₂O.⁶⁹ 

This anharmonicity has its origins in the hydrogen-bonding interaction, which significantly modifies the NPES. An example of this is seen in the 2D IR spectrum of hydrogen-bonded N-methylacetamide (NMA) dimers in dichloromethane, reported in Ref. 61. In the non-hydrogen-bonded monomer form, the N–H stretching mode has a narrow ESA that is consistent with pure dephasing. On the other hand, the ESA of the hydrogen-bonded N–H mode is several times broader than the GSB. This shows that the formation of a single hydrogen-bond results in a deformation of the NPES significant enough to greatly alter the ESA. If this effect is evident even for a single hydrogen bond, certainly it can explain the extreme broadening in H₂O, which has an extended hydrogen-bond network. We therefore interpret the shape of the 2D IR spectrum of water in terms of the anharmonicity of the NPES induced by the strong hydrogen-bonding interactions.

The changes in the O–H stretching line shape between S_∥ and S_⊥ are telling as well. The increased off-diagonal intensity in the perpendicular spectrum suggests the presence of cross peaks beneath the O–H stretching band between overlapping transitions. The origin of the structure observed within the broad O–H stretching resonance has been interpreted alternately as arising from distinct hydrogen bonding environments within the liquid,^70–72 symmetric and asymmetric vibrational modes,^73,74 coherent vibrational energy transfer,⁷⁵ and Fermi resonance between the stretch and bend.^37,76 The presence of cross peaks in short waiting time spectra is inconsistent with the possibility that the substructure comes from distinct liquid environments. That the off-diagonal intensity grows in the crossed-polarization scheme suggests that the substructure under the O–H stretch resonance comes from the same group of molecules. We therefore conclude that the most red-shifted feature under the O–H stretching band is indeed due to the overtone of the bending mode. Our recent analysis of the polarization dependence of the 2D IR spectrum of liquid D₂O revealed symmetric and antisymmetric O–D stretching modes as well as a Fermi resonance with the bending mode.⁶⁹ In the case of H₂O, there is no evidence for symmetric and anti-symmetric stretching, since the depolarization does not take on values greater than 1.

Although the notion of intermolecular delocalized O–H stretching vibrations is rather old,^1,76–79 the language of vibrational excitons to describe liquid water is relatively new.^{18,21,23,25,80–83} Indeed, the exciton language seems to have arisen of a need to describe dynamical effects in the vibrations of water, which are determined by the superposition of eigenstates in the liquid. As such, complexity in interpreting experiments arises from the very nature of the initially prepared state. The IR spectroscopy of H₂O cannot be described in terms of dilute non-interacting chromophores. Rather, all highly coupled O–H oscillators within the excitation volume are simultaneously driven by a coherent external electromagnetic field. Thus the specific temporal and spectral profile will influence the nature of the excitation that is prepared. Broadband, short-pulse lasers generate a coherent superposition of the ground state and many excited O–H stretching eigenstates, each with an amplitude and phase reflecting that of the driving laser field.⁸⁴ Therefore, the excited vibration can really be anything (in principle, even localized to a single bond) and the exact nature of the initial state is going to depend intimately on the driving field. With increasing time, the initial coherent superposition of eigenstates evolves, resulting in rapid propagation of the excitation. Given the density of eigenstates, the odds of any recurrence of the initial superposition are exceedingly small. We therefore imagine that, over short time scales, the vibrations change significantly and irreversibly in character. This has ramifications in how one should compare the results of experiments taken with different laser sources that vary in spectrum and pulse length. Experiments taken with longer pulses^15,17,19 or pulses with a nonlinear phase necessarily generate initial states that are distributed over a narrower range of eigenstates, which evolve differently than superpositions over a broad range of eigenstates.

This picture of the molecular vibrations in water raises the interesting question of how to interpret delocalized excitations that evolve over stochastically fluctuating molecules. If the fluctuations of the water molecules are much slower than the time scale defined by the energy of the exciton, then the evolution ought to be adiabatic and the excitonic vibration evolves continuously with the low-frequency motion of the molecules.⁸⁵ If the fluctuations were much faster than the frequency of the excitation, this would localize the vibrations and result in incoherent dynamics.^86,87 This has been one of the central debates in the vibrational spectroscopy of water, albeit the question is usually framed somewhat differently. Often, vibrations are described in one of two pictures: excitations in water are delocalized and evolve coherently as water molecules move^20,21,25 or excitations are entirely localized and evolve incoherently.¹² Here, we emphasize that these two interpretations are actually at the extremes of a spectrum and one transfers smoothly into the other. The first interpretation assumes that the vibrations localize on time scales long compared to the experiment, while the second assumes that the vibrations localize on time scales that are short compared to the instrument response function. In other words, the question is not whether the dynamics are incoherent or coherent, but rather on what time scale they change from the latter to the former. Considering the time scales associated with low-frequency motions in water (i.e., ∼30 fs for high-frequency librations to ∼300 fs for hydrogen-bonding modes), we conclude that nonlinear IR experiments are in an intermediate regime which is neither completely coherent nor incoherent. An initial excitation evolves coherently until it eventually localizes and evolves incoherently. This picture provides an intuitive explanation for why the dynamics in water are so rapid and why a non-adiabatic description is necessary for understanding the relaxation.

Above it was noted that the red side of the O–H stretching band relaxes slightly faster than does the blue side. As observed in other hydrogen-bonding systems,^15,88–90 oscillators in strongly hydrogen-bonded configurations relax more quickly than those in weak or non-hydrogen-bonded configurations due to the increased density of LFMs. Though the correlation is not as strong as in other systems, it is expected that the red side of the O–H stretching band is going to correspond to more strongly hydrogen-bonded configurations on average, as well as excitons that are strongly mixed with the bend overtone. The overtone of the bend provides, in principle, yet another efficient path for relaxation for the most red-shifted oscillators. This dispersive relaxation was recently attributed⁹¹ to long-lived heterogeneity in liquid water; however, this is inconsistent with all previous measurements of spectral diffusion in H₂O. To understand the frequency-dependent relaxation rate, we must consider that the lifetime reports on a combination of the fluctuations of a particular configuration as well as spectral diffusion into and out of that configuration. Therefore, rapid spectral diffusion is entirely consistent with dispersive relaxation without the need to invoke persistent heterogeneity. That the integrated value of the O–H stretch relaxation, τ_ESA, is closer to $τESAR$ than $τESAB$ can be understood in terms of the vibrational non-Condon effect.^35,36 The transition dipoles for red oscillators are on average larger, meaning that these oscillators contribute more significantly to the measured signal than blue oscillators and so the average relaxation is closer to that of the red side than the blue side.

The fact that the HGS grows in more than twice as slowly as the O–H stretch relaxes suggests that there is some state through which the system relaxes. If this state was the bend, relaxation would proceed via a stretch to bend to LFM mechanism characterized by time scales τ_ESA and τ_HGS, respectively, in which population would build up in the bend state as evidenced by an increasing bleach of the bending-mode cross peak. Previous TA measurements, using narrower bandwidth pulses, concluded that this is the mechanism of relaxation.^19,20 However, we do not observe this in broad bandwidth experiments. Rather, the bend cross-peak bleach disappears uniformly on a time scale similar to τ_ESA as the HGS grows in. Due to the strong overlap of the early- and late-time spectra in this region, it is difficult to extract a precise time scale. The data therefore suggest that any intermediate state is outside the window of direct detection. This interpretation is corroborated by the fact that population relaxation out of the bending mode occurs on a 170 fs time scale,^92,93 which is significantly faster than τ_HGS. That is, the bend cannot provide the observed delay between the relaxation of the O–H stretch and growth of the HGS.

The polarization anisotropy measurement shows that on the fastest time scale that molecules can move at room temperature, orientational correlation of the O–H stretch is completely lost. Clearly, the correlation loss is not due to physical reorientation of the molecules but rather to the evolution of the initial excitation. A comparison of H₂O to isotopically dilute HOD shows that while HOD does have a small amplitude 70 fs component in the polarization anisotropy decay, it is largely lost on a 2.7 ps time scale.⁴⁴ The former component is attributed to rapid librational modes and the latter component is attributed to large-scale reorientation of the liquid, which brings the system to an isotropic distribution. This component is absent in H₂O, which means that the rapid librations in the liquid suffice to remove all orientational correlations, consistent with the notion of librations dramatically changing the initial excitation. In other words, the first inertial librational motions of the molecules making up an exciton result in a rapid and significant spatial and orientational evolution of the initial excitation. In fact, the measured anisotropy decay time scale is consistent with the predicted loss of anisotropy of ∼115 fs in the high-frequency librational modes.⁹⁴ 

The cross-anisotropy between the stretch and the bend is more telling still. For a cross peak in the narrowband pump—broadband probe scheme, the fact that the cross-anisotropy between the stretch and bend is essentially zero for all times can be explained in two ways: either the angle between the stretch and the bend transition dipoles is rigidly fixed at Θ = 54.7° which corresponds to r = 0 or there are no correlations between the stretch and bend transition dipoles.⁴³ Given the highly uncorrelated transition-dipole angles in liquid water,^49,69,95 we are led to conclude that there is no correlation between the orientation of the stretching mode and that of the bending mode. This is also consistent with the notion of vibrational excitons. If the initially excited vibrations are delocalized over several molecules, giving rise to a collective O–H stretching mode, then there is no expected correlation between it and the local bending mode of the individual molecules.

Intuitively, we expect that increasing temperature increases molecular fluctuations and results in faster decays. Therefore, it is unexpected that the O–H stretch lifetime and CLS decay time scales do not strongly depend on temperature. A comparison to dilute HOD in H₂O shows that the CLS decay changes by at least a factor of two (from 2.4 to 1.1 ps) over a similar temperature range⁹⁶ while water remains constant (at ∼170 fs). Therefore, we conclude that this lack of temperature dependence is an excitonic effect in which the intermolecular couplings are too strong and the exciton dynamics too fast to be affected significantly by changes in thermal occupation of the intermolecular modes. The dynamics of the exciton itself, dictated by its collective nature and the strong coupling between high- and low-frequency modes, would result in the rapid decay of the ESA and CLS.

Increasing the temperature results in the polarization anisotropy decaying more slowly, which is contrary to typical systems where orientational correlations decay more quickly due to increased fluctuations. However, excitons are expected to be more localized at higher temperature due to weaker intermolecular coupling, and so orientational correlations persist longer, resulting in an increase in the anisotropy decay time. This is a clear signature of excitonic vibrations. It is worth noting that even at our highest temperature, we do not observe large-scale reorientation of the hydrogen-bonding network, suggesting that the vibrations are still quite delocalized.

The only time scale that varies appreciably with temperature is the growth of the HGS, which grows in 100–200 fs slower at low temperatures compared to room temperature. As temperature is increased, the fluctuations of LFMs become larger and as such are easier to populate. This results in a speedup of the relaxation. The ostensible slowdown of the HGS growth at high temperatures can be understood in terms of the weakening of hydrogen bonds. Since the hydrogen bonds are weaker at high temperature, the density of LFMs is reduced, and the relaxation takes longer.

The HGS, characteristic of the nonlinear spectroscopy of all hydrogen-bonding liquids, is normally attributed to heating of the system due to the initial excitation.^{15,19,20,34,42} The conventional wisdom is that the energy, initially in high-frequency modes, relaxes to LFMs that thermalize on the picosecond time scale, resulting in a system at higher temperature. This interpretation is corroborated by the fact that the HGS spectrum qualitatively resembles a linear thermal-difference spectrum. Subsequently, thermal diffusion brings the excitation volume back to the initial temperature in an estimated time of ∼20 ns based on the heat equation.

Typically, when we refer to the equilibrium temperature of a fully thermalized system, we mean that the population of the eigenstates of the system has reached a Boltzmann distribution. Clearly, the system cannot reach an equilibrium distribution on the time scale of our experiment and temperature is a poorly defined concept. Naively, the growth of the HGS on a time scale of τ_HGS = 700 fs would correspond to a heating rate of ∼10¹³ K/s. This rate of heating, being reminiscent of temperature changes in the early universe, is clearly unphysical. Furthermore, given the speed of sound and rate of thermal diffusion in water, it is in no way possible for hydrodynamic modes to equilibrate on these time scales. We must necessarily adhere to the notion that in the 1 to ∼100 ps after the pump pulse passes through the sample, it is in a highly non-equilibrium state in which the population of eigenstates is not a Boltzmann distribution and as such requires careful consideration.

This implies that the HGS arises entirely due to populating LFMs and that we are spectroscopically insensitive to the other processes that accompany the march to equilibrium, such as changes in density. Indeed, no change in the trend of the linear thermal difference spectrum is observed as the temperature is scanned through the density maximum at 277 K.^70,97 We conclude that only the population of the LFMs affects the spectrum of O–H stretch and that we observe “heating” when the system is out of equilibrium is a result of these LFMs being populated, regardless of whether the energy is delivered thermally or via laser pulse.

This behavior is consistent with the recently developed notion of prethermalization,^98,99 which is a universal feature of non-equilibrium systems. In prethermalization, the equation of state and kinetic temperature reach constant values on time scales much shorter than the equilibration time scale in an approximate equipartition of kinetic and potential energy. Nonetheless, several important quantities (e.g., occupation number for LFMs) may resemble their equilibrium values despite being in a non-equilibrium state that can persist for long times. It is therefore perhaps more reasonable to associate the observed 15 K temperature increase with a kinetic temperature, related to the distribution of RMS momentum of the mass weighted intermolecular coordinates, and not a thermal one. That is, as the O–H stretch relaxes, the kinetic energy of the LFMs that dictate the spectrum of the high-frequency modes almost instantly reaches a value associated with a temperature increase of 15 K; however, the kinetic and potential energy of other modes is still far from their equilibrium distribution. It is interesting to note that our observations on water are similar to those that are encountered in non-thermal phase transitions in semiconductors induced by femtosecond pulse excitation in which the phase transition appears to occur long before phonon emission can heat the lattice.^100,101

Given that within less than 200 fs, essentially the fastest time scale on which the liquid can move at room temperature, signatures of the HGS are evident, we conclude that the excitation of the O–H stretching mode directly results in a change in the LFMs. If the language of a Born-Oppenheimer (BO) separation between O–H stretch and LFMs were to be used, this would correspond to a non-Condon transition (i.e., a non-vertical transitions between adiabatic potential energy curves).⁸⁵ Since the signatures of the HGS, present almost immediately in the spectrum, imply a weakening of hydrogen bonds, we expect that the excitation of the O–H stretch transition results in an increase of the hydrogen-bond length. This is a large deviation from how spectra are normally interpreted, where the fact that the environment does not change upon excitation is usually an implicit assumption. In reality, this is an indication that the BO basis is a poor approximation to the true eigenstates and that the high- and low-frequency modes are strongly coupled. Therefore, great care must be taken when modeling the vibrations of water as it is inherently a multidimensional problem in which a useful separation of time scales is not possible.

The funneling of energy from the excited state to the ground state within a few cycles of the low-frequency tuning coordinates is indicative of vibrational conical intersections,¹⁰² where a true crossing of adiabatic potential energy curves results in ultrafast relaxation. Indeed, it has been theoretically shown that relaxation rates as fast as 60 fs are possible at a conical intersection in a H₂O–HCO₂⁻ dimer including only three high-frequency modes and three LFMs.¹⁰³ Furthermore, such vibrational conical intersections have been shown to play an important role in the relaxation of the water dimer.¹⁰⁴ Given the enormous number of modes in pure liquid water, there is no reason to believe that such vibrational conical intersections do not play a significant role in the relaxation dynamics of liquid water.

We emphasize that this mechanism of relaxation is notably different from that which has been well accepted in the literature. Previously, studies concluded that relaxation occurred through a cascade of weakly coupled normal modes;^19,20,37,91 however, the evidence for non-adiabatic relaxation presented here is in direct contradiction to these interpretations. Recent computational works are also finding that a simple cascaded mechanism does not hold.¹⁰⁵ The ability to look beyond a narrow spectral window shows that the dynamics are far more complex than a simple explanation would suggest and necessarily forces us to go beyond our standard interpretation when trying to understand how liquid water evolves in time.

Though there are already numerous theoretical studies of collective vibrational excitations in H₂O, to describe the dynamics of vibrational excitons we must scrutinize simplifying theoretical assumptions in addition to following several avenues of new work. The development and testing of models is quite advanced, and the O–H stretch Frenkel exciton Hamiltonian, with transition dipole coupling, appears to perform reasonably well at reproducing the IR spectrum at one quantum of excitation;^{21–23,49,65,78,82,95} however, they inherently decouple the O–H stretching vibration from other degrees of freedom. Additional effort is needed to treat anharmonic effects such as Fermi resonance⁷⁶ or higher level excitations in a manner consistent with experimental observations and to test the validity of the transition dipole approximation.¹⁰⁶ 

Even within existing exciton models, there has been limited analysis of the form and dynamics of the predicted exciton states. Beyond the static inverse participation ratio, there are still many questions about the nature of vibrational excitations in water that could benefit from a time-dependent analysis.^107,108 Are there unique structural characteristics to the excitations that depend on hydrogen bond topology? How does their form evolve with time? Are the dynamics diffusive and how far do they propagate? Models that explicitly incorporate additional inter- and intra-molecular degrees of freedom would be desirable to properly address fast energy redistribution processes^94,105 and will provide the opportunity to address the validity of the weak coupling and adiabatic assumptions that underlie most vibrational population relaxation studies. Ab initio molecular dynamics and approximate quantum dynamics methods also can largely reproduce the full IR absorption^{107,109–115} but are not easily extended to calculate nonlinear and time-resolved IR spectroscopy. However, extending these simulation methods to characterize the dynamics, nonequilibrium excitations similar to those created by experiments could provide the fully atomistic insight into the dynamics of vibrations in water.

Finally, we note that most spectroscopic calculations assume the linear response relationship that relates the calculation of the IR spectrum to dipole correlation functions. The fast dissipation of energy and the unusual variation of vibrational dynamics at 0, 1, and 2 quanta of excitation argue that this may be a poor approximation. Although this may be a minor effect in linear IR spectra, the rapid energy relaxation that is currently handled in an ad hoc manner has a large influence on nonlinear IR experiments. Linear and nonlinear response theory assumptions could be tested by direct comparisons to non-equilibrium spectral calculations in which the IR fields are explicitly included.^94,116 Such calculations could also address how the form of the electromagnetic field influences the nature of the vibrational excitation.

To summarize, the vibrational dynamics of liquid water are determined by the anharmonic NPES, which results in intermolecular couplings comparable to intramolecular couplings and enables strong mixing between high- and low-frequency modes. We emphasize that the strong deformation of the NPES is an effect of the extended hydrogen-bond network in liquid water. This strong coupling results in eigenstates at 3400 cm⁻¹ that are excitonic in nature and are delocalized over multiple water molecules; while they are predominantly O–H stretch in character, they contain a significant amount of bending, librational, and hydrogen-bond-stretching motion. This strong mode mixing results in an extremely broadened ESA.

In an experiment, a coherent superposition of exciton eigenstates is excited, and the nature of the excited state depends intimately on the spectrum and phase of the excitation pulse. This superposition evolves coherently for some time during which the initial excitation is distributed over many eigenstates. As the molecules constituting the delocalized exciton fluctuate, the vibration is localized after which it evolves incoherently. This combination of coherent and incoherent dynamics results in ultrafast relaxation of the various correlation functions that nonlinear IR spectroscopy is sensitive to on time scales of 100–250 fs —the fastest thermal motions in the liquid. Since the relaxation rates are comparable to the rates of thermal fluctuations, they are almost entirely insensitive to the temperature.

After the initial relaxation on a 250 fs time scale, there is a growth of a HGS signal on a 700 fs time scale that results from fast non-adiabatic dynamics, which we believe are effected by vibrational conical intersections expected to occur for a large number of hydrogen-bond configurations. The resulting prethermalized state qualitatively resembles a linear thermal difference spectrum with ΔT = 15 K, although this is a non-equilibrium state, better described by an effective kinetic temperature increase as a result of populating LFMs.

This picture of vibrational dynamics in H₂O has important consequences for solutes in water that should influence our understanding of aqueous chemical reactivity and the role of water in biological processes. Given that most aqueous solutes present hydrogen-bond donor and acceptor sites to water, and that hydrogen bonding is the origin for the observed intermolecular coupling, our perspective suggests that water’s influence on the high frequency solute vibrations and the energy transport between solute and solvent should be considered collectively. Polar solvation in water has been described as such for many years, though it has primarily been in the context of the dielectric response of the liquid and librations. We conclude that these collective dynamics, which result from the highly anharmonic non-adiabatic nature of the donor–acceptor interaction potential, extend to high-frequency modes and may be collective enough as to blur the lines between solute and solvent.

See supplementary material for More Information on Dynamics of Alcohols and Examples of Temperature-Dependent 2D IR Spectra.

This work was funded by the Office of Basic Energy Sciences, U.S. Department of Energy (Grant No. DE-FOA-0000995). L.D.M. thanks NSERC for a fellowship.