Nucleation barriers in tetrahedral liquids spanning glassy and crystallizing regimes
aa r X i v : . [ c ond - m a t . s o f t ] J un Nucleation barriers in tetrahedral liquids spanning glassy and crystallizing regimes
Ivan Saika-Voivod, Flavio Romano, and Francesco Sciortino Department of Physics and Physical Oceanography,Memorial University of Newfoundland, St. John’s, NL, A1B 3X7, Canada Dipartimento di Fisica, Sapienza – Universit`a di Roma, P.le A. Moro 5, 00185 Roma, Italy Dipartimento di Fisica and CNR–ISC, Sapienza – Universit`a di Roma, P.le A. Moro 5, 00185 Roma, Italy (Dated: September 7, 2018)Crystallization and vitrification of tetrahedral liquids are important both from a fundamentaland a technological point of view. Here, we study via extensive umbrella sampling Monte Carlocomputer simulations the nucleation barriers for a simple model for tetrahedral patchy particles inthe regime where open tetrahedral crystal structures (namely cubic and hexagonal diamond andtheir stacking hybrids) are thermodynamically stable. We show that by changing the angular bondwidth, it is possible to move from a glass-forming model to a readily crystallizing model. From theshape of the barrier we infer the role of surface tension in the formation of the crystalline clusters.Studying the trends of the nucleation barriers with the temperature and the patch width, we areable to identify an optimal value of the patch size that leads to easy nucleation. Finally, we findthat the nucleation barrier is the same, within our numerical precision, for both diamond crystalsand for their stacking forms.
I. INTRODUCTION
Crystallization is central to several fields, from materi-als research to biological science, not to mention its tech-nological relevance [1, 2]. Several human pathologies arealso caused by crystal nucleation in protein solutions [3].Understanding crystal nucleation requires both the eval-uation of the stability fields of the fluid and crystal phases(i.e. knowledge of the chemical potentials of all possiblephases), as well as the evaluation of the thermodynamicbarriers controlling the formation of the stable phase andof the kinetic pre-factors fixing the timescale of the dif-fusional processes.In recent years, several numerical methodologies havebeen developed for accurately evaluating phase diagramsfrom the free energies of the fluid and crystal phases. Werefer the reader to the review by Vega and coworkers [4].Also for crystallization, various methods are now avail-able for calculating free energy barriers and nucleationrates [5–9], making it possible to generate accurate datafor model potentials and, more importantly, to comparethe numerical results with theoretical predictions, mostlybased on classical nucleation theory (CNT) [10–15], aswell as, when possible, with experimental data.One of the main motivations for the development ofthe special methods for studying nucleation is the propersampling of the equilibrium cluster size distribution N ( n )within the metastable liquid, i.e. the number of crystal-like clusters composed of n particles. The work of forminga cluster of size n from the liquid is given in terms of N ( n )by [8], β ∆ G ( n ) = − ln (cid:20) N ( n ) N p (cid:21) , (1)where N p is the number of particles in the system and β ≡ ( k B T ) − , where T is the temperature and k B isthe Boltzmann constant. Below the melting tempera-ture, β ∆ G ( n ) has a maximum value β ∆ G ∗ at a critical nucleus size n ∗ . β ∆ G ∗ enters into the expression for therate of nucleation exponentially, and is therefore the mainconsideration in determining the rate from a thermody-namic perspective. Very generally, the larger the differ-ence in chemical potential between liquid and crystal ∆ µ ,often referred to as the driving force for crystallization,the smaller β ∆ G ∗ is. Conversely, the surface tension γ between crystallite and surrounding liquid always actsagainst the growth of crystallites, and so β ∆ G ∗ increaseswith γ . Hence, the ability of a system to crystallize isgoverned by the interplay between ∆ µ and γ .Recently, simulation studies have begun to addresscrystallization of tetrahedral liquids [16–21]. This classof interesting liquids includes biological and technologi-cally important molecular and atomic materials such aswater, silica, silicon and carbon. One common featurein these systems is the formation of open, low densitystructures such as the diamond cubic (DC) crystal. Forcarbon, simulations have shown the importance of the liq-uid structure in governing nucleation barriers [22, 23]. Ina generalized form of the Stillinger-Weber model for sili-con [24], the degree of tetrahedrality was shown to havea strong impact on DC nucleation and its interrelationwith the appearance of a metastable liquid–liquid criti-cal point [25]. Further, for silicon and germanium, thelower density of the DC crystal with respect to the liquidwas shown to give rise to a preference for nucleation nearthe surface, with similar implications for water [26]. Theself-assembly of a DC colloidal crystals is also of interestin the field of photonics, since such an ordered structureof dielectric spheres is expected to exhibit a band struc-ture with a large gap within the visible wavelengths oflight [27].Tetrahedral liquids are interesting from another per-spective as well. Often, in these systems crystal forma-tion is inhibited by the onset of a glassy behavior whichdramatically slows down the microscopic dynamics, thuspreventing the transition to the ordered structure. Theglass forming ability of silica and the ease with whichwater crystallizes stand in stark contrast to each other.Motivated by these issues of scientific and technolog-ical interest, recently, some of us [28] have thoroughlyinvestigated the question of how to best design tetrahe-dral patchy colloidal particles that spontaneously assem-ble into open crystal structures through computer sim-ulations of the Kern-Frenkel (KF) model [29]. It wasshown that the angular width of the patch plays a keyrole in determining if the system, at low temperatures,forms a disordered glass structure or a crystal, a resultwhich may also be of interest for interpreting the glass-forming ability of atomic and molecular systems. Whenthe patch width is small, the crystal coexists with a fluidand β | ∆ µ | (the driving force for crystallization) increasesquickly with undercooling, while when the patch widthis large, the crystal coexists with a fully bonded liquidstate and β | ∆ µ | grows only moderately with undercool-ing. It was also shown that the model spontaneouslyforms a crystal composed of a mixture of two open tetra-hedral structures, the DC and the diamond hexagonal(DH) crystals, and that the chemical potentials for thetwo polymorphs were indistinguishable within the preci-sion of the calculations.In this article we proceed one step further by investi-gating, for the same tetrahedral model, the patch widthdependence of the free energy barrier to nucleating theDC and/or DH crystals. Interestingly, we find that, atcomparable undercooling, the barrier height β ∆ G ∗ doesnot monotonically decrease with increasing patch widthas one would expect from a consideration of β | ∆ µ | alone.Indeed, comparing the barrier shape with CNT predic-tions, we find that the surface tension γ increases on de-creasing patch width, a result that we tentatively connectto the larger structural and density difference betweenfluid and crystal. We do find that for narrow patches β ∆ G ∗ generally shows a strong T dependence, rapidlydecreasing to a homogenous nucleation limit, while forwide patches β ∆ G ∗ decreases more slowly with decreas-ing T , allowing the system to reach the glassy regime.Additionally, at the widest patch studied, the barriersremain very large (of the order of 50 k B T ) even for signifi-cant undercooling. This confirms that despite the benefitof a lower surface tension, the lack of a buildup of a largedifference in chemical potential is mainly responsible fordisfavoring crystallization of particles with wide patches.The presented results confirm that DC/DH will spon-taneously form when the angular width of the patch issufficiently small.The remainder of this article is organized as follows. InSection II, we outline our criteria for defining crystal-likeclusters and the umbrella sampling Monte Carlo proce-dure for calculating nucleation barriers. In Section IIIwe present our results, including a fairly detailed exami-nation of the robustness of β ∆ G ∗ to varying the cluster-defining criteria and the dependence of n ∗ and clustercomposition (e.g. percentage of DC particles) on them.We then present our summary and conclusions in Section BCC FluidDC / DH (d) cos θ = 0.920 0.1 0.2 T θ = 0.94 (b) cos θ = 0.9600.020.040.06 P (a) cos θ = 0.98 FIG. 1: Equilibrium phase diagrams for all the models stud-ied in the low P , low T region, where the open crystal struc-tures are stable. The short blue lines correspond to the tem-peratures for which the nucleation barriers have been investi-gated. The circles mark the gas–liquid critical point. For thenarrow patch models [(a) and (b)], nucleation is so effectivethat the location of the critical point can not be properly es-timated. However, we have checked that the studied isobar islocated significantly above the critical pressure. IV.
II. METHODS
We perform biased Monte Carlo (MC) simulations [33]at constant T and pressure P of the Kern-Frenkelmodel [29] with N p = 1000 particles, each with hardsphere diameter σ , and each having four tetrahedrallyarranged attractive patches of range δ = 0 . σ , strength u and angular width 2 θ [28, 34]. We report T and P as dimensionless quantities, after rescaling by u /k B and u /σ , respectively. We investigate four different valuesof θ , namely those corresponding to cos θ = 0 .
98, 0.96,0.94 and 0.92 ( θ = 11 . ◦ , 16.3 ◦ , 19.9 ◦ and 23.1 ◦ ). For thevalue of δ we use, the condition cos θ > . − u ). Thelow density crystal is an open tetrahedral structure andit can exist in two polymorphs, DC and DH. As notedin Ref. [28], the chemical potentials for the DH and DCstructures are largely indistinguishable, and indeed it hasbeen observed that when the liquid crystallizes sponta-neously to the open tetrahedral structure, it does so byforming a stacking of DC and DH planes [28]. Indeed, DHand DC crystalline layers can stack in an analogous wayto the FCC/HCP stacking in hard-spheres. The densephase, composed of two interpenetrating fully bondedDC structures, is a body centered cubic (BCC) crystal.The fluid separates into gas and liquid phases below thecritical temperature. Since the range of the potentialis short compared to the particle size, the critical pointis almost always metastable with respect to the crystalphase.In the following, we study the crystallization barriersto the DC/DH crystal at one selected pressure. Specifi-cally, we choose P = 0 .
03 in order to avoid interferencewith gas-liquid phase separation. For all the models stud-ied, with the exception of cos θ = 0 .
98, the most stablephase at that pressure is the open tetrahedral crystal.In the case of cos θ = 0 .
98, BCC is the stable phase atthis pressure, but as previous studies have shown, spon-taneous crystallization results always in the open struc-ture [28], an example of the Ostwald step rule [30–32].The nucleation to tetrahedral crystals has been stud-ied previously in simulations [22, 23, 26] and we followthe established methodology to evaluate the free energybarriers. A novel issue arises from the presence of twocrystals with different symmetries in the same stabilityfield, as we discuss in the following.We use Steinhardt bond order parameters [35] basedon spherical harmonics of order l = 3. For each particlewe define the complex vector q lm ( i ) = 1 N b ( i ) N b ( i ) X j =1 Y lm (ˆ r ij ) , (2)where the sum is over the N b ( i ) neighbors of particle i ,defined as those particles within a distance of (1 + δ ) σ =1 . σ of particle i . The dot product c ij = l X m = − l ˆ q lm ( i )ˆ q ∗ lm ( j ) , (3)where ˆ q lm ( i ) = q lm ( i ) / l X m = − l | q lm ( i ) | ! / (4)and ˆ q ∗ lm ( i ) is its complex conjugate, determines the de-gree of orientational correlation between neighboring par-ticles i and j . Fig. 2 shows the distributions for thefluid and the DC and DH crystals. The DC distributionis peaked around c ij = − c ij = − .
1. In the DH crystal,each particle has three neighbors with c ij ≈ − c ij ≈ − .
1, while in the DC crystal all four bondedneighbors have c ij ≈ −
1. This provides a local basis fordistinguishing particles as being DC or DH. The fluid dis-tributions are very wide and show sharp peaks for large values of cos θ (small bonding angles). The peaks becomemore intense on heating, which we attribute to a lowerdensity and higher energy; we identify the peaks as sig-nals of specific geometrical assemblies in the fluid withunfilled bonds. For example, a dimer has c ij = −
1, whileboth bonds in a trimer have c ij = − . c ij ranges in whichcrystal-like or fluid-like particles are mostly contributingoffers a way of associating a value of c ij with a localstructure (crystal or fluid). Usually a threshold numberof crystal-like connections is selected to distinguish par-ticles as being fluid-like or crystal-like. In the following,we start by defining a crystal-like connection betweenneighbors as one with c ij < q u , with q u = − .
87 anda crystal-like particle as one which has three or morecrystal-like connections. This definition does not allowone to discriminate between DC and DH and hence ap-pears to be a reasonable first step in the investigation ofthe nucleation barriers [22, 23, 26]. We complement thisstudy with an additional investigation where we differen-tiate between DC and DH by requiring that a solid-likeparticle have four neighbors with c ij < − .
87 in the DCcase or three neighbors with c ij < − .
87 and one with − . < c ij < . φ = κ ( n max − n ) , (5)where κ is a suitably chosen constant that controls therange of sampled crystal cluster sizes, centered near n ,and n max is the size of the largest crystalline cluster inthe system. Depending on the steepness of ∆ G ( n ), weadjust the value of κ to obtain good sampling, but foralmost all simulations, κ = 0 . n , we perform a trajectoryof 20 Metropolis MC steps, where one such step is on av-erage N p − n . The new configuration is acceptedwith probability min(1 , exp {− β [ φ ( n ) − φ ( n )] } ).If the new configuration is accepted, it becomes the start-ing point for the next MC trajectory. If it is rejected, theentire trajectory is discarded and the n configurationis kept and used as the starting point for generating an-other MC trajectory. For the slowest state points, weperform 10 US attempts. We perform several US sim-ulations for different values of n . Generally, two adja- T = 0.140 T = 0.150 T = 0.140 DC T = 0.140 DH T = 0.158 T = 0.165 T = 0.158 DH-1 -0.5 0 c ij P ( c ij ) T = 0.154 T = 0.168 T = 0.154 DH -1 -0.5 0 c ij T = 0.158 T = 0.165 T = 0.158 DH (a) cos θ = 0.98 (b) cos θ = 0.96(d) cos θ = 0.92(c) cos θ = 0.94 FIG. 2: Probability distributions of the dot product c ij in theliquid and crystal structures. The distribution for DC, shownonly in panel (a), is a single peak near -1 and is nearly indis-tinguishable from the DH peak at -1 on the scale of the plots.The distributions for the crystals do not vary significantlyover our T range. cent US windows are spaced by ∆ n = 3 to ensure goodsampling of the reaction path. For each US window wethen evaluate the solid cluster size distribution ˜ N ( n ) andP max ( n ), the probability that the largest cluster in thesystem is of size n .The cluster size distribution N ( n ) in the N P T ensem-ble for each window is worked back from the biased en-semble as in Ref. [8], with N ( n ) = D exp [ βφ ( n max )] ˜ N ( n ) E biased (6)up to a multiplicative constant, where h . . . i indicate anensemble average. Portions of β ∆ G ( n ) are determinedfrom simulations at different values of n up to additiveconstants, in agreement with Eq. 1. These pieces arematched by minimizing the difference between overlap-ping portions after discarding data for which P max ( n ) isless than 0.01 (to ensure good sampling). Alternatively,we find the free energy difference ∆ G ( n ) − ∆ G ( n −
1) as aweighted average of the values obtained with simulationsat different n in which clusters of size n and n − − k ( n − n )]. As a result ∆ G ( n ) = P nj =1 [∆ G ( j ) − ∆ G ( j − G (0) = 0. While this is not formallycorrect, since the number of liquid-like particles [ N (0)] isnot precisely equal to N p , the error is negligible (of theorder of 0.01 k B T or less). We stress here that at this n β ∆ G ( n ) n (a) cos θ = 0.98 (b) cos θ = 0.96(c) cos θ = 0.94 (d) cos θ = 0.92 TTT T
FIG. 3: Nucleation barriers β ∆ G ( n ) for various bond angularwidths and temperatures. Panel a) cos θ = 0 .
98, from bottomto top T = 0 . θ = 0 .
96, from bottom to top T = 0 . θ = 0 .
94, from bottom to top T = 0 . θ = 0 . T = 0 . stage we do not differentiate between DH and DC parti-cles. We have checked that indeed the largest crystallinecluster is a mixture of the two crystal local structures.We will address this point in more detail later on.To ensure independence of the initial conditions, weconstruct our barriers starting from two different setsof data. In the first set we use the biasing potentialto grow clusters in the fluid, starting the US simulationat n from a configuration extracted from the previouswindows. In the second set, we seed the biased simu-lations at cos θ = 0 .
98 with crystallite-containing con-figurations from homogeneously nucleating unbiased MCruns at cos θ = 0 .
98. Simulations at smaller values ofcos θ are started from the cos θ = 0 .
98 data set, afterequilibration at higher values of T . III. RESULTSA. Barrier profiles
In Fig. 3 we show our results for the ∆ G ( n ) profiles forthe four models of different patch widths at various T .The range of T where the barrier can be calculated withthe present methodology is restricted both from aboveand below for different reasons. At large T (small su-percooling) the critical nucleus is large compared to thesystem size. At low T , different reasons conspire againsta proper evaluation of the barrier height. For example,in the case cos θ = 0 .
92, the slow kinetics associated withthe proximity of a glass transition prevent proper equi-libration already at temperatures where the critical nu-cleus is of comparable size to the studied system. In thissame case, it has been suggested on the basis of the ex-plicit evaluation of the difference in chemical potentialbetween the crystal and the fluid that, at odds with thestandard behavior, barriers do not grow with further su-percooling. The difficulty of evaluating barriers for thewide patch model is consistent with this proposition. Inthe case of narrow patches, cos θ = 0 .
98, we are limitedto a barrier height of the order of 30 k B T . For lower T ,despite slow dynamics not being an issue, the results ofthe calculations become more and more affected by thepresence of secondary crystal clusters. While in theorythe presence of more than one cluster is not a problem, inour opinion this effect highlights an incorrect definitionof crystallinity at the local level, which artificially breaksa single cluster into two (or more) pieces. Hence, at thelowest T , the barrier calculations become difficult owingto the appearance of clusters of comparable size to thelargest one in the system. A visual inspection of con-figurations confirms the formation of secondary clustersoccurring next to the primary one, separated by particlesthat fall outside the cutoffs defining crystal-like particles.We recall that for a proper evaluation of the barriers inthe low T limit, where homogeneous nucleation is tak-ing place, one can apply the methodology based on meanfirst-passage times [36–40]. Finally, we note that we havechecked for state points where the barrier is of the orderof 20 k B T the quality of our calculations by comparingthe US results with N ( n ) evaluated in unconstrained sim-ulations. B. Fits
Within the phenomenological framework of CNT, thework of forming an n -sized cluster can be written as, β ∆ G ( n ) = − β | ∆ µ | n + βγA, (7)where A is the surface area of the cluster. We thereforefit our profiles to β ∆ G ( n ) = − a n + b n / . (8)where a = β | ∆ µ | and b ∼ βγ and we have assumed acompact cluster. Below, we compare our results for a against β | ∆ µ | , as some studies suggest a very good cor-respondence [8, 9, 41], and justify our assumption that b is equal to βγ up to a constant factor which depends onthe shape and density of the crystallites. Recently, it hasbeen shown for a soft-core colloidal model that the pro-cedure we follow here for determining n may not be suf-ficient for describing crystal-like structures, but that the n β ∆ G ( n ) cos θ = 0.98 T = 0.143cos θ = 0.92 T = 0.156CNT fit FIG. 4: Fit of two of the barrier profiles studied to the CNTform (Eq. 8). While there is a deviation at small n , the overallrepresentation of the curve is fairly good, especially in termsof the barrier height ∆ G ∗ and size of the critical nucleus n ∗ .Inset: low- n region of the same data, hilighting the deviationfrom the CNT form. scaling represented by Eq. 8 may still hold, i.e., changingthe parameters used to define a crystalline cluster willchange the numerical values of a and b [42]. We explorethis in some detail below.We show two representative β ∆ G ( n ) curves along fitsto Eq. 8 in Fig. 4. While there is a deviation between dataand fit at small n , the overall description is rather good.Moreover, the values for n ∗ and ∆ G ∗ extracted from thefits coincide with those obtained directly from the bar-riers. Improvement of the small n description could beaccomplished along the lines suggested in Ref. [43] but itis outside the scope of the present work.The fits allow us to plausibly extrapolate our ∆ G ( n )profiles to values of n larger than what we can simulatein our N p = 1000 system, allowing us to estimate n ∗ and β ∆ G ∗ for all the T that we study. While we do not relystrongly on the accuracy of the extrapolations, they doprovide a stronger sense of the trends in the data. C. T dependence of ∆ G ∗ and n ∗ Fig. 5(a) shows β ∆ G ∗ as a function of T /T m , where T m is the melting temperature for the four models [45][(cos θ , T m ): (0.98, 0.153), (0.96, 0.169), (0.94, 0.172),(0.92, 0.174)]. For small angles, the T dependence of β ∆ G ∗ is significant, and can be modeled rather well inthe present range of barriers with an exponential func-tion. The fast decrease of β ∆ G ∗ provides evidence forthe inevitability of crystallization. For larger angles, weobserve both a significant change in slope, i.e., a muchslower decrease of the barrier height with supercooling,as well as a progressive increase of the barrier height atfixed supercooling with increasing angle. The trend is ac-counted for by the results reported in [28], namely that β | ∆ µ | becomes a weaker function of temperature as θ in-creases. On the basis of the results presented in Fig. 5(a), T / T m β ∆ G * cos θ = 0.92cos θ = 0.94cos θ = 0.96cos θ = 0.98cos θ = 0.98 q u = - T / T m n* (b) FIG. 5: (a) Barrier height β ∆ G ∗ as a function of T /T m . Forhigh values of cos θ , the barrier decreases significantly withdecreasing temperature, while the slope is much more shal-low for low values of cos θ . The two points at lowest T forcos θ = 0 .
98 are obtained from US simulations employing amore relaxed c ij cutoff of q u = − .
65 after checking for consis-tency at higher T . (b) Critical cluster size n ∗ as a function ofthe reduced temperature. Also in this case the trend with T isstronger in narrow-patch models than in wide-patch models. the model with cos θ = 0 .
96 appears to be the optimalcandidate for crystallization from a thermodynamic per-spective, since modest supercooling is sufficient to inducebarriers of height of the order of 10 k B T .Fig. 5(b) shows n ∗ as a function of T /T m , showingsimilar trends as for β ∆ G ∗ ( T ). The cos θ = 0 .
98 and0.96 models are similar, attaining small n ∗ for a relativelysmall degree of supercooling, while the models with widerpatches seem to require larger sizes of critical nuclei as T decreases. The T dependence of n ∗ is much flatter,suggesting that the critical nucleus remains significantlylarge even under deep supercooling. D. Surface tension
To illustrate the dependence of the surface tension onpatch width, motivated by the relation 3 β ∆ G ∗ = bn ∗ / that follows from Eq. 8, we plot in Fig. 6(a) β ∆ G ∗ asa function of n ∗ / . For all models, a reasonable lineardependence is observed, with a slope that progressively ( n* ) β ∆ G * cos θ = 0.92cos θ = 0.94cos θ = 0.96cos θ = 0.98(a) T / T m b (b) T ρ σ (c) FIG. 6: (a) Barrier height as a function of n / (proportionalto the surface area of the cluster). All models show a linearbehaviour. (b) Surface coefficient b of the CNT fit as a func-tion of the reduced temperature. The surface term does notchange significantly with T , as compared to the barrier height[see Fig. 5(a)]. (c) Reduced number density ρσ as a functionof T at P = 0 .
03 for the models studied. Large open symbolsindicate the range of state points used in the calculation of b .Curves are a guide to the eye. increases on with cos θ beyond 0.94. This results in alarger critical size for a given barrier height as the patchwidth increases. The near linearity also suggests that thesurface tension does not have a strong dependence on su-percooling. In Fig. 6(b) we show b as a function of T /T m as obtained from the fitting procedure for the differentstate points studied. The surface tension increases sig-nificantly on decreasing the angular patch width. Once T β | ∆ µ | , . a cos θ = 0.92cos θ = 0.94cos θ = 0.96cos θ = 0.98 FIG. 7: CNT bulk fit parameter a as a function of temper-ature (points) compared with the difference in chemical po-tential β | ∆ µ | (lines). While the two quantities show similartrends, to achieve quantitative agreement we multiply valuesof a in the plot by a universal scaling factor of 0.4. again, however, the two models with the largest patchesare quite similar in their behavior.It is worthwhile noting that the densities at our simu-lated P for these two wide patch models (at our lowest T ) are very similar with values near ρσ = 0 .
47, almostmatching the crystal density (where ρ is the number den-sity). The density is smaller for cos θ = 0 .
96, and smallestfor cos θ = 0 .
98. Thus, the values of b for the differentmodels appear to reflect the differences in density, witha better match in density between fluid and crystal giv-ing rise to a smaller surface tension. However, a plotof the T dependence of the densities in Fig. 6(c) showsthat density does not solely determine b , as indicated bythe variation of density without a correspondingly largevariation in b within each model. E. a vs β ∆ µ In Fig. 7 we plot β | ∆ µ | [45] and a (obtained from theCNT fits) as functions of T . For comparison purposes,all the values of a appearing in the figure have been mul-tiplied by a factor of 0.4, i.e. the values of a are signif-icantly higher than expected from independent calcula-tion of β | ∆ µ | . However, it is quite comforting that therescaled a matches β ∆ µ quite well across T and θ , asthis strengthens our assumption that the fit parameter b is proportional to βγ with the same proportionality con-stant across all the models.On the other hand, it is somewhat perplexing that sucha large rescaling is required at all. In a few previous stud-ies, comparison between a and β ∆ µ have shown that insome cases a close correspondence is observed [8, 9, 41],while in principle the value of a should vary with thedefinition crystal-like particles. There are several po-tential reasons why such a disagreement can be found.First of all, the assumption of CNT may not be fullysatisfied, including those concerning the structure of the n β ∆ G ( n ) T = 0.146 T = 0.144 T = 0.142 T = 0.140 -0.9 -0.85 -0.8 -0.75 -0.7 -0.65 q u a T = 0.138 T = 0.136 FIG. 8: Comparison of the nucleation barriers obtained withtwo different order parameters, q u = − .
65 (dashed lines)and q u = − .
87 for the cos θ = 0 .
98 model. As expected, thebarrier heights are consistent (see Ref. [44]) but the criticalcluster sizes differ significantly. The more relaxed order pa-rameter leads to a larger critical cluster, possibly due to thefact that it includes more particles on the surface. In the insetwe also show the dependence of the bulk CNT term a as afunction of the order parameter, showing that it increases bymaking the order parameters more strict. crystallites and the sharpness of the interface [42]. Sec-ond, an imperfect classification of particles as solid-likeleads to an imperfect evaluation of the cluster size and aclassification-dependence of a and b .As a test of this second hypothesis and to illustratethe effect of the parameters chosen to define crystallineclusters, we plot in Fig. 8 sets of barrier profiles obtainedfor cos θ = 0 .
98 using q u = − .
65 compared against using q u = − .
87. While ∆ G ∗ remains invariant to within er-ror with respect to q u , at all investigated T , n ∗ increases,as expected, for smaller values of q u . In other words, thesize of the critical nucleus is strongly dependent on thecriteria chosen to define solid-like particles. It is not asurprise that ∆ G ∗ is a much more robust value, since thework required to produce the critical nucleus is indepen-dent on how the cluster is described. In other words, the“real” critical nucleus (i.e. configurations which 50 percent of the time will crystallize and 50 per cent will meltagain into a fluid state) is what needs to be sampled in thesimulation. The number of particles that compose thiscluster depends on the definition, affecting the perceivedsize, but not the barrier height. These conclusions are inagreement with the work of Filion et al [44] on sphericalparticles.In the inset of Fig. 8, we show a as a function of q u for T = 0 .
146 and cos θ = 0 .
98 and see that a varies byroughly a factor of 2 for the reasonable range of q u wehave explored. Coincidentally, the value for β | ∆ µ | forthis state point is 0.56. Thus, for an appropriate choiceof q u , Eqs. 7 and 8 describe the nucleation barrier profilesboth in scaling and in the numerical value of ∆ µ , but howto choose the appropriate q u a priori is not clear.The more inclusive choice of q u = − .
65 allows one tomore easily explore the regime of low barrier heights, less -0.92 -0.90-0.88-0.88-0.91-0.95 -0.10 DCDH
FIG. 9: Schematic description of the criteria that requirefour connections for defining solid-like particles. The valuesof c ij are displayed next to the bonds that the two centralparticles form. The particle on the left has three bonds with c ij < − .
87 and one bond with − . < c ij < .
1, registering itas a DH particle. The particle on the right has four bonds with c ij < − .
87, registering it as a DC particle. In the mixed four-connections case, both particles would be considered solid-like. n β ∆ G ( n ) T = 0.144 T = 0.144 pure DC T = 0.144 pure DH T = 0.144 mixed DH-DC T = 0.142 T = 0.142 mixed DH-DC FIG. 10: Nucleation barriers for pure DC, pure DH, andmixtures of the two. For T = 0 .
144 the four-connectionscurves (pure DC, pure DH and mixed DH-DC) are all quitesimilar to each other and yield the same barrier height asthe three-connections case, but smaller n ∗ . Barriers at T =0 .
142 compare the mixed three- and four-connections casesand confirm the invariance of ∆ G ∗ and the change in n ∗ withdifferent criteria for defining clusters. hampered by the formation of secondary clusters. Using q u = − .
65 we are able to add two additional points tothe curve for cos θ = 0 .
98 in Fig. 5(a).
F. DH vs DC
As we have alluded to before, the clusters are formedby a mixture of local DC and DH particles. Indeed,the criterion of having at least three connected neigh-bors ( c ij < q u ) for defining solid-like particles, accepts either a locally DC or DH particle as solid. However, thisdefinition is more likely to misidentify a DH-like particleas liquid-like, since particles in the DC crystal have foursuch connections, while DH particles only have three.The biasing potential may therefore tend to grow clustersricher in DC than would occur naturally. In the follow-ing we therefore compare this criterion with other criteriawhich enforce the growth of pure DC, pure DH or an un-biased mixture of the two. We do this at cos θ = 0 . q u = − .
87 and simply re-quire that the number of connections be exactly four (asopposed to at least three). With this more restrictivecriterion, the solid-like particles must have a DC envi-ronment. For pure DH, we require exactly three connec-tions with c ij < q u = − .
87 and one connection with − . < c ij < .
1. In this way, DC particles are excludedand can not participate in crystalline clusters. For themixed DH/DC case, a particle is solid-like if it is deter-mined again by exactly four connections, but this timethese four connection must satisfy either the criterion fora DC particle or the criterion for a DH particle. Thecartoon in Fig. 9 provides a graphical explanation of thecriteria for identifying solid-like particles.In Fig. 10, we show the barrier profile results at T = 0 .
144 with cos θ = 0 .
98 for the four different defini-tions of solid-like particles (three-connections, pure DC,pure DH and unbiased four-connections mixture), in or-der to asses if the energy barriers of the pure crystalsare different from that of the mixture and to ascertainthe difference between using three or four connections inthe mixed case. We can not discern any real differencein terms of the barrier height between the four criteria.This interesting observation supports the view for thisshort-range model, that in addition to the DH and DCfree energies being essentially identical, the surface ten-sions are also similar. Even more, if we restrict ourselvesto the comparison of all criteria in which four connec-tions are required, then not only is the barrier heightidentical, but also (within our numerical precision) the n dependence of the barrier profile. This suggests that thepathways to crystallizing DC and DH are quite similaras well. The pure DH case appears to produce a slightlylarger n ∗ , but this will also likely depend on the boundson c ij near − .
1. Comparing the three-connections cri-terion with the more restrictive four-connections cases,we see that the four-connections criterion simply resultsin an apparently smaller critical cluster size, providingone additional piece of evidence for the sensitivity of theprofile on the solid-like definition accompanied by the in-sensitivity of the barrier height.We repeat the comparison at the lower T = 0 . n , the presence of a pure DC cluster ap-pears to act as a template for DH growth (and viceversa).Essentially, the DH particles, which are not counted assolid-like, become part of the growing crystal, providinga bridge between only apparently different DC crystalclusters. Thus, the system crystallizes while still regis-tering a small largest cluster in the system and our orderparameter no longer describes nucleation.This issue is greatly reduced for the four-connectionmixed case and we are able to calculate the barrier at T = 0 . T , also for the four-connections mixed case, it be-comes impossible to properly evaluate the barrier, sincethe overly restrictive criterion will misidentify as liquid-like particles sufficiently crystalline to participate in crys-tal growth.While the barrier heights are the same across the dif-ferent cluster criteria, the DC/DH composition of theclusters varies. To quantify this in a basic way, we ana-lyze an ensemble of largest clusters extracted from a setof US simulations employing the mixed four-connectionssolid-like criterion in a wide range of n values. Wefind that 54% of the solid particles are DC (standarddeviation 17), which is the same as the composition forclusters appearing in spontaneous nucleation studies asreported in Ref. [28]. Performing the same evaluationstarting from US simulations, employing this time themixed three-connections solid-like criterion ( q u = − . q u = − .
65. This confirms that the choice of at leastthree connections with c ij close to minus one does in-deed introduce a bias toward the DC structure. For thismodel, however, this bias does not measurably affect thenucleation barrier. G. BCC
For all the patch widths considered except cos θ = 0 . θ = 0 .
98, the stable crystal phase at P = 0 . θ = 0 .
98, we evaluate the bar-rier at T = 0 .
142 using the same procedure described forthe tetrahedral crystals. The definition of solid-like BCCparticle is based on the l = 6 spherical harmonics [8],defining neighbors to be connected when the scalar prod-uct of the l = 6 harmonics is greater than 0.5. A particleis classified as a solid-like if it has at least six connections.We note that Ref. [8] does not use normalized q lm ( i ) vec-tors for calculating c ij , so we have determined cutoffsbased on our distributions of c ij and of the subsequentnumber of connections a particle has in the liquid andin the crystal. We do find that the barrier to nucleatingBCC is at least 70 k B T larger than for tetrahedral crys-tals, supporting the lack of observation of spontaneousBCC nucleation at the pressure we consider here. cos θ T X / T m T / T m β ∆ G * FIG. 11: Estimates of homogeneous nucleation temperature T (crosses) and T (circles) as functions of cos θ . Thesolid vertical line indicates the value of cos θ = 0 . δ = 0 .
24. The upper horizontal dashed line marks the ap-proximate value of T g /T m = 0 . T g = 2 T m /
3. Inset shows thedetermination of T and T from the crossing of exponen-tial fits to the four lowest T points for each value of cos θ with β ∆ G ∗ = 10 and 20, respectively. IV. SUMMARY AND CONCLUSIONS
Previous work for this tetrahedral patchy modelshowed that the driving force for nucleation (quantifiedas difference in the crystal and fluid chemical potentials)decreases as the patch width increases. Fig. 7 summa-rizes this result by showing β | ∆ µ | increasing rapidly tolarge values as T decreases below T m for narrow patches,while increasing only slowly and remaining fairly smallfor wide patches.We find that despite this reduced driving force, thebarriers to nucleation actually are smallest for the cos θ =0 .
96 case [as shown in Fig. 5(a)], in the sense thatthis model reaches the homogeneous nucleation limit,marked for example by the barrier reaching a value of β ∆ G ∗ = 10, at the highest T /T m amongst the stud-ied models. The increased similarity between liquid andcrystal in terms of energy and density as patch width in-creases not only brings the chemical potentials of liquidand crystal closer in value (tending to increase the nucle-ation barrier), but also reduces the surface tension (tend-ing to lower the barrier). Thus, in the range of narrowangles where crystallization is readily observed, compe-tition between | ∆ µ | and γ leads to optimal nucleation atan intermediate patch width.Increasing the patch width beyond cos θ = 0 .
94 nolonger significantly reduces γ , while | ∆ µ | continues to de-crease, causing an increase in nucleation barrier heights.For cos θ = 0 .
92, this resulting increase is quite large,with β ∆ G ∗ estimated to be of order 50 at the lowest T that we can simulate. Moreover, the rate of decrease of β ∆ G ∗ with T appears to be quite slow, requiring signif-icant supercooling to reach accessible nucleation barrier0heights.The evaluation of the barriers requires a definition ofsolid-like particles. We have checked that the barrierheight is essentially insensitive to the exact choice ofthe cutoff used to define solid-like connections, consis-tent with the observation in Ref. [44] for hard-spheres.The size of the critical nucleus is instead significantly de-pendent on the definition of solid-like particles, again inagreement with the conclusions in Ref. [44]. We havealso compared several definitions of solid-like particles toprobe the crystallization of a pure DC, a pure DH and amixed DC-DH structure. Interestingly, we find the samevalue of β ∆ G ∗ for all structures and, in this case, similar n ∗ . The fact that the entropy gain in creating a mixedstructure does not appreciable lower the work of forminga critical nucleus, but is sufficient to generate a preva-lence of mixed structures when spontaneous nucleationtakes place [28], warrants further investigation.At this point, we would be remiss if we were not tocomment on the interplay between dynamics and ther-modynamics in controlling the rate of nucleation, andindeed governing the glass-forming abilities of the sys-tem. We already see an example of this within our data.The barriers at the lowest T studied for cos θ = 0 .
98 andcos θ = 0 .
94 are the same within error, and have a value of β ∆ G ∗ ≈
17. However, as mentioned previously, the dy-namics (in terms of the diffusion coefficient) are 40 timesslower for the wider patch case. For simulations, this isa significant number. The slow decrease in β ∆ G ∗ with T , combined with the expected slow-down in dynamics,suggest that a search for spontaneous nucleation in un-constrained simulations would target T not far from thelowest T for which we have calculated the barrier. Thecase of cos θ = 0 .
92 is more obvious since its dynamicsare even slower at a given T than those of cos θ = 0 . β ∆ G ∗ to find estimates for the temperature atwhich the homogeneous nucleation limit is reached, T X .We define two such estimates using β ∆ G ( T ) = 10 and β ∆ G ( T ) = 20, and determine them by fitting the low-est four points in T from Fig. 5(a) for each model toexponentials and extrapolate where necessary. In Fig. 11we present the resulting estimates of T X alongside a verti-cal line showing the widest patch that guarantees a singlebond per patch, i.e. four bonds per particle. Also shownare two horizontal lines indicating the approximate valueof T g /T m = 0 .
80 for silica [46] (where T g is the glass tran-sition temperature), as well as the general rule of thumb T g /T m = 2 / T X /T m below T g /T m = 2 /
3, a model may be consid-ered a glass-former; our estimates suggest that we maybe approaching this regime with our widest patch model.The presence of a maximum in T X confirms that decreas-ing the angular range for bonding favors glass formation.Tetrahedral colloidal particles will thus form crystals onlyif the bonding angular width is small.As a final remark, we point out that our findingsmay shed some light on the glass-forming and crystalliz-ing abilities of molecular or atomic tetrahedral network-forming liquids, contributing insight as to why, for exam-ple, water crystallizes while silica more readily forms aglass. V. ACKNOWLEDGMENTS
We thank ACEnet (Canada) for computational re-sources. I.S.-V. acknowledges financial support fromNSERC (Canada). FS and FR acknowledge sup-port from ERC-226207-PATCHYCOLLOIDS and ITN-COMPLOIDS. [1] D. Kashchiev,
Nucleation (Butterworth-Heinemann, Ox-ford, 2000).[2] K.F. Kelton and A.L. Greer,
Nucleation in CondensedMatter: Applications in materials and biology (Elsevier,2010).[3] P.G. Vekilov, Soft Matter , 5254 (2010).[4] C. Vega, E. Sanz, J.L.F. Abascal, and E.G. Noya, J.Phys.: Condens. Matter , 153101 (2008).[5] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J.Chem. Phys. , 9932 (1996).[6] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J.Chem. Phys. , 1591 (1999).[7] S. Auer and D. Frenkel, Nature (London) , 1020(2001).[8] S. Auer and D. Frenkel, J. Chem. Phys. , 3015 (2004).[9] C. Valeriani, E. Sanz, and D. Frenkel, J. Chem. Phys. , 194501 (2005). [10] J.W. Gibbs, The Scientific Papers of J. Willard Gibbs (Dover, New York, 1961).[11] M. Volmer and A. Weber, Z. Phys. Chem. , 227(1926).[12] L. Farkas, Z. Phys. Chem. , 236 (1927).[13] R. Becker and W. D¨oring, Ann. Phys. , 719 (1935).[14] K.F. Kelton, Crystal Nucleation in Liquids and Glasses (Academic, Boston, 1991), Vol. 45, pp. 75-177.[15] P.G. Debenedetti,
Metastable Liquids. Concepts andPrinciples (Princeton University Press, Princeton, NewJersey, 1996).[16] M. Matsumoto, S. Saito, and Iwao Ohmine, Nature ,409 (2002).[17] M. Yamada, S. Mossa, H.E. Stanley, and F. Sciortino,Phys. Rev. Lett. , 195701 (2002).[18] T. Motooka and S. Munetoh, Phys. Rev. B , 073307(2004). [19] E. Sanz, C. Vega, J.L.F. Abascal, and L.G. MacDowell,Phys. Rev. Lett. , 255701 (2004)[20] P. Beaucage and N. Mousseau, Phys. Rev. B , 094102(2005).[21] E.G. Noya, C. Vega, J.P.K. Doye, and A.A. Louis, J.Chem. Phys. , 234511 (2010).[22] L.M. Ghiringhelli, C. Valeriani, E.J. Meijer and D.Frenkel, Phys. Rev. Lett. , 055702 (2007).[23] L.M. Ghiringhelli, C. Valeriani, J.H. Los, E.J. Meijer, A.Fasolino, and D. Frenkel, Mol. Phys. , 2011 (2008).[24] F.H. Stillinger and T.A. Weber, Phys Rev B , 5262(1985).[25] V. Molinero, S. Sastry, and C.A. Angell, Phys. Rev. Lett. , 075701 (2006).[26] T. Li, D. Donadio, L.M. Ghiringhelli, and G. Galli, Na-ture Mater. , 726 (2009).[27] M. Maldovan and E.L. Thomas,. Nature Mater. , 593(2004).[28] F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys. , 174502 (2011).[29] N. Kern and D. Frenkel, J. Chem. Phys. , 9882(2003).[30] W. Ostwald, Z. Phys. Chem. , 289 (1897).[31] R.A. van Santen, J. Phys. Chem. , 5768 (1984).[32] P.R. ten Wolde and D. Frenkel, Phys. Chem. Chem.Phys. , 2191 (1999).[33] D. Frenkel and B. Smit, Understanding Molecular Sim-ulation: From Algorithms to Applications (AcademicPress, San Diego, 1996).[34] F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys. , 184501 (2010).[35] P.J. Steinhardt, D.R. Nelson, and M. Ronchetti, Phys.Rev. B 28, 784 (1983). [36] J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. , 134103 (2007).[37] J. Wedekind and D. Reguera, 17th International Confer-ence on Kinetic Reconstruction of the Nucleation FreeEnergy Landscape in Nucleation and Atmoshperic Ae-orsols, Galway, Ireland, 2007 (Springer, Netherlands,2008).[38] J. Wedekind and D. Reguera, J. Phys. Chem. B ,11060 (2008).[39] G. Chkonia, J. Wlk, R. Strey, J. Wedekind, and D.Reguera, J. Phys. Chem. B , 064505 (2009).[40] S.E.M. Lundrigan and I. Saika-Voivod, J. Chem. Phys. , 104503 (2009).[41] I. Saika-Voivod, Peter H. Poole, and R.K. Bowles, J.Chem. Phys. , 224709 (2006).[42] W. Lechner, C. Dellago, and P.G. Bolhuis, Phys. rev.Lett. , 085701 (2011).[43] S. Ryu and W. Cai, Phys. Rev. E , 030601(R) (2010).[44] Laura Filion, M. Hermes, R. Ni, and M. Dijkstra, J.Chem. Phys. , 244115 (2010).[45] Here, ∆ µ refers specifically to the difference in chemi-cal potential between the DC crystal and the fluid, ascalculated in Ref. [28], even in the case of cos θ = 0 . T m is determined fromthe condition ∆ µ = 0.[46] We have used T g = 1450 K and T m = 1823 K for quartz.[47] P.G. Debenedetti and F.H. Stillinger, Nature , 259(2001).[48] S. Sakka and J. D. Mackenzie, J. Non-Crys. Solids2