Understanding the coherent dynamic structure factor of liquid water measured by neutron spectroscopy with polarization analysis: a Molecular Dynamics simulations study

. This work is focused on atomistic molecular dynamics (MD) simulations of water carried out at 300 K. The main goal is to better understand the experimental results of the coherent dynamic structure factor S ( Q , n ) of D 2 O that were obtained by means of neutron scattering with polarization analysis and previously reported by us [A. Arbe et al. Phys. Rev. Res. 2 , 022015 (2020)]. From the simulations, we have calculated the coherent dynamic structure factor in the time domain S ( Q , t ) as well as its self-and distinct-contributions. We have also calculated S ( Q , t ) corresponding to a H 2 O sample. The main results obtained are: (i) The Q -independent relaxation process identified in S ( Q , n ) in the mesoscopic range ( Q 0 -mode) is the responsible of the restructuring of the hydrogen bond (HB) network at times shorter than that corresponding to the molecular diffusion; (ii) the vibrational contribution identified at high frequency in S ( Q , n ) corresponds to a hydrodynamic-like mode propagating in an elastic medium (fixed HB bonding pattern); (iii) in the crossover range from mesoscopic to intermolecular scales, diffusion also progressively contributes to the decay of density fluctuations; (iv) MD-simulations suggest that it would be basically impossible to measure S ( Q , n ) of H 2 O in the mesoscopic range with the current neutron scattering capabilities.


Introduction
Nowadays there is not any doubt that water is a liquid that plays a prominent role in nature.Moreover, water is also one of the most interesting liquids concerning its physical and chemical properties.It is believed the peculiar behaviour of water is mainly due to the fact that water is an associated liquid, showing a network of hydrogen bonds (HB) that breaks and reforms continuously.
The molecular dynamics in general of liquid water has been deeply investigated by means of many different techniques, including neutron scattering.In particular, the collective dynamics has been previously investigated by inelastic neutron scattering (INS) [1][2][3][4][5], quasielastic neutron scattering (QENS) [6], inelastic Xray scattering (IXS) [5,[7][8][9] and molecular dynamics (MD) simulations [10][11][12][13][14].However, in spite of this effort, the collective dynamics of water is still poorly understood, in particular, at the mesoscale, i.e., at length scales larger than the intermolecular ones but not yet in those corresponding to the hydrodynamic range.The coherent structure factor (, ) ( being the frequency) measured by IXS was usually described in terms of a viscoelastic model that involves a 'structural' relaxation process, usually called 'a-relaxation'.However, due to the current resolution of this technique (≳ 1.5 meV) and to the shape of the resolution function, a direct * juan.colmenero@ehu.eusobservation of this relaxation process has remained elusive.In principle, QENS is the ideal technique for directly observing the 'structural' collective relaxation of water in a wide Q-range, due to the available high energy resolution.However, the problem is that the measured intensity always contains a combination of coherent (collective) and incoherent (self) contributions.In samples like water with a high content of hydrogen atoms, it is generally assumed that replacing hydrogen atoms by deuterons (i.e., measuring D 0 O) the measured neutron intensity would be basically coherent.Although this is true for D 0 O in the Q-range of the first maximum of the static structure factor S(Q), this is not the case at the mesoscale.In that Q-range (0.1Å 34 ≲  ≲ 0.7Å 34 ) it has been demonstrated that the incoherent intensity is of the order of the coherent one [15].The only way to unambiguously separate coherent and incoherent neutron scattering contributions is by using scattering of polarized neutron beams and analysing the polarization of the scattered neutrons.However, this technique has also many technical difficulties, in particular, for neutron spectroscopy with sub-meV resolution.The first measurements of quasielastic neutron scattering with polarization analysis on D 0 O, in a wide Q-range from mesoscopic to intermolecular scales, were recently reported by us in Ref. [15].They were obtained by means of LET direct geometry time of flight TOF spectrometer at the ISIS Neutron and Muon Source, Oxfordshire, UK.In that paper, it was clearly proved the separation between coherent and incoherent contributions to the measured scattering.The obtained coherent and incoherent data were analyzed in terms of the imaginary part of the corresponding susceptibility, "(, ).In general, "(, ) shows a low-frequency relaxational-like process and a high frequency vibrational-like contribution.Fig. 1 displays the Qdependence of the relaxation times obtained by fitting the low-frequency peak of "(, ) by a simple Lorentzian-like expression --equivalent to a simple exponential in the time domain.As can be seen, the incoherent relaxation time follows the expected behaviour for a diffusion-like Q-dependence (~ 30 ) in the low-Q range.However, the coherent relaxation time hardly depends on Q (within the experimental uncertainties) in the low-Q mesoscopic range.This implies a non-diffusive decay of density fluctuations in such a range.In this short paper, we will focus on MDsimulations results trying to contribute to the understanding of the experimental results mentioned above.In particular, we would like to answer the following questions: (i) why the coherent relaxation in the mesoscopic range is not diffusive?(ii) what is the physical origin of such a process?On the other hand, the experimental results we are dealing with correspond to D 0 O.We note that we tried to measure the coherent contribution in the mesoscale of a H 0 O sample without success.Then, the questions emerging are: (iii) what should be observed if we were able to measure separately the coherent scattering on H 0 O sample? and (iv) are this kind of measurements possible with the current neutron scattering capabilities?
After the Introduction (Section 1), the paper is organized as follows.First (Section 2) we will give a summary of the details of MD-simulations carried out.In Section 3 we will describe the relevant neutron scattering magnitudes to be calculated from the MDsimulations atomic trajectories.After that, in Section 4, we will describe and discuss the main results obtained.Finally, the main conclusions will be summarized in Section 5.

Molecular Dynamics Simulations
MD-simulations described in this work were carried out by means of the facilities implemented in the GROMACS module.The simulated cubic cell constructed was large enough (100555 water molecules) to facilitate the calculation of collective properties in the mesoscopic range.After the equilibration procedure and the needed NPT (fixed number of molecules, pressure value and temperature) dynamics, the side of the equilibrated cubic cell at 300K and 1 atm was L = 144.56Å.The production runs were carried out under the canonical ensemble, NVT, conditions (fixed number of molecules, cell side value and temperature).The water model used was that of TIP4P-EW.In order to evaluate the electrostatic-like terms of the potential, we used the Ewald algorithm.Different production runs were performed allowing to extend the simulation time up to 10ns.The simulated cell allowed us to calculate the coherent dynamic structure factor, S(Q,t), for Q-values in the range 0.1 ≲  ≲  <=> ≈ 2 Å 34 where  <=> means the Q-value corresponding to the first maximum of the static structure factor of D 0 O.More details about the simulation procedure can be found in Ref. [14].

Calculation of Magnitudes Measured in Neutron Scattering Experiments
From the atomic trajectories obtained in the MDsimulation runs we have calculated the partial dynamic structure factors defined as: where we have assumed isotropic behaviour.We note that in the case of water, ,  ∈ {H, O}.  BC (, ) in expression ( 1) is the van Hove correlation function that in the classical limit and for an isotropic sample is defined as where r is the radial distance from a given atom and  ⃗ `B e ⃗ aC f is the position of the  hi ( hi ) atom of type (). B ( C ) is the number of nuclei of kind () and N is the total number of atoms.The brackets in (2) mean the ensemble average that considers a large number of atomic configurations.
In neutron scattering experiments, the partial dynamic structure factors are weighted with the neutron scattering lengths of the involved isotopes  m B ,  m C .Thereby, the coherent dynamic structure factor in the time domain, as should be experimentally accessed, is obtained as: For t=0, this expression gives the static structure factor S(Q).In the case of water (,  ∈ {H, O}) expression (3) reads as: where  m n =5.803fm and  m p =-3.741fm.
As it has already been mentioned, the MDsimulations were carried out with the model corresponding to light water (H 0 O).Then, in order to calculate S(Q,t) for deuterated water (D 0 O) we have just used the same partial dynamic structure factors in expression (4) but changing  m p by  m q =6.671fm.We have shown that this approximation works rather well [14].By considering a shift factor of 0.74 in the frequency scale of the simulations, an excellent agreement was found between simulated data of the imaginary part of the susceptibility for D 0 O and those measured by neutron scattering with polarization analysis on an actual D 0 O sample [14].
It is worthy of remark that the van Hove correlation function  BC (, ) contains the so-called self-part of the van Hove correlation function.This is obtained restricting the correlations included in expression (2) to those relating the positions of a single particle of type () at different times.As consequence, the partial dynamic structure factor,  BC (, ), for  =  can also be decomposed in 'distinct' and 'self' parts.For example,  BB (, ) =  BB rstu (, ) +  BB v`rh (, ).Then, for water, S(Q,t) given by expression (4) can be rearranged as (, ) =  rstu (, ) +  v`rh (, ) where, for light water,  rstu (, ) =  m n 0  nn rstu (, ) +  m p 0  pp rstu (, ) and  v`rh (, ) =  m n 0  nn v`rh (, ) +  m p 0  pp v`rh (, ) + 2 m n  m p  np (, ).In the framework of our approximation, for heavy water it would be the same but replacing  m p by  m q as it was mentioned above.

Results and Discussion
First of all, we would like to point out that the MDsimulations considered in this work were previously validated by direct comparison with the neutron scattering results [S(Q) and the imaginary part of the susceptibility "(, )] as it has been mentioned above and described in a previous publication [14].Before focusing on the different partial dynamic structure factors,  BC (, ), Fig. 2 shows both, coherent and incoherent dynamic structure factors, as they are calculated from the simulations for Q-values in the range 0.1Å 34 ≲ Q ≲ 0.5 Å 34 .For the incoherent contribution we observe the expected behaviour in this low-Q regime: only one time-decay function (basically exponential), which characteristic relaxation time, , strongly depends on Q (~ 30 ).The coherent dynamic structure factor is very different.It shows a two-step decay, being the second one basically Q-independent in the low-Q range considered in the plot.In Ref. [14] we demonstrated that the first --Q-dependent--decay can be explained in terms of a longitudinal hydrodynamic-like component.We note that similar results were recently reported by Sciortino et al. [13] from massive MDsimulations of water but corresponding to the supercooled regime of water (T~250K) and to a rather low Q-range (0.03Å 34 ≤ Q ≤ 0.3 Å 34 ).Now we obtain similar results but for liquid water (300K) and for a larger Q-range.Even though these differences, in the rest of the manuscript we will follow the terminology introduced in Ref. [13] and we will call the Qindependent relaxation process  w -mode.We note that the experimental relaxation time also hardly depends on Q in the low Q-range (see Fig. 1).Following the procedure carried out for obtaining the relaxation times reported in Fig. 1 that was described in the introduction, we have fitted the low-Q time-domain coherent simulation data by the expression: Here, Γ is the acoustic attenuation coefficient and c the adiabatic sound velocity.Both parameters were allowed to vary with Q (molecular hydrodynamics approximation).The values obtained for the coherent relaxation time  } in the Q-range from Q=0.15 Å 34 until Q=0.8 Å 34 are shown in Fig. 1.They nicely agree with the experimental times obtained by a similar procedure but in the frequency domain to fit the susceptibility "(, ) obtained by neutron scattering with polarization analysis [15].This agreement gives additional support to our MD-simulations.Now, taking advantage of MD-simulations, we can try to understand why the relaxation process observed does not depend on Q in the low-Q mesoscopic range.To do that, we have calculated from the simulated trajectories the self and distinct parts of the partial dynamic structure factors  BC (, ) (,  ∈ {H, O}).They are shown in Fig. 3 for a representative low-Q value Q = 0.2 Å 34 .As can be seen in the figure, the distinct part of both oxygen and hydrogen atoms are negative for this Q-value.The short-time range ( ≲ 3 ps) of the hydrogen-hydrogen distinct part shows clear signatures of two processes.They are also evident in the oxygen-hydrogen cross correlations that are positive.If we magnify the function corresponding to the oxygenoxygen distinct correlations, we realize that these shorttime processes are also there.On the other hand, both oxygen-oxygen and hydrogen-hydrogen correlations show a long-time decay to zero.This long-time decay is the only process in the self-part of both oxygen and hydrogen atoms, which are positive.Doing the same calculations but for other values in the low-Q range, we realized that the long-time decay depends on Q as a diffusive process.Interestingly enough, this long-time decay is the same -apart from the sign-for both selfand distinct contributions.Thereby, they cancel each other when we make the summations to calculate  pp (0.2, ) and  nn (0.2, ).These two partial dynamic structure factors are represented in Fig. 4.This figure also contains the partial dynamic structure factor corresponding to oxygen-hydrogen correlations that, obviously, only have distinct component.As can be seen, the three partial dynamic structure factors show the same behaviour apart from the amplitude: (i) they decay to zero at about 10 344 s; (ii) they show a first decay, that can be identified with a hydrodynamic-like longitudinal mode, and a second relaxation decay; the amplitude of the  pp (0.2, ) is the same as that corresponding to 2 np (0.2, ).The fact that all these partial dynamic structure factors display the same behaviour tells us that all atoms are participating in a similar way in the two processes involved -something that could be expected taking into account the low-Q range considered.Based on these results, we understand why the coherent dynamic structure factor S(Q,t) at the mesoscale does not show any signature of diffusive processes.The reason is that the diffusive component of the self and distinct parts of the partial dynamic structure factor cancel each other.
Starting from the partial dynamic structure factors shown in Fig. 4 for Q = 0.2 Å 34 and taking into account the suitable coherent scattering lengths ( m n =5.803fm,  m p =-3.741fm and  m q =6.671fm) we can calculate the dynamic structure factor S(Q,t) [eq.( 3)].In Fig. 5 we show the behaviour obtained for heavy water (D 0 O) and  light water (H 0 O).The intensity of S(Q,t) for Q = 0.2 Å 34 and  ≅ 10 34' s results to be in the case of D 0 O about 125 times larger than that corresponding to H 0 O.We note that in the low-Q thermodynamic limit ( ⟶ 0), the coherent intensity measured in a neutron experiment on a molecular fluid can be expressed as Where n is the molecular number density and the summation runs over all the nuclei of a molecule. š is the Boltzmann constant and  oe the isothermal compressibility, which is basically the same for heavy and light water (see, e.g.[16]).Then, considering the values of the coherent scattering lengths about mentioned, we obtain that  }wi q ‡ n  }wi pn • ≃ 130 in rather good agreement with the value of about 125 obtained from the simulations (see above).This, not only gives additional support to our simulations, but also tells us that it is basically impossible to measure S(Q,t) of H 0 O by neutron scattering, even using methods based on polarization analysis.1 and 2) this implies that, in the time range of the fast relaxation decay ( } ≈2ps) the water molecules are trapped in the hydrogen bond (HB) network.Thereby, the decay of the density fluctuations by this process has to be related with the restructuring process of the HB network, likely mediated by a sequence of elementary bond-switch processes.We note that experimental results, corresponding to different temperatures in the range 280K-350K, deliver an activation energy of ≈4 Kcal/mol for  } [15].This energy is larger than the values usually associated with the enthalpy for breaking a hydrogen bond in water (≈2 Kcal/mol) [17].It is worth mentioning that signatures of this HB network restructuring processes have also been observed in the short-time behaviour of single-particle correlations as incoherent neutron scattering of H 0 O [6,15] or single water molecule reorientation [18].On the other hand, at times shorter than  } (the HB-network restructuring time) the HB network has, thereby, a fixed bonding pattern.Then, the hydrodynamic-like modes mentioned above [first decay of S(Q,t)] are propagating in an elastic medium.Since the hydrodynamic-like decay of S(Q,t) depends on Q (see Fig. 2) there will be a low Q-value at which the vibrational process takes place on a time scale comparable with the network restructuring.Only at Q-values below this crossover-Q, the hydrodynamic predictions are expected to properly model the decay of density fluctuations in the system.From our MD-simulation results we can estimate (from extrapolation) a value of the crossover Q in the range 0.03 --0.1 Å 34 , which was not directly covered by the simulations [14].
Now, once we have clarified why the diffusion processes do not play any role in the relaxation of density fluctuations in the mesoscopic Q-range and what is the physical origin of the remaining short-time relaxation process, the question is in which Q-range the diffusion processes start to be an effective mechanism for the decay of density fluctuations.To address this question we have calculated the self and distinct part of S(Q,t), at different Q-values from the mesoscopic to the intermolecular range, for both D 0 O and H 0 O, and the static part of these magnitudes as well.The results obtained for D 0 O are shown in Figs. 6 and 7. Fig. 6 displays the self-part of S(Q,t) and the distinct part but with the sign changed to facilitate the comparison.The resulting S(Q,t) is also shown in the same scale.As it was previously discussed, at Q = 0.2 Å 34 the diffusive contribution to self and distinct parts is, basically, the same but with opposite sign.Thereby, they cancel each other.Fig. 6 shows that the same happens at Q = 0.5 Å 34 and Q = 0.7 Å 34 , although, as can be seen, the characteristic time of the long-time diffusive part decreases as expected ( v ~30 ).However, for Q > 0.7 Å 34 , the absolute value of the distinct part diminishes and, as consequence, (i) the diffusive contribution to S(Q,t) starts to be more and more important as soon as Q increases (see the case of Q = 1.2 Å 34 in Fig. 6 (d)), and (ii) this contribution merges with the fast process identified at low Q.Thereby, S(Q,t) cannot be described anymore by only a hydrodynamiclike term plus a single-exponential fast process [eq.( 6)].In Ref. [15] we proposed a "convolution model" for describing the relaxation contribution to the decay of density fluctuations in this crossover Q-range from mesoscopic to intermolecular scales.In that framework, the relaxational contribution can be expressed as a convolution product: where  v (, ) = [−   v ⁄ ] is a diffusive process ( v ∝  30 ) and  } (, ) is identified with the  Q -mode above mentioned. } (, ) is parametrized as: with amplitude [1 − ()] and Q-independent relaxation time,  } Q .We note that when ()~0 and  v ≫  } Q expression (7) reduces to the simple exponential used in the low Q-range, where the  Qmode dominates the relaxation decay of density fluctuations.As soon as Q increases, atomic diffusion enters in the game and this kind of processes dominate in the Q-range of the maximum of S(Q).The idea behind the convolution model is that in the crossover range S(Q) as Q = 2.4 Å 34 , represented in Fig. 8 (b), the self and distinct parts are clearly different and thereby S(Q,t) becomes 'visible'.These MD-simulation results suggest that the only possibility of measuring the coherent dynamic structure factor S(Q,t) of H 0 O by neutron scattering -even using polarization analysis techniques--would be in the Q-range of the first maximum of S(Q), which for H 0 O takes place at about 3 Å 34 .However, in this Q-range the decay of the density fluctuations would be dominated by diffusive processes.Then, it would be basically impossible to reach experimental evidence by neutron scattering of the so-called  Q − mode in the case of H 0 O.

Summary and Conclusions
In this paper we have first given a short overview of the behavior of the collective density fluctuations in liquid water (300K) in the mesoscale and in the crossover range towards intermolecular scales.The work was focussed on MD-simulations results that were carried out in order to better understand the first neutron scattering measurements of the coherent dynamic structure factor of liquid water (D 0 O) with neutron polarization analysis that were previously performed by us [15].After the validation of the MD-simulations, we have exploited the advantages of this method to calculate the neutron scattering magnitudes in the time domain and in a broader frequency range than the experimental counterpart.We have also been able to separate the self and distinct contributions to the coherent dynamic structure factor as well as the different atomic pair-correlations involved.We have also calculated the neutron scattering dynamic coherent structure factor corresponding to a fully protonated H 0 O sample.The main conclusions obtained can be summarized as follows: -In the mesoscopic range e0.1 Å 34 ≲  ≲ 0.7Å 34 f the diffusive contribution to the self-and distinct-part of the coherent dynamic structure factor cancel each other and thereby the decay of the density fluctuations take place only by hydrodynamic-like modes and a basically Qindependent relaxation process ( Q − mode).-In the time scale of the  Q -mode, the water molecules are trapped in the HB-network (diffusion is not possible).Then, the physical nature of this relaxation mode has to be related with the restructuring of the HBnetwork, likely mediated by a sequence of elementary bond-switch processes.
-At times shorter than the time scale of the  Q -mode the HB-network has to show a fixed bonding pattern.Thereby, the hydrodynamic modes are propagating in an elastic medium.
-The density fluctuations will decay only by hydrodynamic modes in the Q-range at which the time scale of the  Q -mode is shorter than that characterizing the hydrodynamic mode decay, which is Q-dependent.It can be estimated that this crossover-Q should be in the range 0.03 Å 34 -0.1 Å 34 .-In the crossover range from mesoscopic to intermolecular scales e0.7 Å 34 ≲  ≲ 2Å 34 f, the diffusion processes do not cancel and they progressively also contribute to the decay of density fluctuations.A convolution model for the relaxational part of the decay of density fluctuations seems to be a rather good approximation in this crossover range.
-The MD-simulations results show that the amplitude of the intensity measured by neutron scattering in a H 0 O sample, is extremely low in the Q-range from ~0.7 Å 34 to ~2 Å 34 .This suggests that it would be basically impossible to measure this coherent dynamic structure factor in this Q-range with the current neutron scattering capabilities.
We acknowledge the Grant PID2021-123438NB-I00 funded by MCIN/AEI/10.13039/501100011033and by "ERDF A way of making Europe".We also acknowledge financial support of Eusko Jaurlaritza, code: and IT1566-22, as well as from the IKUR Strategy under the collaboration agreement between Ikerbasque Foundation and the Materials Physics Center on behalf of the Department of Education of the Basque Government.

Fig. 1 .
Fig. 1.Characteristic relaxation times (blue symbols: coherent; red symbols: incoherent) obtained by the procedure described in the text from the neutron scattering data with polarization analysis corresponding to a D 0 O sample at 300 K [15].Different solid symbols correspond to data obtained with different incident neutron wavelengths.Open blue squares correspond to the coherent relaxation times obtained from the MD-simulations described in the text.

Fig. 5 .
Fig. 5. Coherent dynamic structure factor (, ) corresponding to neutron scattering measurements calculated from the simulations for both D 0 O and H 0 O at Q = 0.2 Å 34 .Now we can ask what is the physical origin of the two-steps decay shown by the partial dynamic structure factors and S(Q,t) as well.Since in the mesoscopic scale the diffusion process takes place at times longer than  } (see Figs.1 and 2) this implies that, in the time range of the fast relaxation decay ( } ≈2ps) the water molecules are trapped in the hydrogen bond (HB) network.Thereby, the decay of the density fluctuations by this process has to be related with the restructuring process of the HB network, likely mediated by a sequence of elementary bond-switch processes.We note that experimental results, corresponding to different temperatures in the range 280K-350K, deliver an activation energy of ≈4 Kcal/mol for  }[15].This energy is larger than the values usually associated with the enthalpy for breaking a hydrogen bond in water (≈2 Kcal/mol)[17].It is worth mentioning that signatures of this HB network restructuring processes have also been observed in the short-time behaviour of single-particle correlations as incoherent neutron scattering of H 0 O[6,15] or single water molecule reorientation[18].On the other hand, at times shorter than  } (the HB-network restructuring time) the HB network has, thereby, a fixed bonding pattern.Then, the hydrodynamic-like modes mentioned above [first decay of S(Q,t)] are propagating in an elastic medium.Since the hydrodynamic-like decay of S(Q,t) depends on Q (see Fig.2) there will be a low Q-value at which the vibrational process takes place on a time scale comparable with the network