Nuclear quantum effects on the vibrational dynamics of the water-air interface

We have applied path integral molecular dynamics simulations to investigate nuclear quantum effects on the vibrational dynamics of water molecules at the water-air interface. The instantaneous fluctuations in the frequencies of the O-H stretch modes are calculated using the wavelet method of time series analysis, while the time scales of vibrational spectral diffusion are determined from frequency-time correlation functions and joint probability distributions, as well as the hydrogen bond number correlation functions. We find that the inclusion of nuclear quantum effects leads not only to a redshift in the vibrational frequency distribution by about 120~cm$^{-1}$ for both the bulk and interfacial water molecules, but also to an acceleration of the vibrational dynamics on the water-air interface by as much as 60$\%$. In addition, a blueshift of about 30 ~cm$^{-1}$ is seen in the vibrational frequency distribution of interfacial water molecules compared to that of the bulk. Furthermore, the dynamics of water molecules beyond the topmost two molecular layers was found to be essentially similar to that of bulk water.


I. INTRODUCTION
Understanding the behavior of water at interfaces is a primary interest of the biological, chemical and earth sciences [1][2][3].An asymmetrical solvent environment exists at the interface, which leads to variation in the dynamics, structure, reactivity and selective adsorption at the water-air interface [4][5][6][7][8][9][10].This asymmetrical local solvent structure at the interface is characterized by the hydrogen bonding of the interfacial water molecules.The temporal and spatial reorganization of the local hydrogen bond network is intimately connected with the instantaneous fluctuations in vibrational frequencies of the OH/OD modes in liquid water [11,12].Recent experimental techniques like vibrational sum-frequency generation (vSFG) spectroscopy measure the spectral contributions from the non-centrosymmetric region of a solution and thus has enabled a direct surface-specific probe of the molecular structure and dynamics of interfacial water molecules [13][14][15][16][17][18].However, the interpretation of the observed vSFG spectrum is complicated, as there is no direct way for uniquely decomposing the observed peaks.In this context, several molecular dynamics (MD) studies have been performed to interpret the relationship between the vibrational spectrum and the vibrational dynamics of interfacial water molecules [19][20][21][22][23][24][25][26][27][28][29][30][31][32].Previous simulation studies have shown that the hydrogen-bonded network becomes less structured at the interface, with an enhanced population of single donor (SD) configurations, which is basically attributed to the sharp peak around 3700 cm −1 [8,33].Hynes and coworkers employed a mixed quantum/classical approach for describing vibrations and reported the first simulated vSFG spectrum of an water-air interface [19,20].Along similar lines, Skinner and coworkers used an empirical map-based quantum/classical approach to determine the vSFG spectrum of the water-air interface [24][25][26][27][28].The surface-specific velocity-velocity correlation function approach is also extensively used to obtain the vSFG spectrum of the aqueous solutions [29,34].Quite recently, fully DFT-based ab-initio MD studies have also reported the layer-wise signal contribution to vSFG for the water-air interface [35,36].However, because water consists of light hydrogen atoms, the inclusion of zero-point energy and quantum tunneling effects is important for understanding the microscopic behavior of aqueous solutions and interfaces [37].
The path-integral molecular dynamics (PIMD) scheme is known to incorporate nuclear quantum effects (NQEs) within atomistic simulations [38,39].Earlier studies have shown that the inclusion of NQEs leads to a less structured first coordination shell accompanied by an enhanced order in the second coordination shell of liquid water [37].Similarly, the diffusion coefficients were reported to show an increase of 50 % owing to NQEs [40].However, recent PIMD-based studies have shown that a competition between the inter-and intramolecular NQEs exists, so the overall effects can be significantly less [41,42].Interestingly, the contribution of NQEs also varies, depending on the dynamical property being investigated.While the reorientational relaxation times see an increment of 15 % [43], the vibrational frequency correlation decay is accelerated by nearly 30 % [44].In the present work, we have investigated the role of NQEs on the vibrational dynamics of water molecules at the water-air interface.We have also determined the extent to which the effect of interfacial effects persists by determining the layer-specific vibrational spectral diffusion.The timescales of decay of the frequency-time correlation function (FTCF) are correlated with the hydrogen bond reorganization using the temporal decay of R O•••O fluctuation decay and hydrogen bond number correlation functions.Moreover, the time-dependent frequency joint probability distributions, frequency structure correlation and vSFG spectrum are also computed.

II. COMPUTATIONAL DETAILS
In the present study, we have conducted PIMD simulations in the canonical NVT ensemble consisting of 216 water molecules using the flexible q-TIP4P/F water model developed by Habershon and coworkers [41].Earlier studies have shown the flexible q-TIP4P-F model to be well suited to determine the impact of NQEs on the structural, thermodynamical, and dynamical, as well as spectroscopic properties of liquid water [41,42,44].The water-air interface was generated by creating a rectangular simulation box of 18.64 Å in x-and y-direction and 93.22 Å in z-direction.The system was replicated periodically along all three dimensions using the minimum image convention.Short-range interactions were truncated at 9 Å and the Ewald summation mehtod was employed to determine the long-range electrostatic interactions.The ring-polymer contraction scheme with a cutoff value of σ = 5 Å was used to reduce the computationally expensive part of electrostatic forces calculation to a single Ewald sum [41].While a p = 32 ring polymer bead was employed, the computationally expensive electrostatic calculations were contracted to the centroid only.In contrast to the original PIMD scheme, the partially adiabatic centroid MD (PACMD) technique supports the evaluation of dynamical quantities within the PIMD framework [45].The effective mass of the ringpolymer beads is adjusted by modifying the elements of the Parrinello-Rahman mass matrix so as to recover the correct dynamics of the centroids and allow for an integration time-step close to the ionic resonance limit.The temporal evolution of the ring-polymer was performed analytically in the normal mode representation by a multiple time-step algorithm using a discretized time-step of 0.5 fs for the intermolecular and 0.1 fs for the intramolecular interactions [46].The systems were initially equili-Figure 1. Definition of the different layers using the ITIM method [47,48].The surface or topmost layer is first identified.Then, the procedure is repeated for the second and third layers.Finally, the inner molecules (transparent) are shown.The definition used in this work uses the topmost layer as L1, surface and second layers as L2, first three layers as L3 and, finally, the whole system is considered as bulk water.This figure was produced using VMD [49].brated for 20 ps, followed by production runs of 100 ps each.For comparison, an additional simulation with classical nuclei, i.e. p = 1 was performed.

A. Structure and reorientational dynamics
The identification of different layers of the water/air system was done using a scheme for the identification of truly interfacial molecules (ITIM) [47,48].This algorithm uses a probe sphere to detect the molecules at the surface.The radius of the probe sphere was set as 0.2 nm, which has been proven to be a good value for water [48].A distance-based cluster search was also performed using cutoff value of 0.35 nm, which corresponds to the first minimum within the oxygen-oxygen radial distribution function.This search prevents the identifying of evaporated molecules as surface ones.Once the topmost layer is identified, the whole procedure is repeated, thus being capable of labeling each molecule belonging to the successive layers.We have used three layers in this work.Fig. 1 shows the definition of the layers for subsequent analyses: L1 corresponds to the surface molecules, L2 considers both the surface and the immediate second molecular layer below the surface, whereas L3 includes surface molecules and molecules in two consecutive molecular layers.Finally, the whole system is considered as bulk water.Once the surface molecules are identified for every frame of the simulation, the intrinsic mass density profile of water is obtained according to the equation where z i is the position of the i-th particle and ξ(x i , y i ) contains the location of the water-air interface in the zdirection.Fig. 2 shows the resulting intrinsic profile including the ones of the first three consecutive molecular layers.The first layer shows a Dirac delta peak at z = 0, given that all molecules in this layer are at the surface.This figure shows rich features up to about 4 Å , which would have been spread out in an averaged density profile.We note that the inclusion of NQEs does not alter the intrinsic density profile at the interface, which is in agreement with our previous work [33].
We have analyzed the local structure at the water-air interface and bulk using Voronoi polyhedra analysis [50].The Voronoi region of an atom is the region in space to which all the points are closer than to any other atom.The distribution of Voronoi polyhedra volume for water molecules at layer L1 and for the bulk are shown in Fig. 3 for both classical MD and PIMD simulated systems.It is evident for the distributions that there is a significant increase in the volume of Voronoi polyhedra at the interface in comparison to the bulk, which in principle implies that water molecules are less compactly packed at the water-air interface than in bulk.Furthermore, the distribution volume of Voronoi polyhedra of interfacial water molecules obtained from PIMD is also relatively wider than from classical MD simulations.The reorientational dynamics of water molecules at the surface, i.e. layer L1, is analyzed using the orientational time correlation function, which is given as where µ(t) is the orientation of the molecular dipole vector at time t and P 2 is the second-order Legendre polynomial.The integrated time scale of the orientational time correlation function, which is shown in Fig. 4, gives the average time of the loss of correlation within the reorientational dynamics.For the interfacial molecules in layer L1, the time scales of the orientational correlation function of the molecular dipole vectors are found to be 1.93 ps for the classical MD and 1.11 ps for the PIMD simulations.A recent study has shown that the water molecules on the water-air interface with free OH modes have a reorientational relaxation timescale of 0.81 ps [51], which is in agreement with our PIMD results.

B. Frequency fluctuation dynamics
The time-dependent fluctuations in ground state frequencies of OH modes were calculated from PACMD trajectories using the wavelet method of time-series analysis.The details for the calculation of OH frequencies have been described elsewhere [52,53].Here, we present a brief overview of the method.The underlying principle of the wavelet method employed here is that a timedependent function can be expressed in terms of translations and dilations of a mother wavelet where the coefficients of the wavelet expansion are given by the wavelet transform of f (t), i.e.
with a and b being real quantities with a> 0.
We have employed the so-called Morlet-Grossman form for the mother wavelet [56].The inverse of the scale factor a is proportional to the frequency over a time window around b.The time-dependent function f (t) is constructed as a complex function, with its real and imaginary parts corresponding to instantaneous OH bond distance and momentum of an OH mode projected along the OH bond.This method has been applied to all the OH modes present in a given system.
The distribution of OH frequencies as a function of position within the water-air interface has been analyzed for both systems and shown in Fig. 5 (a) and (b), respectively.We observe that the distribution of OH stretch frequencies corresponding to the layer L1 shows a distinct shoulder in the high-frequency regime of 3800-3900 cm −1 from classical simulations.The shoulder is less prominent in the case of PIMD simulations, although an overall blueshift is still observable.The mean frequency of OH stretch modes in layer L1 obtained from classical MD and quantum mechanical PIMD simulation are 3640 cm −1 and 3518 cm −1 , respectively.Similarly, the frequency of OH stretch modes averaged over molecules in layer L2, comprising the two topmost layers of the waterair interface, as obtained from the classical and PIMD simulations are 3620 cm −1 and 3502 cm −1 .The distribution of OH stretch modes in the L3 layer is nearly indistinguishable with that of the bulk for both of the cases, with the average frequencies for L3 being 3612 and 3496 cm −1 , wheras the mean frequencies for all the OH modes are 3608 and 3493 cm −1 , respectively.Purely from a static perspective, this implies that the water-air interface is characterized by the two topmost layers of water molecules.Moreover, a blueshift of ∼ 25-30 cm −1 is observed in the mean frequency of OH modes present in layer L1 as compared to bulk water.A persistent redshift of 122 cm −1 is seen for the interfacial, as well as bulk water molecules upon inclusion of NQEs.
Water molecules in liquid state are known to exhibit strong electronic asymmetry, which is a consequence of difference between the strength of first and second strongest hydrogen bond formed by a solvated water molecule [54,55].Since, the hydrogen bond strength is directly related to the vibrational frequencies of the OH modes in liquid water, we have analyzed the asymmetry between the intramolecular OH modes of water molecules in the bulk, as well as in layer L1 of the interface.The vibrational asymmetry is defined as the difference between the two vibrational frequency of the intramolecular OH modes of a given water molecule.We note that a vibrational asymmetry of 125 cm −1 exists in bulk as obtained from our classical MD simulation.As the water molecules directly at the water-air interface are exposed to a more non-uniform solvent environment, the observed asymmetry gets amplified to 150 cm −1 .Interestingly, including NQEs by means of quantum mechanical PIMD simulations, the observed asymmetry between the intramolecular OH modes in bulk is with 140 cm −1 slightly enhanced with respect to the corresponding classical simulations.Similarly, with the inclusion of NQEs, the interfacial water molecules experiences only a marginal increase in asymmetry of just 10 cm −1 as compared to its classical counterpart with 150 cm −1 , as shown in Fig. 6.Moreover, the joint probability distribution functions of the stronger versus the weaker intramolecular OH mode for bulk and interfacial water, as obtained by the present classical MD simulations are also in Fig. 7.The observed asymmetric probability distribution further supports that the intramolecular OH mode asymmetry is more prominent at the water-air interface as compared to bulk.
Similarly, the dynamical evolution of vibrational fre-  quencies of OH modes present in different layers of the water-air interface is determined using FTCFs.The timedependent decay of FTCFs for different layers and the bulk are shown in Fig. 8 (a) and (b), respectively.The FTCFs for all cases show an initial fast decay which extends up to the first 100 fs, followed by a slow inertial decay which decays within a few picoseconds.The timescales of correlation decay were obtained from the raw data by using the following biexponential fitting function The timescale of decay of the long-tail component for layer L1, L2, L3 and bulk obtained from the classical simulations are 2.37, 1.98, 1.70 and 1.50 ps, respectively.In the case of PIMD, the timescale of decay of the long-tail component for layer L1, L2, L3, and bulk are 1.45, 1.31, 1.20 and 1.02 ps, respectively.It is evident that the dynamics of OH modes present in layer L1 with reference to bulk is slower by nearly 38 and 58 % for classical and PIMD simulations, respectively.Furthermore, the decay of layer L3 and bulk are found to be similar for both the initial fast and long-tail components in both simulations.
It is indicative of the fact that from both a static and dynamical perspective, the water molecules beyond the top two layers are indistinguishable from bulk molecules.With reference to the NQEs, it is known that the vibrational dynamics of liquid water become accelerated by 30 % on inclusion of NQEs.However, for interfacial water molecules, these effects can lead to nearly 63 % acceleration.Along similar lines, we have also calculated the joint probability distribution of vibrational frequencies for OH modes of layer L1 and bulk obtained from the PIMD simulation for different waiting times of 250, 750 and 1500 fs, as shown in Fig. 9.These distributions relax qualitatively at the timescale corresponding to the timescale of loss of frequency correlation.It is evident from the plot that for bulk corresponding to the waiting time of 1500 fs frequency, the distribution becomes completely delocalized and spherical along both of the frequency axes.However, for the layer L1, the distribution still shows some asymmetric alignment along the diagonal, indicating slower vibrational spectral diffusion.
The instantaneous fluctuations in the vibrational frequencies of the OH modes in liquid water are strongly connected to the modulations in the local hydrogen bond structure.Accordingly, we have calculated the timedependent decay of fluctuations in , for the OH modes in the layers L1, L2, L3 and bulk for the classical and PIMD simulations, as shown in Fig. 10.The timescales are calculated by the least-squares fit to a biexponential function in analogy to Eq. 5.The longtime decay constant τ 1 is attributed to the timescale of loss of local structural correlation.The timescales of loss of structural correlation for the layers L1, L2, L3 and bulk obtained from classical simulations are 3.50, 2.70, 2.29 and 1.95 ps, respectively.Explicitly including NQEs leads to an overall accelerated dynamics, with resultant timescales for L1, L2, L3 and bulk being 2.40, 2.06, 1.90 and 1.64 ps, respectively.We observe that although the dynamics are slow compared to the decay of frequency correlation, the overall trends remain similar.Furthermore, the NQE contribution to the acceleration of vibrational dynamics of interfacial water molecules in layer L1 is 45 %.
While demonstrating that the temporal evolution of vibrational frequencies of OH modes is correlated with the instantaneous fluctuations in the R H−O•••H distance both in the bulk and at the interface, we also analyzed the spatial correlation between the modulations in vibrational frequencies as a function of water molecules in bulk and at interface, as shown in Fig. 11.Deviations of the OH modes at the water-air interface from the bulk behavior are easily seen in the region corresponding to higher vibrational frequencies.The increased spread of hydrogen bond network that is less correlated with the fluctuations of the OH stretch modes in comparison to the bulk.
The frequency correlation decay of water molecules is strongly influenced by the reorganization of the local hydrogen bond network.In order to correlate the vibrational correlation loss with the hydrogen bond network, we have calculated the hydrogen bond number correlation function, which shown in Fig. 12, and defined as where δn(t) = n(t)− n and n(t) is the number of hydrogen bonds that a water molecule forms at time instant t.
The timescales of decay of the layer-specific correlation functions were determined using a biexponential leastsquares fitting function of Eq. 5.The longer timescale τ 1 is attributed to the loss of correlation in the average hydrogen bond number.The timescales for the layer L1, L2, L3 and the bulk from our classical MD simulations were found to be 3.40, 2.73, 2.40 and 1.90 ps, respectively.The inclusion of NQEs leads to an overall accelerated dynamics, with timescales for L1, L2, L3 and bulk water of 2.12, 1.66, 1.45 and 1.27 ps, respectively.Furthermore, here we also observe an acceleration by a factor of nearly 60 % with respect to the vibrational dynamics at the interface.
While the vibrational correlation loss is strongly correlated with the reorganization of the hydrogen bond network, the increased timescales of vibrational spectral diffusion at the water-air interface in comparison to bulk indicates significant differences in time-averaged hydrogen bonding patterns at the interface.We have calculated the population of water molecules in different hydrogen bonding states, namely the (i) DD state (with both OH modes hydrogen-bonded), (ii) SD state (with one of the OH modes hydrogen-bonded and the other dangling), and (iii) ND state (with free OH modes).The population distribution of DD, SD and ND states in the bulk from the classical simulation was found to be 0.74,  0.23 and 0.01, respectively.With NQEs, the population distribution for the bulk is modified to 0.61, 0.34 and 0.04, respectively.We note that explicitly accounting for NQEs leads to a similar population distribution, but an evidently less structured hydrogen bond network, as seen through the increase in the population of the SD state.The population distribution of DD, SD and ND states at the water-air interface obtained from our classical MD simulations is 0.56, 0.40 and 0.03, respectively.Upon considering NQEs, the population of the three states becomes 0.43, 0.47 and 0.09, respectively.Thus, we note that the increased population of the SD state is responsible for the overall slowed down decay of frequency correlation at the interface.With reference to the inclusion of NQEs, we find that the majority of water molecules are predicted to be in the SD state.The population distribution for the bulk and the interface are both shown in Fig. 13.
As the population distribution is increased at the interface, we calculated the configuration-wise frequency correlation decay for all OH modes present in the system.The OH modes were classified into two classes: free OH (fOH), i.e.SD and ND configurations, and bonded OH (bOH), which correspondins to the DD state.
The timescale of decay of FTCFs were obtained using the biexponential fitting function, as shown in Eq. 5.
The timescale of decay of fOH and bOH for the PIMD simulated system were 1.20 and 0.74 ps, respectively.Similarly, for the classical case, the timescales obtained for fOH and bOH modes were found to be 2.78 and 1.00 ps, respectively.We note a striking similarity in the timescales of decay of fOH modes and the timescales of decay of the FTCF for layer L1, which evidently implies that the slowdown of vibrational relaxation is predominantly due to the increased population of fOH modes on the surface that characteristically have a slower decay of frequency correlation.All results for the fOH and bOH modes are shown in Fig. 14.The correlation decay of all modes is also shown for comparison.
The vSFG spectrum of the water-air interface is calculated for the layers L1, L2 and L3 for both of the systems, as shown in Fig. 15.The vSFG was calculated using the surface-specific velocity-velocity correlation function approach [29][30][31], i.e.where χ abc 2 is the resonant component of susceptibility, while r OH j and ṙOH j refer to the intramolecular distance and velocity of a given OH mode, respectively.With respect to the contributions of the individual layers to the vSFG intensity, it is seen that as we move below the water-air surface, the contributions become negligible, which is in agreement with the decay trends of FTCF for the three layers.Furthermore, the explicit inclusion of NQEs leads to an apparent redshift in the peak positions and an overall broadening of the peaks corresponding to the free and bonded OH modes in comparison to the classical MD simulations.

IV. SUMMARY
To summarize, the impact of NQEs on the vibrational dynamics of water molecules directly at the liquid-vapor interface is investigated in detail by means of PIMD simulations.We find that the explicit inclusion of NQEs leads to a reduction of frequency correlation timescales by as much as 60 %.While an overall reduction in timescales is easily discerned, a significantly faster decay in the ini-tial diffusive regime is also noticeable.Furthermore, it is also found that the static vibrational distribution, as well as the temporal decay of frequency correlation of water molecules beyond the two topmost layers of molecules is essentially identical to that of bulk water.A redshift of 120 cm −1 is noticed for both the interfacial and bulk water molecules upon the inclusion of NQEs.The water molecules on the surface show a blueshift due to the increased population of SD configurations at the interface.The timescales of decay of the hydrogen bond number correlation and local structure correlation also follow the observed decay pattern of frequency correlation loss.The computed vSFG spectra leads to similar inferences.

Figure 2 .
Figure 2. Intrinsic density profile for the molecular layers on the surface as obtained from our classical MD (top panel) and PIMD (bottom panel) simulations.

Figure 3 .Figure 4 .
Figure 3.The Voronoi polyhedra distribution for water molecules present in (a) layer L1 obtained from PIMD, (b) layer L1 obtained from classical MD, (c) bulk obtained from PIMD, and (d) bulk obtained from classical MD simulations.

Figure 5 .
Figure 5. Vibrational frequency distributions of the O-H stretching modes in layer L1, L2, L3 and the bulk as obtained from (a) classical MD and (b) PIMD simulations.

Figure 6 .
Figure 6.Vibrational frequency distributions of the stronger (OH1)) and weaker (OH2) intramolecular OH modes of water molecules in (a) layer L1 and (b) bulk, as obtained by classical MD simulations, as well as (c) layer L1 and (d) bulk from quantum mechanical PIMD simulations.

Figure 7 .
Figure 7. Joint probability distribution functions of the stronger (OH1)) and weaker (OH2) intramolecular OH modes of water molecules in (a) layer L1 and (b) in bulk, as obtained by classical MD simulations.

Figure 8 .
Figure 8. Time-dependent decay of the frequency-time correlation function of OH modes present in layer L1, L2, L3 and bulk, as obtained by (a) classical MD and (b) quantum mechanical PIMD simulations.

Figure 9 .Figure 10 .
Figure 9. Joint probability distributions of finding the O-H stretching frequencies ω1 and ω3 at a time delay of t2 for water molecules in layer L1 (top) and bulk (bottom) obtained from PIMD simulations.

Figure 11 .
Figure 11.Frequency structure correlation between the instantaneous vibrational frequency of the O-H stretch mode as a function of RH−O•••H , RO•••O and the angle θO−H•••O for layer L1 (top panel) and bulk (bottom panel) obtained from PIMD simulations.

Figure 12 .
Figure 12.The time-dependent decay of the hydrogen bond number correlation function of water molecules in layer L1, L2, L3 and bulk, as obtained from (a) classical MD and (b) PIMD simulations.

Figure 13 .
Figure 13.The population distribution of water molecules in different hydrogen bonding states in (a) layer L1 and (b) bulk obtained from classical MD simulation, as well as (c) layer L1 and (d) bulk obtained from PIMD simulations.

Figure 14 .
Figure 14.Time-dependent decay of frequency correlation of free OH (fOH) and bonded OH (bOH) modes, as obtained from (a) classical MD and (b) PIMD simulations.The timeaveraged decay of all OH modes (in black) is for comparison.

Figure 15 .
Figure 15.vSFG of OH modes present in layer L1, L2 and L3 determined from (a) PIMD and (b) classical MD simulations.