Published on

28-Jul-2016View

214Download

2

Transcript

ARTICLESFrustration on the way to crystallization inglassHIROSHI SHINTANI AND HAJIME TANAKA*Institute of Industrial Science, University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan*e-mail: tanaka@iis.u-tokyo.ac.jpPublished online: 19 February 2006; doi:10.1038/nphys235Some liquids do not crystallize below the melting point,but instead enter into a supercooled state and oncooling eventually become a glass at the glass-transitiontemperature. During this process, the liquid dynamics notonly drastically slow down, but also become progressivelymore heterogeneous. The relationship between the kineticslowing down and growing dynamic heterogeneity is a keyproblem of the liquidglass transition. Here, we study thisproblem by using a liquid model, with a crystalline groundstate, for which we can systematically control frustrationagainst crystallization. We found that slow regions havinga high degree of crystalline order emerge below themelting point, and their characteristic size and lifetimeincrease steeply on cooling. These crystalline regions leadto dynamic heterogeneity, suggesting a connection to thecomplex free-energy landscape and the resulting slowdynamics. These ndings point towards an intrinsic linkbetween the glass transition and crystallization.The mechanism of the drastic slowing down of the structuralrelaxation of a liquid on cooling is one of the central issues ofthe physics of the liquidglass transition. In relation to this,experimental14 and numerical5,6 evidence has been accumulatedfor the existence of dynamic heterogeneity in a supercooled stateof liquid, whose length scale increases on cooling. The link betweenthis growing length scale and the slowing down of the structuralrelaxation has been actively discussed on the basis of the AdamGibbs theory7, the spinglass model8,9, the mode-coupling theory10and the kinetically constrained lattice model11, as it may providea clue to our understanding of the liquidglass transition. Despiteintensive eorts, however, the physical factors that control dynamicheterogeneity remain elusive.Many dierent ideas have been proposed on the origin of theliquidglass transition112. Among them, frustration has often beenconsidered to play a key role in the glass transition in conjunctionwith other glassy systems such as spin glass and dipolar glass12.The most popular scenario is geometrical frustration associatedwith the fact that locally favoured structures such as icosahedralorder cannot ll up the space1216. In other words, liquids tendtowards global icosahedral order, but they cannot achieve it,which is the origin of frustration in this scenario. There is also adierent possibility17,18 (see ref. 19 on the dierence in the physicalmechanism between the two scenarios): liquids tend to order intothe equilibrium crystal, but frustration eects of locally favouredshort-range ordering on long-range crystalline ordering preventcrystallization and help vitrication. Even a simple one-componentliquid may suer from such frustration if the local symmetry of theinteraction potential does not perfectly match the symmetry of theequilibrium crystal.Motivated by this frustration scenario, we have developed adierent type of simulation model, where we can systematicallychange the degree of frustration against crystallization. Moleculardynamics simulation is a useful means for investigating thedynamics of liquids microscopically6. We can also directly controlthe intermolecular potential in molecular dynamics, whereas itis dicult to do so in experiments. We modify a sphericallysymmetric interaction potential by including an anisotropic part;more specically, we introduce a spin director ui into particle isuch that the interaction between particle i and j depends not onlyon the interparticle distance |rij|, but also on the angles among200 nature physics VOL 2 MARCH 2006 www.nature.com/naturephysicsUntitled-3 1 2/18/06, 11:26:27 AMNature Publishing Group 2006ARTICLES = 0.2 = 0.6= 0.5 = 0.6= 0.65= 0.4= 0.5= 0.6LiquidAFMcrystal AFMcrystalPlasticcrystal1.11.21.31.4 3.5 2.5 3.0 4.0 4.0 3.8 3.6 3.4 3.2 4.0 4.2 3.8 3.6 3.4 3.2 3.00.10.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.60.2 0.3 0.1 0.2 0.3VUUUT TT Ta bc dFigure 1 Phase change on cooling and heating. a,b, T-dependence of V (a)and U (b) for various on cooling. The arrow indicates the location of Tg. For 0.5, we see a step in the volume change, reecting crystallization into a plasticcrystal. For= 0.6, on the other hand, there is no step and only the slope ofdV /dT changes at Tg. c,d, T-dependence of U on heating are shown for= 0.5 (c) and= 0.6 (d). For= 0.5, we see two steps, reecting theAFM-to-plastic crystal transformation and the melting of the plastic crystal. For= 0.6, on the other hand, we see only one step, which reects the melting of theAFM crystal.ui, uj , and rij (see the Methods section). The modulation strengthis given by a parameter . This potential locally favours ve-fold symmetry (see Supplementary Information, Fig. S1), but hasa crystalline state with long-range antiferromagnetic (AFM) spinorder as the ground state. As ve-fold symmetry is not consistentwith any crystallographic symmetry13,14, this can be regarded asthe strength of frustration against crystallization. To the best of ourknowledge, this is the rst two-dimensional glass-forming modelliquid system that consists of single-component particles and has acrystalline state as the ground state. By using this lattice-free spin-glass model, we investigate the roles of frustration in vitrication aswell as the relationship between crystallization and vitrication.We investigated the phase behaviour of this system as a functionof by studying the temperature (T) dependence of the averagevolume per particle V and the average potential energy perparticle U at pressure P = 0.5 for a xed cooling rate Q = 104(see Fig. 1a,b). All results are given in reduced units (see Methods).For weak frustration (that is, for < 0.6), the system crystallizes(into a plastic crystal) on cooling with a distinct jump in thevolume. On the other hand, for strong frustration (that is, for 0.6), there is no jump in the volume itself, but there arerather sharp changes in both dV /dT and dU /dT aroundthe glass-transition temperature (Tg). This indicates that thesystem is vitried without crystallization under strong frustrationeven for the same cooling rate Q: energetic frustration againstPlastic crystalAFMcrystalGlassLiquidT Supercooledliquid0.00.0 0.2 0.4 0.6 0.8 1.00.10.20.30.40.5Figure 2 The phase diagram of our model system. The horizontal axis represents, which is a measure of the strength of frustration against crystallization. The lledcircles represent the melting point of a plastic crystal. The open squares representthe temperature where an AFM crystal transforms into a plastic crystal or melts onheating. The lled squares represent the Tg, which is dened as (Tg) = 106. Thelled triangles represents the VogelFulcher temperature T0. T0 = 0.099, 0.076 and0.026 for= 0.6, 0.7 and 0.8, respectively. The increase in Tg T0 withmeans that a liquid becomes stronger (that is, less fragile) with an increase in(see the text). We also show a simulated structure corresponding to each state ofour system. The arrows indicate the spin direction.crystallization prevents crystal nucleation. To understand whatkinds of crystalline ordering take place in this system, we studiedcrystalline structures that can be formed at P = 0.5 for various. We found that the lowest-energy ground state of our systemat > 0 is a crystal with AFM spin order (see Fig. 1c). Thiscrystal has translational symmetry concerning both position andorientation of particles. Note that this crystalline lattice is notperfectly hexagonal, but slightly distorted uniaxially because ofstrong anisotropic interactions. We also studied the change ofthe potential energy on heating to work out the phase behaviourfor various values of (see Fig. 1c,d). For weak frustration (for< 0.6), the crystal with AFM spin order rst transforms into aplastic crystal to gain the rotational entropy, and then the plasticcrystal melts into a liquid at Tm(), where Tm is the meltingpoint. This plastic crystal has hexagonal positional order but noorientational order. For strong frustration (for 0.6), on theother hand, no plastic crystal exists and the crystal with AFM spinorder melts directly into a liquid. Interestingly, in the region where,on heating, the AFM crystal directly melts into a liquid withoutpassing through the plastic crystal (namely, for 0.6), a liquidvitries without crystallization on cooling. Figure 2 summarizes thephase behaviour of this system including the non-equilibrium stateas a function of the degree of frustration.We calculated the T-dependence of the radial distributionfunction g(r). Figure 3a shows such an example for = 0.6.There is no long-range correlation even at low temperatures, whichindicates that the system forms a glassy state for this value of. Thus, we conrm that a liquid vitries on cooling withoutcrystallization for 0.6. As explained in the Methods section,our system tends to form locally favoured structures of ve-foldsymmetry. From the T-dependence of g(r), we can clearly see thatboth clusters with crystallographic symmetry and locally favouredstructures progressively develop with a decrease in T . We conrmedthat the T-dependence of the height of a peak corresponding to thenature physics VOL 2 MARCH 2006 www.nature.com/naturephysics 201Untitled-3 2 2/18/06, 11:26:28 AMNature Publishing Group 2006ARTICLES1 1.5 2 2.5T*m T*m1/T 1/TLogg (r)g (r)024685012341.00.50.60.70.80.91 2 3 4 5 61 2 3 4 5 6r012340 5 10 15rab cFigure 3 Change in structures and dynamics on cooling. a, T-dependence of theradial distribution function g(r ) for= 0.6. T is 0.18, 0.22, 0.25 and 0.4 from thetop to the bottom. g(r ) in a wider r range for T = 0.18 is shown in the inset. Thearrows at r 1.13 and r 1.63 correspond to a peak of the AFM crystal and apeak of the locally favoured structure of ve-fold symmetry, respectively (see theblack lines in the diagrams of the two structures). This can also be conrmed in realspace in Fig. 4. b, T-dependence of for= 0.6. The dashed line is theArrhenius T-dependence tted to the data above T m, whereas the solid line is theVogelFulcher behaviour tted to all of the data. c, T-dependence of for= 0.6.Owing to strong anisotropic interactions, the relaxation function is not a singleexponential ( 0.95) even at high temperatures. Such behaviour is known forhydrogen-bonding liquids such as glycerol.locally favoured structure is consistent with the prediction of ourtwo-order-parameter model of liquid (see equation (2) in ref. 20)that assumes that a liquid consists of only two types of structure;namely, locally favoured structures (blue particles in Fig. 2) andnormal-liquid structures (green particles in Fig. 2).We also calculated the rotational autocorrelation functionCR(t) = (1/N )i ui(t) ui(0), where t is time and N is thenumber of particles. Then we tted the stretched exponentialfunction to CR(t) as CR(t) exp((t/)). Figure 3b and cshows the T-dependence of the characteristic rotational relaxationtime and the Kohlausch exponent obtained by the tting,respectively, for = 0.6. can be well-tted by the VogelFulcher law = 0 exp(DT0/(T T0)) with 0 = 0.61, D = 7.4and the VogelFulcher temperature T0 = 0.099. The increase of on cooling is often classied between strong and fragile extremesusing Tg as a scaling parameter3,5. Liquids whose obeys theArrhenius law are called strong, whereas fragile liquids havesuper-Arrhenius behaviour. D is an indicator of the strength ofa liquid: the larger D is, the stronger the liquid is. Here D = 7.4corresponds to a fragile liquid. The T-dependence of starts todeviate from the Arrhenius behaviour below Tm (=0.46), which wedene as Tm of a crystal free from frustration (that is, at = 0)17. is almost constant above Tm, whereas it monotonically decreaseswith a decrease in T below Tm. Both results indicate that a liquidshows complex cooperative behaviour only below Tm. We alsocalculated the mean square displacement r2(t) (detailed later)and the intermediate scattering function and found that there is atwo-step relaxation (fast and relaxation) at temperatures lowerthan Tm. We conrmed dynamic heterogeneity of a supercooledliquid by means of a time-dependent four-point density correlationfunction26 (see the Supplementary Information) and found thatdynamic heterogeneity becomes more pronounced with a decreasein T (we provide further details later). Thus, we may say that thissystem reproduces the most essential features of the liquidglasstransition well, which are known from previous studies16.Here we focus on the relationship between the dynamics andthe structure; specically, the medium-range crystalline order.It was pointed out16,21 that revealing the relationship betweendynamics and structure is key to the understanding of the slowdynamics associated with the liquidglass transition. This problemis also related to the problem of the complex energy landscape5,2224.The most popular method is to focus on an inherent structure25or an ideal glass structure, which is believed to form local minima(basins) in the energy landscape. In contrast to this popular model,hereafter we propose a possibility that a crystalline structure is a keystructure of a supercooled state.First, we show evidence for the existence of long-lived clustersof medium-range crystalline order in a supercooled liquid state.In our system, we know the equilibrium crystalline structures (seeFig. 2). Note that only one type of crystal, which has AFM spinorder, exists in the glass-forming region ( 0.6). The orderparameter i characterizing it can then be dened as i(t) =(1/6)i,j |ui(t) uj(t)|, where i, j represents the summationover the nearest neighbours of particle i. When i 0.8, we judgethat particle i is involved in clusters of high crystalline order. Wejudge that particle j is also involved in the same cluster if the spin ofparticle j is parallel or antiparallel to that of particle i; specically,if |ui uj| 0.9. To avoid the eects of thermal noise and detectregions with slow dynamics, we averaged the order parameter overthe duration of asi(t) = (1/) t +/2t /2dt i(t ).We emphasize that this orientational order parameter is muchmore sensitive to the crystalline order than the positional orderparameter, which is one of the merits of our model system.In Fig. 4ac, we present snapshots of the spatial distributionof the orientational order parameter for three temperatures.At T = 0.35, there are few clusters of high crystalline order.With a decrease in T , clusters of high crystalline order appearand grow in size and lifetime. We also present rotationaldynamics in Fig. 4df and translational dynamics (particletrajectories) in Fig. 4gi corresponding to Fig. 4ac. Despitethe fact that Fig. 4ac are snapshots, Fig. 4df show the time-averaged information and Fig. 4gi show trajectories, theyare well-correlated: particles are less mobile both rotationallyand translationally in ordered regions than in disorderedregions. To demonstrate this correlation more quantitatively, wealso calculated the rotational correlation function CR( , t)and the mean square displacement r2( , t) for regionshaving a certain range of the orientational order parameter(i ) at t = 0. Here CR( , t) = (1/N )iui(t) ui(0), where N is the number of particles belonging toi and r2( , t) = (1/N )i [ri(t) ri(0)]2. As shownin Fig. 5a, we reveal that the higher the degree of crystalline order,the slower the local rotational relaxation. The correlation betweentranslational motion and the degree of the crystalline order can also202 nature physics VOL 2 MARCH 2006 www.nature.com/naturephysicsUntitled-3 3 2/18/06, 11:26:29 AMNature Publishing Group 2006ARTICLESa b cd e fg h iFigure 4 Medium-range crystalline order and dynamic heterogeneity. ac, Snapshots of liquid structures (at a certain time ts) for three different temperatures:a, T = 0.35, b, T = 0.18 and c, T = 0.17 (see also Supplementary Information, Video S1, which corresponds to b). Red particles havei greater than 0.75. Blue particles,on the other hand, represent locally favoured structures with ve-fold symmetry. The arrows indicate the spin direction or the direction of ui. df, The slowness of therotational dynamics. The red and blue particles denote those slowly rotating during ts < t < ts + /2: the particles belonging to the locally favoured structures are colouredblue. gi, The trajectories of corresponding to ac. The trajectories in gi are drawn from ts to ts + , ts +2 and ts +4 , respectively. Note that ac are snapshots at ts.be seen in Fig. 5b. In Fig. 5c and d we also show the distributions ofCR( , t1) (t1 = 16,400) and r2( , t2) (t2 = 7,900), respectively,where t1 roughly corresponds to the rotational relaxation time, (more precisely, t1 = 0.87), and t2 is the time when theheterogeneity of the translational motion, which is characterizedby 4 (see the Supplementary Information), becomes maximum.These results clearly indicate a positive correlation between thedegree of crystalline order and the slowness of the structuralrelaxation, that is, a link between crystalline order and dynamicheterogeneity.The above results show that clusters of high crystalline orderemerge below Tm and their average size and lifetime increaseon cooling. The lifetime is longer than (see SupplementaryInformation, Fig. S2) and reaches the order of 10 near Tg. Thisis consistent with experimental ndings2. Particles do not rotate fora long time in these clusters, as can be seen from SupplementaryInformation, Video S1. As this heterogeneity has a nite lifetime,it can be regarded as dynamic heterogeneity. The clusters aretemporally created and annihilated, and their average size remainsconstant with time. Thus, we speculate that the size of the clustersis still smaller than the size of the critical crystal nucleus, but thisneeds to be checked carefully in the future.Figure 6a shows the T-dependence of the characteristic sizeof clusters, . We dene this size (T) as (T) =N , whereN is the average number of particles contained in a clusterwith 0.75. We also plot the dynamic coherence lengthobtained by the analysis using the four-point time-correlationfunction26, 4 (see the Supplementary Information), in Fig. 6a.Note that 4 characterizes the heterogeneity of translational motion(see also Fig. 4gi). We found that 4 is proportional to ,and both and 4 can be tted by a function of (4) (T) =A(4) (T T0)2/d (d is the spatial dimensionality; d = 2 in ourcase) with T0 = 0.099, which is obtained by the tting of theVogelFulcher relation to the T-dependence of (see Fig. 3b).The behaviour of and 4 above is consistent with the scalingargument8,9, which connects the increase of and 4 to that of .nature physics VOL 2 MARCH 2006 www.nature.com/naturephysics 203Untitled-3 4 2/18/06, 11:26:31 AMNature Publishing Group 2006ARTICLESLog t Log tSlowSlowFrequencyFrequency1.010110010 110 20.80.60.40.20.0 1 1 10 0 0.2 0.4 0.6 0.8 10 2 3 4 5 61 1 0 2 3 4 5 61024681002468Total= 0.65= 0.70= 0.75= 0.80= 0.85Total= 0.40 (LFS)= 0.65= 0.70= 0.75 = 0.85 = 0.80Total= 0.40 (LFS)= 0.65 = 0.70= 0.75= 0.85 = 0.80Total= 0.65= 0.70= 0.75= 0.80= 0.85CR(,t)CR( ,t1)r2 (,t)r 2( ,t2)a bc dFigure 5 Relation of rotational and translational dynamics to the degree of local crystalline order. a, CR( , t ) for T = 0.17 and= 0.6. Note that the dynamicheterogeneity smears out for large t. b, r2( , t ) for T = 0.17 and= 0.6. c, The distribution of CR( , t1) (t1 = 16,400). Here values are 0.02. The particles withlarger are slower in rotation, but the particles belonging to the locally favoured structures (LFS, 0.4) are also slow. d, The distribution of r2( , t2) (t2 = 7,900). Theparticles with larger are slower in translation. Unlike rotational motion, the particles belonging to the locally favoured structures are rather fast in translational motion. Notethat these distribution functions are normalized such that the integral becomes one.In addition to the growth of clusters, the increase in the fractionof ordered clusters and the resulting percolation and/or theincrease in the frictional contacts between clusters may also causeslow dynamics. Although we cannot specify which is the mostrelevant mechanism for slow dynamics, we may say that dynamicheterogeneity and the slow dynamics are intrinsically related tomedium-range crystalline ordering in a supercooled liquid. Nextwe focus on another important structure, that is, locally favouredstructure, which induces frustration against crystallization. Wefound that the rotational dynamics of particles in locally favouredstructures is slow (see Figs 4df and 5c), although their translationaldynamics are not so slow (see Fig. 5d). Thus, both medium-rangecrystalline ordering and short-range bond ordering contributeto slow rotational dynamics near Tg. However, as the spatialdistribution of locally favoured structures is rather random, theirconnection to dynamic heterogeneity (or the growing length scale)remains the future problem.We note that medium-range crystalline ordering can also beconrmed in the radial distribution function g(r) shown in Fig. 3a.We can clearly see that the height of the peak correspondingto the crystalline order increases with decreasing T . We do notnd any signicant change in the structure factor S(q) (q is thewavenumber) in the low-q region (q 1/). This is consistent withthe well-known fact that dynamic heterogeneity cannot be detectedin small-angle scattering experiments. The reason may be due tothe fact that clusters with high crystalline order appear randomly inspace without spatial correlation. Our simulation results indicatethat there is a possibility to detect dynamic heterogeneity, ormedium-range ordering, in g(r) and S(q) not in a low-q regionwhere q 1/, but rather in a high-q region where q 1/ ( isthe particle diameter) through medium-range crystalline ordering.In relation to this, it is worth mentioning the experimentalstudy27 on the T-dependence of S(q) in propylene carbonateand phenyl salicilate (salol). There it was found that shouldersin S(q), which were identied as peaks from clusters of highcrystalline order, progressively develop in the supercooled liquidsbelow Tm on cooling. Our simulation results are consistent withthis experimental nding, which may imply that this behaviour isgeneric to any glass-forming liquids. Further experimental studiesalong this line are highly desirable.Our study indicates that glassy slow dynamics can beinduced by frustration between long-range density orderingtowards crystallization and short-range bond-ordering towards theformation of locally favoured structures. According to this1719, aliquid suering from stronger frustration should be stronger. In oursimulations, we found that D increases with an increase in , asshown in Fig. 6b: D = 7.4, 10.9, 17.2, 29.7 and 84.3 for = 0.6,0.65, 0.7, 0.75 and 0.8, respectively. We stress that this range of D204 nature physics VOL 2 MARCH 2006 www.nature.com/naturephysicsUntitled-3 5 2/18/06, 11:26:34 AMNature Publishing Group 2006ARTICLES864265432100.2 0.3 0.4 0.2 0.4 0.6 0.8 1.0Log ( T)T Tg/T4= 0.60= 0.65= 0.70= 0.75= 0.80a bFigure 6 Growth of dynamic heterogeneity on cooling and the relationshipbetween and the fragility. a, T-dependence of the characteristic size ofcrystal-like clusters and the correlation length of translational motion 4 obtainedfrom the analysis of the four-point time-correlation function (see the SupplementaryInformation) for= 0.6. The solid and dashed lines are the ttings of A/(T T0)with T0 = 0.099 for and 4, respectively. We found that 4. The error barindicate the range of observed during 10 . b, The Angell plot of against Tg/T.Tg = 0.150, 0.158, 0.163, 0.167 and 0.166 for= 0.6, 0.65, 0.7, 0.75 and0.8, respectively.almost covers from the fragile (for example, o-terphenyl, D 5) tothe strong extremes (for example, silica, D 60). This is anotherindication that the frustration against crystallization controls thenature of the liquidglass transition including the fragility of liquid.Network-forming tendency seems not to be a prerequisite formaking liquids strong, in contrast to common beliefs. Our studysuggests that the model where the slow dynamics are caused byjamming owing to dense packing may not be enough to understandthe glass transition. It may be important to consider why a systemcan remain in a disordered state and be jammed without ordering;namely, frustration adversely aects crystallization.METHODSMODELHere we use molecular dynamics simulation to investigate the roles ofcrystallization in glass transition in the spirit of our two-order-parametermodel of the liquidglass transition1719. For this purpose, we introduce apotential that can directly and systematically control the strength of frustrationagainst crystallization.U (rij ,ij) = U (rij)+U (rij ,ij),where rij is the distance between two particles i and j, andij expresses theorientation between these particles. The rst term of the potential is anisotropic potential that tends to maximize the density of a liquid and encouragecrystallization. The second term is an anisotropic potential that leads to theformation of locally favoured structures. In this paper, we adopted theLennardJones potential as the isotropic potential U (rij):U (rij) = 4((rij)12(rij)6),where is the well depth of the potential and is the diameter of the particle.The anisotropic part of the potential is given byU (rij ,ij) =4(rij)6[h(i 0c)+h(j 0c) 6435c],h(x) = 13x2 +3x4 x6 (1 < x < 1); h(x) = 0 (x 1, x 1).i is an angle between the relative vector rji = rj ri and the unit vector ui ,which represents the orientation of the axis of particle i. j is an angle betweenthe relative vector rij = ri rj and the unit vector uj (see SupplementaryInformation, Fig. S1). The function h(( 0)/c ) (0 = 126 and c = 53.1,where c is the width of the angular part of the potential) has a maximum at = 126, and 0 is an energetically favoured value of i , and is the angle whereh(( 0)/c ) becomes maximum. Thus, this term stabilizes the locallyfavoured structure of ve-fold symmetry, as shown in SupplementaryInformation, Fig. S1. controls the stability of the locally favoured structure as well as thestrength of frustration against crystallization. Our system has frustration that isrelated not only to the orientation of the particle axis, which may be regarded asspin, but also to the position of the particle. In this sense, our model can beconsidered as a lattice-free spin system with translational degrees of freedom(lattice-free spin glass).SIMULATION METHODThe number of particles in our simulations is N = 1,024. We used periodicboundary conditions. We cut o and smoothed both the potential and itsderivative at 3.0 to save computational time. Here we show all results inreduced units. We used particle mass m, and as the basic unit of mass,length and energy, respectively. Thus, the moment of inertia I , the temperatureT , the pressure P, the distance r and the time t are scaled as I I/m2,T kBT/ (kB is Boltzmanns constant), P P2/, r r/ and t t/( =m2/).To control T and P of a system, we used the NosePoincareAndersen(NPA) thermostat28 and extended the NPA method to include the rotationaldegrees of freedom of particles. Our modied NPA hamiltonian (HNPA) isgiven byHNPA(t) = (HNA(t)HNA(0))s,HNA = p22mVs2+ p22I s2+U (qV 1/2,{i})+ gkBT lns+ 2V2MV+ 2s2Ms+PV ,where NA is the NoseAnderson hamiltonian, p is the momentum, p is theangular momentum, i is the angle between x axis and ui , g is the total numberof the degrees of freedom of this system, s is the variable to controltemperature, V is the two-dimensional volume of the system,V is themomentum conjugate to V ,s is the momentum conjugate to s, MV is themass of V , Ms is the mass of s and q is a scaled position vector. Here we usedI = 0.04, MV = 0.0001 and Ms = 1.0. We conrmed the equipartition of theenergy into the translational and the rotational degrees of freedom of particles.According to the theoretical estimate, the variance of the internal temperaturedistribution should be T2/N and 2T2/N for translation and rotation,respectively, which is also conrmed in our simulations within an accuracy of2%. We also proved that the above NPA method extended to a system with therotational degrees of freedom constitutes the NPT ensemble. The resultsshown in this article are for P = 0.5 and a cooling rate of Q dT/dt = 104. Asystem with the LennardJones potential U (r) melts at Tm = 0.46 at P = 0.5.After quenching the system into a target temperature, we annealed it for a longtime (during 30200 ) before making any analysis.Received 21 September 2005; accepted 18 January 2006; published 19 February 2006.References1. Sillescu, H. Heterogeneity at the glass transition: a review. J. Non-Cryst. Solids 243, 81108 (1999).2. Ediger, M. D. Spatially heterogeneous dynamics in supercooled liquids. Annu. Rev. Phys. Chem. 51,99128 (2000).3. Angell, C. A., Ngai, K., McKenna, G. B., McMillan, P. F. & Martin, S. W. Relaxation in glassformingliquids and amorphous solids. J. Appl. Phys. 88, 31133157 (2000).4. Richert, R. Heterogeneous dynamics in liquids: uctuations in space and time. J. Phys. Condens.Matter 14, R703R738 (2002).5. Debenedetti, P. G. & Stillinger, F. H. Supercooled liquids and the glass transition. Nature 410,259267 (2001).6. Andersen, H. C. Molecular dynamics studies of heterogeneous dynamics and dynamic crossover insupercooled atomic liquids. Proc. Natl Acad. Sci. USA 102, 66866691 (2005).7. Adam, G. & Gibbs, J. H. On the temperature dependence of cooperative relaxation properties inglass-forming liquids. J. Chem. Phys. 43, 139146 (1965).8. Kirkpatrick, T. R., Thirumalai, D. & Wolynes, P. G. Scaling concepts for the dynamics of viscousliquids near an ideal glassy state. Phys. Rev. A 40, 10451054 (1989).9. Xia, X. & Wolynes, P. G. Microscopic theory of heterogeneity and nonexponential relaxations insupercooled liquids. Phys. Rev. Lett. 86, 55265529 (2001).nature physics VOL 2 MARCH 2006 www.nature.com/naturephysics 205Untitled-3 6 2/18/06, 11:26:36 AMNature Publishing Group 2006ARTICLES10. Toninelli, C., Wyart, M., Berthier, L., Biroli, G. & Bouchaud, J.-P. Dynamical susceptibility of glassformers: Contrasting the predictions of theoretical scenarios. Phys. Rev. E 71, 041505 (2005).11. Garrahan, J. P. & Chandler, D. Coarse-grained microscopic model of glass formers. Proc. Natl Acad.Sci. USA 100, 97109714 (2003).12. Tarjus, G., Kivelson, S. A., Nussinov, Z. & Viot, P. The frustration-based approach of supercooledliquids and the glass transition: a review and critical assessment. J. Phys. Condens. Matter 17,R1143R1182 (2005).13. Frank, F. C. Supercooling of liquids. Proc. R. Soc. Lond. A 215, 4346 (1952).14. Steinhardt, P. J., Nelson, D. R. & Ronchetti, M. Bond-orientational order in liquids and glasses. Phys.Rev. B 28, 784805 (1983).15. Dzugutov, M. Glass formation in a simple monatomic liquid with icosahedral inherent local order.Phys. Rev. A 46, R2984R2987 (1992).16. Doye, J. P., Wales, D. J., Zetterling, F. H. & Dzugutov, M. The favored cluster structures of model glassformers. J. Chem. Phys. 118, 27922799 (2003).17. Tanaka, H. Two-order-parameter description of liquids: I. A general model of glass transitioncovering its strong to fragile limit. J. Chem. Phys. 111, 31633174 (1999).18. Tanaka, H. Roles of local icosahedral chemical ordering in glass and quasicrystal formation inmetallic glass formers. J. Phys. Condens. Matter 15, L491L498 (2003).19. Tanaka, H. Two-order-parameter model of the liquid-glass transition: I. Relation between glasstransition and crystallization. J. Non-Cryst. Solids 351, 33713384 (2005).20. Tanaka, H. Simple physical model of liquid water. J. Chem. Phys. 112, 799809 (2000).21. Widmer-Cooper, A., Harrowell, P. & Fynewever, H. How reproducible are dynamic heterogeneities ina supercooled liquid? Phys. Rev. Lett. 93, 135701 (2004).22. Goldstein, M. Viscous liquids and the glass transition: A potential energy barrier picture. J. Chem.Phys. 51, 37283739 (1969).23. Sciortino, F. Potential energy landscape description of supercooled liquids and glasses. J. Stat. Mech.P05015 (2005).24. Sastry, S., Debenedetti, P. G. & Stillinger, F. H. Signatures of distinct dynamical regimes in the energylandscape of a glass-forming liquid. Nature 393, 554557 (1998).25. Stillinger, F. H. & Weber, T. A. Hidden structure in liquids. Phys. Rev. A 25, 978989 (1982).26. Lacevic, N., Starr, F. W., Schrder, T. B. & Glotzer, S. C. Spatially heterogeneous dynamics investigatedvia a time-dependent four-point density correlation function. J. Chem. Phys. 119, 73727387 (2003).27. Eckstein, E. et al. X-ray scattering study and molecular simulation of glass forming liquids: Propylenecarbonate and salol. J. Chem. Phys. 113, 47514762 (2000).28. Sturgeon, J. B. & Laird, B. B. Symplectic algorithm for constant-pressure molecular dynamics using aNosePoincare thermostat. J. Chem. Phys. 112, 34743482 (2000).AcknowledgementsThe authors are grateful to C. P. Royall for a critical reading of our manuscript. This work was partiallysupported by a grand-in-aid from the Ministry of Education, Culture, Sports, Science andTechnology, Japan.Correspondence and requests for materials should be addressed to H.T.Supplementary Information accompanies this paper on www.nature.com/naturephysics.Competing nancial interestsThe authors declare that they have no competing nancial interests.Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/206 nature physics VOL 2 MARCH 2006 www.nature.com/naturephysicsUntitled-3 7 2/18/06, 11:26:38 AMNature Publishing Group 2006 /ColorImageDict > /JPEG2000ColorACSImageDict > /JPEG2000ColorImageDict > /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 450 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.00000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict > /GrayImageDict > /JPEG2000GrayACSImageDict > /JPEG2000GrayImageDict > /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 2400 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.00000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict > /AllowPSXObjects false /CheckCompliance [ /PDFX1a:2001 ] /PDFX1aCheck true /PDFX3Check false /PDFXCompliantPDFOnly true /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox false /PDFXBleedBoxToTrimBoxOffset [ 0.30000 0.30000 0.30000 0.30000 ] /PDFXOutputIntentProfile (OFCOM_PO_P1_F60) /PDFXOutputConditionIdentifier () /PDFXOutputCondition (OFCOM_PO_P1_F60) /PDFXRegistryName (http://www.color.org) /PDFXTrapped /False /SyntheticBoldness 1.000000 /Description >>> setdistillerparams> setpagedevice