Finite-size scaling investigation of the liquid-liquid critical point in ST2 water and its stability with respect to crystallization

The liquid-liquid critical point scenario of water hypothesizes the existence of two metastable liquid phases---low-density liquid (LDL) and high-density liquid (HDL)---deep within the supercooled region. The hypothesis originates from computer simulations of the ST2 water model, but the stability of the LDL phase with respect to the crystal is still being debated. We simulate supercooled ST2 water at constant pressure, constant temperature and constant number of molecules N for N<=729 and times up to 1000 ns. We observe clear differences between the two liquids, both structural and dynamical. Using several methods, including finite-size scaling, we confirm the presence of a liquid-liquid phase transition ending in a critical point. We find that the LDL is stable with respect to the crystal in 98% of our runs (we perform 372 runs for LDL or LDL-like states), and in 100% of our runs for the two largest system sizes (N=512 and 729, for which we perform 136 runs for LDL or LDL-like states). In all these runs tiny crystallites grow and then melt within 1000 ns. Only for N<=343 we observe six events (over 236 runs for LDL or LDL-like states) of spontaneous crystallization after crystallites reach an estimated critical size of about 70+/-10 molecules.


I. INTRODUCTION
For many centuries, water and its anomalies have been of much interest to scientists. A particular rise of interest occurred in the late 1970s after experiments done by Angell and Speedy seemed to imply some kind of critical phenomenon in supercooled liquid water at very low temperatures [1][2][3][4]. Even though liquid water experiments are limited by spontaneous crystallization below the homogenous nucleation temperature (T H ≈ 233 K at 1 bar), it is possible to further explore the phase diagram by quenching water to far lower temperatures [5][6][7]. The result of these experiments is an amorphous solid, i.e. a glassy ice, corresponding to an out-of equilibrium state that is very stable with respect to the equilibrium crystalline ice phase. The amorphous depends on the applied pressure: at low pressure, below ≈ 0.2 GPa, the low density amorphous ice (LDA) is formed, while at higher pressure the high density amorphous ice (HDA) is observed [8]. It has been shown by Mishima et al. that these two amorphous ices are separated by a reversible abrupt change in density that resembles in all its respects an equilibrium first order phase transition [9][10][11][12].
Raising the temperature of either LDA or HDA does not turn the sample into a liquid, but leads once again to spontaneous crystallization (around T X ≈ 150 K). In fact, between T H and T X , often called the "no man's land" of bulk water, crystallization occurs at a time scale that is too short for current experimental methods, although a new technique is possibly succeeding in the task of measuring the metastable liquid phase [13].
Computer simulations of water, however, involve time scales small enough to witness spontaneous crystallization and are therefore able to explore liquid water in the "no man's land". In 1992 Poole et al. [14] performed a series of molecular dynamics simulations using the ST2 water model [15], using the reaction field method for the long-range interactions (ST2-RF), and discovered a liquid-liquid phase transition ending in a critical point, separating a low density liquid (LDL) and a high density liquid (HDL). These two liquids can be considered to be the liquid counterparts of the LDA and HDA, respectively.
Recently Limmer and Chandler used Monte Carlo umbrella sampling to investigate the ST2-Ew model, but claimed to have found only one liquid metastable phase (HDL) rather than two [63]. They therefore concluded that LDL does not exist because it is unstable with respect to either the crystal or the HDL phase. The emphasis in their work is about the difference between a metastable phase, i.e. separated from the stable phase by a finite free-energy barrier, and an unstable state, where the free-energy barrier is absent and the state does not belong to a different phase.
Shortly after, Poole et al. [64] and Kesselring et al. [65] presented results using standard molecular dynamics for ST2-RF showing the occurrence of the LLCP with both HDL and LDL phases metastable with respect to the crystal, but with the LDL not unstable with respect to either the crystal or the HDL. This result was confirmed, using the same method as Limmer and Chandler, by Sciortino et al. [66] and Poole et al. [67] in ST2-RF water and by Liu el al. in ST2-Ew water [68].
The aim of this paper is to confirm the presence of a liquid-liquid critical point in water in the thermodynamic limit using finite size scaling techniques, and confirm that LDL is a bona fide metastable liquid. We use the ST2-RF model because it has been well-studied in the supercooled region, making it easier to compare and verify our data. In the supercooled phase it has a relatively large self-diffusion compared to other water models, therefore suffers less from the slowing down of the dynamics at extremely low temperatures. We explore a large region of the phase diagram of supercooled liquid ST2-RF water ( Fig. 1) using molecular dynamics simulations with four different system sizes by keeping constant the number N of molecules, the pressure P and the temperature T (N P T ensemble).
Within the explored region we find both LDL and HDL, separated at high pressures by a LLPT, ending in a LLCP estimated at P C ≈ 208 MPa and T C ≈ 246 K. This phase transition is particularly clear in Fig. 2 where one can see from the density how the system continuously flips between the two states. However, due to finite size effects this phase flipping also occurs below the critical point along the Widom line (the locus of correlation length maxima) [69,70]. For this reason it is necessary to apply finite size scaling methods to establish the exact location of the critical point.
For six state points and for small system size N ≤ 343 we observe, in only one over the (on average) seven simulations we performed for each state point, irreversible crystal growth, indicated as red stars in Fig. 1. Each of these crystallization events occurred within the LDL (or LDL-like) region. Analysis of these crystals revealed them to have a diamond cubic crystal structure. As we will discuss later, because these events disappears for larger systems, we ascribe these crystallization to finitesize effects.
We start in Sec. II with a description of the model this fact is subsequently used in Sec. VI where we discuss the growth and melting of crystals within the LDL liquid. In Sec. VII, by defining the appropriate order parameter, we show that the LLCP in ST2-RF belongs to the same universality class as the 3D Ising model. We accurately determine where the LLCP is located in the phase diagram in the thermodynamic limit by applying finite size scaling on the Challa-Landau-Binder parameter. We discuss our results and present our conclusions in Sec. VIII.

II. SIMULATION DETAILS
In the ST2 model [15] each water molecule is represented by a rigid tetrahedral structure of five particles. The central particle carries no charge and represents the oxygen atom of water. It interacts with all other oxygen atoms via a Lennard-Jones (LJ) potential, U LJ (r ij ) ≡ 4ε (σ/r ij ) 12 − (σ/r ij ) 6 with ε ≡ 0.31694 kJ/mol and σ ≡ 3.10Å. Two of the outer particles represent the hydrogen atoms. Each of them carries a charge of +0.2357 e, and is located a distance 1Å away from the central oxygen atom. The two remaining particles carry a negative charge of −0.2357 e, are positioned 0.8Å from the oxygen, and represent the lone pairs of a water molecule.
The electrostatic potential in ST2 is treated in a special way. To prevent charges a and b from overlapping, the Coulomb potential is reduced to zero at small distances: where S(r ij ) is a function that smoothly changes from one to zero as the distance between the molecules decreases, with R L ≡ 2.0160Å, R U ≡ 3.1287Å, and where r ij is the distance between the oxygen atoms of the interacting molecules. In the original model a simple cutoff was used for the electrostatic interactions. In this paper, however, we apply the reaction field method [71] which changes the ST2 Coulomb potential to where T (r ij ) is another smoothing function: We use a reaction field cutoff R c ≡ 7.8Å together with R T ≡ 0.95R c . These parameters define our ST2-RF water model and are the same that were used in previous ST2-RF simulations.
For the LJ interaction we use a simple cutoff at the same distance of 7.8Å. We do not adjust the pressure to correct for the effects of the LJ cutoff [72,73], since these adjustments come from mean field calculations which become increasingly weak as one approaches a critical point.
We use the SHAKE algorithm [74] to keep the relative position of each particle within a ST2 molecule fixed. The temperature and pressure are held constant using a Nosé-Hoover thermostat [73,75,76] together with a Berendsen barostat [77]. In all simulations periodic boundary conditions are applied.
Our code is validated by simulating the same state points as those published by Poole et al., see Fig. 1b in [56], where pressure corrections for the LJ cutoff were applied in the N V T (constant N , T and volume V ) ensemble. Averaging at each state point over 10 simulations with different initial conditions allows us to estimate the error bars. In Fig. 3 we compare our results for N = 216 molecules and density 0.83 g/cm 3 , and find that our data, after pressure correction, matches that of Ref. [56] well.
For each of the simulations done in the N P T ensemble, we use the following protocol. We first create a box of N molecules at n different initial densities (with n up to 21) ranging from 0.85 to 1.05 g/cm 3 . We then perform a 1 ns N V T simulation at T = 300 K. In this way we obtain n independent configurations all at T = 300 K in the prefixed range of densities. Next, we use these independent configurations as starting points for N P T To validate our code, we compare our simulation results with those from Poole et al. [56] at density ρ = 0.83 g/cm 3 and for N = 216 molecules. We performed simualtions in the N V T ensemble applying pressure corrections and find the same results as Ref. [56] within the error bars. At this density the pressure correction due to the LJ cutoff (proportional to ρ 2 ) is equal to −12.66 MPa. The variation of P with T along this isochore shows the occurence of both a density maximum at 300 K and a density minimum near 265 K, as at these state points (∂ρ/∂T )P = −ρKT (∂P/∂T )V = 0 with KT > 0 the isothermal compressibility.
simulations at T = 265 K and different pressures ranging from 190 to 240 MPa, and continue the simulation for an additional 1 ns. This results in n independent configurations at T = 265 K and the given pressure. For all pressures considered here, this will lead the system into the HDL phase. Finally the system is quenched to the desired temperature at the given pressure, followed by 100-200 ns of equilibration time. In Sec. IV it will be shown that this provides enough time for the system to reach equilibrium for the state points above the line marked with the label T g in Fig. 1

III. INTERMEDIATE SCATTERING FUNCTION
The intermediate scattering function S(k, t) plays an essential role in the analysis of liquid structure, since it is frequently measured in experiments as well as easily calculated from simulation data. It describes the time evolution of the spatial correlation at the wave vector k, and can be used to distinguish between phases of different structure, such as LDL and HDL or crystal. It is defined as where ... t denotes averaging over simulation time t , and r (t ) the position of particle at time t . For sim-plicity we only apply the intermediate scattering function to the oxygen atoms, which we denote as S OO (k, t).
Since the system has periodic boundary conditions, the components of k have discrete values 2πn/L, where L is the length of the simulation box and n = 1, 2, 3, . . . . We define S OO (k, t) ≡ S OO (k, t) n where the average is taken over all vectors k with magnitude k belonging to the nth spherical bin π(n − 1 2 )/L ≤ k < π(n + 1 2 )/L for n = 2, 3, . . . , 300. Similarly, we define the structure factor S OO (k) ≡ S OO (k, t) t as the time-averaged intermediate scattering function, with (unless indicated otherwise) the average taken over the whole duration of the run. We study S OO (k) above and below our estimate for the LLCP pressure. At P = 210 MPa > P C (Fig. 4a,b) we observe a discontinuous change in the first two peaks of S OO (k) as T changes between 245 and 246 K, and a continuous change above and below this temperatures. This is the expected behavior for a first order phase transition occurring at 245 K < ∼ T < ∼ 246 K and P = 210 MPa between two phases with different structure, consistent with our results in Fig. 1. The fact that for both phases S OO (k) ∼ O(1) for all k shows that both phases are fluid. Indeed, for a crystal-like configuration, with a long-range order, there would be at least one wave vector such that S OO (k) ∼ O(N ) [78]. Furthermore, the fact that at lower T the first peak increases and the other peaks only have minor changes indicates that the lower-T liquid has a smaller density than the higher-T liquid. Therefore, this result show a first-order phase transition between the LDL at lower-T and HDL at higher-T . This transition occurs at the same temperature at which we observe the phase flipping in density (Fig. 2) and corresponds to the purple region at P > P C in Fig. 1.
The fact that the peaks of S OO (k) are sharper in LDL than HDL is an indication that the LDL phase is more structured. We can also observe that the major structural changes in S OO (k) between LDL and HDL are for k 1.8 and 2.8Å −1 , corresponding to r = 4π/k 7 and 4.5Å, respectively, i.e. are for the third and the second neighbor water molecules. This change in the structure is consistent with a marked shift inwards of the second shell of water with increased density, and almost no change in the first shell (at k 4.6Å −1 and r 2.75Å), as seen in structural experimental data for supercooled heavy water interpreted with Reverse Monte Carlo method [79]. This changes are visible also in the OO radial distribution function g OO (r) (Fig. 5a,b). For P < P C (Fig. 4c,d) by increasing T we observe that the first peak of S OO (k) merges with the second, transforming continuously in a shoulder. Same qualitative behavior is observed for g OO (r) (Fig. 5c,d). These quantities show us also that the lower-T structure is LDLlike, while the higher-T structure is HDL-like. However, the absence of any discontinuous change in the structure implies the absence of a first-order phase transition in the structure of the liquid. This is consistent with the occurrence of a LLCP at the end of the first-order phase transition somewhere between 200 and 210 MPa, at a temperature between 245 and 250 K. In Sec. VII we shall apply a different method to locate the LLCP with more precision.
At P < P C , in the one-phase region, we expect to find the Widom line emanating from the LLCP. The Widom line is by definition the locus of maxima of the correlation length, therefore, for general thermodynamic considerations [70] near the LLCP it must be also the locus of maxima of the response functions. In particular, it must be the locus where the isobaric heat capacity C P ≡ T (∂S/∂T ) P , where S is the entropy of the system, has its maximum along a constant-P path. This maximum occurs where the entropy variation with T is maximum, expected where the structural variation of the liquid is maximum, i.e. where the derivatives of the values of S OO (k) (Fig. 4d) and g OO (r) (Fig. 5d) with T are maximum. The interval of temperatures for each P where this occurs corresponds to the purple region at P < P C in Fig. 1, indicated as the Widom line.
It is actually possible to follow the structural changes during the simulation. An example is given in Fig. 6 where we focus on a 30 ns time period of a simulation at 200 MPa and 248 K. We divide this time period into six 5 ns intervals and for each interval we calculate the intermediate scattering function, time-averaged over those 5 ns. We observe that the liquid is LDL-like for the first and third interval, having low density and LDL-like S OO (k) (first peak near 2Å −1 , separated from the second). On the contrary, for the fifth and sixth interval the density is high and S OO (k) is HDL-like (the first peak is merely a shoulder of the second peak), indicating that the liquid is HDL-like. For the second and fourth interval, the liquid has an intermediate values of density and S OO (k), indicating that it is a mix of LDL-like and HDL-like structures.

IV. CORRELATION TIME
Apart from its use in structure analysis, the intermediate scattering function S OO (k, t) can also be used to define a correlation time τ , i.e. the time it takes for a system to lose most of its memory about its initial configuration [80,81].
In Fig. 7 we show how S OO (k, t) decays with time for a fixed value of k. Its decay is characterized by two relaxation times, the α-relaxation time τ α and the β-relaxation time τ β . On very short time scales, the molecules do not move around much and each molecule is essentially stuck in a cage formed by its neighbors. The β-relaxation time τ β is of the order of picoseconds. On longer time scales, the molecule can escape from its cage and diffuse away from its initial position. The time τ α is the relaxation time of this structural process.
Mode-coupling theory of supercooled simple liquids predicts that [82] F OO (k, t) ≡ S OO (k, t)/S OO (k, 0) The factor A(k) is the Debye-Waller factor arising from the cage effect, which is independent of the temperature and follows A(k) = exp(−a 2 k 2 /3) with a the radius of the cage. We are able to fit Eq. (5) remarkably well to all our data, as for example in Fig. 7. Data in Fig. 7 was collected every 10 fs for simulations of 1 ns. This rate of sampling results in a large amounts of data and is unfeasible for our runs up to 1000 ns. Therefore, for the 1000 ns runs we collect data at 10 ps intervals. At this rate of sampling it is no longer possible to estimate τ β or the cage size a, but it is still possible to determine τ α accurately, utilizing the fact that S OO (k, t) reaches a plateau near t ≈ τ β . One can therefore define which is S(k, t) normalized by its value at the plateau (Fig. 8). A good estimate of τ α is then the time for which C OO (k, τ α ) = 1/e ≈ 0.37. From the shorter 1 ns runs (which were mostly done in the HDL regime) we find that the cage radius is a = 0.35 ± 0.09Å with a stretching exponent of b = 0.63 ± 0.09. Both parameters a and b do not show a significant dependence on the state point within the studied range of temperature and pressure.
As shown in Fig. 7, different k result in slightly different values for τ α . We use as the correlation time τ the largest value of τ α which is usually found at k = k 1 , the first maximum in S OO (k) (inset Fig. 7). As is to be expected, the correlation time does not seem to depend on the box size (Fig. 9). It does however depend strongly on the phase, which is evident from Fig. 10. At high temperatures the system is in the HDL phase, and has a correlation time τ on the order of 10-100 ps. As we decrease the temperature at fixed pressure, the value of τ has a large increase when we cross the phase transition line or the Widom line, depending if P is above or below P C , respectively. Apparently, the LDL states evolve nearly four orders of magnitude slower than HDL states, with correlation times in the nanosecond range.
If we lower the temperature further, the correlation time slowly increases until the system becomes a glass rather than a liquid, and we are no longer able to fully equilibrate the system. As we can only run simulations up to 1000 ns, we consider the state points with a correlation time above 100 ns to be beyond our reach. We therefore designate the effective glass transition temperature T g as the temperature for which τ > 100 ns (see Fig. 1).

V. STRUCTURAL PARAMETERS
Apart from the intermediate scattering function, there are other ways to quantify the structure of a liquid. In this section we shall examine several structural parameters, and determine which of those are the most effective in distinguishing between LDL, HDL, and the crystal. For simplicity, we approximate the center of mass of a water molecule with the center of its oxygen atom.
The structural parameters are designed to distinguish between different phases by analyzing the geometrical structure. This is typically done by evaluating the spherical harmonics Y m (ϕ, ϑ) for a particular set of neighboring atoms, with ϕ and ϑ the polar angles between each pair of oxygen atoms in that set. In this paper we consider two different sets: we define the first coordination shell n 1 (i) to be the four nearest neighbors of molecule i, and define the second coordination shell n 2 (i) as the fifth to sixteenth nearest neighbors (the sixteenth nearest neighbors minus those in the first shell).
Different values of are sensitive to different symmetries. The spherical harmonics with = 3, for example, are sensitive to a diamond structure. Those with = 6 are more sensitive to the hexagonal closest packing (hcp) structure. Since we expect the liquid and crystal structures to be hcp, diamond, or a mix of these, we focus primarily on = 3 and = 6.

A. Parameters q3 and q6
All parameters defined in this section are based on q (s) ,m (i) which quantifies the local symmetry around molecule i. It is defined as where and m are integers, s = 1, 2 indicates the shell we are considering, with N s the number of molecules within that shell (i.e. N 1 ≡ 4 for the first coordination shell, and N 2 ≡ 12 for the second). Y m is normalized according to |Y m | 2 sin(ϑ)dϕdϑ = 1. We can consider q This means that we can define an inner product ,m (j)) + Im(q (s) ,m (i)) Im(q (s) ,m (j)) (8) and a magnitude The local parameter q (s) (i) is one way to distinguish between different structures, and can be used to label individual molecules as LDL-like or HDL-like. We can convert it into a global parameter by averaging over all molecules, In Fig. 11 we see that all global q (s) are sensitive to the difference between LDL and HDL, especially q (1) 3 and q (2) 6 . We conclude that the structural difference is visible in both the first and second shell, and that LDL and HDL differ mostly in the amount of diamond structure of the first shell and the amount of hcp structure in the second shell. This is confirmed by the histograms in Fig. 12, in which the largest difference between LDL and HDL is seen in q (1) 3 and, next, in q (2) 6 . The latter is the parameter that better discriminate with respect to the crystal structure.

B. Global parameters Q3 and Q6
An alternative approach, as used by Steinhardt et al. [83], is to first average q (s) ,m (i) over all molecules, defining ,m (i), then calculate the magnitude Our calculations show that the parameters Q and Q (s) 6 , with s = 1, 2, are not efficient in discriminating between LDL and HDL (Fig. 11), although Q 6 ≡ Q (1) 6 has been proposed recently as a good parameter to this goal [63] and consequently has been used by several authors [66][67][68]. In particular, we observe that there is not much correlation between the fluctuations of Q (s) and those of the density, except for Q 3 . However, we confirm that Q   approximately 10 times larger for the crystal than it is for the liquids (Fig. 13). This large increase of Q (s) 6 for crystal-like structures might be related to the few instances in Fig. 11 where an increase in Q (s) 6 corresponds to a decrease of density (such as within interval t = 230-237 ns), consistent with the observation that the crystal-like structures have a density comparable to the LDL structure and smaller than the HDL structure.
To confirm that LDL remains a liquid in the thermodynamic limit, we look at how Q 6 changes with the system size. For liquids Q 6 scales like N −1/2 while for crystals the value Q 6 remains finite as N → ∞. We find that the probability distribution functions of Q 6 N 1/2 for N = 216, 343, 512, and 729 overlap, which means that Q 6 N 1/2 is independent of the system size, therefore Q 6 ∼ N −1/2 (Fig. 14). We conclude that the metastable LDL is not transforming into the stable crystal in the thermodynamic limit. This implies that the LDL and the crystal phase are separate by a free-energy barrier that is higher than k B T at the temperatures we consider here and that the system equilibrates to the stable (crys-  tal) phase only in a time scale that is infinite with respect to our simulation time (1000 ns), as occur in experiments for metastable phases. Therefore, the LDL is a bona fide metastable state. Our conclusion is consistent with recent calculations by other authors [66][67][68].

C. Bond parameters d3 and ψ3
We define the bond order parameter d (s) similar to that defined by Ghiringhelli et al. in Ref. [84], where the quantity d (1) 3 (i, j) characterizes the bond between molecules i and j, and is designed to distinguish between a fluid and a diamond structure. The local parameter d (s) (i, j) is defined as the cosine of the angle between the vectors q (s) (i) and q (s) (j): with the inner product and magnitude as defined in Eqs. (8) and (9). A crystal with a perfect diamond structure has d the bonds within the same layer (three out of four) have d (1) 3 (i, j) = −1, while the bonds connecting atoms in different layers (one out of four) have d (1) 3 (i, j) = −1/9. We find that the parameters d (s) for = 3, 6 and s = 1, 2 do not distinguish well between the two different liquidlike structures, but that d for both s = 1 and 2 are suitable to discriminate between the crystallike structure and the liquids (Fig. 15). In particular, for the crystal, most molecules have d (1) 3 < −0.87, and we therefore consider a molecule to be part of a crystal if at least three out of its four bonds with its nearest neighbors have d (1) 3 < −0.87. This is the same cutoff used by Ghiringhelli et al. in [84].
The global parameter associated to d (s) (i, j) is defined as where FIG. 14: Finite size scaling of parameter Q6 in the LDL phase (210 MPa, 243 K). The probability distribution function of Q6N 1/2 is independent of the system size N , which means LDL scales like a liquid in the thermodynamic limit: is the average of d (s) (i, j) over the first four nearest neighbors of the molecule i. We observe that each ψ (s) (i) has the same features of the corresponding 3 (i) discriminating well between the crystal-like and the liquids-like structures (Fig. 16). We observe that ψ (1) 3 discriminates well between LDL-like and HDL-like structures (Fig. 11), while ψ (s) 6 for s = 1 and 2 might be able to emphasize the temporary appearance of crystal-like structures, as noted for Q (s) 6 .

VI. GROWTH AND MELTING OF CRYSTAL NUCLEI
In a small percentage of our simulations, the system was found to spontaneously crystallize. These are interesting events because spontaneous crystallization of water in molecular dynamics is extremely rare; only recently Matsumoto et al. were the first to successfully simulate the freezing of water on a computer [85]. Crystallization events in supercooled ST2 water are particularly important to study, as it has been proposed that LDL is unstable against crystallization [63].
Following the discussion in Sec. V, we define a crystal as a cluster of molecules which has three out of four bonds with d 3 < −0.87 and belong to the first coordination shell of each other. In this section we shall study the growth and melting of these crystal nuclei, and estimate the critical nucleus size needed to overcome the free energy barrier. The existence of this barrier allows us to conclude that LDL is in fact a bona fide metastable state with respect to the crystal.
In Fig. 17 we show the density evolution for 11 different configurations, each with 343 molecules and at 205 MPa and 246 K. Each of these runs started at a dif- in (a), these parameters do not distinguish well between the two different liquid-like structures, but d ferent initial density (between 0.85 and 0.95 g/cm 3 ) and was subsequently equilibrated to the final temperature and pressure using the procedure described in Sec. II. Because this state point lies close to the LLPT, we see phase flipping in all of them. However, the two configurations C and F display a sudden jump to a stable low density plateau. This is a hallmark of crystalization. We confirm this by calculating the size of the largest crystal as a function of time (Fig. 18). During most runs the largest crystal continuously grows and shrinks, but never reaches a size larger than 30 molecules. On the other hand, configurations C and F show a jump in crystal size exactly matching the jump in density. Run F ends up partially crystalized, while for C we find that over 90% of the box is crystallized in a diamond structure with a density of about 0.92 g/cm 3 (Fig. 19).
The correlation time increases dramatically if crystals appear with a size comparable to the system size, as is evident from Fig. 20. The correlation functions of C and F decay very slowly, leading to correlation times of 200-400 ns, while the other configurations have a correlation time of less than 4 ns.
For spontaneous crystallization to occur, a sufficiently large crystal nucleus needs to form within the liquid. According to classical nucleation theory, this nucleus needs to reach a minimum size to prevent it from melting. We observed in many simulations that a small nucleus grows and melts, and a few runs in which the nucleus grows further or remains stable. Therefore, we can make an estimate of the critical nucleus size. The two largest crystals that formed and subsequently melted, both reached a size of about 50-60 molecules ( Fig. 21a and 21b). The smallest crystal that formed and remained stable, had a size of about 50-80 molecules (Fig. 21c). We therefore conclude that the critical nucleus size is approximately 70 ± 10 molecules. A similar value of 85 molecules was found by Reinhardt and Doye [86] for ice nucleation in the monatomic water model [87].
For a more accurate estimate it is necessary to run longer simulations, as the crystal nuclei can survive for hundreds of nanoseconds (e.g., Fig. 21d in which a small crystal lasts for 700 ns).

VII. LOCATING THE CRITICAL POINT
In Sec. III we used the intermediate scattering function S OO (k) to estimate the position of the liquid-liquid critical point, and found it to lie near 200-210 MPa and 244-247 K. It is commonly believed that the LLCP falls in the same universality class as the three-dimensional Ising model [57]. At the critical point the order parameter distribution function (OPDF) of a system has the same bimodal shape as all other systems that belong to the same universality class. Therefore we can locate the LLCP accurately by fitting our data to the OPDF of the 3D Ising model (Fig. 22).

FIG. 22:
Fitting the order parameter distribution function (obtained from the simulation data) to that of the 3D Ising model at the critical point (black curve, from Ref. [88]). The shape of the order parameter distribution function depends on temperature and pressure. Here (for N = 343) we find an excellent fit for PC = 206 MPa and TC = 246 K, and we therefore confirm that the LLCP indeed belongs to the same universality class as the 3D Ising model. Based on our fit, we locate the LLCP to be at PC = 206 ± 3 MPa and TC = 246 ± 1 K.
In the 3D Ising model the order parameter M is simply the spontaneous magnetization, but for liquids the order parameter turns out to be a linear combination of two independent quantities such as the density and the potential energy [89,90]. We therefore define M ≡ ρ+sE, with s a constant known as the field mixing parameter. The value of s depends only on the model and should therefore be independent of the number of molecules. Our fits indeed confirm this; we find s = 0.0362 (g/cm 3 )/(kJ/mol) for all values of N .
Only the shape of the OPDF is dictated by the theory, which means we are free to move and stretch our OPDF to acquire an accurate fit. So, instead of fitting the order parameter M , we actually fit x ≡ A(M − M C ) to the 3D Ising model. The critical order parameter M C is chosen such that the mean value of x is zero, and the amplitude A has been chosen such that the variance equals unity.
To calculate the OPDF for a particular pressure and temperature, we create a two-dimensional histogram of the density and energy. An example of such a 2D histogram is shown in Fig. 23  To fit our OPDF to that of the 3D Ising model, we need to calculate our OPDF at different pressures and temperatures, until we find the (P C , T C ) that gives us the best fit. The state point (P C , T C ) is then our best estimate of the location of the LLCP. We only have simulation data for a finite number of state points, therefore some kind of interpolation is necessary. The method of choice here is histogram reweighting [91]; we use the algorithm as described by Panagiotopoulos in Ref. [92].
The results of fitting our data to the 3D Ising model are shown in Fig. 24. Tables I and II indicate which data was used by the histogram reweighting method to obtain these fits. For N = 343, 512, and 729, we are able to fit our data very accurately to the OPDF of the 3D Ising model, and find the critical point to be located at T C = 246 ± 1 K, P C = 206 ± 3 MPa for N = 343, and at T C = 246 ± 1 K, P C = 208 ± 3 MPa for N = 529 and 729. Theory predicts that the location of the critical point depends on N , and these findings agree with that prediction. In particular, the 3D Ising model predicts that the amplitude A should scale with box size L as A ∼ L β/ν ∝ N β/3ν with β/ν = 0.52 [89,93], in agreement with the slope of A(N ) in Fig. 25. This figure   FIG. 24: Best OPDF fits for N = 343, 512, and 729 molecules.
In all cases we obtain s = 0.0362 (g/cm 3 )/(kJ/mol) for the field mixing parameter. We find the critical point to be located at TC = 246 ± 1 K, PC = 206 ± 3 MPa for N = 343, and at TC = 246 ± 1 K, PC = 208 ± 3 MPa for N = 529 and 729. also indicates that N = 216 cannot provide an accurate estimate of the location of the LLCP.
the Widom line, while there remain two at the phase transition line. Therefore, Π → 2/3 at the Widom line, while Π < 2/3 at the LLPT even in the limit N → ∞. Hence, the finite-size scaling of Π allows us to distinguish whether an isobar crosses the LLPT or the Widom line, and is yet another method of estimating the location of the critical point.
We study Π versus temperature T and system size N for different pressures, finding minima Π min at specific temperature for each pressure (Fig. 26). The finite-size dependence of Π min (P ) reveals if P < P C or P > P C (Fig. 27).
For P < P C the mimimum Π min approaches 2/3 linearly with 1/N , while for P ≤ P C it approaches the limit [94] This limiting value is also approached linearly with 1/N . Here ρ H ≡ ρ H (P ) and ρ L ≡ ρ L (P ) are the densities of the two phases LDL and HDL [96]. Above the critical pressure the limiting value of Π min decreases as P increases, i.e. the two peaks of the bimodal D(ρ) move further apart. This happens because ρ H − ρ L increases at coexistence as (P − P C ) β where β ≈ 0.3 is the critical exponent of the 3D Ising universality class [100,101]. From this analysis (Fig. 27) we conclude that our results agree with theory and that the critical pressure P C ≈ 190-210 MPa, in agreement with the estimate of Sec. VII. Furthermore, as Π remains less than 2/3 for P > P C even in the limit N → ∞, we conclude that the LLPT does not vanish in the thermodynamic limit.

VIII. CONCLUSIONS
We performed molecular dynamic simulations in the N P T ensemble for ST2-RF water in the supercooled region of the phase diagram for different system sizes with simulation times of up to 1000 ns. Using several different techniques we confirmed the existence of two liquid phases, LDL and HDL, separated by a liquid-liquid phase transition line. Near the LLPT line the system continuously flips between the two phases. Because of finite size effects this phenomenon also occurs near the Widom line, but by fitting the order parameter distribution function to that of the 3D Ising model, we were able to accurately FIG. 27: Minima of the Challa-Landau-Binder parameter Π as a function of system size N for different pressures. The minimum Πmin occurs at the pressures and temperatures of the LLPT and the Widom line, and is always less than 2/3 for a finite system because of the bimodality of the density histogram. As N → ∞ the bimodality disappears in the onephase region but remains at the LLPT, and therefore Πmin → 2/3 at the Widom line while Πmin < 2/3 on the LLPT, even in the thermodynamic limit. We conclude that the critical point survives in the thermodynamic limit, and that it is located between P = 200 and 210 MPa (in agreement with previous results of this paper).
determine the location of the liquid-liquid critical point (at T C = 246 ± 1 K, P C = 208 ± 3 MPa). Finite size scaling of the Challa-Landau-Binder parameter indicates that the critical point does not disappear in the thermodynamic limit.
Both phases have been confirmed to be bona fide metastable liquids that differ substantially in structural as well as dynamical properties. It is found that the LDL phase is a more "structured" liquid, and that it has a correlation time of almost four orders of magnitude larger than that of HDL, with LDL correlation time of the order of 100-1000 ns. We show that Q 6 structural parameter is not able to discriminate between HDL and LDL, but can discriminate well between liquids and crystal. Finite size scaling of the Q 6 parameter confirms that LDL scales as a liquid and not as a crystal.
The different structures of LDL and HDL are better discriminated by structural parameters like q (1) 3 and q (2) 6 . These parameters show that LDL and HDL differ mostly in the amount of diamond structure of the first shell and the amount of hcp structure in the second shell.
For small box sizes (N = 343) there were a few simulation runs that resulted in spontaneous crystallization, always within the LDL region of the phase diagram. Further analysis revealed that during all simulations small crystals grow and melt within the liquid, a clear indication that LDL is metastable with respect to the crystal. From the few crystalization events that occurred, we were able to conclude that the critical nucleus size is approximately 70 ± 10 molecules.