Unravelling the contribution of local structures to the anomalies of water: the synergistic action of several factors

TL;DR

This study uses molecular dynamics simulations to explore how local structural changes in water's hydrogen bond network contribute to its anomalous properties, especially near the Widom line and under different pressures.

Contribution

It provides a detailed microscopic analysis of water's anomalies, highlighting the role of LDL and HDL environments and their spatial organization, revealing the synergistic factors involved.

Findings

HDL-like clusters fragment at a critical LDL concentration of 20-30%.

The hydrogen bond network in LDL environments shows maximal frustration with pentagonal and hexagonal rings.

Global order decreases near the Widom line, supporting the liquid-liquid critical point hypothesis.

Abstract

We investigate the microscopic origin of water's anomalies by inspecting the hydrogen bond network (HBN) and the spatial organization of low-density-liquid (LDL) like and high-density-liquid (HDL) like environments. Specifically, we simulate --via classical molecular dynamics simulations-- the isobaric cooling of a sample composed by 512 water molecules from ambient to deeply undercooled conditions at three pressures, namely 1 bar, 400 bar and 1000 bar. \emph{In correspondence with the Widom line}, (i) the HDL-like dominating cluster undergoes fragmentation caused by the percolation of LDL-like aggregates following a spinodal-like kinetics; (ii) such fragmentation always occur s at a "critical" concentration of $\sim20-30\%$ in LDL; (iii) the HBN within LDL-like environments is characterized by an equal number of pentagonal and hexagonal rings that create a state of maximal frustration between a configuration that promotes crystallization (hexagonal ring) and a configuration that hinders it (pentagonal ring); (iv) the spatial organization of HDL-like environments shows a marked variation. Moreover, the inspection of the global symmetry shows that the intermediate-range order decreases in correspondence with the Wid om line, and such decrease becomes more pronounced upon increasing the pressure, hence supporting the hypothesis of a liquid-liquid critical point. Our results reveal and rationalize the complex microscopic origin of water's anomalies as the cooperative effect of several fac tors acting synergistically. Beyond implications for water, our findings may be extended to other materials displaying anomalous behaviours.

Equations

$$S(j)=\underset{\theta ,\varphi ,\psi ;\mathcal{P}}{\max }\left\{\prod _{i=1}^M\exp \left(-\frac{{\left|{\mathbf{P}}_{i\mathcal{P}}^j-{\mathbf{A}}^j{\mathbf{R}}_i^j\right|}^2}{2{\sigma }^{2}M}\right)\right\}$$

$$S=\frac{1}{N}\sum _{j=1}^{N}S(j)$$

$$I(i)=\frac{1}{n(i)}\sum _{j=1}^{n(i)}{\left[\Delta (j;i)-\bar{\Delta }(i)\right]}^2$$

Full text

Supercooled water, water anomalies, local structures, hydrogen bond network

††preprint: AIP/123-QED

I Introduction

In 1933, Bernal and Fowler provided the first microscopic description of the structure of liquid water and identified the presence of an underlying hydrogen-bond network (HBN) connecting neighbouring molecules in a tetrahedral-like configuration Bernal and Fowler (1933). In order to account for the observed changes of water upon cooling, Bernal and Fowler suggested that water is characterized by a continuous transformation of the HBN from a denser quartz-like structure four-fold coordinated, to a less denser tridymite structure four-fold coordinated also, where the tridymite structure is the structure of hexagonal ice (Ih). Only a few years later, in 1939, Pauling recognized that the HBN is a characteristic that differentiate water from simpler liquids Pauling (1960), a notion that is nowadays widely accepted. In 1951, People further developed the concept of HBN framing it in a random network picture, showing that the HBN can accommodate arrangements in which molecules in the second or third shells of neighbours can collapse towards the first coordination shell of a central molecule People (1951). To further enrich this emerging complex picture, in 1976, Speedy and Angell observed that, at variance with other liquids, the isothermal compressibility of supercooled water increases upon cooling, suggesting the presence of a thermodynamic singularity at -45*∘*C Angell and Angell (1976). This hypothesis was further supported by the evidence that other thermodynamic response functions also tend to diverge upon supercooling. In 1992, based on the results of classical molecular dynamics simulations, Poole et al hypothesized that a liquid-liquid critical point (LLCP) located at deeply undercooled conditions and low/intermediate pressures could be the source of the observed thermodynamic anomalies Poole et al. (1992, 1993a, 1993b). In correspondence with the LLCP, two liquids (an high-density liquid –HDL– and a low-density liquid –LDL–) are in metastable equilibrium with each other. The existence of the LLCP has been proven via numerical simulations for a molecular model of water Palmer et al. (2014), and several other computational Poole et al. (1992, 1993a, 1993b); Xu et al. (2005); Abascal and Vega (2010); Liu et al. (2010); Abascal and Vega (2011); Corradini and Gallo (2011); Li, Li, and Wang (2013); Palmer, Car, and Debenedetti (2013); Yagasaki, Matsumoto, and Tanaka (2014); Holten et al. (2014); Smallenburg, Filion, and Sciortino (2014); Smallenburg and Sciortino (2015); Ni and Skinner (2016); Pathak et al. (2016); Biddle et al. (2017); Palmer et al. (2018) and experimental Mishima and Suzuki (2002); Mishima and Stanley (1998); Mishima (2000); Fuentevilla and Anisimov (2006); Kanno and Miyata (2006); Mallamace et al. (2008); Mishima (2010); Bertrand and Anisimov (2011); Holten et al. (2012); Sellberg et al. (2014a, b); Fanetti et al. (2014a); Kim et al. (2017) investigations point toward the same confirmation. A remarkable evidence in favour of the LLCP is the polyamorphic character that water acquires at low temperatures and low/intermediate pressures. Two amorphous states (a low-density amorphous –LDA– ice, and an high-density amorphous –HDA– ice) are separated by a first-order phase boundary Mishima, Calvert, and Whalley (1985); Mishima (1994); Mishima and Suzuki (2002); Klotz et al. (2005); Winkel, Mayer, and Loerting (2011); Martelli et al. (2017, 2018a) and correspond to the glassy state of LDL and HDL. A very high-density amorphous state has also been hypothesized Loerting et al. (2001). Signatures of LDA and HDA have been found in the inherent potential energy surface (IPES, a collection of local potential energy minima resulting from systematic quenches along a liquid trajectory allowing to separate packing from thermal motion effects Stillinger and Weber (1984a, b)) of ab initio liquid water at ambient conditions, where the LDA-like and the HDA-like environments are characterized by having distinct HBNs Santra et al. (2015). The distinct HBNs in the IPES of ab initio liquid water indicate that there is an energy barrier associated with the different HBNs, and this barrier could be large enough at low temperatures and pressures to account for the first-order nature of the LDA-HDA phase boundary. Two distinct local structures have also been observed experimentally in liquid water from ambient to supercooled conditions, and interpreted as LDL-like and HDL-like environments Taschin et al. (2013) as well as under pressure Fanetti et al. (2014b). Water can then be regarded as a dynamical mixture of two different local structures –or two states– whose fraction changes with the thermodynamic conditions, accounting for the thermodynamic anomalies Tanaka (2000a, b); Russo and Tanaka (2014); Holten and Anisimov (2012); Holten et al. (2014); Nilsson and Pettersson (2015); Biddle et al. (2017); Anisimov (2018). Evidence supporting the existence of two states in liquid water has been found by Raman spectroscopy Walrafen (1964, 1967); Walrafen et al. (1986), femtosecond mid-IR pump-probe spectroscopy Woutersen, Emmerichs, and Bakker (1997), time-resolved optical Kerr effect spectroscopy Taschin et al. (2013), X-ray absorption Sellberg et al. (2014a) and emission spectroscopies Huang et al. (2009).

Within the LLCP scenario, a line of maxima of the correlation length emanates from the LLCP, at which the correlation lengths diverge. This line, known as the Widom line (WL), is the locus of maxima of the thermodynamic response functions, and where a strong to fragile crossover occurs Starr, Sciortino, and Stanley (1999); Faraone et al. (2004); Xu et al. (2005); Liu et al. (2005); Mallamace et al. (2006); Zheng et al. (2009); Gallo, Rovere, and Chen (2010); Gallo and Rovere (2012); Gallo, Rovere, and Chen (tter); Wang et al. (2015); Xu et al. (2016); De Marzio et al. (2016, 2017), suggesting that thermodynamic anomalies are intimately linked with dynamic anomalies. In this work, we define TWL the temperature range in which, at a given pressure, the isobaric heat capacity $C_{P}$ and the isothermal compressibility $\kappa_{T}$ show a maximum. In fig. 1 we report $C_{P}$ (panel (a)) and $\kappa_{T}$ (panel (b)) for the three pressures here investigated. We can observe that TWL is $\sim 220-230$K for 1 bar, $\sim 210-220$K for 400 bar and $\sim 200-205$ for 1000 bar. In water, the WL is negatively sloped, so the Clapeyron equation implies that the high-temperature phase is HDL dominated, and the low-temperature phase is LDL dominated. Recently, Shi et al provided a microscopic basis for a two-state description of water, and proposed that the emergent fragile-to strong transition is actually a crossover between two Arrhenius regimes with different activation energies, bringing into question the glass-transition scenario Shi, Russo, and Tanaka (2018a, b).

The outstanding importance of water anomalies is nowadays recognized to play a central role in processes in several fields, from physics and chemistry, to biology, material science, geology and climate modeling Eisenberg and Kauzmann (1969); Angell (1995); Mishima and Stanley (1998); Debenedetti (2003). For this reason, a lot of effort has been devoted to rationalizing water’s anomalous behaviours, but a coherent and simple picture has yet to be found.

In this article, we have performed equilibrium classical molecular dynamics (MD) simulations crossing the WL at three pressures, namely 1 bar, 400 bar and 1000 bar. Therefore, we concern ourselves only with the lines of maxima determined along isobars for the remainder of the current work111The Widom line can also be defined, e.g., on the isothermal pathway fomin_2015; schienbein_2018. We have observed an hitherto unknown intimate connection between the percentage of LDL-like molecules and their HBN, and between HDL-like molecules and their spatial organization occurring in correspondence with the WL. Our results indicate that water’s anomalies occur under specific conditions of composition and connectivity, and enlighten the role of the underlying HBN. Therefore, our results elucidate the basic microscopic physical picture hidden behind water’s anomalies, and may represent a generic feature of anomalous liquids at large.

The article is organized as follows: in Section II we describes the molecular dynamics set up and two order parameters adopted in the present work. In Section III we present our results, while final remarks are reported in Section IV.

II Simulation details

In this Section we describe the protocol we have implemented in our simulations and the two order parameters we have employed to investigate the temperature-dependence of the global order while crossing the WL and to distinguish local environments.

II.1 Classical molecular dynamics simulations

We have performed classical MD simulations of a sample of 512 water molecules interacting with the TIP4P/2005 potential Abascal and Vega (2005) in the NPT ensemble. We have employed Nosé-Hoover thermostat Nos/’e (1984); Hoover (1985) with 0.2 ps relaxation time to maintain constant temperature, and Parrinello-Rahman barostat Parrinello and Rahman (1981) with 2 ps relaxation time to maintain constant pressure. We have explored three pressures, namely 1 bar, 400 bar and 1000 bar, and isobarically scanned temperatures from 300 K down to 205-190 K depending on the pressure. We have performed simulations with the GROMACS 5.0.1 package Pronk et al. (2013). We have truncated short-range interactions at 9.5 Å, and we have computed long range electrostatic terms using particle mesh Ewald with a grid spacing of 1.2 Å. At each state point, we have computed and carefully monitored the decay of the self-part of the intermediate scattering function (ISF) with time Hansen and McDonald (2006). All reported trajectories are at least 500 times longer that the structural relaxation time as computed from the ISF. Depending on the thermodynamic conditions, production runs vary between 50 ns and 12 $\mu$s. No sign of crystallization have been observed.

II.2 The order parameters

In this section we describe the two order parameters employed in this work.

II.2.1 The score function

The local environment of an atomic site $j$ in a snapshot of a molecular dynamics or Monte Carlo simulation defines a local pattern formed by $M$ neighboring sites. Typically these include the first and/or the second neighbors of the site $j$. There are $N$ local patterns, one for each atomic site $j$ in the system. The local reference structure is the set of the same $M$ neighboring sites in an ideal lattice of choice, the spatial scale of which is fixed by setting its nearest neighbor distance equal to $d$, the average equilibrium value in the system of interest. For a given orientation of the reference structure and a given permutation $\mathcal{P}$ of the pattern indices, we define the LOM $S(j)$ as the maximum overlap between pattern and reference structure in the $j$ neighborhood by:

[TABLE]

Where $\theta,\phi$ and $\psi$ are Euler angles, $\mathbf{P}^{j}_{i\mathcal{P}}$ and $\mathbf{R}_{i}^{j}$ are the pattern and the reference position vectors in the laboratory frame of the $M$ neighbors of site $j$, respectively, and $\mathbf{A}^{j}$ is an arbitrary rotation matrix about the pattern centroid. The parameter $\sigma$ controls the spread of the Gaussian functions ($\sigma=d/4$ in this work, where $d$ is the characteristic length of the local pattern). The LOM satisfies the inequalities $0\lesssim S(j)\leq 1$. The two limits correspond, respectively, to a local pattern with randomly distributed points ($S(j)\rightarrow 0$) and to an ordered local pattern matching perfectly the reference ($S(j)\rightarrow 1$). We also define a global order parameters based on $S(j)$, as the average score function $S$:

[TABLE]

We have recently employed the LOM and the score function to enlighten structural properties of water and its homogeneous crystallization at various conditions Martelli et al. (2016, 2018a, 2018b); Samatas et al. (2018), and -as a collective variable- to drive the formation of boron-nitride nanotubes Santra et al. (2018). The details of the numerical algorithm can be found in Ref Martelli et al. (2016).

In this work, we have employed eq. 2 to measure the global symmetry of our samples upon isobarically cooling and crossing the Widom line at the three pressures here investigated.

II.2.2 The local structure index

We order the set of radial oxygen-oxygen distances ${r_{j}}$ corresponding to the $N$ neighbouring molecules that are within a cut-off distance of $3.7$ Åfrom the reference molecule as follows: $r_{1}<r_{2}<...<r_{j}<r_{j+1}...<r_{N}<3.7<r_{N+1}$. The local structure index (LSI) value $I$ is then defined as the inhomogeneity in this distribution of radial distances, i.e.,

[TABLE]

where $\Delta(j;i)=r_{j+1}-r_{j}$ and $\bar{\Delta}(i)$ is the average over all neighbours $j$ of a molecule $i$ within a given cutoff. Hence, $I$ provides a convenient quantitative measure of the fluctuations in the distance distribution surrounding a given water molecule within a sphere defined by a radius of $\sim 3.7$ Å Shiratani and Sasai (1996, 1998). In doing so, the index $I$ measures the extent to which a given water molecule is surrounded by well-defined first and second coordination shells.

In this work, we employ the LSI to characterize local environments in our samples as a function of the temperature and of the pressure.

III Results

In this section we present and discuss the main results of our work. Computational and experimental studies on the structural changes of liquid water have shown that the short range order (SRO) is only marginally affected by the change in pressure and temperature Soper and Ricci (2000); Svishchev and Kusalik (1993); Schwegler, Galli, and Gygi (2000); Saitta and Datchi (2003); Sciortino, Geiger, and Stanley (1990, 1991); Kumar et al. (2006); Franzese and Stanley (2007); Lapini et al. (2016). Rather, relevant effects occur beyond the SRO, at the level of the intermediate range order (IRO) Martelli et al. (2016). In particular, in 2000, Soper and Ricci showed that part of the second shell collapses towards the first shell increasing the local density and pointing out the importance of the second shell Soper and Ricci (2000). The role of the second shell has consequently been enlightened via several local descriptors such as, e.g., the distance of the fifth neighbour $g_{5}$ Cuthbertson and Poole (2011), the LOM described in eq. 1 Martelli et al. (2016), and the LSI Shiratani and Sasai (1996, 1998) that measures the degree of inhomogeneity between the first and the second shell. Wiktfeldt et al showed that, at the level of the IPES, the LSI distribution in liquid water has a bimodal distribution that, in correspondence with the WL, gives a 1:1 ratio between HDL-like and LDL-like environments Wikfeldt, Nilsson, and Pettersson (2011). It is worthy to remark that, with respect to Wiktfeldt et al, this work inspects the behaviour of water in the presence of thermal noise. Consequently, Russo and Tanaka recognized the existence of locally favoured structures by measuring the degree of translational order on the second shell, enlightening their role in the anomalies of water Russo and Tanaka (2014).

III.1 Structural order

III.1.1 Global order from the score function

In this paragraph, we inspect the local order of liquid water by separating the contribution of the first shell from the contribution of the second shell employing the score function described above. In Fig. 2 (a) we report $S_{\text{SRO}}$, the score function (eq. 2) for the SRO measured using, as a reference, a perfect tetrahedron. As expected, for a given temperature the $S_{\text{SRO}}$ decreases upon increasing the pressure, while it increases upon fixing the pressure and isobarically cooling. The increment of the $S_{\text{SRO}}$ upon isobarically cooling the sample is an effect caused by the reduction of the thermal energy, as also measured by the increase of the intensity of the first peak in the oxygen-oxygen radial distribution function.

In Fig. 2 (b) we report $S_{\text{IRO}}$, the score function (eq. 2) for the IRO computed using, as a reference, the anticuboctahedron that describes the spatial configuration of the oxygen atoms at the level of the second shell in Ih. With respect to $S_{\text{SRO}}$, the $S_{\text{IRO}}$ decreases upon increasing the pressure at a given temperature. Since in our maximization algorithm we let the reference structure to rescale its dimensions with respect to the local pattern Martelli et al. (2016), the increase of $S_{\text{IRO}}$ upon increasing the pressure indicates that the second shell tends to increase the local packing Martelli et al. (2018a). Remarkably, upon cooling isobarically our samples, the $S_{\text{IRO}}$ shows a minimum in correspondence with the location of the WL. We have fitted the loci of each minimum along with the two closest state points on the left and on the right, with a quadratic function at all pressures (dashed lines in fig. 2 (b)). Notably, the amplitude of the quadratic function shrinks from $\sim 8\times 10^{-6}$ at 1 bar, to $\sim 10^{-5}$ at 400 bar and $\sim 2\times 10^{-3}$ at 1000 bar. The minima in $S_{\text{IRO}}$ is indicative of the decrease in the local order caused by the density fluctuations that become more pronounced approaching the hypothesized LLCP Abascal and Vega (2010). Therefore, the shrinkage of the amplitude supports the LLCP scenario for this water model and confirm that the thermodynamic and dynamic water anomalies have a structural counterpart.

III.1.2 High-density and low-density environments from the local structure index

In order to enlighten the role played by local environments in the emergence of the anomalies of water, we have characterized LDL-like and HDL-like molecules based on the corresponding value of the LSI. As reported by Wikfeldt et al Wikfeldt, Nilsson, and Pettersson (2011) and Appignanesi et al Appignanesi, Rodriguez, and Sciortino (2009); Montes de Oca et al. (2016), the LSI distribution in the IPES is a bimodal function with an isosbestic point $I_{is}$ mostly unaffected by the thermodynamic conditions and located at $I_{is}\sim 0.13$ Å2. Thus, it provides a well defined criteria to distinguish LDL-like and HDL-like environments at supercooled conditions, as recently shown by Shi and Tanaka Shi and Tanaka (2018). In this work, since water’s anomalies occur in the presence of thermal energy, we have borrowed the $I_{is}$ from the IPES to define LDL-like and HDL-like environments at the level of the dynamical trajectories, i.e., in the presence of thermal energy 222The position of the isosbestic point in the presence of thermal energy fluctuates within the range $0.12-0.14$ Å2. Any change within this range would produce very minor quantitative changes in our results. The position of the isosbestic point also depends on the cutoff employed in the evaluation of eq. 3 accordino_2011.

The advantages of investigating local structures in the presence of thermal energy have been also recently enlightened by Russo and Tanaka Russo and Tanaka (2014) and by Saito et al Saito, Bagchi, and Ohmine (2018).

In fig. 3 we report the probability density distribution of the LSI order parameter ($P(I)$) and we begin our discussion by analyzing how $P(I)$ depends on the thermodynamic conditions. Panel (a) reports $P(I)$ computed at 1 bar in the temperature range $205\leq T\leq 235$K, panel (b) at 400 bar and in the temperature range $195\leq T\leq 230$K and panel (c) at 1000 bar and in the temperature range $190\leq T\leq 230$K. The bimodal character of $P(I)$ become more pronounced upon cooling the samples. This finding is indicative of the fact that, at high temperatures, the majority of the molecules are situated in HDL-like locally disordered environments. Upon cooling isobarically, we observe a systematic decrease in the relative population of water molecules with lower LSI values (i.e., molecules in locally disordered environments) coupled with an increase in the relative population of water molecules with higher LSI values (i.e., molecules in locally ordered environments). In the remainder of this work we will focus only on the aforementioned temperature intervals.

In fig. 4 we report the percentage of LDL-like molecules and of HDL-like molecules as a function of the temperature for the three pressures here investigated. As expected, the percentage of LDL-like environments increases upon cooling the samples at all pressures, indicating that the space between the first and the second shell of neighbours become less and less populated and the local environments progressively more tetrahedral, as a result of the reduction of thermal energy. This trend can also be visualized graphically by observing the three representative snapshots respectively obtained at $230$, $215$ and $195$ K and at 400 bar (similar pictures hold for the other pressures inspected in this work). Red spheres represent HDL-like molecules, while blue spheres represent LDL-like molecules. Remarkably, in correspondence with the WL (filled simbols), the composition of the sample relative to LDL acquire the same value in the range of $\sim 20-30\%$ for all pressures (colored stripe). Since this behaviour occur at the same concentrations of LDL-like particles at all pressures, we infer that this concentration represents a kind of ”critical mass” needed to sustain changes in the dynamical and thermodynamic behaviour. We therefore infer that this relatively narrow range of composition could be a key ingredient for a system like water to display anomalous behaviours, and a common feature for other liquids also displaying anomalous behaviours. We will provide further proofs of the ”criticality” of such concentration in the remainder of the article.

III.2 Clustering of LDL-like and of HDL-like environments

In order to delve deeper into the role of LDL-like and HDL-like environments in the emergence of the anomalies of water, we have inspected the tendency of LDL-like and of HDL-like environments to aggregate and form clusters. We here define a cluster an aggregate of at least two water molecules belonging to the same class -LDL or HDL- as measured by the LSI, and separated by an oxygen-oxygen distance shorter than $3.3$ Å.

The number of LDL-like clusters is reported in fig. 5 (a) and shows a non monotonic behaviour with maxima at temperatures slightly higher than the WL and a decreasing trend starting, upon cooling, in correspondence with the TWL. This observation indicates that, in correspondence with the TWL, the LDL-like clusters start forming an extended network further expanding upon cooling. As a confirmation, in fig. 5 (b) we report the size of the largest LDL-like cluster. We can observe that the pace of growth of the largest cluster shows two distinct regimes: a slow growth at temperatures above TWL, ad a faster growth starting in correspondence with TWL upon cooling the sample. This trend is indicative of the fact that smaller clusters are absorbed into the largest expanding cluster. At high temperatures, the low number of LDL-like clusters (fig. 5 (a)) and the low percentage of LDL-like molecules (fig. 4) indicate that the sample is composed by many LDL-like small, scattered clusters. On the other hand, in correspondence with the WL, the percentage of LDL-like molecules reaches the critical concentration of $\sim 20-30\%$, high enough for LDL-like clusters to percolate and form a large aggregate, further proving that this amount of LDL-like molecules represent a critical concentration.

In fig. 6 (a) we report the number of HDL-like clusters, while in fig. 6 (b) we report the size of the largest HDL-like cluster. At temperatures above the TWL, the system is composed by one large HDL-like aggregate spanning the entire simulation box and enclosing smaller LDL-like clusters. Upon cooling, the largest HDL-like cluster gets continuously depleted (panel (b)) at the expenses of LDL-like environments (fig. 4). Remarkably, the largest HDL-like cluster gets scattered into 2-3 smaller clusters in correspondence with the TWL (panel (a)), confirming that LDL-like clusters percolate inside the HDL-like large cluster, fragmenting it into smaller domains. Therefore, the term ”critical” referred to the $\sim 20-30\%$ amount of LDL-like molecules marks the limiting value of LDL-like environments at which these molecules form an aggregate large enough to percolate within the HDL-like network and to split it. Upon further cooling isobarically the samples, the number of HDL-like clusters increases (panel (a)) while the dimension of the HDL-like cluster keep decreasing (panel (b)), indicating that the LDL-like network is now further percolating inside the HDL-like network, creating several HDL-like patches.

In summary, by inspecting the evolution of LDL-like and of HDL-like clusters as a function of the temperature and of the pressure, we observe that the system is composed by a large HDL-like cluster at temperatures above TWL. Upon cooling, the percentage of LDL-like molecules increases and, in correspondence with TWL, it reaches a ”critical” concentration of $\sim 20-30\%$ that allows the LDL-like clusters to connect each other, percolating and fragmenting the large HDL-like cluster. Upon further cooling isobarically the sample, the LDL-like network keep expanding generating a fast-growing large LDL-like cluster. This picture is consistent with a spinodal-like decomposition, in which phase separation occur in correspondence with TWL. We have recently observed a similar behaviour during the simulated isothermal compression of LDA to produce HDA via a mild first-order phase transition Martelli et al. (2018a).

In fig. 7 we report a two-dimensional 4 Å-thick slice of our sample at 400 bar and at three temperatures: T=230K –(a)–, T$\sim$TWL=215K –(b)– and T=195K –(c)–. Blue spheres represent LDL-like environments, while red spheres represent HDL-like environments. At T=230K, the sample is permeated by the large HDL-like cluster and its underlying HBN (red sticks), embracing LDL-like molecules and clusters. Upon cooling the sample to T$\sim$TWL=215K, the main HDL-like network no longer span the entire simulation box; rather, it is now fragmented and alternated with the likewise fragmented LDL-like network (blue sticks). Finally, at T=195K we can observe that the HDL-like and the LDL-like networks are now almost independent and separated. The formation of two almost independent networks at low temperatures, allows us to speculate that the tendency of the system is to effectively unmix and separate the two HBNs, a mechanism that occur in its entirety in correspondence with the LLCP. Moreover, since local environments characterized by different structural properties diffuse with different velocities Xu et al. (2009); Mallamace, Corsaro, and Stanley (2013), the drastic changes in terms of composition, clustering and cluster growth occurring in correspondence with TWL may explain the fragile-to-strong crossover occurring in correspondence with the TWL.

III.3 Radial distribution functions and the hydrogen bond network

In the previous paragraph we have observed that the clustering of HDL-like and LDL-like networks plays a fundamental role in the anomalies of water in correspondence with the WL. We now investigate the details of the spatial organization and of the HBN for both LDL-like and HDL-like environments. The definition of HB follows Ref. Luzar and Chandler (1996). In this regard, any quantitative measure of HBs in liquid water is somewhat ambiguous, since the notion of a HB itself is not uniquely defined. However, qualitative agreement between many proposed definitions have been deemed satisfactory over a wide range of thermodynamic conditions Prada-Gracia, Shevchuk, and Rao (2013); Shi, Russo, and Tanaka (2018b).

We have inspected the spatial correlation via the two-body radial distribution functions (RDFs) computed among the LDL-like and the HDL-like water molecules only, and we have inspected the topology of the HBN using the ring statistics analysis, a tool which has been instrumental in theoretically characterizing the amorphous states LDA and HDA Marton̆ák, Donadio, and Parrinello (2004, 2005), and the HBN in the IPES of ab initio liquid water at ambient conditions Santra et al. (2015).

In fig. 8 we report gll(r), the RDFs computed among LDL-like oxygens only at 1 bar, in the temperature range $205\lesssim T\lesssim 235$K. Similar distributions characterize higher pressures. We can observe that the first and the second peak are very well separated, indicating that the space between the first and the second shells is mostly empty, hence reflecting the high tetrahedrality of LDL-like environments. Overall, the RDFs are qualitatively similar to the experimental RDF for LDA obtained using the empirical potential structure refinement of neutron diffraction data reported in Ref. Finney et al. (2002).

We now inspect the topology of the HBN connecting LDL-like molecules. In fig. 9 we report P(n), the normalized probability of having an n-folded ring, n$\in[3,10]$, for LDL-like environments only. Panel (a) shows the distributions at all temperatures at 1 bar, panel (b) the distributions at 400 bar and panel (c) the distributions at 1000 bar. Interestingly, at high temperatures LDL-like molecules are mostly connected via pentagonal rings and the P(n)’s show a maximum in correspondence with n5. We have encountered a similar distribution in the IPES of ab initio liquid water at ambient conditions Santra et al. (2015) and we have concluded that such distribution, that describes a network deviating from the hexagonal one characteristic of the amorphous state LDA, is an artifact due to the quenching at the level of the IPES. On the other hand, our results here indicate that a distribution maximized at n5 is a genuine effect, and we provide here a justification supporting this evidence: at high temperatures, LDL-like environments are low in number and they tend to form small clusters rather than arranging in open wires. Therefore, the pentagon is the spatial arrangement with hydrogen bond distances and angular strains the closest to the hexagonal geometry. Hence, the pentagon represents the most energetically favourable configuration at high temperatures. Upon cooling isobarically the samples, P(n) shifts from a distribution maximized at n5 to a distribution maximized at n6. This observation indicates that LDL-like environments can develop the expected hexagonal network characteristic of the amorphous state LDA only when their number overcome the critical concentration of $\sim 20-30\%$. This observation further stresses the importance of such critical concentration. It is worthy to remark that, although the gll(r) shown in fig. 8 are mostly independent with respect to the thermodynamic conditions and resemble –at least qualitatively– the experimental RDF of LDA, the underlying HBNs are drastically different not only with respect to LDA, but also upon minor changes in the temperature in the liquid phase. It has been recently shown that, upon rapid quenching equilibrated liquid water at different conditions, one obtains LDAs that are indistinguishable in terms of RDFs, but that are characterized by different potential energy landscapes (PELs) Giovambattista et al. (2016); Handle and Sciortino (2018). Our results suggest that the different underlying HBNs within LDL-like environments in the equilibrated liquid phases may represent the source of the different PELs in the amorphous states.

Remarkably, in correspondence with the WL, P(n) for the three pressures here investigated show an inflection point at which n5=n6. This inflection point represents a state of maximal frustration in a system composed by an equal number of hexagonal rings, known to be the precursor of crystallization, and pentagonal rings, known to frustrate against crystallization. Therefore, we infer that the inflection point at which n5=n6 in the LDL-like HBN represents another key feature, essential for the bulk system to display anomalous behaviours, and may be a general feature of liquids endowed with a LLCP displaying anomalous behaviours.

The inflection point is emphasized by the orange line in fig. 10, where we report the ratio between the number of hexagonal and pentagonal rings for the three pressures here examined in the temperature range $190\lesssim T\lesssim 240$K. It is possible to observe that the inflection point occurs in correspondence with the TWL (filled symbols) at all pressures.

We now focus our attention on the HDL-like environments. At variance with the LDL-like environments, the ghh(r) show a marked temperature dependence that propagates beyond the SRO (fig. 11). In fig. 11 we report the RDFs computed among HDL-like oxygens only at 1 bar, in the temperature range $205\lesssim T\lesssim 235$K. Similar distributions characterize higher pressures. Upon cooling the sample, we observe an increase in the intensity of the first peak accompanied with a deepening of the first minimum located at $\sim 3.2$ Å. Moving beyond the SRO, we observe a marked peak at $\sim 3.7$ Å, whose intensity increases on cooling and whose position is mostly independent on the temperature. On the other hand, the consequent minimum only slightly decreases on cooling, but its position shifts from $\sim 3.8$ Åat $T=235$K to $\sim 3.9$ Åat $T=205$K. This peak accounts for the presence of interstitial molecules that populate the space between the first and the second shell of neighbours, and becomes more pronounced upon cooling, suggesting that the interstitial configuration become favourable upon reduction of the thermal energy. Interestingly, the intensity of this interstitial peak is comparable with the intensity of the main second peak at the lowest temperatures. The main second peak also shows a marked temperature dependence, becoming more pronounced upon cooling. The ghh(r) become indistinguishable at radial distances of $r>6$ Å, confirming that the HDL-like environments play and active role in the appearance of water anomalies.

In order to investigate in deeper detail the role of HDL-like environments, we have focused our attention on the interstitial molecules. While an inspection of the HBN suggests that the interstitial molecules are bonded to the central molecules at all thermodynamic conditions 333It is worthy to mention that in LDA, the interstitial molecules are not bonded to the central one Finney et al. (2002); Martelli et al. (2017), suggesting that the HBN of LDA and LDL are substantially different, as reported in Ref. Martelli et al. (2017), we can acquire interesting information by computing the ratio, as a function of the temperature, between the intensity of the maximum of the interstitial peak (max${}_{\textit{int}}$) and the corresponding minimum (min${}_{\textit{int}}$). In fig. 12 we report the ratio max${}_{\textit{int}}$/min${}_{\textit{int}}$ computed at the three pressures and in the temperature range $190\lesssim T\lesssim 235$K.

It is possible to recognize a change in slope in max${}_{\textit{int}}$/min${}_{\textit{int}}$ for all pressures as a function of the temperature. This change in slope is emphasized with dashed lines that serve as a visual guide. Remarkably, the change in slope occurs for all pressures in correspondence with the TWL, suggesting that the spatial organization of HDL-like molecules plays an active role in water’s anomalies. It is worthy to emphasize here that our analysis are limited by the relatively small samples, and by the small number of HDL-like molecules located at the interstitial configuration. Further inspections with larger statistics may help in better quantifying and framing the -hitherto unknown- role of interstitial molecules in the anomalies of water.

In fig. 13 we report P(n), the

normalized probability of having an n-folded ring connecting HDL-like environments only via their HBN. Panel (a) shows the temperature dependent distributions at 1 bar, panel (b) at 400 bar and panel (c) at 1000 bar. At variance with the LDL-like network, the P(n) for the HDL-like network does not show any interesting feature in correspondence with the WL. At T$>$TWL, the distributions show a maximum at n6 and are broad, indicating the existence of longer rings to accommodate the higher density. Upon cooling the samples, the corresponding P(n) show a decrease in the number of longer rings and a corresponding increase in n5, reflecting the decrease in the thermal energy and the progressive fragmentation of the large HDL-like cluster into smaller domains.

IV Conclusions

We have explored the structural properties and the HBN of equilibrated supercooled liquid water described by the TIP4P/2005 interaction potential. We have inspected the behaviour of the SRO and of the IRO by employing the score function, and we have characterized local environments based on their LSI value.

We have inspected the composition, clustering and HBN of LDL-like and of HDL-like environments and we have shown that, in correspondence with the WL:

(i) the percentage of LDL-like molecules is in the narrow window of $\sim 20-30\%$ at all the inspected pressures, indicating that this window may represent a ”critical” concentration at which water shows anomalous behaviours;

(ii) a large, dominating HDL-like cluster undergoes fragmentation caused by the increasing number of LDL-like clusters that percolate and connect each other, in a spinodal-like decomposition scenario. We speculate here that such kinetics could be modeled with a modified Cahn-Hilliard equation Cahn and Hilliard (1958), describing a phase separation occurring in correspondence with TWL. We are currently making progresses in this direction;

(iii) the topology of the HBN within the LDL-like molecules is characterized by an equal number of pentagonal and hexagonal rings, creating a state of maximal frustration;

(iv) the ratio max${}_{\textit{int}}$/min${}_{\textit{int}}$ shows, upon cooling, a rapid change in slope, suggesting that the interstitial arrangement in HDL-like environments plays a central role in the emergence of water’s anomalies;

(v) the IRO -as measured by the score function- shows a minimum in correspondence with TWL that decrease upon increasing the pressure, reflecting the enhanced density and entropy fluctuations upon approaching the hypothesized LLCP.

Our results indicate that water’s anomalies are the combined effect of a delicate balance between HDL-like and LDL-like environments, how they cluster, how they form networks and how they organize in space. We have unbundled these effects, shedding light on the complex physics of water’s anomalies. Beyond implications for water, our results may help in understanding the anomalous behaviour of other materials.

Acknowledgements.

This work was supported by the STFC Hartree Centre’s Innovation Return on Research programme, funded by the Department for Business, Energy and Industrial Strategy.