On the nature of the glass transition in atomistic models of glass formers
OOn the nature of the glass transition in atomistic models of glass formers
Alexander Hudson ∗ and Kranthi K. Mandadapu ∗ Department of Chemical and Biomolecular Engineering, University of California, Berkeley, andChemical Sciences Division, Lawrence Berkeley National Laboratory
We study the nature of the glass transition by cooling model atomistic glass formers at con-stant rate from a temperature above the onset of glassy dynamics to T = 0. Motivated by theEast model, a kinetically constrained lattice model with hierarchical relaxation, we make severalpredictions about the behavior of the supercooled liquid as it passes through the glass transition.We then compare those predictions to the results of our atomistic simulations. Consistent withour predictions, our results show that the relaxation time τ of the material undergoes a crossoverfrom super-Arrhenius to Arrhenius behavior at a cooling-rate-dependent glass transition tempera-ture T g . The slope of ln τ with respect to inverse temperature exhibits a peak near T g that growsmore pronounced with slower cooling, matching our expectations qualitatively. Additionally, thelimiting value of this slope at low temperature shows remarkable quantitative agreement with ourpredictions. Our results also show that the rate of short-time particle displacements deviates fromthe equilibrium linear scaling around T g , asymptotically approaching a different linear scaling. Toour surprise, these short-time displacements, the dynamic indicators of the underlying excitationsresponsible for structural relaxation, show no spatial correlations beyond a few particle diameters,both above and below T g . This final result is contrary to our expectation, based on previous resultsfor East model glasses formed by cooling, that inter-excitation correlations should emerge as theliquid vitrifies. I. INTRODUCTION
A glass may be understood as a liquid whose struc-tural relaxation time, τ , the timescale over whichits constituent particles can rearrange, exceeds theaccessible timescale of observation. On timescalesshort compared to the structural relaxation time,i.e., t (cid:28) τ , the liquid does not flow, instead demon-strating the rigidity of a solid despite lacking obvi-ous structural ordering. Vitrification, the process bywhich a flowing liquid transforms into a rigid glass,can be attributed to the explosive growth in the re-laxation time that accompanies decreasing temper-ature in supercooled liquids. Specifically, it is wellknown that the relaxation time of a liquid exhibitssuper-Arrhenius growth at temperatures below amaterial-dependent onset temperature T o [1, 2]. For T < T o , τ ( T ) grows super-exponentially with theinverse temperature. For historical overviews of thefield, covering both theory and experiment, see ref-erences [3–7]. ∗ We are deeply indebted to David Chandler, without whomthis work would not have been possible. David was inti-mately involved with this work from the beginning and weregret that it could not be completed before he passed away.We wish to emphasize that the opinions and perspective ex-pressed in this paper belong solely to the listed authors, AHand KKM. Although David’s intellectual contributions tothis work warrant his inclusion as an author on this paper,we felt it inappropriate to add his name without his permis-sion. If not for David’s untimely passing, he too would havecontributed to the writing of the manuscript and likely sug-gested significant revisions to meet his exacting standards.
Although there is disagreement about the preciseform of τ ( T ), in this paper we take the perspectiveof the East model [8], a kinetically constrained lat-tice model of glass formers, for which the relaxationtime has been shown theoretically and numericallyto obey parabolic scaling, [9–13]ln τ ( T ) ∼ J (cid:18) T − T o (cid:19) for T (cid:28) T o , (1)where J is a material property and is related to thefree-energetic cost of creating an excitation — a lo-calized region of space where particles can undergonontrivial short-time displacements — in the super-cooled liquid [14]. Equation 1 is noteworthy in thatit remains finite at all nonzero temperatures, in con-trast to other functional forms proposed for τ ( T ),such as the Vogel-Fulcher-Tammann equation com-mon in the literature [3], that diverge at finite tem-perature.The origin of the parabolic scaling in equation 1 isa relaxation barrier that grows linearly with inversetemperature, E ( T ) ∝ /T . In the East model, thislinear scaling is a consequence of two factors. Thefirst is hierarchical dynamics, which dictates that thefree-energetic barrier to relaxing an inter-excitationdomain of length (cid:96) grows logarithmically with (cid:96) , i.e., E (cid:96) = γJ ln (cid:96), (2)where J is the energy associated with an excita-tion and γ measures the entropy of relaxation path-ways [9, 10, 13]. The second factor is the growinglengthscale between excitations as temperature is a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r lowered. Because excitations cost energy, their con-centration will decrease with temperature accordingto c ( T ) ≈ exp ( − J /T ), so that the mean inter-excitation spacing (cid:96) = 1 /c will obey (cid:96) ( T ) ≈ e J /T . (3)Combining equations 2 and 3, we find that E ( T ) ≈ γJ /T . The relaxation time is then given by τ ( T ) = τ e E ( T ) /T ≈ τ e γJ /T . (4)This form for the relaxation time of the East modelwas first derived by Sollich and Evans in reference [9]and later confirmed analytically [12, 13]. Comparingequation 4 to equation 1, we see that the curvatureparameter in the parabolic law is related to the ex-citation energy scale via J = γJ .Regardless of the exact form of τ ( T ), the super-Arrhenius growth of the relaxation time ensures thata modest lowering of temperature can transform anergodic liquid into a non-ergodic glass, the latter un-able to relax on the finite timescale of interest. Al-though it is possible, in principle, to prepare an equi-librium liquid at any temperature T >
0, in practicethe finite timescale implicit in any preparation pro-cess, whether performed in the laboratory or on acomputer, will be dwarfed by the relaxation time ofthe liquid at sufficiently low temperatures.To make the preceding discussion concrete, welimit our consideration to protocols in which a bulkliquid is in contact with a heat bath whose temper-ature follows a schedule T ( t ). Although T ( t ) couldbe a complicated function, in this paper we focusexclusively on constant-rate cooling protocols, suchthat T ( t ) decreases linearly from an initial temper-ature T i to a final temperature T f in a time t c . As-suming T i is a high temperature, say T i (cid:38) T o , theprotocol is entirely characterized by just two param-eters: the final temperature T f and the cooling rate ν c ≡ | dT ( t ) /dt | = ( T i − T f ) /t c .The inverse cooling rate, ν − , is the timescale overwhich the protocol produces significant changes intemperature. If this timescale is long compared tochanges in the relaxation time of the liquid as tem-perature is lowered, then the liquid will remain inequilibrium. Mathematically, we can express thiscondition as ν − (cid:29) | dτ /dT | . If the opposite is true, ν − (cid:28) | dτ /dT | , then the timescale set by the pro-tocol is short compared to changes in the relaxationtime and the liquid will fall out of equilibrium.If we choose the final temperature of the protocolto be T f = 0, then there will be a crossover at an in-termediate temperature T g where the two timescalesare equal, ν − = (cid:12)(cid:12)(cid:12)(cid:12) dτ ( T g ) dT (cid:12)(cid:12)(cid:12)(cid:12) . (5) For T > T g , the liquid remains in equilibrium; for T < T g , the liquid is out of equilibrium. For T (cid:28) T g ,the liquid has fully transformed into a glass. Thus,equation 5 formally defines the glass transition tem-perature T g [15, 16]. Although other definitions arepossible, equation 5 clearly expresses the concept ofthe glass transition as a dynamic process in whichan equilibrium supercooled liquid gradually trans-forms into a non-ergodic material. Equation 5 alsomakes explicit the protocol dependence of the glasstransition. Note that if τ ( T ) is non-singular for all T >
0, as it is for the parabolic law, then in the limitof infinitely slow cooling, ν c →
0, equation 5 impliesthere is no finite-temperature glass transition.In this paper, we subject two atomistic models ofglass formers commonly used in simulation studies ofglassy dynamics, the Wahnstr¨om and Kob-Andersenmodels [17, 18], to constant-rate cooling as describedabove. Using equation 5 and the East model rela-tions, equations 1–4, we make predictions about theglass transition and compare those predictions to theresults of our atomistic cooling protocols.Consistent with our prior expectations, we find acrossover in the relaxation time from the parabolicscaling of equation 1 to linear (Arrhenius) scaling.This crossover is characterized by a peak in the slopeof ln τ ( T ) versus 1 /T , which grows more pronouncedwith slower cooling. This behavior agrees qualita-tively with our predictions, which are based on ap-proximate but intuitive arguments. Specifically, theslope exhibits its peak near T g given by equation 5.Furthermore, our prediction for the slope at T (cid:28) T g ,which we obtain by evaluating the energy barrier inequation 2 at T g , agrees quantitatively with the sim-ulation results for the Kob-Andersen model.In addition to the super-Arrhenius-to-Arrheniuscrossover, we also observe a crossover in the rate atwhich particles displace over short times t (cid:28) τ , fromthe equilibrium linear scaling established in reference[14] to a different linear scaling below T g . Becausethe rate of short-time displacements is closely re-lated to the concentration of excitations in the Eastmodel, we interpret this result as indicating that thenumber of excitations in the liquid approaches a con-stant, limiting value for T (cid:28) T g , with the residualslope corresponding to a temperature-independentfree-energetic barrier to short-time dynamics withinan excitation.To our surprise, the spatial distribution of short-time particle displacements remains ideal-gas-like inthe atomistic models, even for T < T g . To the ex-tent that these short-time displacements are faith-ful indicators of underlying excitations (which areotherwise undetectable), this contradicts East modelresults showing pronounced inter-excitation correla-tions emerging below T g [14, 19, 20]. These correla-tions are a natural consequence of the East model’shierarchical dynamics, and their absence in the spa-tial distribution of particle displacements below T g is puzzling. II. MODELS AND METHODSA. Models
We consider binary mixtures of Lennard-Jonesparticles, which are governed by pairwise potentialsof the form u αβ ( r ) = − (cid:15) αβ (cid:20)(cid:16) σ αβ r (cid:17) − (cid:16) σ αβ r (cid:17) (cid:21) , (6)where subscripts α, β denote the particle type, A or B . The models chosen for this study are character-ized by the Lennard-Jones parameters (cid:15) αβ and σ αβ ,the dimension d , and the mole fraction f A of A -typeparticles. Additionally, for each model we fix thetotal density, ρ ≡ N/V , choosing values common inthe literature. With the model and density specified,the behavior of the supercooled liquid — character-ized by parameters such as the onset temperature T o and energy scale J — is determined [1, 2, 14].Binary mixtures are common choices in studies ofglass-forming liquids due to their resistance to crys-tallization, which enables study of the liquid into thesupercooled regime. We consider two binary mixturemodels in this paper: the Wahnstr¨om model, a 50-50 mixture, and the widely studied Kob-Andersenmodel, an 80-20 mixture with non-additive parame-ters [17, 18]. Table I lists the parameters for thesemodels.At sufficiently low temperatures, both models willsuccumb to crystallization. In our equilibrium simu-lations, we observe the Wahnstr¨om model crystalliz-ing at temperatures satisfying T − ≡ β (cid:38) . β (cid:38) .
3. Given our estimates of T o for these twomodels, the temperature ranges over which we cansafely study the supercooled liquid are 1 . (cid:46) β (cid:46) . . (cid:46) β (cid:46) . B. Simulation details
Unless otherwise noted, we express all distancesin units of σ ≡ σ AA and all times in units of τ LJ ≡ (cid:0) mσ /(cid:15) (cid:1) / , where m ≡ m AA . We performed all simulations us-ing the HOOMD-blue molecular dynamics softwarepackage [25, 26], running primarily on GPUs. Wecompiled HOOMD-blue in single precision and typ-ically used a time step of δt = 0 . . C. Calculating non-equilibrium averages
In this paper, we sometimes refer to averages inand out of equilibrium with some degree of ambigu-ity. Here, we clarify the meaning of an average andoffer some comments about how these averages are
Model ρ f A σ AA σ AB σ BB (cid:15) AA (cid:15) AB (cid:15) BB m A m B T o W 1.296 0.5 1 11 /
12 5 / . .
88 1 1.5 0.5 1 1 0.72TABLE I. Parameters of the models considered in this work. The parameters listed between the double columnbars fully define the models. Definitions of these parameters are as follows: ρ is the number density N/V ; f A is thefraction of particles that are type A ; σ ij and (cid:15) ij are the Lennard-Jones parameters for ij pairs; and m i is the mass oftype- i particles. The last parameter in the table, T o , is onset temperature as determined by fitting a parabolic lawto equilibrium relaxation time data. At T < T o , the models exhibit glassy dynamics. computed in practice. To begin, consider a quantity x that is a function of the configuration C of the sys-tem, x = x ( C ). Of interest is the ensemble average of x , (cid:104) x (cid:105) ≡ (cid:88) C w ( C ) x ( C ) , (7)where w ( C ) is a weight function that determinesthe contribution of configuration C to the aver-age. Choosing an ensemble means choosing thefunction w . The micro-canonical ensemble weightsall allowed configurations equally, while the canon-ical ensemble weights configurations by the Boltz-mann factor exp( − E ( C ) /T ), where E ( C ) is the en-ergy of configuration C and T is the temperature[28]. Non-equilibrium ensembles, corresponding tosystems driven out of equilibrium by a protocol oftime-varying external parameters, are also possible,though an analytical expression for w may be knownonly in special cases.In practice, the ensemble average of x is computedby a simple average of x over configurations drawnfrom the corresponding distribution w . If we sample M configurations, then we estimate (cid:104) x (cid:105) with x ≡ M M (cid:88) i =1 x ( C i ) . (8)Provided the set of configurations {C i } is represen-tative of the ensemble, x will be a good estimate ofthe true ensemble average. For equilibrium ensem-bles, this representative subset is typically generatedby periodically sampling configurations from a suffi-ciently long time series following a generous periodof equilibration.Non-equilibrium ensembles require specification ofan initial condition , i.e., the ensemble that prevailsat t = 0, as well as the time dependence of allexternal parameters for t >
0. This time depen-dence means that time averaging is not an option outof equilibrium. Instead, we generate representativesubsets by independently preparing a large numberof systems according to the protocol, drawing initialconfigurations from the known ensemble at t = 0 andthen subjecting those configurations to the protocol of time-varying external parameters. The configura-tions obtained at the end of the protocol comprisea representative subset of the non-equilibrium en-semble. Because the cooling protocols we considerbegin with the system in equilibrium at a high ini-tial temperature T i , the starting configurations aredrawn from the canonical ensemble at temperature T i . All parameters besides the temperature are fixedand T decreases linearly with time, so the rate ν c issufficient to fully specify a cooling protocol.So far we have implicitly limited ourselves to cool-ing protocols that take the temperature all the waydown to zero. However, we could imagine cool-ing to an intermediate temperature T such that T i > T >
0. The rate ν c then corresponds to afamily of protocols and associated non-equilibriumensembles, differentiated by the final temperature T to which the system is cooled. Note that the ensem-ble corresponding to ν c and T could also be producedby pausing the full cooling protocol when it reachestemperature T and sampling the configurations atthat point, which is the approach we take in prac-tice. If x = x ( C ) is our quantity of interest, then wedenote the average over this non-equilibrium ensem-ble by (cid:104) x ( T ) (cid:105) ν c . For a fixed rate ν c , this quantity is afunction of temperature indicating how the averageof x changes as the system is cooled to lower tem-peratures at rate ν c . As the rate changes, so doesthe function.For any ν c >
0, there will be a glass transitiontemperature, T g , satisfying equation 5, below whichthis function will deviate substantially from the equi-librium result. In the limit of infinitely slow cooling, ν c →
0, we expect (cid:104) x ( T ) (cid:105) ν c to approach the equilib-rium ensemble average. In the work that follows, wewill typically drop the subscript ν c for ease of no-tation, but it should be clear from context whetherthe average is equilibrium or non-equilibrium, andfor the latter, what the corresponding rate is. D. Iso-configurational averaging
In the preceding discussion of non-equilibrium av-erages we restricted our attention to static observ-ables that depend on a single configuration. How-ever, in much of the work that follows, we will beconcerned with dynamic observables that depend onthe configuration at multiple points in a trajectory,e.g., x = x ( C ( t ) , C ( t ) , . . . , C ( t N )) , where t < t < · · · < t N . Typically we will focuson observables that depend only on the endpoints ofthe trajectory, so that the above simplifies to x = x ( C (0) , C ( t )) , where t is the total length of the trajectory, or ob-servation time.The so-called iso-configurational ensemble [29] as-sociated with an initial configuration C (0), a tem-perature T , and an observation time t is gener-ated by drawing initial velocities from the Maxwell-Boltzmann distribution at temperature T and thenintegrating Newton’s equations of motion forwardto time t to produce a new configuration C ( t ). Theaverage of a dynamical observable x over this newensemble is termed the iso-configurational average ,and is itself a static observable that depends only onthe initial configuration C (0). For convenience, wedenote this observable by x iso ( C ).We will often write (cid:104) x ( C (0) , C ( t )) (cid:105) when we meanto refer to the ensemble average of x iso , so that thefollowing expressions are equivalent: (cid:104) x ( C (0) , C ( t )) (cid:105) ≡ (cid:104) x iso (cid:105) . That is, when dealing with a dynamical observable,the angle brackets denote an ensemble average aswell as an iso-configurational average. In practice,this average is computed via an extension of equa-tion 8: x = 1 M M (cid:88) i =1 N iso N iso (cid:88) j =1 x ( C i , C ij ( t )) . In this expression, N iso is the number of trajecto-ries run per initial configuration C i to sample its iso-configurational ensemble, and C ij ( t ) is the configu-ration produced at time t in the j th trajectory thatstarts with configuration C i .For x to produce a good estimate of (cid:104) x ( t ) (cid:105) , it isimportant that M is large. By contrast, it is not particularly important that N iso be large, provided M is large. The iso-configurational average is essen-tially just an average over the initial velocities of thetrajectories, and for a sufficiently large value of M ,this velocity averaging is satisfied even for N iso = 1.In practice, however, the value of M will be limited(by computational resources), and a larger value of N iso will sometimes be necessary. III. RESULTS AND DISCUSSION
In this section, we present the results of ouratomistic simulations and compare them, in detail,to expectations based on the one-dimensional Eastmodel. For each model listed in table I, we per-formed cooling protocols as described in Models andMethods, for various cooling rates. For each modeland cooling rate, at various values of the externaltemperature T , we compute the following propertiesof interest:(i) the relaxation time of the liquid;(ii) the rate at which particles make significant dis-placements over short observation times;(iii) the spatial distribution of particles that makethese short-time displacements.Each of these properties obeys East model scalingin equilibrium and is expected to deviate from thatscaling as T drops below T g .In the remainder of this section, we discuss each ofthese properties in greater detail, our expectationsfor them (in and out of equilibrium), and how theresults of our atomistic simulations compare to theseexpectations. A. Relaxation time
Following reference [14], we define the relaxationtime of the liquid as the time τ satisfying F s ( q , τ ) =1 /e , where F s ( k, t ) ≡ (cid:104) exp [ i k · ( r i ( t ) − r i (0))] (cid:105) (9)is the self-correlation function and q is the wavevec-tor that maximizes the static structure factor. Theangle brackets in equation 9 may denote either anequilibrium or non-equilibrium average. In the non-equilibrium case, the average is performed over anensemble of constant-energy trajectories of length t , with initial momenta drawn from the Maxwell-Boltzmann distribution and the initial configurationdrawn from a specified non-equilibrium ensemble.See Models and Methods for more detail.At high temperatures, T > T o , the relaxation timeas defined above is Arrhenius, obeying the linear re-lationship ln τ ( T ) = ln τ + E (cid:18) T (cid:19) , (10)where τ and E are material properties governingthe dynamics of the high-temperature liquid. Theseparameters can easily be determined by fitting equa-tion 10 to relaxation time data computed at high ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● Wahnström modelKob-Andersen model J T T o J T T o l n ( ⌧ / ⌧ h i ) l n ( ⌧ / ⌧ h i ) equilibriumequ fit ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ equilibriumequ fit ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ FIG. 1. Logarithm of the relaxation time (with the high-temperature Arrhenius contribution removed) versus inversetemperature for the Kob-Andersen (left) and Wahnstr¨om (right) models, for various cooling rates. Rates have unitsof temperature/time, where time is expressed in units of τ LJ . Relaxation times are computed as described in themain text. The dashed black curves show the result of fitting equilibrium data (gray points) to the equilibrium formin equation 12. Solid lines are guides to the eye. temperatures. The energy scale E may be inter-preted as an intrinsic energy barrier that particlesmust overcome in order to undergo a local reorgani-zation. We expect this energy barrier to be presentat all temperatures, so that even for T < T o , equa-tion 10 contributes to the relaxation time.As discussed in the Introduction, for T (cid:28) T o ,the relaxation time obeys the parabolic scaling ofequation 1. Combining this parabolic contributionwith the high-temperature contribution in equation10 yieldsln τ ( T ) = ln τ + E (cid:18) T (cid:19) + J (cid:18) T − T o (cid:19) , (11)for T < T o . We can easily combine the high- andlow-temperature forms into a single formula valid atall temperatures:ln τ ( T ) = ln τ o + E (cid:18) T (cid:19) ++ J (cid:18) T − T o (cid:19) Θ (cid:0) T − − T − (cid:1) , (12)where Θ( x ) is the Heaviside step function and en-sures that the parabolic contribution only applies attemperatures below T o .The above equations for the relaxation time arevalid for glass-forming liquids in equilibrium. How-ever, as discussed in the Introduction, a liquid sub-jected to constant-rate cooling will eventually fallout of equilibrium at the glass transition tempera-ture, T g , given by equation 5. At temperatures be-low T g , the structure of the liquid is arrested and energetic barriers to relaxation are fixed. Conse-quently, we expect the relaxation time to exhibit Ar-rhenius behavior for T < T g . This super-Arrhenius-to-Arrhenius crossover is frequently observed in ex-periments and is a key signature of the glass transi-tion [5, 30–34].For the models listed in table I, for each coolingrate considered, we computed the relaxation time τ at various temperatures T . The results are shownin figure 1, where we plot ln ( τ /τ hi ) versus 1 /T foreach combination of model and cooling rate, where τ hi ( T ) ≡ ln τ o + E/T denotes the high-temperaturecontribution in equation 10. At high temperatures, T ∼ T o , where τ changes only modestly with T , therelaxation time shows little dependence on the cool-ing rate. However, the results for each cooling ratebegin to separate noticeably as temperature is low-ered and τ begins to grow more rapidly with decreas-ing T . The effect of the cooling rate on the relaxationtime out of equilibrium is as expected: faster coolingleads τ ( T ) to quickly deviate from the equilibriumscaling, whereas slower cooling leads τ ( T ) to followthe equilibrium scaling to lower temperatures, con-sistent with the lower glass transition temperatureimplied by equation 5. Furthermore, for all coolingrates, for both models, ln τ ( T ) versus 1 /T appears toapproach a straight line, consistent with our expec-tation that the relaxation time exhibits Arrheniusbehavior below T g .The slope of ln τ ( T ) versus 1 /T for T < T g canbe estimated from simple arguments based on theEast model. To start, we write the low-temperaturecontribution to the relaxation time, which we denote τ lo ( T ), as the product of a (temperature-dependent) d l n ⌧ / d /T /T Wahnström modelKob-Andersen model (fast) ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ (medium) (slow) (fast) (slow) FIG. 2. Slope of ln τ ( β ), where β = 1 /T is the inverse temperature, for the Kob-Andersen model (left) and theWahnstr¨om model (right) cooled at varying rates. Open circles show numerical estimates of the slope from the datain figure 1, while dashed lines show predictions from East-like scaling, assuming an abrupt glass transition at T g . free energy barrier and an explicit temperature de-pendence: ln τ lo ( T ) = F ( T ) (cid:18) T − T o (cid:19) . As discussed in the Introduction, super-Arrheniusrelaxation is a consequence of F ( T ) growing with1 /T . Arrhenius behavior below T g , such as that ob-served in figure 1, implies that this barrier has be-come temperature independent, a consequence of thestructural arrest of the liquid as T passes through T g .Assuming the liquid abruptly falls out of equilibriumat T g , the relaxation time will obeyln τ lo ( T ) = F ( T g ) (cid:18) T − T o (cid:19) , for T < T g . That is, the equilibrium barrier at T g is frozen into the material for all T < T g . In theEast model, this barrier is related, through equation2, to the mean inter-excitation spacing, which growsas temperature is lowered according to equation 3.The inter-excitation spacing at the glass transition,denoted (cid:96) ne , is given by (cid:96) ne = e J (1 /T g − /T o ) . This lengthscale is fixed below T g [15, 16], leadingto a fixed energy barrier F ( T g ) = γJ ln (cid:96) ne = γJ (cid:18) T g − T o (cid:19) that prevails below T g . This barrier is also preciselythe out-of-equilibrium slope of ln τ lo ( T ) with respect to 1 /T . Defining J ≡ γJ , the relaxation time foran East-like system with glass transition tempera-ture T g is given byln τ lo ( T ) = J (cid:16) T − T o (cid:17) , T o > T > T g ; J (cid:16) T g − T o (cid:17)(cid:16) T − T o (cid:17) , T > T g . (13)From equation 13, we can see that the derivative ofln τ lo ( T ) with respect to 1 /T drops discontinuouslyby a factor of two, from 2 J (1 /T g − /T o ) immedi-ately before the glass transition to J (1 /T g − /T o )immediately after. This factor-two change in slopeis a corollary of the parabolic law [15].Figure 2 compares the prediction of the precedingarguments with the results of our atomistic simula-tions. The dashed lines show the predicted slope,including the discontinuity at T g . The open cir-cles show numerical estimates of the slope using thedata from figure 1. For the Kob-Andersen model,for the two slower cooling rates (second and thirdpanes), the computational results show fairly goodagreement with the prediction, at least over the lim-ited range of available data. The major differencebetween our predictions and the simulation results,which is clearly visible for the moderate cooling rate(second pane), is the smoothness of the slope around T g . However, this is not unexpected; the disconti-nuity in the predicted slope is a consequence of as-suming that the system falls abruptly out of equi-librium at T g . To the extent that this process isgradual rather than abrupt – perhaps reflecting thenontrivial distribution of length and energy scalesin the system – the slope will exhibit a peak that issmoothed over a finite temperature range around T g .This is indeed what we see in figure 2. Interestingly,for the fastest cooling rate, the Kob-Andersen modelexhibits almost no peak at all in the slope. This maybe a result of cooling too quickly, such that the liq-uid begins falling out of equilibrium before it clearsthe crossover region around T o . Nonetheless, for allthree cooling rates, the limiting slope for T < T g seems to agree fairly well with the predicted slope,at least within the limits of the data we could obtainwith our computational resources.The results for the Wahnstr¨om model are consid-erably worse than for the Kob-Andersen model. Al-though the predicted values for T g seem reasonable,the peaks in the data fall far short of the predictedpeaks and, even worse, there is substantial quan-titative disagreement in the limiting slope. Thesediscrepancies are likely due to the great difficulty inaccurately measuring the parameter J for the Wahn-str¨om model. To obtain a reliable estimate of J requires computing the equilibrium relaxation timeat deeply supercooled conditions, far away from thecrossover to heterogeneous dynamics at T o . Unfor-tunately, the Wahnstr¨om model crystallizes at onlymodest supercooling, so that only a limited range ofdata is available to determine this parameter. TheKob-Andersen model resists crystallization to muchlower temperatures, so that its J value may be ob-tained more reliably. B. Rate of short-time displacements
In supercooled liquids, non-trivial particle dynam-ics is confined to localized excitations or “soft spots”where particles are locally unjammed and can makesmall ( ∼ σ -sized) displacements over short times [14].In analogy with the East model, these localized exci-tations carry a free-energetic cost, which we denote J , and can facilitate motion in neighboring regionsthrough the creation (and destruction) of neighbor-ing excitations.Soft spots in supercooled liquids are structural inorigin, depending only on the set of particle coordi-nates. Nonetheless, the construction of a local, staticorder parameter capable of reliably distinguishingmobile from immobile regions has proven elusive.For this reason, and because non-trivial dynamicsoccurs only near excitations, we identify soft spotsby looking for significant particle displacements overshort observation times. We associate particle i withan excitation if, in a short observation time t , it dis-places a specified minimum distance a ∼ σ . Math-ematically, this is accomplished with the following dynamic order parameter: h i ( t, a ) ≡ Θ ( | r i ( t ) − r i (0) | − a ) . (14) In the above equation, Θ( x ) is the unit step func-tion, and the overlines above the position vectorsindicate inherent structure coordinates. The latteris necessary to remove vibrational motion and can beobtained by steepest descent to the nearest potentialenergy minimum or by a suitable coarse-graining intime [14]. If at time t particle i is at least a distance a away from its initial position, then h i = 1 and weassociate particle i with an excitation.The order parameter h i is dynamic, so that it de-pends on a trajectory rather than a single config-uration. Thus, it depends on both the initial posi-tions and the initial momenta. To construct an orderparameter suitable for detecting excitations (whichdepend on coordinates and not momenta) requiresperforming an iso-configurational average over theinitial momenta, as discussed in the Methods sec-tion. The resulting order parameter, (cid:104) h i (cid:105) iso , givesthe propensity of particle i to move from the speci-fied initial configuration. Particles with (cid:104) h i (cid:105) iso ≈ c ( T ) = e − J /T (15)until the glass transition, at which point the concen-tration plateaus to a cooling-rate-dependent value.We expect an analogous result to hold for atomisticmodels. To test this expectation, we consider thequantity ρ σ ≡ (cid:104) Θ( | r i ( t ) − r i (0) | − σ ) (cid:105) , (16)which is just the ensemble average of h i with a = σ .Note that while the concentration c is a direct mea-surement of the number of excitations in the Eastmodel, ρ σ is merely a proxy for the number of under-lying excitations in an atomistic system. The quan-tity ρ σ gives the average fraction of particles thatmove in an observation time t . In the East modelpicture, particle motion occurs only near excitations,so that ρ σ should be proportional to the number ofunderlying (but unobservable) excitations present inthe liquid. To make this concrete, we write ρ σ ( t ) = r ( t ) N exc , (17)where N exc is the number of particles in the systemassociated with an underlying excitation and thatcould potentially move within a short time interval.In practice, only a fraction of the N exc particles thatcan move will do so, such that ρ σ is proportional tothe total fraction of particles eligible to move.We write the proportionality constant r ( t ) as afunction of time to emphasize that it tends to grownwith the observation time t , as more eligible parti-cles displace. However, the relationship in equation ●●●●●●●●● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● Wahnström modelKob-Andersen model
10 20 30 40 5020 40 60 80 J T T o J T T o equilibriumequ fit ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ equilibriumequ fit ⌫ c = 8 . ⇥ ⌫ c = 8 . ⇥ l n ⇢ l n ⇢ FIG. 3. Fraction ρ σ of particles that have displaced at least σ in an observation time t obs . For the Kob-Andersenmodel (left), we choose t obs = 300. For the Wahnstr¨om model (right), we choose t obs = 10. The shape of the curveschanges little for T < T g over this range of observation times. The gray dashed lines show the expected equilibriumscaling; solid lines are guides to the eye.
17 is meaningful only as long as t is short enoughthat N exc , which depends on the instantaneous con-figuration of the system, may be approximated asconstant. For longer times, the system configura-tion changes appreciably and ρ σ is no longer relatedto an instantaneous set of “excited” or mobile par-ticles. Although this timescale is unknown, it mustbe at least as long as the mean or typical instantontime ∆ t [14], which is the time required for a parti-cle to commit to a new position once it has begun adisplacement, and it must be shorter than the struc-tural relaxation time τ . In practice, we typicallychoose t ≈ ∆ t . At conditions of deep supercooling,under which the system evolves very slowly and therate of displacements is low even within excitations(due to the presence of intrinsic energy barriers),we sometimes choose t (cid:29) ∆ t in order to generatea larger set of observations. In all cases, however,we are careful to ensure that the chosen observationtime satisfies t (cid:28) τ .With a suitable analog for the concentration nowdefined for atomistic systems, we can investigatewhether the crossover in c ( T ) described above for theEast model has a counterpart in atomistic systemscooled below the glass transition. Figure 3 shows ρ σ over temperatures ranging from T ∼ T o to β (cid:29) β g for the set of models and cooling rates considered inthe previous section. The results look qualitativelysimilar to those in figure 1 for the relaxation time:at high temperatures, the data are identical acrosscooling rates, but at T < T g , different cooling ratesproduce different results.At high temperatures, the liquid remains in equi-librium and ln ρ − σ follows the linear scaling expectedfrom equation 15. As T approaches T g , ln ρ − σ begins to deviate from this equilibrium scaling. Consistentwith the effect of cooling rate on T g , slower cool-ing leads to results that begin their deviation fromequilibrium at lower temperatures, so that the equi-librium scaling prevails over a larger temperaturerange. Additionally, at deep supercooling, β (cid:29) β g ,ln ρ − σ increases as the cooling rate is lowered, imply-ing that slower cooling leads to smaller values of ρ σ and lower density of excitations out of equilibrium.This result is expected from the East model and isalso supported by figures 1 and 2, which show thatslower cooling leads to slower relaxation times andlarger free-energetic barriers to relaxation.Also noteworthy in figure 3 is that ln ρ − σ slowlyapproaches linear scaling again at very low temper-atures, with the asymptotic slope noticeably smallerthan the slope in equilibrium. Interestingly, the out-of-equilibrium slope seems to be independent of thecooling rate; indeed, for β (cid:29) β g , the results for eachmodel in figure 3 seem only to differ by their verti-cal offset. This suggests the existence of an intrinsicenergy barrier that controls the rate of short-timedisplacements within excitations, a contribution tothe overall rate of displacements ρ σ that would bepresent in equilibrium but only be clearly manifestedat T < T g , where the number of underlying excita-tions becomes fixed. Although the East model hasno concept of intra-excitation dynamics, it exhibitsanalogous behavior in the instantaneous rate of spinflips, following different linear scalings above and be-low T g .0 C. Spatial distribution of soft spots
A distinguishing feature of the East model is theemergence of inter-excitation correlations as the lat-tice is cooled through the glass transition. In ref-erence [19], Keys, Garrahan, and Chandler demon-strate the emergence of correlations between excita-tions frozen into the lattice by cooling, and relatethe cooling-rate-dependent lengthscale of these cor-relations to the response of the heat capacity uponsubsequent heating.Because the East model Hamiltonian lacks any ex-plicit inter-excitation coupling, excitations should bespatially uncorrelated in equilibrium. However, theEast model’s nontrivial dynamics can lead to spa-tial correlations emerging when the lattice is forcedout of equilibrium by a time-dependent perturba-tion such as cooling. Such perturbations impose fi-nite timescales on the system, preventing its returnto equilibrium and turning the dynamic correlationsintroduced by the kinetic constraint into static cor-relations between excitations.The nonequilibrium correlation length reported inreference [19] follows directly from the hierarchicaldynamics of the East model, and we expect thata similar correlation length should emerge betweenthe dynamical indicators of excitations in atomisticmodels subjected to cooling.In East model glasses, inter-excitation correlationscan be detected by computing P ( (cid:96) ), the distributionof inter-excitation distances, which has the formalmathematical definition P ( (cid:96) ) ≡ (cid:42) n i n i + (cid:96) (cid:96) − (cid:89) j =1 (1 − n i + j ) (cid:43)(cid:30) (cid:104) n i (cid:105) , (18)where the angle brackets denote, as usual, an ensem-ble average. The distribution P ( (cid:96) ) is the probabilitythat a randomly selected excitation is separated bya distance (cid:96) from the closest neighboring excitation.These inter-excitation domains, which lack excita- tions and are thus inactive, are the physical originof the barriers to relaxation in the East model. Ifexcitations are energetically uncoupled, as they arein the East model, then it can be easily shown that P ( (cid:96) ) is exponential, and the distribution of excita-tions on the lattice is described by a Poisson processin space. In this case, the spatial distribution ofexcitations is ideal-gas-like and the location of eachexcitation unrelated to and uncorrelated with thelocations of all other excitations.The exponential distribution of inter-excitationdomain sizes that holds at equilibrium gives way to anon-exponential distribution as the lattice is cooledbelow T g . This is a non-equilibrium effect that owesto the hierarchical dynamics of the East model: as T decreases, smaller inter-excitation domains – whichrelax more quickly – are annihilated, while larger do-mains – which relax more slowly – persist. Domainslonger than a cooling-rate-dependent lengthscale (cid:96) ne are unable to relax on the finite timescale of the pro-tocol and will be frozen into the material, while do-mains shorter than this lengthscale are eliminated.The exponential distribution that prevails at equilib-rium thus gives way to a non-monotonic distributionthat peaks at (cid:96) ne (cid:54) = 0. Naturally, (cid:96) ne increases as thecooling rate decreases, as shown in reference [19].As before, we consider particle displacements, thedynamical indicators of excitations, when workingwith atomistic systems. To probe the spatial dis-tribution of excitations in atomistic glass formers,we need an analogue of the nearest-neighbor distri-bution function defined in equation 18 for the Eastmodel. For convenience, we define the quantities θ ij ( r ) ≡ Θ ( | r i ( t ) − r j (0) | − r )and δ ij ( r ) ≡ δ ( | r i ( t ) − r j (0) | − r ) , which are indicator functions that register whetherparticles i and j are separated at t = 0 by at least the distance r or exactly the distance r , respectively.With these, we can define p ( (cid:96) ) ≡ (cid:42) h i (cid:89) j (cid:54) = i (cid:16) (1 − θ ij ( (cid:96) ))(1 − h j ) + θ ij ( (cid:96) ) + δ ij ( (cid:96) ) h j (cid:17)(cid:43)(cid:30) (cid:104) h i (cid:105) . (19)Recall from equation 14 that h i indicates whetherparticle i has displaced by a minimum distance a inan observation time t and is thus associated with anunderlying excitation.The physical meaning of the distribution p ( (cid:96) ) de- fined in equation 19 can be understood by consider-ing the set of all particles that displace a minimumdistance a , as measured by inherent structure coor-dinates, over the course of a trajectory (i.e., particlessatisfying h i = 1). The probability that one of these1 ` / □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □ □□□□□□□□□□□□□□□□□□□□□□□□ □□□ ` / ⌦ ` ↵ l n P ( ` / ⌦ ` ↵ ) l n P ( ` ) T = 0 . . . . . . . T = FIG. 4. Distribution of (cid:96) for the Kob-Andersen model. Results are shown for the slowest cooling protocol thatwe performed. The main figure shows P ( (cid:96) ) at four temperatures near or below the glass transition temperature T g = 0 .
412 for this choice of model and cooling rate. For small separations, (cid:96) (cid:46) σ , there is a noticeable enhancementin the probability, which we attribute to the collective nature of displacements in a supercooled liquid. Beyond thissmall-lengthscale enhancement, the distributions are remarkably exponential, as demonstrated by the linearity of thedata. The solid curves show linear best fits to the data, excluding the outlier point at (cid:96) = 0 and any points thatrepresent five or fewer samples. Inset: When the distributions are scaled and shifted by the mean value of (cid:96) , theycollapse onto the exponential form ln P ( x ) = − x (dashed line). For all distributions in this figure, we used theindicator function defined in equation 14 with a displacement cutoff a = 0 . σ and observation time t obs = 10. mobile or “active” particles, selected at random, willbe separated from the nearest neighboring activeparticle by a distance x , satisfying (cid:96) ≤ x ≤ (cid:96) + d(cid:96) , isgiven by p ( (cid:96) ) d(cid:96) .For an ideal gas in d dimensions, it can be shownthat this nearest neighbor distribution function hasthe form p ( (cid:96) ) ∝ (cid:96) d − exp (cid:16) − ( (cid:96)/(cid:96) ) d (cid:17) , where (cid:96) is the average distance between particles inthe gas. Although p ( (cid:96) ) for an ideal gas is exponen-tially distributed only when d = 1, a straightforwardchange of variables from (cid:96) to (cid:96) d recovers an exponen-tial distribution: p (cid:0) (cid:96) d (cid:1) ∝ exp (cid:16) − ( (cid:96)/(cid:96) ) d (cid:17) . (20)If the excitations in a supercooled liquid are spatiallyuncorrelated, like the particles in an ideal gas, thenthey will obey the exponential scaling in equation 20.It follows that particles displacing over short time in-tervals, whose positions are strongly correlated withthe locations of the excitations, should also obey this scaling. If inter-excitation correlations emerge out ofequilibrium, as they do in the East model, then p (cid:0) (cid:96) d (cid:1) should deviate from this scaling.Figure 4 shows p (cid:0) (cid:96) (cid:1) at several temperatures T (cid:46) T g for the Kob-Andersen model subjected to ourslowest cooling protocol. For (cid:96) (cid:38) σ , all of the distri-butions appear remarkably (and unexpectedly) ex-ponential. The only deviation from this exponentialscaling is an enhancement of probability at small (cid:96) ,a feature also observed in references [14, 20]. Thissmall-lengthscale feature can be understood as man-ifesting the finite size and localized nature of an ex-citation, which facilitates dynamics for a collectionof particles in a contiguous region of space no morethan a few particle diameters in extent. The dis-placement of a particle signifies the presence of anexcitation and, by extension, the proximity of otherparticles that may readily displace. As a result, par-ticles that displace over short times tend to be highlycorrelated over short distances.The preceding observations make clear that p (cid:0) (cid:96) (cid:1) combines the spatial distribution of underlying exci-tations with that of displacing particles within exci-2tations. As discussed above, the latter contributionis short-ranged and decays on the lengthscale of anexcitation, leaving only the inter-excitation distri-bution at large lengthscales. The distributions plot-ted in figure 4, which appear to be exponential for (cid:96) > σ , suggest that the spatial distribution of exci-tations is ideal-gas-like at all temperatures and un-der all cooling protocols studied. To confirm theexponential scaling at large values of (cid:96) , we removethe small-lengthscale ( (cid:96) < σ ) portion of the distri-butions shown in figure 4, then shift and normalizewhat remains. The modified distributions, plottedin the inset of figure 4, collapse onto the exponentialform. Results similar to those in figure 4 were ob-tained for the Wahnstr¨om model and are presentedin the Supplemental Material.That excitations would be spatially uncorrelatedat high temperatures, where the liquid is at or nearequilibrium, is expected. That they would remain so at temperatures well below T g is a surprise giventhe hierarchical nature of dynamics in supercooledliquids. As discussed previously for the East model,only excitations within ≈ (cid:96) ne of other excitations canbe eliminated in the finite time afforded by the cool-ing protocol; consequently, the out-of-equilibriumliquid should display negative inter-excitation corre-lations similar to those observed for the East model[19]. The lack of such correlations in our results, atconditions where the liquid is clearly out of equilib-rium, is unexpected. To verify the validity of our re-sults, we performed additional calculations, varyingthe definition of the indicator function, the choicesof distance cutoff and observation time, and even thechoice of spatial distribution function. The results ofmany of these calculations are presented in the Sup-plemental Material. In all cases, at all temperatures,for all models and cooling rates considered, there isno clear evidence of correlations between displacingparticles, except for the short-range correlations dis-cussed above.The most straightforward interpretation of the ap-parent contradiction between the above result andour prior expectation is that short-time particle dis-placements are not reliable indicators of excitationsin atomistic systems, contrary to our stated assump-tion. This could be the case if many of the observeddisplacements comprise relatively unimportant (butstill non-vibrational) motions that do not contributeto structural relaxation. If so, inter-excitation cor-relations could be present in the glasses that we pre-pared, but unobservable using any of the indicatorfunctions that we tested. Direct observation of thisphenomenon would require developing indicators ca-pable of faithfully distinguishing between trivial mo-tion unrelated to and uncorrelated with underlyingexcitations and non-trivial motion facilitated by ex- citations and responsible for structural relaxation.Whether better indicators or a more complete eval-uation of the applicability of the East model is nec-essary is left to future work. IV. CONCLUSIONS
The results in the previous section complicate thesimple picture of glassy dynamics offered by theone-dimensional East model. Although the out-of-equilibrium behavior of the relaxation time and dis-placement rate largely agree with our East-model-based expectations, the absence of spatial correla-tions between displacing particles at temperatureswell below T g stood in stark contrast to them. Rec-onciling this surprising result with our understand-ing of the East model would represent an importantstep forward in the study of glassy phenomena andthe glassy state.In the previous section, we suggested the possibil-ity that a failure of our selected indicator functionsto faithfully detect underlying excitations may bethe source of the discrepancy. If this is the case anda better indicator may be constructed, then the ap-parent contradiction is easy to resolve and our East-model-based understanding of glassy materials re-quires at most minor modifications.However, if more sophisticated indicators confirmthe results of this work, then our current picture ofglassy phenomena requires more significant revision.Specifically, the lack of spatial correlations betweenexcitations implies that the energy barriers govern-ing relaxation and aging dynamics cannot be solelyrelated to the spatial distribution of an underlyingset of excitations, as is the case for the East andother kinetically constrained models. This possibil-ity is supported by the observation in reference [35]that the kinetic stability of vapor-deposited glassescannot be encoded in the configuration of excitationson the lattice, instead requiring careful tuning of afree parameter.Given the many successes of the East model inpredicting the behavior of glass-forming materials,including in this paper, we suspect that revisions toour theory will preserve East model scaling and sim-ilarly feature a hierarchical relaxation mechanism,but include important physical details that clarifythe observed absence of inter-excitation correlationsbelow T g . ACKNOWLEDGMENTS
This work was supported by Department of En-ergy Contract No. DE-AC0205CH11231, FWP No.3CH-PHYS02. We would like to thank Shachi Katira for helpful comments on the manuscript. [1] Y. S. Elmatad, D. Chandler, and J. P. Garrahan,Journal of Physical Chemistry B , 5563 (2009).[2] Y. S. Elmatad, D. Chandler, and J. P. Garrahan,Journal of Physical Chemistry B , 17113 (2010).[3] M. D. Ediger, C. A. Angell, and S. R. Nagel, Jour-nal of Physical Chemistry , 13200 (1996).[4] C. A. Angell, Current Opinion in Solid State andMaterials Science , 578 (1996).[5] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F.McMillan, and S. W. Martin, Journal of AppliedPhysics , 3113 (2000).[6] P. G. Debenedetti and F. H. Stillinger, Nature ,259 (2001).[7] G. Biroli and J. P. Garrahan, Journal of ChemicalPhysics , 12A301 (2013).[8] J. J¨ackle and S. Z. Eisinger, Physik B – CondensedMatter , 115 (1991).[9] P. Sollich and M. R. Evans, Physical Review Letters , 3238 (1999).[10] P. Sollich and M. R. Evans, Physical Review E ,031504 (2003).[11] S. Kim, D. G. Thorpe, C. Noh, J. P. Garrahan,D. Chandler, and Y. Jung, Journal of ChemicalPhysics , 084504 (2017).[12] D. Aldous and P. Diaconis, Journal of StatisticalPhysics , 945 (2002).[13] P. Chleboun, A. Faggionato, and F. Martinelli,Journal of Statistical Mechanics: Theory and Ex-periment , L04001 (2013).[14] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C.Glotzer, and D. Chandler, Physical Review X ,021013 (2011).[15] D. T. Limmer and D. Chandler, Proceedings of theNational Academy of Sciences , 9413 (2014).[16] D. T. Limmer, Journal of Chemical Physics ,214509 (2014).[17] G. Wahnstr¨om, Physical Review A , 3752 (1991).[18] W. Kob and H. C. Andersen, Physical Review E ,4626 (1995). [19] A. S. Keys, J. P. Garrahan, and D. Chandler, Pro-ceedings of the National Academy of Sciences ,4482 (2013).[20] A. S. Keys, D. Chandler, and J. P. Garrahan, Phys-ical Review E , 022304 (2015).[21] S. Toxvaerd, U. R. Pedersen, T. B. Schrøder, andJ. C. Dyre, Journal of Chemical Physics , 224501(2009).[22] U. R. Pedersen, T. B. Schrøder, J. C. Dyre, andP. Harrowell, Physical Review Letters , 105701(2010).[23] U. R. Pedersen, N. P. Bailey, J. C. Dyre, and T. B.Schrøder, arXiv:0706.0813.[24] A. Ninarello, L. Berthier, and D. Coslovich, Physi-cal Review X , 021039 (2017).[25] J. A. Anderson, C. D. Lorenz, and A. Travesset,Journal of Computational Physics , 5342 (2008).[26] J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui,F. Spiga, J. A. Millan, D. C. Morse, and S. C.Glotzer, Computer Physics Communications ,97 (2015).[27] J. A. Anderson and S. C. Glotzer, arXiv:1308.5587.[28] D. Chandler, Introduction to Modern Statistical Me-chanics (Oxford University Press, 1987).[29] A. Widmer-Cooper, P. Harrowell, and H. Fynew-ever, Physical Review Letters , 135701 (2004).[30] O. V. Mazurin, Y. K. Startsev, and S. V. Stoljar,Journal of Non-Crystalline Solids , 105 (1982).[31] A. Alegr´ıa, E. Guerrica-Echevarr´ıa, L. Goitiand´ıa,I. Teller´ıa, and J. Colmenero, Macromolecules ,1516 (1995).[32] D. J. Plazek and J. H. Magill, Journal of ChemicalPhysics , 3038 (1966).[33] D. J. Plazek and J. H. Magill, Journal of ChemicalPhysics , 3678 (1968).[34] Z. Wojnarowska, K. L. Ngai, and M. Paluch, TheJournal of Chemical Physics , 174502 (2014).[35] R. Guti´errez and J. P. Garrahan, Journal of Sta-tistical Mechanics: Theory and Experiment2016