Correlation between thermodynamic anomalies and pathways of ice nucleation in supercooled water
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n Correlation between thermodynamic anomalies and pathways of ice nucleation insupercooled water
Rakesh S. Singh and Biman Bagchi ∗ Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore 560012, India (Dated: April 11, 2018)The well-known classical nucleation theory (CNT) for the free energy barrier towards formationof a nucleus of critical size of the new stable phase within the parent metastable phase fails to takeinto account the influence of other metastable phases having density/order intermediate betweenthe parent metastable phase and the final stable phase. This lacuna can be more serious thancapillary approximation or spherical shape assumption made in CNT. This issue is particularlysignificant in ice nucleation because liquid water shows rich phase diagram consisting of two (highand low density) liquid phases in supercooled state. The explanations of thermodynamic and dy-namic anomalies of supercooled water often invoke the possible influence of a liquid-liquid transitionbetween two metastable liquid phases. To investigate both the role of thermodynamic anomaliesand presence of distinct metastable liquid phases in supercooled water on ice nucleation, we employdensity functional theoretical approach to find nucleation free energy barrier in different regions ofphase diagram. The theory makes a number of striking predictions, such as a dramatic lowering ofnucleation barrier due to presence of a metastable intermediate phase and crossover in the depen-dence of free energy barrier on temperature near liquid-liquid critical point. These predictions canbe tested by computer simulations as well as by controlled experiments.
I. INTRODUCTION
Liquid water shows pronounced thermodynamic(isothermal compressibility, isobaric heat capacity, nega-tive thermal expansion coefficient below 4 o C etc.) as wellas dynamic anomalies in its supercooled state [1–10]. Onepopular, yet controversial, interpretation of the observedanomalies invokes the concept of liquid-liquid transitionand Widom line in supercooled water [11–13]. Poole etal. [11] observed the existence of two liquid phases – highdensity liquid (HDL) and low density liquid (LDL) in su-percooled ST2 model of water. These two liquid phasesundergo first order liquid-liquid phase transition (LLPT)on supercooling. The HDL-LDL coexistence line ends atliquid-liquid critical point (LLCP). We must note that insimulations one observes LLPT at thermodynamic con-ditions far from ambient conditions [11]. Several recentexperiments also suggest the existence of liquid-liquidtransition in water [14–16] as well as in other molecu-lar liquids [17–19]. However, the unambiguous confirma-tion of existence of LLPT in bulk supercooled water ishampered by fast ice nucleation at thermodynamic con-ditions predicted for the existence of LLPT. Alternativeinterpretations of the increase in response functions uponsupercooling also exist that does not invoke the conceptof LLPT [20, 21].Despite water being the most ubiquitous as well as themost studied liquid on earth, the origin of anomalies insupercooled water is still debated [22–27]. Most surpris-ingly, while much discussion has focused on the pos-sible presence (or absence) of HDL/LDL phases andLLCP [22–27], relatively fewer studies have focused on ∗ Corresponding author: [email protected] crystallization of ice and the effects of pronounced ther-modynamic and dynamic anomalies observed in super-cooled water on the pathway of ice nucleation [22, 23,28, 29].Nucleation of a new phase in correlated molecular sys-tems (such as water) can be a highly complex process.There are a multitude of factors (both thermodynamicand dynamic) for this difficulty. The interplay betweendifferent length and energy scales present in molecularsystems gives rise to unusually slow relaxation of the sys-tem and also rich phase diagram (polymorphism) [30] insupercooled state. Rich phase behavior can offer diverse(non-classical) nucleation pathways [31–33]. Both thesefactors along with observed thermodynamic and dynamicanomalies in supercooled water enhance the complexityof ice nucleation from supercooled water and poses fun-damental limitations to the classical nucleation theory(CNT) [34–36].
II. DENSITY FUNCTIONAL THEORETICALFORMULATION OF ICE NUCLEATION
Two major assumptions of CNT [34–36] are (i) thecapillary approximation that alllows us to write the freeenergy of a growing nucleus in terms of a sum of the sur-face energy and the bulk free energy terms, and (ii) thespherical shape of the nucleus. It is virtually impossibleto remove these assumptions in any self-consistent waywithout making use of substantial numerical work whichdestroys the simplicity of CNT. However, these approxi-mations might not introduce errors of many orders unlesswe approach the spinodal line.Density functional theory (DFT) [37–41] allows us toaddress the problem of the free energy of growing nu-cleus without making the capillary approximation. In
P T W i do m li n ec o e x i s t e n ce HDLLDL L D L s p i nod a l H D L s p i n o d a l (a) LLCPIIIIII ρ ω icek l = 1 x 10 k l = 5 x 10 k l = 1 x 10 ∆ l (b) ρ ω HDLLDL LLCP
ICE ∆ l (c) ρ ω ICE LDL HDL ∆ LDL ∆ HDL (d)
FIG. 1. (a) A schematic phase diagram of supercooled water. LLCP stands for liquid-liquid critical point. (b) The free energysurfaces along the green dotted lines indicated by I, II and III are shown in (b), (c) and (d), respectively. Note that for simplicitywe have only shown the dependence of the curvature of supercooled liquid phase free energy basin on changing thermodynamicconditions.
DFT one gets directly the (unstable) equilibrium densityprofile by minimizing the grand potential of the inhomo-geneous system,Ω[ ρ ( r )] = Z d r (cid:20) ω ( ρ ( r )) + 12 k ( ∇ ρ ( r )) (cid:21) , (1)with respect to density profile ( ρ ( r )) and correspondingfree energy cost for formation of the density profile ofcritical cluster. In the above equation, k is density cor-relation length and ω is the grand potential density. Thenon-local effects in the system due to inhomogeneity inthe density are accounted for in the square gradient term. A. Phase diagram of supercooled water
As mentioned, ice nucleation in water could be sub-stantially different due to the possible presence of multi-ple intermediate thermodynamic phases and also fromthe fact that the participating free energy surfaces(FES’s) of the parent and daughter phases can be quite different from each other. Thus, the FES of supercooledwater can become soft (as we approach ≈ K ) due tomultiple reasons. In Fig. 1(a), we have shown a schematicphase diagram of supercooled water. Liquid-liquid coex-istence line ends at LLCP. Beyond LLCP, the red dot-ted line shows the Widom line (locus of the maxima ofthermodynamic response functions). As evident from thefigure, depending on the change of controlled thermo-dynamic variables nucleation scenario will differ. TheFES’s (emphasizing the change in the curvature of liquidfree energy basin) on changing thermodynamic variablesalong the dotted green arrows indicated by I, II and III inFig. 1(a) are shown in Fig. 1(b), Fig. 1(c) and Fig. 1(d),respectively.In this article, we have used one order parameter DFTapproach to explore the effects of thermodynamic anoma-lies and distinct metastable liquid phases in supercooledwater on ice nucleation by studying ice nucleation in dif-ferent regions of phase diagram and made a number ofstriking predictions that can be tested either by computersimulation or controlled experiments. The present analy- r ρ k l = 1 x 10 k l = 3 x 10 k l = 5 x 10 k l = 7 x 10 k l = 10 x 10 (a) k l ∆ Ω ∗ (b) CNT
10 15
30 1000 2000 k l ∆ l ∆ Ω * (c) FIG. 2. (a) Density profiles of the critical cluster for differentvalues of the curvature of liquid free energy basin ( k l ). Wehave chosen k i = 2 × , density of ice, ρ i = 0 .
90, densityof liquid, ρ l = 1 . l = 25 . k l . Thedotted line (with diamonds) indicates the nucleation barrierpredicted by CNT. (c) The dependence of free energy barrieron both curvature of liquid free energy basin ( k l ) and freeenergy gap between supercooled liquid and ice basins (∆ l ).Note the linear decrease of the free energy barrier as well asflattening of the interface on decreasing k l . sis offers an explanation of the proximity of homogeneousnucleation temperature to the apparent LLCP. This DFTbased approach can also be extended to multiple orderparameters [39, 42]. B. Effects of softening of free energy surface on icenucleation
The anomalous increase in thermodynamic responsefunctions on approaching towards LLCP (dotted greenarrow I) indicates the softening of the metastable liq-uid FES. As illustrated in Fig. 1(a), on moving alongthe Widom line towards the LLCP, two thermodynamicvariables (pressure and temperature) change. Increasein pressure has an opposite effect (favors high densityphase) to that of decrease in temperature (favors ice-likelow density phase) on density as well as free energy gapbetween supercooled water and ice. Considering theseopposite dependencies and sake of simplicity we have ne-glected the thermodynamic condition dependence of boththe density of supercooled water as well as the relativefree energy gap between two phases.We include the softening of in FES through Landautype free energy expansions in order parameter(s). Thegrand potential densities of ice and metastable liquidphases are assumed as, ω i ( ρ ) = 12 k i ( ρ − ρ i ) ω l ( ρ ) = 12 k l ( P, T ) ( ρ − ρ l ( P, T )) + ∆ l ( P, T ) , (2)where k i is the curvature of ice free energy basin, ρ i isthe equilibrium density of ice phase and ρ l is the equi-librium density of metastable liquid phase. ∆ l is the su-percooling parameter. The values of these parameters arementioned in the caption of Fig. 2. At a particular super-saturation we can evaluate the density profile of the crit-ical nucleus by solving analytically the Euler-Lagrangeequation, δ Ω [ ρ ( r )] /δρ ( r ) = 0, with appropriate bound-ary conditions (see Appendix A for details). Ω [ ρ ( r )] isdefined by Eq. 1 with correlation length, k = 1 × and ω ( ρ ) = min [ ω i ( ρ ) , ω l ( ρ )]. The advantage of such a sim-plified, yet representative, approach is that one can solveEuler-Lagrange equation analytically (still-non-trivial).The expression for density profile of critical cluster is, ρ = ρ i + ( ρ c − ρ i ) r c r × " exp (cid:0) r √ c i (cid:1) − exp (cid:0) − r √ c i (cid:1) exp (cid:0) r c √ c i (cid:1) − exp (cid:0) − r c √ c i (cid:1) , ρ < ρ c ρ = ρ l + ( ρ c − ρ l ) r c r exp ( − ( r − r c ) √ c l ) , ρ > ρ c (3)where c i = k i /k and c l = k l /k . ρ c is the density of thesystem at r = r c . Relationship between ρ c and r c canbe established by equating the derivatives of density pro-files for ρ < ρ c and ρ > ρ c at ρ = ρ c . This conditionis necessary for the smoothness of the composite densityprofile at r = r c . The supercooling parameter (∆ l ) is re-lated with ρ c as, 2∆ l = k i ( ρ c − ρ i ) − k l ( ρ c − ρ l ) . Thedensity profiles given by Eq. 3 for different curvatures ofliquid FES are shown in Fig. 2(a). We note the broad-ening of the interface as well as decrease of the criticalcluster size on softening of supercooled liquid FES. Thenucleation barrier is the extra energy cost (with respectto metastable homogeneous liquid phase) for the forma-tion of (unstable) equilibrium density profile of criticalcluster and is given as, ∆Ω ∗ = Ω ( ρ ( r )) − Ω ( ρ l ).The final expression for the nucleation barrier (∆Ω ∗ ) caneasily be derived using the analytical expression of den-sity profile from Eq. 3 and is∆Ω ∗ = 2 πkr c (1 + r c √ c l ) ( ρ c − ρ l ) ( ρ i − ρ l ) − π r c ∆ l . (4)The effect of softening of the metastable liquid phaseFES is also reflected in the plot of ∆Ω ∗ vs . k l , shownin Fig. 2(b). Fig. 2(c) consists a two dimensional plotshowing dependence of ∆Ω ∗ on both ∆ l and k l . Notethe linear increase of ∆Ω ∗ on increasing k l as well asstronger dependence of ∆Ω ∗ on k l at low supersatura-tion compared to high supersaturation (see Fig. 2(c)).The dotted line in Fig. 2(b) indicates the nucleation bar-rier predicted by CNT, ∆Ω ∗ CNT = 16 πγ / l , where γ is the surface tension at coexistence. The surface ten-sion at coexistence (∆ l = 0) is the extra energy costper unit area for the formation of planner interface andis given by, γ = (cid:0) Ω ( ρ ( z )) − Ω (cid:0) ρ l/ice (cid:1)(cid:1) /A , where A isthe surface area of the planner interface. The coexis-tence density profile can be easily calculated by solvingthe Euler-Lagrange equation, δ Ω [ ρ ( z )] /δρ ( z ) = 0, withappropriate boundary conditions (see Appendix B for de-tails). The expression for surface tension is, γ = √ k i k ρ c − ρ i )( ρ l − ρ i ) . (5)As, inclusion of the effects of curvature increases thesurface tension, DFT prediction of nucleation barrier islarger than the nucleation barrier predicted by CNT withcoexistence (planner interface) surface tension.It should be noted that on increasing supercoolingthe relaxation time of the system also increases and onewould expect that the non-equilibrium effects (especiallyin correlated molecular systems) will also become morepronounced. These non-equilibrium effects are not ac-counted in our study. C. Beyond the harmonic approximation: role ofanharmonicity of free energy surface on nucleationphenomena
The majority of theoretical studies of phase transi-tions assume the participating free energy surfaces as ρ ω icea p = 0a p = 1 x 10 a p = 2 x 10 a p = 3 x 10 (a) ∆ l r ρ a p = 0a p = 1 x 10 a p = 2 x 10 a p = 3 x 10 p (x 10 )2.42.62.83.03.2 ∆ Ω ∗ ( x ) (b) FIG. 3. (a) Free energy surfaces as a function of densityfor ice and metastable liquid phases for different values ofanharmonicity parameter ( a p ). The chosen parameters are - k i = 2 × , k l = 1 × , ρ i = 0 . ρ l = 1 . l =25 . a p . In inset, the dependence of nucleation barrier (∆Ω ∗ )on a p is shown. The points are fitted with a straight line,∆Ω ∗ = ma p + c , where m = − . c = 3091 .
78. Notethe linear decrease in ∆Ω ∗ on increasing anharmonicity in theliquid phase free energy surface. harmonic. However, in computer simulation studies oneoften observes an anharmonic softening of the metastableliquid phase FES on supercooling. As on increasing an-harmonicity, the bulk free energy barrier (which in turnrelated to surface tension) decreases, one would expectthat the enhance anharmonic fluctuations on increasingsupercooling has pronounced effect on nucleation. In thissection, we quantify the effects of anharmonic softeningof the metastable liquid FES on nucleation barrier anddensity profile of the critical cluster. The grand potentialdensities of ice and metastable liquid phases are assumedas, ω i ( ρ ) = 12 k i ( ρ − ρ i ) ω l ( ρ ) = 12 k l ( ρ − ρ l ) + a p ( ρ − ρ l ) + ∆ l , (6)where a p is the anharmonicity parameter and is a mea-sure of the anharmonicity of the liquid phase FES (grandpotential density).For a fixed supersaturation (∆ l ), the grand potentialdensities for different values of anharmonic parameter areshown in Fig. 3(a). The density profile is computed us-ing relaxation method [43]. In Fig. 3(b), we have shownthe dependence of the density profile of critical cluster aswell as nucleation barrier (see the inset of Fig. 3(b)) onanharmonicity parameter. Note that on increasing an-harmonicity of the metastable liquid phase FES, the sizeof the critical cluster as well as the nucleation free en-ergy barrier decreases. This decrease can be attributedto the decrease in the bulk free energy barrier (shown inFig. 3(a)) on increasing a p . Thus, anharmonic softeningof the FES on increasing supercooling has an importantrole in decreasing the free energy barrier of nucleationand hence enhancing the rate of nucleation. D. Nucleation near liquid-liquid critical point
Now, we discuss the nucleation of ice from supercooledwater near LLCP (indicated by dotted green arrow IIin Fig. 1(a)). In order to study the effects of metastableLLCP on the ice nucleation scenario we have assumed thegrand potential densities of ice and supercooled water as, ω i ( ρ ) = 12 k i ( ρ − ρ i ) ω l ( ρ ) = 12 k l ( T ) ( ρ − ρ l ( T )) + ∆ l ( T ) . (7)The temperature dependence of k l is assumed as, k l ( T ) = (cid:16) a l ( T − T c ) + 1 (cid:17) k c , with T c = 5 . k c = 1000. T c is thecritical temperature and k c is the curvature of metastableliquid free energy basin at T c . The temperature depen-dence of density of metastable liquid, ( ρ l ) is assumed as, ρ l ( T ) = ρ c + ( T − T c )∆ ρ l , with critical density, ρ c = 1 . ρ l = 0 .
02. The grand potential densities corre-sponding to different phases are shown in Fig. 1(c).The calculated nucleation barrier and density profiles areshown in Fig. 4. The temperature is scaled with T c , T r = T /T c . The blue line indicates the dependence ofnucleation barrier on temperature when ∆ l is assumed tobe independent of temperature (∆ l = 25 .
0) and a l = 5for LDL basin and a l = 1 for HDL basin and the blackline indicates the same when the curvatures of metastableHDL and LDL basins are same ( a l = 1). The red lineindicates the case same as blue line, however, the tem-perature dependence of ∆ l ( T ) is also taken into account.In this case, ∆ l is varied in such as way that ∆ l = 25 . T r = 1 . l = 40 . T r = 0 . there is a crossover in thetemperature dependence of free energy barrier for ice nu-cleation near LLCP . Similar crossover behavior is alsoobserved in the temperature dependence of critical clus-ter size (see inset of the figure). We also note that the T r ∆ Ω ∗ r r * (a) T c r ρ T r = 0.7T r = 0.8T r = 0.9T r = 1.0T r = 1.1 (b) FIG. 4. (a) Dependence of nucleation barrier on tempera-ture near liquid-liquid critical point (LLCP). The black lineindicates the case when curvatures of free energy basins cor-responding to LDL and HDL phases are same ( a l = 1). Thered line indicates the case when the curvature of LDL basinis larger ( a l = 5) than the HDL basin ( a l = 1) and the blueline indicates the case when the free energy gap between themetastable liquid and ice is fixed corresponding to its valueat T = T c and LDL basin curvature is larger than the HDLbasin. In inset, we have shown the dependence of criticalcluster size on supercooling for the conditions indicated bycorresponding colors in ∆Ω ∗ vs. T r plot. (b) Density profilesof the critical cluster near LLCP at different supersaturationsfor the case when the free energy gap between the metastableliquid and ice is fixed corresponding to its value at T = T c and LDL basin curvature is larger than the HDL basin. Notethat temperature is scaled with T c . free energy barrier is relatively insensitive to the super-cooling near LLCP.In Fig. 4(b), we have shown the temperature dependenceof the density profile of critical ice nucleus for the casewhen the free energy gap between the metastable liquidand ice is fixed corresponding to its value at T = T c , andLDL basin curvature is larger than the HDL basin. Notethe increase in the stiffness of density profiles as well asdecreases in the density difference between bulk phases(ice and supercooled water) on supercooling. Decrease indensity difference between bulk phases reduces the sur-face tension cost for formation of ice nucleus; however, atthe same time increase in the stiffness of density profile(or, curvature of LDL basin) leads to an increase in thesurface tension. This delicate balance leads to an inter-esting crossover behavior in the temperature dependenceof nucleation barrier near T c . This observed crossoverbehavior of ice nucleation barrier can be tested in com-puter simulation studies. We must note that, dependingon the stiffness of LDL basin with respect to HDL basin,one might observe a significant increase in the free en-ergy barrier of nucleation (in place of being insensitiveor weakly sensitive) on increasing supercooling just below T c . The exact nature of crossover can only be quantifiedby inserting more realistic parameters in our theoreticalformalism which is not available at the moment.The nucleation scenario near LLCP has striking sim-ilarity to the (non-classical) pathway of protein crystal-lization near metastable gas-liquid critical point. Themetastable critical point enhances density fluctuationand thus decreases the nucleation barrier [44]. One mightalso observe a similar non-classical nucleation pathwaynear LLCP by invoking a two order parameter descrip-tion (density and order) where ice nuclei will grow insidelow density ice-like domains formed due to large scaledensity fluctuations in the system. E. Wetting mediated nucleation pathway
In the last section, we discuss the ice nucleation sce-nario at thermodynamic conditions where one observesdistinct metastable minima for LDL and HDL phases(indicated by dotted arrow III in Fig. 1(a)). The grandpotential densities of ice, LDL and HDL phases are as-sumed as, ω i ( ρ ) = 12 k i ( ρ − ρ i ) ω LDL ( ρ ) = 12 k LDL ( ρ − ρ LDL ) + ∆ LDL ω HDL ( ρ ) = 12 k HDL ( ρ − ρ HDL ) + ∆ HDL , (8)where ∆ LDL and ∆
HDL are the (meta)stabilities of LDLand HDL phases with respect to ice phase. ρ i , ρ LDL and ρ HDL are the densities of bulk ice, LDL and HDLphases, respectively. For numerical computation we haveassumed the bulk densities and curvatures of grand po-tential densities of ice, LDL and HDL phases as indepen-dent of temperature. The grand potential densities ofdifferent phases at thermodynamic condition when HDLis metastable with respect to both LDL and ice phases areshown in Fig. 1(d). In order to get the density profile ofcritical cluster, we have solved the corresponding Euler-Lagrange equation using relaxation method [43]. Thedependence of density profile of critical ice nucleus on(meta)stability of HDL phase with respect to ice (∆
HDL )is shown in Fig. 5(a). Note that we have neglected the r ρ ∆ HDL = 50.0 ∆ HDL = 60.0 ∆ HDL = 70.0 ∆ HDL = 80.0 ∆ HDL = 90.0 (a)
LDLICE (b)HDL
FIG. 5. (a) The dependence of density profile of critical icenucleus on the stability of ice phase with respect to HDLphase ∆
HDL . Note the crossover from wetting mediated tran-sition to one step Ostwald step rule type of scenario. For thecomputation of density profile, we have chosen ρ i = 0 . ρ LDL = 0 . ρ HDL = 1 .
05 and ∆
LDL = 25 .
0. The cur-vatures of grand potential densities for different phases are k i = 2 × , k LDL = 5 × and k HDL = 1 × . (b) Aschematic diagram showing the wetting of ice nucleus by anintermediate metastable LDL phase. supercooling dependence of ∆ LDL with respect to ∆
HDL (∆ LDL is fixed at 25 . ρ i is relatively closer to ρ LDL than ρ HDL . We have also neglected the non-equilibriumeffects arising due to rapid increase of relaxation time ofLDL phase on decreasing temperature.As evident from Fig. 5(a), when the HDL phase hasminimal metastability with respect to the ice phase, aone step density profile (without pronounced wiggle) in-dicates the absence of (or negligible) wetting of ice nu-cleus by intermediate metastable LDL phase. On gradu-ally increasing the stabilities of ice and LDL phases withrespect to HDL phase, we observe a significant deviationin the density profile of the critical cluster. This indicatesa change in the composition of the critical cluster of ice byan intermediate LDL phase. In Fig. 5(b), we have showna schematic diagram of the ice nucleus wetted by an in-termediate metastable LDL phase. Note the extent ofwetting (width of the metastable LDL region) depends onthe stability of LDL phase. On further increasing ∆
HDL we observe a transition where a critical cluster of inter-mediate LDL phase appears inside the bulk metastableHDL phase. This is the Ostwald step rule scenario, wheretransition from metastable HDL phase to final stable icephase occurs via sequential transformations.Thus, on increasing supercooling we observe acrossover from the wetting enhanced one step transitionto a sequential two step (following Ostwald step rule)transition. The existence of a metastable LDL/LDA-like phase, with order/density intermediate between HDLand ice can greatly facilitate the nucleation of ice fromHDL metastable phase by decreasing the surface tensionbetween stable ice and metastable HDL phases. Recentcomputer simulation studies of freezing of water by Mat-sumoto et al. [28] as well as Moore et al. [29] indicatethe wetting of ice nucleus by a metastable low densityphase. This type of wetting mediated transition has alsobeen predicted in the nucleation of liquid from a glassyphase [45], crystallization of simple (such as hard sphere)and complex (protein crystallization) systems [42, 44–49].The application of CNT in the case of wetting mediatedtransitions is flawed as it does not take into account thedecrease in surface tension due to wetting (indirect par-ticipation of intermediate phase(s)) of the nucleus.
III. CONCLUSION
To summarize, we have used DFT with phenomeno-logical free energy surfaces to explore the thermodynamic condition dependent diverse plausible pathways of ice nu-cleation from supercooled water. We show that both thesoftening of the free energy surface (due to even a distantpresence of LLCP in the phase plane) and the presenceof distinct metastable liquid phases (LDL and HDL) canlower significantly the nucleation free energy barrier. Asboth the two situations lower the free energy barrier, weneed to look for the features that can allow us to uniquelyidentify the actual pathway, at least in computer simu-lations. The main discernible difference occurs in theorder parameter profiles as we move from the core of thegrowing ice nucleus to the surface. However, a detailedanalysis of the order parameter profile around a criticalis yet to be carried out. We also observe an interestingcrossover in the temperature dependence of nucleationbarrier near LLCP.It should be noted that, in computer simulation stud-ies, the phase diagram as well as formation of ice fromsupercooled water is quite sensitive to the force field used.This often leads to uncertainty and even conflicting re-sults.
ACKNOWLEDGMENTS
We thank Dr. Mantu Santra, Prof. Shinji Saito andProf. Iwao Ohmine for help and many stimulating dis-cussions. We thank the Department of Science and Tech-nology (DST) and the Board of Research in Nuclear Sci-ences (BRNS), India, for partial financial support for thiswork. B.B. thanks DST for J. C. Bose fellowship. [1] F. Franks,
Water: A Matrix for Life (Royal Society ofChemistry, Cambridge, 2000)[2] A. G. Smart, Phys. Today , 16 (2013)[3] R. J. Speedy and C. A. Angel, J. Chem. Phys. , 851(1976)[4] C. A. Angel, Annu. Rev. Phys. Chem. , 593 (1983)[5] O. Mishima, J. Chem. Phys. , 144503 (2010)[6] C. A. Angell, J. Shuppert, and J. C. Tucker, J. Phys.Chem. , 3092 (1973)[7] D. E. Hare and C. M. Sorensen, J. Chem. Phys. , 4840(1987)[8] V. Holten, C. E. Bertrand, and J. V. Anisimov,M. A. Sengers, J. Chem. Phys. , 094507 (2012)[9] S. Saito, I. Ohmine, and B. Bagchi, J. Chem. Phys. ,094503 (2013)[10] B. Jana, R. S. Singh, and B. Bagchi, Phys. Chem. Chem.Phys. , 16220 (2011)[11] P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley,Nature , 324 (1992)[12] O. Mishima, J. Chem. Phys. , 5910 (1994)[13] T. Loerting and N. Giovambattista, J. Phys.: Condens.Matter , R919 (2006)[14] O. Mishima and H. E. Stanley, Nature , 164 (1998)[15] O. Mishima, Phys. Rev. Lett. , 334 (2000)[16] K. Amann-Winkel, C. Gainaru, P. H. Handle, M. Seidl, H. Nelson, R. Bhmer, and T. Loerting, Proc. Nat. Acad.Sci. , 17720 (2013)[17] M. Durandurdu and D. A. Drabold, Phys. Rev. B ,041201 (2002)[18] R. Kurita and H. Tanaka, Science , 845 (2004)[19] M. Beye, F. Sorgenfrei, W. F. Schlotter, W. Wurth, andA. Fhlisch, Proc. Nat. Acad. Sci. , 16772 (2010)[20] R. J. Speedy, J. Phys. Chem. , 982 (1982)[21] S. Sastry, P. G. Debenedetti, F. Sciortino, and H. E. Stan-ley, Phys. Rev. E , 6144 (1996)[22] D. Limmer and D. Chandler, J. Chem. Phys. , 134503(2011)[23] D. Limmer and D. Chandler, J. Chem. Phys. , 214504(2013)[24] F. Sciortino, I. Saika-Voivod, and P. H. Poole, Phys.Chem. Chem. Phys. , 19759 (2011)[25] P. H. Poole, R. K. Bowles, I. Saika-Voivod, andF. Sciortino, J. Chem. Phys. , 034505 (2013)[26] Y. Liu, J. C. Palmer, A. Z. Panagiotopoulos, and P. G.Debenedetti, J. Chem. Phys. , 214505 (2012)[27] T. A. Kesselring, E. Lascaris, G. Franzese, S. V.Buldyrev, H. J. Herrmann, and H. E. Stanley, J. Chem.Phys. , 244506 (2013)[28] M. Matsumoto, S. Saito, and I. Ohmine, Nature , 409(2002) [29] E. B. Moore and V. Molinero, Nature , 506 (2011)[30] R. S. Singh, M. Santra, and B. Bagchi, J. Chem. Phys. , 184507 (2013)[31] S. Whitelam, Phys. Rev. Lett. , 088102 (2010)[32] S. Whitelam, J. Chem. Phys. , 194901 (2010)[33] L. O. Hedges and S. Whitelam, J. Chem. Phys. ,164902 (2011)[34] R. Becker and W. Doring, Ann. Phys. , 719 (1935)[35] A. C. Zettlemoyer, Nucleation (Dekker: New York, 1969)[36] P. G. Debenedetti,
Metastable Liquids: Concepts andPrinciples (Princeton University Press, 1996)[37] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. , 258(1958)[38] Y. C. Shen and D. W. Oxtoby, Phys. Rev. Lett. , 3585(1996)[39] V. Talanquer and D. W. Oxtoby, J. Chem. Phys. ,223 (1998)[40] L. Granasy and D. W. Oxtoby, J. Chem. Phys. , 2410(2000)[41] C. K. Bagdassarian and D. W. Oxtoby, J. Chem. Phys. , 2139 (1994)[42] M. Santra, R. S. Singh, and B. Bagchi, J. Phys. Chem.B , 13154 (2013)[43] W. H. Press, S. A. Teukolsky, W. T. Vitterling, andB. P. Flannery, Numerical Recipes (Cambridge Univer-sity Press, 1986)[44] P. R. ten Wolde and D. Frenkel, Science , 1975 (1997)[45] X. Xia and P. G. Wolynes, Proc. Nat. Acad. Sci. , 2990(2000)[46] T. Kawasaki and H. Tanaka, Proc. Nat. Acad. Sci. ,14036 (2010)[47] J. Russo and H. Tanaka, Scientific Reports , 505 (2012)[48] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel,Phys. Rev. Lett. , 2714 (1995)[49] P. R. ten Wolde and D. Frenkel, Phys. Chem. Chem.Phys. , 2191 (1999) Appendix A: Analytical expression for nucleationfree energy barrier
At a particular supersaturation, density profile of thecritical nucleus can be evaluated by solving the Euler-Lagrange equation, δ Ω [ ρ ( r )] /δρ ( r ) = 0. The resultingequation is k ∇ ρ ( r ) = dω ( ρ ( r )) dρ ( r ) . (A1)The above equation in spherical coordinate can be rewrit-ten as, d ρdr + 2 r dρdr = dωdρ , (A2)where ω ( ρ ) = min [ ω i ( ρ ) , ω l ( ρ )]. This equation asclose resemblance with free particle time independentSchrdinger equation in spherical coordinate. One caneasily derive an analytical expansion for density profileby substituting ρ − ρ i/l = y and y = u/r . The boundaryconditions are, dρ/dr = 0 at r = 0 and ρ = ρ l at r → ∞ . The expression for density profile is, ρ = ρ i + ( ρ c − ρ i ) r c r × " exp (cid:0) r √ c i (cid:1) − exp (cid:0) − r √ c i (cid:1) exp (cid:0) r c √ c i (cid:1) − exp (cid:0) − r c √ c i (cid:1) , ρ < ρ c ρ = ρ l + ( ρ c − ρ l ) r c r exp ( − ( r − r c ) √ c l ) , ρ > ρ c (A3)The relationship between ρ c and r c can be established byequating, dρdr (cid:12)(cid:12)(cid:12)(cid:12) ρ<ρ c r = r c = dρdr (cid:12)(cid:12)(cid:12)(cid:12) ρ>ρ c r = r c . (A4)For each ρ c , one can get r c by solving numerically thenon-linear algebraic equation given by Eq. A4. The freeenergy barrier of nucleation is the grand potential for for-mation of critical cluster in bulk metastable liquid phase,∆Ω ∗ = Ω ( ρ ( r )) − Ω ( ρ l )= 4 π Z ∞ drr (cid:20) ω ( ρ ( r )) + 12 k ( ∇ ρ ( r )) − ω ( ρ l ) (cid:21) . (A5)The above equation can be splitted as, ∆Ω ∗ = ∆Ω ∗ i +∆Ω ∗ l , where∆Ω ∗ i = 4 π Z r c drr (cid:20) ω i ( ρ ( r )) + 12 k ( ∇ ρ ( r )) − ω ( ρ l ) (cid:21) , (A6)and∆Ω ∗ l = 4 π Z ∞ r c drr (cid:20) ω l ( ρ ( r )) + 12 k ( ∇ ρ ( r )) − ω ( ρ l ) (cid:21) . (A7)Now, on integrating the above equation and on furthersimplification using the relationship between ρ c and r c (Eq. A4), the expression for nucleation barrier reducesto∆Ω ∗ = 2 πkr c (1 + r c √ c l ) ( ρ c − ρ l ) ( ρ i − ρ l ) − π r c ∆ l . (A8) Appendix B: Analytical expression for surfacetension at coexistence
Surface tension is the extra energy cost for theformation of planner interface between two coexist-ing phases at equilibrium and is defined as, γ = (cid:0) Ω ( ρ ( z )) − Ω (cid:0) ρ l/ice (cid:1)(cid:1) /A , where A is area of the plan-ner interface and Ω ( ρ ( z )) is grand potential for inhomo-geneous density profile ρ ( z ) and Ω (cid:0) ρ l/ice (cid:1) is the grandpotential of the bulk phases at coexistence. Followingsimilar procedure one can derive an expression for sur-face tension, γ = √ k i k ρ c − ρ i )( ρ l − ρ i ) . (B1)Now, when the curvature of the grand potential densi-ties corresponding to ice and liquid phases are same, i.e. k i = k l = k ′ , then ρ c = ( ρ i + ρ l ) / γ = √ k ′ k ρ ) , (B2) where ∆ ρ = ( ρ l − ρ i ) is the order parameter differencebetween two coexisting phases and kk