Mixture-like behavior near a liquid-liquid phase transition in simulations of supercooled water

Mixture-like behavior near a liquid-liquid phase transition in simulations of supercooled water


In simulations of a water-like model (ST2) that exhibits a liquid-liquid phase transition, we test for the occurrence of a thermodynamic region in which the liquid can be modelled as a two-component mixture. We assign each molecule to one of two species based on the distance to its fifth-nearest neighbor, and evaluate the concentration of each species over a wide range of temperature and density. Our concentration data compare well with mixture-model predictions in a region between the liquid-liquid critical temperature and the temperature of maximum density. Fits of the model to the data in this region yield accurate estimates for the location of the critical point. We also show that the liquid outside the region of density anomalies is poorly modelled as a simple mixture.

64.30.-t, 64.70.Ja, 64.60.De

The possibility that a liquid-liquid phase transition (LLPT) occurs in supercooled water and other tetrahedral liquids (e.g. silicon) continues to be a subject of investigation and debate (1); (2). In the LLPT proposed for water, two phases, a low-density liquid (LDL) and high-density liquid (HDL), become distinct below a critical temperature located in the supercooled regime. While a LLPT for water has yet to be confirmed experimentally, simulation studies have identified unambiguous LLPTs in several model tetrahedral liquids, including water (4), silicon (3), and tetrahedrally coordinated colloids (5).

Long before any discussion of LLPTs, there were recurring proposals that the thermodynamic anomalies of water, such as the density maximum, could be understood if the liquid is modelled as a mixture of two “species” differing in local molecular structure: one of lower density and disorder, and the other of higher density and disorder (6). Following the emergence of evidence for the continuum nature of the local structure and bonding in the liquid under ambient conditions, mixture models of water faded from prominence (7). However, the proposal of a LLPT has renewed interest in mixture models, in which spontaneous LDL-like and HDL-like structural fluctuations play the role of the mixed species (8); (9). Recent experiments have been interpreted as evidence for such mixture-like fluctuations in water (10), although this interpretation is disputed (11).

One commonly discussed mixture model for water (8); (9) is based on the Gibbs free energy of a binary regular solution, given by,


where is the temperature, is the concentration of component B, and are the respective free energies of the pure A and pure B liquids, and is the gas constant. The energy of mixing is quantified by . It is further assumed that the two species can interconvert (), and that the free energy difference between the two pure phases is given by . Here , and are respectively the difference in energy, entropy and molar volume of the two pure phases (e.g. ), and for simplicity are assumed to be constant with respect to and pressure . In this modified regular solution (MRS) model, the equilibrium value of at fixed and is determined by minimizing Eq. 1 with respect to . For , a liquid-liquid critical point occurs at , , and . This MRS model was originally developed by Rapaport to account for melting line maxima in pure liquids (12), and is a member of a large group of two-state models that have been applied to thermodynamic and relaxation phenomena in one-component liquids (13); (14) and crystals (15); (16).

MRS models have illustrated how a LLPT and the thermodynamic anomalies of water may be interrelated. These models have also been reported to be in quantitative agreement with experimental data for water (8); (9); (13), but such comparisons are complicated by two significant challenges. First, the four model parameters , , , and are not known for real water, and so they must be estimated indirectly, e.g. from the properties of the amorphous ices. Second, in order to compute thermodynamic properties from Eq. 1, the free energy function of a reference state [e.g. ] is required. By themselves, the four model parameters determine only the “anomalous” contribution to thermodynamic properties arising from the variation of . Hence, to fit the model to experimental data without knowledge of , a “normal” contribution must first be estimated and subtracted from the data. This process introduces additional and difficult-to-estimate parameters. Consequently, the regime of validity of MRS models for describing behavior near a LLPT, if such a regime exists at all, remains uncertain.

In this Letter, we use simulation data to test if a thermodynamic regime exists in which a MRS model can describe a water-like liquid near a LLPT. We study the ST2 model of water (17), as it provides a context in which the two difficulties described above can be avoided. First, ST2 exhibits a well-characterized LLPT (18); (4). We can therefore determine the four model parameters of the MRS model directly. Second, we identify a property of the local structure in ST2 that allows us to estimate the concentrations of LDL-like and HDL-like species. We compare these concentrations directly to the predictions of the MRS model, thus avoiding the need to decompose thermodynamic properties into normal and anomalous contributions.

Figure 1: (a) - and (b) - projections of the properties of ST2 water. From Ref. (18) we show the line of density maxima (thick black line); density minima (thin black line); maxima (thick blue line); minima (thin blue line); the liquid spinodal (diamonds); the HDL spinodal (down triangles); and the LDL spinodal (up triangles). We also show the locus along which (open red squares); the location of the critical point (open black circle); and the state points falling inside the region defined in the text (crosses). The green line in (a) is a line of slope , and in (b) is the coexistence curve predicted by the MRS model.

Figure 2: (a) at  K, for to  g/cm in steps of  g/cm. (b) at  g/cm, for to  K in steps of  K. For clarity, curves for  K are successively shifted upward by 0.05. (c) A snapshot of a system of ST2 molecules at  K and  g/cm; at this state we find  MPa. Blue molecules have  nm; red molecules have  nm.

Our results are based on molecular dynamics simulations of a system of ST2 molecules. Our runs are conducted at fixed volume , and is controlled using Berendsen’s method (19). Long-range contributions to electrostatic interactions are approximated using the reaction field method. We study a wide range of states: from to 400 K, in 5 K steps; and from density to  g/cm, in steps of  g/cm. Complete details of our simulation procedure are as described in Ref. (18). Fig. 1 summarizes the known phase behavior of ST2. As reported in Ref. (18), a liquid-liquid critical point occurs in the vicinity of  K,  MPa and  g/cm.

First, we seek a criterion for assigning molecules to LDL-like and HDL-like species. In Fig. 2, we analyze the liquid structure near the critical point in terms of the distance from the O atom of each molecule to its fifth-nearest neighbor. Following Ref. (20), we define such that is the average density of fifth-nearest neighbors of an O atom at the origin, as found in a volume element at a distance . So defined, the conventional pair correlation function . In the range of studied here, fifth-nearest neighbors are located over a range of distances that span the first minimum in the O-O pair correlation function, and thus is an indicator of the degree to which the tetrahedral structure of the first coordination shell is disrupted by additional neighbors. At we find that is typically greater than 0.35 nm in the LDL phase (), while in the HDL phase () is typically less than 0.35 nm [Fig. 2(a)]. Near for , has a bimodal shape indicating the presence of distinct populations of LDL-like and HDL-like coordination environments [Fig. 2(b)]. We therefore adopt as a local order parameter for assigning molecules to two species: “A” molecules are LDL-like and have  nm; “B” molecules are HDL-like and have  nm. Fig. 2(c) shows an equilibrium configuration from a separate simulation of ST2 molecules at a state within error of the critical temperature and density. This image, in which A and B molecules are shown in different colors, confirms the presence of large, spatially correlated clusters of each species, as expected near a critical point.

We evaluate , the equilibrium number concentration of B molecules in the system, at each state point as a time average over the equilibrium configurations generated during the run. We also evaluate , to allow us to analyze as well as . Our results for are shown as isotherms as a function of in Fig. 3(a), and as isobars as a function of in Fig. 3(b). Approaching the critical point, the qualitative behavior of is as expected for a liquid mixture approaching a LLPT: Both isotherms (as ) and isobars (as ) become more steeply sloped.

Figure 3: (a) Isotherms of as a function of . (b) Isobars of as a function of . Solid curves are the corresponding predictions of Eq. 1.

Next we evaluate the parameters of the MRS model appropriate for ST2. We confirm the estimates of Ref. (18) for and by examining two - isotherms straddling (inset, Fig. 4). Along the  K isotherm, is monotonic in , while the first sign of a “van der Waals loop” occurs along the  K isotherm. Given the scatter of the data along these isotherms, we estimate  K and  MPa.

We estimate by noting that in the MRS model , where is the “normal” or non-singular contribution to , and is the “anomalous” contribution due to the variation of  (9). in general depends on both and . However, along the critical isotherm near , both and are constant (inset, Fig. 4), and therefore is constant. Accordingly, we estimate from the slope of versus along the critical isotherm in the interval . As we will see below, this range of spans the states near . By averaging the slopes (and their errors) obtained for and  K, we estimate  cm/mol (Fig. 4).

As from above, the MRS model predicts that the value of the isothermal compressibility is increasingly dominated by a term containing , which diverges at the critical point (9). In the model, the maxima of isotherms of occur at for all . As a consequence, the locus along which is a maximum, which converges to the “Widom line” as  (21), should also converge with the locus in the region approaching the critical point. The set of points satisfying is plotted in Figs. 1(a) and (b), along with the locus of maxima reported in Ref. (18). We find that the locus is in excellent agreement with the locus of maxima for  K, and also that the location of the critical point is consistent with the model prediction of in both the - and - planes. Computing the value of at which at predicts  g/cm.

We also find that the locus, both above and below , closely approximates a straight line in the - plane for  K. In the MRS model, both the locus for , as well as the coexistence curve for , follow a straight line in the - plane given by  (9). From a linear fit to the locus from 250 to 280 K, we obtain the slope  MPa/K. We note that the resulting prediction for the coexistence line (with a Clapeyron slope of ) as expected lies between the LDL and HDL spinodal lines in Fig. 1(a).

Figure 4: Isotherms near of as a function of . The lines are linear fits to the data for . For clarity, the data for  K have been shifted downward by 0.5 cm/mol. Inset: Isotherms of versus straddling .

Our estimates for , , , and completely determine the four model parameters (, , , and ) required for Eq. 1. We use these values to obtain the MRS model prediction for , and compare the results to the data for from our simulations (Fig. 3). We find that inside a region , defined approximately as and , the MRS model is in good agreement with the values of computed from our simulations. At the boundaries of and beyond, the agreement rapidly degrades.

To test the robustness of the agreement between the MRS model and our data in the region , we carry out a least-squares fit of the model to our data for in this region, allowing all four model parameters to vary. We select all distinct pairs of isotherms of in that are at least 15 K apart, and fit the model to each of the 21 data subsets so defined. This gives 21 separate estimates for each fit parameter, from which we compute the mean and the standard deviation. The results are in excellent agreement with the values obtained above:  K,  MPa,  cm/mol, and  MPa/K.

Our results thus demonstrate that a mixture model can indeed provide a quantitatively accurate description of a water-like liquid, in the specific case that the liquid exhibits a LLPT. Our MRS model successfully predicts the concentrations of LDL-like and HDL-like structural fluctuations in the region , which lies inside the locus of density extrema, above , and spans the range centered on the Widom line (Fig. 1). The signatures for the onset of this “mixture-model regime” are the merging of the Widom line with the locus, and the observation of a linear Widom line in the - plane. These signatures may be useful for assessing other water-like liquids (either in simulations or experiments) for mixture-like behavior.

Outside , the MRS model does a poor job estimating . Fig. 1(a) shows that the high- boundary of is in the vicinity of  K, which is of the highest ( K) reached by the line of density maxima for ST2 water. For real water, of the temperature of maximum density at ambient pressure (277 K) gives an estimate of  K ( C) for the highest at which mixture-like behavior might be observed experimentally. This estimate supports the view that mixture models are not appropriate for interpreting the behavior of real water at ambient  K (11).

We also note that the MRS model fails for . Fig. 1(b) shows that the coexistence curve predicted by the MRS model (9) is too narrow since, in violation of thermodynamics, it lies inside the estimate of the LDL and HDL spinodal lines. The MRS model is a mean-field theory, and thus the shape of the coexistence curve near obeys , with . However, the LLPT in ST2 water has been shown to belong to the 3D Ising universality class, for which  (4). Hence it is to be expected that the MRS model will underestimate the density difference between the coexisting phases as decreases below . In Fig. 3 we also note that there are significant deviations between the model and the data in the limits of large and small for , highlighting that the nearly pure A and B phases are poorly described by the model. This may also contribute to the discrepancy between the model and the data for the coexistence curve for .

However, within the mixture-model regime defined by , our work shows that the examination of a local structural property as a function of and can yield accurate information concerning the location of the Widom line and the critical point of a LLPT. In our case, the local structure is quantified in terms of , which is determined by the values of individual molecules. A number of simulation studies have used the behavior of and related measures as evidence for a LLPT, in both tetrahedral (3); (20) and non-tetrahedral liquids (22). Our results validate this approach, and further, confirm that the structures relevant to the LLPT in water-like liquids are highly localized, extending no farther than the second coordination shell.

If a mixture-model regime exists for real water, our results suggest that it will be found in the vicinity of the Widom line. States on the Widom line have yet to be studied in experiments on bulk supercooled water, due to the onset of rapid ice crystallization. However, for tetrahedral liquids in which the region of the Widom line is accessible, which may be the case for nanoconfined water (23), our results demonstrate that probes of local molecular structure, in concert with mixture-model concepts, can be used to elucidate the properties of a LLPT and its associated critical point.

We thank K. Fraser and C. Creelman for useful discussions; ACEnet for computing resources; and NSERC, AIF and the CRC program for financial support.


  1. P.H. Poole, F. Sciortino, U. Essmann, and H.E. Stanley, Nature 360, 324 (1992).
  2. O. Mishima and H.E. Stanley, Nature 396, 329 (1998).
  3. S. Sastry and C.A. Angell, Nature Materials 2, 739 (2003).
  4. Y. Liu, A.Z. Panagiotopoulos, and P.G. Debenedetti, J. Chem. Phys. 131, 104508 (2009).
  5. C.W. Hsu, et al., Proc. Natl. Acad. Sci. USA 105, 13711 (2008).
  6. W.C. Röntgen, Ann. Phys. Chem. N.F. XLV 91, 1891.
  7. G. Malenkov, J. Phys.: Condens. Matter 21, 283101 (2009).
  8. C.T. Moynihan, Mater. Res. Soc. Symp. Proc. 455, 411 (1997).
  9. E.G. Ponyatovsky, V.V. Sinitsyn, and T.A. Pozdnyakova, J. Chem. Phys. 109, 2413 (1998).
  10. C. Huang, et al., Proc. Natl. Acad. Sci. USA 106, 15214 (2009).
  11. G.N.I. Clark, et al., Proc. Natl. Acad. Sci. USA 107, 14003 (2010).
  12. E.J. Rapoport, J. Chem. Phys. 46, 2891 (1967); J. Chem. Phys. 48, 1433 (1968).
  13. H. Tanaka, Europhys. Lett. 50, 340 (2000); J. Chem. Phys. 112, 799 (2000).
  14. C.A. Angell, C.T. Moynihan, M. Hemmati, J. Non-Cryst. Solids 274, 319 (2000).
  15. S. Strässler and C. Kittel, Phys. Rev. A 139, 758 (1965).
  16. I.L. AptekarÕ and E.G. Ponyatovsky, Fiz. Met. Metalloved. 25, 777 (1968); Fiz. Met. Metalloved. 25, 1049 (1968).
  17. F.H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 (1974).
  18. P.H. Poole, I. Saika-Voivod, and F. Sciortino, J. Phys.: Condens. Matter 17, L431 (2005).
  19. H.J.C. Berendsen, et al., J. Chem. Phys. 81 3684 (1984).
  20. I. Saika-Voivod, F. Sciortino, and P.H. Poole, Phys. Rev. E 63, 011202 (2001).
  21. L. Xu, et al., Proc. Natl. Acad. Sci. USA 102, 16558 (2005).
  22. B. Boates and S.A. Bonev, Phys. Rev. Lett. 102, 015701 (2009).
  23. D. Liu et al., Proc. Natl. Acad. Sci. USA 104, 9570 (2007).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description