The hydrogen-bond network of water supports propagating optical phonon-like modes

The local structure of liquid water as a function of temperature is a source of intense research. This structure is intimately linked to the dynamics of water molecules, which can be measured using Raman and infrared spectroscopies. The assignment of spectral peaks depends on whether they are collective modes or single-molecule motions. Vibrational modes in liquids are usually considered to be associated to the motions of single molecules or small clusters. Using molecular dynamics simulations, here we find dispersive optical phonon-like modes in the librational and OH-stretching bands. We argue that on subpicosecond time scales these modes propagate through water's hydrogen-bond network over distances of up to 2 nm. In the long wavelength limit these optical modes exhibit longitudinal–transverse splitting, indicating the presence of coherent long-range dipole–dipole interactions, as in ice. Our results indicate the dynamics of liquid water have more similarities to ice than previously thought.

) Supplementary Figure 4: Longitudinal (left) and transverse (right) dielectric susceptibility for a system of 1,000 MeOH molecules. The longitudinal librational peak at ≈ 700 cm −1 clearly disperses with k, while the transverse peak at ≈ 600 cm −1 disperses slightly with k. The higher frequency peaks exhibit no dispersion. The static dielectric function ε(k, 0) has not converged properly in the transverse case, so the magnitude of the peaks is not converged. The broad band which peaks at 100 cm −1 exhibits dispersion. We hypothesize this dispersion is due entirely to the translational modes, however we cannot say for sure since the librational and translational modes overlap in this region. The peak at ≈ 500 cm −1 is due to CCN bending. The static dielectric function ε(k, 0) has not converged properly, so the magnitude of the transverse peaks is not converged correctly, but the position of the peaks and dispersion can be seen. appear in the residual -the lower frequency peak is dispersive, having the same dispersion relation as the fitted peak, suggesting that it is actually part of the dispersive peak lineshape that is not captured by our lineshape function. The higher frequency peak in the residual is non-dispersive and is in the same location for both the transverse and longitudinal susceptibility.
type ω L1 ω L2 ω L3 ref. There are several different ways to decompose a spectra into contributions from molecules separated by distance R: Kirkwood dipole-sphere method This is the method we choose, which is a modification of the sphere-sphere method (see below). We start with the time-correlation function of interest : Here µ can be replaced with any dynamical variable of interest, for instance p T (k, t) or j(t). We omit the k dependence for simplicity. The most straightforward way is to limit the molecules around each molecule i to those in a sphere of radius R: This is similar to the method employed by Bopp & Kornyshev. During the the course of a simulation molecules enter and leave each sphere, which creates noise, requiring longer averaging times. This can be improved by utilizing a smooth cutoff function: where Here D is a sharpness parameter determining the relative sharpness of the cutoff. We choose not to use smoothing however, finding it to be unnecessary. The result is a spectra χ(k, ω, R) showing contributions from molecules up to radius R. The resulting function exhibits the expected R → 0 limit, yielding only the self contribution. In the R → ∞ limit, the original full response function is recovered. This function can then be numerically differentiated to show the contributions from shells of thickness ∆R centered at distance R.
Sphere-sphere method Another method discussed by Heyden, et. al. (2010) is to study the autocorrelation of the total dipole moment of a sphere of radius R centered around a reference molecule, and then average over each molecule in the system. [11] Heyden, et. al. recommend the normalization factor N i (t) = (1 + j∈Ri P 2 ij ) −1/2 to normalize for number of molecules in each sphere. This normalization factor is chosen so that in the bulk limit (R → ∞) the original full response function is obtained (in that limit N i = 1/ √ N mol ). In the limit R → 0 only the self-term contributes. Results from this method must be interpreted with a bit of care since the calculation includes all cross-correlations between molecules within the sphere centered around the reference molecule. We found that this method is more sensitive to intermolecular correlations, in particular the H-bond stretching at ≈ 250 cm −1 (not shown). Altogether though we found the results from this method are complementary with our results from the dipole-sphere method.
Spatial grid method To achieve higher resolution, Heyden, et. al. also introduce a spatial grid method. [11] The method works by bining the molecular dipoles into grid cells. To reduce noise caused by moleules moving in and out of bins the binning is Gaussian, meaning the dipoles are smeared with a Gaussian function. Unlike the other methods the spatial grid method does not yield the self part as R → 0 so this limit requires special interpretation.
Bopp & Kornyshev show that to get accurate results in k space it is important to use the polarization vectors for each molecule rather than just the dipole moment. To calculate the polarization vector we use the method of Raineri & Friedman. [12] We utilize the defining relation for the polarization: When transformed into Fourier space this becomes: We introduce polarization vectors for each molecule p i (k) so that we have where The molecules are indexed by i and the atoms on each molecule are indexed by α. r iα = r i (t) − r α (t) is the distance from each atomic site to the center of mass of molecule i. Following Raineri & Friedman, we use the identity and taking into account the charge neutrality of each molecule we obtain The transverse part is then calculated as P T =k × P , while the longitudinal component is P L =k · P . The longitudinal component can also be calculated more directly from: This yields the following Kubo formula for the longitudinal part of the response: For a system composed of point charges, the charge density is : Again, the index i runs over the molecules while α runs over the atomic sites on each molecule. The charge density in k-space becomes: Note that this can be Taylor expanded as: (−ik · r iα (t)) n n! = M (k, t) + Q(k, t) + O(k, t) + · · · Here M (k, t), Q(k, t), O(k, t) are contributions due to the molecular dipoles, quadrupoles and octupoles. In the limit k → 0 it from supplementary equation 15 it can be seen that only the dipole term contributes to the susceptibility. In the k → 0 limit one obtains This type of expression has been used previously as an approximate expression at small k. [13] However, Bopp & Kornyshev show quite convincingly that for water the higher order multipole terms are very important, even at the smallest k available in computer simulation. [14] Neglect of the higher order terms leads to severe consequences at large k, and one will not recover the physical limit lim k→∞ χ L/T (k) = 1 unless higher order terms are included.