The Putative Liquid-Liquid Transition is a Liquid-Solid Transition in Atomistic Models of Water, Part II

This paper extends our earlier studies of free energy functions of density and crystalline order parameters for models of supercooled water, which allows us to examine the possibility of two distinct metastable liquid phases [J. Chem. Phys. 135, 134503 (2011) and arXiv:1107.0337v2]. Low-temperature reversible free energy surfaces of several different atomistic models are computed: mW water, TIP4P/2005 water, SW silicon and ST2 water, the last of these comparing three different treatments of long-ranged forces. In each case, we show that there is one stable or metastable liquid phase, and there is an ice-like crystal phase. The time scales for crystallization in these systems far exceed those of structural relaxation in the supercooled metastable liquid. We show how this wide separation in time scales produces an illusion of a low-temperature liquid-liquid transition. The phenomenon suggesting metastability of two distinct liquid phases is actually coarsening of the ordered ice-like phase, which we elucidate using both analytical theory and computer simulation. For the latter, we describe robust methods for computing reversible free energy surfaces, and we consider effects of electrostatic boundary conditions. We show that sensible alterations of models and boundary conditions produce no qualitative changes in low-temperature phase behaviors of these systems, only marginal changes in equations of state. On the other hand, we show that altering sampling time scales can produce large and qualitative nonequilibrium effects. Recent reports of evidence of a liquid-liquid critical point in computer simulations of supercooled water are considered in this light.


I. INTRODUCTION
This is our second paper examining whether molecular simulation provides support for the hypothesis that supercooled water possesses two distinct liquid phases with a reversible coexistence line ending at a critical point. 1 In the first (Paper I), 2 we described our results for pertinent free energy functions of three different models: the mW model of water, 3 a variant of the ST2 model for water, 4 and the Stillinger-Weber model for Si. 5 Each of these models has an equilibrium liquid phase with pair distributions and thermodynamic anomalies like those of water, and each has an equilibrium phase transition like that of water-ice freezing.For each, others have claimed numerical evidence of liquid-liquid transitions at supercooled conditions, 1,[6][7][8][9][10][11][12][13][14][15][16][17][18][19] but our calculations described in Paper I reveal no such behavior.Rather, for each system we found that there can be at most one stable or metastable liquid phase, and this liquid can coexist with crystalline ice.Here, we establish that results contrary to our findings derive from lack of equilibration of slow fluctuations in long range order.We also present new calculations for several additional models reaching consistent conclusions in each case.Specifically, for computersimulation models of water that exhibit liquid and ice-like phases, there is no liquid-liquid transition, but there is non-equilibrium coarsening of ice that others have misinterpreted as evidence of a liquid-liquid transition.
a) Electronic mail: chandler@berkeley.edu Figure 1 shows the relevant part of water's phase diagram and corresponding free energy surfaces.The liquid is the stable equilibrium phase for temperatures above the melting temperature, i.e., T > T m , and it is unstable below a stability temperature, i.e., T < T s .In between, for T s < T < T m , the liquid is metastable with respect to crystal ice.Throughout much of that intermediate region, structural reorganization of water is slow, and it becomes slower in a super-Arrhenius fashion as temperature is lowered. 20This sluggishness can present problems for straightforward molecular simulation, as noted below, but it is not so sluggish to prevent certain crystallization of water when the liquid is cooled close to or below T s .Coarsening of water in that regime occurs on the time scale of microseconds -fast for experiment, but slow for simulation. 21All speculations on the existence of a liquid-liquid phase transition in water locate that transition near or below T s , the so-called "no-man's land" for liquid water.As such, it is difficult for experiments to prove or disprove the liquid-liquid hypothesis.It is left to simulation, which can reversibly control crystallization, to see if such an idea could be correct within the purview of statistical mechanics for plausible models of water or water-like systems.
Calculations of free energy functions of relevant order parameters are required when using simulation to establish phase behavior. 22As noted, such calculations can be difficult, especially for supercooled water because fluctuations in this regime are collective and slow.To address this difficulty and sort out the phase behavior of supercooled water, we have found it convenient to consider two order parameters.One is molecular density, ρ, that dis-FIG.1. Phase diagram, free energy surfaces and typical configurations of cold water.Typical of all systems considered in this paper, the specific pictures render results from molecular simulations of one particular model.Quantitative scales of temperature, pressure and free energy depend upon simulation model, and these scales are omitted here because this figure serves a qualitative purpose only.For experimental water, the phase diagram covers pressures p ranging from 1 bar to 3 kbar, and temperatures T ranging from 150 K to 300 K. Tm and Ts stand for temperatures at the melting and liquidstability lines, respectively.Blue region in the p-T plane is where liquid is stable, red region is where liquid is metastable with respect to ice I, and grey region is where liquid is unstable (i.e., it is the liquid's "no-man's land").Corresponding free energy surfaces are shown above as functions of density ρ and crystal-order variable Q6.Corresponding molecular configurations shown below are cuts through a simulation box at the ends of trajectories that are initiated in a liquid configuration and that run for times shorter than required to crystalize the entire sample.Molecules located in crystal-like regions are colored red.
tinguishes different amorphous phases.The other can be the Steinhardt-Nelson-Ronchetti Q 6 that distinguishes an amorphous phase from a symmetry-broken crystalline phase. 23,24The two variables fluctuate on significantly different time scales.For example, the liquid structural relaxation time (i.e., the equilibration time for ρ) around T ≈ T s is of order 10 −8 s or shorter, whereas the relevant equilibration time for Q 6 in this regime is the time to form a crystal, 10 −6 s or longer.This wide separation of time scales is typical of systems undergoing crystallization transitions. 25In the case of water, we will see in this paper that it is a principal source of confusion in simulation studies that claim evidence for a liquid-liquid phase transition.
To foreshadow this point, consider the equilibrium joint distribution function for the order parameters, P (ρ, Q 6 ).It is related to the free energy (or reversible work) surface for these variables in the usual way: where k B is Boltzmann's constant.This is the free energy function illustrated in Fig. 1.Over time scales large compared to those of liquid relaxation but possibly not large compared to those of crystal formation, the joint distribution is in general a non-equilibrium distribution, where P (ρ|Q 6 ) is the equilibrium distribution for ρ given a specific value for Q 6 , and P ne (Q 6 , t) is the nonequilibrium distribution for Q 6 .The non-equilibrium distribution depends upon the protocol with which the system is prepared, and its time dependence is irreversible.For large enough t, presuming ergodicity, P ne (Q 6 , t) approaches the equilibrium P (Q 6 ).But this limit can require simulation times thousands of times longer than those needed to equilibrate ρ.Not accounting for this behavior can give the illusion of a reversible polyamorphism because the non-equilibrium free energy, −k B T ln[P (ρ|Q 6 ) P ne (Q 6 , t)], can have a low-Q 6 basin for times shorter than those required for Q 6 to diffuse towards its equilibrium crystal value at high Q 6 .This possibility, which we refer to as "artificial polyamorphism," can be appreciated by comparing the free energy surfaces shown in Fig. 1.In particular, imagine studying the system on time scales where Q 6 can diffuse over no more than the left halves of the pictured free-energy panels.P ne (Q 6 , t) would then be peaked at a low value of Q 6 , even for cases where a high Q 6 value would be the correct equilibrium value.Thus, if Q 6 is limited in this way to small values, the low-temperature (i.e., left-most) panel would then yield a pseudo free energy, −k B T ln P ne (ρ, Q 6 , t), with an illusory "amorphous basin" at a density lower than that of the metastable liquid.This irreversible behavior was discussed in the Supplement to Paper I. 26 Specifically, we showed that for liquid water at pressures and temperatures in or close to "no man's land," small values of Q 6 will survive while the crystal phase begins to coarsen.The bottom left of Fig. 1 shows a configuration of water in that regime.
Section II of the current paper provides a quantitative theoretical analysis of this behavior.It shows specifically how the polyamorphism of Refs. 10, 11, 15, 16, and 19 is an irreversible effect reflecting the time-scale separation between fluctuations in ρ and fluctuations in Q 6 .During the time of coarsening, the faster order parameter, ρ, fluctuates between typical crystal and liquid values.References 15 and 16 report this type of behavior, which they call "phase flipping."On the time scale of the flipping, the drift in P ne (Q 6 , t) can be almost imperceptible, but drift it does.The authors of Refs.15 and 16 describe the flipping as evidence of a second metastable liquid.The analysis of Section II shows that this flipping is not a consequence of such metastability, but rather the coarsening of the crystal from the unstable or nearly unstable liquid, occurring steadily and irreversibly on a time scale long compared to those considered in Refs.15 and 16.This finding could already be anticipated from the Supplement to Paper I 26 and from Moore and Molinero's previous study of crystallization of mW water. 27he specific pathways by which simulated models coarsen depend upon system size.For example, free energy barriers separating coexisting equilibrium basins scale as N 2/3 , manifesting the presence of an interface. 28imilarly, with increasing N , the thermally accessible widths of the two basins decrease as N −1/2 , and the width of the barrier grows. 28Due to the growing barrier height and barrier width, the frequency of spontaneous events carrying a system between stable phases is therefore vanishingly small in the limit of large N .Similarly, in or near the region of liquid instability (i.e., T close to or lower than the crossover temperature T s ), the slope towards the crystal basin will increase in magnitude as N increases.This behavior will affect the rate of "phase flipping" during coarsening, giving the transient impression of a non-equilibrium barrier between low-density and high-density states that changes with N .
These finite-size effects are fundamental to the nature of phase transitions.Establishing the existence of a phase transition requires studying system-size dependence, for example, by computing changes in free energy barriers with respect to changing N .No such computations have yet been performed for putative liquid-liquid transitions in models of water that exhibit water-like structure of the liquid and crystal phases.To do so requires algorithms that can attend to the collective nature of systems undergoing phase transitions.Free energy methods are among the tools that are suitable for the task, provided they are combined with trajectory algorithms that are appropriately efficient and reversible. 29n Section III, we detail how pertinent free energies can be computed for supercooled water, and we consider different variants of the ST2 model as applications.Juxtaposition of free energy surfaces for three different variants indicates that reasonable changes in electrostatic boundary conditions do not change general phase behaviors.Section IV presents free energy surfaces obtained for other systems: the mW of water, the TIP4P/2005 model of water, 30 and the SW model of Si.In every case, the models are found to exhibit one stable or metastable liquid phase plus ice-like crystal phases.Coexistence between two distinct liquid phases does not occur.A summary of our findings is given in Section IV, and Appendices A, B and C present further details and results.
Reversibility is particularly important to the issues addressed here and in Paper I. Distinct reversible phases can be interconverted, with properties that are independent of the paths by which they are prepared.Reversible liquid phases are thus not the same as amorphous solids or glasses.The former are reversible and ergodic, so their measured stationary behaviors are independent of history.The latter, like high-density or low-density amorphous ices (HDA and LDA), are not ergodic, so their behaviors depend much on history (i.e., preparation protocols).Observed transitions between HDA and LDA phases, 31 therefore, are necessarily different than reversible liquid-liquid transitions.Melting amorphous ice to produce a non-equilibrium liquid that then crystalizes is different too. 32rystallization following the melting of glass 33 and crystallization following the rapid quench of water into the liquid's "no-man's land" 34 are much like nonequilibrium dynamics evolving from low to high Q 6 on the middle and left free energy surfaces pictured in Fig. 1, an observation worthy of future study.But these interesting non-equilibrium processes and the transitions between different amorphous solids of water are not our focus in this work.Rather, we are concerned with whether water-like systems when constrained to not freeze can exhibit two distinct liquid phases.6][37][38] If, instead, reversible molecular simulation models exhibit only ice and one liquid, then the symmetry differences between ice and liquid exclude the possibility of an associated critical point.We believe the systematic evidence provided herein and in Paper I indicates that there is only one liquid and no low-temperature critical point.

II. THEORY OF COARSENING AND ARTIFICIAL POLYAMORPHISM IN COMPUTER SIMULATION OF WATER
This section provides a quantitative theoretical analysis showing the difficulty in obtaining correct reversible free energy surfaces of supercooled water.We do so by examining the effects of time-scale separation for dynamics on a reversible free energy surface.The particular surface we employ is the free energy F (ρ, Q 6 ) derived in Paper I for a variant of the ST2 model.This free energy Negative logarithm of the non-equilibrium distribution for crystal order, Q6, as it relaxes from the liquid state.It is computed from the Fokker-Planck equation with the free energy surface given in Panel (a) under the assumption that the density, ρ, remains at equilibrium with the instantaneous value of Q6.(c) and (d) Non-equilibrium pseudo free energy surfaces computed from Eq. 1 at two intermediate stages of relaxation, t = 10 τQ 6 and t = 1, 000 τQ 6 .The unit of time, τQ 6 , is the autocorrelation time for Q6 fluctuations in the liquid basin (i.e., at small Q6).The reduced density is ρ = (ρ − ρ xtl )/∆ρ, where ρ xtl is the mean density of the crystal basin (i.e., at large Q6), and ∆ρ is the difference between the mean densities of the liquid and crystal basins.Contour lines are separated by 1 kBT and statistical uncertainties are about 1 kBT is shown in Fig. 2(a).The methods used to obtain that surface are the subject of the next section, but here we only need to assume that there is such a surface, and that it is qualitatively like the surface shown in Fig. 2(a).
Free energy surfaces for several models are derived in Sections III and IV.The generic features of the free energy surfaces are the same for all the models.There is a liquid basin at small Q 6 and large ρ, and a crystal basin at large Q 6 and small ρ.For a given temperature, T , the relatively stabilities of the basins are controlled by the pressure.The free energy at pressure p is related to the free energy at pressure p + ∆p in the usual way, Accordingly, lowering pressure tips the surface towards the crystal basin, and raising pressure tips it towards the liquid basin.
In cases where the crystal is stable but the system is prepared in the liquid, an irreversible drift towards the crystal will occur.To the extent that ρ and Q 6 are the principal slow variables, this coarsening can be described in terms of motion on the F (ρ, Q 6 ) surface.By using this perspective, and specifically by adopting the free energy surface pictured in Fig. 2(a) , we illustrate here the generic behavior of early-stage coarsening of ice.The behavior is not specific to the particular free energy surface.Rather, it is general consequence of a separation of time sales, where the density ρ equilibrates on time scales that are at least two orders of magnitude shorter than the time scales on which Q 6 fluctuates.Such separations of time scales are typical in natural and computer simulated supercooled water.See Table I.
A. Time dependence of Pne(Q6, t) Due to the separation in time scales, relaxation of P ne (Q 6 , t) can be accurately estimated by assuming density ρ is always in equilibrium with the current value of Q 6 .An appropriate Fokker-Planck equation 40 is therefore where The quantity F (Q 6 ) is the equilibrium free energy for the crystal-order parameter, β = 1/k B T , and D = (δQ 6 ) 2 /τ Q6 is the diffusion constant projected along the Q 6 direction.The quantity (δQ 6 ) 2 ≈ 0.01 is the mean-square fluctuation of Q 6 in the liquid basin for the 216-molecule system considered in Fig. 2(a).The longtime limit is set by the diffusion constant, For quantitative treatments of the ultimate equilibration (i.e., of the final stages of crystal coarsening), Eq. 3 could be generalized to include Q 6 -dependence and memory effects in D. Such generalization could account for the complexity of pathways by which multiple ordered domains reorganize and connect and would be expected to increase the timescales for equilibration.These refinements are unnecessary for the current analysis of earlystage coarsening, where Q 6 does not progress far from its values in the liquid.
We have integrated Eq. 3 using a first-order finite difference approach with a small enough discretization of  shows how the distribution evolves in time from an initial Gaussian distribution centered in the liquid region of Q 6 .What is notable is that the relaxation to equilibrium takes orders magnitude longer than the basic timescale, τ Q6 .
With this time evolved probability distribution, we have used Eq. 1 to estimate a non-equilibrium joint free energy, This pseudo free energy function for two intermediate times is shown in Figs. 2 (c) and (d).At the first of these intermediate times, t = 10 τ Q6 , F ne (ρ, Q 6 , t) exhibits two minima at low values of Q 6 , and these minima are separated by a small barrier of a few k B T .The low density basin is centered at the mean density of the crystal, and the high density basin is centered at the mean density of the liquid.We use the reduced density variable, ρ = (ρ − ρ xtl )/∆ρ, to emphasize these connections to the crystal and liquid basins.This behavior shown in Fig. 2 (c) is precisely the behavior found in Refs.10, 11, 15, and 19 -both the bistability and the length of time allowed for equilibration.Those workers find τ Q6 ≈ 10 8 Monte Carlo sweeps, and they use 10 9 sweeps to estimate free energies.The lowdensity liquid minimum eventually disappears, but the time scale for that to occur is orders of magnitude longer than considered in Refs.10, 11, 15, and 19.To further illustrate the connection between the nonequilibrium calculation shown here with the finite-time sampling results of Ref. 10, 11, 15, and 19, Fig. 3 shows the effects of pressure variation on the pseudo free energy surfaces.Upon re-weighting to lower pressure, Fig. 3 (a), or to higher pressure, Fig. 3 (c), one of the disordered minima disappears.The specific dependence on re-weighting and the relative locations of the minima are also consistent with the results of Refs. 10, 11, 15, and 19.We have thus reproduced the principal results of those papers by identifying the limited time over which the system was allowed to equilibrate.See for example, Fig. 2 of Ref. 15, as we have purposely used a similar color code in our Fig. 3 to emphasize the similarities of our finite-time results with the free energies reported in that paper.
References 10, 11, 15, and 19 employ slightly different variants of the ST2 model and slightly different temperatures.In Appendix B, we consider another variant and another temperature to illustrate that the results of our analysis are generic.The analysis uses the Fokker-Planck equation in order to elucidate the generality of the phenomena.An explicit simulation calculation establishing the same conclusion is also provided in Appendix B.

B. Phase flipping
The time-dependent pseudo free energies illustrated above also shed light on previous reports of phase flipping between two seemingly distinct liquids. 15,16 mittently between smaller amplitude motion while the system is globally liquid-like.This behavior is expected for trajectories of ρ when driven by the pseudo free energy surface graphed in Fig. 3 (b).Indeed, Fig. 4 shows representative trajectories obtained by running overdamped Langevin dynamics 40 on the free energy surface of Fig. 2(a) with the time-scale separation τ Q6 = 100 τ ρ .The trajectories were initiated in the liquid basin at temperature T = 235 K. 42 The trajectories look exactly like those presented as phase flipping in Refs.15 and 16.(The structural relaxation time in those studies is τ ρ 1 ns.)Thus, so-called "phase flipping" is not a flipping between distinct liquid phases.It emerges in a single liquid phase because there exists a separation of timescales between density fluctuations and long range order fluctuations.Large density fluctuations occur as the system is attempting to crystallize, accompanied by the formation of sub-critical nuclei formation that shrink and cause further large density fluctuations.
This behavior is transient, as the pseudo free energy is ever-changing with time.An analysis of time series would elucidate this nature of the phenomenon.Specifically, the mean transition time will show a dependence on trajectory length, reflecting the non-stationarity of the system.Consequently, the distribution of transition times will deviate from simple Poisson statistics.Those offering phase flipping as evidence for two distinct liquids have not provided a quantitative analysis supporting such statistics.Rather, studies illustrated in Ref. 17 catalogue many independent trajectories some of which exhibit anomalously long quiescent periods between transitions lasting several times 100 ns ≈ 100 τ ρ .Data supplied for the ST2 model in Ref. 17 also illustrate that phase flipping occurs at pressures far below a proposed critical pressure, in apparent contradiction to the authors' claims that a critical point exists.The model behavior we discuss is a consequence of liquid-crystal coexistence, which by symmetry must not have a critical point and would explain the insensitivity of this behavior to pressure.
The only quantitive analysis provided in studies of phase flipping is the calculation of bi-modal density distributions, which are then fit to a universal Ising form, from which supposed critical parameters are extracted. 16,18The barostat used in those simulations is well known to not reproduce correct equilibrium density fluctuations, 43 making such a calculation difficult to interpret.

III. CALCULATION OF REVERSIBLE FREE ENERGY SURFACES
The previous section emphasizes the significance of a time-scale separation and the pitfalls for simulation that result from it.In this section, we detail procedures that can overcome the problems inherent in that time-scale separation, and we apply the procedures to the ST2 model to evaluate the liquid-liquid hypothesis in that case.In light of recent speculations, 10,11,15,19 we also analyze how these equilibrium surfaces are affected by changes in model parameters and boundary conditions.

A. Methodology
We have carried out Monte Carlo simulations with constant number of molecules, N , pressure, p, and temperature, T .
This ensemble is appropriate for determining conditions of phase coexistences and is well suited for sampling dense liquid and crystalline states.Density, ρ = N/V , fluctuates with p and N fixed because volume, V , fluctuates.The alternative Grand Canonical calculation for sampling global density fluctuations 22 can suffer from poor acceptance ratios (e.g., acceptance ratios are reported 10 to be as low as 10 −5 ) and can also be ill-defined in applications to crystals. 44ajectory algorithm designed for the task.The move set we employ for our Monte Carlo calculations is chosen to mitigate long correlation times expected for supercooled liquids and coarsening crystals.Specifically, to allow for collective reorganizations, we use a hybrid Monte Carlo algorithm 45 that propagates an initial configuration with Boltzmann distributed velocities under symplectic, norm preserving, molecular dynamics.For the models with internal degrees of freedom we use the SETTLE integrator while for single-site models we use a velocity Verlet integrator. 46The configuration is integrated for a time n δt, where δt is the integration timestep.Each move is accepted with a Metropolis criterion, 22 so energy need not be conserved and consequently δt need not be small.
In practice, we generally range from δt ≈ 5 − 30 fs and the number of steps in a trial trajectory, n, to vary between 1-20 depending of the steepness of the free energy landscape.The choices are made systematically to minimize correlation times.Volume moves are used at a ratio of 2 hybrid Monte Carlo moves to 1 trial volume displacement.The relatively large value of δt serves to swiftly propagate dynamics over long time scales.We have found that at supercooled conditions for N ≈ 200, it significantly reduces correlation times for structural relaxation relative to energy conserving dynamics.For instance, in the case of the ST2 model discussed in detail below, the characteristic structural relaxation times under these moves are between 10 2 to 10 3 Monte Carlo steps, depending on the specific value of density and temperature.In contrast, single particle Monte Carlo moves reported previously 10,11,15,19 yield structural relaxation times that are between 10 5 to 10 7 Monte Carlo steps, and molecular dynamics 16 yields structural relaxation times that are between 10 6 to 10 8 integration steps.Accounting for the factor n = O(10), we see that our choice of hybrid Monte Carlo moves is computationally more efficient by 1-3 orders of magnitude over single particle moves and by 2 orders of magnitude over molecular dynamics.This remarkable speed up must in part reflect the highly non-linear and correlated nature of dynamics at supercooled conditions.Controlled biasing methods.Fluctuations that result in phase transformations are exponentially rare at conditions of coexistence or modest supercooling.Standard methods of umbrella sampling get around the rareevent problem by adding biasing potentials, W (x N ), to the Hamiltonian in order to enhance occurrences of otherwise improbable fluctuations.Re-weighting configurations correct for the biasing.Here, x N stands for a point in the configuration space for the system.The form of an added energy function can be chosen for convenience, but in order to guarantee the system will reach a stationary state it must be time-independent.
Free energy methods that do not strictly adhere to this condition, such as meta-dynamics 47 and Wang-Landau sampling, 48 converge only conditionally in the limit that the biasing degree of freedom is the slowest mode or when the change in the biasing term asymptotes to zero.This issue is particularly relevant for supercooled water because pathways beyond the early stages of coarsening generally involve several slow variables in addition to the global crystal order parameter, Q 6 .An incorrect free energy estimate will be obtained if one or more of those slow variables are not controlled in meta-dynamics or Wang-Landau algorithms.We do not exclude the possibility of inventing adaptive methods that could account for the physical behaviors of supercooled water and be more efficient than the approach we employ, but such adaptive methods are not currently standard.
The order parameters, ρ and Q 6 , are chosen to distinguish phases of broken symmetries expected to result in water-like models at low temperatures.Our previous study demonstrated that Q 6 is a sufficiently sensitive order parameter to distinguish globally ordered from disordered states accompanying a freezing transition. 2 This virtue noted, it must also be appreciated that Q 6 deviates from its disordered value only after a substantial amount of orientational order has developed in the system.This fact is demonstrated below and also in the Supplement to Paper I. 26 The umbrella biasing potentials we employ are of the form where ρ(x N ) and Q 6 (x N ) are the density and crystalsymmetry order parameters, respectively for configuration x N .The biasing potentials, with its force constants κ and k, keep these order parameters close to their target values ρ * and Q * 6 .Each pair of target values defines a specific sub-ensemble or so-called "window" in configuration space.After collecting statistics in one window, the window is moved by changing the pair of target values, ρ * and Q * 6 , whereupon statistics in the new window are collected.The procedure is carried out throughout the ρ-Q 6 plane, making sure that passage from one region to another is fully reversible, and that adjacent regions have sufficient overlap of statistics to enable further analysis.
For this purpose, we find that κ in the range of 500 to 10,000 k B T and k in the range of 1,000 to 2,000 k B T cm 6 g −2 is satisfactory.Statistics gathered in these biased ensembles are then unweighted and the free energy differences between each ensemble are estimated using MBAR. 49titching together data from different sub-ensembles yields the joint probability and thus free energy F (ρ, Q 6 ; p, T ), where p and T , respectively, denote the pressure and temperature at which the Monte Carlo trajectory is carried out.Assuming this function is determined accurately over the relevant range of densities, free energies at other relevant pressures are determined from Eq. 2, i.e., by straightforward re-weighting. 22,50are applied to ensure reversibility.In each window, we initially equilibrated for up to 100 structural relaxation times as evaluated for the larger of either the initial liquid density or the equilibrium of bias window density.Production runs were between 50 and 1000 structural relaxation times, depending upon the length of time required to obtain reliable statistics as judged from cumulative averages for Q 6 and potential energy.Here, structural relaxation time refers to the simulation time, t, required for mean square fluctuations in structure factors to decay from their initial to 90% of their relaxed value.
Three different techniques for generating initial conditions were used.Initial seeds were created by cooling an equilibrium liquid initially prepared at T = 330 K at a rate of 10 K/ns until it reached the target temperature.Seeds from this procedure were biased into different windows in steps between adjacent windows, by gradually changing parameters of the biasing potential W (x N ), Eq. 7, and with re-equilibration runs in between each step.At high Q 6 , the crystal that was spontaneously formed using this procedure in all cases was a defected Ice Ic.For second generation seeds, we assumed that the spontaneously formed crystal was the relevant solid phase, so we prepared a perfect Ice Ic configuration, which was used to sample intermediate and high Q 6 states as well as bias them into low Q 6 regions to sample liquid states.Third generation seeds were obtained by melting an Ice Ic configuration and then using states along the melting trajectory to seed intermediate Q 6 windows.These configurations were subsequently biased into the high and low Q 6 regions, again by gradually changing the parameters of W (x N ), and again with new re-equilibration runs.In total, of the order of 10 3 independent biasing windows are used in the calculation of an individual free energy surface.Specifically, we generated about 800 for each SW and mW free energy surface and about 2500 for each ST2 and TIP4P/2005 model free energy surface.
One issue to keep in mind if one chooses to use parallel tempering and replica exchange 22 is that the shape of the 5kBT and statistical errors over the surfaces average to less then 1 kBT .Quantitative features will change with system size.For example, as N grows, the mean value of Q6 in the liquid basin will vanish as 1/N 1/2 , while in the crystal basin it will remain finite.
free energy surface changes with temperature.A liquid basin exists for T > T s , but it does not exist for T T s .Moving between replicas that traverse this crossover will produce a transient shadow of the higher temperature liquid basin in the lower temperature replica.Appendix C illustrates how overlooking this behavior can lead to poor free energy estimates and false impressions of a second liquid basin.By using a control variable other than temperature, related methods might be designed to increase the relevant diversity of sampled configurations, but the underlying physics that makes a particular variable either applicable or inapplicable must always be considered.
In principle, free energies can be correctly computed by any number of different methods, provided the procedures are truly reversible.This property, reversibility, was explicitly checked in our calculations by constructing plots of all two-dimensional histograms and checking for hysteresis.Estimates of errors in free energy differences were made by computing overlaps and gradients of the distributions obtained by various routes.These steps, including bidirectional biasing to and from the crystalline phase and creating many independent realizations of initial conditions, follow standard practices for computing free energies articulated in reviews such as Ref. 29.

B. Results for different variants of the ST2 model
Figure 5 shows free energy surfaces we have computed for three different versions of the ST2 model, variants that differ only in the manner by which long-ranged forces are computed.The phase behaviors in each case are similar, with one liquid basin and one crystal basin.Indeed, the existence of singularities in a partition sum is usually not sensitive to subtle changes in potential energy function.This is true because the existence of a phase transition is mostly dictated by dimensionality and the general form and symmetry of the potential energy function. 51In contrast, the locations of the singularities (e.g., temperatures and pressures of coexistence) are often sensitive to subtle changes. 52Recent reports on variants of the ST2 model 15,19 have hypothesized that differences in electrostatic boundary condition and non-electrostatic cutoff parameters can account for why our previously published results 2 found no liquid-liquid phase transition where others suggest it does exist.The results shown in Fig. 5 dismiss this hypothesis.
The ST2a model.Panel (a) is for the model we used in Ref. 2, but at a temperature slightly lower than that considered in our earlier work.Our prior reported calculations were for T = 235 K, whereas Fig. 5(a) is for T = 230 K.The results for T = 230 K differ very little from those at T = 235 K.The model employs a modification of the original ST2 potential for water. 5The modification includes forces from the long ranged electrostatics that were neglected in the original model.The inclusion uses an Ewald summation with conducting boundary conditions.The non-electrostatic Lennard-Jones potential is truncated and shifted at 7.5 Å.This model was referred to as the mST2 model in Ref. 2. Here, we identify it as the ST2a model, and distinguish it from a modified ST2 model with insulating boundary conditions.
The ST2b model.Panel (b) shows the free energy we have computed at the same temperature as considered in Panel (a), but now with insulating boundary conditions.In addition, a tail correction has been added to account for the portion of the Lennard-Jones potential neglected in truncation.Panel (b), therefore, shows our results for the equilibrated free energy surface of precisely the variant of the ST2 model used in Refs. 10 and 15.Here, we identify it as the ST2b model.As in Panel (a), the surface exhibits a single crystal basin at large Q 6 and a single liquid basin at small Q 6 .The only significant differences between the two surfaces in Panels (a) and (b) are in the location of the liquid basin and the relative stability of the crystal.The density for the liquid basin in Panel (b) is higher than that in Panel (a), and the crystal stability in Panel (b) is reduced from that in Panel (a).Reasons for these differences will be discussed shortly.
The ST2c model.Finally, Panel (c) shows our results for the ST2 model using the variant described in Refs.11 and 19, which we identify as the ST2c model.It uses a reaction-field treatment of electrostatic interactions and a Lennard-Jones tail-correction.This Panel (c), plus Panels (a) and (b) together with the results of Paper I show that small changes in temperature and variation in boundary condition have only marginal effects on the phase behavior of ST2 water.
When effects of changing potential energy or temperature are marginal, the nature of those effects are easily interpreted and computed with the identity 28 Here, ∆F (ρ, Q 6 ) is the change in free energy surface due to changing the temperature or Hamiltonian by the amount ∆H/k B T .The angle brackets refer to the ensemble average over configurations with the original temperature and Hamiltonian, and with the order parameters fixed at the values indicated by the subscripts.When ∆H/k B T is large, or when its effects are large, calculations with this formula will yield poor estimates of ∆F (ρ, Q 6 ) because exponential averages are slowly converging. 29On the other hand, if ∆H is small, or if its effects are small, averages converge reasonably quickly.In that case, Eq. 8 becomes a computationally convenient estimator.We have checked explicitly that this measure of marginal behavior is satisfied with respect to the above cited variations for the ST2 model at T ≈ 230 K and N ≈ 200.For example, the typical size of ∆H/k B T for the comparison between Figs. 2 (a) and 5 (b) is 30% of the root-mean-square fluctuations in the net energy per k B T .The conducting boundary condition for Ewald sums, used for Fig. 5(a), strictly cancels electrostatic surface potentials in the energy function.These surface potentials arise from instantaneous polarization fluctuations and transient dipole moments of the total system.The conducting boundary condition is generally chosen for use with molecular dynamics simulations because the alternative, insulating boundary condition, typically introduces discontinuities in the potential energy whenever a molecule crosses the the periodic boundary. 53For a dipole disordered system, like liquid water and ice Ih, this boundary condition should be irrelevant as its energy will average to zero.
Assuming the Ewald parameters are chosen such that energy is well converged in both cases, the pressure of the system with insulating boundary conditions is larger than that with conducting boundary conditions by an amount ∆p surf .This pressure difference is given by 54 Here, is the dielectric constant of the surrounding medium, M is the total system dipole, and V is the volume.For a disordered system, the average M is zero in the thermodynamic limit, and M 2 fluctuates with values that are a extensive in system size.Accordingly, Eq. 5 shows that this strictly positive contribution to the pressure vanishes as ∼ 1/V .More discussion on subtleties involved in the implementation of the Ewald summation is given in Appendix A.
The non-electrostatic part of the ST2 potential is a 6-12 Lennard-Jones interaction.The simulations yielding Panel (a) truncate and shift this potential to zero beyond the oxygen-oxygen distance r c = 7.5 Å.This truncation produces a net potential energy that is slightly higher than that with no truncation.An accurate correction to the potential energy that accounts for the neglected tail is the mean-field estimate −16π LJ σ 6 ρN/3r 3 c , where the Lennard-Jones energy and length parameters are LJ and σ, respectively.Differentiation with respect to volume thus gives a tail correction for the pressure, 4 Thus, while the two simulations yielding Panels (a) and (b) in Fig. 2 are both carried out at p = 2.2 kbar, the effective pressure of the former differs from the latter by the amount ∆p = ∆p surf + ∆p tail , so that the mean density of the latter will be higher than that of the former by the amount where the term in square brackets is the compressibility, (∂ ρ /∂p) T .Evaluating ∆p from the formulas above for the surface and tail corrections, and estimating the meansquare density fluctuations from the widths of the liquid basins in the free energy surfaces, we find ∆ρ ≈ 0.1 g/cc, in harmony with the differences seen in Panels (a) and (b) of Fig. 5.The system is much less compressible in the crystal basin than in the liquid basin (i.e., the density fluctuations are smaller in the crystal than in the liquid), and as a result, the shift in position of that basin between Panels (a) and (b) is much less than that found for the liquid.
The relative stability of the crystal basin in Panels (b) is notably less than that in Panel (a).This juxtaposition is another manifestation of the fact that with a given external pressure p, the effective pressure of the ST2b model is higher than that of the ST2a model.
In the reaction-field treatment used to compute F (ρ, Q 6 ) of Fig. 5(c), Coulomb interactions are summed directly up to a cutoff distance, R c , and contributions from larger separations are approximated as those from an ideal polarizable continuum.With this approximation, and assuming the medium has a large dielectric constant, the term must be added to the potential energy evaluated with the truncated direct Coulomb sums.Here, µ i is the dipole of molecule i, r ij is the distance between molecules i and j, and the Θ-function is unity for positive arguments and zero for negative arguments.This reaction-field approximation is reasonable for a homogeneous system, 55 and the asymptotic large dielectric assumption is isomorphic to the conducting boundary condition used with Ewald sums to construct Panel (a).While, the potential energy computed with a reaction field method is not guaranteed to be the same as that computed by an Ewald summation, judicious choice of cutoff in this instance results in reasonable agreement.Indeed, there is closer correspondence between Figs. 5 (a) and (c) than between Figs.

(a) and (b).
The free energy surfaces at pressures other than those considered in Fig. 5 are easily produced by re-weighting, Eq. 2. Such re-weighted surfaces are shown in Appendix D.

C. Global order contraction 26
Here, we explicitly demonstrate the importance of Q 6 for analyzing phase behavior of supercooled water by showing the effects of controlling the range of accessible Q 6 fluctuations.This discussion follows closely that provided in Ref. 26, and it augments the results of Sec.II.
In particular, we define a contraction of a constrained free energy, We compute these functions from our estimates of the unconstrained reversible free energy surface.Functions so obtained are shown in Fig. 6(a).The unconstrained free energy from which they are derived is F (ρ, Q 6 ) graphed in Fig. 2(a) and re-weighted to the pressure 2.7 kbar.It is the reversible free energy surface for the ST2a model at a pressure that puts the system close to coexistence between the liquid and the crystal at the temperature T = 235 K.At lower pressures, the liquid will be supercooled; at higher pressures, the liquid will be stable with respect to the crystal.The surfaces for those different pressures are obtained from the pictured surface by applying Eq. 2. An upper limit of Q max 6 = 0.13 encompasses the liquid basin.The contracted free energy in that case is unimodal, and it does not exhibit statistically meaningful changes in convexity.That is to say, no re-weighting with Eq. 2 in that case will produce bi-modality.Thus, there is not a second liquid for all densities (and corresponding pressures) in the range considered.
Notice, however, that this contracted free energy function is skewed in a fashion where fluctuations towards low density are more probable than fluctuations towards higher density.In the limit of large system size, fluctuations within stable or metastable basins are Gaussian.The skewed behavior is therefore a finite system-size effect.Its physical origin can be resolved by increasing Q max 6 .For Q max 6 = 0.25, a shoulder and change in convexity appears.For larger values of Q max 6 , there is systematic growth of the shoulder into a basin.This low density basin is the crystal.The mean value of Q 6 as a function of , is shown in Fig. 6(c).This mean value remains at its liquid state value for Q max 6 < 0.45.Thus, by sampling the full surface pictured in Fig. 2(a), it is possible to identify the low-density basin as the crystal phase.
When limiting the range of Q 6 to Q max 6 < 0.45, the features shown in Fig. 6 could be easily misinterpreted as indicative of liquid-liquid coexistence.Indeed, by applying an external pressure to re-weight the curves in Fig. 6 (a), with Eq. 2 as if F (ρ, Q max 6 ) was an equilibrium free energy for 0.13 < Q max 6 < 0.43, a bistable density distribution can be obtained with a low mean value of Q 6 .This type of construction is illustrated in Panel (b) of Fig. 6.The bi-stability of precisely the sort reported in Refs.10, 11, 15, and 19 is thereby found.But contrary to the interpretation expressed in Ref. 19, the behavior is not reflective of liquid-liquid transition.Rather, the unconstrained free energy shows that the appearance of convexity loss is associated with moving towards and then over the barrier separating liquid and crystal basins.
As the geometry of the total F (ρ, Q 6 ) changes with system size, N , the re-weighting (i.e., pseudo Maxwell construction) illustrated in Panel (b) and alluded to in Ref. 19 will depend upon system size in ways that are inconsistent with two-phase coexistence.The free energy barrier separating basins of truly coexisting phases scales as N 2/3 .But the low density phase with Q 6 confined to low values cannot equilibrate and thus cannot coexist.

IV. STUDIES OF ADDITIONAL COMPUTER SIMULATION MODELS
The potential existence of a second liquid phase has been conjectured for a number of different models of water and other tetrahedral liquids.In this section we summarize our work on applying robust free energy calculations to a selection of these models.
Larger phase space for the mW model In Paper I, 2 we reported results for the mW model over a range of temperatures spanning 160 K to 320 K and pressures spanning 1 bar to 3 kbar.A corresponding states argument 39 indicates that this pressure range for mW water is a factor of 3 smaller than that accessed in experiments.As a consequence, the range of pressures relevant to the existence of a liquid-liquid transition in the mW model extends up to 10 kbar.We have thus expanded the domain of our free energy calculations and ruled out a liquid-liquid transition up to 10 kbar and 160 K< T < 320 K.The free energy surfaces throughout are consistent with those shown in Paper I.

Quantitative water model (TIP4P/2005)
The recently parameterized TIP4P/2005 water model 30 has had success in quantitatively reproducing many essential properties of water and ice, including the density versus temperature line at p = 1 bar.Previous studies have extrapolated equation of state data for this model to estimate the location of a putative liquid-liquid critical point and first order transition line. 56Using free energy calculations we can test the validity of this extrapolation.7. The free energy in Fig. 7 is computed below this purported critical point temperature, yet the distribution is clearly monostable at low values of Q 6 .High values of Q 6 were sampled, but not shown for clarity.On re-weighting this free energy, Eq. 2, bi-stability does not appear throughout the pressure range, 1 kbar < p < 3 kbar.We can conclude, therefore, that for this model of water, a second liquid does not exist at conditions studied.Here too, the putative liquid-liquid transition seems to be an artifact of finite-time sampling, reflecting a liquid-to-ice transition where coarsening is incomplete.

General tetrahedral model (SW)
We have also studied the behavior of the Stillinger-Weber (SW) model of silicon.This model has been the subject of numerous studies 6,9,14 that have used equation of state data to propose the existence of a second liquid phase.Using free energy methods, we can test this proposal by examining conditions below the putative liquid-liquid critical temperature.
Figures 8 (a)-(c) shows F (ρ, Q 6 ) calculation for a system of 512 particles at T = 1050 K and a range of pressures, 0 < p < 10 kbar.Based on the phase diagram proposed in Ref.  pressure range should traverse the first-order liquid-liquid phase boundary.Rather than finding two liquid basins upon decreasing pressure, we find the system crosses the line of liquid stability, T s (p), for pressures lower than p = 1 kbar.Put into context with our discussion in Sec.II, the geometry of the free energy in Fig. 8 (a) illustrates how finite-time sampling can produce the illusion of liquidliquid bistability.Specifically, the facile equilibration in density would result in an abrupt change between the liquid with density at ρ = 2.45 g/cc and an amorphous material with ρ = 2.35 g/cc, while over short timescales the slow diffusion in Q 6 would keep the value small.In fact this point was noted by the authors of Ref. 14 when they state that out of 10 to 50 trajectories, "Non-crystallizing samples (an average of 5) were run for up to 10 relaxation times when possible."

V. SUMMARY OF CALCULATIONS ON PUTATIVE LIQUID-LIQUID TRANSITIONS IN SUPERCOOLED TETRAHEDRAL FLUIDS
Table II summarizes all that is now known about the putative liquid-liquid transition in supercooled water and related systems.We see that the low-temperature behaviors of several different models of water-like liquids are similar.The reversible phase behavior in all cases is that of one liquid, a liquid that can coexist with and transform to a lower density ice-like phase.In all cases, the time scales for density fluctuations in the supercooled liquid are several orders of magnitude shorter than those for long-ranged order fluctuations, and both time scales are strongly temperature dependent.
These features lead to a rich non-equilibrium behavior, some of which we have illustrated here in our treatment of the early stages of coarsening, and previously in our treatment of confined supercooled water. 39This non-equilibrium behavior is clearly responsible for the numerous reports of polyamorphism in computer simulations of water, and it is likely important in processes that lead to the formation of glassy phases of water.While the former represent artifacts of finite-time sampling, the behaviors of glasses represent an important class of phenomena worthy of future study.Such far-from equilibrium behavior, however, is beyond the scope of what we have considered here.
The generic nature of what is found in so many models makes it seem unlikely that plausible models of water will exhibit liquid-liquid coexistence and a second lowtemperature critical point.It also suggests that any model exhibiting a crystal phase with lower density than the liquid phase will also exhibit transient behavior that looks like polyamorphism when viewed on the time scales no longer than those of early-stage coarsening.The reason why this non-equilibrium behavior is not observed in a recent study of the SPC/E model 57 is because that study 58 quenches from temperatures significantly higher than T s for that model.
The generic nature also implies that simplified models, like mW water, can reliably serve as useful descriptors of water.Hypotheses on the nature of low temperature water could often be studied in that way.Further, computational efforts that struggle with overcoming the separation of time scales inherent in low-temperature water could often be tested with such models.In our own work (see Appendix C below), it seems clear that errors in prior announcements of liquid-liquid transitions in supercooled water could have been detected by first examining the behavior of mW and SW models.
In documenting the necessity of attending to time scales and relaxation, we have outlined robust methods by which equilibration and reversibility can be achieved.There is clearly need for independent assessment of this growing body of work, which we look forward to seeing in the future.   of the long range part of the potential to 1 part in 10 4 using a standard error estimator and a spherical wavevector cutoff. 61In order to maintain this level of accuracy with a changing box size the parameters are updated over the course of the trajectory.In the unstable trajectory a cubic wave-vector cutoff is used, and the parameters are held fixed which reduces the accuracy of the potential estimate.There may be other ways to induce this pathological behavior, and these ways will depend upon the type of dynamics used.The principal point is that it is a finite-size effect that is sensitive to technical but unphysical details in implementing Ewald sums, and the effect can be avoided when sufficient care is applied to the algorithm.In addition to the Fokker-Plank analysis, we have computed P ne (Q 6 , t) directly from a molecular simulation of the ST2b model, starting from an ensemble of liquid configurations and running for t sim = 10τ Q6 .We have then convoluted that non-equilibrium distribution with the equilibrium P (ρ|Q 6 ) obtained from the reversible free energy in Fig. 5

FIG. 2 .
FIG. 2. Slow relaxation behavior and its consequences for free energy calculations.a) The reversible free energy surface for 216 molecules with the ST2a potential energy function at temperature T = 235 K and pressure p = 2.2 kbar.(See text for definition of the ST2a potential.)Contour lines are separated by 1.5 kBT , and statistical uncertainties are about 1 kBT .b) Negative logarithm of the non-equilibrium distribution for crystal order, Q6, as it relaxes from the liquid state.It is computed from the Fokker-Planck equation with the free energy surface given in Panel (a) under the assumption that the density, ρ, remains at equilibrium with the instantaneous value of Q6.(c) and (d) Non-equilibrium pseudo free energy surfaces computed from Eq. 1 at two intermediate stages of relaxation, t = 10 τQ 6 and t = 1, 000 τQ 6 .The unit of time, τQ 6 , is the autocorrelation time for Q6 fluctuations in the liquid basin (i.e., at small Q6).The reduced density is ρ = (ρ − ρ xtl )/∆ρ, where ρ xtl is the mean density of the crystal basin (i.e., at large Q6), and ∆ρ is the difference between the mean densities of the liquid and crystal basins.Contour lines are separated by 1 kBT and statistical uncertainties are about 1 kBT

FIG. 3 .
FIG.3.Non-equilibrium pseudo free energy surfaces, Fne(ρ, Q6, t), at three different pressures, illustrating how artificial polyamorphism arises as a finite-time effect.All three surfaces are evaluated by propagating from an initial liquid distribution for a time t = 10τQ 6 .Computed as in Fig.2, with the same notation as used in that figure.Contour lines are separated by 1 kBT and statistical uncertainties are about 1 kBT .

Q 6
and time to ensure numerical stability.41Figure 2 (b)

FIG. 5 .
FIG. 5. Free energy surfaces, F (ρ, Q6; p, T ), for three variants of the ST2 model at temperatures where others report evidence of liquid-liquid coexistence for the ST2 model.Possibilities of two-phase coexistence require changes in convexity, as re-weighting through Eq. 2 in that case can produce two basins of equal statistical weight.A coexistence pressure is then the value p + ∆p at which there is equal statistical weight.For the three variants considered, the only changes in convexity are associated with coexistence between a liquid (low Q6) and a crystal (high Q6).a) Free energy for the ST2a variant at T = 230 K and p = 2.2 kbar, with N = 216.b) Free energy for the ST2b variant at T = 230 K and p = 2.2 kbar with N = 216.c) Free energy for the ST2c variant at T = 235 K and p = 2.05 kbar with N = 216.See text for definitions of the different variants.Contour lines are separated by 1.5kBT and statistical errors over the surfaces average to less then 1 kBT .Quantitative features will change with system size.For example, as N grows, the mean value of Q6 in the liquid basin will vanish as 1/N 1/2 , while in the crystal basin it will remain finite.

FIG. 6 . 6 .
FIG. 6. Free energy functions and mean order parameters computed for the ST2a model at T = 235 K, illustrating artificial polyamorphism arising from incomplete knowledge of Q6-dependence in F (ρ, Q6).(a) The restricted contracted free energy function, F (ρ; Q max 6 ), for N = 216 at several indicated choices of Q max 6 .The functions are computed from integrating the reversible free energy function in Fig. 2 (a) re-weighted to pressure p = 2.7 kbar.Error bars indicate one standard deviation.Higher or lower pressures shift the free energy to favor the liquid or crystal, respectively.See Eq. 2. (b) Re-weighted contracted free energy function, with Q max 6 = 0.4, as if the re-weighting coincided with a Maxwell construction for two coexisting liquid phases.(c) The mean value of Q6 as a function of the maximum order-parameter value for T = 235 K and p = 2.7 kbar.

FIG. 7 .
FIG. 7. Reversible free energies and equation of state for liquid TIP4P/2005 at T = 185K and N = 216 for pressures between 1 and 3 kbar.a) Free energy, F (ρ, Q6).Contour lines are separated by 1 kBT and error estimates are less than 1 kBT .b) Contracted free energy as a function of density when restricting Q6 to the liquid basin, i.e., Q6 < 0.14.c) Liquid phase equation of state accessible from the free energy surfaces shown in Panels (a) and (b)

Figure 7
shows F (ρ, Q 6 ) for TIP4P/2005 computed at T = 185K and p = 1.8 kbar for N = 216 molecules.Reference 56 estimates the location of the critical point to be T = 193 K and p = 1.35 kbar, as shown in Panels (b) and (c) of Fig.

FIG. 8 .
FIG. 8. Reversible free energy, F (ρ, Q6), computed for the SW model for N = 512, T = 1050 K and three pressures.Contour lines are separated by 2 kBT and error estimates are less than 1 kBT .
1.8) d (240, 2.0) d (235, 2.2) d (230, 2.4) d (230-240, 1.0-3.0)i ST2a -(230-240, 1.0-3.0)i SW (1070, 0.) f (950, 7.5) f (920, 11.3) f (1050, 0-10.0)j TIP4P/2005 (195, 1.45) g (180-190, 1.0-2.6)i a Conditions where artificial polyamorphism has been reported.b Conditions where Q 6 -equilibrated free energy calculations rule out a liquid-liquid transition.c Short Grand Canonical simulations from Ref. 10, N ≈ 200.d Q 6 -unequilibrated free energy calculations from Ref. 15, N ≈ 200.e Q 6 -unequilibrated free energy calculations from Ref. 19, N ≈ 200 f Short molecular dynamics trajectories from Ref. 14, N = 512.g Short molecular dynamics trajectories from Ref. 56, N = 500.h This paper and Paper I, N = 216 to N = 1000.i This paper, N = 216.j This paper and Paper I, N = 512.version of the manuscript.We have also benefitted from discussion with the authors of Refs.15 and 19.Early work on this project was supported by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division and Chemical Sciences, Geosciences, and Biosciences Division of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.While for the later stages the authors gratefully acknowledge the Helios Solar Energy Research Center, which is supported by the Director, Office of Science, Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

FIG. 10 .
FIG. 10.Nonequilibrium pseudo free energy surfaces, Fne(ρ, Q6, t), at three different pressures, illustrating how artificial polyamorphism arises as a finite-time effect.All three surfaces are evaluated by propagating from an initial liquid distribution for a time t = 10τQ 6 .Computed for the ST2b model from the Fokker-Planck analysis with the underlying reversible free energy surface shown in Fig. 5b.The appearance of artificial polyamorphism is like that found in Ref. 15 for similar conditions.The time t = 10 τQ 6 corresponds to the time used for averaging in Ref. 15. Contour lines are separated by 1 kBT and statistical uncertainties are less than 1 kBT .The correct, reversible free energy surface for p = 2.2 kbar is shown in Fig. 5(b).

FIG. 11 .
FIG. 11.Nonequilibrium pseudo free energy surfaces, Fne(ρ, Q6, t) illustrating how artificial polyamorphism arises as a finite-time effect.The correct reversible surface is shown in Fig.5(b).The incorrect surface shown here is computed for the ST2b model from a non equilibrium distribution for Q6, where the distribution is obtained from simulations initiated in the liquid basin and propagated for a time 10 τQ 6 .The surface shows the appearance of artificial polyamorphism like that found in Ref. 15 for similar conditions.The simulation time tsim = 10 τQ 6 corresponds to the time used for averaging in Ref.15.Contour lines are separated by 1 kBT and statistical uncertainties are less than 1 kBT .
Appendix B: Early stage coarsening and artificial polyamorphism in the ST2b model The calculations carried out in Section II for the ST2a model have been also done for the ST2b model at T = 230 K with N = 216 and pressures ranging from 2 kbar to 2.4 kbar.These are the conditions and model considered in Ref. 15.The results are shown in Fig. 10.The theoretical results agree with those found in Ref. 15, thus indicating that the free energies reported in that work suffer from finite-time effects.
(b) according to Eq. 1.The results are shown in Fig. 11.Again, behavior like that reported in Ref. 15 is obtained, thus indicating that the results of that paper suffer from finite-time effects.

TABLE I
25timated for T = 220K and p = 1 bar from analysis in Ref.39.Estimated from the critical cooling rate of 10 6 K/s needed to form amorphous ice.25 g 14, also calculated with 512 particles and citing agreement with Ref. 6, this temperature and