of 10
Article
https://doi.org/10.1038/s41467-023-37667-7
Two-dimensional infrared-Raman spectro-
scopy as a probe of water
stetrahedrality
Tomislav Begu
š
i
ć
1
&GeoffreyA.Blake
1,2
Two-dimensional spectroscopic techni
ques combining terahertz (THz), infra-
red (IR), and visible pulses offer a weal
th of information about coupling among
vibrational modes in molecular liquids, thus providing a promising probe of
their local structure. Howe
ver, the capabilities of these spectroscopies are still
largely unexplored due to experimental limitations and inherently weak non-
linear 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) spe
ctrum. The structure-spectrum rela-
tionship can explain the te
mperature dependence of the spectral features
corresponding to the anharmonic coupling between low-frequency inter-
molecular and high-frequ
ency intramolecular vibrational modes of water. In
light of these results, we propose new e
xperiments and discuss the implica-
tions for the study of tetrahedrality of liquid water.
Understanding the dynamical structure of liquid water is paramount to
a number of chemical and biological processes. In particular, the tet-
rahedral ordering of molecules, stemming from the directionality of
hydrogen bonds, has been proposed as the origin of water
sanom-
alous behavior in the liquid phase
1
. However, the tetrahedral structure
of water has also been contested both computationally
2
and
experimentally
3
. To date, most of our microscopic, structural infor-
mation about liquid water comes from molecular dynamics (MD)
simulations, which depend strongly on the choice of electronic
structure theory or force
fi
eld parametrization. In contrast, experi-
mental tools capable of studying the local arrangement of water
molecules remain scarce
4
. For example, common techniques that
probe the structure of liquids, such as X-ray and neutron scattering,
typically report on highly time-averaged quantities and can be accu-
rately reproduced with very different structural motifs
2
,
5
7
.Vibrational
spectroscopy, such as IR absorption and Raman scattering, offers
complementary information about the strength of intermolecular
hydrogen bonds, which depends on instantaneous, local arrangement
of water molecules
8
. Different spectral regions have been probed to
investigate intermolecular hydrogen-bond bending and stretching
modes (up to 300 cm
1
), frustrated rotational (librational) modes
(400
1000 cm
1
), intramolecular bending (1650 cm
1
) and stretching
(3000
3800 cm
1
) modes, as well as combination bands between
intermolecular and intramolecular modes. Reported experiments and
simulations range from high-resolution spectroscopy on molecular
clusters
9
12
to studies of interfacial
13
15
or bulk water
16
25
.Recently,a
joint experimental and computational study
26
,
27
of the temperature
dependence of the Raman spectrum of liquid water revealed a clear
structure-spectrum relationship between the local tetrahedral order
parameter
28
, a measure of structuring of liquid water, and frequency
shifts and intensities of spectral peaks. In fact, the temperature
dependence of the spectrum was, in this way, fully explained as a
change in the thermal distribution of the tetrahedral order parameter.
Yet, conventional, steady-state spectroscopy of liquids typically pro-
duces broad, unresolved spectral features and misses dynamical
information.
Time-resolved and two-dimensional vibrational spectroscopies
have emerged in the past years as probes sensitive to local water
coordination. Even here, however, the long-standing tool of two-
dimensional infrared spectroscopy (2DIR)
4
,
29
,
30
, which measures the
interaction of light with intramolecular, high-frequency modes, pro-
vides only an indirect probe of intermolecular dynamics. 2DIR has
been successfully applied to analyze coupling among high-frequency
modes in ice
31
and in proteins
32
, but also to understand the vibrational
Received: 31 October 2022
Accepted: 22 March 2023
Check for updates
1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA.
2
Division of Geological and Planetary Sciences,
California Institute of Technology, Pasadena, CA 91125, USA.
e-mail:
tbegusic@caltech.edu
;
gab@caltech.edu
Nature Communications
| (2023) 14:1950
1
1234567890():,;
1234567890():,;
relaxation mechanisms in liquids
33
. For example, several 2DIR studies
on liquid water and ice speculated that mechanical anharmonic cou-
pling to intermolecular modes contributes to the relaxation of the OH
stretch
34
36
. In addition, non-Condon effects, which are related to
electrical anharmonic coupling between high- and low-frequency
modes, were shown to affect the 2DIR spectra of OH stretching mode
in liquid water
37
. To target the low-frequency modes directly, a number
of hybrid spectroscopic techniques have been proposed, involving
different sequences of THz, IR, and visible pulses, such as the THz-THz-
Raman
38
41
, THz-Raman-THz, Raman-THz-THz
42
45
, and THz-IR-Raman
(also called THz-IR-visible
46
,
47
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
48
. The THz-THz-
Raman and TIRV methods are related to other two-dimensional IR-
Raman techniques, namely the two-dimensional IR doubly vibration-
ally enhanced and IR-IR-visible sum-frequency generation
spectroscopies
49
54
. For example, TIRV experiments have revealed
unambiguous spectral signatures of coupling between the intramole-
cular O-H stretch and intermolecular hydrogen-bond bending and
stretching modes
46
. Theoretical simulations by Ito and Tanimura
55
predicted such spectral features and assigned them to both mechan-
ical 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
56
. Finally, Raman-THz-THz and THz-Raman-THz spectro-
scopieshavebeenusedtostudytheinhomogeneityofliquidwaterand
aqueous solutions
57
,
58
, as well as the coupling among intermolecular
and intramolecular modes in liquid and solid bromoform
59
,
60
.Evenso,
due to limited availability of ef
fi
cient THz emitter materials, not all
frequencies have been covered by the reported techniques. Speci
fi
-
cally, 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 experi-
ments. To date, only the TIRV frequency region has been experimen-
tally 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 temperatures, we analyze
whether nuclear quantum effects are discernible in the two-
dimensional IIR spectrum.
Results
Theoretical model
We simulated the IIR response function
40
,
55
R
ð
t
1
,
t
2
Þ
=

1
_
2
Tr
f½½
^
Π
ð
t
1
+
t
2
Þ
,
^
μ
ð
t
1
Þ
,
^
μ
ð
0
Þ
^
ρ
g
ð
1
Þ
using the equilibrium-nonequilibrium MD approach, in which the
quantum-mechanical trace is replaced by a classical average
61
64
R
MD
ð
t
1
,
t
2
Þ
=
β
ε
Π
ð
q
+,
t
2
Þ
Π
ð
q

,
t
2
Þ
_
μ
ð
q

t
1
Þi
,
ð
2
Þ
where
^
ρ
is the thermal density operator,
^
μ
=
μ
ð
^
q
Þ
is the dipole moment
operator, and
^
Π
=
Π
ð
^
q
Þ
is the polarizability.
q
t
denotes the position
vector of a classical trajectory initiated at (
q
0
,
q
0
), whereas
q
±,
t
corre-
sponds to the initial conditions
ð
q
±,0
=
q
0
,
p
±,0
=
p
0
±
εμ
0
ð
q
0
Þ
=
2
Þ
after
an instantaneous interaction with the electric
fi
eld.
ε
is a free para-
meter in the calculations and corresponds to the magnitude of the
external electric
fi
eld integrated over the short interaction time. For
suf
fi
ciently small
ε
, the thermal average of Eq. (
2
) is linear in
ε
, meaning
that the time-dependent response function
R
MD
(
t
1
,
t
2
) is, as expected,
independent of its exact value
63
,
65
. The two-dimensional spectra are
computed through a double sine transform
42
,
55
R
ð
ω
1
,
ω
2
Þ
=
Z
1
0
Z
1
0
R
ð
t
1
,
t
2
Þ
sin
ð
ω
1
t
1
Þ
sin
ð
ω
2
t
2
Þ
dt
1
dt
2
:
ð
3
Þ
To model water, we used a
fl
exible, point-charge qTIP4P/F force
fi
eld
66
, which has been well studied for spectroscopic simulations
25
and
benchmarked against a number of experimental thermodynamic
properties of water, including the radial distribution functions,
dielectric constant, density, and melting point. As a point-charge force
fi
eld, the model is computationally ef
fi
cient, which is needed for sta-
tistically averaging two-time response functions. In addition, because
it is
fl
exible, we could study all vibrational degrees of freedom,
including high-frequency intramolecular modes. To allow for non-
linear dependence of the dipole moments and polarizabilities on
nuclear coordinates, we employed the truncated dipole-induced-
dipole (DID) model. Following Hamm
44
, each water molecule was
amended with 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.
44
,
here the permanent polarizability was an explicit function of intra-
molecular degrees of freedom, according to ref.
67
.Theintroduced
induced dipole and polarizability effects were not used in the evalua-
tion of the forces, which were computed according to the original,
non-polarizable qTIP4P/F model. 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
,
44
,
68
. However, while some of them
were not readily available, others considered only rigid water mole-
cules or were
fi
tted to experiments using classical MD simulations,
which would prevent us from accurately exploring nuclear quantum
Fig. 1 | Steady-state spectra of water.
IR absorption (
a
) and anisotropic Raman (
b
)
spectra of liquid water at 300 K simulated with the classical MD and TRPMD,
compared with the experiments of refs.
16
,
27
,
91
. The inset in panel
b
compares the
low-frequency part of the simulated and experimental anisotropic Raman spectra.
Article
https://doi.org/10.1038/s41467-023-37667-7
Nature Communications
| (2023) 14:1950
2
effects. In Fig.
1
, we show that the employed force
fi
eld combined with
our induced dipole and polarizability models can reproduce the main
features of the experimental IR absorption and anisotropic Raman
spectra of liquid water. Additional details about the model and MD
simulations can be found in the Methods section.
Most of the results presented below rely on the validity of the
classical MD approach, which neglects the quantum-mechanical
properties of atomic nuclei. However, due to the presence of light
hydrogen atoms, their effect might not be negligible
66
. 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 equilibrium-nonequilibrium MD have
been available to study nuclear quantum effects on two-dimensional
IR-Raman spectra
69
71
. Recently, we have developed a new ring-
polymer MD (RPMD) approach
65
, which can simulate, at least
approximately, the nuclear quantum effects on the two-dimensional
IIR spectra. Brie
fl
y, the RPMD method
72
replaces the original quantum-
mechanical problem by an extended classical system consisting of
N
replicas (beads) of the original system connected by harmonic springs.
For a given potential energy surface, the extended classical system in
the limit of
N
reproduces the exact quantum-mechanical thermal
distribution of nuclear degrees of freedom, while if
N
=1, RPMD
reduces to classical MD. In our simulations, we used
N
= 32, which is
suf
fi
ciently 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 arti
fi
cial harmonic springs appear in the
spectra, we employed its thermostatted version (TRPMD)
73
,
74
.IR
absorption and anisotropic Raman spectra simulated with TRPMD are
presented in Fig.
1
.
Two-dimensional IIR spectrum of liquid water
Figure
2
shows the full two-dimensional IIR spectrum of liquid water
simulated at 300 K. Apart from the spectral features along the
ω
1
=
ω
2
diagonal, which appear due to mechanical or electrical anharmonicity
of individual vibrational modes, the spectrum contains off-diagonal
peaks that correspond to coupling among different vibrations. Overall,
the spectrum qualitatively agrees with the simulation of ref.
55
,with
the main difference in the relative intensities of different spectral
regions, which can depend strong
ly on the details of the potential
energy, dipole, and polarizability surfaces. In the following, we focus
on the two frequency regions indicated by the grey rectangles. One
region targets the coupling between intermolecular modes
(0 cm
1
<
ω
1
<1200cm
1
) and the intramolecular O-H stretch mode
(2900 cm
1
<
ω
2
<4100cm
1
), whereas the other covers all low-fre-
quency, intermolecular modes, and the intramolecular bend mode.
Model simulations of ref.
55
. (Supplementary Discussion 1 and Sup-
plementary Fig. 3) imply that the complex spectral lineshape of the
former, which we will refer to as the TIRV region
46
, are a product of
interplay between mechanical and electrical anharmonicity. Namely,
mechanical anharmonicity leads to a shape comprising a positive and a
negative lobe above and below the central
ω
2
frequency, whereas
electrical anharmonicity produces an approximately symmetric fea-
ture. To verify this interpretation, we would ideally construct
approximate dipole and polarizability models that depend linearly on
atomic coordinates, which would allow us to study the mechanical
anharmonic coupling pathways directly. This can be easily done for the
dipole moments by neglecting the induced part because the perma-
nent molecular dipole is a linear function of coordinates (see Meth-
ods). Unfortunately, the same cannot be achieved for polarizability
because both permanent and induced parts are nonlinear in atomic
coordinates. For this reason, we consider an alternative, Raman-IR-IR
(RII) pulse sequence, in which the polarizability is responsible for the
fi
rst interaction with the external electric
fi
eld. Assuming weak elec-
trical and mechanical anharmonic coupling, it can be shown that the
nonlinearity in the
fi
rst interaction does not contribute to the RII
spectrum (Supplementary Discussion 2). Therefore, the RII spectrum
with permanent (i.e., linear) dipole moments (shown in Fig.
3
a) con-
sists only of mechanical anharmonic coupling pathways. Because the
Raman response in the THz frequency range is weak, the corre-
sponding RII spectrum exhibits roughly an order of magnitude lower
intensity than the TIRV spectrum. For comparison with the full TIRV
spectrum, the simulated RII spectrum must be appropriately scaled
(Fig.
3
b) along the frequency axes (Supplementary Eq. (9)). The result
agrees with the proposed interpretation that the nodal shape of the
spectrum results from mechanical anharmonic coupling pathways.
In the low-frequency re
gion of interest, we observe a strong peak
at about (600 cm
1
,1650cm
1
) due to anharmonic coupling between
librations (400
1000 cm
1
) and the intramolecular bending mode
75
.
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 combina-
tion band is not captured in our model (Fig.
1
), even though the cor-
responding two-dimensional spectral feature clearly appears in the IIR
spectrum. This discrepancy can be explained by the fact that the
fundamentals of the one-dimensional spectra follow harmonic selec-
tion rules, while the peaks in the two-dimensional spectra appear solely
due to the anharmonic excitation pathways. Therefore, the combina-
tion bands in the one-dimensional spectra can be orders of magnitude
weaker than the fundamentals and still exhibit strong off-diagonal
peaks in the two-dimensional spectra. Indeed, we see that the diagonal
Fig. 2 | Two-dimensional IIR spectrum of liquid water simulated with classical
MD at 300 K.
Grey squares indicate the regions of the spectrum that we studied in
more detail, namely the region of TIRV spectroscopy and the low-frequency region
that covers intermolecular modes and the intramolecular bending mode.
Fig. 3 | Analysis of the mechanical anharmonic coupling. a
RII spectrum with
μ
ind
= 0 simulated using classical MD at 300K (multiplied by 10 so that the spectral
features can be seen on the same scale as in panel
b
).
b
The same spectrum mul-
tiplied by the frequency-dependent scaling factor derived in Supplementary Eq. (11),
which approximates the mechanical anharmonic coupling contribution to the TIRV
spectrum.
Article
https://doi.org/10.1038/s41467-023-37667-7
Nature Communications
| (2023) 14:1950
3
peak at around (1650 cm
1
,1650cm
1
) is much lower in intensity than
the off-diagonal libration-bending peak, implying that the anharmo-
nicity within the bending mode is weaker than its anharmonic coupling
to the intermolecular libration
s. Let us note that other force
fi
elds and
dipole/polarizability models also heavily underestimate or fail to
reproduce the combination band at 2150 cm
1
19
,
21
,
55
,
76
,unliketheab
initio approaches, which appear to systematically reproduce it
26
,
77
,
78
.
Since the force
fi
eldweusedreproducesthestructuralanddynamical
properties of liquid water rather accurately
66
, we tentatively assign the
absence of this combination band in the simulated spectra (Fig.
1
)to
the limitations of our dipole and polarizability models.
Temperature dependence of the IIR spectrum
We now turn to the temperature dependence of the two IIR
spectral regions (Fig.
4
). The TIRV part of the spectrum experi-
ences 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
,3600cm
1
). This indicates
that the electrical anharmonici
ty contribution to the spectrum,
marked by a strongly symmetric lineshape, increases compared to
that of the mechanical anharmonicity. Unlike the frequency shift,
this feature cannot be accessed t
hrough steady-state spectro-
scopy. Finally, we note that a peak around (400 cm
1
,3700cm
1
)
appears at elevated temperature
, while the other peak at about
(700 cm
1
,3400cm
1
) almost disappears. Interestingly, both of
these coupling terms have been di
scussed as possible origins of
the libration-stretch combination band appearing in the aniso-
tropic Raman spectrum of liquid water at 4100 cm
1
26
,
79
.Thesame
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
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 (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 anharmonicity of the hydrogen-bond
stretching modes. These seem to progressively disappear at elevated
temperatures. Importantly, this change in intensity could not be sim-
ply 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 arti
fi
cial broadening.
Longer simulations are possible with rigid water models, which have
been used to study explicitly the intermolecular modes and corre-
sponding two-dimensional THz-Raman signals in the time domain
44
,
62
.
The TIRV spectral region has been studied experimentally
46
,
47
.
However, the experiments could only measure an absolute value
Fourier transform spectrum, which is quite different from the real sine
transform discussed here and in ref.
55
.Speci
fi
cally, the experimental
spectrum reported in Fig. 8 of ref.
47
is related to the time-dependent
Fig. 4 | Temperature dependence of the IIR spectrum.
Two-dimensional IIR spectra of liquid water simulated with classical MD at different temperatures.
a
TIRV part of
the spectrum. The grey squares highlight the peak around (250 cm
1
,3600cm
1
).
b
Low-frequency part of the spectrum.
Fig. 5 | Temperature dependence of the TIRV spectrum related to the experiment of ref.
47
.
Two-dimensional TIRV spectra of liquid water at different temperatures, as
calculated using Eq. (
4
) and time-dependent response functions
R
(
t
1
,
t
2
) equal to those used in Fig.
4
.
Article
https://doi.org/10.1038/s41467-023-37667-7
Nature Communications
| (2023) 14:1950
4
response function by
~
S
ð
ω
1
,
ω
2
Þ
=
S
ð
ω
1
,
ω
2
Þ
+
S
ð
ω
2

ω
1
,
ω
2
Þ
E
THz
ð
ω
1
Þ
E
IR
ð
ω
2

ω
1
Þj
,
ð
4
Þ
S
ð
ω
1
,
ω
2
Þ
=
Z
1
0
Z
1
0
R
ð
t
1
,
t
2
Þ
e

i
ω
1
t
1
e

i
ω
2
t
2
dt
1
dt
2
,
ð
5
Þ
where
E
THz
(
ω
)and
E
IR
(
ω
) are the THz and IR electric
fi
elds obtained as
square roots of the pulse intensity spectra presented in Fig.
2
b, c of
ref.
46
. The main differences between our simulation and the experi-
ment of ref.
47
arise due to the limitations of our model, namely the
narrow and blue-shifted OH stretch (Fig.
1
). Furthermore, the
experimental spectrum exhibits additional spectral features that
cannot be explained with Eq. (
4
) and the pulse shapes we used
because they fall outside of the bandwidth of the instrument response
function. Unlike the sine-transform spectra shown in Fig.
4
, Fourier-
transform spectra convolved with external electric
fi
elds (Fig.
5
)only
become less intense with increasing temperature but show no
interesting spectral change. Therefore, an accurate determination of
the full response function will be needed to experimentally measure
the spectral features appearing in the sine-transform TIRV spectra.
Fortunately, a scheme that could achieve this goal has been recently
applied to the TIRV spectrum of liquid dimethyl sulfoxide
80
,
demonstrating that similar experiments on water are within reach.
Spectral signatures of water
s tetrahedrality in the IIR spectrum
To understand the temperature-dependent spectral features, we fol-
low 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 quanti
fi
ed by the local tetrahedral order para-
meter
Q
81
83
, which measures how the arrangement of the four neigh-
boring 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 (Fig.
6
a). In the studied tem-
perature range, the bimodal distribution can be decomposed into two
temperature-independent components whose populations change
with temperature
26
.However,asshownbyGeissler
84
, the appearance
of an isosbestic point does not necessarily imply that water is a het-
erogeneous mixture. More recently, the increased population of low-
Q
water molecules at higher temperatures has been assigned to the
appearance of neighboring molecules in the interstitial position
between the
fi
rst and second solvation shells
83
.
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
(Supplementary Discussion 3). Then, we simulated the spectra
(Fig.
6
b
e) originating predominantly from the molecules with either
high (
Q
> 0.72) or low (
Q
< 0.62) local tetrahedral order parameters.
The two IIR spectra align remarkably well with the observed
temperature-dependent changes. Speci
fi
cally, the high-order TIRV
spectrum (Fig.
6
c) 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-stretching peak at
(400 cm
1
,3700cm
1
) appears almost exclusively in molecules with
low tetrahedral order. Similarly, the libration-bending peak in the low-
frequency part the spectrum (Fig.
6
d, e) appears at lower
ω
1
fre-
quencies for low
Q
. Overall, the temperature dependence of the IIR
spectra can be almost exclusively assigned to the changes in the dis-
tribution 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
cov-
ering up to
800 cm
1
, could provide an alternative probe of the local
molecular structure in liquid water and aqueous solutions. Impor-
tantly, 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.
This result is consistent with other experimental observations and
theoretical models. Namely, it is known that as the tetrahedral order
parameter increases, the hydrogen-bond stretching mode frequency
increases and the OH stretching frequency decreases
26
, which can be
interpreted in terms of a growing mechanical anharmonic coupling
strength between the two modes. A more detailed analysis is possible
based on the work of Auer and Skinner
85
, who studied the dependence
of the intramolecular stretching frequency and the corresponding
dipole moment derivative on the electric
fi
eld
E
generated by the
surrounding water molecules at the hydrogen atom of the central
molecule and projected along the OH bond. They found that the fre-
quency can be
fi
t to a quadratic function of this electric
fi
eld, while the
dipole moment derivative is approximately linear in
E
:
ω
OH
=
ω
ð
0
Þ
OH

ω
ð
1
Þ
OH
E

ω
ð
2
Þ
OH
E
2
,
ð
6
Þ
μ
0
OH
=
μ
0
Þ
OH
+
μ
1
Þ
OH
E
,
ð
7
Þ
where all parameters
ω
ð
α
Þ
OH
and
μ
α
Þ
OH
are positive and de
fi
ned in Table I.
of ref.
85
. More recent models
86
,
87
included a very weak quadratic term
for the dipole moment as well, which can be neglected within the
typical range of values of the electric
fi
eld. For us, the relevant
anharmonic coupling quantities are the derivatives of
ω
OH
and
μ
0
OH
Fig. 6 | Effect of tetrahedral order on the IIR spectra. a
Distribution of the local
tetrahedral order parameter of liquid water as a function of temperature. The
isosbestic value is indicated by the solid, grey line. Two dashed, grey lines specify
the region excluded from subsequent two-dimensional IIR spectra simulations
(0.62 <
Q
< 0.72). Two-dimensional IIR spectra of water molecules with high
(
Q
>0.72,
c
,
e
)orlow(
Q
<0.62,
b
,
d
) tetrahedral order parameter, simulated with
the classical equilibrium-nonequilibrium MD approach at 320 K.
Article
https://doi.org/10.1038/s41467-023-37667-7
Nature Communications
| (2023) 14:1950
5
with respect to an intermolecular hydrogen-bond mode
q
HB
,
ω
OH
=
q
HB
=
ω
ð
1
Þ
OH
+2
ω
ð
2
Þ
OH
E
Þ
E
=
q
HB
and
μ
0
OH
=
q
HB
=
μ
1
Þ
OH
E
=
q
HB
.
From these equations, we see that the mechanical anharmonic
coupling, determined by
ω
OH
/
q
HB
, grows with the electric
fi
eld,
whereas the electrical anharmonic coupling, proportional to
μ
0
OH
=
q
HB
, has no explicit dependence on
E
.Inref.
83
, it has been
shown that
E
correlates positively with the tetrahedral order parameter
Q
, which further validates the interpretation of our two-dimensional
spectroscopy simulations, i.e., stronger mechanical anharmonic
coupling in high-
Q
molecules.
Nuclear quantum effects on the two-dimensional IIR spectrum
Finally, we report the
fi
rst TRPMD simulation of the two-dimensional
IIR spectrum of liquid water and compare it with the MD simulation in
Fig.
7
. Some nuclear quantum effects can be easily explained in terms
of the differences between one-dimensional spectra (Supplementary
Fig. 7). Namely, the red-shifted O-H stretch band in the TRPMD steady-
state spectra re
fl
ects in a red shift along
ω
2
of the spectral features in
the TIRV spectrum (Fig.
7
a, b). 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 (Fig.
7
c,
d), 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 fea-
ture at (
ω
1
750 cm
1
,
ω
2
50 cm
1
)inTRPMDorthechangeinthe
spectral lineshapes in the TIRV region at
ω
1
>400cm
1
. Nevertheless,
these differences are not suf
fi
ciently large to alter the above conclu-
sions 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 tem-
perature dependence of the spectrum, we have identi
fi
ed features
that report on the degree of local tetrahedral ordering in the
fi
rst
coordination shell. Further computational work is needed to con-
fi
rm whether IIR spectra of related systems, such as aqueous solu-
tions, alcohols, or ice, 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 con-
tributions of mechanical and electrical anharmonic coupling change as
a function of local tetrahedrality, which cannot be studied with con-
ventional, one-dimensional spectroscopy. Speci
fi
cally, the electrical
anharmonicity dominates at higher temperatures, where the tetra-
hedral order parameter is low, whereas the mechanical coupling, due
to the anharmonic terms in the potential energy surface, is pro-
nounced at lower temperature and higher tetrahedrality. This
fi
nding
was related to the electric
fi
eld caused by the surrounding molecules
and its positive correlation to the tetrahedral order parameter. Overall,
the observed change in the mechanical anharmonic coupling could
impact our understanding of the OH stretch relaxation mechanism.
Although this relaxation has been largely assigned to the coupling with
the intramolecular bending overtone, the intermolecular, hydrogen-
bond modes are also believed to play an important role
34
36
,
88
.Strong
mechanical anharmonic coupling between intramolecular stretching
and intermolecular hydrogen-bond stretching modes could shorten
the excited-state lifetime of the OH stretch. Such correlation would
also agree with the experimentally observed trend of shorter OH
stretch excitation lifetime at a lower temperature.
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 physi-
cal chemists for decades.
Methods
Dipole moment and polarizability models
The induced dipole moments and polarizabilities were modeled as
44
μ
ind
=
X
i
,
j
i
j
Π
i

E
ij
,
ð
8
Þ
Π
ind
=

X
i
,
j
i
j
Π
i

T
ij

Π
j
,
ð
9
Þ
where
E
ij
=
P
a
2f
M,H1,H2
g
Q
ja
r
ja
,
i
O
=
r
3
ja
,
i
O
is the electric
fi
eld produced by
molecule
j
on the oxygen atom of molecule
i
,
T
ij
=
T
(
r
j
O,
i
O
),
T
(
r
)=
(
r
2
3
r
r
)/
r
5
is the 3 × 3 dipole-dipole interaction tensor,
r
ja
,
i
O
=
r
ja
r
i
O
,
r
ja
denotes the three-dimensional position vector of
atom
a
{M, 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. M denotes the M site of th
e TIP4P model, whose position is
related to the molecular geometry of molecule
i
by
r
i
M
=
γ
r
i
O
+(1
γ
)
(
r
i
H1
+
r
i
H2
)/2, where
γ
=0.73612
66
. Permanent dipole moment of
molecule
i
was
μ
i
=
a
{M,H1,H2}
Q
ia
r
ia
/
c
,where
c
=1.3isascalingfactor
that reduces the dipole moment of water to its gas-phase value
44
,
63
.The
permanent polarizability of each molecule
i
,
Π
i
Π
i
(
r
i
O
,
r
i
H1
,
r
i
H2
), was
computed according to ref.
67
inthemolecularreferenceframeand
rotated into the laboratory frame. Then, the total dipole moment and
polarizability are given by
μ
=
X
i
μ
i
+
μ
ind
,
ð
10
Þ
Π
=
X
i
Π
i
+
Π
ind
:
ð
11
Þ
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 in
fi
nity
Fig. 7 | Nuclear quantum effects on the IIR spectrum.
Two-dimensional IIR
spectra of liquid water at 280 K, simulated with equilibrium-nonequilibrium MD
(
a
,
c
) and TRPMD (
b
,
d
).
Article
https://doi.org/10.1038/s41467-023-37667-7
Nature Communications
| (2023) 14:1950
6