---
abstract: '[Liquids generally become more ordered upon cooling. However, it has been a long-standing debate on whether such structural ordering in liquid water takes place continuously or discontinuosly: continuum vs. mixture models. Here, by computer simulations of three popular water models and analysis of recent scattering experiment data, we show that, in the structure factor of water, there are two overlapped peaks hidden in the apparent “first diffraction peak”, one of which corresponds to the neighboring O-O distance as in ordinary liquids and the other to the longest periodicity of density waves in a tetrahedral structure. This unambiguously proves the coexistence of two local structural motifs. Our findings not only provide key clues to settle long-standing controversy on the water structure but also allow experimental access to the degree and range of structural ordering in liquid water. ]{}'
author:
- Rui Shi
- Hajime Tanaka
bibliography:
- 'fsdpWater.bib'
title: |
    Distinct signature of two local structural motifs of liquid water\
    in the scattering function
---

Water is ubiquitous in our planet and plays vital roles in many biological, geological, meteorological, and technological processes. Despite its simple molecular structure, water shows many unique thermodynamic and dynamic properties in the liquid state, such as a density maximum at 4 $^\circ$C, a rapid increase of isothermal compressibility and a dynamic fragile-to-strong transition upon cooling [@Debenedetti2003; @gallo2016water]. These unusual properties, which are absent in ordinary liquids, are well-known as “water’s anomalies”. It is widely believed that the anomalies are linked to water’s structural ordering towards tetrahedral structures stabilized by four hydrogen bonds (H-bonds). Even after intensive studies for more than a century, however, how such structural ordering takes place is still a matter of hot debate without convergence.

Two conflicting different scenarios have continued to exist until now: ‘continuum models’ based on a broad unimodal distribution of structural components and ‘mixture models’ based on a bimodal distribution of structural components reflecting the coexistence of two (or more) types of local structures [@narten1969observed; @eisenberg2005structure; @handle2017supercooled]. The mixture model dates back to Wilhelm Röntgen, who proposed in 1892 that water can be regarded as ’icebergs’ in a fluid ’sea’ [@rontgen1892ueber]. Later various mixture models have been developed. One famous example is the mixture model of Linus Pauling, who proposed in 1959 that water is mixture of clathrate-like structure and interstitial molecules [@pauling1959structure]. These mixture models, however, have been continuously challenged by the continuum model dating back to John Pople, who proposed in 1951 that water’s structure can be described by a continuously distorted H-bond network [@pople1951molecular].

These two types of models lead to fundamentally different understandings of water structure. Despite such a clear difference in the physical picture, there has been no convergence of this debate over a century. The main reason is the lack of [*experimental*]{} evidence exclusively supporting either of the two models. In this Letter, we provide such clear evidence that liquid water is indeed a mixture of two types of local structural motifs, from simulations of three popular water models and detailed analysis of recent scattering experiments.

First we need to explain the precise nature of our two-state model to clarify essential differences from ‘continuum models’ and other types of mixture models that regard water as a mixture of two types of liquids, i.e. low-density (LDL) and high-density liquids (HDL). Our two-state model that regards water as a mixture of ordered ($S$) and less ordered local structural motifs ($\rho$-state) [@Tanaka_review; @Tanaka2000] is characterized by the following five features: (1) $S$- and $\rho$-states in liquid water are defined by [*local structures around a central molecule*]{} and characterized by low and high local symmetry, energy, density, and entropy, respectively. In the one-phase regime of water far from the second-critical point (if it exists), the two structural motifs can have only short coherence lengths. Thus, our model is essentially different from a mixture model of LDL and HDL. We stress that they are macroscopic phases of water that can exist only below the second critical point; (2) Reflecting the presence of the two states, the distribution of a proper structural descriptor should have a bimodal distribution composed of two Gaussians (not necessarily two delta functions; note that there is no unique configuration for each state under thermal fluctuations); (3) The two-state model effectively transforms to a continuum-like model at high temperatures/high pressures where there exists only $\rho$ state, because the ordered $S$-structure can hardly survive due to the entropy/volume penalties; (4) The $T, P$-dependence of the fraction of the two structural motifs should obey the thermodynamic two-state equations [@Tanaka_review]. (5) The existence of a second critical point is a sufficient but not necessary condition for the two-state model.

Recently we have shown that we can detect two structure motifs by a microscopic structural descriptor $\zeta$ (see Methods), which characterizes the translational order in the second shell, and confirmed the above five features on a microscopic level by computer simulations of several popular water models [@Russo2014; @shi2018impact; @shi2018Microscopic; @shi2018origin; @shi2018common]. These studies have clearly indicated that water is a dynamic mixture of the two states [@Tanaka_review; @Tanaka2000; @holten2013nature; @Russo2014; @singh2016two; @Biddle2017; @Singh2017; @de2018viscosity; @russo2018water; @shi2018impact; @shi2018common; @shi2018origin]—$S$-state \[locally favored tetrahedral structure (LFTS)\] and $\rho$-state \[disordered normal-liquid structure (DNLS)\]. The former stabilized by four H-bonds has lower symmetry, density, energy and entropy than the latter. A typical snapshot of LFTS and DNLS is shown in Fig. \[fig:fnfsZeta\]a. We have also found that the fraction $s$ of the LFTS, serving as an order parameter characterizing the degree of structural ordering, changes with temperature $T$ and pressure $P$, obeying the prediction of the thermodynamic two-state model [@Tanaka_review; @Tanaka2000; @Singh2017; @de2018viscosity; @russo2018water; @shi2018common; @shi2018origin]. These results provide strong computational support for the two-state model.

The shape of the distribution function of a physical quantity is a key to judge which of mixture and continuum models is relevant, because the former predicts a bimodal distribution at a certain range of $T$ and $P$ whereas the latter always predicts a unimodal Gaussian distribution. Our structural descriptor $\zeta$ clearly shows the bimodality. Then, a key question is whether a quantity directly related to local density shows such bimodality or not. The answer is yes. We show the distribution $P(N_\mathrm{fs})$ of the coordination number $N_\mathrm{fs}$ in Fig. \[fig:fnfsZeta\]b. We note that $N_\mathrm{fs}$ is the number of water molecules in the spherical first-shell volume $V_{\rm fs}$ of radius of 3.5 Å, and thus proportional to the local number density, $N_{\rm fs}/V_{\rm fs}$ (see Methods (Characterization of local density) for the relevance of this estimation of local density). We can see that $P(N_\mathrm{fs})$ has clear bimodality. Furthermore, $P(N_\mathrm{fs})$ and $P(\zeta)$ both can be properly characterized by two Gaussian functions (see Methods) with the same fraction $s$ (see Fig. \[fig:fnfsZeta\]b), following the prediction of the thermodynamic two-state model (Fig. \[fig:fnfsZeta\]c-e and Fig. \[fig:fnfsT\]-\[fig:fnfs2Gau-small\]). This clearly indicates the anti-correlation between $\zeta$ and local density \[see the above feature (1)\]. This result strongly contradicts with the prediction of the continuum models that $P(N_\mathrm{fs})$ should be unimodal Gaussian, which is the case for simple liquids such as Lennard-Jones liquids (Fig. \[fig:fnfs-lj\]). We note that three popular water models all show the bimodal distributions of $P(N_\mathrm{fs})$ (Fig. \[fig:fnfsT\]). We can see that $P(N_\mathrm{fs})$ exhibits a unimodal distribution at very high $T$, but it transforms to a bimodal one (composed of two Gaussians) upon cooling for all the three water models. This clearly indicates the failure of continuum models, and supports the two-state model (see the above Features (1)-(4)).

![Structural bimodality in the coordination number distribution of liquid water. [**a**]{}, A snapshot of liquid TIP4P/2005 water at 1 bar and 240 K, where $s\approx 0.5$, i.e., near the Schottky temperature $T_{s=1/2}$. Two structural motifs, LFTS and DNLS, are highlighted by blue and red colors, respectively. LFTS has four H-bonded nearest neighbors with tetrahedral coordination, whereas DNLS typically has five or more nearest neighbors (typically three of them being H-bonded) with broken tetrahedral symmetry. [**b**]{}, Correlation between the number of water molecules in the spherical first-shell volume of radius of 3.5 Å, $N_\mathrm{fs}$ (i.e., the coordination number), and the structural descriptor $\zeta$ at 1 bar and 240 K. The distributions of $N_\mathrm{fs}$ and $\zeta$ are shown in the right and top sides, respectively. Both distributions show clear bimodal features, which can be properly described by the sum of two Gaussian functions (black lines), with blue and red shades corresponding to LFTS and DNLS respectively. [**c**]{}, Distribution of coordination number $N_\mathrm{fs}$ as a function of temperature, $P \left(N_\mathrm{fs}, T\right)$, is shown by colored balls with red and blue for higher and lower temperatures, respectively. The colored surface is the fit to two Gaussian functions by Eqs. \[eq:pnfs\] - \[eq:SDNLS\]. [**d**]{}-[**e**]{}, $T$-dependence of the two Gaussian components of $P \left(N_\mathrm{fs}, T\right)$ in side ([**d**]{}) and top ([**e**]{}) view. One Gaussian component locating at $N_\mathrm{fs} \approx 4$ corresponds to LFTS, whereas the other one at $N_\mathrm{fs} \simeq 5$ corresponds to DNLS, in agreement with the snapshot in [**a**]{}. The fraction of the two Gaussian components is consistent with the fraction determined by $\zeta$ distribution in [**b**]{} and the theoretical two-state model (Eq. \[eq:s\]).[]{data-label="fig:fnfsZeta"}](fnfsZeta-small){width="14.5cm"}

It has sometimes been argued that the unimodal Gaussian distribution of density fluctuations is a signature against mixture models. However, we point out that it is not the case: Under thermal fluctuations, any thermodynamic order parameters, e.g. density $\rho$ and local structural order $s$, should have unimodal Gaussian distributions. This is because the free energy of a system, $f(\rho,s)$, can be expressed by quadratic terms in the one-phase homogeneous region (see, e.g., Ref. [@tanaka1998simple]). Although theoretically obvious, we can confirm it from the fact that the macroscopic density distributions in liquid water and other single- or two-component liquids commonly show Gaussian distributions (see the results of Lennard-Jones (LJ) liquid, SiO$_2$, and Cu$_{64}$Zr$_{36}$ in Fig. \[fig:fdensity\]) irrespective of whether the local density distribution is unimodal or bimodal (unimodal for LJ liquid whereas bimodal for H$_2$O, SiO$_2$ and Cu$_{64}$Zr$_{36}$). The same is applied for the distribution of another macroscopic order parameter $s$ (estimated from $\zeta$): $P(s)$ has a unimodal Gaussian distribution, even when the underlying microscopic structural descriptor $\zeta$ has a bimodal distribution composed of two Gaussians (Fig. \[fig:fs\_zeta\]). This difference between macroscopic and microscopic distributions clearly indicates that the bimodal structural ordering in liquid water is highly localized, in agreement with Feature (1) in the introduction. Here we note that a mixture model of LDL and HDL should result in the bimodal distributions of macroscopic order parameters $\rho$ and $s$, contrary to the above results.

So far we show that computer simulations of classical water models allow us to directly access the distributions of $\zeta$ and $N_\mathrm{fs}$ and provide strong evidence for the presence of the two types of structural motifs. Unfortunately, however, we cannot access such microscopic molecular-level information by experiments. So an experimentally accessible structural descriptor is highly desirable to close a long-standing debate on the structure of liquid water.

The most powerful experimental methods to access the local structures of materials are x-ray and neutron scatterings, by which we can measure the structure factor, i.e. the density-density correlation in reciprocal space: $$S(\bm{k}) = \frac{1}{N}\langle\rho_{\bm{k}} \rho_{-\bm{k}}\rangle
\label{eq:sk}$$ where $\langle \cdots \rangle$ denotes the ensemble average, $N$ is the number of particles, $\rho_{\bm{k}} = \sum_{i=1}^{N} \exp \left(-i\bm{k}\cdot \bm{r}_{i}\right)$ is the number density, $\bm{r}_i$ is the position vector of particle $i$, and $\bm{k}$ is the wave vector. In a crystal, the density $\rho_{\bm{k}}$ has components only at particular wave vectors $\bm{k}$’s because of the periodic arrangement of particles, leading to sharp diffraction spots at those wave vectors in the structure factor. These spots provide a complete description of a crystal structure. On the other hand, liquids and amorphous solids do not possess long-range translational order, and, as a result, only board isotropic amorphous halos are usually observed, which makes their structural characterization extremely difficult. This has also been the case for liquid water. So far no evidence of the coexistence of two types of structural motifs has been detected in $S(k)$ ($k=|\bm k|$), which has been a main cause of a continuous doubt on the two-state model. In this Letter, however, we report a new analysis of $S(k)$ focusing on the first few peaks, which provides direct experimental evidence for the coexistence of the two types of structural motifs and thus supports the two-state model.

To do so, we focus on the lowest wave number peak in liquid water. In simple liquids such as hard spheres and LJ liquids, the first diffraction peak usually appears at the wave number $k \hat{r}/2\pi \simeq 1$ corresponding to the average nearest neighbor distance $\hat{r}$, or the average interparticle distance. However, it has been reported that $S(k)$ of a wide class of materials has a peak at a lower wave number whose corresponding length is longer than the average nearest neighbor distance [@elliott1991medium]. Such a peak has been widely observed in the so-called tetrahedral liquids such as SiO$_2$, GeO$_2$, BeF$_2$, Si, Ge and C, and widely known as the first sharp diffraction peak (FSDP) [@shi2019distinct]. The emergence of FSDP has been considered as a signature of intermediate-range structural ordering in liquids and amorphous states. Recently we have discovered [@shi2019distinct] that FSDP of these liquids originates from the scattering from the density wave characteristic of a tetrahedral unit in LFTS, which is the fundamental structural motif of tetrahedral materials. More precisely, a density wave whose wave vector corresponds to the height $H$ of the tetrahedral structure (e.g., along the $Z$ direction in Fig. \[fig:sqooTetraZeta\]a) generates a sharp diffraction peak specifically at $k_\mathrm{T1}=k \hat{r}/2\pi \simeq 3/4$, i.e., FSDP. We note that a tetrahedral unit produces four peaks in the range $0.5 \leq k \hat{r}/2\pi \leq 3.0$ with peak wave numbers labelled as $k_{\mathrm{T}i}$ ($i=1\sim4$) from low to high $k$ [@shi2019distinct] (Fig. \[fig:sqooTetraZeta\]b). If the two structural motifs revealed by $\zeta$ for water models are also relevant to real water, there should be the corresponding distinct signatures in the experimentally measured structure factor. Such a signature is indeed seen from the locations of the first diffracton peak in the structure factor of low $T$ and high $P$ water (Fig. \[fig:sqoo\_lowT\_highP\]). We can see a more distinct signature in simulated model waters, for which we are able to access both much lower temperatures (predominantly composed of LFTS), and higher pressures (predominantly composed of DNLS) than for experiments, without suffering from ice crystallisation. Figures \[fig:sqooTetraZeta\]b and c show the partial O-O structure factor of TIP4P/2005 water at low $T$ (LFTS-dominant) and high $P$ (DNLS-dominant), respectively, together with those of typical amorphous tetrahedral materials, C, Si and Ge. We can clearly see that low-$T$ water shows the structure factor very similar to the typical amorphous tetrahedral materials and its FSDP is exactly located at $k_\mathrm{T1}=k \hat{r}/2\pi \simeq 3/4$, as expected [@shi2019distinct]. On the other hand, high-$P$ water has a first diffraction peak at $k_\mathrm{D1}=k \hat{r}/2\pi \simeq 1$, as simple liquids do, reflecting its (partially) disordered nature. In the two-state regime lying between the two extreme conditions, where LFTS and DNLS coexist with comparable populations, distinct signatures from the two structural motifs are expected to appear in the structure factor of liquid water.

To reveal local structural characteristics in the wave-number space, we employ what is called “the Debye scattering function” [@debye1915zerstreuung] (see Eqs. \[eq:skdebye\] - \[eq:skzeta\] and Fig. \[fig:sqdebye\]). This allows us to access the correlation between the local structure of each structural motif characterized by $\zeta$ and its local structure factor on the firm theoretical basis. Figures \[fig:sqooTetraZeta\]d and e show the $\zeta$ dependent O-O structure factor $S(k,\zeta)$ of TIP4P/2005 water at 1 bar and 240 K, where water has equal amount of LFTS and DNLS (or, $s=1/2$). We refer this particular temperature to the Schottky temperature [@shi2018common] and denote it as $T_{s=1/2}$. Strikingly, we can see two distinct peaks at $k_\mathrm{T1}$ and $k_\mathrm{D1}$ in different $\zeta$ domains, which are nicely characterized by the two Gaussian components in the distribution of $\zeta$. Thus, we may conclude that the two peaks at $k_\mathrm{T1}$ and $k_\mathrm{D1}$ in the structure factor should correspond to LFTS and DNLS, respectively (Figs. \[fig:sqooTetraZeta\]e and f). We have also confirmed the same feature for TIP5P and ST2 water (Fig. \[fig:sqooZeta\]). Together with the bimodality of the structural descriptor $\zeta$ and coordination number $N_\mathrm{fs}$, this result further supports the two-state model. We emphasize that our finding indicates that we can now access the two-state signature [*experimentally*]{} by analysing the structure factor of real water.

Unfortunately, however, because the $k_\mathrm{T1}$ and $k_\mathrm{D1}$ peaks are close to each other, they are heavily overlapped under substantial thermal fluctuations, which makes a clear separation difficult. This difficulty is a source of the long-standing controversy. Thanks to the strong two-state nature in liquid silica—a tetrahedral liquid structurally similar to water [@shi2018impact]—and large scattering cross sections of the atoms, we recently found that the apparent “first diffraction peak” in the Si-Si partial structure factor of silica is indeed a doublet: A Lorentzian peaked at $k_\mathrm{T1}$ and a Gaussian peaked at $k_\mathrm{D1}$ are necessary to properly describe the apparent “first diffraction peak”. Moreover, the integrated intensity of the $k_\mathrm{T1}$ component is proportional to the fraction $s$ of LFTS, which is determined independently from a microscopic structural descriptor $z$ [@shi2018impact]; namely, it obeys the prediction of the thermodynamic two-state model [@shi2019distinct].

![Structural bimodality in the structure factor of liquid water. [**a**]{}, A schematic representation of a regular tetrahedron formed by five water molecules with the nearest O-O distance $r_{\rm OO}$, height $H$, width $D$, edge length $L$. [**b**]{}, The O-O partial structure factor $S_\mathrm{AA}(k)$ ($\mathrm{AA}=\mathrm{OO}$) of a regular tetrahedron (formed by five oxygen atoms) and simulated TIP4P/2005 liquid water at 1 bar, 194 K, together with the structure factor of typical amorphous tetrahedral materials, C [@gilkes1995], Si [@laaziri1999] and Ge [@etherington1982]. Deeply supercooled liquid water clearly shows four characteristic peaks ($k_{Ti}$ ($i=1 \sim 4$)) common to tetrahedral materials. The peak position of FSDP, $k_{T1}=kr_\mathrm{OO}/2\pi \simeq 3/4 $, indicated by the arrow cooresponds to the height $H$ of an LFTS [@shi2019distinct]. [**c**]{}, The O-O partial structure factor of TIP5P/2005 water at 10000 bar, 250 K (see Fig. \[fig:sqooP\] for results of real and TIP5P water). It shows a characteristic peak of normal disordered systems at $k_{D1}=kr_\mathrm{OO}/2\pi \simeq 1 $ (see the arrow). [**d**]{},[**e**]{}, The $\zeta$-dependent partial O-O structure factor $S(k,\zeta)$ of TIP4P/2005 water at 1 bar, 240 K, calculated by Debye’s scattering equation (Eqs. \[eq:ski\] - \[eq:skzeta\]), in side ([**d**]{}) and top ([**e**]{}) view (see Fig. \[fig:sqooZeta\] for TIP5P and ST2 models). The characteristic peaks, $k_\mathrm{T1}$ and $k_\mathrm{D1}$, are highlighted by blue and red circles, respectively. [**f**]{}, The distribution of $\zeta$ shows two Gaussian components, corresponding to LFTS (blue shade) and DNLS (red shade) respectively. The wave numbers in [**b**]{}, [**c**]{}, [**d**]{}, [**e**]{} are scaled by the nearest neighbor distance $r_\mathrm{AA}$ for all cases (A = C, Si, Ge and O in water).[]{data-label="fig:sqooTetraZeta"}](sqooTetraZeta-small){width="16cm"}

![Analysis of O-O partial structure factors of real water and model waters at ambient pressure. [**a**]{}-[**d**]{}, X-ray scattering data of real water [@skinner2014structure; @pathak2019intermediate]. [**e**]{}-[**h**]{}, TIP4P/2005 model. [**i**]{}-[**l**]{}, TIP5P model. [**m**]{}-[**p**]{}, ST2 model. The O-O partial structure factors are shown by spheres in 3D plots ([**a**]{}, [**e**]{}, [**i**]{}, [**m**]{}) and circles in 2D plots ([**b**]{}, [**f**]{}, [**j**]{}, [**n**]{}) with more blue and red color for lower and higher temperatures respectively. The colored surfaces in ([**a**]{}, [**e**]{}, [**i**]{}, [**m**]{}) are the fits of our model (two Lorentzian (L1 + L2) and two Gaussian (G1 + G2) functions; see Eq. \[eq:sk2\] for the details) to the structure factors. In [**b**]{}, [**f**]{}, [**j**]{}, [**n**]{}, four characteristic peaks obtained from the fit are displayed by broken lines and assigned as indicated by the arrows in [**b**]{}. The colored surfaces in the right two columns show the temperature dependences of the four characteristic peaks from the fit in side view ([**c**]{}, [**g**]{}, [**k**]{}, [**o**]{}) and top view ([**d**]{}, [**h**]{}, [**l**]{}, [**p**]{}). The color bar is shown in [**c**]{} and the number in it denotes the ratio of the peak height to the maximum height of each peak over the temperature and wave number ranges shown in each image. The wave number is scaled by the nearest neighbor O-O distance $r_\mathrm{OO}$.[]{data-label="fig:sqooT"}](sqooT-small){width="16.6cm"}

Recent progress of x-ray scattering techniques enables to measure structure factors of liquid water with high accuracy down to 254.1 K [@skinner2014structure; @skinner2013benchmark], which makes a detailed structural analysis possible even for real water, as for silica. Here we analyse the O-O partial structure factors of real water as well as TIP4P/2005, TIP5P, and ST2 waters by using four peak functions for fitting (see Fig. \[fig:sqooT\]). Indeed, we find that the apparent “first diffraction peak” in the O-O partial structure factors of real water as well as model waters can be nicely described by the sum of a Lorentzian (L1) and a Guassian (G1) over a wide temperature range (Figs. \[fig:sqooT\] and \[fig:sqoofit\]). We call this fitting scheme ‘scheme II’ (see Methods). In particular, the Lorentzian and Gaussian functions have peaks at $k_\mathrm{T1}=kr_{\mathrm{OO}}/2\pi \sim 3/4$ and $k_\mathrm{D1}=kr_{\mathrm{OO}}/2\pi \sim 1$, corresponding to LFTS and DNLS respectively, in agreement with the above-mentioned silica case and the Debye scattering function shown in Fig. \[fig:sqooTetraZeta\]. The integrated intensity of the Lorentzian peak follows the prediction of the two-state model and agrees well with the fraction of LFTS, $s$, determined independently by $\zeta$ and $N_\mathrm{fs}$. Here we note that the Lorentzian and Gaussian shapes reflect the different nature of the two structural motifs: LFTS has rather unique local tetrahedral order, whereas DNLS intrinsically has high structural fluctuations. The presence of the bimodality in the experimental structure factor of liquid water (Fig. \[fig:sqooT\]), as well as in $\zeta$ [@Russo2014; @shi2018Microscopic; @shi2018impact; @shi2018common] and $N_\mathrm{fs}$, together with their inter-consistency, unambiguously show the existence of the two types of local structures (LFTS and DNLS) in liquid water and thus support the two-state description of liquid water.

                     Real water   TIP4P/2005    TIP5P      ST2
  ----------------- ------------ ------------ --------- ---------
   $\Delta E$ (K)     -1929.0      -1802.0     -3355.9   -4612.5
   $\Delta \sigma$    -8.2845      -7.5779     -13.134   -16.106
     $T_{s=1/2}$       232.85       237.80     255.51    286.39

  : Two-state parameters for real water and model waters.

\[Tablepara\]

Unlike at low $T$, we find that at high $T$ only one Gaussian function is enough to properly describe the apparent “first diffraction peak” in the experimental and simulated O-O structure factors. We call this fitting scheme ‘Scheme I’. One might think that Scheme I might work even at any temperatures, which is expected for continuum models. Thus, to rationalise the relevance of Scheme II at low $T$, or to confirm the bimodality of the apparently first diffraction peak in an unambiguous manner, we show in Fig. \[fig:fitError\] the difference in the mean squared residual, which measures the deviation of the fit from the data, between Schemes I and II as a function of the scaled temperature $T/T_{s=1/2}$. We can clearly see a tendency common to both real water and simulated model waters: At temperatures above 1.1$T_{s=1/2}$ a single Gaussian (Scheme I) can describe the apparent “first diffraction peak” in the structure factor equally well as a Gaussian plus a Lorentzian function (Scheme II). Below 1.1$T_{s=1/2}$, on the other hand, Scheme I starts to seriously fail in describing the data, reflecting the rapid growth of the fraction of LFTS below that temperature. The failure of Scheme I at low temperatures not only supports the emergence of the bimodality in the apparent “first diffraction peak” there, but also explain why the two-state feature can hardly be observed in liquid water at ambient condition [@smith2005unified; @clark2010small; @niskanen2019compatibility].

Moreover, our two-state description (Scheme II) of the structure factor provides a direct experimental access to the degree and range of local structural ordering in real water. The fraction $s$ of LFTS, which is propotional to the integrated intensity of FSDP at $k_\mathrm{T1}$, increases rapidly towards the LDL limit ($s\simeq 1$) upon cooling, as shown in Fig. \[fig:sl\]a. The increase is faster for TIP5P and ST2 water than for TIP4P/2005 and the real water, indicating the “over-structured" tendency in the former two models. In the two-state language, TIP5P and ST2 models overestimate the energy gain and entropy loss upon the formation of LFTS, as shown by the two-state-model parameters in Table \[Tablepara\].

Figure \[fig:sl\]b shows the increase of the coherence length estimated from the width of FSDP (Eq. \[eq:L\]) upon cooling. Below $T_{s=1/2}$, the coherence lengths estimated from the experimental and simulated structure factors increase and converge towards the $s \rightarrow 1$ limit upon cooling. Above $T_{s=1/2}$, on the other hand, the fraction of LFTS, i.e. the integrated intensity of FSDP, is rather small and thus the data suffer from large uncertainty. In any case, the coherence length is very short, bounded between $\sim2$ Å of a single tetrahedron and $\sim6.5$ Å of LDA ice (see Fig. \[fig:sqooTetraLDA\] for the detail), in agreement with the previous measurements of structural correlation length [@xie1993noncritical; @huang2009inhomogeneous; @kim2017maxima] in real water and dynamic correlation length in TIP5P water [@shi2018origin; @shi2018common] \[Feature (1) in the introduction\].

![Degree and range of local tetrahedral ordering in liquid water at ambient pressure. [**a**]{}, The integrated intensity of FSDP ($k_\mathrm{T1}$ peak in Fig. \[fig:sqooT\]b) as a measure of the degree of local tetrahedral ordering, or the fraction $s$ of LFTS, which monotonically increases upon cooling. The horizantal dash line indicates the upper limit of tetrahedral ordering in liquid water. [**b**]{}, The coherence length of FSDP as a measure of the range of local tetrahedral ordering, which monotonically increases with decreasing temperature. The high and low temperature limits of the coherence length from a single tetrahedron and LDA ice [@mariedahl2018] are shown by violet dash dot line and navy dot line, respectively (see Fig. \[fig:sqooTetraLDA\]). The correlation lengths determined by the Ornstein-Zernike analysis of small-angle X-ray scattering data (typically $k<0.5$ Å$^{-1}$) by different groups [@xie1993noncritical; @huang2009inhomogeneous; @kim2017maxima] are shown by dotted lines. The maximum correlation length of $\xi\simeq 4.1$ Å at $T_{s=1/2} = 229.2$ K estimated from recent small-angle X-ray scattering measurements of liquid water droplets [@kim2017maxima] is shown by the magenta star symbol. []{data-label="fig:sl"}](sl){width="14cm"}

We show the first clear experimentally accessible evidence in the structure factor for the dynamical coexistence of the two types of structural motifs, LFTS and DNLS, supporting the two-state description of liquid water. We reveal that liquid water exhibits the so-called FSDP in the structure factor as other tetrahedral liquids do. The FSDP provides crucial information on the fraction of LFTS (degree of structural ordering, i.e., the order parameter of the two-state model (Eq. \[eq:s\])) and its coherence length (range of structural ordering). We hope that these findings will contribute to the convergence of long-standing debates on the structure of water.

[**Acknowledgements**]{} The authors are grateful to R. Evans for his valuable suggestion on the analysis regarding the criticality. This study was partly supported by Scientific Research (A) and Specially Promoted Research (KAKENHI Grants No. JP18H03675 and No. JP25000002 respectively) from the Japan Society for the Promotion of Science (JSPS) and the Mitsubishi Foundation.

Methods {#methods .unnumbered}
=======

Simulation of water {#simulation-of-water .unnumbered}
-------------------

Classical molecular dynamics simulations were performed in a periodic cubic box containing 1000 TIP4P/2005 [@abascal2005general] water molecules by using the Gromacs package [@Hess2008] with a time step of 2 fs. Intermolecular van der Waals forces and Coulomb interactions in real space were truncated at 9 Å, and long-range Coulomb interactions were treated by the particle-mesh Ewald method. Long-range dispersion corrections for energy and pressure were applied. All simulations were performed in *NPT* ensemble with temperature and pressure kept constant by Nos[é]{}-Hoover thermostat and Parrinello-Rahman barostat, respectively. All the bonds are constrained by using the LINCS algorithm. Long-time simulations (typically longer than 100 times molecular reorientation time) were performed after equilibration in a wide temperature range from 194 to 300 K at 1 bar, and in a wide pressure range from 1 to 10000 bar at 250 K. At 194, 197 and 200 K, two independent trajectories were generated to enhance the statistics. The simulation times for production runs are summarized in Tables \[TabletimeT\] and \[TabletimeP\]. The simulation details of TIP5P and ST2 model can be found in Refs. [@shi2018origin; @shi2018common]. Ice nucleation has not been observed at any temperature and pressure studied in this work for all the water models.

   $T$ (K)               194                         197                         200              210   220   230   240
  ---------- --------------------------- --------------------------- --------------------------- ----- ----- ----- -----
   $t$ (ns)   $\mathbf{2} \times 38000$   $\mathbf{2} \times 41000$   $\mathbf{2} \times 23000$   600   40    10     5
   $T$ (K)               250                         260                         270              280   290   300  
   $t$ (ns)               3                           2                          1.2              1.2   1.0   1.0  

  : Simulation times $t$ used for production runs for TIP4P/2005 water at 1 bar. The bold multiplers indicate the number of independent runs.

\[TabletimeT\]

   $P$ (bar)   1   1000   1800   4000   6000   8000   10000
  ----------- --- ------ ------ ------ ------ ------ -------
   $t$ (ns)    3    2     2.6     20     40     80     100

  : Simulation times $t$ used for production runs for TIP4P/2005 water at 250 K.

\[TabletimeP\]

Simulation of silica {#simulation-of-silica .unnumbered}
--------------------

A two-component system of 3456 BKS [@Beest1990; @Saika2000] silica (3456 silicon and 6912 oxygen ions) in a periodic cubic box was simulated with a time step of 0.5 fs by using the LAMMPS package [@Plimpton1995]. Short range interactions were truncated at 5.5 Å, and long-range electrostatic interactions were treated by using the particle-particle particle-mesh method. All simulations were performed in *NPT* ensemble with temperature and pressure controlled by Nos[é]{}-Hoover thermostat and barostat, respectively. Production runs were carried out at ambient pressure for 200 ps each at 5500 K and 6000 K, after equilibration runs of the same length.

Simulation of metallic glass Cu$_{64}$Zr$_{36}$ {#simulation-of-metallic-glass-cu_64zr_36 .unnumbered}
-----------------------------------------------

A binary metallic glass Cu$_{64}$Zr$_{36}$, containing 6400 copper and 3600 zirconium atoms was simulated by using the embedded atom method potential [@mendelev2009]. Molecular dynamics simulations were carried out with a time step of 1.0 fs by using the LAMMPS package [@Plimpton1995]. Periodic boundary condition was applied to all directions of the cubic box. Nos[é]{}-Hoover thermostat and barostat were employed to keep the temperature at 1800 K and pressure at 1 bar, respectively. A production run of 400 ps was performed at 1 bar, 1800 K after a 300 ps equilibration run.

Simulation of Lennard-Jones (LJ) liquid {#simulation-of-lennard-jones-lj-liquid .unnumbered}
---------------------------------------

A one-component LJ system of 6912 particle interacting via the standard 12-6 LJ potential was simulated with a time step of 0.005 by using the LAMMPS package [@Plimpton1995]. The interaction was truncated and force-shifted at $2.5\sigma$, so that both the potential and force smoothly go to zero at $2.5\sigma$. *NVT* and *NPT* simulations were performed at $\rho=0.7$, $t=0.75$ and $p=0.05$, $t=0.75$ (in reduced unit) for 4000000 steps, respectively.

Characterization of local density {#characterization-of-local-density .unnumbered}
---------------------------------

In this work the local density is characterized by the coordination number that measures the number of neighboring water molecules in the first coordination shell of a center molecule. The first coordination shell is defined as a sphere with a radius corresponding to the position of first minimum in the oxygen-oxygen radial distribution function, which is typically $\sim3.5$ Å, in comparable with the coherence length for liquid water (Fig. 4b).

Other approaches such as Voronoi tessellation [@duboue2015characterization] and density in grids [@soper2010recent; @english2011density] have also been used to measure the local density of liquid water. The former estimates local density by using the Voronoi volume of each molecule, whereas the latter calculates the number of molecules in a small cubic box with different sizes (typically $>9$ Å). Both of these two methods report unimodal density distributions, which have often been taken as direct evidence against the two-state model. However, we argue that neither of the two methods is a proper measure of the local density. For the Voronoi tessellation method, it has been shown that its application to tetrahedral materials such as amorphous silicon suffers from a serious deficiency because of the low coordination number [@tsumuraya1993statistics]. For the grid method, on the other hand, because of the small coherence length of the structural motifs (Fig. 4b), a box of $9 \times 9 \times 9$ Å is too large to detect the local density fluctuation associated with them: For a box significantly larger than the size of the local structural motifs, the density distribution is expected to be unimodal and Gaussian, as shown in Fig. \[fig:fdensity\].

Fitting formula for the coordination number distribution {#fitting-formula-for-the-coordination-number-distribution .unnumbered}
--------------------------------------------------------

At high enough temperatures, the distribution of coordination number, $P\left(N_\mathrm{fs}\right)$, of liquid water shows a unimodal Gaussian distribution as simple LJ liquid (see Fig. \[fig:fnfs-lj\] and Fig. \[fig:fnfsT\]a,d,g). However, $P\left(N_\mathrm{fs}\right)$ of liquid water significantly deviates from a single Gaussian distribution and instead displays a bimodal distribution upon cooling, which strongly suggests the development of two structural motifs in supercooled water (see Fig. 1c and Fig. \[fig:fnfs2Gau-small\]). As a result, we find that $P\left(N_\mathrm{fs}\right)$ can be properly described by the sum of two Gaussian functions, whose integrated intensity corresponds to the fraction of LFTS and DNLS respectively, $$\begin{split}
P\left(N_\mathrm{fs}\right) & = \frac{s}{\sigma_\mathrm{LFTS} \sqrt{2\pi}} \exp \left[-\frac{\left(N_\mathrm{fs}-N_\mathrm{LFTS}\right)^2}{2\sigma_\mathrm{LFTS}^2}\right] \\ & + \frac{1-s}{\sigma_\mathrm{DNLS} \sqrt{2\pi}} \exp \left[-\frac{\left(N_\mathrm{fs}-N_\mathrm{DNLS}\right)^2}{2\sigma_\mathrm{DNLS}^2}\right]
\label{eq:pnfs}
\end{split}$$ In this equation, $s$ is defined by Eq. \[eq:s\] and other parameters can be described as follows, $$N_\mathrm{LFTS} = a_{00} + a_{01}T + a_{02}T^2
\label{eq:NLFTS}$$ $$N_\mathrm{DNLS} = a_{10} + a_{11}T + a_{12}T^2
\label{eq:NDNLS}$$ $$\sigma_\mathrm{LFTS} = b_{00} + b_{01}T + b_{02}T^2
\label{eq:SLFTS}$$ $$\sigma_\mathrm{DNLS} = b_{10} + b_{11}T + b_{12}T^2
\label{eq:SDNLS}$$

Calculation of the structure factor {#calculation-of-the-structure-factor .unnumbered}
-----------------------------------

The structure factor $S(\bm{k})$, defined as the density-density correlation in reciprocal space, can be obtained by $$S(\bm{k}) = \frac{1}{N}\langle\rho_{\bm{k}} \rho_{-\bm{k}}\rangle
\label{eq:skS}$$ where $\langle \cdots \rangle$ represents the ensemble average, $\bm{k}$ is the wave vector, $N$ is the number of particles and $\rho_{\bm{k}}$ is the Fourier component of the number density $\rho$, which is given by $$\rho_{\bm{k}} = \stackrel[i=1]{N}{\sum}\exp \left(-i\bm{k}\cdot \bm{r}_{i}\right)
\label{eq:rhok}$$ where $\bm{r}_i$ is the coordinates of atom $i$. For an isotropic system, the structure factor is a function of only the magnitude of the wave vector, $k$: $S(k)$.

Debye scattering function {#debye-scattering-function .unnumbered}
-------------------------

The structure factor can also be calculated by the Debye scattering function [@debye1915zerstreuung]: $$S(k) = 1 + \frac{1}{N}  \stackrel[i=1]{N}{\sum}  \stackrel[j\neq i]{N}{\sum} \frac{\sin (kr_{ij})}{kr_{ij}} W(r_{ij})
\label{eq:skdebye}$$ where $W(r_{ij}) = \frac{\sin (\pi r_{ij}/r_{c})}{\pi r_{ij}/r_{c}}$ is the window function [@zhang2019improving] and $r_{c}$ is a cutoff distance.

Debye scattering function allows for a local structural characterization by the molecular structure factor: $$S_{i}(k) = 1 + \stackrel[j\neq i]{N}{\sum} \frac{\sin (kr_{ij})}{kr_{ij}} W(r_{ij})
\label{eq:ski}$$

Then, the correlation between molecular structure factor $S_{i}(k)$ and local structure descriptor $\zeta$ can be evaluated by the $\zeta$-dependent structure factor: $$S(k,\zeta) = \frac{\stackrel[i]{N}{\sum} S_{i}(k)\delta(\zeta-\zeta(i))}{\stackrel[i]{N}{\sum} \delta(\zeta-\zeta(i))}
\label{eq:skzeta}$$

Fitting formula for the structure factor {#fitting-formula-for-the-structure-factor .unnumbered}
----------------------------------------

The first three peaks, the first of which is actually a doublet, in the structure factor of liquid water can be well described by two Lorentzian and two Gaussian functions as $$\begin{aligned}
S\left(k\right) = \frac{f_\mathrm{T1}}{\pi}\frac{\Gamma_\mathrm{T1}}{\left(k-k_\mathrm{T1}\right)^2 + \Gamma_\mathrm{T1}^2} + \frac{f_\mathrm{D1}}{\sigma_\mathrm{D1} \sqrt{2\pi}} \exp \left[-\frac{\left(k-k_\mathrm{D1}\right)^2}{2\sigma_\mathrm{D1}^2}\right] \\+ \frac{f_\mathrm{T2}}{\pi}\frac{\Gamma_\mathrm{T2}}{\left(k-k_\mathrm{T2}\right)^2 + \Gamma_\mathrm{T2}^2} + \frac{f_\mathrm{T3}}{\sigma_\mathrm{T3} \sqrt{2\pi}} \exp \left[-\frac{\left(k-k_\mathrm{T3}\right)^2}{2\sigma_\mathrm{T3}^2}\right]
\label{eq:sk2}
\end{aligned}$$ where the subscripts denote the peaks in the O-O partial structure factor as shown in Fig. 3b. In this equation, all the parameters depend on temperature and pressure, and therefore a large set of parameters are needed to describe the structure factor of liquid water at different thermodynamic conditions (12 parameters for each temperature and pressure).

However, thanks to the weak temperature dependence of the fitting parameters (except for $f_\mathrm{T1}$ and $f_\mathrm{D1}$), we found that they can be well described by a set of polynomial functions up to the second order: $$k_{x} (T) = \tilde{k}_{x0} + \tilde{k}_{x1} \hat{T} + \tilde{k}_{x2} \hat{T}^2
\label{eq:polyk}$$ $$\Gamma_{x} (T) = \tilde{\Gamma}_{x0} + \tilde{\Gamma}_{x1} \hat{T} + \tilde{\Gamma}_{x2} \hat{T}^2
\label{eq:polyg}$$ $$\sigma_{x} (T) = \tilde{\sigma}_{x0} + \tilde{\sigma}_{x1} \hat{T} + \tilde{\sigma}_{x2} \hat{T}^2
\label{eq:polys}$$ $$f_{x} (T) = \tilde{f}_{x0} + \tilde{f}_{x1} \hat{T} + \tilde{f}_{x2} \hat{T}^2
\label{eq:polyf}$$ where the subscript $x$ ($x = \textrm{T1, D1, T2 or T3}$) denotes parameters for each peak and $\hat{T}=T/T_\mathrm{ref}$ with $T_\mathrm{ref} = 373.15$ K. This procedure allows for a simultaneous fitting of a large set of structure factors measured in a wide temperature range, which largely reduces the number of fitting parameters. Moreover, we found in practice that the fitting accuracy will not be affected if we set $\tilde{k}_\mathrm{x2}=0$ (for $x = \textrm{T1, D1, T2 or T3}$), and fix the value of $k_\mathrm{T3}$ and $\tilde{k}_\mathrm{T11}/\tilde{k}_\mathrm{D11}$ properly. Although the $k_\mathrm{T3}$ peak is only partially included in the fitting, its position is read from the data and fixed in the fitting precedure.

On the other hand, the intensities of T1 and D1 peaks vary significantly with temperature and pressure, corresponding to the change in the fractions of the two structural motifs. In our previous study [@shi2019distinct], we found that the integrated intensity $f_{T1}$ of FSDP is proportional to the fraction $s$ of LFTS in liquid silica as $$f_\mathrm{T1} (T,P) = a\cdot s
\label{eq:ft1}$$ where $a$ is a positive constant. This knowledge is directly applied to liquid water, since they are both characterized by the same two-state features [@shi2019distinct].

Here $s$ can further be described by the two-state model with negligible cooperativity as [@Tanaka_review; @Tanaka2000; @russo2018water; @shi2018impact; @shi2018Microscopic; @shi2018origin; @shi2018common]: $$s(T,P) = \frac{1}{1+\exp \left(\frac{\Delta E - T \Delta \sigma + P \Delta V}{k_\mathrm{B}T}\right)}
\label{eq:s}$$ where $k_\mathrm{B}$ is the Boltzmann constant, $\Delta E$, $\Delta \sigma$ and $\Delta V$ are the energy, entropy, and volume differences between LFTS and DNLS in the two-state model. At ambient pressure, the term $P\Delta V$ is negligibly small. The parameters $\Delta E$ and $\Delta \sigma$ for TIP5P and ST2 waters have already been determined in our previous work [@shi2018common]. For TIP4P/2005 and real water, $\Delta E$ and $\Delta \sigma$ were independently determined by applying the two-state model (Eq. \[eq:s\]) to the fraction $s(T)$ of LFTS that can be obtained from the $g_\mathrm{OO}(r)$ by $s(T)=1-g_\mathrm{OO}(r=r_\mathrm{HB})$ [@Russo2014], where $r_\mathrm{HB}=3.5$ Å, according the Luzar-Chandler definition of H-bond [@Luzar1996hydrogen]. The parameters $\Delta E$ and $\Delta \sigma$ for real water, TIP4P/2005, TIP5P and ST2 waters are summarised in Table 1 in the main text.

Since $D1$ peak is exclusively from DNLS, whose fraction is given by $(1-s)$, we can formulate the temperature and pressure dependence of $f_{D1}$ as $$f_{D1} (T,P) = b(1-s)=b \left[1- \frac{1}{1+\exp \left( \frac{\Delta E - T \Delta \sigma + P \Delta V}{k_\mathrm{B}T}\right)}\right]
\label{eq:fd1}$$ where $b$ is a positive constant. After all the above considerations, only 25 free fitting parameters are necessary to fit O-O partial structure factors of liquid water at all the temperatures studied in this work.

The Gaussian function represents the scattering peak coming from the interatomic correlation of DNLS (see main text). Simple liquids such as LJ and hard spheres liquids usually have this peak at the wave number corresponding to the neighboring interactomic distance $r$, and thus we constrained its position $k_{0}r_\mathrm{OO}/2\pi$ to be close to 1. The Fourier transform of a Lorentzian function is an exponentially decaying function in real space. The coherence length $\lambda$ of FSDP, which characterizes the range of coherent tetrahedral ordering, can be estimated by $$\lambda=1/\Gamma_\mathrm{T1}
\label{eq:L}$$ where $\Gamma_\mathrm{T1}$ is the half width of FSDP.

![An example of the typical distribution of coordination number $N_\mathrm{fs}$, $P\left(N_\mathrm{fs}\right)$, for a simple liquid. Here we show $P\left(N_\mathrm{fs}\right)$ of Lennard-Jones liquid at $\rho = 0.7$ and $T=0.75$. $P\left(N_\mathrm{fs}\right)$ can be well described by a Gaussian distribution as indicated by the red curve.[]{data-label="fig:fnfs-lj"}](fnfs-lj){width="12cm"}

![The distribution of coordination number $N_\mathrm{fs}$, $P\left(N_\mathrm{fs}\right)$, for simulated water. Here we show $P\left(N_\mathrm{fs}\right)$ for TIP4P/2005 water at 380 K ([**a**]{}), 280 K ([**b**]{}), 240 K ([**c**]{}), for TIP5P water at 360 K ([**d**]{}), 280 K ([**e**]{}), 260 K ([**f**]{}), and for ST2 water at 360 K ([**g**]{}), 300 K ([**h**]{}), 285 K ([**i**]{}), at ambient pressure. At high temperatures, $P\left(N_\mathrm{fs}\right)$ shows a Gaussian distribution as simple Lennard-Jones liquid does (Fig. \[fig:fnfs-lj\]), whereas at low temperatures it changes to bimodal distributions, which can be properly described by the sum of two Gaussian functions (with blue and red shades). The fraction of the two Gaussian components agrees well with the fraction independently determined by $\zeta$ distribution and the prediction of the theoretical two-state model (Eq. \[eq:s\]). One component (with blue shade) corresponds to LFTS, in which the central water typically has $\sim 4$ H-bonded nearest neighbors, whereas the other (with red shade) corresponds to DNLS, in which the central water has $\sim 5$ nearest neighbors (three of them being H-bonded typically) on average at ambient pressure.[]{data-label="fig:fnfsT"}](fnfsT){width="16cm"}

![$T$-dependence of the distribution of coordination number $N_\mathrm{fs}$, $P\left(N_\mathrm{fs}, T\right)$, for three water models. [**a**]{}, [**d**]{}, [**g**]{}, $P \left(N_\mathrm{fs}, T\right)$ for simulated liquid TIP4P/2005 ([**a**]{}), TIP5P ([**d**]{}) and ST2 ([**g**]{}) water at ambient pressure. The colored surfaces in [**a**]{}, [**d**]{}, [**g**]{} are the fits to two Gaussian functions by Eqs. \[eq:pnfs\] - \[eq:SDNLS\]. The two Gaussian components of $P \left(N_\mathrm{fs}, T\right)$ are shown in side ([**b**]{}, [**e**]{}, [**h**]{}) and top ([**c**]{}, [**f**]{}, [**i**]{}) view for TIP4P/2005 ([**b**]{}, [**c**]{}), TIP5P ([**e**]{}, [**f**]{}) and ST2 ([**h**]{}, [**i**]{}) water.[]{data-label="fig:fnfs2Gau-small"}](fnfs2Gau-small){width="16cm"}

![Macroscopic density distribution for single- and two-component systems. Macroscopic density distribution $P(\rho)$ for two-component BKS silica at $T=6000$ K, $P=1$ bar ([**a**]{}), two-component metallic glass Cu$_{64}$Zr$_{36}$ at $T=1800$ K, $P=1$ bar ([**b**]{}), single-component Lennard-Jones liquid at $T=0.75$, $P=0.05$ ([**c**]{}) and single-component TIP4P/2005 water at $T=240$ K, $P=1$ bar ([**d**]{}), TIP5P water at $T=260$ K, $P=1$ bar ([**e**]{}) and ST2 water at $T=285$ K, $P=1$ bar ([**f**]{}). The macroscopic density distributions in panels [**d**]{}, [**e**]{} and [**f**]{} were measured at the same conditions as in panels [**c**]{}, [**f**]{} and [**i**]{} in Fig. \[fig:fnfsT\] for TIP4P/2005, TIP5P and ST2 water, respectively. The macroscopic density distribution always remains unimodal and Gaussian for both single- and two-component liquids, under thermal fluctuations, as it should be.[]{data-label="fig:fdensity"}](fdensity){width="14cm"}

![Distribution of local structural descriptor $\zeta$ and macroscopic order parameter $s$. [**a**]{}-[**c**]{}, Distribution of local structural descriptor $\zeta$, $P(\zeta)$, for TIP4P/2005 ([**a**]{}), TIP5P ([**b**]{}) and ST2 water ([**c**]{}) at ambient pressure. [**d**]{}-[**f**]{}, Distribution of macroscopic structural order parameter $s$, $P(s)$, for TIP4P/2005 ([**d**]{}), TIP5P ([**e**]{}) and ST2 water ([**f**]{}) at ambient pressure. All the distributions are shown at $T\approx T_{s=1/2}$ (240 K for TIP4P/2005 water, 260 K for TIP5P water and 285 K for ST2). The distribution of the macroscopic order parameter, $P(s)$, always remains unimodal and Gaussian for all the models, whereas the distribution of local order, $P(\zeta)$, is clearly bimodal with two Gaussian components (blue and red shades correspond to LFTS and DNLS, respectively). The macroscopic order parameter $s$ is estimated as the fraction of molecules whose $\zeta > \zeta_c$ at each time frame. A threshold value $\zeta_c$ ($\simeq 0.5$ Å) is chosen to satisfy $s\approx 0.5$ after time averaging.[]{data-label="fig:fs_zeta"}](fs_zeta){width="14cm"}

![Partial O-O structure factors of liquid water measured by x-ray scattering experiments at 0.1 MPa, 254.1 K [@skinner2014structure] (blue circle) and 362 MPa, 300 K [@skinner2016structure] (red square). Distinct characteristic peaks $k_\mathrm{T1}$ and $k_\mathrm{D1}$ corresponding to LFTS and DNLS respectively, are clearly located at different wave numbers for these two conditions, whereas the other peaks are located at almost the same positions.[]{data-label="fig:sqoo_lowT_highP"}](sqoo_lowT_highP){width="12cm"}

![The validity of the Debye’s scattering function analysis. The parital O-O structure factors obtained by the two-body density correlation (Eq. \[eq:skS\], black circle) and Debye scattering equation (Eq. \[eq:skdebye\], red curve) for TIP4P/2005 ([**a**]{}), TIP5P ([**b**]{}) and ST2 ([**c**]{}) water at 240 K, 260 K and 285 K, respectively, at ambient pressure. The structure factors calculated from the two methods agree well with each other. A small deviation at the second peak $kr_\mathrm{OO}/2\pi\approx 1.3$ mainly comes from the effect of the Window function we employed (see Methods).[]{data-label="fig:sqdebye"}](sqdebye){width="10cm"}

![The $\zeta$ dependent partial O-O structure factor $S(k,\zeta)$ obtained by Debye’s scattering equation (Eqs. \[eq:ski\] and \[eq:skzeta\]). Side ([**a**]{}, [**d**]{}, [**g**]{}) and top ([**b**]{}, [**e**]{}, [**h**]{}) views for TIP4P/2005 ([**a**]{}, [**b**]{}) TIP5P ([**d**]{}, [**e**]{}) and ST2 ([**g**]{}, [**h**]{}) water at 240 K, 260 K and 285 K, respectively, at ambient pressure. In ([**b**]{}, [**e**]{}, [**h**]{}) the characteristic peaks, $k_\mathrm{T1}$ and $k_\mathrm{D1}$, are highlighted by blue and red circles, respectively. [**c**]{},[**f**]{},[**i**]{}, The distribution of $\zeta$ shows two Gaussian components, corresponding to LFTS (blue shade) and DNLS (red shade) respectively for TIP4P/2005 ([**c**]{}), TIP5P ([**f**]{}) and ST2 ([**i**]{}) water. The wave numbers are scaled by the nearest neighbor O-O distance $r_\mathrm{OO}$.[]{data-label="fig:sqooZeta"}](sqooZeta-small){width="14cm"}

![Partial structure factor $S_\mathrm{AA}(k)$ of simulated water (A=O) and silica (A=Si) at $T\approx T_{s=1/2}$. At ambient pressure, $T_{s=1/2} \approx $ 240, 260, 285 and 5000 K for TIP4P/2005, TIP5P, ST2 water and BKS silica, respectively. Silica shows very similar anomalous behaviors as water does, although no critical point has been found in silica [@lascaris2014search]. Recently strong evidence supporting the existence of two structural motifs—LFTS and DNLS, in liquid silica has been revealed for BKS silica [@shi2018impact; @shi2019distinct]. The similarity between the structure factor of liquid water and liquid silica at the same fraction of LFTS further supports the two-state feature in both liquids.[]{data-label="fig:sqWaterSilicaTw"}](sqWaterSilicaTw){width="12cm"}

![Experimental O-O partial structure factor of water at high pressure. [**a**]{}, X-ray scattering data [@skinner2016structure]. [**b**]{}, TIP4P/2005 model. [**c**]{}, TIP5P model. The arrow denotes the direction of pressure increase. The wave number is scaled by the nearest neighbor O-O distance $r_\mathrm{OO}$.[]{data-label="fig:sqooP"}](sqooP){width="10cm"}

![Decomposition of the O-O partial structure factor of real water and model waters at ambient pressure. [**a**]{}-[**c**]{}, X-ray scattering data [@skinner2014structure; @skinner2013benchmark; @pathak2019intermediate]. [**d**]{}-[**f**]{}, TIP4P/2005 water. [**g**]{}-[**i**]{}, TIP5P water. [**j**]{}-[**l**]{}, ST2 water. The O-O partial structure factors are shown by black squares. The red lines are the cumulative fits of two Lorentzian (L1 + L2) and two Gaussian (G1 + G2) functions to the structure factors. The four characteristic peaks obtained from the fits are displayed by green, blue, orange and magenta curves. The FSDP is highlighted by green shading. The wave number is scaled by the nearest neighbor O-O distance $r_\mathrm{OO}$.[]{data-label="fig:sqoofit"}](sqoofit){width="16cm"}

![The fitting quality of structure factors. Two schemes were applied to the analysis of the partial O-O structure factor $S_\mathrm{OO}(k)$ of liquid water. In Scheme I, the “apparent” first three peaks (from low to high wave number) in $S_\mathrm{OO}(k)$ are fitted to one Gaussian (1st peak), one Lorentizan (2nd peak), and one Gaussian (3rd peak) function, respectively. All the fitting parameters are free, and thus we need a new set of 9 parameters for each temperature. In Scheme II, the “apparent” first three peaks (from low to high wave number) in $S_\mathrm{OO}(k)$ are fitted to one Lorentzian + one Gaussian (1st peak), one Lorentizan (2nd peak), and one Gaussian (3rd peak) function, respectively. We constrain the temperature dependence of the fitting parameters. We need only up to 25 free fitting parameters for simultaneous fitting of $S_\mathrm{OO}(k)$ at all the temperatures (see Methods). Scheme I is chosen for comparison, since it reasonably well describes $S_\mathrm{OO}(k)$ of liquid water at high temperatures with only three functions (8 fitting parameters at each temperature with position of $k_\mathrm{T3}$ peak being fixed). The mean squared residual $\chi^2$ measures the deviation of the fits from the data. The difference in the mean squared residual between Schemes I and II, $\Delta \chi^2 = \chi^2_\mathrm{I} - \chi^2_\mathrm{II}$, is shown for real water, simulated TIP4P/2005, TIP5P and ST2 waters as a function of scaled temperature $T/T_{s=1/2}$. All the curves collapse into a master curve, suggesting a universal behavior independent of models. For all the systems, $\Delta \chi^2$ remains almost zero above $1.1T_{s=1/2}$, whereas sharply increases below $1.1T_{s=1/2}$, reflecting the fast growth of LFTS upon cooling (Fig. 4a). The failure of Scheme I at lower temperatures strongly suggests that the structure of liquid water fundamentally changes with temperature. At low enough temperatures (below $\sim 1.1T_{s=1/2}$), liquid water can no longer be treated to be microscopically homogeneous unlike at high temperatures, and should be regarded as a dynamical mixture of two local structure motifs. The two structure motifs have different contributions to the “apparent” first diffraction peak in $S_\mathrm{OO}(k)$ (Figs. 2-3), leading to the failure of Scheme I that assumes only a single type of contribution to the peak. For real water, there are three data points significantly deviating from the master curve, which are attributed to the experimental errors in the structure factor data (see Figure \[fig:noise\])[]{data-label="fig:fitError"}](fitError){width="12cm"}

![Noise analysis in the experimental structural factor of supercooled water. [**a**]{}, The height of the 1st, 2nd and 3rd peak in the experimental partial O-O structure factor as a function of temperature. The height of the 3rd peak is multiplied by a factor of 1.16 for clarification. [**b**]{}, The height of the 1st and 2nd trough in the experimental partial O-O structure factor as a function of temperature. In [**a**]{} and [**b**]{} the filled and open symbols are data from Ref. [@skinner2014structure] and [@pathak2019intermediate], respectively. [**c**]{}, The experimental partial O-O structure factor of supercooled water at 234.8, 240.9 and 243.8 K (Ref. [@pathak2019intermediate]). [**d**]{}, The experimental partial O-O structure factor of supercooled water at 254.1, 263.1, 264.0, 268.1 and 277.1 K (solid line from Ref. [@skinner2014structure] and dashed line from Ref. [@pathak2019intermediate]). The inset shows the enlarged plot of the 1st peak inside the black square. Clearly, the structure factors measured at 234.8 K, 264.0 K and 268.1 K suffer from large uncertainty, because they significantly deviate from the normal temperature dependence, as highlighted by the black and red circles in [**a**]{} and [**b**]{}. The large errors at these three temperatures result in the three outliers in the $\Delta \chi^2$ value for real water, deviating from the master curve. For example, at 234.8 K the 3rd peak and the 2nd trough of the structure factor is substantially smaller than the expectation from the temperature trend, which is the major source of the deviation of the fits from the data by our fitting scheme II (see Figure \[fig:sqoofit\]a). We note that scheme I independently fits structure factors at different temperatures, whereas scheme II simultaneously fits structure factors at all the temperatures. Because scheme II uses much less fitting parameters, it is more robust, but at the same time, more sensitive to the errors in the data than scheme I.[]{data-label="fig:noise"}](noise){width="16cm"}

![Analysis of the O-O partial structure factor of a single regular tetrahedron and experimental LDA ice. [**a,**]{} Decomposition of the structure factor of a single regular tetrahedron (Fig. 2b) by two Lorentzian (L1 + L2) and one Gaussian (G) functions. [**b,**]{} Decomposition of the structure factor of experimental LDA ice [@mariedahl2018] by the same functions. The three characteristic peaks obtained from the fit are displayed by blue, orange and magenta curves. The coherence lengths of a single regular tetrahedron and experimental LDA ice were thus estimated to be $\sim2$ Å and $\sim6.5$ Å respectively, from the widths of the corresponding FSDP’s (blue curves). The wave number is scaled by the nearest neighbor O-O distance $r_\mathrm{OO}$.[]{data-label="fig:sqooTetraLDA"}](sqooTetraLDA){width="10cm"}
