Role of bond orientational order in the crystallization of hard spheres
With computer simulations of the hard sphere model, we examine in detail the microscopic pathway connecting the metastable melt to the emergence of crystalline clusters. In particular we will show that the nucleation of the solid phase does not follow a two-step mechanism, where crystals form inside dense precursor regions. On the contrary, we will show that nucleation is driven by fluctuations of orientational order, and not by the density fluctuations. By considering the development of the pair-excess entropy inside crystalline nuclei, we confirm that orientational order precedes positional order. These results are at odd with the idea of a two-step nucleation mechanism for fluids without a metastable liquid-liquid phase separation. Our study suggests the pivotal role of bond orientational ordering in triggering crystal nucleation.
After almost 150 years from the pioneering works of Gibbs, a satisfactory understanding of the crystallization process is still lacking. Classical Nucleation Theory (CNT) describes nucleation as an activated process in which the free energy gain of the solid phase over the fluid one is in competition with the free energy penalty of forming an interface for nuclei of size smaller than a characteristic length, the critical nucleus size. The two main limitations of CNT have been recognized as the following: i) the capillarity approximation, i.e. the assumption that the thermodynamic properties of the small crystal nuclei are the same as the bulk crystal, such as specific volume or surface tension; ii) all the order parameters involved in the transition proceed simultaneously, and the transition can be described effectively by one order parameter (as density for the nucleation of liquid droplets from a supersaturated gas). Many efforts have been done in order to relax some of the assumptions of CNT, showing in fact that nuclei do not form at bulk condition and that the process is non-classical, i.e. the order parameters involved in the transition do not change simultaneously Oxtoby (2003).
On these grounds, the two-step nucleation mechanism has emerged as a candidate description for many of the discrepancies between CNT and experiments. The idea is the following. The transition from the fluid phase to the solid phase involves at least two order parameters, one associated with the breaking of translational symmetry and the other with the breaking of rotational symmetry. In the crystallization of spherical particles, usually these order parameters are identified with density and bond orientational order. In the two-step mechanism scenario, the first fluctuation to trigger the transition is density fluctuation, which leads to the formation of a dense droplet in the melt. Then the droplet starts to be structured and a crystal is formed. Large-amplitude density fluctuations can emerge from the criticality of a nearby liquid-liquid (L-L) phase transition, which is the case of globular proteins Galkin and Vekilov (2000) or colloids with small range attractions Savage and Dinsmore (2009). This also led to the speculation that critical fluctuations could enhance the crystallization rate ten Wolde and Frenkel (1997); Xu et al. (2012). Rather surprisingly, the same mechanism was claimed for the nucleation pathway of many systems for which a dense liquid phase is not known Garetz et al. (2002); Pouget et al. (2009) or even does not exist at all Schöpe et al. (2006); Schilling et al. (2010).
Recently we have provided a detailed microscopic study of the nucleation pathway in ultrasoft potential systems Russo and Tanaka (2012a) and in hard sphere systems Russo and Tanaka (2012b), both having purely repulsive potentials and thus without a L-L transition. These works highlighted the role played by bond orientational order in the crystallization process Kawasaki and Tanaka (2010), finding no evidence of the dense precursors upon which the two-step assumption is built. On the contrary, the transition resembles a two-step process where the first step is the formation of extended structured regions of high orientational order, which progressively densify. The bulk density is in fact reached only at a very late stage, when the nucleus becomes several times the critical nucleus size. The first step, structuring at (almost) constant density is achieved due to the wetting of bond-orientational regions to small crystalline nuclei. The identification of the bond orientational order as the appropriate coordinate to describe the nucleation stage has unveiled that the polymorph selection stage starts already in the metastable fluid phase before nucleation occurs Russo and Tanaka (2012a, b). With bond orientational order we have also studied the competition between crystalline structures and icosahedral (or five-fold symmetric) structures Russo and Tanaka (2012b); Leocmach et al. (2012), i.e. the mechanism by which crystallization is prevented and that paves the way to out-of-equilibrium states (i.e. glasses).
In the present contribution we will focus on the nucleation pathway in monodisperse hard spheres. We will show that, while density fluctuations have a very short lengthscale, bond orientational order has an extended lengthscale that increases with supersaturation. The nucleation events are localized inside the regions of high bond orientational order. The process is thus non-classical, with orientational order preceding translational order as the nuclei grow.
We run Umbrella Sampling Monte Carlo simulations of monodisperse hard spheres in the isothermal-isobaric (NpT) ensemble. Hard spheres are ideal for studying crystallization and have already provided tremendous contributions to our basic understanding of crystal nucleation Gasser et al. (2001); Zaccarelli et al. (2009); Schilling et al. (2010); Sanz et al. (2011); Valeriani et al. (2012); Taffs et al. (2012). In the following, lengths are given in unit of the diameter , and pressure in units of , where is the Boltzmann constant. We place the spheres randomly in a simulation box at packing fraction and equilibrate the system at reduced pressure . At this pressure the liquid is metastable with respect to crystallization, with a difference in chemical potential between the liquid and solid state of Filion et al. (2010). We run Umbrella Sampling simulations, where a biasing potential is added to the system Hamiltonian to sample crystalline clusters of different sizes, in order to extract information on the free energy barrier and the critical cluster size. Details on the implementation can be found in Ref. Auer and Frenkel (2005). Configurations for clusters at different sizes are extracted from the simulations and analyzed. The Umbrella Sampling biasing is used to gain better statistics by stabilizing nuclei at the desired size. But we confirm that the same results are obtained if the configurations are extracted from direct simulations (still possible at this pressure) without biasing potential. At the free energy barrier is and the size of the critical nucleus is .
The identification of crystalline particles follows the usual procedure Auer and Frenkel (2004). A particle is identified as crystal if its orientational order is coherent (in symmetry and in orientation) with that of its neighbors. The -fold symmetry of a neighborhood around each particle is characterized by a dimensional complex vector () as , where is a free integer parameter, and is an integer that runs from to . The functions are the spherical harmonics and is the vector from particle to particle . The sum goes over all neighboring particles of particle . Usually is defined by all particles within a cutoff distance, but in an inhomogeneous system the cutoff distance would have to change according to the local density. Instead we fix which is the number of nearest neighbors in icosahedra and close packed crystals (like hcp and fcc) which are known to be the only relevant crystalline structures for hard spheres. The scalar product quantifies the similarity of the two environments. If it exceeds between two neighbors, they are deemed connected. We then identify a particle as crystalline if it is connected with at least neighbors Auer and Frenkel (2004).
The liquid-to-solid transition is characterized by the symmetry breaking of both translational and orientational order Truskett et al. (2000). Translational order expresses the relative spacing between particles in the system, as shown in the top row of Fig. 1. A good order parameter for translational order is the local density , which is calculated by constructing the Voronoi diagram and measuring the volume of the cell associated with each particle. More generally, translational order can be obtained from two-body correlation functions, such as the pair correlation function Truskett et al. (2000). Since we are dealing with hard sphere systems, where there is no energy involved, the most meaningful definition of translational order comes from the two-body excess entropy, defined as:
Orientational order, instead, expresses the relative orientation between the neighbors around each particle, as shown in the bottom row of Fig. 1. Unlike translational order, which is obtained from two-body correlation functions, orientational order is obtained by considering many-body correlations. For hard spheres it is well known that the local bond-order analysis introduced by Steinhardt Steinhardt et al. (1983) is an adequate measure of orientational order. This order parameter is obtained by constructing a rotational invariant from the quantity previously introduced. In the present work we are going to focus on the coarse-grained quantity , which can be obtained in the following way. First is spatially coarse-grained Lechner and Dellago (2008)
where are the neighbors of particle . Then the following invariant is obtained
where we focus on the case (corresponding to the six-fold symmetry of the crystal structures). The effect of coarse graining has proven effective not only in reducing the signal-to-noise ratio Lechner and Dellago (2008), but also eliminates the signal from five-fold symmetric structures Russo and Tanaka (2012b); Leocmach and Tanaka (2012) which do not add coherently. In this way the order parameter goes continuously from low values in the fluid phase to high values in the crystal phase.
We start by looking at the two-dimensional probability distribution of density and bond orientational order , for a metastable fluid state at pressure (before the appearance of the critical nucleus), reported in Fig. 2. The probability distribution is related to the Landau free energy, . The free energy can be well fitted with a full cubic polynomial, for which the most important term is of the form . This term is responsible for the shape of contours lines (dashed line in Fig. 2): because the interaction is linear in and quadratic in , the system can increase its orientational order without an increase of its translational order, but the contrary is not true, and an increase in density also increases the average . A similar term was proposed in order to describe quasi-crystals formation, where the ordering of the orientational field does not imply an ordering of the translational field Jarić (1985). Note also that a small linear coupling between and exists at high , seen for example in the small slope of the steepest descent path (white dashed arrow) in Fig. 2. As we will discuss later, we believe this small coupling term to be responsible for previous observations of two-step nucleation processes in hard spheres.
Next we address the correlation length of these order parameters. In the two-step scenario, crystal nucleation occurs inside dense droplets that form because of fluctuations in the density field. In a recent study we have compared correlation lengths for both translational and orientational order Leocmach et al. (2012) and found that, while the correlation length associated with density (and any two-body quantity) is always very short and does not grow with increasing supersaturation, the correlation length of orientational order is an increasing function of supersaturation. This already suggests that the mechanism of formation of dense droplets is unlikely in hard spheres, while a viable mechanism would be the nucleation inside localized regions of high orientational order Russo and Tanaka (2012b); Kawasaki and Tanaka (2010). We show this in Fig. 3, by simply looking at a typical configuration in the metastable state. On the left panel we plot a subset of particles ( of the total) having the highest value of . It is immediately clear that the position of particles is approximately uncorrelated, without any apparent lengthscale. The colors are given according to the bond orientational order field, with high particles being colored in red (dark gray), and low particles colored in green (light gray). We see that there is a lack of correlation of particles in dense environments and the value of . Going back to Fig. 2 this is confirmed by observing that, at high density, the value (while increasing on average) is equally likely to have low or high values. On the right panel of Fig. 3 we plot instead the subset of particles ( of the total) having the highest value of . Unlike the density field, bond orientational order clearly displays a strong correlation between the particles, with structures resembling droplets. In this case, the colors are given according to the value of of each particle: red (dark gray) particles having a high value of and green (light gray) particles a low value of . This time there is a clear correlation between high values of and high values of , so clusters of high orientational order will appear denser. This correlation comes from the linear term that we noted in the distribution function of Fig. 1.
As we will show soon, since crystals are born from fluctuations of the orientational order field, a measure of the density of the regions from which the crystal originates would give a value higher than the average melt density. But this does not mean that the crystal originated due to a fluctuation of the density field, but it is simply a consequence of this small linear coupling between and .
We now proceed to study the nucleation stage, where a crystal is nucleated up to the critical size, eventually overcoming the free energy barrier and starting the crystallization process. We first focus on nuclei of a size comparable to the critical nucleus size (so we average over many independent configurations where the biggest nucleus is of size between and particles). In Fig. 4 we plot the radial profile for three different fields, , and the two-body excess entropy , as a function of the distance from the center of the nucleus (). The fields are normalized as to be in the fluid phase and in the bulk crystal phase. First we note that the process is non-classical, where all fields have still not reached their bulk value, and the radial distribution is different for each field. In particular, going from the fluid phase (high value of ) to the center of the crystal nucleus (), the field is the first to develop, while density is lagging behind the development of the orientational field.
This again proves that nucleation is driven by orientational order and not by density. The same result is confirmed by looking at other definitions of translational order, which can be extracted from two-body correlation functions. For hard spheres systems, where there is no energy term in the free energy, a physically motivated definition of translational order is given by the two-body excess entropy defined in Eq. 1. To obtain the radial dependence of we introduce the functions which measure the radial distribution function averaged over particles at distance from the critical nucleus. In formula
where the ensemble average is over Umbrella Sampling configurations with a crystal nucleus of critical size. The two-body excess entropy function is then
where the integral is carried up to a distance corresponding to the second nearest-neighbors shell (as nuclei are still small). The result shown in Fig. 4 clearly indicates that the translational order (expressed by ) is lagging behind orientational order (expressed by ).
The slow positional ordering inside the critical nucleus can be seen in Fig. 5, where the pair correlation function is plotted at different distances () from the center of the critical nucleus. The at the center is the one closest to the dashed line, which is the pair correlation function of a bulk crystal. Even at the center of the critical nucleus (), the pair correlation function is far from that of a bulk crystal.
Already from Fig. 4 it is clear that the two-step scenario with dense precursors is unlikely, as orientational order precedes the densification process. To have a direct confirmation of the absence of a two-step process, in Fig. 6 we plot the density at the center of the nucleus as a function of the nucleus size. The increase in density is linear with the size of the nucleus, and extrapolates to the fluid density in the limit of vanishing size of the nucleus. This result is in contrast to the expectation of the two-step nucleation mechanism, which instead predicts a non-linear dependency of the density on the nucleus size: first an enrichment at constant size, followed by a size growth with little density change. Figure 6 thus confirms that there is nothing anomalous with the density increase during the nucleation process.
In the present article we have examined in detail the possibility of a two-step nucleation process from the supercooled state in hard spheres. The two-step nucleation is based on a specific behavior of two order parameters: density and structure. The fluid phase first forms a dense precursor, which is then structured to form a crystalline nucleus. Nucleation is then initiated by density fluctuations, and followed by the superposition of a structural fluctuation. Unlike this scenario, we stress that hard spheres take another route to crystal nucleation. Nucleation is promoted by extensive fluctuations in bond orientational order. These fluctuations act as precursors for the formation of crystalline nuclei. The increase of density inside these regions lags behind the development of orientational order, so that one can effectively have small crystalline nuclei with densities still very far from bulk conditions. The reason why the transition from a fluid structure to a solid structure occurs with little density change is due to the fact that nuclei form in regions of high orientational order, and so are wetted by fluid regions with a high value of this field. The structuring of the nuclei then requires just small adjustments of particle positions, without big density changes Russo and Tanaka (2012b).
This decoupling between the orientational order parameter and the translational order parameter has also important consequences for the glass transition. Once the positional order transition is avoided, orientational order keeps growing with increasing supersaturation. It has been recently shown that this increase of orientational order is linked with the slowing down of the dynamics, displaying critical-like behavior Tanaka et al. (2010); Leocmach and Tanaka (2012).
Interestingly, the mechanism by which the positional ordering transition is avoided seems to be also correlated with the development of orientational order. We have recently shown Russo and Tanaka (2012b) that crystalline structures are not the only structures that are stabilized by an increase of orientational order. There is in fact a competition between crystalline structures and icosahedral structures (five-fold symmetric structures), and we have shown that both an increase in density and polydispersity Leocmach et al. (2012) stabilize the icosahedral structure, effectively suppressing crystallization. This again confirms the effect of frustration on crystallization Tanaka (1998), as in the case of 2D spin liquids Shintani and Tanaka (2006, 2008), and as recently observed in metallic glasses Jakse and Pasturel (2008); Hwang et al. (2012).
We stress that the two-step scenario is instead fully consistent for systems which display a L-L transition. In this case, the presence of a minima (even if metastable) in the system’s free energy corresponding to the liquid phase can start the nucleation process without a macroscopic phase separation ten Wolde and Frenkel (1997); Galkin and Vekilov (2000); Savage and Dinsmore (2009); Desgranges and Delhommelle (2011); Murata and Tanaka (2012).
This study was partly supported by a grant-in-aid from the Ministry of Education, Culture, Sports, Science and Technology, Japan (Kakenhi) and by the Japan Society for the Promotion of Science (JSPS) through its “Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST Program)” and a JSPS Postdoctoral Fellowship.
- D. Oxtoby, Philos. T. R. Soc. A 361, 419–428 (2003).
- O. Galkin, and P. Vekilov, Proc. Nat. Acad. Sci. U.S.A. 97, 6277–6281 (2000).
- J. R. Savage, and A. D. Dinsmore, Phys. Rev. Lett. 102, 198302 (2009).
- P. ten Wolde, and D. Frenkel, Science 277, 1975–1978 (1997).
- L. Xu, S. V. Buldyrev, H. E. Stanley, and G. Franzese, Phys. Rev. Lett. 109, 095702 (2012).
- B. A. Garetz, J. Matic, and A. S. Myerson, Phys. Rev. Lett. 89, 175501 (2002).
- E. Pouget, P. Bomans, J. Goos, P. Frederik, G. de With, and N. Sommerdijk, Science 323, 1455–1458 (2009).
- H. Schöpe, G. Bryant, and W. van Megen, Phys. Rev. Lett. 96, 175701 (2006).
- T. Schilling, H. J. Schöpe, M. Oettel, G. Opletal, and I. Snook, Phys. Rev. Lett. 105, 025701 (2010).
- J. Russo, and H. Tanaka, Soft Matter 8, 4206–4215 (2012a).
- J. Russo, and H. Tanaka, Scientific Reports 2 (2012b).
- T. Kawasaki, and H. Tanaka, Proc. Nat. Acad. Sci. U.S.A. 107, 14036 (2010), ISSN 0027-8424.
- M. Leocmach, J. Russo, and H. Tanaka, Importance of many-body correlations in glass transition: an example from polydisperse hard spheres. Accepted in J. Chem. Phys. (2012).
- U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey, and D. A. Weitz, Science 292, 258 (2001).
- E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates, and P. N. Pusey, Phys. Rev. Lett. 103, 135704 (2009).
- E. Sanz, C. Valeriani, E. Zaccarelli, W. C. K. Poon, P. N. Pusey, and M. E. Cates, Phys. Rev. Lett. 106, 215701 (2011).
- C. Valeriani, E. Sanz, P. Pusey, W. Poon, M. Cates, and E. Zaccarelli, Soft Matter 8, 4960 (2012).
- J. Taffs, S. Williams, H. Tanaka, and C. Royall, arXiv preprint arXiv:1206.5526 (2012).
- L. Filion, M. Hermes, R. Ni, and M. Dijkstra, J. Chem. Phys. 133, 244115 (2010).
- S. Auer, and D. Frenkel, Adv. Polym. Sci. 173, 149–207 (2005).
- S. Auer, and D. Frenkel, J. Chem. Phys. 120, 3015–3029 (2004).
- T. M. Truskett, S. Torquato, and P. G. Debenedetti, Phys. Rev. E 62, 993 (2000).
- P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784–805 (1983).
- W. Lechner, and C. Dellago, J. Chem. Phys. 129, 114707 (2008).
- M. Leocmach, and H. Tanaka, Nature Commun. 3, 974 (2012), ISSN 2041-1723.
- M. V. Jarić, Phys. Rev. Lett. 55, 607–610 (1985).
- H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe, Nature Mater. 9, 324–331 (2010).
- H. Tanaka, J. Phys.: Condens. Matter 10, L207–L214 (1998).
- H. Shintani, and H. Tanaka, Nature Phys. 2, 200–206 (2006).
- H. Shintani, and H. Tanaka, Nature Mater. 7, 870–877 (2008).
- N. Jakse, and A. Pasturel, Phys. Rev. B 78, 214204 (2008).
- J. Hwang, Z. Melgarejo, Y. Kalay, I. Kalay, M. Kramer, D. Stone, and P. Voyles, Phys. Rev. Lett. 108, 195505 (2012), ISSN 0031-9007.
- C. Desgranges, and J. Delhommelle, J. Am. Chem. Soc. 133, 2872–2874 (2011).
- K. Murata, and H. Tanaka, Nature Mater. 11, 436–443 (2012).