Combinational Vibration Modes in H2O/HDO/D2O Mixtures Detected Thanks to the Superior Sensitivity of Femtosecond Stimulated Raman Scattering

Overtones and combinational modes frequently play essential roles in ultrafast vibrational energy relaxation in liquid water. However, these modes are very weak and often overlap with fundamental modes, particularly in isotopologues mixtures. We measured VV and HV Raman spectra of H2O and D2O mixtures with femtosecond stimulated Raman scattering (FSRS) and compared the results with calculated spectra. Specifically, we observed the mode at around 1850 cm–1 and assigned it to H–O–D bend + rocking libration. Second, we found that the H–O–D bend overtone band and the OD stretch + rocking libration combination band contribute to the band located between 2850 and 3050 cm–1. Furthermore, we assigned the broad band located between 4000 and 4200 cm–1 to be composed of combinational modes of high-frequency OH stretching modes with predominantly twisting and rocking librations. These results should help in a proper interpretation of Raman spectra of aqueous systems as well as in the identification of vibrational relaxation pathways in isotopically diluted water.

Polarized (VV) and depolarized (HV) Raman spectra were calculated by Fourier transform of the corresponding response function C(t) [1] where C(t) is a time-correlation function of a quantummechanical operatorÔ: C(t) = ⟨Ô(0)Ô(t)⟩, wherê O(t) = e iĤt/ℏÔ e −iĤt/ℏ is the Heisenberg representation of the operatorÔ,Ĥ is the Hamiltonian of the system, * mpastorczak@ichf.edu.pl † akanane@udel.edu and ⟨. . .⟩ denotes a quantum-mechanical equilibrium statistical average. The particular form of the operatorÔ depends on the type of line shape of interest [1]. The following time-correlation functions correspond to VV and HV Raman spectra where C iso (t) and C aniso (t) are isotropic and anisotropic time-correlation functions, correspondingly. [1] They are given by the time-correlation functions of the polarizability tensor α. The isotropic Raman response function is the auto-correlation function of the spherical part of the polarizability tensor A(t) [1] where A(t) = 1 3 Tr[α(t)], and Tr[. . .] denotes a trace over three-dimensional space. The anisotropic Raman S2 response function is given by the autocorrelation function of the anisotropic polarizability B(t) [1] where B(t) = α(t) − IA(t), and I is the unit tensor. Full quantum-mechanical calculation of Raman line shapes for realistic condensed-matter systems is computationally infeasible due to the prohibitively large size of the Hilbert space required for the description of such systems. Mixed quantum-classical methods often are the only practical ways to calculate Raman spectra for large realistic systems. The basic idea of a mixed quantumclassical approach used in this work is to divide all degrees of freedom in the system into the subsystem comprising vibrations of interest and the so-called environment or "bath", which is composed of slow degrees of freedom [2,3]. Within this approach vibrations of interest, which are those being probed spectroscopically and typically have frequencies ℏω >> k B T , are treated quantum mechanically (here k B is the Boltzmann constant). Low-frequency degrees of freedom associated with translations, rotations, and low-frequency vibrations constitute the "bath" and are treated classically. These degrees of freedom ideally (but not necessarily in practice) have frequencies less than k B T . Finally, there are other degrees of freedom, typically high-frequency vibrations that are thought not to contribute to the spectroscopy of interest, and therefore are ignored. All bath degrees of freedom are treated classically via classical molecular dynamics. High-frequency OH (OD) stretch and HOH (DOD, HDO) bend vibrations are treated quantum-mechanically as explained below.
We used the one-exciton model to describe molecular vibrations (chromophores). In the local-mode basis the vibrational exciton Hamiltonian is given by This Hamiltonian describes a system with the ground state |0⟩ where all chromophores are in their ground states with H 0 being the corresponding Hamiltonian.
The system is also composed of a manifold of N singly excited states, {|n⟩ : n = 1, . . . , N }, corresponding to a single excitation on the nth chromophore with all other chromophores in their ground states: For n = m the Hamiltonian H nm describes the local excitation on chromophore n; otherwise it is the coupling between n and m chromphores. Within the mixed quantum-classical theory based on the one-exciton model above, for a system of coupled chromophores Raman line shape can be written as [4] whereC(t) is the corresponding (e.g., VV or HV) mixed quantum-classical time-correlation functions [1] where a jpq is the p, q ∈ {x, y, z} component of the transition polarizability tensor for chromophore j,p is the unit vector in the direction of polarization of the incident light beam,q is the unit vector in the direction of polarization of the scattered beam, T 1 is the average lifetime of the vibrationally excited states which is taken to be 260 ps [5]. F jk (t) are the matrix elements of matrix F(t), which satisfies the following equation of motion where the matrix elements of K(t) are given by where the site frequencies are denoted by ω n , and {ω nm : n ̸ = m} are the vibrational chromophore-chromophore couplings. If n and m belong to two chromophores on the same molecule, then ω nm is the intramolecular coupling element and will be denoted by ω intra nm , while if n and m belong to different molecules, ω nm is an intermolecular coupling element and will be denoted as ω inter nm . Site frequencies as well as inter-and intramolecular couplings generally change with time due to fluctuations of the bath variables, giving rise to time-dependence of the matrix K(t).
Numerical solution of Eq. (S10), subject to the initial condition F nm (0) = δ nm , produces a set of time-evolved matrices F(t), which are used to calculate Raman line shapes according to Eqs. (S8)-(S9).
In addition to the time-propagation matrix F(t), knowledge of the components of the time-evolved transition polarizability tensor α for each chromophore is required in order to calculate Raman line shapes.
Intramolecular and intermolecular chromophorechromophore couplings are also required to construct S3 matrix K. The intramolecular coupling between two OH (OD) stretch chromophores n and m is given by [8,9] where p n is the 1-0 momentum matrix elements of OH (OD) stretching chromophore n, ϕ is the angle between OH (OD) stretch chromophores n and m, which is fixed at 104.52 • for the rigid water model used in this work, m O is the mass of oxygen atom and k intra nm is the second derivative of the potential energy with respect to O-H(D) bond lengths, calculated at r eq OH . Due to frequency mismatch intramolecular couplings in HOD molecules were set to 0.
The intermolecular coupling between two OH (OD) stretch chromophores n and m is described using transition dipole coupling (TDC) scheme where k inter nm is given by (S15) wheren nm is a unit vector along the line connecting the locations of OH (OD) stretch transition dipoles n and m, and r nm is the distance between these locations. The transition dipole of an OH (OD) stretch chromophore is located along the OH (OD) bond 0.67Å away from the oxygen atom [9].
The nature of the intermolecular coupling between bend fundamental vibrations is less clear [10,11]. However, because in this work we are only concerned with the bend overtone vibration and, since its 2-0 transition matrix element is smaller than the 1-0 transition matrix element, the intermolecular coupling between bend overtone modes is neglected.
The intramolecular coupling between the fundamental OH stretching and overtone of the HOH bending vibrations is known as Fermi resonance coupling. As in the previous work it is assumed constant [12,13]. Its value 25 cm −1 (in the local-mode basis) is determined by fitting the calculated polarized Raman spectra to its experimental counterpart [12]. This value is in a good agreement with earlier works based on somewhat similar fitting strategies [14][15][16][17][18]. Fermi resonance in D 2 O and HDO was ignored because, in contrast to H 2 O, the value of the Fermi resonance coupling in such systems has not been determined.
It remains to describe how transition frequencies ω n , vibrational couplings ω nm , dipole moment derivatives µ ′ n , coordinate x n , and momentum p n matrix elements change with time. As done previously [3,19,21,22], we correlate vibrational frequencies with the effective electric field computed for each configuration of water molecules at each time step from the point charges utilized in the molecular dynamics force field. The OH (OD) stretch frequency is determined from the electric field E at the location of the hydrogen (deuterium) atom along the direction of the OH (OD) bond. The mapping relationships between E and OH (OD) stretch transition frequencies ω n are determined so as to reproduce the correct distribution of vibrational OH (OD) stretch frequencies determined by the density functional theory calculations. For example, in order to calculate transition frequencies and dipole derivatives of the OH-stretch, a MD trajectory with E3B2 force field was generated first. Then clusters of water molecules centered on a randomly chosen water hydrogen atom were selected from the MD trajectory and electronic structure calculations are performed to determine transition frequency and transition dipole moment of the OH stretch of the central molecule [19]. All water molecules whose oxygen atom was within a cutoff of 4.0Å of the central hydrogen atom were explicitly included in the calculation, while all waters between 4.0 and 7.831Å were implicitly in-cluded only as point charges according to E3B2 water model (+0.52e for hydrogen and -1.04e for the dummy atom M site [23]). The OH bond of the central water molecule was systematically stretched and contracted to map out a one-dimensional potential energy surface using the Gaussian software package [24]. The selected OH bond is varied between 0.6572Å and 1.7572Å with a grid spacing of 0.02Å.The coordinates of all other atoms in the cluster are kept fixed as the bond length is changed. At each bond length, both the total energy and total dipole moment of the system are calculated at the B3LYP/6-311++G(d,p) level of theory and basis set [25][26][27] using the gaussian09 quantum chemistry program. It should be noted that the same procedure performed for the HOD monomer leads to OH-stretch frequency of 3699 cm −1 [19] which can be compared to the experimental value of 3707 cm −1 [28]. To correct for this small difference, we multiply all calculated transition frequencies by the scale factor of 1.0021. Similar approach was also used to determine OD stretching frequency.
The OH and OD dipole derivative µ ′ map is determined similarly as a function of E by a quadratic leastsquares fit to the data from the B3LYP/6-311++G** calculations. Coordinate x and momentum p matrix elements were found to depend linearly on the transition frequency. The intramolecular coupling ω intra nm was fit to a linear function of the sum of the electric fields for chromophores n and m.
Both the HOH bend fundamental transition frequency and first overtone frequency are correlated with the sum of the components of electric fields perpendicular to the two OH bonds and pointing towards the HOH bisector, evaluated at the hydrogen atoms [10]. HOH bending maps for 1-0 and 2-1 transitions were determined previously [10]. Since in this work the 2-0 transition (first overtone) frequency is needed, it was obtained as the sum of 1-0 and 2-1 transition frequencies from Ref. 10. Note that the 2-1 frequency map includes an anharmonicity of ∼30 cm −1 . Because DOD and HOD bending maps have not been developed, in this work, we set their overtone frequencies to 2380 and 2950 cm −1 correspondingly as estimated from experimental spectra [29,30].
A summary of all spectroscopic maps and other parameters used in this work is illustrated in Table S1.

S3. BOOK HEXAMER GEOMETRY
Cartesian coordinates of all atoms of the water hexamer used in VPT2 calculations are shown in Tab. S2.