Slow Dynamics of the High Density Gaussian Core Model
aa r X i v : . [ c ond - m a t . s o f t ] M a y Slow Dynamics of the High Density Gaussian Core Model
Atsushi Ikeda and Kunimasa Miyazaki Institute of Physics, University of Tsukuba, Tennodai 1-1-1, Tsukuba 305-8571, Japan (Dated: October 18, 2018)We numerically study crystal nucleation and glassy slow dynamics of the one-component Gaussiancore model (GCM) at high densities. The nucleation rate at a fixed supersaturation is found todecrease as the density increases. At very high densities, the nucleation is not observed at allin the time window accessed by long molecular dynamics (MD) simulation. Concomitantly, thesystem exhibits typical slow dynamics of the supercooled fluids near the glass transition point.We compare the simulation results of the supercooled GCM with the predictions of mode-couplingtheory (MCT) and find that the agreement between them is better than any other model glassformersstudied numerically in the past. Furthermore, we find that a violation of the Stokes-Einstein relationis weaker and the non-Gaussian parameter is smaller than canonical glassformers. Analysis of theprobability distribution of the particle displacement clearly reveals that the hopping effect is stronglysuppressed in the high density GCM. We conclude from these observations that the GCM is moreamenable to the mean-field picture of the glass transition than other models. This is attributedto the long-ranged nature of the interaction potential of the GCM in the high density regime.Finally, the intermediate scattering function at small wavevectors is found to decay much fasterthan its self part, indicating that dynamics of the large-scale density fluctuations decouples withthe shorter-ranged caging motion.
PACS numbers:
I. INTRODUCTION
Essential aspects of the glass transition of the super-cooled liquids remain still elusive despite of decades ofstudy. Many theories and scenarios have been proposedto explain the dramatic slow down of the systems andthe associated growing cooperative length scales near theglass transition point [1–4]. They can explain the exper-imental results equally well or equally poorly but noneof them have been proved to be decisively better thanother. Even a satisfactory mean-field picture of the glasstransition has not been established [5, 6]. Numerical sim-ulation of simple model fluids is an ideal route to exam-ine the competing theories. Considerable efforts havebeen put forward to gain insight from the dynamical be-haviors of simple model glassformers in silico , but com-pelling answers are still lacking. There are several reasonswhy the simulation studies are not successful in sortingout numerous scenarios and theories. First, the modelsystems are more or less similar; the pair potentials ofcanonical glassformers studied in the past are exclusivelycharacterized by short-ranged strong repulsions. Exam-ples are Lennard-Jones, its WCA counterpart, soft-core,and the hard sphere potentials. Since the strong repul-sion dominates thermodynamic and dynamic propertiesof dense fluids, it is hardly surprising that the resultsfor these models are qualitatively similar [7, 8]. Studiesof a completely different class of potential systems maypotentially diversify our views and perspectives on theglass transition within the limited accessible time win-dows of the simulations. Secondly, the model systemsare not clean enough. Even the simplest class of modelglassformers (with a few exceptions [9, 10]) are inevitablybidisperse or polydisperse in order to avert the nucleation to the crystalline phase [7]. This complicates quantita-tive assessment of the simulation results. Finally, we stilllack a realistic model glassformer which conforms to themean-field picture in finite dimensions. Concept of themean-field scenario of the structural glass transition isbasically borrowed from the mean-field theory developedin the spin glass communities [3, 4, 11 ? ]. The replicatheory [12, 13] and mode-coupling theory (MCT) [14]are believed to be the static and dynamic versions of themean-field theory of the glass transition, simply becauseof their apparent resemblance to the spin-glass counter-parts. The mosaic pictures of the random first ordertransition theory has been developed as the finite di-mension version of this mean field pictures [3, 11, 15].Accumulated simulation data are not inconsistent qual-itatively from the prediction of the mean field theoriesbut the quantitative agreement between simulation re-sults and theoretical predictions are far from compelling.The best way to verify the mean-field scenario would beto take the mean-field limit by either going to higherdimensions or making the system’s interactions longer-ranged. Recently, simulations for four dimensional sys-tems have been performed [10, 16]. Results therein hintthat the dynamic heterogeneities are suppressed com-pared with three dimensional systems and agreementwith MCT moderately improves [10]. However, consider-ing the current computational abilities, it would be hardto simulate the system beyond four dimension, whereasthe upper critical dimension of the glass transition is ar-gued to be eight [17, 18]. On the other hand, few studieshave been done for realistic liquids with long-ranged par-ticle interactions [19–21].The Gaussian core model (GCM) is a candidate to dis-pel all of the above-mentioned concerns and could be anideal and clean bench to compare with various glass theo-ries. The GCM consists of the point particles interactingwith a Gaussian shaped repulsive potential [22–31]; v ( r ) = ǫ exp[ − ( r/σ ) ] , (1)where r is the interparticle separation, ǫ and σ are the pa-rameters which characterize the energy and length scales,respectively. The GCM is one of the simplest models ofthe so-called ultrasoft potential systems which are char-acterized by the bounded and long-tailed repulsive po-tential [32]. Recently, we have reported that the one-component GCM vitrifies at very high densities [33]. TheGCM or the ultrasoft particles in general have very dis-tinct and exotic properties both thermodynamically anddynamically [22–31, 33, 34], such as the re-entrant melt-ing at high densities, negative thermal expansion coef-ficient, and anomalous density dependence of the diffu-sion coefficient. There are several studies on the glasstransition of the ultrasoft particles [35–38] and it wasfound that they exhibit rich dynamical behaviors differ-ent from conventional model glassformers [35, 36]. Oneof the advantages to study the glass transition of the ul-trasoft particles is that, due to the mild repulsion tailof the potential, the density as well as the temperaturecan be used as a parameter to control the system. Ex-ploring the wide range of density–temperature param-eter space makes it easier to establish various scalinglaws, to bridge the gaps between temperature-driven or-dinary glasses and density-driven colloidal glasses, andto help unifying the concepts of the finite-temperatureglass transition and the zero-temperature jamming tran-sition [37, 38]. However, most studies in the past focusedon the relatively low density regime, where the genericnature of the glass transition is not extremely differentfrom that of the conventional model glassformers. Thesystems at low densities, including the GCM, also hadto be either polydisperse or bidisperse in order to avoidcrystallization.The GCM at very high densities is very differ-ent [33]. First of all, the system vitrifies withoutpoly(bi)dispersity. The nucleation rate systematicallydecreases as the density increases and the system startsexhibiting typical slow dynamics observed in supercooledfluids near the glass transition point. Furthermore, thedynamics is quantitatively well-described by MCT. Es-pecially, the MCT nonergodic transition point extractedfrom the simulation unprecedentedly matches with thetheoretical prediction. Besides, the violation of theStokes-Einstein relation and the amplitude of the non-Gaussian parameter, both of which is the manifestationof the heterogeneous fluctuations of dynamics, are sup-pressed. We conjecture that these facts can be attributedto the long-ranged nature of the interaction potential atthe high densities where particles are overlapped. Theseresults suggest that the high density GCM is not onlyone of the cleanest model glassformers in silico , but alsothe closest to the mean-field model.In this paper, we present thorough and complete nu-merical analysis of the nucleation and glassy dynamics of the high-density and one-component GCM. We not onlypresent the exhaustive set of the numerical results butalso provide with the new evidence which bolsters the va-lidity of MCT. Detailed analysis of thermodynamic andstructural properties of the high density GCM, such asthe phase diagram and the static structure factors arediscussed in Ref. [34]. In the previous study [33], wehave attributed the weak violation of the SE relation andsmaller non-Gaussian parameter to the suppression of thedynamic heterogeneities. We provide stronger and moredirect evidence that intermittent heterogeneous motionis suppressed by monitoring the distribution of the parti-cle displacement as a function of time. We also evaluatethe correlation functions of single and collective densityfluctuations. Surprisingly we find that dynamics of thecollective density decouple from the single particle den-sity at large length scales, where the former relaxes muchfaster than the latter. This is in stark contrast with theordinary model glassformers for which the slow glassydynamics set in over the whole length scales for bothcollective and single particle densities alike. We com-pare these simulation results with MCT predictions andfind that MCT beautifully captures the decoupling of dy-namics at the large length scales. However, we also find asubtle but noticeable disagreement of MCT from the sim-ulation results at intermediate length scales, where thenonergodic parameter (the plateau height of the two steprelaxation in the density correlators) predicted by MCTshows a weak shoulder which tends to grow as the densityincreases. This shoulder is reminiscent of those found forthe d -dimensional hard sphere glasses at large d evalu-ated from MCT [5, 6] and may be a signal of breakdownof MCT at the mean field limit.This paper is organized as follows. In Sec. II, we sum-marize the simulation method, theoretical background,and the setting of the system. The nucleation dynamicsfrom fluid to crystalline phase is discussed in Sec. III. InSec. IV, we present all simulation results on various staticand dynamical observables. Detailed analysis and care-ful comparison of the simulation results with the MCTpredictions are made. Suppression of the dynamic het-erogeneities are also discussed. Finally, Sec. IV concludesthe paper with a summary. II. PRELIMINARIESA. Simulation Methods
We investigate the dynamics of the one-componentGCM using a molecular dynamics (MD) simulation inthe
N V T ensemble with a Nos´e thermostat. The sys-tem is a cubic cell and a periodic boundary conditionis imposed. A time-reversible integrator, similar to thevelocity-Verlet method, is used with a potential cut-off at r = 5 σ [39]. Hereafter, σ , ǫ/k B , and σ ( m/ǫ ) / are takenas the units of the length, temperature, and time, respec-tively. The time step is fixed at 0.2, which is sufficiently fcc bccfluid FIG. 1: State points at which MD simulations were per-formed (crosses). Squares with solid line and filled circleswith dotted line are the solid-fluid phase boundary obtainednumerically by us [33, 34] and Prestipino et al. [26], respec-tively. The melting and freezing lines are indistinguishable atthis scale. small to conserve the Nos´e Hamiltonian during the longsimulation runs. We focus on the four densities, ρ = 0 . .
0, 1 .
5, and 2 .
0, and perform the MD simulations forvarious temperatures in the vicinity of the melting tem-perature T m . The state points which we performed simu-lation are shown in Fig. 1 along with the solid-fluid phaseboundary line [26, 33, 34]. As discussed in detail in theprevious study [34], the melting temperature, T m , at thehigh density regime ρ & T m ∝ − ρ / which was originally conjectured by Still-inger [22]. For all densities which we study, thermody-namically stable crystalline structure is bcc [33, 34]. Werun the simulations for the total run time always 50 timeslonger than the structural relaxation time. For example,the simulation time was t sim = 10 for the lowest tem-perature at ρ = 2 .
0. This was confirmed to be sufficientlylong to neglect aging effect. The first half of the simu-lation run was used for the equilibration and we usedthe trajectories of the second half for the analysis of thestationary dynamics. For each state point, five indepen-dent runs are performed and the results are obtained byaveraging over those trajectories in order to improve thestatistics. Configurations obtained from the high temper-ature simulation were used as the initial configurations.The system size is fixed at N = 3456. The simulations for N = 2000 and 9826 confirmed that the finite-size effectis negligible. B. Mode Coupling Theory
In this work, we compare our simulation results for dy-namics of the high density GCM in the supercooled statewith the prediction of MCT. In the context of the glass transition, MCT is commonly expressed as a set of theself-consistent nonlinear equations for correlation func-tions. These correlation functions are the intermediatescattering function (the correlation of the collective den-sity), F ( k, t ) ≡ h δρ ( ~k, δρ ( − ~k, t ) i /N , where δρ ( ~k, t ) isthe k -dependent density fluctuation, and the self inter-mediate scattering function or the correlation of the sin-gle particle density, F s ( k, t ) ≡ h ρ s ( ~k, ρ s ( − ~k, t ) i , where ρ s ( ~k, t ) is the density of a single particle. The time evolu-tion of F ( k, t ) is given by the generalized Langevin equa-tionΩ − ( k ) ¨ F ( k, t ) + F ( k, t ) + Z t ds M ( k, t − s ) ˙ F ( k, s ) = 0 , (2)where Ω( k ) = p k B T k /mS ( k ) is the frequency term. S ( k ) = F ( k, t = 0) is the static structure factor. M ( k, t )is the memory kernel which, according to MCT, is ap-proximated as M ( k, t ) = ρS ( k )2 k Z d~q (2 π ) V ~k ( ~q, ~k − ~q ) F ( q, t ) F ( | ~k − ~q | , t ) . (3)Here V ~k ( ~q, ~p ) ≡ { ~k · ~qc ( q )+ ~k · ~pc ( p ) } /k is the vertex, where c ( k ) = { − /S ( k ) } /ρ is the direct correlation function.In Eq. (3), we neglect the short time contribution for thememory kernel, which does not affect the slow dynam-ics. MCT predicts that F ( k, t ) undergoes the ergodic-nonergodic transition at a finite temperature, T c , belowwhich lim t →∞ F ( k, t ) = F ∞ ( k ) remains finite. F ∞ ( k ) isreferred to as the nonergodic parameter. The nonergodicparameter and T c can be evaluated by taking t → ∞ ofEqs. (2) and (3), which is expressed as F ∞ ( k ) /S ( k )1 − F ∞ ( k ) /S ( k ) = M ∞ ( k ) , (4)where M ∞ ( k ) is the long time limit of the memory kernel.As the temperature approaches to T c from above, MCTfirst predicts that F ( k, t ) exhibits the two-step relaxationbehavior characterized by a finite plateau and the slowstructural relaxation. The height of the plateau is iden-tical to F ∞ ( k ) at T = T c . The structural relaxation orthe alpha relaxation time, τ α , increases and eventuallydiverges at T c . MCT predicts that the increase of τ α is given by a power law τ α ∼ | T − T c | − γ , where γ isa system-dependent parameter which can be evaluatedfrom the MCT equation.Likewise, the MCT equation for the self intermedi-ate scattering function, F s ( k, t ), is written in the sameform as Eq. (2), but with the frequency term Ω s ( k ) = p k B T k /m instead of Ω( k ) and the self memory kernel M s ( k, t ) = ρ k Z d~q (2 π ) ( ~k · ~qk c ( q ) ) F s ( q, t ) F ( | ~k − ~q | , t )(5)instead of M ( k, t ) in Eq. (3). The MCT equation for F s ( k, t ) undergoes the nonergodic transition exactly atthe same temperature, T c , as for F ( k, t ), at least formost model systems studied in the past (see Ref. [40]for exceptions). By taking the small k -limit of theMCT equation for F s ( k, t ), we can also construct theself-consistent equation for the mean square displace-ment h R ( t ) i . MCT predicts that the self-diffusion co-efficient D ≡ lim t →∞ h R ( t ) i / t follows the power law D ∼ | T − T c | γ and vanishes at T c . Note that the powerlaw exponent γ is identical with that for τ α . In addi-tion to the MCT nonergodic transition and power law ofthe transport coefficients, MCT predicts many importantdynamical properties such as the dynamic scaling knownas von Schweidler’s law at the plateau regime (the betaregime) and the time-temperature superposition at thealpha relaxation regime [41].In order to solve the MCT equations, the static struc-ture factor, S ( k ), is required as an input. We used S ( k )obtained directly from simulations. For the numerical in-tegration of Eq. (3) and (5), we employed equally spaced400 grids with the grid spacing ∆ k = 0 . III. CRYSTALLIZATION
Ordinary simple atomic fluids nucleate to form crystalsquickly as the temperature is lowered below the meltingpoint. In this section, we analyze the crystal nucleationdynamics of the high density GCM and show that thenucleation rate systematically decreases as the densityincreases. In order to monitor the crystallization fromthe homogeneous fluid phase, we use the potential energy U and the bond order parameter q [42]. The bond orderparameter is defined by q ≡ N N X i =1 q ( i ) , (6)where q l ( i ) is the l -th bond order parameter of the i -theparticle defined by q l ( i ) = vuut π l + 1 l X m = − l | q lm ( i ) | . (7)Here q lm ( i ) is the complex bond parameter of the i -thparticle given by q lm ( i ) = 1 N b ( i ) N b ( i ) X j =1 Y lm ( ~R i − ~R j ) , (8)where ~R i is the position of the i -th particle, N b ( i ) isthe number of nearest neighbor particles around the i -thparticle, and Y lm ( ~r ) is the spherical harmonic function ofthe degree l and the order m . q is zero in the fluid phase and q ≈ . q and U of thefive representative trajectories as a function of the lapse (a)(b)(c)(d) FIG. 2: The time dependence of the bond order parameter q and potential energy U of the representative trajectoriesmeasured from the time when the system is prepared. (a) ρ = 0 . T = 2 . × − , (b) ρ = 1 . T = 2 . × − , (c) ρ = 1 . T = 2 . × − , and (d) ρ = 2 . T = 2 . × − .The short bold line in each figure indicates the time scale of τ α . of time measured from the moment when the system isprepared. At a relatively low density ρ = 0 . T = 2 . × − (Fig. 2 (a)), one observes that q ’s of all five trajectoriesabruptly increase from zero to a finite value and con-comitantly U ’s decrease. These behaviors are the hall-mark of the crystal nucleation. This figure shows thatthe nucleation initiates only after the lapse of time sev-eral times longer than the structural relaxation time τ α which is indicated by the short bold lines in the figures(the precise definition and compiled data set of τ α aregiven in Sec. IV). The degree of supersaturation definedby ∆ = 1 − T /T m at this state point is 0.43. Next, welook at the higher density ρ = 1 .
5. Five runs of q and U at T = 2 . × − are shown in Fig. 2 (b). Despite of thedeeper supersaturation (∆ = 0 .
55) and much longer sim- liquid liquidbcc hcp fcc bcc hcp fcc (a) (b)
FIG. 3: The ¯ q -¯ q correlation map for the configurationsobtained at the end of all the five simulation runs at ( ρ, T ) =(1 . , . × − ) (left panel) and ( ρ, T ) = (2 . , . × − )(right panel). Four circles are the characteristic distributionfor the bcc, hcp, fcc crystal, and fluid phase. ulation runs (over 40 τ α ) than Fig. 2 (a), q and U do notshow any sign of nucleation. Decreasing the temperaturefurther to T = 2 . × − where ∆ = 0 . ρ = 2 .
0, allfive trajectories fail to nucleate even at a very low tem-perature T = 2 . × − with the similar degree of thesupersaturation, ∆ = 0 .
6, for the whole simulation runs.In order to ensure that the nucleated samples are un-ambiguously the bcc crystal and that samples whichfailed to nucleate remain in the homogeneous fluid phase,we evaluate new parameters which were recently intro-duced by Lechner et al. [43]. They have used the two av-eraged bond order parameters ¯ q ( i ) and ¯ q ( i ) and demon-strated that the correlation map of them improves abilityto determine the crystalline structures [43, 44]. The aver-aged bond order parameter is defined by replacing q lm ( i )in Eq. (7) with the averaged value ¯ q lm ( i ) defined by¯ q lm ( i ) = 1˜ N b ( i ) ˜ N b ( i ) X k =0 q lm ( k ) , (9)where q lm ( k ) is given by Eq. (8) and the sum runs from k over all ˜ N b ( i ) neighbors of the i -th particles, including the i -th particle itself ( k = 0 in the sum). In Fig. 3, we placedall ¯ q ( i ) and ¯ q ( i ) ( i = 1 , , · · · , N ) in the correlation mapfor the configurations obtained at the end of simulationruns of the two state points ( ρ, T ) = (1 . , . × − ) and( ρ, T ) = (2 . , . × − ). The four circles represent thecharacteristic areas for the bcc, hcp, fcc crystals, andfluid phase [43]. The results for ( ρ, T ) = (1 . , . × − )show that the two trajectories remain in the fluid phasewhereas the rest formed the bcc crystal. It is clear thatno other structures are formed in the course of the sim-ulations. Note that the results for the three trajecto-ries which nucleated slightly deviate from the bcc region,which we presume is due to defects or imperfectness of (a) (b)(c) (d) FIG. 4: The radial distribution function g ( r ) (left panels)and the static structure factors S ( k ) (right panels). (a) and(b) are for ρ = 1 . T = 7 . × − (dashed line) and T = 2 . × − (solid line). (c) and (d) are for ρ = 2 . T = 7 . × − (dashed line) and T = 2 . × − (solid line).The insets of (b) and (d) are the closeup of S ( k ) at small k ’sin the semilog plot. Dotted lines in (a) and (c) are the barepotential v ( r ). the obtained crystalline structures. On the other hand,all the five trajectories for ( ρ, T ) = (2 . , . × − ) donot show any hint of crystal nucleation and the configura-tions remain completely disordered. Hereafter, we focuson the densities ρ = 1 . IV. GLASSY DYNAMICSA. Structural functions
Before discussing the slow dynamics in the supercooledstate, we summarize the fluid structures of the high den-sity GCM to demonstrate the difference from those ofconventional model glassformers. In Fig. 4, we plot theradial distribution functions g ( r ) and static structure fac-tors S ( k ) of the GCM for ρ = 1 . g ( r ) and S ( k )show typical behaviors of dense fluids characterized bythe prominent peaks near the position and the wavevec-tor corresponding to the first coordination shell. Theirpeak heights increase as the temperature decreases. Asdensity increases from ρ = 1 . g ( r ) shifts from r = 0 .
94 to 0.85 and for S ( k ) from k = 7 . v ( r ) stretchesbeyond the first coordination shell, as demonstrated inFigs. 4 (a) and (c). This considerable overlap of parti-cles imparts the character of the long-ranged interactionsystems to the high density GCM. The long-ranged na-ture also appears as the anomalously small S ( k ) at smallwavevectors. The insets of Figs. 4 (b) and (d) show that S ( k ≈ S ( k ) vanishes at k →
0. More detailed analysis of thesimulation results for the structural functions and com-parisons with the predictions of the liquid state theoryhave been reported in Ref. [34]. S ( k )’s obtained here areused in the MCT analysis discussed below. B. Mean square displacement and self intermediatescattering function
In this subsection, we evaluate various dynamic quan-tities and observe their slow dynamics, focusing on thetrajectories which did not crystallize even when deeplysupercooled. The mean square displacement h R ( t ) i ≡ N − P Ni =1 h| ~R i ( t ) − ~R i (0) | i , the self intermediate corre-lation function F s ( k, t ), and the intermediate correlationfunction F ( k, t ) are evaluated for the densities ρ = 1 . h R ( t ) i and F s ( k, t ) at severaltemperatures well below the melting temperature. Thesefigures clearly display the canonical behaviors of the su-percooled liquids near the glass transition point. For ρ = 1 .
5, we could not observe the glassy dynamics below T = 2 . × − because the crystallization intervened.At ρ = 2 .
0, all trajectories did not crystallize down tothe lowest temperature which we accessed. In Figs. 5 (a)and (c), one observes that, as the temperature is low-ered, h R ( t ) i develops the long plateau regimes followedby the usual diffusive behaviors h R ( t ) i ∝ t at longertimes. The appearance of the plateau signals the for-mation of a cage of a particle surrounded by its neigh-bors and is the hallmark of the supercooled fluid nearthe glass transition point. The value of p h R ( t ) i atthe plateau region is a measure of the sizes of the cages.They are about p h R ( t ) i ≈ ρ = 1 . ρ = 2 .
0. These values are slightly smaller than the val-ues for conventional model glassformers. For example, p h R ( t ) i ≈ . F s ( k = k max , t ) for sev-eral temperatures, where k max is the wavevector where S ( k ) show the maximum peak. F s ( k max , t ) relaxes ex-ponentially at high temperatures. As the temperaturedecreases, a plateau with a finite height appears and itstretches over longer times as the temperature decreasesfurther, while the plateau height remains almost con-stant. This two-step relaxation behavior is another hall-mark of the slow dynamic near the glass transition point.The terminal relaxation following the plateau is called the structural or alpha relaxation. We define the structuralrelaxation time τ α by F s ( k max , t = τ α ) = e − . In Fig. 6,we plot F s ( k max , t ) against the time scaled by τ α . Theresult shows that relaxation curves are collapsed at thealpha relaxation regime. This is the universal property ofthe glassy systems known as the time-temperature super-position (TTS) [41]. Furthermore, the all curves whereTTS holds are fitted by a stretched exponential function e − ( t/τ α ) β with the exponent β ≈ .
8. This value is com-parable with that for the KA model [46] and for the hardsphere mixture [47].In Fig. 7, the structural relaxation time τ α and the selfdiffusion constant defined by D ≡ lim t →∞ h R ( t ) i / t areplotted against the inverse temperature. We plotted D − and adjusted its ordinate so that the data collapses with τ α at high temperatures. For both densities, ρ = 1 . . τ α and D − drastically increase as the temperatureis lowered. Both data almost collapse to each other forthe whole temperatures except for the slight deviation atthe lowest temperature. As we shall discuss later, thisis the direct reflection of a weak violation of the Stokes-Einstein relation.So far, all simulation data show no sign of peculiar-ity in the slow dynamics of the high density GCM atthe qualitative level. They are all similar to conventionalmodel glassformers. In order to assess the properties ofthe high density GCM more quantitatively, we comparethe simulation results with the predictions of MCT. Forthis purpose, we solve the MCT equations Eqs. (2)–(5) bynumerically integrating the equations in a self-consistentmanner. As inputs, we used S ( k ) obtained numericallyin the previous subsection. First, we compute the MCTtransition temperature T c by solving Eq. (4). The resultsare T (theory) c = 2 . × − and 3 . × − for ρ = 1 . T (theory) c in order to emphasize that they areobtained by solving the MCT equations. The exponent γ ≈ . D − , τ α ∝ | T − T c | − γ with the same parameters γ and T c . We fitted D − and τ α obtained by simulation withthis MCT power law, using T c as a fitting parameter. Wedenote it as T (sim) c . By plotting D γ and τ − γα against T − ,we found that they both vanish at the same temperatureand we identified T (sim) c = 2 . × − and 2 . × − for ρ = 1 . τ α in Fig. 7 using ε ≡ − T /T (sim) c instead of 1 /T . The results for the KA model [46] arealso plotted. These data are scaled by a time unit, t ,defined by a relaxation time at the short time scale, F s ( k max , t = t ) = 0 .
95. This figure shows that the re-laxation times for both the GCM and KA model ride onthe MCT power law for the range of temperatures whichthe simulation can access. Collapse of the data of twosystems on a single power law is a reflection that the val-ues of γ ’s of both systems are close ( γ ≈ . (a) ρ = 1.5 (b) ρ = 1.5(c) ρ = 2.0 (d) ρ = 2.0 FIG. 5: h R ( t ) i ((a) and (c)) and F s ( k max , t ) ((b) and (d)). The filled circles are simulation results for ρ = 1 . T × = 7, 4, 3, 2 . . ρ = 2 . T × = 10, 7, 5, 4, 3 .
4, 3 . .
93 (lower panel). The dashed line in (c) is the mean square displacement of the KA model at T = 0 .
475 [46] shiftedto fit with the GCM’s result at the lowest temperature at long times (see text). The dashed lines in (b) and (d) are the MCTsolutions obtained using the same reduced temperatures, ε , as those for the simulation data. model [48]). This figure also demonstrates that ε is agood parameter to measure the distance from the onsetof the glassy slow dynamics for different systems. Here-after, we refer to ε as the reduced temperature. In Fig. 5(c), we plotted the simulation data of h R ( t ) i for the KAmodel at T = 0 .
475 by shifting the time unit in such away that the long time diffusive regime collapses with thedata for the GCM at T = 2 . × − and ρ = 2 . T (sim) c , where our MD simulation can access.However, there are two noticeable differences betweenthe high density GCM and conventional model glass-formers. First, the MCT transition temperature ob-tained from fitting the simulation data, T (sim) c , is un-precedentedly close to the theoretical prediction T (theory) c for the GCM. The agreement improves as the density in-creases. The deviation of T (sim) c from T (theory) c are only32 % for ρ = 1 . ρ = 2 .
0. It is in starkcontrast with the KA model for which T (sim) c = 0 .
435 and T (theory) c = 0 .
92 with the deviation of more than100% [49, 50]. The KA model at T (theory) c is still a high-temperature fluid and F s ( k, t ) decays exponentially with-out a sign of two-step relaxation. Contrarily, the GCM at T (theory) c already lies deep in the region where the plateauof F s ( k, t ) is well developed (see Fig. 5 (d)). Consider-able deviation of T (sim) c from T (theory) c for conventionalmodel glassformers is known as one of serious drawbacksof MCT. These deviations have been attributed to the ef-fect of the activated processes in the ragged energy land-scapes, which smears out the clear-cut dynamical tran-sition [51–54]. Second, the MCT parameters T c and γ obtained from fitting simulation data for τ α match verywell with that obtained from the data of D − . This is alsoin contrast with the model glassformers such as the KAmodel [46, 50] and poly(bi)disperse hard spheres [47, 55],for which T (sim) c (or the transition density ρ (sim) c ) and γ obtained from fitting the simulation data vary dependingon the observables ( τ α or D − ) and also on the compo-nents (large or small particles components of the binarysystems). These variances are partly attributed to thepresence of strong dynamic heterogeneities which decou-ple the diffusion from the structural relaxation time, aswe shall discuss in the next subsection. FIG. 6: Same as Fig. 5 (d) but plotted against t scaled by τ α . Filled circles are the simulation results for ρ = 2 . T × = 5, 4, 3 .
4, 3 .
2, 3, 2 .
93 from right to left. The solidline is a fit by a stretched exponential function.
The direct evidence that MCT works better for theGCM than any other model glassformers is the remark-able agreement of the simulated F s ( k, t ) with the MCTprediction. In Fig. 5 (b) and (d), we plotted the solutionsof MCT for exactly the same reduced temperatures ε asthe simulation data. Only free parameter is the time unit,which is determined solely from the short time dynam-ics. Long time behaviors of the MCT solution agree verywell with the simulation results. MCT also correctly pre-dicts the exponent of the stretched exponential relaxation β . The agreement is striking given that for other modelglassformers, ε (and sometimes the wavevectors as well)needs to be adjusted at each temperature to obtain a rea-sonable fit [49, 56] (an exception is the four-dimensionalsystem [10]). C. Intermediate scattering function
Next, we look at the intermediate scattering function F ( k, t ). For conventional model glassformers, it is knownthat behavior of F ( k, t ) is qualitatively the same as thatof F s ( k, t ), except for the wiggly k -dependence of the non-ergodic parameter for the former, reflecting the wigglyprofiles of the static structure factor (see the discussionbelow). Contrarily, for the high density GCM, F ( k, t )and F s ( k, t ) differ from each other considerably. Fig. 9shows F ( k, t ) at two wavevectors. Fig. 9 (a) is the re-sult at k = k max ( ≈ .
4) which is the peak position of S ( k ). There, the relaxation behavior of F ( k, t ) is verysimilar to that of F s ( k, t ), suggesting the relaxations ofboth functions at the interparticle length scales are dic-tated by the same relaxation mechanism. Fig. 9 (b) is (a)(b) FIG. 7: The temperature dependence of the structural re-laxation time (filled circles) and the inverse of the diffusioncoefficient (empty squares) for (a) ρ = 1 . ρ = 2 . τ − /γα and D /γ as a function of inverse temperature,where γ is fixed to 2.7. the result at k = 6 .
4, which corresponds to a slightlylonger length scale than the interparticle distance. Therelaxation of F ( k, t ) is very fast and shows no sign of twostep relaxation. F ( k, t ) almost fully relaxed at t ∼ h R ( t ) i appears (see Fig. 5). Thequick decays are followed by the phonon-like oscillationsand very weak tails persisting over the time scale of thestructural relaxation time. This tail vanishes at smaller k ’s. This behavior is in sharp contrast with the KAmodel, where the relaxation time at small wavevectorsis comparable with that at the interparticle distance andthe plateau heights remains finite down to very smallwavevectors [57]. These results indicate that, in the highdensity GCM, the large scale density fluctuations are de-coupled from the slow structural relaxation processes at FIG. 8: τ α /t as a function of the reduced temperature ε forthe GCM and KA model. t is the short-time relaxation timedefined by F s ( k max , t ) = 0 . the shorter length scales.In order to see this qualitative difference of F ( k, t )of the GCM more clearly, we plot the k -dependence ofthe plateau heights, or the nonergodic parameter, F ∞ ( k )and F s, ∞ ( k ) together with the MCT predictions obtainedfrom Eq. (4). In Fig. 10, we show F ∞ ( k ) /S ( k ) and F s, ∞ ( k ) at ρ = 2 . F ∞ ( k ) /S ( k ) above k max remains compatible with that of F s, ∞ ( k ), while keeping a wiggly behavior characteristicof the collective density fluctuations. Absence of slowdynamics at small k ’s is a consequence of the anomalousstructural properties inherent in the high density GCM.In the previous subsection, we discussed that the staticstructure factor at the small wavevectors, or the com-pressibility, is extremely small compared with those ofordinary model glassformers. This makes the amplitudeof the memory kernel at small k ’s negligibly small (seeEq. (3)). Consequently the large scale fluctuations decou-ple from the fluctuations at the length scales of the inter-particle distance which trigger the glassy slow dynamics.We argue that this decoupling between small and longlength scales should be commonly observed for the sys-tems with small compressibilities which are an universalfeature of the dense and long ranged interaction systemsincluding the Coulomb interaction systems as predictedin the framework of MCT [58].The nonergodic parameters in Fig. 10 exhibit anothersubtle but noticeable feature which may have relevance tofundamental problems of MCT as the mean field descrip-tion of the glass transition. Although MCT reproducesthe overall behaviors of the nonergodic parameters forboth F ∞ ( k ) /S ( k ) and F s, ∞ ( k ), its prediction systemat- (a) k = 8.4(b) k = 6.2 FIG. 9: The intermediate scattering function at (a) k = 8 . k = 6 .
2. For both panels, ρ = 2 . T × = 10, 7, 5, 4, 3.4, 3.2, 3and 2.93. The inset in (b) shows a closeup of the weak andlong tails of the main panel. ically overestimates the simulation results at the inter-mediate wavevectors (in the range of, say, 5 . k . F s, ∞ ( k ) is well fitted by a Gaussian function,whereas the MCT nonergodic parameter has a small butnon-negligible shoulder which the Gaussian function cannot fit. This shoulder is reminiscent of those observed inthe MCT solution for hard sphere glasses in large spatialdimensions [5, 6]. There, we have found that the devia-tion from the Gaussian function for F s, ∞ ( k ) increases asthe dimension d increases. This observation has lead usto conclude that MCT is not rigorously a bona fide meanfield theory [5]. This glitch of MCT which we found inone of the mean field limits, i.e. , the high d limit, couldalso show up in another mean field limit, that is, the long-ranged interaction limit, which is realized in the highdensity limit of the ultrasoft potential systems such as0 s , ∞∞ (a)(b) FIG. 10: Nonergodic parameters for the collective part F ∞ ( k ) /S ( k ) (upper panel) and self part F s, ∞ ( k ) (lower panel)of the intermediate scattering functions. Filled circles are thesimulation data and solid lines are the MCT solutions. Thedotted line in the lower panel is a fit by a Gaussian function. the GCM. This may explain the shoulder of the F s, ∞ ( k )in Fig. 10 (b). Remember that the anomalously small S ( k ) at small k ’s is also due to the long-ranged inter-action. Interestingly, this small S ( k ) may explain theanomalous shoulder of the MCT solution. By artificiallyenhancing the amplitude of S ( k ) at small k ’s by a minuteamount and plugging the modified S ( k ) into the MCTequation, we find that the nonergodic parameter F ∞ ( k )at small k ’s jumps from zero to finite values. At thesame time, the shoulder of F s, ∞ ( k ) at the intermediatewavevectors disappears and MCT’s F s, ∞ ( k ) gets closerto the simulation results. This observation implies thatthe long range interaction affects the static properties of the large length scales, which eventually amplifies the pu-tative non-Gaussian behaviors of the MCT solution. Asubtle interplay between the long and short length fluc-tuations may be quite common for the glass or/and jam-ming transition: For example, the the hyper-uniformity(vanishing S ( k ) at small k ) and diverging radial distri-bution function at the contact length r = σ are known tobe the two facets of a universal character of the jammingtransition [59]. D. Violation of Stokes-Einstein relation
For many glassformers, the Stokes-Einstein (SE) rela-tion D ≈ T /η , where η is the shear viscosity, is violatednear the glass transition point and the violation is be-lieved to be the manifestation of spatially heterogeneousdynamics which grows as the temperature is lowered [60].Indeed, MCT can not capture the SE violation due to itsmean field character. In this section, we show that theSE violation for the high density GCM is suppressed. InFig. 11 (a), we plot Dτ α for ρ = 1 . ε .Note that Dτ α instead of Dη has been plotted, because η and τ α are roughly proportional to each other. In thesame figure, we have also plotted the results for the largeand small particles for the KA model [46]. It is obviousthat the variations of Dτ α for the GCM is much weakerthan that of the KA model. Similar suppression of theSE violation was observed in the four-dimensional hardsphere system [10]. τ α was defined by F s ( k, τ α ) = e − at k = k max . Inorder to study the length scales which are relevant to theSE violation, we generalize the structural relaxation timeto the k -dependent form, τ ( k ), defined by F s ( k, τ ( k )) = e − . Note that τ α = τ ( k max ). In the small wavevectorlimit, the self intermediate scattering function behaves as F s ( k, t ) = e − Dk t . Therefore, τ ( k ) ∼ /Dk as k →
0. Inthe opposite limit, the system should behave as an idealgas, so that F s ( k, t ) = e − k B T k t /m . Thus, τ ( k ) ∝ /k as k → ∞ [61]. Fig. 11 (b) shows Dk τ ( k ) as a func-tion of k for ρ = 2 . etal . [50]. At a high temperature T = 7 . × − wherethe two-step relaxation of F s ( k, t ) is set off (see Fig. 5(d)), Dk τ ( k ) is nearly constant and almost 1 at smallwavevectors up to k max . It then decreases as k increasesfurther, followed by a turn over to a mildly increasingfunction. The decrease is a reflection of the vanishingof the cages at length scales shorter than the interparti-cle distance. The increase at larger k is a crossover tothe ideal gas limit where Dk τ ( k ) ∝ k . The qualita-tive behavior remains unchanged at T = 4 . × − , butthe drop at k & k max is more pronounced, reflecting thestronger cage effect at lower temperatures. At the low-est temperature T = 2 . × − which corresponds toabout ε ≈ . k & k max is more dramatic.Furthermore, slight positive bump at 3 . k . k max is1 k max (b)(a) FIG. 11: (a) Reduced-temperature dependence of Dτ α at ρ = 1 . Dτ α ) ref . (b) Dk τ ( k )for T × = 7 . ρ = 2 .
0. The arrow indicates k max , the first peak of S ( k ). observed. This deviation corresponds to a weak SE vio-lation observed in Fig. 11 (a). This behavior is noticeablydifferent from that for the KA model for which Dk τ ( k )significantly increases as k increase before dropping near k max [50]. E. Non Gaussian dynamics
Another good measure to monitor the extent ofthe departure from the mean field behavior is the (a)(b)(a)
FIG. 12: (a) The non-Gaussian parameter α ( t ) for T × =10, 7, 5, 4, 3.4, 3.2 3 and 2.93 at ρ = 2 .
0. (b) The temperaturedependence of the maximum value of α ( t ) at ρ = 1 . non-Gaussianity of the dynamics. At high tempera-tures, F s ( k, t ) or its real space expression, G s ( r, t ) ≡ P i h δ ( | ~R i ( t ) − ~R i (0) | − r ) i , also known as the van Hovefunction, becomes almost a Gaussian function. However,as the temperature is lowered to the supercooled regime,these function substantially deviates from the Gaussian.This deviation is also considered to be a manifestation ofdynamic heterogeneities. To quantify this, it is common2to introduce the non-Gaussian parameter defined by α ( t ) ≡ h R ( t ) i h R ( t ) i − , (10)where h R ( t ) i = N − P i h| ~R i ( t ) − ~R i (0) | i . In Fig. 12(a), we plot α ( t ) for ρ = 2 . t near or slightly before τ α whose heights increase as the temperature decreases.However, the heights of the peaks are considerably lowerthan other model glassformers at the comparable reducedtemperatures ε . Fig. 12 (b) shows the temperature de-pendence of the maximum value of the non-Gaussian pa-rameter α max for both ρ = 1 . α max of the GCM is far smallerthan that of the KA model. Furthermore, one observesthat α max for ρ = 2 . ρ = 1 .
5. These results suggest that the dynamic hetero-geneities are suppressed for the GCM and the suppressionis stronger at higher densities. This is another collateralsupport that the high density GCM is more “mean-field-like” than other glassformers.More direct evidence that the dynamics of the highdensity GCM is closer to a Gaussian process and dynamicheterogeneities are weaker can be obtained by monitoringthe probability distribution of the particle displacement r , denoted as P (log r ; t ). P (log r ; t ) is related to thevan Hove function G s ( r, t ) by [50, 62, 63] P (log r ; t ) = (ln 10)4 πr G s ( r, t ) . (11)If the dynamics is purely a Gaussian process, G s ( r, t ) alsobecomes a Gaussian function, G s ( r, t ) = (cid:18) π h R ( t ) i (cid:19) / e − r / h R ( t ) i . (12)From Eqs. (11) and (12), P (log r ; t ) becomes a functionof solely r/ p h R ( t ) i ; P (log r ; t ) = (ln 10)4 π (cid:18) r π h R ( t ) i (cid:19) / e − r / h R ( t ) i . (13)Thus, the shape of P (log r ; t ) for a Gaussian processshould be unchanged as t is varied, but only shifted ifplotted as a function of log r . The peak height shouldbe a constant value of 2 .
13. In Fig. 13, we plotted thesimulated P (log r ; t ) for ρ = 2 . T = 7 . × − ( ε ≈ .
2) and T = 2 . × − ( ε ≈ . P (log r ; t ) is almost given by Eq. (13); theshape of the function is almost Gaussian and the peakheight remains very close to 2 .
13 over the long time. Onthe other hand, the low temperature result in Fig. 13(b) shows that the peak height of the function becomeslower and the width becomes slightly larger at t ∼ τ α . (a)(b) l og l og FIG. 13: The probability distribution of the logarithm of theparticle displacements. (a) The results for T = 7 . × − and ρ = 2 .
0. From left to right, t = 44, 180, 512, 1024, 2896and 5792. (b) The results for T = 2 . × − and ρ = 2 . t = 500, 32000, 181000, 362000, 1448000,and 4096000. This non-Gaussian behavior at the beta to alpha relax-ation time regime is a common properties of P (log r ; t )at a mildly supercooled state. Note that, however, the ex-tent of the non-Gaussianity shown in Fig. 13 (b) is muchweaker than that of other glassformers such as the KAmodel [50]. P (log r ; t ) for typical model glassformersis known to split into the binodal shape at low temper-atures, corresponding to a separation of the constituentparticles into the mobile and immobile ones. This is oneof the most salient feature of the dynamic heterogeneities.The peak of P (log r ; t ) in Fig. 13 (b) does not show anyhint to split into the binodal shape. P (log r ; t ) of theKA model at ε = 0 .
08 ( T = 0 .
47 in the LJ unit), a com-3parable reduced temperature as that of Fig. 13 (b), iscompletely separated to the two peaks, corresponding tothe distribution of mobile and immobile particles. Thedecrease of the peak height of P (log r ; t ) in Fig. 13 (b)is compatible with that of the KA model at much highertemperature, ε = 0 .
38 ( T = 0 . V. SUMMARY AND OUTLOOK
In this paper, we presented the detailed analysis ofdynamics of the high density GCM. The results are sum-marized below.(i) The crystal nucleation becomes slower as the den-sity increases. Analysis of the two orientational bond or-der parameters, ¯ q and ¯ q , reveals that the crystal struc-ture is bcc at all densities beyond the reentrant point.(ii) The system which failed to crystallize shows cleartwo-step and stretched exponential relaxation in the(both self and collective) intermediate scattering func-tions, which is the hallmarks of the supercooled fluidnear the glass transition point. All dynamical proper-ties which we have analyzed are well described by MCT.First, the temperature dependence of the diffusion coef-ficient and the structural relaxation time is well fittedby the MCT power law. The parameter T (sim) c used tofit the simulation data is unprecedentedly close to thetheoretical prediction. The time dependence of the selfintermediate scattering function F s ( k, t ) is well fitted byMCT, using the reduced temperature ε as a sole pa-rameter. Furthermore, the nonergodic parameters forboth collective and self intermediate scattering functions, F ∞ ( k ) and F s, ∞ ( k ), are well described MCT. Here wefind two noticeable differences from the typical glassform-ers. First, the shape of F ∞ ( k ) is qualitatively differentfrom F s, ∞ ( k ) at small wavevectors regime, where F ( k, t )decays very fast and the nonergodic parameter vanishes,whereas F s ( k, t ) decays very slowly and its nonergodicparameter remains finite down to k = 0. We conjec-ture that this decoupling of the collective density dy-namics from the single particle dynamics is universal forthe systems with the long-ranged interactions. This indi-cates that the large-scale density fluctuation is decoupledto the slow structural relaxation processes. Similar de-coupling has been predicted from the MCT analysis ofthe systems with the power law interactions v ( r ) ∼ /r n with small n [58]. Second, the agreement between MCTand simulation for F s, ∞ ( k ) is satisfactory but conceiv-ably worse than those for other model glassformers suchas the KA model [47, 57]. We found a weak shoulderat the intermediate wavevectors. This shoulder is remi- niscent of those found in the MCT analysis of the hardsphere glasses at the high dimensions [5]. We conjecturethat the anomalous shoulders are the deficiency of MCTwhich appears only at the mean-field limit.(iii) Dynamic heterogeneities are suppressed in thehigh density GCM. The SE violation is very weak andthe peak height of the non-Gaussian parameter is muchlower than the conventional model glassformers at thecomparable reduced temperatures. The weak dynamicheterogeneities of the high density GCM was most obvi-ous from the observation of the probability distribution ofthe particle displacement P (log r ; t ). We find no obvi-ous change in the shape of P (log r ; t ) which remains al-most Gaussian, though the width slightly widens aroundthe beta to alpha relaxation time regimes. Even at thelowest reduced temperature, at which the typical modelglassformers exhibit the very clear binodal distributionof mobile and immobile particles, due to the growing dy-namic heterogeneities, the probability distribution of theGCM remains to be a single peak function.We conclude that the high density GCM is the idealmodel system to study the glass transition. It is not onlythe cleanest glass model in that it is the one-componentsystem. But it is also the closest to the “mean-field” model in that dynamic heterogeneities are strongly sup-pressed and the way how MCT predicts simulation resultsis synchronized with the way it does for the high dimen-sional systems. The mean-field nature comes from thelong-range nature of the interaction potential, which iscaused by the overlapping of the particles at the high den-sities. Both the excellent agreement with MCT and smalldeviation from MCT (the shoulder of F s, ∞ ( k )) also leadus to reconsider the validity of MCT as the the mean fieldtheory of the glass transition. Mean-field models of theglass transition have been proposed and analyzed by tak-ing the long-range limit of the interactions [20] but it hasnever been realized in the simulation box. The anothermean field limit, i.e. , the high dimension limit, is anotherinteresting challenge but given the current CPU power,going beyond d = 5 would be unrealistic. In this sense,the high density GCM might be the first realistic fluidmodel which may be able to bridge the gap between thefinite dimensional system with the mean-field limit. It istempting to consider the high density limit of the GCMwhere the small parameter 1 /ρ may make the analyti-cal treatment of especially the static/thermodynamic pa-rameters tractable and leads us the exact mode-couplingtheory (or alike). Acknowledgments
This work is partially supported by Grant-in-Aid forJSPS Fellows (AI), KAKENHI; [1] P. G. Debenedetti and F. H. Stillinger, Nature , 259(2001).[2] A. Cavagna, Phys. Rep. , 51 (2009).[3] G. Biroli and J. P. Bouchaud, arXiv:0912.2542.[4] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti,and W. van Saarloos, eds.,
Dynamical Heterogeneities inGlasses, Colloids, and Granular Media (Oxford Univer-sity Press, Oxford, 2011).[5] A. Ikeda and K. Miyazaki, Phys. Rev. Lett. , 255704(2010); ibid. , 049602 (2011).[6] B. Schmid and R. Schilling, Phys. Rev. E , 041502(2010); R. Schilling and B. Schmid, Phys. Rev. Lett. ,049601 (2011).[7] H. C. Andersen, Proc. Natl. Acad. Sci. U. S. A. , 6686(2005).[8] L. Berthier and G. Tarjus, Phys. Rev. Lett. , 170601(2009).[9] F. Sausset and G. Tarjus, Phys. Rev. Lett. , 065701(2010).[10] P. Charbonneau, A. Ikeda, J. A. van Meel, andK. Miyazaki, Phys. Rev. E , 040501(R) (2010).[11] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes,Phys. Rev. A , 1045 (1989).[12] M. M´ezard and G. Parisi, Phys. Rev. Lett. , 747(1999).[13] G. Parisi and F. Zamponi, Rev. Mod. Phys. , 789(2010).[14] W. G¨otze, ”Complex Dynamics of Glass-Forming Liq-uids” (Oxford University Press, Oxford, 2009).[15] V. Lubchenko and P. G. Wolynes, Annu. Rev. Phys.Chem. , 235 (2007).[16] J. D. Eaves and D. R. Reichman, Proc. Natl. Acad. Sci.U. S. A. , 15171 (2009).[17] G. Biroli and J.-P. Bouchaud, J. Phys.: Condens. Matter , 205101 (2007).[18] G. Biroli, J.-P. Bouchaud, K. Miyazaki, and D. R. Re-ichman, Phys. Rev. Lett. , 195701 (2006).[19] E. Zaccarelli, F. Andreev, Stefan Sciortino, and D. R.Reichman, Phys. Rev. Lett. , 195701 (2008).[20] V. S. Dotsenko, J. Stat. Phys. , 823 (2004); V. S. Dot-senko and G. Blatter, Phys. Rev. E , 021502 (2005).[21] R. Mari and J. Kurchan, arXiv:1104.3420.[22] F. H. Stillinger, J. Chem. Phys. , 3968 (1976); F. H.Stillinger and T. A. Weber, ibid. , 3837 (1978); F. H.Stillinger and T. A. Weber, ibid. , 4879 (1979); F. H.Stillinger, Phys. Rev. B , 299 (1979).[23] F. H. Stillinger and D. K. Stillinger, Physica A , 358(1997).[24] A. Lang, C. N. Likos, M. Watzlawek, and H. L¨owen, J.Phys.: Condens. Matter , 5087 (2000).[25] A. A. Louis, P. G. Bolhuis, and J. P. Hansen, Phys. Rev.E , 7961 (2000).[26] S. Prestipino, F. Saija, and P. V. Giaquinta, Phys. Rev. E , 050102(R) (2005); S. Prestipino, F. Saija, and P. V.Giaquinta, J. Chem. Phys. , 144110 (2005).[27] B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, andC. N. Likos, Phys. Rev. Lett. , 045701 (2006).[28] P. Mausbach and H.-O. May, Fluid Phase Equilibria ,17 (2006).[29] C. E. Zachary, F. H. Stillinger, and S. Torquato, J. Chem.Phys. , 224505 (2008). [30] W. P. Krekelberg, T. Kumar, J. Mittal, J. R. Errington,and T. M. Truskett, Phys. Rev. E , 031203 (2009);W. P. Krekelberg, M. J. Pond, G. Goel, V. K. Shen, J. R.Errington, and T. M. Truskett, ibid. , 061205 (2009);M. J. Pond, W. P. Krekelberg, V. K. Shen, J. R. Erring-ton, and T. M. Truskett, J. Chem. Phys. , 161101(2009);M. J. Pond, J. R. Errington, and T. M. Truskett,arXiv:1101.1982.[31] L. A. Shall and S. A. Egorov, J. Chem. Phys. ,184504131 (2010).[32] C. N. Likos, Phys. Rep. , 267 (2001); Soft Matter ,478 (2006).[33] A. Ikeda and K. Miyazaki, Phys. Rev. Lett. , 015701(2011).[34] A. Ikeda and K. Miyazaki, (unpublished).[35] G. Foffi, F. Sciortino, P. Tartaglia, E. Zaccarelli, F. L.Verso, L. Reatto, K. A. Dawson, and C. N. Likos, Phys.Rev. Lett. , 238301 (2003).[36] E. Zaccarelli, C. Mayer, A. Asteriadi, C. N. Likos,F. Sciortino, J. Roovers, H. Iatrou, N. Hadjichristidis,P. Tartaglia, H. L¨owen, et al., Phys. Rev. Lett. ,268301 (2005); C. Mayer, E. Zaccarelli, E. Stiakakis,C. N. Likos, F. Sciortino, A. Munam, M. Gauthier,N. Hadjichristidis, H. Iatrou, P. Tartaglia, et al., NatureMaterials , 780 (2008).[37] L. Berthier and T. A. Witten, Europhys. Lett. , 10001(2009); L. Berthier and T. A. Witten, Phys. Rev. E ,021502 (2009).[38] L. Berthier, A. J. Moreno, and G. Szamel, Phys. Rev. E , 060501(R) (2010).[39] D. Frenkel and B. Smit, ”Understanding Molecular Sim-ulation” (Academic Press, 2001).[40] T. Voigtmann, arXiv:1010.0440.[41] K. Binder and W. Kob, ”Glassy Materials and DisorderedSolids” (World Scientific, Singapore, 2005).[42] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys.Rev. B , 784 (1983).[43] W. Lechner and C. Dellago, J. Chem. Phys. , 114707(2008).[44] T. Kawasaki and H. Tanaka, Proc. Natl. Acad. Sci. U. S.A. , 14036 (2010).[45] S. Ichimaru, Rev. Mod. Phys. , 1017 (1982).[46] W. Kob and H. C. Andersen, Phys. Rev. Lett. , 1376(1994); Phys. Rev. E , 4626 (1995); ibid. , 4134(1995).[47] G. Foffi, W. G¨otze, F. Sciortino, P. Tartaglia, andT. Voigtmann, Phys. Rev. E , 011505 (2004).[48] E. Flenner and G. Szamel, Phys. Rev. E , 031508(2005).[49] W. Kob, M. Nauroth, and F. Sciortino, J. Non-Cryst.Solids , 181 (2002), .[50] E. Flenner and G. Szamel, Phys. Rev. E , 011205(2005).[51] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature , 554 (1998).[52] Y. Brumer and D. R. Reichman, Phys. Rev. E , 041202(2004).[53] P. Mayer, K. Miyazaki, and D. R. Reichman, Phys. Rev.Lett. , 095702 (2006).[54] S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes,Proc. Natl. Acad. Sci. U. S. A. , 16077 (2008). [55] S. K. Kumar, G. Szamel, and J. F. Douglas, J. Chem.Phys. , 214501 (2006).[56] T. Voigtmann, A. M. Puertas, and M. Fuchs, Phys. Rev.E , 061506 (2004).[57] T. Gleim, W. Kob, and K. Binder, Phys. Rev. Lett. ,4404 (1998).[58] S. Shiroiwa, A. Ikeda, and K. Miyazaki, unpublished .[59] S. Torquato and F. H. Stillinger, Phys. Rev. E , 041113(2003); A. Donev, F. H. Stillinger, and S. Torquato, Phys.Rev. Lett. , 090604 (2005). [60] M. D. Ediger, Annu. Rev. Phys. Chem. , 99 (2000).[61] J. P. Hansen and I. R. McDonald, ”Theory of simpleliquids” (Academic Press, 1986).[62] M. E. Cates, M. Fuchs, K. Kroy, W. C. K. Poon, andA. M. Puertas, J. Phys.: Condens. Matter16