Temperature dependent anomalous fluctuations in water: shift of ≈1 kbar between experiment and classical force field simulations

Here we report on the temperature dependence of the anomalous behaviour of water in terms of (i) its growth in tetrahedral structures, (ii) instantaneous spatial correlations from small angle x-ray scattering (SAXS) data, (iii) estimates of thermodynamic response functions of isothermal compressibility and (iv) thermal expansion coefficient. Water’s thermal expansion coefficient is estimated for the first time at supercooled conditions from liquid water’s structure factor. We used previously published data from classical force-fields of TIP4P/2005 and iAMOEBA to compare experimental data with molecular dynamics simulations and observe that these force-fields underestimate water’s anomalous behaviour but perform better upon increasing pressure. We demonstrate that the molecular dynamics simulations can describe better the temperature dependent anomalous behaviour of ambient pressure water if simulated at 1 kbar. The deviation in anomalous fluctuations in the simulations is not restricted to ≈228 K but extends all the way to ambient temperatures. GRAPHICAL ABSTRACT


Introduction
Water is one of the most anomalous liquids at conditions where essential physical, chemical, biological, and geological processes occur in nature [1]. It is a key component of life and a lot of effort has been made to understand it's peculiar properties [1][2][3][4][5][6][7][8]. In particular, water exhibits intermolecular structure of water -the nature of the Hbonding network and its collective fluctuations [10][11][12][13][14].
There are indications that the structure of water is related to structural fluctuations around two local configurations, high-density liquid (HDL) and low-density liquid (LDL) [1,8,11,15,16]. The latter corresponds to a tetrahedral structure similar to what is seen locally in hexagonal ice, whereas the former is characterised by a collapse of the tetrahedral structure allowing for higher density [17] featuring distortions in the H-bonding of the first coordination shell. LDL and HDL local configurations have also been denoted by S and ρ states, and these were shown to have a temperature and pressure dependent population in a way that described many of the water properties [12,18,19]. The anomalous behaviour of water is explained by the transient and competing nature of these structures where LDL is preferred at supercooled temperatures and HDL at higher temperatures [11,19]. Upon cooling water from ambient conditions, Hbonding becomes stronger and LDL population increases at the expense of HDL resulting in higher volume and entropy fluctuations until these fluctuations reach a maximum. Soper et al. [7] extracted the pair distribution functions for pure phases of HDL and LDL using trends in temperature and pressure from neutron scattering experiments.
What ignited a large interest is that the power-law fits of thermodynamic response functions appear to diverge towards a temperature (T S ) of about 228 K [5,20]. The scenario that best describes most of the existing experimental and theoretical data are a hypothesis that at high pressure and low temperature the LDL and HDL exist as pure phases, separated through a liquid-liquid transition (LLT) line that ends in a liquid-liquid critical point (LLCP) at positive pressure [21] or very close to ambient pressure [19]. In this scenario, the Widom line, defined as the locus of correlation length maxima in the P-T plane, emanates from the LLCP as a continuation of the LLT line into the one-phase region [22]. Maxima in isothermal compressibility for supercooled water has been recently observed experimentally by Kim et al. [23] as well as Holten et al. [24]. The locus of points based on maxima in volume or entropy fluctuations may deviate slightly at regions away from a critical point, but they eventually coincide at a critical point, as also shown for simulated supercritical water [25].
In the present study, we discuss at what pressure the LLCP would be located if it truly exists, based on a comparison of recent x-ray scattering data against different force field molecular dynamic (MD) simulations. We will base our estimate on measurements of the tetrahedral structure, small angle x-ray scattering (SAXS) data and the thermodynamic response functions of isothermal compressibility and thermal expansion coefficient as function of temperature. The LLCP has not been seen based on experiments but some water models show clear evidence of its existence, in particular TIP4P/2005 [26,27], iAMOEBA [23,28], E3B3 [29], WAIL [30] and ST2 [31]. The response functions become higher in magnitude and more sensitive to pressure and temperature close to the critical point. We investigate the iAMOEBA and TIP4P/2005 force-fields of water since they accurately predict the temperature dependence of tetrahedrality of real water [28]. Coincidentally, they also have almost similar LLCP location of around 1700 bars and 182 K for TIP4P/2005 [27] and around 1750 bar and 185 K for iAMOEBA [28]. We demonstrate that the simulated results give a qualitative close resemblance to the experimental data at ambient pressures if the simulations are conducted at ≈ 1 kbar. This deviation is not only with respect to the enhancement close to the Widom line but extends all the way up to ambient temperatures.

Tetrahedral structure
Kim et al. [23] measured small and wide angle x-ray scattering (SAXS and WAXS) of supercooled water droplets at PAL-XFEL in a similar manner as the previous work done at LCLS [32] and observed the maxima in κ T and correlation length for both H 2 O and D 2 O. It is well known [33][34][35] that the strongest feature of the structure factor in water is split into two peaks denoted by S 1 at peak position, q 1 ≈ 2 Å −1 and S 2 at peak position, q 2 ≈ 3 Å −1 where q is the magnitude of the scattering vector. The distance between the two peaks, q 2 -q 1 = q is correlated to water's tetrahedrality, defined by the second peak at r ≈ 4.5 Å of the O-O pair correlation function, where a high value of q corresponds to a high degree of tetrahedrality [28,32].
The measurement at PAL-XFEL [23] had a limited q-range where it was possible to measure only the first peak, S 1 but not the second peak. This was because the experiment was optimised for small angle scattering and the detector covered a q-range from 0.12 to 1.92 Å −1 . As shown earlier, q increases on supercooling due to a decrease in q 1 as well as an increase in q 2 . However, the decrease in q 1 is the major contributing factor for q for temperatures T < 300 K [32]. A plot of q 1 vs T is shown in Figure 1 where we can see that q 1 decreases with lowering of the temperature. The data in Figure 1 summarises various measurements done at different x-ray facilities such as Linac Coherent Light Source (LCLS) [32,37], Stanford Synchrotron Radiation Lightsource (SSRL) [32], Pohang Accelerator Lab XFEL (PAL-XFEL) [23] and  [32], the PAL-XFEL data are taken from ref. [23] and APS data are taken from ref. [36]. The figure in the inset is S(q) simulated from the iAMOEBA force-field model of water for the temperature range of 218 K (blue) to 298 K (red) in steps of 20 K and it clearly shows the sharp variation of peak position of S 1 with temperature.
Advanced Photon Source (APS) [36]. The measurements at these facilities are slightly offset from each other. This is due to errors in absolute q-calibration. The SSRL data have been obtained on a static liquid sample, experiments at PAL-XFEL and LCLS were performed with a microdroplet train whereas the APS data are obtained on different sample environments, which include acoustically levitated droplets, flowing water stream and water in thin glass capillaries. Although the absolute value of q 1 may be slightly inaccurate, the relative change of q 1 position and hence, its derivative with respect to temperature, is accurate within the same dataset.
The experimental data are compared to simulation results from the iAMOEBA and TIP4P/2005 water model. Details of the simulation are given in ref [28]. There is a good agreement in q 1 between experiments and TIP4P/2005 but an offset when compared to data from iAMOEBA at T > 255 K. At lower temperatures, the experimental data are closer to that from iAMOEBA rather than TIP4P/2005. Overall, we observe that q 1 is more sensitive to temperature at supercooled temperatures rather than at ambient temperatures. A more accurate way to quantify this is by measuring the slope, (∂q 1 /∂T) P . Figure 2 compares the measured temperature dependent (∂q 1 /∂T) P from the PAL-XFEL data at 1 bar in comparison to results from molecular dynamics simulations at different pressures. We observe that (∂q 1 /∂T) P increases as temperature is decreased until a temperature of 229.2 K where it reaches a maximum of 0.0093 Å −1 K −1 and then decreases with decrease in temperature [23]. On the other hand, simulation data from TIP4P/2005 and iAMOEBA at 1 bar have a much broader maximum in a temperature range of 221-235 K and the maximum value is about half compared to experiments. For the experimental data, the existence of a maxima in (∂q 1 /∂T) P and the temperature where it occurs is similar to the existence of maxima in κ T [23,24]. With q 1 being correlated to water's tetrahedrality, the maximum of (∂q 1 /∂T) P indicates a maximum in temperature dependence of growth of tetrahedral structures. For the iAMOEBA model at higher pressures, we see that the maxima of (∂q 1 /∂T) P moves to lower temperatures and the value of the maximum increases as we increase the pressure. We observe that the maximum for the isobar of 1000 bar has a value of around 0.01 Å −1 K −1 and is located at ≈ 203 K. This value of maxima is similar to the experimental data measured at 1 bar. The rise in (∂q 1 /∂T) P becomes sharper with increase in pressure. Experimental data have a resolution of ≈ 0.5 K near the maximum of (∂q 1 /∂T) P but the simulation data have only a resolution of 5 K (iAMOEBA) and 10 K (TIP4P/2005).

Small angle X-ray scattering
Small angle x-ray scattering (SAXS) measures the structure factor at small q values. For a liquid, it has a sensitivity to the instantaneous density fluctuations [15,20,39]. SAXS is a volumetric probe and a slope towards low q is obtained if the instantaneous density over a distance d is non-uniform, where d = 2π/q. Figure 3(a) shows temperature dependent SAXS data [23] and it is observed that the structure factor starts to increase at low values of q as the temperature decreases and becomes rather high at 231 K. What is most interesting is that all the structure factor curves cross each other around q ≈ 0.4 Å −1 indicating an isosbestic point. The question is if the simulations would capture such an isosbestic point and thereby describe the long-range structure of the fluctuating liquid as the temperature varies.
In order to calculate the low q structure factor from simulations, a large box size is necessary. This was done in the past with TIP4P/2005 simulations using a box size of 45000 molecules [38]. Figure 3(b,c) show the temperature dependent SAXS curve for TIP4P/2005 simulations at 1 bar and 1 kbar, respectively [38]. In all SAXS curves there is a decrease in the structure factor at q values above 0.5 Å −1 with decreasing temperature. However, if there is a structural heterogeneity that is increasing upon cooling, there will be a stronger enhancement at lower q. This structural heterogeneity is due to a result of tetrahedral structures surrounded by patches of Hbond distorted structures. The average length-scale of such structures is around 2 nm near the Widom line [23], which is close to the length scales of 2π/q probed at low q. We see an isosbestic point when the loss of contrast at high q and increase of contrast at low q exactly balance each other allowing for a q where the intensity is constant. For the 1 bar simulation as shown in Figure 3(b), clearly the development of structural heterogeneity is not strong enough and thereby no isosbestic point is seen. On the other hand, the 1 kbar simulation illustrated in Figure 3(c) has a narrow q-range where the S(q) curves cross close to the experimental isosbestic point indicating that the structural heterogeneities develop close to real water at 1 bar. The fact that there is a narrow q-range in 1 kbar simulations rather than an isosbestic point would indicate that the loss of contrast at high q and increase of contrast at low q closely rather than exactly balance each other. Again, this represents a shift of ≈ 1 kbar between SAXS experiments and force field simulations to obtain a good agreement.

Thermodynamic response functions
Isothermal compressibility (κ T ): Thermodynamic response functions represent how a system reacts to changes in temperature and pressure. κ T is a thermodynamic function that measures the relative change in density on application of pressure at constant temperature. It is proportional to local density fluctuations. Kim et al. recently observed a maxima in κ T [23] in supercooled water along an isobar of 1 bar. Although some aspects of this study were not unanimously accepted by other researchers [40,41] in this field, such a maximum value in κ T has also been observed in other experiments of stretched water by Holten et al. [24] and has been seen in molecular dynamics simulations of the iAMOEBA [28,42], E3B3 [29] and TIP4P/2005 [38,43,44] and ST2 [45] water models for different pressures. It was recently pointed out that when talking about maxima in the P-T plane, it is important to define the path taken. This aspect is not generally discussed in supercooled water literature and was recently pointed out in recent studies [46][47][48]. In all of the cases discussed here, the maxima are always determined along the isobar at a given pressure.
κ T increases with decrease in temperature from ambient conditions until it reaches a maximum, interpreted as the Widom line shown in Figure 4. Note that the maximum in the experimental data are seen in the inset of Figure 4. Clearly we observe that the simulations at 1 bar underestimate the sharpness of the maximum and its absolute value whereas the 1 kbar simulation data are close to the experimental data. Späh et al. [49] quantified the slope of κ T with respect to temperature and fitted the experimental data as well as the iAMOEBA water model with a power law where = (T − T S )/T S is the reduced temperature, T S is the singularity temperature and it is close to the temperature of the observed maxima of isothermal compressibility. γ in Equation (1) correlates with the slope of κ T with respect to temperature where a steep slope corresponds to a high value of γ . The behaviour of κ T with respect to temperature gets sharper and the value of the maximum is increased as we approach a critical point from a single-phase region, which in this case, means increase in pressure. We observe in Figure 5 that γ increases as a function of pressure. When we fit experimental data to Equation (1), we get γ = 0.40 ± 0.01 and it is shown as a red dashed line in Figure 5. This value lies between the value for isobars of 500 and 1000 bar for the iAMOEBA model [49].
The value of the compressibility maximum (κ T ,max ) increases as a function of pressure for different force-fields of water [23,29,43] and is shown for iAMOEBA water model in Figure 5. The experimental value of κ T ,max for H 2 O is ≈ 105 × 10 −6 bar −1 and is shown as a black dashed line in Figure 5 and it corresponds to pressures between 700 and 1000 bar in simulations [46]. So, we conclude that there is a good agreement between experiment and simulations if simulations are shifted by ≈ 0.7-1 kbar.
Thermal expansion coefficient (α): Here we discuss this property based on a new analysis where we estimate the thermal expansion coefficient (α). Thermal expansion coefficient is proportional to a combined effect of density and entropy fluctuations. Although for most liquids, density and entropy fluctuations are positively correlated, but for liquid water at 1 bar and colder than 277 K, these are anti-correlated. This results in a negative thermal expansion coefficient for supercooled water. Here we attempt to relate the structure factor of water to its thermal expansion coefficient since the decrease in density originates from the growth of tetrahedral structures with its open network. In the section 'Tetrahedral structure', we discussed about the observation of a maxima in (∂q 1 /∂T) P .
(∂q 1 /∂T) P can be related to thermal expansion coefficient (α) through Equation (2) 2) (∂q 1 /∂T) P is an experimentally measured quantity that represents the change in structure with respect to temperature. (1/ρ) × (∂ρ/∂q 1 ) P represents the relative change in bulk density with respect to structure and has not been measured experimentally but its trend can be extracted from MD simulations.
The density of supercooled water and α has been measured before down to a temperature of 240 K by Hare and Sorenson [50] in glass capillaries. This density data were fit with a 6th-degree polynomial with a very low standard deviation of 2 × 10 −4 g/cc. However, extrapolation of the fit to lower temperatures was not recommended by Hare and Sorenson [50]. Wölk et al. [51] conducted her studies on water nucleation at temperatures less than 223 K and found that such extrapolation of the density at T < 223 K results in densities lower than 0.91 g/cc, which is lower than densities of amorphous and hexagonal ice. Wölk et al. [51] found these density estimates unreasonable and found a way to extrapolate density from Kell's data [52] to 218 K without the predicted density being lower than that of amorphous ice by using her own fit function. These different extrapolations of density for our temperature range of interest of T > 227 K from Wölk et al. [51] and Sorenson et al. [50] give density values differing by no more than 3% but differ almost by a factor of 5 in α. For best estimates of α at 227 K < T < 240 K, we rely on estimates of (1/ρ) × (∂ρ/∂q 1 ) P rather than extrapolated or estimated density fits. Figure 6 shows (1/ρ) × (∂ρ/∂q 1 ) P for iAMOEBA and TIP4P/2005 simulations as a function of temperature as well as q 1 . Figure 6(a) shows that there is a trend of continuous rise in (1/ρ) × (∂ρ/∂q 1 ) P with decrease in temperature. There are no sharp changes in (1/ρ) × (∂ρ/∂q 1 ) P corresponding to the maxima in (∂q 1 /∂T) P . This means that the temperature dependence of α is mostly determined by the temperature dependence of (∂q 1 /∂T) P , which is an important observation. One can also plot (1/ρ) × (∂ρ/∂q 1 ) P with respect to q 1 as shown in Figure 6(b). Figure 6(b) shows that (1/ρ) × (∂ρ/∂q 1 ) P increases with decrease in q 1 and for a given value of q 1 , there is a general trend of higher (1/ρ) × (∂ρ/∂q 1 ) P for higher pressures.
A linear fit to (1/ρ) × (∂ρ/∂q 1 ) P for the 1000 bar isobar from Figure 6(b) is used together with Equation (2) in order to calculate α. We choose the 1 kbar isobar since it shows a close similarity to experimental (∂q 1 /∂T) P as shown in Figure 2. The resultant α is shown as 'esti-mate1000bar' in Figure 7. For comparison, α estimated from a linear fit from the isobar of 500 and 1200 bar from Figure 6(b) is also shown as 'estimate500bar' and 'esti-mate1200bar' respectively. In Figure 7, we illustrate that α shows a clear minimum for isobars of molecular dynamics simulations for the TIP4P/2005 38 and iAMOEBA force fields of water at ambient pressure. For iAMOEBA force field, the minimum shifts to lower temperatures as the pressure is increased and its location in terms of temperature coincides with the maxima observed in isothermal compressibility [23]. The estimates of α from 500, 1000 and 1200 bar are close to each other and we observe that α rapidly decreases with decrease in temperature until 229.2 K. Although this is an approximate way of determining α, the temperature dependence of α is accurately captured by these estimates. Figure 7 also illustrates the α determined by Hare and Sorenson [50] who fitted their density data to a 6th-degree polynomial equation down to 239.14 K. The temperature dependence of α is much steeper at colder temperatures and Figure 7. Thermal expansion coefficient (α) of water. Simulation data for iAMOEBA is shown for isobars of 1, 500, 1000 and 1200 bar respectively. TIP4P/2005 data are from ref. [38] (pink). The dashed line is the measurement from Hare and Sorenson [50]. 'estimate500bar', 'estimate1000bar' and 'estimate1200bar' are our estimations of α based on measured x-ray scattering spectra and estimated (1/ρ) × (∂ρ/∂q 1 ) P from iAMOEBA curves in Figure 6 it is not accurately captured by ambient pressure isobars of iAMOEBA or TIP4P/2005 force-fields. However, the temperature dependence of α for 1 kbar iAMOEBA is close to our estimates. We see a rise in α for the two coldest temperatures where T < 229.2 K. However, we are aware that the possible error in our estimates of α is too high to definitely say if there exists a minimum at 229.2 K. Figure 8 summarises the hypothesis that is consistent with all current experimental data [23]. There could exist two separate liquid phases, HDL and LDL, with a phase coexistence line in the P-T diagram deep in the supercooled regime and at elevated pressure [2,21]. This liquid-liquid transition (LLT) line terminates in the LLCP and its equivalent in the one-phase region corresponds to the Widom line [53]. At the Widom line, the density fluctuations would reach a maximum consistent with equal population of molecules in HDL and LDL [54]. Note that HDL is on the ambient temperature side of the Widom-line, whereas LDL is on the low temperature side. This explains why the high density, H-bond distorted species dominate at ambient conditions. The fluctuating anomalous region around the Widom line becomes wider on decreasing the pressure thereby, moving further away from the LLCP. The closer we are to the critical point, the narrower the region becomes and the stronger the anomalous behaviour will be, leading to steeper rise in the response functions. What is unique with water is that at ambient pressure, the location of the LLCP is such that the anomalous region with fluctuations extends up to around 320 K. At 320 K, water's isothermal compressibility starts to increase with decreasing temperature.

Discussion
All experimental data presented here shows evidence of stronger anomalous fluctuations in comparison to the two force fields iAMOEBA and TIP4P/2005 used in the current study. Similar arguments have also been given regarding the E3B3 model [46]. What we note is that the deviation is not only for the coldest temperatures but all the way to ambient temperatures. The two vertical arrows in Figure 8 indicate the pressure relative to the LLCP of real water and simulated water. This means that real water is more anomalous than most MD water models because it is closer to the LLCP. Most likely this means that the fluctuating regions are less well defined in the simulations in comparison to real water [11]. There are further evidences for this in some recent Xray spectroscopy investigations of ambient water with other comparisons between experiment and simulations [55] as well in recent intermediate scattering function measurements using speckle visibility spectroscopy [56]. It was observed that the timescale for the cage effects in real water becomes much shorter than what is seen in the TIP4P/2005 model, which is also reflected in the temperature dependence of the isothermal compressibility. Since the anomalous fluctuations are observed to be correlated with the growth of tetrahedral structures that form patches of around 1-2 nm size [11,20,23,39] depending on temperature, we assume that collective many body forces are essential to improve the simulation models [11].
In order to quantify the differences in anomalous fluctuations observed in simulations as compared to experiments, we have used the fact that the simulations show an enhancement in anomalous behaviour with increasing pressures towards the LLCP. The closer water is to a critical point, the stronger critical fluctuations are to be expected. Here we observe that the simulations qualitatively agree with the experimental data when shifted by ≈ 1 kbar. iAMOEBA and TIP4P/2005 have their critical point at 1750 [28] and 1700 bar [27] respectively. This means that for real water the LLCP, if it exists, would be located at ≈ 700 bar. This is consistent with recent estimates of LLCP, which were also derived by comparing experiments to simulations [23,46,49].

Conclusion
We have looked at experimental data of supercooled water and compared it to molecular simulations of two different force-fields of water, namely TIP4P/2005 and iAMOEBA. We have focused on the temperature dependent anomalous fluctuations related to the growth of tetrahedral structure, spatial correlation obtained from SAXS, thermodynamic response functions including isothermal compressibility and thermal expansion coefficients. Water's thermal expansion coefficient has been estimated for the first time at deeply supercooled conditions from liquid water's structure factor. Based on our observations, we conclude that ambient pressure water is best represented by molecular dynamics forcefields when simulated at ≈ 1 kbar pressure. This means that classical force field simulations underestimate the anomalous fluctuations in water extending from the Widom line to ambient temperatures.