Connecting real glasses to mean-field models
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov Connecting real glasses to mean-field models
Ujjwal Kumar Nandi,
1, 2
Walter Kob, a) and Sarika Maitra Bhattacharyya
1, 2, b) Polymer Science and Engineering Division, CSIR-National Chemical Laboratory, Pune-411008,India Academy of Scientific and Innovative Research (AcSIR), Ghaziabad 201002,India Laboratoire Charles Coulomb and CNRS, University of Montpellier, Montpellier F-34095,France (Dated: 1 December 2020)
We propose a novel model for a glass-forming liquid which allows to switch in a continuous manner froma standard three-dimensional liquid to a fully connected mean-field model. This is achieved by introducing k additional particle-particle interactions which thus augments the effective number of neighbors of eachparticle. Our computer simulations of this system show that the structure of the liquid does not changewith the introduction of these pseudo neighbours and by means of analytical calculations, we determine thestructural properties related to these additional neighbors. We show that the relaxation dynamics of thesystem slows down very quickly with increasing k and that the onset and the mode-coupling temperaturesincrease. The systems with high values of k follow the MCT power law behaviour for a larger temperaturerange compared to the ones with lower values of k . The dynamic susceptibility indicates that the dynamicheterogeneity decreases with increasing k whereas the non-Gaussian parameter is independent of it. Thus weconclude that with the increase in the number of pseudo neighbours the system becomes more mean-field like.By comparing our results with previous studies on mean-field like system we come to the conclusion that thedetails of how the mean-field limit is approached are important since they can lead to different dynamicalbehavior in this limit. I. INTRODUCTION:
The details of the relaxation dynamics of glassy sys-tem and the properties of the glass has been and con-tinues to be in the focus of an intense research activ-ity . These investigations are motivated by the fact thatglasses are not only important for many daily and tech-nological applications but are also an intellectual chal-lenge for fundamental studies since so far there is notheoretical framework that is able to give a satisfactorydescription of the unusual properties of glassy systemsand glasses. Although there are sophisticated mean-fieldtheories, like the mode-coupling theory (MCT) of theglass transition , or the random first order transitiontheory , that are able to give in some cases a surpris-ingly good description of real glass former , these ap-proaches still have many flaws since they fail to give areliable description of many features of glass-forming sys-tems opening thus the door to other approaches that at-tempt to describe glassy systems . Note that thesetheories are mean-field in nature, whereas the experi-ments and computer simulation studies are three or lowerdimensional systems. Moreover, it has been found thatMCT, although expected to be mean-field in nature, doesnot become exact even at high dimensions , a flawwhich might, however, be related to the approximationsused to describe the structure of the liquid in high di-mensions. Thus it is important to understand how these a) Electronic mail: [email protected] b) Electronic mail: [email protected] theories are connected to real glass-forming systems andhow the properties change as the mean-field character ofthe system is modified. To establish such a connection itis useful to study systems whereby varying a parameterone can go from d dimensional system to mean-field (MF)system. In the past various possibilities have been pro-posed to take this limit, see Ref. 23 for an overview, butmost of them do have some drawbacks that prevent toreach a solid understanding how three-dimensional (3d)and MF systems are related to each other .One interesting model that allows approaching the MFlimit in a continuous manner has been proposed by Mariand Kurchan (MK) . The MK-model is a hard-spheresystem in which the interaction range between two par-ticles i and j is a random variable with a variance thatallows switching from a standard three-dimensional sys-tem to MF like system. For this model, it is found thatwith increasing interaction range the Stokes-Einstein re-lation holds down to lower temperatures and that thedynamic heterogeneity of the system, measured by thefour-point susceptibility and non-Gaussian parameter,decreases. The increase in interaction range also makesthe system follow MCT like behaviour for a larger rangein temperature. Although all these results indicate thatthe MK model can indeed be used to study the transitionfrom 3d to MF, there are certain features of the modelthat are disturbing. First of all, the structural proper-ties of the system becomes very different from the one ofa normal liquid if the MF limit is approached in that,e.g., the radial distribution function becomes gas-like.Related to this is the fact that the three-point correlationfunctions vanish. As a consequence one looses the prop-erty that nearest neighbors can cage a tagged particle, anotion that is fundamental for the slowing down of thedynamics in real glass-forming systems . Secondly, themaximum attainable packing fraction diverges in the MFlimit, a behavior that is very different from the one foundin finite dimensions. Some of these oddities are avoidedif one considers models on a lattice . However, latticemodels, notably kinetic Ising models with non-conservedparticle density, do have the drawback that it is not ob-vious to what extent their relaxation dynamics is relatedto any off-lattice systems. As a consequence one has tobe cautious when applying results from lattice models todescribe the dynamics of real systems.Another approach to connect the properties of 3d sys-tems with the MF behavior has been proposed in a seriesof papers by Miyazaki and coworkers who have studiedthe properties of the Gaussian-Core-Model (GCM) .Due to the long interaction range, each particle has alarge number of neighbours, and hence the system can beexpected to be MF like. These authors showed that com-pared to the (short-ranged) Kob-Andersen (KA) model ,in the GCM the Stokes-Einstein relation is followed till alower temperature regime and that the relaxation dynam-ics shows a qualitatively better agreement with the MCTpredictions . Furthermore, it was found that the GCMshows less dynamic fluctuation and that activated pro-cesses are suppressed , in agreement with recent studiesof the thermodynamic properties of this system .A further possibility to connect the properties of lowdimensional systems with the MF predictions is to con-sider systems with increasingly higher dimensions. Sen-gupta et al. have studied the properties of some stan-dard glass formers in 2, 3, and 4 dimensions and foundthat with increasing dimensionality the breakdown ofthe Stokes-Einstein relation becomes less pronounced andthat the dynamical heterogeneity decrease . Charbon-neau et al. have studied systems up to 6 dimensionsand found that the shape of the cage does not becomeGaussian-like, as expected from MF , showing that theapproach to this limit might be more complex than ex-pected.In the present paper we introduce a simple approachthat allows crossing over in a continuous manner from anormal 3d liquid to a MF system. In practice we do thisby increasing for each particle the number of particles itcan interact with, thus increasing the effective interactionof the particle with the rest of the system. In contrast tothe studies discussed above, our method does not modifyin a significant manner the local structure of the liquideven when the MF limit is reached, i.e. the structureis always similar to the one of the 3d system. So thisallows us to study how increasing connectivity affects therelaxation dynamics, without modifying in a noticeablemanner the structure, and hence to probe the dynamicsupon approaching the MF limit.The rest of the paper is organized as follows: The sys-tem and simulation details are described in Sec. II. InSec. III, we present the result while in Sec. IV we sum-marize and conclude. II. DETAILS OF SYSTEM AND SIMULATIONS
As mentioned in the Introduction, our system is givenby N particles that interact with each other via a stan-dard short-range potential. In addition, each particleinteracts also with “pseudo neighbors”, i.e. particles thatare not necessarily close in space. Hence the total inter-action potential of the system is given by U tot ( r , ...r N ) = N X i =1 N X j>i u ( r ij ) + N X i =1 k X j =1 u pseudo ( r ij ) (1)= U + U pseudo k . (2)The first term on the right-hand side is the regular in-teraction between particles while the second term is theinteraction each particle has with its pseudo neighbours.Here we consider the case that the regular interaction de-scribes a binary Lennard-Jones (LJ) system, with 80% ofthe particles of type A and 20% of the particles of typeB. Thus the interaction between the particles i and j isgiven by u ( r ij ) = 4 ǫ ij h(cid:16) σ ij r ij (cid:17) − (cid:16) σ ij r ij (cid:17) i , (3)where r ij is the distance between the particles, σ ij is theeffective diameter of the particle and ǫ ij is the interac-tion strength. We use σ AA and ǫ AA as the unit of lengthand energy, setting the Boltzmann constant k B = 1.The values of the other parameters are given in Ref. 9,i.e. σ AB = 0 . σ BB =0.88, ǫ AB =1.5, and ǫ BB =0.5, achoice which makes this binary system to be a good glass-former. This potential is cut and shifted at r c = 2 . σ ij .The masses are m A = m B = 1 and time is expressed inunits of p m A σ AA /ǫ AA .The interaction potential with the pseudo neighboursis modelled in terms of a modified LJ potential, u pseudo ( r ij ) = u ( r ij − L ij ) (4)= 4 ǫ ij h(cid:16) σ ij r ij − L ij (cid:17) − (cid:16) σ ij r ij − L ij (cid:17) i , (5)where L ij is a random variable defined below. In oursimulations we impose the restriction that any two par-ticles interact either via u ( r ij ) or via u pseudo ( r ij ). Thiscondition determines how for a given configuration equi-librated with the potential u the pseudo neighbors andthe values L ij are chosen: Taking this configurations weselect for each particle, i , k random numbers L ij in therange r c ≤ L ij ≤ L max , where L max ≤ L box / − r c , with L box the size of the simulation box. (The distribution ofthese random variables will be denoted by P ( L ij ) andin the following, we will consider the case that the dis-tribution is uniform.) Subsequently we choose k distinct r g (r) g AA (r) for k=0g AB (r)for k=0g BB (r) for k=0g AA (r) for k=28g AB (r) for k=28g BB (r) for k=28T=0.9 FIG. 1: The partial radial distribution functions for k = 0 and k = 28 at T = 0 .
9. The structure remainsinvariant under the introduction of the pseudoneighbours.particles j with r ij > r c and use the L ij to fix perma-nently the interaction between particles i and j . Thisprocedure thus makes that each particle i interacts notonly with the particles that are within the cutoff distancebut in addition to k particles that can be far away. Thesystem, as defined here, can then be simulated using astandard simulation algorithms.The molecular dynamics (MD) simulation have beendone using N = 2744 particles. We have performed con-stant volume, constant temperature simulations (velocityrescaling) at density ρ = 1 .
2, thus L box = 13 . t = 0 . L max we havetaken 4.0, slightly below the maximum value of 4.09. Wehave simulated four different systems with the number ofpseudo neighbours, k = 0 , , , and 28. III. RESULTSA. Structure of the liquid
To start, we discuss the effect of the pseudo neighbourson the structure of the liquid. In Fig. 1 we show the threepartial radial distribution function, g αβ ( r ) with α, β ∈{ A, B } , for the k = 0 and the k = 28 systems. Thetemperature is T = 0 .
9, which for the k = 0 system isslightly above the onset temperature, see Ref. 9, whilefor the k = 28 system it corresponds to a state at whichthe system is already rather viscous (see below). Thegraph shows that the radial distribution functions for thetwo systems overlap perfectly well, i.e. the structure isindependent of k for this value of k . Thus this indicatesthat the interactions due to the pseudo neighbours donot affect the local structure of the system, one of thereasons for our choice of the interactions of the model.To probe whether the structure of the liquid on alarge scale is influenced by the introduction of the pseudo q -0.500.511.522.53 S ( q ) S AA (q) for k=0S AB (q) for k=0S BB (q) for k=0S AA (q) for k=28S AB (q) for k=28S BB (q) for k=28 T=1.0
FIG. 2: The partial structure factors for k = 0 and k = 28 at T = 1 .
0. Similar to what we have obtained inthe radial distribution function, the structure remainsinvariant under the introduction of the pseudoneighbours.neighbors we have calculated the partial static structurefactors and show them in Fig. 2 for the case of k = 0 and k = 28. Since the two sets of curves match each otherperfectly well, we can conclude that also the large scalestructure is not influenced by the additional neighbors. B. Static properties of the pseudo neighbors
In this subsection, we characterize some of the struc-tural properties of the pseudo neighbors with respect toa tagged particle.To start, we first calculate the probability P L that agiven pseudo neighbor j interacts with the tagged particle i , where L = L ij . Neglecting the indirect interactions(via the direct neighbors) between the tagged particleand the pseudo neighbor one can express P L as P L = R V acc d r e − βu ( r − L ) y ( r ) R V acc d r e − βu ( r − L ) . (6)Here β = 1 /k B T , V acc is the volume accessible to thepseudo neighbor, and y ( r ) is a step function that takesinto account that the potential is cut off at 2.5 σ αβ ,i.e. y ( r ) = 1 if L ≤ r ≤ L + 2 . σ αβ and y ( r ) = 0 forall other values of r . The volume integrals in Eq. (6) canbe decomposed into a spherical part that is contained in-side the cubic box, and the rest. The latter volume isgiven by ∆ V = L − π (cid:16) L box (cid:17) (7)= L (1 − π . (8)A spherical integration in Eq. (6) gives then T P , k e / k k=4k=12k=28from Eqs.(11) and (12) P , k e / k FIG. 3: Probability that a pseudo neighbour is withinthe interaction range as a function of temperature. Thepink line is the theoretical prediction from Eqs. (11)and (12). Inset: Same quantities extending thetemperature range to T = 0. The theoretical curveshows a sigmoidal shape. P L = R L + r c L dr r e − βu ( r − L ) R L box / L dr r e − βu ( r − L ) + ∆ V . (9)Note that in the above expression, L = L ij is fixed.Hence for a distribution of L , the probability of findinga pseudo neighbour within the interaction range of thetagged particle is given by P = Z L max r c dL P ( L ) R L + r c L dr r e − βu ( r − L ) R L box / L dr r e − βu ( r − L ) + ∆ V . (10)In the numerator we make the substitution r ′ = r − L which allows to interchange the two integrals: P = Z r c dr ′ Z L max r c dL P ( L ) ( r ′ + L ) e − βu ( r ′ ) R L box / L dr r e − βu ( r − L ) + ∆ V . (11)We thus find that this probability is independent of k ,a result that is reasonable since we have neglected anycorrelations between the pseudo neighbors. Also notethat P depends on the interaction potential via u ( r ) and r c . For a binary system, we can generalize this calcula-tion to obtain the partial probabilities P αβ and then thetotal probability is given by P = x A P AA + 2 x A x B P AB + x B P BB , (12)where x α is the concentration of species α . In the simula-tion, this probability can be obtained by calculating theratio k e /k , where k e is the number of pseudo neighborsthat have a non-zero interaction with the tagged particle. r´ -0.500.511.522.53 g AA p s e udo (r´) k=4k=12k=28Eq. (16)exp(- β u(r´)) T=1.0
FIG. 4: Radial distribution function for pseudoneighbours from simulations at T = 1 . k = 4 , k . The solid line is theresult from the theoretical expression given by Eq. (16).The dashed line is the theoretical prediction from thebare potential.In Fig. 3 we show the temperature dependence of P asobtained from Eqs. (11) and (12) (solid line) and compareit with the corresponding quantity k e /k determined fromthe simulations (symbols). One recognizes that k e /k is asexpected independent of k and that the simulation datamatches perfectly well the theoretical prediction given byEqs. (11) and (12). Note that at the lowest temperaturesat which we could equilibrate the systems for the differ-ent value of k the probability is around 0.3, i.e. for theglassy dynamics we will discuss below only a relativelysmall part of the pseudo neighbors are actually interact-ing with the tagged particle. The inset of the figure showsthat P becomes 0.5 at around T = 0 .
4, a temperature atwhich already the k = 0 system is very viscous , and for T → j with respect to a tagged particle i we canconsider the corresponding radial distribution function g pseudo ( r ′ ) = ρ k πr N X i =1 k X j ( i ) h δ ( r ′ −| r i − r j | + L ij ) i , (13)where in the second sum the index runs over the pseudoneighbors of the tagged particle i and ρ k is the averagepseudo neighbour density, ρ k = Z L max r c k P ( L ) V − πL dL , (14)where V is the total volume of the system.To calculate g pseudo ( r ) analytically we can make use ofour result for P given by Eqs. (11) and (12). The number k e of pseudo neighbours within the interaction range canbe expressed in terms of g pseudo ( r ′ ) as k e = ρ k Z r c dr ′ g pseudo ( r ′ ) Z L max r c dL P ( L )4 π ( r ′ + L ) . (15)Since k e can also be written as k e = k × P we get,using Eq. (11) and Eq. (15) g pseudo ( r ′ ) ρ k Z L max r c dL P ( L )4 π ( r ′ + L ) = k Z L max r c dL P ( L ) ( r ′ + L ) e − βu ( r ′ ) R L box / L drr e − βu ( r − L ) + ∆ V (16)from which one obtains directly g pseudo ( r ′ ). Note that g pseudo ( r ′ ) is independent of k , since ρ k is directly pro-portional to k , see Eq. (14).Fig. 4 shows the radial distribution function g pseudo ( r ′ )from the simulations of three different values of k (sym-bols) and we recognize that, as predicted by Eq. (16)the function is indeed independent of k . We have alsoincluded the analytical result from Eq. (16) and we seethat the theory describes perfectly well the simulationdata, thus demonstrating that the approximation thatthe structure of the pseudo neighbors can be obtainedwell by the bare interaction with the tagged particle isvery accurate, at least for the k values considered in thepresent work. We also note that since one has the re-lation g pseudo ( r ′ ) = exp( − βu ( r ′ )), which can be derivedfrom Eq. (16), the function g pseudo ( r ′ ) can also be ob-tained directly from the bare interaction potential u ( r ′ )as shown in Fig.4.Within the standard theory of liquids, the radial distri-bution function allows to obtain the potential energy .Due to the presence of the pseudo neighbors this is nolonger possible, and thus the usual expression has to bemodified as follows. (Note that in the following we givethe expressions for a one-component system. For thebinary system considered here, one will have to do thesum over the various partials.) Since the potential en-ergy of the system has two contributions, one is the reg-ular neighbour and the other the pseudo neighbour (seeEq. (1)), the total potential energy U tot is given by, U tot N = ρ Z ∞ u ( r ) g ( r )4 πr dr + ρ k Z ∞ u ( r ) g pseudo ( r ) Z L max r c P ( L )4 π ( r + L ) dLdr. (17)At this stage it is useful to introduce an “effective ra-dial distribution” function g eff ( r ) by defining ρ eff g eff ( r ) = ρg ( r ) + ρ k g pseudo ( r ) R L max r c P ( L )( r + L ) dLr , (18) where the effective particle density is given by ρ eff = ρ + ρ k . (19)Note that since ρ k increases linearly with k , for large k the density ρ eff is dominated by ρ k and hence in thatlimit g eff will be directly proportional to g pseudo ( r ).Using g eff ( r ) we now can express the total potential en-ergy of the system as a function of the radial distributionfunction g eff ( r ): U tot N = ρ eff Z ∞ u ( r ) g eff ( r )4 πr dr . (20)In Fig. 5 we present g eff ( r ) for the A-A correlation fordifferent values of k . Since the regular radial distribu-tion function g ( r ) is independent of k (see Fig. 1) and g pseudo ( r ) can be calculated analytically from Eq. (16) itis possible to obtain g eff for arbitrary values of k . Thegraph shows that with increasing k , the radial distribu-tion function loses its characteristic structure with themultiple peaks and converges toward a distribution thathas a single peak at r = 1. This result can be understooddirectly from Eq. (18) since for large k the first term onthe right-hand side vanishes (if divided by ρ eff ) while thesecond term is g pseudo ( r ) multiplied by an r − dependentfactor that is independent of k . So we see that in thelarge k limit the effective radial distribution function de-velops a dominant sharp peak at a finite distance. Withdecreasing temperature, this peak increases since mostof the pseudo neighbors will condensate at the optimaldistance L ij . It is this growing peak that signals theincreasing number of constraints in the system which in-duce the slowing down of the relaxation dynamics. Thisloss of structure of the radial distribution function is atypical signature of mean-field-like systems, such as thehard-sphere system of Ref. [23]. (However, unlike theresults in the present study, in the hard-sphere systemthere is no peak at r = 1.) C. Relaxation dynamics
We now analyze how the presence of the pseudo neigh-bours affects the relaxation dynamics. To characterizethis dynamics we consider the self part of the overlapfunction Q ( t ) and the mean squared displacement (MSD)of a tagged particle, ∆ r ( t ). The former observable is de-fined as Q ( t ) = 1 N N X i =1 h ω ( | r i ( t ) − r i (0) | ) i , (21)where the function ω ( x ) is 1 if 0 ≤ x ≤ a and ω ( x ) = 0otherwise. The parameter a is chosen to be 0.3, a value r g AA e ff (r) k=0k=28k=125k=1250 r g AA e ff (r) / y T=0.82
FIG. 5: The effective A-A particle radial distributionfunction for k = 0 , , k the multi-peak structure disappears. Inset: g eff AA /y ( k ) vs r where y ( k ) is the height of the main peak. Thesmoothing of the undulation with increasing k is clearlyseen.that is slightly larger than the size of the cage (deter-mined from the height of the plateau in the MSD at in-termediate times .) Thus the quantity Q ( t ) tells whetheror not at time t a tagged particle is still inside the cageit occupied at t = 0. -1 t Q ( t ) k=0k=4k=12k=28 T=0.9
FIG. 6: Time dependence of the self part of overlapfunction Q ( t ) for systems with different values of k at T = 0 .
9. With increasing k the relaxation dynamicsquickly slows down.In Fig. 6 we show the time dependence of Q ( t ) for dif-ferent values of k . The temperature is T = 0 . k = 0 to a T that is around the onsettemperature . The graph demonstrates that with in-creasing k , the relaxation dynamics slows down quickly,in that the correlator for k = 28 decays on a time scalethat is about two orders of magnitude larger than theone for k = 0. Also note that for the largest k we clearlysee a two-step relaxation, i.e., the hallmark of glassy dy-namics in which the particles are temporally trapped by their neighbors , while for k = 0 one has just a sim-ple one-step relaxation, i.e., a normal liquid state re-laxation. These results demonstrate that the presenceof the pseudo neighbors does have the sought after ef-fect of strongly slowing down the relaxation dynamics ofthe system, although, as demonstrated above, the overallstructure of the liquid is not changed. Interestingly theshape of the time correlation function in the α -relaxationregime does not seem to have a noticeable dependenceon k , indicating that the relaxation mechanism is weaklydependent on k . However, this conclusion only holds forlength scales on the order of ′ a ′ while it could be that onlarger scales differences become noticeable. Here we alsonote that for other mean-field like models, such as theone introduced by Mari and Kurchan , an increase ofthe interaction range leads to an acceleration of the dy-namics, i.e. the hoped for slowing down of the dynamicsis not necessarily guaranteed.Next, we compare the time dependence of the meansquared displacement, averaged over all the particles, oftwo systems, k = 0 and k = 28, Fig. 7. For the k = 0system we show the MSD for T = 0 .
82, i.e., a tempera-ture close to the onset T and as a consequence one seesthat the curve shows between the ballistic regime at shorttimes, ∆ r ( t ) ∝ t , and the diffusive regime at long times,∆ r ( t ) ∝ t , a weak shoulder. Qualitatively the sametime-dependence is found for the k = 28 system, but thistime at the higher temperature, T = 1 .
5, indicating thatthe increase of k leads to an increase of the onset tem-perature. If for the k = 0 system the temperature islowered to 0.445, the MSD shows at intermediate timesa very pronounced plateau that is due to the temporarycaging of the particles . The same behavior is found in -3 -2 -1 t -5 -4 -3 -2 -1 M S D k=0, T=0.445k=28, T=0.82k=0, T=0.82k=28, T=1.5 FIG. 7: Time dependence of the mean squareddisplacement for the k = 0 and k = 28 systems in thehigh and low temperature regimes. The curves are forsimilar value of relaxation time. The k = 28 systemshows a weak sub-diffusive behaviour at high and lowtemperature. -2 -1 t d l og ( M S D A ) / d l og ( t ) T=5.0T=1.0T=0.8T=0.445k=0(a) 10 -2 -1 t d l og ( M S D A ) / d l og ( t ) T=5.0T=2.0T=1.0T=0.65k=12(b)10 -2 -1 t d l og ( M S D A ) / d l og ( t ) T=10.0T=5.0T=2.0T=1.0T=0.82k=28(c) t -2 -1 M S D A t w =0t w =10000t w =20000t w =30000t w =40000 (d) k=28 T=1.0 FIG. 8: Double-logarithmic derivative of the MSD of the A particles as a function of time. (a) System for k = 0. Iftemperature is decreased the derivative shows at low T a local minimum, indicating the presence of caging. (b)System for k = 12. Qualitatively the same time dependence as in panel (a) but now at higher temperatures. (c)System for k = 28. One sees that the curves show at intermediate times a plateau that is due to the caging causedby the pseudo neighbors. The arrows pointing upward [downward] in panels (a)-(c) indicate τ [ τ ], the location ofthe peak in the non-Gaussian parameter α ( t ) [in the dynamic susceptibility χ ( t )]. (d) MSD of the A particles fordifferent waiting times t w (see legend). No waiting time dependence is noticeable.the k = 28 system at T = 0 .
82 with a plateau height andlength that is very close to the one of the k = 0 system.(This similarity is due to our choice of the temperature T = 0 . k , see Fig. 1, the pronounced caging for the k = 28system (at T=0.82) is thus due to the pseudo neighbors,i.e., the non-local interactions. From these curves wehence can conclude that the presence of the additionalinteractions leads to a substantial slowing down of therelaxation dynamics while the details of the MSD, suchas the height of the plateau or its width, at the same ef-fective temperature (discussed below) are modified onlymildly, at least in the parameter regime probed here.At sufficiently long times the motion of the particlesis expected to be diffusive, and hence the MSD shouldincrease linearly in time. Fig. 7 shows that for the k = 0system, this is indeed the case and that this diffusion setsin once the MSD has reached a value around 1.0. Inter-estingly one observes for the k = 28 system even at thelongest times a sub-diffusive behavior, with an exponentthat is around 0.8, and this even for values of the MSDthat are on the order of 10. This behavior can be noticedbetter by calculating the slope of the MSD in the log-logpresentation, see Fig. 8. For k = 0, panel (a), we see that E I S ( T ) k=0k=4k=12k=280 0.5 1 1.5 2 2.5 3 3.5 4 T -0.6-0.4-0.20 E I S ( T )- E I S ( T = . ) (a)(b) FIG. 9: (a) Inherent structure energy, E IS , as a functionof temperature for the k = 0 , , , and 28 systems. (b)Shifted (by E IS ( T = 4 . T . Near T onset the energy starts to deviate from its hightemperature value allowing to determine T onset . Withincreasing k , T onset moves to higher temperatures.at short times the slope is 2.0, as expected for a ballisticmotion. At high temperatures the slope crosses over to1.0 at around t = 3, i.e. the system becomes diffusive.If T is lowered, the slope starts to show a dip with adepth and width that increase rapidly with decreasingtemperature. For long times we see, however, that thecurves again attain the value of 1.0, i.e. the system isdiffusive. Qualitatively the same behavior is found for k = 4 (not shown) and k = 12, panel (b). However, acloser inspection of the curve for T = 2 . ≤ t ≤ k = 28, panel (c),but now the mentioned plateau at intermediate times be-comes more visible since its height has decreased to 0.8,i.e. the deviation from the diffusive regime become morepronounced. We now clearly see that if the tempera-ture is lowered the curves reach this second plateau ata later time, but its height is unchanged (see the curvesfor T = 1 . r c , and, because of geomet-rical reasons (the volume of the spherical cap increaseswith L ij ) and the fact that L ij > r c , this takes certainlymore time than cutting just the interactions between thetagged particle and its nearest neighbors, which explainsthe long time tail in the MSD. Note, however, that forsufficiently long times the MSD can be expected to be-come diffusive for all values of k , see, e.g., the curve for T = 2 . . In order to distinguish in the followingthe two mentioned processes, we will refer to the one cor-responding to the particles leaving their nearest neighborcage as the “NN- α -process”, while the dynamics in whichthe pseudo-neighbors leave the interaction range of thetagged particle will be referred to as the “PN- α -process”.Note that although Fig. 8 clearly indicates that there aretwo processes, we will see in the following that not all ob-servables reveal this in a direct manner. For example, thetime dependence of Q ( t ), presented in Fig. 6, does not in-dicate an obvious presence of two different α − processes,although the pseudo-neighbors can be expected to affectnot only the relaxation time but also the details of thecorrelator.Since the onset temperature is an important point onthe energy scale of the system, we now have a closer lookat the k -dependence of T onset . As mentioned above, thistemperature can be identified from the first occurrenceof a plateau in the MSD. Alternatively one can study τ k=0k=4k=12k=28(a) T g /T τ k=0k=4k=12k=28(a) T g at τ =1000 FIG. 10: (a) Arrhenius plot of the α -relaxation time forsystems with different values of k . The lines are fitswith the Vogel-Fulcher-Tammann expression, Eq. (22).(b) Same data as in (a) but now as a function of thescaled temperature T g /T , with τ ( T g ) = 10 .the inherent structure energy, E IS , which shows at T onset a marked change in its T -dependence . (We recallthat E IS of a configuration is the potential energy eval-uated at the local minimum of the energy reached fromthe configuration via the steepest descent procedure.) InFig. 9(a) we show E IS as a function of T , with the differ-ent curves corresponding to different values of k . Fromthe graph, one recognizes that with increasing k the en-ergy decreases, an effect that is due to the presence of thepseudo neighbors which can lower the energy by occupy-ing the well in the interaction potential. Less trivial is thefact that the temperature at which the curve starts to de-crease rapidly, i.e. the onset temperature, increases withincreasing k . Thus the increase of T onset with k can beseen directly from this static observable. In order to seebetter the k -dependence of T onset , we plot in Fig. 9(b) theinherent structure energy shifted by E IS ( T = 4 . T = 4 . k , demonstrating the increase of the onsettemperature. Fitting two straight lines to the data for T > T onset and
T < T onset , their intersection point canbe used to determine T onset . As we will show elsewhere ,the so obtained values are compatible with the values ofonset temperature as determined from the entropy . InTable I we list the values of T onset obtained from thesecurves and one sees that for k = 28 this temperature isabout 90% higher than T onset for k = 0.A further important quantity to characterize the re-laxation dynamics of a glass-former is the α -relaxationtime τ . Here we define this time scale via Q ( τ ) = 1 /e .This definition is reasonable since we have seen in Fig. 6that the shape of the time correlation functions is basi-cally independent of k . (Note that with this definition of τ we do not distinguish between the NN- α -process andthe PN- α -process discussed in the context of Fig. 6. Forthe values of k considered here, this is justified since thefinal decay of Q ( t ) involves both processes.) Fig. 10(a)is an Arrhenius plot of τ for the different systems. Oneclearly sees that with increasing k , the dynamics quicklyslows down and that the bending of the curve seems toincrease, i.e. the system becomes more fragile. To quan-tify this trend as a function of k , we have fitted τ ( T, k )at intermediate and low temperatures to a Vogel-Fulcher-Tammann(VFT)-law: τ ( T ) = τ exp h K ( T /T − i . (22)Here T is the so-called VFT temperature at whichthe relaxation time of the system is predicted to diverge.The parameter K describes the curvature of the datain an Arrhenius plot and hence can be considered as ameasure for the fragility of the glass-former. The figuredemonstrates that this functional form gives a good fitto the data (solid lines) and hence allows to estimate T and K .The values of T are included in Tab. I as well andone sees that T changes by about a factor of two if k isincreased from 0 to 28, i.e. a factor that is comparableto the one found for T onset . In contrast to this we findthat the parameter K occurring in the Vogel-Fulcher-Tammann-law, Eq. (22), increases by about 30% in theconsidered k -range, see Tab. I. This indicates that theintroduction of the pseudo neighbors renders the systemincreasingly more fragile. Another way to see this is todefine an effective glass transition temperature T g via τ ( T g ) = 10 and to plot the relaxation time as a functionof T g /T . This is done in Fig. 10(b) and one sees thatthe curves for large k are indeed more bent than the onesfor small k , i.e. the fragility of the system increases with k . This trend is thus qualitatively similar to the obser-vation of Ref. 29 in which it was found that increasingthe dimensionality of a glass-former gives rise to a higherfragility. D. MCT power law
Having presented our findings regarding the relaxationdynamics of the system we now probe whether this dy-namics can be described by means of mode coupling the-ory. MCT predicts that close to the critical temperature T c of the theory the relaxation times show a power lawdivergence: τ ( T ) = τ MCT ( T − T c ) − γ . (23)Using this functional form to fit the temperature de-pendence of the relaxation time we obtain T c ( k ) (valuesare given in Tab. I). In Fig. 11(a) we present a log-logplot of the relaxation time as a function of the normalizedtemperature ( T − T c ) /T c . One recognizes that for k = 0,the increase of τ with decreasing T is described well bya power law (dashed line), in agreement with previoussimulations . However, at the lowest T ’s deviationsare observed, and the increase in τ is weaker than thepower law predicted by MCT. This deviation is usuallyattributed to the existence of “hopping processes”, i.e. acomponent in the relaxation dynamics that is not takeninto account in the idealized version of the MCT. Thetwo arrows in the plot delimit the T -range in which thepower law gives a good description to the data.For the system with k = 28 the temperature depen-dence of τ is qualitatively very similar to the one for the k = 0 system, if one plots the data as a function of thereduced temperature ( T − T c ) /T c . The highest temper-ature at which the data follows the power law (dashedline), marked by an arrow, is around 2 T c , and very closeto the corresponding reduced temperature for the k = 0system. However, the lower (reduced) temperature atwhich τ starts to deviate from this power law, see arrow,is smaller for the k = 28 system than the corresponding T for the k = 0 system, showing that for the former sys-tem the mentioned hopping processes are less important,i.e., the system is more mean-field like. For the k = 28system, this lower limit is about a factor of 3 smallerthan the limit for k = 0; thus the T -range in which theidealized MCT can be expected to be reliable has in-creased significantly by the introduction of the pseudoneighbors. In Tab. I we have also included the value of T c and one recognizes that the critical temperature for k = 28 is about 90% higher than the one for k = 0,i.e. the k -dependence of T c is very similar to the one of T onset .According to the analytical calculations for the mean-field p -spin model, for which there is no activated dy-namics, the onset temperature coincides with the MCTtemperature which is also the temperature at which thedynamics diverges . (Note that this is only true inthe thermodynamic limit while for finite systems one hasvery strong finite-size effects that completely wash outthese transitions, see Ref. 43.) For the GCM it was foundthat the relative distance between the three temperatures T onset , T c , and T , is much smaller than the one we findhere for the k = 0 system . Thus the reduction ofthis relative distance with increasing k , given in Tab. I,can also be taken as a signature of increasing mean-fieldlike behaviour.From Fig. 11(a) we recognize that the relaxation timesfor the k = 28 system are shorter than the ones for the0TABLE I: The value of the characteristic temperatures and the kinetic fragility parameter for systems with differentvalues of k . T onset is the onset temperature at which the inherent structure energy starts to deviate significantlyfrom its high temperature value. T c is the MCT transition temperature. T is the singular temperature of theVogel-Fulcher-Tammann equation, Eq.(22). All characteristic temperatures increase with increasing k . Alsoincluded are the normalized differences between various temperatures. K is the kinetic fragility defined in Eq. (22). x ( k ) is the prefactor needed for the scaling plot shown in Fig. 11(b). k T onset T c T T onset − T c T c T onset − T c T onset T onset − T T T c − T T K x ( k )0 0.74 ± .
04 0.43 0.283 0.72 0.42 1.61 0.52 0.184 1.04 0 . ± .
08 0.51 0.362 0.63 0.38 1.29 0.41 0.237 1.5512 1 . ± .
07 0.62 0.465 0.66 0.40 1.22 0.33 0.286 2.028 1 . ± .
22 0.80 0.610 0.60 0.38 1.10 0.31 0.297 2.1 k = 0 system if compared at the same reduced tempera-ture. In fact, as plotted in Fig. 11(b) on an intermediatetime scale the two data sets can be superimposed withhigh accuracy by applying a multiplicative factor x ( k )(see Tab. I for values). Thus we conclude that the maindifference in the two data sets is the prefactor τ MCT inEq. (23). A decrease in τ MCT implies a faster motion in-side the cage, and this is in fact very reasonable since withincreasing k the tagged particle is interacting with moreparticles, thus making its effective cage stiffer. Anotherway to present this result is to plot the time scale τ · x ( k )as a function of T c /T , see Fig. 11(c). We find that thisrepresentation of the data gives rise to a collapse of thecurves for the different values of k , demonstrating thatthe T -dependence is indeed very similar at intermediatetemperatures. Hence we conclude that the introductionof the pseudo neighbors does not only increase the α -relaxation time strongly but also increase somewhat theattempt frequency with which the particle tries to leavethe cage. E. Dynamic Heterogeneity
One of the hallmarks of glassy dynamics is that timecorrelation functions are stretched in time. The reasonfor this non-Debye relaxation has been a long-standingpuzzle with the contrasting views that each small do-main of the sample shows the same stretched time de-pendence or, alternatively, that the stretching is relatedto dynamical heterogeneities . Experiments and sim-ulations have shown that the homogeneous scenario isnot compatible with the observations, i.e. glass-formingsystems do have a significant amount of dynamical het-erogeneities (DH) . In this final section, we thereforediscuss the k -dependence of these DH and probe whetherwith increasing k one does indeed find a decrease of thesefluctuations, the behavior expected for a mean-field sys-tem.One first step to probe the DH is to look at the so- called non-Gaussian parameter (NGP) α ( t ) which is de-fine by α ( t ) = 3 h r ( t ) i h r ( t ) i − , (24)where r ( t ) is the displacement of a tagged particle withina time t . Thus α ( t ) measures whether or not the distri-bution of the particle displacement is Gaussian .In Fig.12(a) we plot the NGP for the A particles in the k = 28 system. Interestingly one finds that at high tem-peratures α ( t ) has two peaks: A first one at t around0.6 and a second one at t ≈ α -process, in agreement with ear-lier studies . The second peak has so far not been seen inthe glass-forming systems considered before and is likelydue to the breaking of the bonds with the pseudo neigh-bors, i.e. the PN- α -relaxation. Note that the presenceof this second peak is coherent with our findings for theMSD, see Fig. 8(c), for which we observed a plateau inthe slope that, for T = 2 .
0, ended at around t = 10 and we had argued that this is due to the motion of thepseudo neighbors. If T is lowered, the first peak in α ( t )rises quickly and dominates the second peak, i.e. on over-all the time dependence of the NGP becomes again quitesimilar to the one that has been observed in previousstudies of glass-forming systems. The main difference isthat in our case the second peak will make the decay of α ( t ) slow since at long times the dynamics will be in-fluenced by the pseudo neighbors, which decorrelate onlyslowly (see the data for the MSD in Fig. 8).The influence of the pseudo neighbors on α ( t ) is shownin Fig. 12(b) where we plot this function for different val-ues of k but keeping ( T − T c ) /T c constant. One sees thatat short and intermediate times, i.e. around the peak, thecurves are independent of k , which shows that the NN- α -process is not affected by the presence of the pseudo-neighbors. Only at longer times, the curves for large k are1 -2 -1 (T-T c )/T c -1 τ T c =0.429 for k=0T c =0.798 for k=28 γ =2.4 for all k(a)10 -2 -1 (T-T c )/T c -1 τ · x T c =0.429 for k=0T c =0.798 for k=28 γ =2.4 for all k x=1.0 for k=0x=2.1 for k=28(b) T c /T τ · x k=0k=28(c) FIG. 11: (a) The relaxation time obtained from theoverlap function as a function of the scaled temperature( T − T c ) /T c for the k = 0 and the k = 28 systems. (b)Same data as in (a) but now with τ multiplied with ascaling factor x ( k ). (c) Same data as in (b) as afunction of T c ( k ) /T .higher than the ones for small k , showing that the pseudoneighbors affect the NGP only at time scales that are be-yond the time scale of the first maximum in the NGP.Since with decreasing temperature the peak correspond-ing to the NN- α -relaxation grows quicker than the secondpeak we can conclude that the dominant feature in α ( t )is due to the NN- α -process, except if k becomes muchlarger than the values we consider here.In Fig. 13 we show α p , the height of the peak in α ( t ),as a function of the reduced temperature ( T − T c ) /T c .Surprisingly we find that this quantity is completely in-dependent of k , i.e. the strength of the non-Gaussianityof the relaxation dynamics does not depend on whetheror not the system is mean-field like. In other words, the -2 -1 t -2 -1 α ( t ) T=5.0T=3.0T=1.4T=1.2T=1.0T=0.9T=0.85T=0.82 k=28(a)10 -2 -1 t -2 -1 α ( t ) k=0, T=0.75k=4, T=0.9k=12, T=1.1k=28, T=1.4(T-T c )/T c = 0.754(b) FIG. 12: (a) The time dependence of the non-Gaussianparameter, α , at different temperatures for the k = 28system. α ( t ) shows a double peak structure. (b) α ( t )at fixed reduced temperature and different values of k .The peak at short times is independent of k while theone at long times grows with increasing k . -2 -1 (T-T c )/T c -1 α k=0k=4k=12k=28 FIG. 13: The peak height of α as a function of thereduced temperature ( T − T c ) /T c for different values of k .statistics of the displacement of a tagged particle is inde-pendent of the number of pseudo neighbors, if measuredat the same reduced temperature. This result reflectsthe fact that the first peak in α ( t ) is dominated by thedynamics in which the tagged particle leaves the cageformed by its nearest neighbors.Note that α p shows a bend at around ( T − T c ) /T c ≈ . τ ·x -1 α k=0k=4k=12k=28 x=2.1 for k=28x=1 for k=0x=2.0 for k=12x=1.55 for k=4 α ~ τ ζ ζ = 0.362 FIG. 14: The peak height of α as a function of the α -relaxation time τ multiplied by x ( k ) for differentvalues of k . Also included is a fit to the data with apower law. τ τ k=0k=4k=12k=28 κ = 0.70 τ ~ τ κ FIG. 15: The time scale τ at which α ( t ) peaks, as afunction of the α -relaxation time τ . The solid line is apower law with an exponent κ = 0 . T -dependence, we expect it to be the signatureof the onset of the hopping processes mentioned above.The bend indicates that these processes start to becomeprominent at around 10% above T c , a value that seems tobe coherent with the observation from Fig. 11 regardingthe T -dependence of the relaxation times.One might wonder whether the master curve in Fig. 13is just due to the choice of the scaling factor of the tem-peratures, i.e. T c . To test this possibility, we show inFig. 14 the same data as a function of the relaxationtime τ multiplied by the same factor x ( k ) that was usedto obtain a master curve in Fig. 11(b). We recognizethat this representation leads to a very nice collapse ofthe data onto a master curve which, for intermediate andlong relaxation times, can be described well with a powerlaw with an exponent close to 0.36 (see solid line in thefigure). It is remarkable that the hopping processes dis- -2 -1 (T-T c )/T c τ · x κ k=0k=4k=12k=28 τ ·x κ ~ ((T-T c )/T c ) -1.54 FIG. 16: τ x ( k ) κ as a function of the reducedtemperature ( T − T c ) /T c . The solid line is a power lawwith exponent -1.54.cussed above, which lead to the bends in the differentcurves if the temperature approaches T c , do not seem toaffect the validity of the power law. At present, it is notclear up to which value of τ this power law will hold, inparticular, whether it will be observed at temperaturesbelow T c . Future studies on this point will certainly beof interest to understand better the relaxation dynamicsof glass-forming liquids.In Fig. 15 we plot τ , the time at which α ( t ) peaks,as a function of the α -relaxation time τ . Surprisingly wefind that the two quantities show a simple relation witheach other in the form of a power law with an exponent κ = 0 .
70 (solid line). This result can be rationalizedwithin the framework of MCT as follows: α ( t ) is relatedto the shape of the self part of the van Hove function inthat it measures its deviation from a Gaussian . At theend of the caging regime, i.e. the β -relaxation, some ofthe particles will have already left their cage, thus givingrise to a tail to the right of the main peak of the vanHove function. It is this tail that is responsible for thenon-Gaussian shape of the van Hove function and henceleads to an increase of α ( t ). Thus it is reasonable toassume that τ is directly related to the time scale of the β -relaxation τ β . MCT predicts that the latter time scaleincreases like τ β ∝ ( T − T c ) − / (2 a ) . (25)The α -relaxation time τ is instead predicted by MCTto increase like τ ∝ ( T − T c ) − / (2 a ) − / (2 b ) = ( T − T c ) − γ . (26)In Eqs. (25) and (26) the parameters a and b can inprinciple be calculated from the T -dependence of thestatic structure factor or, exploiting Eq. (26), determinedfrom the T -dependence of the relaxation time . Forthe k = 0 system it has been found that a is around 0 . b is around 0 . . Combining these last twoequations gives, under the assumption that τ ∝ τ β , τ ∝ τ b/ ( a + b ) . (27)Thus we find a power law dependence with an exponentof 0.63 (using the mentioned values of a and b ), whichis indeed very close to our exponent κ from the fit (0.7).We mention here that the observed power law extendsover the whole accessible range of τ , i.e. it also includesthe temperature regime in which we expect hopping pro-cesses to be present, in agreement with earlier resultsfor the KA model (our k = 0 system) . The presentstudy shows that this relationship between τ and τ is in-dependent of k , corroborating our earlier statements onthe nature of the NGP.To get Eq. (27) we have made the assumption that τ is proportional to τ β . As argued above, this hypothe-sis is reasonable since it can be expected that the non-Gaussian parameter peaks at a time at which a substan-tial number of particles start to leave their cage and MCTdefines τ β as the time at which the correlator starts todrop below the plateau at intermediate times . Previ-ous studies have therefore made the assumption that τ β can be determined from the minimum in the slope of theMSD . However, we argue that such an identificationmight be misleading: For the case of a system with New-tonian dynamics, the phonons that govern the short-timedynamics mask the critical decay of the time correlationfunctions thus also masking the correlation between theabove-mentioned minimum and τ β . (This effect is, how-ever, absent if the system has a Brownian dynamics .)Therefore we think it is more appropriate to determine τ β from a quantity that is not directly influenced by thesevibrational modes, such as the α ( t ) considered here. InFig. 8(a)-(c) we have also included for the various curvesthe times τ , arrows pointing upward, and one sees thatthey do not correspond to the location of the minimum inthe curves but that they are located at somewhat largertimes, as expected because of the mentioned effect of thephonons. Although at present we do not have any solidproof why τ does indeed correspond to τ β , our findingthat the relation between τ and τ given by Eq. (27) isobeyed by our data does speak in favour of this identifi-cation. More tests on this using a system with Browniandynamics would certainly be useful to clarify this pointfurther.Finally we show in Fig. 16 the time at which α ( t )peaks, τ , as a function of ( T − T c ) /T c . Since we have ar-gued in the context of Fig. 11 that the k -dependence of τ will include a factor x ( k ) that is related to the short timedynamics, and we also showed that τ ∝ τ κ (Fig. 15), weplot directly τ · x ( k ) κ , with the values of x ( k ) obtainedfrom Fig. 11 and κ from Fig. 15. We recognize that thedata for the different values of k fall nicely on a mas-ter curve which follows a power law with an exponentaround -1.54. Also this result can be understood withinthe framework of MCT since Eq. (25) predicts that the -2 -1 t -4 -3 -2 -1 χ ( t ) T=5.0T=3.0T=1.4T=1.2T=1.0T=0.87T=0.82 (a) k=28 -1 t χ ( t ) k=0, T=0.75k=4, T=0.9k=12, T=1.1k=28, T=1.4(T-T c )/T c = 0.754(b) FIG. 17: (a) The time dependence of the dynamicalsusceptibility χ ( t ) for different temperatures for the k = 28 system. χ ( t ) increases with decreasingtemperature. (b) Time dependence of χ at a fixedreduced temperature ( T − T c ) /T c for different values of k .slope should be given by − / (2 a ) which for a = 0 . Q ( t )are related to a dynamic susceptibility which indicateswhether or not the system relaxes in a cooperative man-ner, i.e. shows dynamical heterogeneities . Thusone defines χ ( t ) = 1 N (cid:2) h Q ( t ) i − h Q ( t ) i (cid:3) (28)as a measure to quantify this cooperativity. In Fig. 17(a)we show the time dependence of χ for the system with k = 28 at different temperatures. In agreement withearlier studies, , we find that χ shows a marked peakthe height of which increases with decreasing temper-ature and also its position shifts to larger times upondecreasing T , i.e. the cooperativity becomes more pro-nounced and occurs at later times. In panel (b) of thefigure we present χ for different values of k while keepingthe normalized temperature ( T − T c ) /T c constant. Thegraph demonstrates that with increasing k the height ofthe peak decreases quickly, indicating that the systemdoes indeed become more mean-field like, as expected,4and in agreement with previous simulations of mean-fieldlike models . This k -dependence is thus very differentfrom the one seen for the height of the peak in α , high-lighting the difference between the two quantities, despitetheir (apparently) similar time dependence. We also notethat with increasing k the location of the peak in χ ( t )shifts to shorter times, in qualitative agreement with thefact that, at fixed reduced temperature, the α -relaxationtime decreases somewhat, see Fig. 11(a). -2 -1 (T-T c )/T c χ k=0k=4k=12k=28 χ ~ (T-T c ) -1.2 χ ~ (T-T c ) -2 FIG. 18: Height of the peak in χ ( t ) as a function ofthe reduced temperature for different values of k . Thedashed lines are power laws with exponent -1.2 and thesolid line is a power law with an exponent -2.To probe in more detail how the height of the peakin χ ( t ), χ p , depends on T and k we show in Fig. 18this height as a function of the reduced temperature.We see immediately that this representation of the datadoes not give rise to a master curve. With increasing k ,the curves move downwards, a k -dependence that is incontrast to the one we found for α p shown in Fig. 13.Thus we conclude that with increasing k the dynamicalheterogeneities decrease, i.e. the system becomes moremean-field like. However, we point out that even in themean-field limit these heterogeneities cannot be expectedto vanish completely which shows that this aspect ofthe dynamics is a delicate feature that is highly non-trivial.From the figure, one can conclude that for reducedtemperatures higher than around 0.1 the height of thepeak shows a power law dependence on the reduced tem-perature and we find an exponent of -1.2 that is indepen-dent of k , which implies that the dependence of χ p on thenumber of pseudo neighbors is encoded in the prefactorof the power law.The presence of power laws in χ p can be rationalized bymeans of MCT. This theory predicts that the dynamicalsusceptibility in the N V T ensemble is given by χ NVT4 ( t ) = χ NVE4 ( t ) + T c V (cid:18) dQ ( t ) dT (cid:19) , (29)where c V is the specific heat at constant volume . Evaluating this expression at t = τ , thus giving theheight of the peak, χ p , one finds that the first term on theright-hand side of the equation increases like ( T − T c ) − while the second one is found to be proportional to( T − T c ) − . Hence the power law with exponent -1.2 wefind at intermediate and higher temperatures can be in-terpreted to be due to the power law from the first term,i.e. with an exponent -1.0, which is somewhat augmentedby the presence of the second term, thus giving rise toa power law with an effective exponent smaller than -1.Thus if the mentioned hopping processes would be absentone would expect that at sufficiently low temperatures,the power law crosses over to one with an exponent -2.Whether this is indeed the case will have to be tested forsystems in which one is able to suppress these hoppingprocesses, a work that is left for the future. τ ·x χ k=0k=4k=12k=28 χ ~ ( τ ·x) χ ~ ( τ ·x) χ ~ ( τ ·x) FIG. 19: The height of the peak in χ as a function of τ · x ( k ) for different values of k . The solid line is apower law fit to the data for k = 4. The two dashedlines are power laws with exponents that correspond tothe theoretical upper and lower bounds from Eq. (30).Since the representation of the data in Fig. 18 dependson the choice of T c , it is also useful to look at the k -dependence of χ p in a more direct manner. This is donein Fig. 19 where we plot this quantity as a function of the α -relaxation time τ . (Also here we use τ · x ( t ) as abscissa,in order to take into account the trivial k dependence ofthe relaxation time.) We see that the shape of the curvesfor the different k is basically independent of k , but thatthe absolute value of χ p at fixed τ · x ( k ) decreases withincreasing k . (The same conclusion is reached if one usesjust τ as the abscissa.) Hence we confirm the conclu-sion from Fig. 17(b) that the heterogeneity of the systemdecreases with increasing k . For small and intermediatevalues of τ , the data falls approximately on a straightline, and a power law fit gives an exponent 0.51 (solidline). Expressing the T -dependence on the right handside of Eq. (29) as a function of τ = ( T − T c ) − γ , seeEq. (23), we obtain for the height of the peak χ p = Aτ /γ + Bτ /γ , (30)where A and B are expressions that have only a weak5 T -dependence. Using our value γ = 2 . -1 τ -1 τ k=0k=4k=12k=28 τ = τ τ τ FIG. 20: The location of the peak in χ ( t ) as a functionof the α -relaxation time τ . The symbols are fordifferent values of k and different T , and the solid line isa power law with exponent 1.0. Inset: τ as a functionof τ showing a power law connection between the twoquantities. The straight line has a slope of 0.70.Finally, we note that for large τ we find clear deviationsof our data from the predicted power law in that thegrowth of χ p is weaker than predicted. So in this regime,we can again invoke the argument that hopping processesdecrease the cooperativity of the relaxation dynamics.Fig. 17(a) shows that the location of the peak in χ ( t ), τ , quickly moves to larger times if the temperature islowered. To determine the connection between the α -relaxation time τ and the time scale τ we plot in Fig. 20 τ as a function of τ . Also included in the graph is theline τ = τ (solid line) and one recognizes that all thedata points fall on this line with high accuracy. Hencewe can conclude that the time scale at which the systemshows maximum cooperativity is on the time scale of the α -process, which is in agreement with earlier results .Also note that this conclusion is independent of k , i.e. thestrength of the mean-field character does not play a rolefor this result. This result demonstrates that the α -relaxation process is tightly related to the presence of thedynamical heterogeneities and that hence it is useful tostudy the latter in order to understand the slowing downof the relaxation dynamics. Finally we mention that thedirect proportionality of τ to τ and the power law con-nection between τ and τ , (see Fig. 15) implies that wehave the simple connection τ ∝ τ κ , with an exponent κ given by b/ ( a + b ), see Eq. (27). That this relation worksindeed well is shown in the inset of Fig. 20. Since theexponent κ is less than unity, we see that τ is smallerthan τ , as expected . This can also be concluded fromFig. 8 where we have added in panels (a)-(c) the values of τ (downward arrows), in that one recognizes that at low T , these are indeed to the right of the arrows presenting τ . These graphs also show that, interestingly, the (loga-rithmic) slope of the MSD at t = τ is independent of T but weakly dependent on k . IV. SUMMARY AND CONCLUSION
We have introduced a simple glass-forming systemwhich allows to tune in a smooth manner its mean-fieldcharacter. This is achieved by introducing additional k “pseudo neighbors” with which a particle can interact.These additional interactions are long-ranged and hencewith increasing k , each particle becomes increasingly con-nected with the rest of the system. However, since we alsokeep the original interaction between nearest-neighborparticles, our model has the advantage of maintaining aliquid-like structure even in the mean-field limit, i.e. thenearest neighbor distances are always of the order of theparticle diameter, which is in contrast to other modelsthat allow tuning their mean-field character .We find that the structure of the system, as charac-terized by the radial distribution function or the staticstructure factor, remains unchanged with the additionof the pseudo neighbours, also this in contrast to previ-ous models. Due to the way the model is set up, it ispossible to analytically calculate all the static structuralproperties of the system from the knowledge of the k = 0system. This allows us to understand that the additionalinteractions give rise to an effective potential that in-creases with k , thus influencing the relevant temperaturescale of the system.Due to the presence of the pseudo neighbors, the relax-ation dynamics shows a very strong dependence on k inthat the onset temperature as well as the critical temper-ature of mode-coupling theory increase with increasing k .However, once the relaxation times are expressed in termsof the critical temperature of MCT one finds only a mild k − dependence, indicating that for this class of systems T c is the most relevant parameter for the dynamics, atleast in the T − range investigated here. We note thatthe range in temperature in which MCT seems to givea good description of the relaxation dynamics increasessystematically with increasing k , thus indicating that inthe mean-field limit, the theory becomes exact. This isalso confirmed by the observation that the dynamical het-erogeneities, characterized by the dynamic susceptibility χ ( t ), decrease with increasing k .It is often believed that the fragility of the glass-formeris directly related to the presence of dynamical hetero-geneities (or more precisely to the value of the stretchingparameter β in the Kohlrausch-Williams-Watts function6used to fit the time-correlation functions) . Since wefind that the fragility of the system increases with k whilethe dynamic heterogeneity decreases we conclude thatthere is no such (strict) connection between these twoquantities, although we do not want to exclude the pos-sibility that in practice there might be a certain correla-tion. This result is in qualitative agreement with thefindings in earlier studies . Sengupta et. al. have,e.g., reported that compared to a three-dimensional sys-tem, the corresponding four-dimensional system was lessheterogeneous but more fragile . This is also corrobo-rated by experimental data analyzed by Dyre, which in-dicate that there is no direct connection between fragilityand heterogeneity .The possibility to tune the mean-field character of thesystem without changing the structure also allows elu-cidating the relation between the non-Gaussian param-eter α ( t ) and χ ( t ). While previous studies have oftenconsidered both functions to be indicators for the dy-namical heterogeneities, our analysis shows that this isnot the case at all since their dependence on k is verydifferent. Therefore our work clearly shows that thesetwo observables convey information that is very differ-ent, a conclusion that is in line with previous results thatshowed that the peak in α ( t ) has a temperature depen-dence which differs from the one of χ p . Furthermore,we also recall that for the MK-model, Ref. 23, one findsthat χ p decreases with increasing mean-field character ofthe system, i.e. the same behavior as we have found here,but that also the value of α p decreases, while in our casewe find that α p is independent of k . Also in the caseof the Gaussian core model, it was found that it’s α ( t )peak is lower than the one for the Kob-Andersen model,whereas the χ peak is much higher . The authors ofthese papers justified this results by stating that α pro-vides a measure of the degree of dynamic heterogeneityand thus its peak value should be lower for more mean-field like models and χ provides a measure of the sizeof the domains and systems which have larger domainsshould have higher value of χ . Although this interpre-tation might apply to the Gaussian core model, it is notin agreement for the system studied here and hence notgeneral. This suggests that further studies are requiredto understand the exact information provided by χ and α and if these two quantities are indeed related to eachother.Finally, we also note that the decrease of χ with in-creasing k can be due to the fact that the fluctuationsin the overlap function do indeed decrease, i.e. the re-laxation dynamics of the system becomes more homoge-neous, as expected for a mean-field-like system. However,since with increasing k the characteristic temperatures ofthe system also increase, the fluctuations should decrease.So for the moment, it is not clear which one of the twomechanisms is the main cause for the decrease of χ p thatwe observe in the present work.In an earlier study involving different glass-formers ev-idence was given that the locally preferred structures (LPS) are connected to the dynamics only for systemswhich are not mean field like . The ability of the presentmodel to continuously tune the mean-field behaviourmakes it thus an ideal system to check the validity of thisobservation. Since we find that with increasing number ofpseudo neighbours the LPS remains unchanged whereasthe dynamics slows down, this suggests that with an in-crease in the mean field nature the correlation betweenthe LPS and the dynamics decreases, a result that cor-roborates the earlier findings from Ref. 64.This summary clearly indicates that the details howthe mean-field limit is approached are important andfuture studies are needed to clarify this point. Finally,we note that the approach we propose here on how themean-field character is tuned can be applied to any sys-tem. Hence it will be interesting to study whether othertypes of interaction potentials, such as the Coulombpotentials used to describe oxide glass-formers, willgive qualitatively the same behavior, or in other words,whether the approach to the mean-field limit dependson the nature of the local structure of the system. Acknowledgements
W. K. is member of the Institut universitaire de France.S. M. B thanks SERB for funding. U. K. N thanksCSIR for his fellowship. The authors thank C. Dasgupta,D. Coslovich, M. K. Nandi, and M. Sharma for discus-sions. U. K. N thanks S. Sengupta, A. Banerjee andMd. Alamgir for help with the initial setup of the sys-tem.
V. REFERENCES W. Kob and K. Binder, Glassy Materials and Disordered Solids:An Introduction to Their Statistical Mechanics (World Scientific,2011). W. G¨otze, Journal of Physics: Condensed Matter , A1 (1999). W. G¨otze, Complex Dynamics of Glass-Forming Liquids: AMode-Coupling Theory (Oxford Science Publication, 2008). L. M. C. Janssen, Frontiers in Physics , 97 (2018). S. P. Das, Statistical Physics of Liquids at Freezing and Beyond(Cambridge University Press, 2011). T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A , 3072(1987). X. Xia and P. G. Wolynes, Proceedings of the National Academyof Sciences , 2990 (2000). V. Lubchenko and P. G. Wolynes, The Journal of ChemicalPhysics , 9088 (2003). W. Kob and H. C. Andersen, Phys. Rev. E , 4626 (1995). W. Kob and H. C. Andersen, Phys. Rev. E , 4134 (1995). W. Kob, T. Gleim, and K. Binder, AIP Conference Proceedings , 68 (1999). S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes, Proceedingsof the National Academy of Sciences , 16077 (2008). M. K. Nandi and S. M. Bhattacharyya, Journal of Physics: Con-densed Matter , 064001 (2019). S. H. Chong, Phys. Rev. E , 041501 (2008). E. Flenner and G. Szamel, The Journal of Chemical Physics ,12A523 (2013). M. K. Nandi, A. Banerjee, S. Sengupta, S. Sastry, and S. M.Bhattacharyya, The Journal of Chemical Physics , 174504(2015). M. K. Nandi, A. Banerjee, C. Dasgupta, and S. M. Bhat-tacharyya, Phys. Rev. Lett. , 265502 (2017). M. K. Nandi and S. M. Bhattacharyya, arXiv:2011.02299 (2020). J. P. Garrahan and D. Chandler, Phys. Rev. Lett. , 035704(2002). Y. Jung, J. P. Garrahan, and D. Chandler, The Journal of Chem-ical Physics , 084509 (2005). A. Ikeda and K. Miyazaki, Phys. Rev. Lett. , 255704 (2010). B. Schmid and R. Schilling, Phys. Rev. E , 041502 (2010). R. Mari and J. Kurchan, The Journal of Chemical Physics ,124504 (2011). L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli,Phys. Rev. E , 031502 (2012). D. Coslovich, A. Ikeda, and K. Miyazaki, Phys. Rev. E ,042602 (2016). A. Ikeda and K. Miyazaki, The Journal of Chemical Physics ,054901 (2011). A. Ikeda and K. Miyazaki, Phys. Rev. Lett. , 015701 (2011). M. K. Nandi and S. M. Bhattacharyya, The Journal of ChemicalPhysics , 034504 (2018). S. Sengupta, S. Karmakar, C. Dasgupta, and S. Sastry, TheJournal of Chemical Physics , 12A548 (2013). P. Charbonneau, A. Ikeda, G. Parisi, and F. Zamponi, Proceed-ings of the National Academy of Sciences , 13939 (2012). J. P. Hansen and I. R. McDonald, Theory of Simple Liquids(Elsevier, Amsterdam, 1986). D. Coslovich, M. Ozawa, and W. Kob, The European PhysicalJournal E , 62 (2018). A. Banerjee, M. K. Nandi, S. Sastry, and S. Maitra Bhat-tacharyya, The Journal of Chemical Physics , 024504 (2017). P. Chaudhuri, L. Berthier, P. I. Hurtado, and W. Kob, Phys.Rev. E , 040502 (2010). S. Sastry, Nature , 164 (2001). S. Sastry, PhysChemComm , 79 (2000). U. K. Nandi, S. Sengupta, W. Kob, and S. M. Bhattacharyya,manuscript in preparation (2020). C. A. Angell, Science , 1924 (1995). W. Kob and H. C. Andersen, Phys. Rev. Lett. , 1376 (1994). A. Cavagna, I. Giardina, and G. Parisi, Journal of Physics A:Mathematical and General , 5317 (2001). T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A , 3072(1987). T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. Lett. , 2091(1987). C. Brangian, W. Kob, and K. Binder, Europhysics Letters (EPL) , 756 (2001). M. D. Ediger, Annual Review of Physical Chemistry , 99(2000). W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C.Glotzer, Phys. Rev. Lett. , 2827 (1997). W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C.Glotzer, Phys. Rev. Lett. , 2827 (1997). B. Doliwa and A. Heuer, Phys. Rev. Lett. , 4915 (1998). M. M. Hurley and P. Harrowell, The Journal of Chemical Physics , 10521 (1996). K. Kim and S. Saito, The Journal of Chemical Physics ,12A506 (2013). T. Odagaki and Y. Hiwatari, Phys. Rev. A , 1103 (1991). F. W. Starr, J. F. Douglas, and S. Sastry, The Journal of Chem-ical Physics , 12A541 (2013). G. Buchalla, U. Dersch, W. G¨otze, and L. Sj¨ogren, Journal ofPhysics C: Solid State Physics , 4239 (1988). T. Gleim, W. Kob, and K. Binder, Phys. Rev. Lett. , 4404(1998). M. Nauroth and W. Kob, Phys. Rev. E , 657 (1997). E. Flenner and G. Szamel, Phys. Rev. E , 011205 (2005). R. Das, I. Tah, and S. Karmakar, The Journal of ChemicalPhysics , 024501 (2018). L. Berthier et al. , The Journal of Chemical Physics , 184503(2007). L. Berthier et al. , The Journal of Chemical Physics , 184504(2007). C. Brangian, W. Kob, and K. Binder, Journal of Physics A:Mathematical and General , 191 (2002). R. B¨ohmer, K. L. Ngai, C. A. Angell, and D. J. Plazek, TheJournal of Chemical Physics , 4201 (1993). X. Xia and P. G. Wolynes, Phys. Rev. Lett. , 5526 (2001). K. Niss, C. Dalle-Ferrier, G. Tarjus, and C. Alba-Simionesco,Journal of Physics: Condensed Matter , 076102 (2007). J. C. Dyre, Journal of Physics: Condensed Matter , 205105(2007). G. M. Hocky, D. Coslovich, A. Ikeda, and D. R. Reichman, Phys.Rev. Lett.113