Energy transfer within the hydrogen bonding network of water following resonant terahertz excitation

Energy flow in the hydrogen bonding network of water is traced by resonant terahertz excitation and off-resonant optical probing.

Classical Molecular dynamics simulations. In the FFMD, the atomic pair interactions are parameterized using a combination of Lennard-Jones and Coulomb potentials. We use a rigid three-point water model (SPC/E) optimized to reproduce the density and structure of water (50). We simulate a box with 5360 SPC/E water molecules in the NVE ensemble using the GROMACS version 2018 (51). We use 3D periodic boundary conditions, a Lennard-Jones potential which is shifted by a constant such that it is zero at the cut-off at 0.9 nm. 3D Particle Mesh Ewald summation are applied for the electrostatics beyond 0.9 nm. At time 0 = 5 ps we apply a THz pulse according to ( ) = 1 2 0 −( − 0 ) 2 / 2 cos (2 ( − 0 ) + ) (S1) with = 0.5/ps, = 0.5 ps = 0.5 ps, = 0.628319 and 0 = 1 V/nm, which corresponds to a THz electric filed with a maximum amplitude of ~2 MV/cm (similar to the experimental one).
We calculate the total KE from the atomic velocities and the translational KE trans from the velocity of the molecular centers of mass. The KE rot is the difference between the two. During the simulation, the temperature increases on average less than 1 K. The curves shown in the paper are based on the average of 5000 trajectories of 50 ps, each containing a single pulse. We use a time step of 2 fs and a write-out frequency of 20 frames per ps.
We check for finite-size effects by simulating a water box of 128 molecules using the same procedure (see Fig. S8). As expected, the statistics are insufficient to draw any conclusions when only a small number of pulses is averaged, but for 3200 pulses and above, the curves reproduce the results for the larger box size. This shows that finite-size effects are unimportant.
Kinetic energy decomposition. In the AIMD simulations, at each MD snapshot the decomposition was done for every molecule in the instantaneous molecular internal coordinates defined by the Eckart conditions (the Eckart frame), thus minimizing the vibrational-rotational (Coriolis) cross term (52). In the FFMD, the KE rot is calculated as the difference between the total KE and the KE of the molecular centers of mass.
Calculation of polarizability anisotropy. In order to calculate the relaxation of polarizability anisotropy shown in Figure 6c in the main text, we have based our calculations on the dipole code developed by R. Khatib et al. (53) (code downloaded from github: ttps://github.com/RemiKhatib/dippol). We chose to use a single-point model for water's molecular dipole moment and polarizability, i.e. a point dipole and polarizability are placed at the center of mass of each water molecule, and oriented according to the molecular frame of reference. The polarizability of each molecule at each at each sampled MD time step is selfconsistently calculated using the first-order dipole-induced dipole moment model. In this model, for each molecule in the system one starts with it's gas phase dipole 0 and polarizability 0 , and then the induced polarizability (due to the electric field of the surrounding molecules) is calculated from: where is the dipolar interaction tensor: and a cutoff of 7.5 A˚ has been used to compute the induced polarizabilities. Obviously, the solution of equation (S2) requires knowledge of the induced polarizability of every other molecule in the system within the cutoff range, hence the equation is solved iteratively until self-consistency is achieved.
For 0 we have used values obtained from accurate ab initio gas phase calculations at the coupled cluster level of theory (54). The polarizability anisotropy was calculated by first summing the polarizabilities of all the molecules in the laboratory frame at some time , and then calculating Δ (t) = xx (t) − ( yy (t) + zz (t))/2. The final result is obtained as an average over 5000 FFMD trajectories, each with 5360 water molecules.
As a validation of this model, we have compared the polarizability anisotropy values obtained using it, with those obtained from fully periodic ab initio Berry-phase calculations on snapshots extracted from the AIMD trajectories and averaged over all of them (55), and have found the results to be in agreement within the error bars of the AIMD.

Calculating the power exerted by the field on molecular translations and rotations.
In this section, we follow A. Pasquarello

Translations
The power exerted an electric field on a charge moving with velocity v is simply given by the dot product of the electric force and the velocity F ⋅ v. In the case of a point charge, the electric force is simply given by the electric field vector multiplied by the charge. For an isolated charge distribution with an overall zero charge, like a neutral molecule in vacuum, the electric field thus does not exert any work on the molecule when it translates, and the lowest order interaction with the field would be with torque exerted on the dipole moment of the charge distribution, if the latter is not vanishing by symmetry.
For a molecule in condensed phase, even for a molecule that is formally neutral, things are different. Here, as a molecule translates, it interacts with its neighbors leading to distortions in the electronic distribution and polarization, which can lead to a non-vanishing coupling of the molecular translation with the applied field, as if the molecule is charged. Thus in this case, one can define a so-called dynamical molecular monopole in terms of the coupling between translational molecular modes and electric fields.
This dynamical monopole or Born effective charge, (Z) for some atom labelled i (Z i ) in a periodic system is defined as: where Ω is the volume of the repeating periodical cell, and u i is some Cartesian displacement of the atom along the -axis. Z is thus given by a second-rank Cartesian tensor. An equivalent definition is: where F is the force induced on atomic nucleus i due to the electric field E, and Einstein summation notation is used. These dynamical charges are physical observables and measure the coupling between atomic displacements and electric fields.
In this work, we have used the dynamical monopoles of liquid water as reported by Pasquarello and Resta in the work cited in the previous paragraph. These were based on condensed phase calculations using DFT with the BLYP function, in close agreement with our AIMD setup. The total force exerted on a molecule is simply given by equation (S5) summed over all the atoms in the molecule. With the force at hand, we calculate the power exerted by the field on molecular translations, and integration of the latter quantity from time 0 to t gives the total work done on translations.

Rotations
The power exerted by the field on rotations of each molecules is given by the dot product of the electric torque with the molecular angular momentum. We calculate the molecular angular momentum in the Eckart frame, taking the molecular center of mass as the origin, and the electric torque ( ) is given by: where is the Levi-Civita tensor, and r s are the position vectors of the atomic nuclei in the molecule. Again, integration of the power yields the work done on rotations.

Figure S7 | Survival probability of H-bonds.
Blue: The probability for a hydrogen bond that was present at an initial time t=0 (which we choose to coincide with the first peak in the pulse) to survive at later times. For those hydrogen bonds that are broken we have also calculated the probability that the hydrogen bond is broken due to translational (red, see illustration at bottom left) or rotational (green, see illustration at bottom right) motion of the hydrogen bond acceptor relative to the hydrogen bond donor. In all cases the dotted lines show the same quantities for an equilibrium microcanonical fieldfree trajectory, which shows that, while the H-bond lifetime remains unaffected by the pulse, there is an increase in the contribution of the translational motion to the H-bond breaking events which is compensated for by a decrease in the rotational contribution. In this plot we define a hydrogen bond to have an O-O distance < 3.5 Å, and an angle of < 30°. A hydrogen bond is considered broken once these criteria are violated for a continuous time segment longer than 20 fs. Figure S8 | Impact of the finite size of the MD simulation box. Rotational and translational kinetic energy calculated in FFMD simulations using 128 water molecules, averaged over different numbers of pulses. The result for 10000 pulses agrees with the results obtained using 5360 water molecules (Fig. 5 of the main text), showing that finite-size effects are negligible. Here, the applied THz field has the amplitude of 0 = 2.5 V/nm (see Eq. S1) to gain higher signal to noise ratio.

Bipolarity of the TKE signal.
As mentioned in the main text, the prerequisite for observing a bipolar TKE signal of an isolated molecule is its negative polarizability anisotropy, i.e. the polarizability along the dipole axis is smaller than the polarizability average in the plane perpendicular to this axis. As this criterion isn't fulfilled for single water molecules, we go beyond the single molecule intrinsic polarizability and envisage low-order clusters, which is in the simplest case a dimer, to include the impact of the neighboring molecules to the total polarizability of water.
As highlighted in Fig. S10, the partial O-H group of molecule can be regarded as the permanent dipole of the dimer. By applying a THz torque to the OH dipole, the energy of this coupling is transferred to its neighboring molecule and increases its translational kinetic energy. Consequently, the latter excess kinetic energy leads to the stronger collision of molecule with its adjacent neighbors, resulting in the deformation of its electron cloud. The latter electron cloud deformation is indeed the additional collision-induced polarizability contribution to the total polarizability of the dimer. Indeed, the translational motion of molecule can be decomposed to its components along and perpendicular to the OH dipole. The latter of which corresponds to the motion of molecule in the normal mode coordinate of H-bond bending mode, and the former matches the coordinate of the stretch vibration.
We suggest that molecule gains more translational kinetic energy for its motion in the Hbonding mode coordinate because of the spectral vicinity of this mode and the excited region. Moreover, the H-bond bending has a pronounced amplitude in the Raman response of water at 1-2 THz frequency range and its expected dynamic Kerr response would decay with = ( B 2 ) −1 ≈ 0.57ps ( is the speed of light), very close to the relaxation time of the observed bipolar TKE response of water. This intermediate relaxation time has also been resolved in previous OKE experiments and similarly assigned to the relaxation of the H-bond bending of water (56,57). Moreover, our temperature-dependent TKE measurement does not contradict this assignment. As shown in Fig. S9a and Fig. S9b, the magnitude of the T-dependent TKE decay time has the closest affinity to the T-dependence of the relaxation of the H-bond bending mode (58), relative to the other relaxation processes of water.
Accordingly, the collision-induced polarizability of the dimer perpendicular to the O-H dipole becomes larger than its parallel component. As a result, the polarizability anisotropy of the dimer can be conceived as negative ΔΠ I = Π ∥ I − Π ⊥ I < 0. In this simplified picture, the dimer is the rotating entity whose response Δ rot−dimer ∝ ΔΠ I 〈 2 (cos θ dimer )〉 is resolved in the TKE experiment. Note that the dimer picture tends to be a toy model to provide a simple understanding of the bipolar TKE signal of water. Obviously more rigorous studies are required to shed light on the collective intermolecular dynamics of water (59).  Fig. 1b), which has a large spectral overlap with the THz pump spectrum, is almost temperature independent (25). To explain the sign flip of the TKE signal of water, we suggest a simple dimer model. In the highlighted dimer structure, the permanent dipole is given by the O-H group of molecule A. Upon the coupling of the THz electric-field to the latter dipole, the energy of this interaction is transferred to the molecule B and increases its kinetic energy (facilitated by the H-bonding between the two neighbors). As discussed in the text above, for the excitation frequency ~1 THz, the lateral motion of molecule B perpendicular to the OH dipole (along the normal mode coordinate of H-bond bending mode), is more likely to gain larger kinetic energy rise, relative to the translational motion along the OH dipole. As a result, the electron cloud of molecule B with excess kinetic energy is deformed, because of its collision with adjacent neighbors. The latter deformation of electron cloud would then be the origin of the larger polarizability of the dimer perpendicular to the OH dipole relative to the polarizability component along the OH dipole axis, such that ΔΠ I = Π ∥ I − Π ⊥ I < 0. Accordingly, the TKE response of this dimer Δ rot−dimer ∝ ΔΠ I 〈 2 (cos θ dimer )〉 is negative.