Hopping and the Stokes-Einstein relation breakdown in simple glass formers

One of the most actively debated issues in the study of the glass transition is whether a mean-field description is a reasonable starting point for understanding experimental glass formers. Although the mean-field theory of the glass transition -- like that of other statistical systems -- is exact when the spatial dimension $d\rightarrow\infty$, the evolution of systems properties with $d$ may not be smooth. Finite-dimensional effects could dramatically change what happens in physical dimensions, $d=2,3$. For standard phase transitions finite-dimensional effects are typically captured by renormalization group methods, but for glasses the corrections are much more subtle and only partially understood. Here, we investigate hopping between localized cages formed by neighboring particles in a model that allows to cleanly isolate that effect. By bringing together results from replica theory, cavity reconstruction, void percolation, and molecular dynamics, we obtain insights into how hopping induces a breakdown of the Stokes--Einstein relation and modifies the mean-field scenario in experimental systems. Although hopping is found to supersede the dynamical glass transition, it nonetheless leaves a sizable part of the critical regime untouched. By providing a constructive framework for identifying and quantifying the role of hopping, we thus take an important step towards describing dynamic facilitation in the framework of the mean-field theory of glasses.

Introduction -Glasses are amorphous materials whose rigidity emerges from the mutual caging of their constituent particles -be they atoms, molecules, colloids, grains, or cells.Although glasses are ubiquitous, the microscopic description of their formation, rheology, and other dynamical features is still far from satisfying.Developing a more complete theoretical framework would not only resolve epistemological wrangles [1], but also improve our material control and design capabilities.Yet such a research program remains fraught with challenges.Conventional paradigms based on perturbative expan-sions around the low-density, ideal gas limit (for moderately dense gases and liquids), or on harmonic expansions around an ideal lattice (for crystals) fail badly.Because dense amorphous materials interact strongly, lowdensity expansions are unreliable, while harmonic expansions lack reference equilibrium particle positions.These fundamental difficulties must somehow be surmounted in order to describe the dynamical processes at play in glass formation.
A celebrated strategy for studying phase transitions is to consider first their mean-field description, which becomes exact when the spatial dimension d of the system goes to infinity [2], before including corrections to this description.In that spirit, we open with the d → ∞ "ideal" random first-order transition (iRFOT) scenario, which, based on the analysis of simple models, brings together static- [3][4][5] and dynamics-based (mode-coupling) [6] results for glass formation (see, e.g., [7,8] for reviews) [8][9][10][11].In iRFOT, an infinitely slowly cooled simple liquid (or compressed hard sphere fluid) becomes infinitely viscous, i.e., forms a glass in which particles are completely caged, at the (critical) dynamical transition temperature T d (or packing fraction ϕ d ).Upon approaching this transition, caging makes the diffusivity D vanish as a power-law D ∼ (T − T d ) γ , and the viscosity diverge as η ∼ (T − T d ) −γ .Hence, in the critical regime one ex-pects the Stokes-Einstein relation (SER) between transport coefficients, D ∼ η −1 , to hold.In short, the d → ∞ scenario is characterized by (i) a sharp dynamical glass transition associated with perfect caging, (ii) a power-law divergence of η, and (iii) the SER being obeyed.
Part of the difficulty of clarifying the situation in finite d, where the iRFOT description is only approximate and the dynamical transition is but a crossover, lies in the shear number of different contributions one has to take into account.From a purely field-theoretic point of view, one has to include finite-dimensional corrections to critical fluctuations.A Ginzburg criterion gives d u = 8 as the upper critical dimension for the dynamical transition [17][18][19][20], and hence for d < d u critical fluctuations renormalize the power-law scaling exponents.In principle, these corrections could be captured by a perturbative d u −d expansion, and phenomenological arguments along this direction indicate that they could also induce a SER breakdown [17].A number of non-perturbative processes in 1/d must additionally be considered.(i) In the iR-FOT picture, caging is perfect, hence in the glass phase each particle is forever confined to a finite region of space delimited by its neighbors [6].However, it has been theoretically proven [21] and experimentally observed [22] that in low-dimensional systems the diffusivity is never strictly zero.Single particles can indeed hop between neighboring cages [23][24][25][26], and the free space they leave behind can facilitate the hopping of neighboring particles.Facilitation can thus result in cooperative hopping and avalanche formation [27][28][29].(ii) For some glass formers, activated crystal nucleation cannot be neglected and interferes with the dynamical arrest, leading to a glass composed of microscopic geometrically frustrated crystal domains [30].(iii) In the iRFOT scenario, the dynamical arrest is related to the emergence of a huge number of distinct metastable glass states whose lifetime is infinite.In finite dimensions, however, a complex glassglass nucleation process gives a finite lifetime to these metastable states [5,12,31].The dynamics of glassforming liquids is then profoundly affected.Including glass-glass nucleation into iRFOT leads to the complete RFOT scenario [12], in which the mean-field dynamical glass transition becomes but a crossover [12], and both VTF scaling and facilitation are recovered [32,33].
Because the treatment of these different processes has thus far been mostly qualitative, their relative importance cannot be easily evaluated.A controlled firstprinciple, quantitative treatment is for the moment lim-ited to the exact solution for d → ∞ [10,11,34,35].Its approximate extension to finite d [6,8,36] completely ignores the non-perturbative effects mentioned above.This approach therefore cannot, on its own, cleanly disentangle the various corrections.Systematic studies of glass formation as a function of d have encouragingly shown that these corrections are limited, even down to d = 3 [15,16,[37][38][39][40], provided length and time scales are not too large, as is typical of numerical simulations and experiments with colloids and grains.In particular, with increasing d the distribution of particle displacements (the self-van Hove function) loses its second peak associated with hopping [16], the critical power-law regimes lengthen [41], and the SER breakdown weakens [15,16,40], which motivates investigating corrections to iRFOT in a controlled way.
Here we develop a way to isolate the simplest of these corrections, i.e., hopping, by studying a finitedimensional mean-field model.Through the use of the cavity reconstruction methodology developed in the context of spin glass and information theory [42], we carefully describe caging using self-consistent equations that can be solved numerically.We can thus compute the cage width distribution and isolate hopping processes.Our results provide an unprecedentedly clear view of the impact of hopping on the dynamical transition and on the SER breakdown in simple glass formers.
MK Model -We consider the infinite-range variant of the hard sphere (HS)-based model proposed by Mari and Kurchan (MK) for simple structural glass formers [43][44][45] (see SI Sec.IA for details).The key feature of the MK model is that, even though each sphere has the same diameter σ, pairs of spheres interact via an additional constant shift that is randomly-selected over the full system volume.This explicit quenched disorder eliminates the possibility of a crystal state, suppresses coherent activated barrier crossing that leads to glass-glass nucleation [44], and diminishes the possibility of facilitated hopping (as we discuss below).Yet at finite densities the number of neighbors that interact with a given particle is finite and therefore finite-dimensional corrections related to hopping remain, in principle, possible.
MK liquids have a trivial structure.Even in the dense and strongly interacting regime, the pair correlation in the liquid phase is simply g 2 (r) = θ(r − σ) (where θ(x) is the Heaviside step function), because particles are randomly displaced in space.In addition, even if both particles i and k are nearby particle j they need not be close neighbors, hence all higher-order structural correlations are perfectly factorizable.Because only two-body correlations contribute, the virial series can be truncated at the second virial coefficient [44], hence the equation of state for pressure is trivially βP/ρ = 1 + B 2 ρ, where is the volume of a ddimensional ball of radius R, ρ is the number density (the packing fraction ϕ = ρV d (σ/2)), and the inverse temperature β is set to unity [43][44][45] (see SI Sec.IA).
Note that these structural features hold for the liquid phase of the MK model in all d, and for standard HS liquids in the limit d → ∞ [8,46].The MK model therefore coincides with standard HS in that limit.For a given finite d, however, MK liquids are structurally more similar to their d → ∞ counterparts than HS liquids are.One thus sidesteps having to take into account the nontrivial structure of g 2 (r), which muddles the description of standard finite-dimensional HS [8].
For the MK model, one can easily construct equilibrated liquid configurations at all ϕ, even for ϕ > ϕ d (For standard HS, by contrast, prohibitively long molecular dynamics (MD) simulations are necessary in this regime.)This dramatic speedup is accomplished by adapting the planting technique developed in the context of information theory [47] (see SI Sec.IB).It is thus possible to study MK liquids arbitrarily close to, both above and below, the dynamical glass transition at ϕ d .A systematic study of caging beyond ϕ d is also possible thanks to the cavity reconstruction formalism, a method adapted from the statistical physics of random networks [42].

Caging -
The MK model dynamics is studied by event-driven MD simulations of planted initial configurations with N = 4000 particles (see SI Sec.IB for details) [37,38].The mean square displacement (MSD) determined from time evolution of the particle positions r i (t).At short times, before any collision occurs, ballistic motion gives ∆(t) = dt 2 ; at long times, diffusive motion gives ∆(t) ∼ 2dDt.From φonset onwards, the ballistic and the diffusive regimes are separated by an intermediate caging regime where ∆(t) ≈ ∆ is approximately constant, first appearing as an inflection point and then as a full-fledged plateau (see SI Sec.IC for definition).Simply put, after a few collisions with its neighbors, a particle becomes confined to a small region of space of linear size √ ∆, from which it can only escape, and henceforward diffuse, after a very large number of collisions.
In the d → ∞ iRFOT scenario, a sharp dynamical transition occurs at ϕ d [6,8,10], beyond which complete caging results in an infinitely-long plateau and in the disappearance of the diffusive regime.In finite-dimensional systems, one can use an approximate theory based on a Gaussian assumption for the cage shape, to obtain a prediction for ϕ d and ∆ [8,45] (see SI Sec.IIA).One can also estimate φd from the simulation results by fitting the diffusivity using the mean-field critical form D ∼ (ϕ− φd ) γ , and ∆ = lim t→∞ ∆(t) beyond φd (Fig. 1a).As expected from the suppression of various finite d corrections, the critical power-law regime is much longer for the MK model than for standard finite-dimensional HS (Fig. 1b) [44].Marked qualitative discrepancies from the iRFOT predictions are nonetheless observed.(i) Numerical estimates for φd systematically deviate from the approximate Gaussian result for ϕ d (Fig. 1c), even though the two quantities grow closer with dimension.(ii) The diffusion time τ D = σ 2 /D and the structural relaxation time τ α ∝ η (see SI Sec.IC for definitions and a discus-  where τSER = τD(ϕSER) and τ0 is the microscopic time, i.e., the characteristic time for the decay of the velocity autocorrelation function [48] (see SI Sec.IC2 for details).
In order to clarify the physical origin of the above discrepancies, we first determine whether the mismatch between φd and ϕ d is due to the hypothesis made in computing the latter, i.e., that all the cages have a Gaussian shape of a fixed diameter ∆, by using the cavity reconstruction formalism to relax both assumptions [42].Above φd , we can build the equilibrated neighborhood of particle i in order to self-consistently determine the overall cage size and/or shape distribution P f (∆) (see SI Sec.IIA for details).The process involves placing Poisson-distributed neighbors j that are randomly assigned a cage size ∆ j from a prior guess of Pf (∆), with a fixed function shape f Aj (r) (a Gaussian or a ball function, for instance).Averaging over the vibrational relaxation of each neighboring particle gives the cavity field ψ(r) felt by particle i, which is the probability density of the particle being at position r (Fig 2a It follows that deviations from the d → ∞ scenario ought to be ascribed to an imperfect caging above ϕ d in finite-dimensional systems.Microscopically, these imperfections correspond to particles trapped for a finite time before escaping to another cage through a narrow passage (Fig. 3a).Because the above calculations solely consider single-cage forms, a fixed-point distribution P f (∆) can only be reached by removing these "hopping" segments of the particle trajectories (see SI Sec.IIA for details).Not only does φd then appears at higher densities, but as long as the network of connected cages percolates dynamical arrest is formally impossible.In that context, it is interesting to note that for a prior Pf (∆) = δ(∆), the first iteration of the cavity reconstruction formalism is analogous to the void (Swiss-cheese) percolation setup for a Poisson process [49].In addition, for a non-trivial distribution of cage sizes, thresholding volume exclusion maps cavity reconstruction onto void percolation for polydisperse spheres [50] (see SI Sec.IIC).This equivalence between cavity reconstruction and void percolation sheds light on the single-cage assumption.In the iRFOT description, the MSD of each particle should remain finite when ϕ > ϕ d , but by construction the MSD can only be truly bounded if (minimally) ϕ > ϕ p , the void percolation transition.
From MD simulations of the MK model, we detect the first hopping event of each particle (see SI Sec.IIIA for details).Around ϕ d , mode-coupling and hopping processes mix, but hopping quickly dominates the dynamics upon increasing ϕ.Although the hopping of a particle does not leave an empty void in the MK model, it can nonetheless unblock a channel for a neighboring particle to leave its cage and hence facilitate its hopping.Facilitation is thus present, but weaker than in standard finite-dimensional HS, especially at high densities.Weakened facilitation is notably signaled by the fact that the distribution of hopping times computed from a regular MD simulation largely coincides with the distribution obtained in the cavity procedure, where a single particle hops in an environment where neighboring particles are forbidden to do so (Fig. 3b inset).We find the cumulative distribution of hopping times over the accessible dynamical range to be well described by a power law G h (t) = (t/τ h ) 1−µ (Fig. 3b), with the characteristic hopping time τ h increasing roughly exponentially with ϕ > ϕ d and markedly increasing with d (Fig. 3d).This Arrhenius-like scaling form is consistent with a gradual and uncorrelated narrowing of the hopping channels with ϕ.Note that similar phenomenological power-law distributions have recently been reported for other glassforming systems, such as the bead-spring model for polymer chains [51].We get back to this point in the conclusion.Finite-dimensional phase diagram -A clear scenario for hopping in the MK model follows (Fig. 4).Dynamically, the system becomes increasingly sluggish upon increasing ϕ above φonset .Initially, cages are not well formed and the slowdown exhibits a power-law scaling, according to the iRFOT critical predictions.Hopping cannot be defined because cages are too loose.Upon approaching ϕ d , however, cages become much longer-lived.In this regime, iRFOT predictions give a rapidly growing τ D , but hopping processes allow particles to escape their cages and diffuse, hence providing a cutoff to the critical divergence of τ D .The critical-like behavior of the diffusivity is also pushed to denser systems, and fitting to a power-law gives φd > ϕ d .When τ D is comparable to τ h a mixed regime emerges, characterized by a SER breakdown, as we discuss below.Even beyond φd , however, the dynamics is not fully arrested.Hopping remains possible, which shows that φd has no fundamental meaning and is just a fitting parameter associated to an effective power-law divergence of τ D .In fact, the MK dynamical data are better fitted by a VTF form than by the critical power-law (Fig. 4a), although the fitting parameter ϕ 0 has no direct static interpretation because it is intermediate between ϕ d and ϕ p .The dynamics can also be understood from the organization of cages.The critical density ϕ d of iRFOT corresponds to the emergence of a connected network of cages.Typical networks for ϕ d < ϕ < ϕ p span the system volume.When ϕ > ϕ p , they become finite and the mean network volume Vnet (sum of cage volumes in the network) follows a critical scaling from standard percolation (Fig. 4b).Based on this analysis, in the absence of facilitation the dynamical arrest should take place at ϕ p [53].Note that although above ϕ p the single-particle MSD is bounded, a particle can still explore a finite number of cages.Perfect single-cage trapping can only be found at ϕ → ∞ in finite d.Hopping is then infinitely suppressed because both the width and the number of hopping channels between cages vanish.However, even if hopping interferes with caging, well above ϕ d vibrational relaxation within the cage is sufficiently quick to numerically distinguish it from hopping.This large separation of timescales enables the facile detection of hopping in MD simulations and cavity reconstruction.But upon approaching ϕ d the task becomes acutely sensitive to the arbitrary thresholding inherent to any hopping detection algorithm [22,54] (see SI Sec.IIIA for details).
As expected from the exactness of the iRFOT descrip- tion in d → ∞, φd /ϕ d → 1 with increasing d.Both γ and γ also appear to converge to the d = ∞ value (Fig. 1b) [34].Because ϕ d < ϕ p for all d, the suppression of hopping with increasing d (see Fig. 1d inset) ought to be ascribed either to the narrowing of the hopping channels or to topological changes to the cage network.Because the pressure at the dynamical transition increases only slowly with dimension (p d ∼ d), the typical channel width is expected to stay roughly constant.The topology of the cage network, however, has a larger dimensional dependence.The cage network at percolation, for instance, has a fractal dimension d f d [52], e.g, d f = 4 for d ≥ d u = 6.Although this result is only valid at ϕ p proper, the local network structure persists at smaller ϕ because the loss of the cage network fractality takes place through the single-point inclusion of non-percolating clusters [52].The network topology is therefore such that the hopping channels (even assuming that their cross-section remains constant) cover a vanishingly small fraction of the cage surface as d increases.
The limited number of ways out of a local cage thus entropically suppresses hopping.
SER Breakdown -With hopping events clearly identified, it becomes possible to isolate the pure critical iR-FOT (or mode-coupling) regime.Within this regime, we obtain a power-law scaling that is consistent with ϕ d (see SI Sec.IIB for details), and the SER is followed.Deviations from the extrapolated critical scaling coincide with the SER breakdown in all d.Although φonset occurs at a roughly constant distance from ϕ d , the SER breakdown occurs in systems that are increasingly sluggish with d, ϕ SER → ϕ d , and thus properly converges to the idealized mean-field behavior as d → ∞.In the MK model, the SER breakdown is thus clearly due to hopping.
By modifying the cavity reconstruction analysis, a self-consistent caging determination of ϕ d and ϕ p should also be possible for standard finite-dimensional HS.We do not attempt such a computation here, but instead use the insights gained from the MK model to associate the SER breakdown in HS with hopping.We fit the dynamical data from the regime over which the SER is obeyed to extract ϕ d and γ, and the full dynamical regime to extract φd and γ [38].As for the MK model, the two procedures converge as d increases (Fig. 5), while φonset clearly remains distant, as is observed in many other glass formers [55,56].Interestingly, for HS, ϕ SER and ϕ d are relatively close to begin with.The fairly structured pair correlation function in HS and the much larger pressure at ϕ d lead to smaller interparticle gaps.Particles are thus caged more efficiently, which suppresses hopping.Contrasting Figs.1d and 5a suggests that near ϕ SER the SER breakdown exponent ω is similar for HS and the MK model.In this regime, HS hopping is consistent with MK-like hopping.In HS, however, single-particle hopping leaves an actual structural void that enhances the correlation (and hence the facilitation) of hopping events [27][28][29].As HS become more sluggish, cooperativity plays a growing role.As a result, a pronounced difference between HS and MK hopping for ϕ ϕ SER can be observed.The lack of a notable dimensional dependence of the master curve suggests that if the SER breakdown is also affected by critical fluctuations, as suggested in Ref. [17], that effect may be hard to detect.In contrast to Ref. [16], we now understand the reduction of the measured ω as d increases to a delayed onset of hopping.

Conclusions -
We have numerically and theoretically studied a model glass former in which it is possible to isolate hopping from the critical mode-coupling dynamical slowing down, and in which no other dynamical effects are present besides these two.The results illuminate the key role played by hopping in suppressing the iRFOT dynamical transition in finite-d and in breaking the SER scaling.The MK model gives an example where single-particle hopping is sufficient to cause the SER breakdown, but in HS facilitation likely amplifies the effect, which could explain the dependence of ω on density (Fig. 5) [57].
For standard finite-dimensional HS and other structural glass formers, we expect the situation to be made more complex by the other dynamical processes mentioned in the introduction.One might then conjecture the existence of at least three dynamical regimes for glass formers, upon increasing density.(i) A iRFOT/modecoupling regime below ϕ SER .(ii) A MK-like hopping regime around ϕ SER , where hopping is the dominant correction to the iRFOT description, the mode-coupling critical scaling holds but the apparent mode-coupling transition shifts to higher densities and the effective exponent γ changes, and the SER breakdown is incipient.In this regime the hopping timescale increases (exponentially) quickly with density (Fig. 3d).We expect this increase to be similar for HS and MK liquids, because the probability of finding a neighboring cage is roughly exp(−ϕ) for both models.(iii) At yet higher densities, hopping becomes too slow and other dynamical effects likely become important.If glass-glass nucleation barriers do not grow as quickly as the hopping barriers, then these processes may eventually become the dominant relaxation mechanism, following the RFOT prediction [5,12,31].In this regime (hence in deeply supercooled liquids much below T d ) the VTF law and the associated Adam-Gibbs relation should be reasonably well obeyed.Note that other processes such as cooperative hopping dressed by elasticity might also occur in this regime [26].Note also that these different regimes are probably not separated by sharp boundaries in realistic systems, hence all these relaxation processes might coexist, making their identification quite challenging.
We would also like to stress, in line with previous studies, that VTF fits of the structural relaxation time in regimes (i) and (ii) should not be used to extract the putative Kauzmann transition point.In our opinion it makes no sense to test the Adam-Gibbs relation in these dynamical regimes.In the MK model, although the VTF law can be used to fit the dynamical data, there is indeed no associated Adam-Gibbs relation and ϕ 0 has no thermodynamic meaning.In particular, ϕ 0 is not associated with a Kauzmann transition (which in the MK model only happens at ϕ = ∞ [44]).This observation is particularly important for numerical simulations and experiments on colloids and granular systems, which are most often performed in the vicinity of ϕ d and ϕ SER , and hence are found within the first two regimes.
Finally, we note that the MK model could also serve as a test bench for descriptions of hopping [24,25,58], as well as for relating percolation and glassy physics more broadly [59].These studies may clarify other finitedimensional effects, such as the correlation observed between local structure and dynamics [30].The infinite-range variant of the Mari-Kurchan (MK) model [44] is defined by adding to the distances between pairs of particles an additional quenched random shift that spans the whole system size (Fig. 6a).The Hamiltonian for N hard spheres (HS) is thus where U (r) for |r| = r is the HS potential (e −βU (r) = θ(r − σ)), for spheres of diameter σ, and Λ ij is a uniformly distributed vector within the system volume V , i.e., with a probability distribution P (Λ ij ) = 1/V .Note that the standard HS model corresponds to Λ ij = 0. Note that even if in principle all particles interact with all others, in practice because U (r) is short-ranged, a given particle only interacts directly with a finite number of neighbors (the first coordination shell), as in usual liquids.Hence, the model is akin to a mean-field spin-glass model with finite connectivity with the connectivity depending on the number of neighbors in the first coordination shell, and thus on both density and dimension.
MK liquids have a simple structure in all spatial dimension d, because random shifts eliminate higher-order correlations.For example, consider two particles j and k both near particle i, i.e., Unlike in regular HS, in the MK model particles j and k have a negligible probability of being near each other (|r j − r k + Λ jk | > σ), because their effective distance is shifted by Λ jk , which is of the order of the system size.This argument can also be generalized to interactions between more particles.Particle i thus has hard-core interactions with its neighbors, but with probability one in the thermodynamic limit these neighbors can overlap with each other.The two-point correlation function seen from one particle is simply where ) be the volume of a d-dimensional ball of diameter σ, and . For the MK model, the virial expansion of the equation of state (EOS) terminates at the second-order  Compared to HS, the MK model has several unique features.(i) Although monodisperse HS easily crystallize in low dimensions, the random shifts in the MK model impose a quenched disorder that is incompatible with crystal symmetry and fully suppresses the crystal phase.(ii) Glass-glass nucleation [5,31] is also suppressed, because if a nucleus forms around a particle, the particles inside this nucleus are actually randomly distributed in real space and thus no surface can be formed (Fig. 6).The free energy cost of forming a nucleus hence scales with the system size and diverges in the thermodynamic limit.(iii) Particle hopping is much less correlated.By contrast to HS, where the hopping of a particle increases the chance that one of its neighbors also hops due to the real-space void it leaves behind, facilitation is limited to unblocking an escape channel in the MK model (Fig. 6b).(iv) MK particles are distinguishable because the quenched shifts {Λ ij } are fixed.Besides the lack of structure, the partition function Z of the MK model is therefore different from that of HS by a factor of N !, i.e., Z MK /N !∼ Z HS , and hence S MK liq ∼ S HS liq + ln N .As a result, the density of the Kauzmann transition in MK diverges in the thermodynamic limit [44].
Introducing a set of quenched random shifts brings two key advantages from a methodological point of view.First, in computer simulations, it is convenient to "plant" an equilibrated MK configuration (Sec.A 2 b).Planting avoids the (circular) difficulty encountered in most other glass-forming liquids of equilibrating an initial liquid configuration before studying its equilibrium relaxation dynamics.Second, one can map the model onto a constraint satisfaction problem defined on a random graph (Bethe lattice).It is therefore possible to study its properties with the cavity method (Sec.B 1 b), which is, in principle, exactly solvable.

Simulation details a. Molecular dynamics simulations
We adapt the event-driven molecular dynamics (MD) algorithm of Refs.[37,38,41] for HS to simulate the MK model in dimensions d = 2 − 6 with N = 4000 particles.Periodic boundary conditions with the minimum image convention are implemented on the shifted distances |r i −r j +Λ ij |.For each ϕ, we perform 8 independent realizations, each corresponding to a different set {Λ ij } for a planted initial configuration (see Sec.A 2 b).Simulations are run at constant unit β, for a time t (given in units of βmσ 2 , where the particle mass m is also set to unity) sufficiently long to reach either the diffusive regime in the liquid or the asymptotic plateau in the glass.As described in Refs.[16,41], HS data is obtained from simulations of N = 8000 identical particles in d = 4 − 8, and, in order to prevent the system from crystallizing [60], from a HS binary mixture with diameter ratio σ 2 /σ 1 = 1.2 in d = 3 [16] .

b. Planting
Planting, which here consists of switching the order of determining initial particle positions {r i } and constraints {Λ ij }, is an expedient technique for studying equilibrium ensembles in random constraint satisfaction problems [47].In general, the planted ensemble is different from the annealed ensemble, but for the liquid phase it can be shown that both are equivalent, as we detail below.
In the following, we will be interested in physical observables F that depend on some initial condition {r i } and on their time evolution under deterministic MD dynamics, e.g. the mean square displacement defined in Eq. (A7).In the presence of disorder, the average of physical observables should be measured by the so-called "quenched" average where F and H depend on both {r i } and {Λ ij }.In fact, because the disorder is independent of time for a given sample, one should first perform the thermal ensemble average F for a given realization of disorder, and then repeat this operation for many extractions of {Λ ij } to average over the disorder.In simulations, however, once {Λ ij } is fixed, equilibrating independent configurations at large ϕ is very time consuming, because one should first anneal the system quasi-statically slowly up to the desired density.
Let us define the so-called "annealed" average: This average corresponds to a very different situation, where the averages over {r i } and {Λ ij } are interchangeable.Physically, this describes a situation where both variables and disorder fluctuate together; their timescales are indistinguishable.Mathematically, the last equality in Eq. (A5) shows that the integration measure can be obtained by first extracting a uniformly random configuration {r i }, and next extracing a configuration {Λ ij } from the distribution Because P ({Λ ij }|{r i }) is factorized, each Λ ij must be extracted independently, uniformly in the volume V with the constraint that In summary, we use the following procedure to compute F a : Procedure-Planting-MK 1. Generate N particle positions {r i } according to a Poisson (ideal gas) process.
2. For each pair of particles i and j, randomly and independently draw a vector Λ ij , uniformly in the sub-region of the whole volume V that is compatible with r i and r j , 3. Starting from the state given by {r i }, and for the given {Λ ij }, compute the time evolution {r i (t)} from MD simulations.From this trajectory compute F.
The key to the success of this approach is determining if, and under what conditions, the quenched and the annealed averages over the disorder are the same, F Λ = F a .Equations (A4) and (A5) coincide if the equality log Z Λ = log Z Λ holds, where Z = N i=1 dr i e −βH is the partition function for given {Λ ij } [47,61,62].This situation arises if the fluctuations of Z induced by the fluctuations of quenched disorder {Λ ij } are very weak in the thermodynamic limit.This condition is satisfied in the liquid phase, but is violated in the glass phase away from the equilibrium liquid line [47,61,62].
According to the analysis of Ref. [47], in order to check numerically that the annealed and the quenched average coincide, one should compute the vibrational (internal) entropy of the planted glass state.This can be done for example using the procedure described in Ref. [63].If the internal entropy of the glass turns out to be larger than the liquid entropy given by Eq. (A3), then the annealing average does not coincide with the quenched average [47].Fortunately, in the MK model the liquid entropy per particle diverges proportionally to log N (see Eq. (A3)), while the glass entropy per particle is finite, because particles cannot exchange (at least if one neglects hopping, as discussed below).Therefore, for N → ∞ the liquid entropy is always larger than the glass entropy, and the annealed average is correct.In other words, because the Kauzmann transition for the MK model is located at infinite density [44], the procedure is valid.
Numerical simulations show that the annealed average done using the planting procedure discussed above is in perfect agreement with the liquid EOS Eq. (A3) (Fig. 7).This result is not a surprise, because it can easily be shown that the annealed equation leads to the same liquid EOS in Eq. (A3), but it is a consistency test for the numerical procedure.Note that the pressure remains stable over time, as it should be if one initializes the MD simulation in an equilibrium configuration.

Basic phenomenology of glassy behavior and definitions of physical quantities
Before turning to a more detailed explanation of our results, in this section we summarize the main physical observables that we investigate in this study, with a short account of their definition and of the main results.

a. Mean square displacement (MSD) and cage sizes
Despite its trivial liquid phase, the MK model presents a complex glass-forming and glassy behavior.Above the onset φonset of sluggish dynamics (see Sec.A 3 c for definition), We can distinguish three main regimes in the mean square displacement: (i) a ballistic regime with ∆(t) = dt 2 ; (ii) a caging regime with a plateau ∆(t) ∼ ∆, where ∆ is the mean cage size; and (iii) a diffusive regime with ∆(t) = 2dDt, where D is the diffusivity.According to mode-coupling theory (MCT), the plateau becomes asymptotically stable beyond the dynamical transition ϕ d (see Sec. B 2 a).We can then formally define the mean cage size as the infinite time limit of the MSD and the individual cage size ∆ i of each particle i From this definition and the equilibrium conditions r i (0) = r i (t) and |r i (0)| 2 = |r i (t)| 2 , we obtain another expression for ∆ i The definition of ∆ i in Eq. (A9) can be directly used to measure individual cage sizes in numerical simulations.Equation (A10) also suggests that ∆ i is twice the variance of the distribution of particle positions within a cage.In theoretical calculations, a cage form ansatz f A (r) is usually used for this distribution.Two commonly used functions are the Gaussian and the ball functions Below, we use the Gaussian anzatz in the replica method (Sec.B 1 b), and both ansatz in the cavity method (Sec.B 1 b).The parameter A in these functions can be related to ∆ i using Eq.(A10) for the Gaussian function and for the ball function.

b. Characteristic timescales
In this subsection, we define the characteristic timescales, their physical interpretations, and how they are numerically determined.
• τ 0 -microscopic time.This natural timescale serves as reference to compare the evolution of other timescales with spatial dimension d.Its definition is such that the velocity autocorrelation function d(τ 0 ) = 1/e, where (see Fig. 8c) [48].
• τ D -diffusion time.The characteristic time for diffusion time is defined as τ D = σ 2 /D, such that ∆(t) vs t/τ D collapses in the caging and diffusive regimes (Fig. 8a), as predicted by MCT (see Eq. (B16)).Using this collapse, we can determine τ D without explicitly extracting D, which allows us to estimate τ D close to the dynamical transition, even when the fully diffusive regime itself is beyond numerical reach.
• τ α -structural relaxation time.In standard glass-forming liquids, τ α is typically extracted from the decay of the self-intermediate scattering function where k * is the first particle peak of the structure factor For the MK model, however, this method cannot be directly applied because the trivial structure of g 2 (r) (and hence of S(k)) leaves k * ill defined.Here, we use a slightly different, although consistent, approach to measuring τ α .We first generalize the definition of the MSD to the typical displacement of particles which is the zeroth moment of the self van Hove function G s (r, t), i.e., the displacement of the majority of particles at time t.In practice, to determine r typ (t), we use z = 0.1 which is very close to the limit z → 0. By analogy to τ D , we then determine the relaxation time τ α by ensuring that r 2 typ (t) vs t/τ α collapses in the MSD caging regime (Fig. 8b).Note that τ α is then only defined up to an overall constant that is independent of density.
For HS, this (re)definition of τ α is consistent with the traditional one, because the condition in Eq. ( A16) is equivalent to kr typ (τ α ) ∼ 1.The length scale 1/k * indeed corresponds to that of the maximum density fluctuation, which should be of the order of the typical cage diameter.The scaling r typ (τ α ) ∼ 1/k * ∼ √ Ā shows that r typ is near the caging regime at τ α , and hence should be independent of density.Our estimate of τ α is therefore consistent with the proportionality relation for the viscosity τ α ∼ η observed in very sluggish fluids [16].Note that the above definitions of τ α and τ D give additional weight to slower and faster particles, respectively.In this context, the breakdown of SER is consistent with a proportion of fast particles that is larger than expected [14].
• τ h -hopping time.The typical time for a caged particle to escape (see Sec. C for more details).

c. Characteristic densities
In this subsection, we define the characteristic densities (number density and volume fraction are used interchangeably), their physical interpretations, and how they are numerically and theoretically determined.Results for HS and the MK model are reported in Table I.
• φonset -onset density of the glassy behavior.It corresponds to the lower limit of the caging regime [64].Its choice is such that for ϕ < φonset no inflection point appears in the logarithmic-scale MSD; for ϕ ≥ φonset , the MSD shows an inflection point, i.e., a point where d 2 ln ∆(t) (d ln t) 2 = 0.In this regime a non-Fickian behavior is observed.Hence, φonset also corresponds to the density at which the minimum value of the non-Fickian coefficient ν min is unity (Fig. 8d), where ν min = min t ν(t) and the non-Fickian coefficient ν(t) ≡ d ln ∆(t) d ln t [65].Note, however, that our estimate of φonset likely underestimates the onset calculated from the emergence of a finite configurational entropy in static calculations [66].
• ϕ SER -characteristic density for the breakdown of the Stokes-Einstein relation (SER).Below ϕ SER , hopping is irrelevant because, if present, it is indistinguishable from the regular liquid dynamics, and the MCT scaling relations are satisfied (Sec.B 2 a); above ϕ SER , hopping becomes faster than the characteristic MCT time (τ h < τ D ), and consequently both the MCT scaling and the SER are violated.
• ϕ d -dynamical glass transition threshold.In the MK model, this density is theoretically calculated from the replica method (Sec.B 1 c), and numerically confirmed by testing the MCT scaling ) in the density range over which hopping is negligible ( φonset < ϕ < ϕ SER ).In the HS model, however, we lack reliable theoretical predictions for ϕ d in low dimensions.We therefore determine ϕ d from fitting the simulation results for D in the regime φonset < ϕ < ϕ SER .Our results are consistent with those reported in Ref. [38], where ϕ d was extrapolated from slowly quenching the fluid.Note that ϕ d is only sharply defined when hopping contributions can be separated without ambiguity.Hence, ϕ d is only well defined in the replica calculation, where hopping is excluded by construction.
• φd -effective dynamical glass transition threshold.Empirically, φd is determined by fitting the diffusivity data, as is commonly done in glass formers.In this study we show, however, that φd is systematically shifted with respect to ϕ d ( φd > ϕ d ).Note that because φd is a fitting parameter, it also depends on the density range one chooses (or is available) for the power-law fit.
• φ0 -phenomenological parameter from the Vogel-Tammann-Fulcher (VTF) fit to τ D ∼ e BVTF/( φ0−ϕ) .In the MK model, φ0 is clearly different from the thermodynamic Kauzmann transition point ϕ K = ∞.Recall, however, that the MK model lacks the glass-glass nucleation processes assumed by the Adam-Gibbs (AG) and the random first-order transition (RFOT) theories, in order to associate the divergence of the relaxation timescale with the thermodynamic singularity at ϕ K .
• ϕ p -percolation threshold for the cage network.Below ϕ p , a particle can diffuse by successive hops on the percolating network of cages.Because the infinite time limit of the MSD is only truly bounded above this threshold, ϕ p also provides a upper bound for ϕ d , i.e., ϕ p > ϕ d .
• ϕ K -Kauzmann transition.Density at which the complexity Σ (or configurational entropy) vanishes.As discussed above, because Σ MK ∼ Σ HS + ln N , the density of the Kauzmann transition diverges (ϕ K = ∞) in the thermodynamic limit.References [8,67] used the replica approach to obtain HS results, and it is straightforward to check that these derivations only rely on the pair correlation function in the liquid phase; terms corresponding to third-and higherorder structural correlations are neglected.The treatment of Ref. [8] can therefore be directly applied to the MK model, for which these assumptions are exact.The results from Ref. [67] have also been obtained using the approximation g 2 (r) = y HS liq (ϕ)θ(r − σ), see [67, Eq.( 21)] (the soft-sphere temperature is set to zero to study hard-core systems).Comparing this result with Eq. (A2), we see that for the MK model y HS liq (ϕ) = 1.All the results of Refs.[8,67] can thus be straightforwardly extended to the MK model by setting y HS liq (ϕ) = 1.(Note that the discussion of Ref. [67] was restricted to d = 3, but a general discussion for all d can be found in Ref. [8]).The replicated entropy thus has the form where f G A (r) is the d-dimensional Gaussian cage given in Eq. (A11) and I n (x) is the modified Bessel function.The last expression for q A (r) is obtained using bipolar coordinates to compute the convolution [8].Remarkably, in odd dimensions the integral over u can be computed analytically, which facilitates the numerical evaluation of the replicated entropy.
From Eq. (B1), we can derive the equation for A from the condition ∂S/∂A = 0, which reads The cage radius in the liquid can be obtained by solving this equation in the limit m → 1, and ∆ = 2dA.From Eq. (B2), one sees that the dynamical transition ϕ d corresponds to the point where 2 d ϕ d max A F (1, A) = 1.Note that for ϕ > ϕ d Eq. (B2) admits two solutions, but only the smaller of the two is a stable physical solution [8].
b. Calculation of the cage size distribution: the cavity method More information on the distribution of individual cage shapes and sizes can be obtained from the cavity method [61,68].Its application to the MK model has been developed in Ref. [45], where the cavity equations are derived and discussed.Here, we only present the main steps.
a. Cavity fields and replica symmetric cavity equations -In the cavity approach, the system is described by a set of cavity fields ψ(r).Each cavity field describes the probability of finding a particle at position r, when it is added to a system of N − 1 particles.The replica symmetric cavity equations provide a recurrence equation for determining these cavity fields In this recurrence, the new particle interacts with the N 0 other particles, each described by its own cavity field ψ j (r j ).
The interaction is given by the hard-core constraint χ(r) = e −βU (r) = θ(r − σ).In this equation the quenched random variables Λ 0j are the random shifts that appear in the Hamiltonian, but they should be independently extracted at each cavity iteration.They are independently distributed in the whole volume V with a uniform distribution P (Λ 0j ) = 1/V .Note that in Ref. [45], the cavity equations were obtained for a model defined on a random graph that is locally tree-like, corresponding to a situation where N 0 remains finite as N → ∞.The method, however, is also applicable to the MK model, where N 0 = N − 1, corresponding to the fully connected graph [61].A convenient way to obtain the fully connected graph is to first take the limit N → ∞ and then N 0 → ∞.One can show that this procedure is indeed equivalent to considering N 0 = N − 1 [61].b.Translational invariance and irrelevance of the random shifts -In order to describe the liquid and the glassy states of the MK model, we are interested in solutions of the cavity equation that have statistical translational invariance.To be more precise, the liquid phase is described by uniform fields ψ(r) = 1/V for all particles.Physically, this situation corresponds to particles diffusing everywhere within the system volume, which mathematically reproduces the virial expansion [45].In the glass phase, each individual cavity field has the form ψ(r) = f A (r − R), where f A (r) is a cage function localized around r = 0.The cavity field is thus localized around point R, but the localization centers R themselves must be uniformly distributed in the whole volume because the glass is globally translationally invariant.Hence, in Eq. (B3), when neighbors are picked at random, they are localized around uniformly distributed random positions in space, which makes the random shifts redundant.In the following, we can thus neglect the random shifts and write the replica symmetric cavity equations as and take N 0 → ∞. c.The glass phase and the 1RSB cavity equations -In the glass phase, as mentioned above, ψ(r) are random variables described by a probability distribution Q[ψ].In the regime that is here of interest, the glass is described by the 1RSB cavity equations derived in Refs.[61,68].From these equations, we obtain that the probability distribution Q[ψ] satisfies the self-consistent equation which makes explicit that the N 0 cavity fields describing the neighborhood of the new particle are extracted independently from Q[ψ].The new cavity field is constructed according to Eq. (B4) and weighted according to z m 0 .Note that z 0 is the free volume associated with the new particle, and therefore, by varying the free parameter m, one can select glassy states according to their free volume or, equivalently, their internal entropy.
d. Reconstruction equations -Beyond the dynamical transition, ergodicity is broken in the liquid phase, which corresponds to the liquid splitting into many distinct glassy states.It is well known, however, that if configurations are sampled with the equilibrium Gibbs-Boltzmann measure, then the entropy and pressure are analytic around ϕ d and the equilibrium glass phase is the analytical continuation of the liquid phase.In order to weight glassy states according to the equilibrium Gibbs-Boltzmann measure, one has to weight them proportionally to their free volume, hence one must set m = 1 [68].In the case m = 1, the 1RSB equations greatly simplify thanks to a mapping onto the reconstruction formalism [42].Reconstruction is then done by introducing new fields R r [ψ(r )] ≡ ψ(r)Q[ψ(r )].This change of variable ensures that only the fields that are localized around point Note that the reweighting term z m 0 has now disappeared from the equations.Note also that only R 0 [ψ] enters the equations and therefore all cavity fields are localized around the origin.The r j in Eq. (B6) are random shifts of the cavity fields that are constrained to be outside a sphere of radius σ around the origin.The neighbors j are thus localized outside that sphere, which guarantees that around the origin there exists a void to accommodate an additional particle.
e. Ansatz on the cage shape -As discussed in Ref. [45], numerically solving the cavity equations in Eq. (B6) remains a formidable task.Here, we make a simple ansatz on the cage shape to facilitate this computation.We first assume that the cavity fields all have the form ψ j (r) = f Aj (r − R), where f is a fixed (spherically symmetric) cage shape.We then choose either a Gaussian (Eq.(A11)) or a ball (Eq.(A12)) cage shape, with ∆ i given by Eqs.(A13) and (A14), respectively.We assume that the cage sizes are distributed according to a function P f (A) while the centers R are uniformly distributed within the volume, as discussed above.We therefore obtain the ansatz The above equations show that fields contributing to R 0 [ψ] are localized around a point R that is distributed according to f A (r), and hence R is itself localized close to the origin.Plugging this ansatz in Eq. (B6), we obtain where Note that the factor of 2 is introduced to follow the notational convention of Ref. [8].
f. Reconstruction procedure -The physical interpretation of the reconstruction equation is quite straightforward.In order to construct a new cavity around the origin, one should draw at random N 0 → ∞ particles that are located at random positions r j outside a sphere of radius σ around the origin.These particles are themselves within a cage, whose size A j is extracted from P f (A).The point r j is not the center of the cage, but a point that is typical of the distribution inside the cage.The cage center is therefore at r j +R j , where the shift R j is extracted from the cage shape f Aj (R j ).Each neighbor rattles around its cage center at r j + R j , producing an effective potential that convolutes the HS constraint with the cage shape, e −βv j eff (r0) = q Aj /2 [r 0 − (r j + R j )].The new cavity field is then given by the (normalized) exponential of the sum of all effective potentials, ψ(r 0 ) ∝ j q Aj /2 [r 0 − (r j + R j )] = exp[−β j v j eff (r 0 )].Finally, we note that although the number of neighbors N 0 should be sent to infinity, distant neighbors do not affect the new cavity field because q Aj /2 (r) tends to 1 when r → ∞.We can therefore introduce an arbitrary spatial cutoff and only consider the neighbors (whose number distribution is Poissonian) that are within this cutoff, and then increase the cutoff until the results converge.This approach is expressed by the following recursive procedure for selfconsistently determining the distribution P f (A), which is the only remaining unknown in the cavity reconstruction.Note that once P f (A) has been calculated, one can easily obtain the distribution of mean square displacements in the cage, P f (∆), according to Eq. (A13) or (A14).This observable is also easily measured in numerical simulations (and experiments).

Procedure-Reconstruction-MK
1. Consider a spherical shell σ < r < σ + σ cut of volume V 0 (the upper bound σ cut should be sufficiently large for the results to be independent of it).Consider a number of centers N 0 distributed according to a Poisson law with average N 0 = ρV 0 .Uniformly draw these sphere centers r j within the shell.
2. Independently draw N 0 cage radii A j from P f (A), and N 0 displacements R j from f Aj (R j ).
3. From these N 0 random variables, derive a new cavity field . (B10) 4. Compute the mean square displacement in the new cavity as r 0 = dr 0 r 0 ψ(r 0 ) The value ∆ new is the long-time mean square displacement corresponding to Eq. (A10).It allows one to determine the new cage parameter A new that enters in the new cage shape in Eq. (A11) (or (A12)), according to Eq. (A13) (or (A14)).
5. Repeat steps (1-4) to get N samples A new in order to construct a new distribution P f (A new ).
6. Repeat (5) until the distribution converges P f (A new ) P f (A) within the statistical error.
g. Numerical details -In principle, the above procedure provides a theoretical way to compute P f (A), but practically it must be implemented numerically, with two additional tricks.
First, we note that it is difficult to calculate the normalization of the cavity field ψ(r 0 ) in Eq. (B10), because one has to integrate over the whole space.It is more convenient to compute the variance δr 2 0 in Eq. (B11) using the Metropolis algorithm without explicitly obtaining ψ(r 0 ).We can then write Eq. (B10) as where the non-normalized probability ψ(r 0 ) is From this expression, it is clear that ψ(r 0 ) is analogous to the Boltzmann factor in the Gibbs measure with an effective potential H eff (r 0 ), ψ(r 0 ) = e −βH eff (r0) .We can thus use the standard Monte Carlo (MC) algorithm to sample any average quantity, such as δr 2 0 , with acceptance rate acc(r old 0 → r new 0 ) = min{1, ψ(r new 0 )/ ψ(r old 0 )}.(B14) Interestingly, we actually derived from the cavity formalism a "local" MC simulation.In this local MC sampling, the positions of all the particles, except for the caged particle at r 0 , are fixed and their vibrational contribution to the motion of the caged particle is integrated into the effective potential H eff (r 0 ).In our simulations, we perform 4 × 10 5 MC steps with step size 0.1 √ Ā to calculate each cage size.
Second, we have to remove hopping from the cavity procedure, or otherwise the cavity solution does not properly converge (see Fig. 9).To achieve this task, during the calculation of A new in the local MC simulations, we record the spatial trajectory of r 0 .We then check if any hopping occurs during this trajectory using the detection algorithm described in Sec.C 1. We only include A new in the statistics of P f (A new ) if no hopping is detected, as otherwise the particle is not truly caged.Once hopping is removed, our results indicate that the cavity solution properly converges when ϕ > ϕ d (see Fig. 9).We represent the distribution P f (A) by a number N = 10 4 N 0 of samples

c. Comparing theoretical predictions with simulations
We first show that the mean square displacement ∆ predicted from both the replica and the cavity methods is generally in good agreement with the simulation data (Fig. 10).The simulation ∆ is extracted from the asymptotic time limit of the MSD data, according to the MCT scaling (see Eq. (B21) below).Close to ϕ d , the replica theory predicts a scaling (Fig. 10) where ∆ d = ∆(ϕ d ), that is consistent with the MCT prediction.However, around ϕ d precisely determining ∆ from either simulations or cavity reconstruction requires a careful consideration of hopping.The simulation and the cavity data therefore unsurprisingly deviate from Eq. (B15) in that regime (Fig. 10).
Cavity reconstruction provides a theoretical prediction for the distribution P f (∆) of individual mean square displacements.In order to obtain individual cages from simulation, we use Eq.(A10) at t = 2, which is sufficiently long for the cages to form, but not so long that a large fraction of particles have hopped.Note that at densities well above ϕ d , hopping is so rare that this choice of timescale is irrelevant.As discussed in the main text, our theoretical results agree well with simulations, and are independent of the Gaussian or of the ball ansatz for the shape functional.The replica calculation for ∆ using the Gaussian functions and P f (∆) = δ(∆ − ∆) also agrees with the simulation results (Fig. 10).In summary, we find a basic consistency between our MD simulations and theoretical calculations, including (i) the replica calculation with a Gaussian anzatz for the cage shape and a δ-function approximation for the cage size distribution function P f (∆) ≈ δ(∆ − ∆), and (ii) the cavity method with both Gaussian and ball anzatzs.It has been shown that in the limit d → ∞, the theoretical result (of replica calculation) is independent of the cage shape anzats [10], and we also expect it to be independent of the method we use (replica/cavity).In finite dimensions, weak dependence is expected, but according to our results presented here, it is insignificant compared to the numerical accuracy of the resolution of the cavity equations.

Caging dynamics: mode-coupling theory (MCT) and beyond
In this section, we compare the MD results with the dynamical caging behavior predicted by MCT.The MCT scalings are found to only be consistent with our data when ϕ < ϕ SER .Above ϕ SER , MCT predictions are violated, which is well captured by the breakdown of SER and is a consequence of entangling caging with hopping, as discussed in Sec. C. It is important to note that we here only refer to MCT as the general scaling laws predicted by the schematic MCT equation [6], which can be also independently derived from the static framework [69].The traditional MCT kernel being incorrect for the MK model [44], the numerical MCT predictions are indeed unsuitable for comparison.

a. Testing the mode-coupling theory
We first compile the MCT predictions tested in our study.The derivations of these predictions as well as many important physical interpretations can be found in Ref. [6] and references therein.We denote = ϕ−ϕ d ϕ d as the distance from the dynamical transition, and τ as the characteristic time for the β-relaxation.Note that MCT does not predict any breakdown of the SER, so we do not distinguish between the α-relaxation time τ α and the diffusion time τ D in this analysis (τ D ∼ τ α τ ).Below the dynamical transition ϕ < ϕ d , MCT predicts that the time evolution of the MSD has the form where B and C are density-independent constants, and the exponents a and b are related by the exponent parameter λ as Equation (B16) shows that the relaxation process of ∆ can be divided into three regimes: (i) an early β-relaxation towards the plateau ∆ d , ∆ d − ∆ − (t) ∼ t −a , (ii) a late β-relaxation leaving from the plateau, ∆ − (t) ∼ t b , and (iii) a diffusive process that is linear in time ∆ − (t) ∼ t.We stress that, as described in Sect.A 3 a, before regime (i), there is a ballistic regime characterized by a microscopic time that is much smaller than τ and is not included in Eq. (B16).
One of the most important predictions made by MCT is that, upon approaching ϕ d , a power-law divergence should be observed for and where the exponents are related via These scalings are tested by the following procedure.

Procedure-Testing-MCT
1. Obtain ϕ d and ∆ d from the replica calculation.
2. Fit τ D according to Eq. (B18) (see Fig. 1 of main paper) with the theoretical ϕ d , in order to obtain the exponent γ.A consistent power-law scaling is only observed below some density ϕ SER ; above ϕ SER , τ D becomes smaller than the MCT predictions, implying that an additional relaxation process starts to interfere with the dynamics.This observation suggests that when we fit the diffusivity data, only the data below ϕ SER should be used.If instead we treat ϕ d as a fitting parameter for the entire density range, then we end up with shifted values φd and γ From the analysis presented in the main text, it is clear that MCT actually fails when ϕ > ϕ SER , and thus the apparent power-law fitting of Eq. (B22) is not reliable.Results for ϕ d , φd , γ and γ can be found in Table I.
3. Determine a and b (Table I) from γ using Eqs.(B17) and (B20).A good agreement is found for the entire time regime when ϕ < ϕSER (black solid lines).When ϕ > ϕSER (green dashed lines), we only observe a good agreement for the early β-relaxation regime, which suggests that at later times hopping mixes with the MCT dynamics.τ α .If SER were obeyed, we should obtain τ D ∼ τ α .As shown in Fig. 1 in the main paper, SER breaks down when ϕ > ϕ SER as where the exponent ω = 0.22 is invariant with d.
The breakdown of SER is beyond the MCT description.Interestingly, we observe that three phenomena happen at the same density ϕ SER : (i) violation of MCT scalings, (ii) violation of SER, and (iii) the hopping characteristic time τ h becoming comparable with τ D .Our interpretation of these observations is presented in the main text.

Percolation of the cage network
Because cages can be connected via hopping channels, it is natural to examine how the cages are topologically connected.We find that the network of cages spans the system below density ϕ p (ϕ p > ϕ d ).Above ϕ p , only local cage clusters are formed and particles become strictly confined.We show that this phenomenon can be mapped onto a void percolation transition, which belongs to the same universality class as regular percolation.

a. Mapping the glass transition to a void percolation transition
To do the mapping, we first consider the simplest case, where we assume that all the neighbors of a given particle are frozen, P f (A) = δ(A).We want to know if the caged particle can move to another cage without overlapping with other particles.Equivalently, we can rescale the size of neighbors as σ → 2σ, and look for a hopping path for the point representing the caged particle in the leftover void space (see Fig. 12a).
We next consider the situation where cage sizes are not zero.In this case, particle j is rattling inside a cage with radius A j , whose distribution is a density dependent function P f (A).If a certain channel were closed in the first case, there is now a possibility for it to be open because the particles bounding that channel are now thermally moving.Because we are interested in the upper bound for percolation, i.e., the best case scenario for hopping, we rescale particle sizes as (see Fig. 12a): where A j is drawn from P f (A).If no path in void space is found by this construction, then the particle is confined.Strictly speaking, this procedure only works for cage shapes with sharp boundaries, like the ball function in Eq. (A12).For a Gaussian cage in Eq. (A11), the confined particle always has a finite (but vanishingly small) probability to hop, even if the cage is found to be closed in the percolation mapping.In the following percolation analysis, to avoid any possible confusion, we assume that all cages have ball shapes, which corresponds to assuming that this probability tail is negligible.
Algorithm for radical Voronoi tessellation of polydisperse spheres in any d-Because our systems have a range of cage sizes, i.e., they map onto spheres with a polydisperse diameters, we develop a method to produce the radical Voronoi tessellation for a given configuration.The basic idea is to map the radical Voronoi tessellation in dimension d to a Voronoi tessellation in dimension d + 1, and then use Qhull [72] to compute the Voronoi tessellation.
For the standard Voronoi tessellation, the Voronoi cell for sphere i consists of space points r that satisfy the relation for any j = i.The radical Voronoi tessellation is a generalization of this definition for unequal sized spheres: where R = σ/2 is the particle radius.
In order to map the radical Voronoi tessellation to a Voronoi tessellation, we denote R max the maximum radius, and introduce a set of points in dimension d + 1, ri = (r 1 i , r 2 i , . . ., r d i , R max − R 2 i ), where i = 1 . . .N .The first d coordinates of ri are the same as r i , and the final coordinate is a function of the sphere radius.We further introduce a set of dual points r i = (r 1 i , r 2 i , . . ., r d i , − R max − R 2 i ), as images of ri 's with respect to the last coordinates.For each pair {r i , r i }, we find the d-dimensional polygon P i that is the common Voronoi boundary between these two points (Fig. 12b).According to the definition in Eq. (B25), it is clear that any point r in P i can be written as r = (r, 0), and r satisfies which is equivalent to Because this relation is exactly the definition of the radical Voronoi cell in Eq. (B26), we have proven that the d-dimensional polygon P i is the radical Voronoi cell of sphere i in the original configuration.
Determining the percolation threshold from the scaling theory-The percolation threshold can be determined by finite-size scaling [73].Note that in the void percolation analysis, the variable of interest is the volume fraction of void space η, and not directly the volume fraction ϕ [71,74].Because our planted configuration is essentially a Poisson process of overlapping spheres, however, we have Let η p = e −ϕp be the percolation threshold in the infinite system-size limit.For a system of finite linear size L ∼ V 1/d , the average effective percolation threshold ηp (L) and its variance ∆η p (L) are linearly related [75]  to determine K p (ϕ) (Fig. 12c).We finally compute ϕ p , such that K p (ϕ p ) = 1.Our system is thus percolated at ϕ p , without the extra rescaling factor K p (ϕ) (Fig. 12c).
To check the universality of the percolation transition, we examine the scaling of the mean cluster size Vnet , where V net is the total volume of the cluster of cages connected to the planted central cage.According to percolation theory, Vnet diverges at the percolation threshold as a power-law with exponent γ p : Vnet ∼ |η(ϕ) − η p | −γp . (B33) Figure 12d shows that our results are in agreement with γ p = 1.8 [75] given by lattice percolation, in support of the two problems sharing a same universality class.
Procedure for determining ϕ p -Based on the above discussion, we summarize the procedure for determining ϕ p .
Procedure-ϕ p -determination 1.For a given density ϕ, obtain the distribution P f (A) from the cavity method.
2. Plant a configuration with linear system size L such that the central particle is compatible with all neighbors (this requirement is the same as for the cavity method, see Sec.B 1 b).
4. Find the percolation threshold K p (ϕ, L) for the configuration using a binary search.To determine if the void space is percolated: • Map the d-dimensional configuration of rescaled polydisperse spheres to a (d+1)-dimensional configuration of monodisperse spheres.• Use a Voronoi protocol to calculate the Voronoi tesselation of the (d + 1)-dimensional configuration.
• Map back the (d + 1)-dimensional Voronoi tesselation to the d-dimensional radical Voronoi tesselation.
• From the radical Voronoi tesselation, construct a network.
• Determine if the network is percolated.
4. For particle i, if δ i (t * max ) > ∆ i , hopping is detected, and the process is repeated recursively for each sub-trajectory until δ i (t * max ) < ∆ i in each sub-trajectory.
Following this procedure, we save a sequence of hopping times Note that in this study we are only interested in the time of the first hopping, which is equivalent to the time during which the particle is trapped in a cage before escaping.Because facilitation is reduced in the MK model, especially at high ϕ, we do not specifically distinguish between the first and the subsequent hopping events.The algorithm generally works well at densities ϕ > ϕ d , as shown in Fig. 13.Close to ϕ d , however, hopping is mixed with other relaxation processes, and cages are not clearly defined.Detecting hopping indeed then becomes more sensitive to the specific cutoff thresholds.The hopping between cages are visualized in the time series, with the two detected hopping times (dotted lines) at t = 3948.0and t = 7863.6.

Hopping dynamics
Empirically, we find that the above detected hopping time t follows a power-law distribution (see Fig. 3 of the main paper) p h (t) ∼ t −µ , (C2) with exponent µ < 1.We can write its cumulative distribution function as where τ h is the characteristic hopping time scale, representing the time needed for all particles to hop, G h (τ h ) = 1.As shown in Fig. 3 of the main paper, both µ and τ h depend on ϕ.In particular, τ h is roughly an exponential function of ϕ τ h ∼ e αϕ , (C4) which suggests that there is no diverging density for τ h .

FIG. 1 :
FIG. 1: (a) MSD of the MK model in d = 3 for ϕ =0.40, 1.00, 1.40, 1.65, 1.72, 1.78, 1.84, 1.93, 2.00, 2.20, and 2.50, from top to bottom.The onset of caging φonset (red), the theoretical dynamical transition ϕ d (blue), and its dynamical estimate φd (magenta) are highlighted.Note that at φd and beyond a steady drift of the MSD plateau can be detected.(b) Power-law scaling in d=3 of the characteristic time τD determined by fitting φd =1.93 and γ=4.95 and by using the idealized mean-field result ϕ d =1.78 and by fitting γ=3.27.(inset) Dimensional evolution of γ and γ.The dashed line indicates the d = ∞ result γ=2.33786 [34].Solid lines are guides for the eye.(c) The dimensional scaling of φd , ϕ d , and ϕSER converges as d increases, while the onset of caging at φonset remains clearly distinct.The dashed line is the replica result ϕ d = 4.8d2 −d [8, 10].Solid lines are guides for the eye.(d) Dimensional rescaling of the SER (black) and SER breakdown (red) regimes for the MK model with ω=0.22.(inset)The ratio τSER/τ0 grows exponentially with d (solid line), where τSER = τD(ϕSER) and τ0 is the microscopic time, i.e., the characteristic time for the decay of the velocity autocorrelation function[48] (see SI Sec.IC2 for details).

6 FIG. 2 :
FIG.2:(a) Illustration of a cavity reconstruction in d=2 for a perfectly caged particle at the center.Neighboring particles at their equilibrium positions (circles) provide an effective field ψ(r) that cages the trajectory of the central particle (red line).(b) Examples of P f (∆) in d=3 from the cavity reconstruction formalism for Gaussian (straight lines) and ball (dashed lines) cage shapes compared with MD results (symbols).(c) Rescaled P f (∆) superimposed with a log normal distribution (dashed line).(d) Density evolution of ∆ measured from MD simulations (points) superimposed on the theoretical predictions of Refs.[8,45] (lines).

3 FIG. 3 :
FIG. 3: (a) Illustration of a cavity reconstruction in d = 2 for a hopping particle.In this case the neighboring particles allow the central particle to hop to other cages (red line).(b) Cumulative time probability distribution of hopping events G h (t) for d = 3 systems at densities (from top to bottom) ϕ = 1.78, 1.84, 1.90, 1.97, and 2.10, along with the power-law scaling form (dashed line).(inset) Single-particle hopping from the cavity reconstruction (circles) overlays with the MD simulations at short times (ϕ = 1.90).Phenomenological scaling parameters (c) µ and (d) τ h for the probability distribution of hopping events.Solid lines are a guide for the eye for µ and exponential fits for τ h .

FIG. 4 :
FIG. 4: (a) Dynamical and (b) static phase diagrams for the MK model in d = 3.Early in the critical regime, the relaxation times scale like a power-law, but beyond ϕSER hopping causes large deviations from this scaling.An effective φd is numerically detected instead.A VTF scaling fits the data even better.Statically, cages can be detected from ϕ d onwards by removing hopping.In reality, the fine inter-cage channels that allow hopping, however, result in a cage network.Beyond ϕp the typical network stops percolating and the network volume scales critically Vnet ∼ (ϕ − ϕp) −1.8 with ϕp = 2.40 (dashed blue line) [49, 52].The single-cage limit is reached when Vnet ∼ ∆3/2 .

FIG. 5 :
FIG.5:(a) Dimensional rescaling of the SER (black) and SER breakdown regimes for standard finite-dimensional HS.The early deviation exponent ω is consistent with hopping in the MK model with ω = 0.22 (red line, see Fig.1d), but a growing deviation is observed as ϕ increases.(b) The dimensional scaling of HS results for φd , ϕ d , and ϕSER converge as d increases, while φonset remains distinctly smaller (compare with Fig.1c).Note that in d = 8, φd , ϕ d , and ϕSER are numerically indistinguishable.(inset) Dimensional evolution of γ and γ.Both of which are consistent with the d = ∞ result (dashed line).Solid lines are guides for the eye.
details a.Molecular dynamics simulations b.Planting 3. Basic phenomenology of glassy behavior and definitions of physical quantities a. Mean square displacement (MSD) and cage sizes b.Characteristic timescales c.Characteristic densities B. Caging 1. Thermodynamics: the caging order parameter and the dynamic transition density a. Calculation of the mean caging order parameter: the replica method b.Calculation of the cage size distribution: the cavity method c.Comparing theoretical predictions with simulations 2. Caging dynamics: mode-coupling theory (MCT) and beyond a.Testing the mode-coupling theory b.Breakdown of the Stokes-Einstein relation 3. Percolation of the cage network a. Mapping the glass transition to a void percolation transition b.Methodology for determining the void percolation threshold A3)where p = βP/ρ is the reduced pressure with β the inverse temperature and ρ = N/V the number density, ϕ = ρV d (σ/2) = ρV d 2 −d is the packing fraction (we set σ = 1), S MK liq is the liquid entropy per particle, λ is the thermal de Broglie wavelength, and B 2 = V d /2 is the HS second virial coefficient.

FIG. 6 :
FIG. 6: (a) Illustration of the MK model.Left: particles in the original space.Right: particles in the shifted space with respect to particle i (red).Although neighbors cannot overlap with particle i, they are allowed to overlap with each other because of the random shifts.(b) Comparing hopping in the HS and the MK models.Left: in HS, removal (hopping) of a neighbor (blue) creates an open channel for the caged particle (red) to hop.Right: In MK, removing a neighbor is much less likely to open a channel because the other neighbors are allowed to overlap.

5 FIG. 7 :
FIG.7: Comparison between the reduced pressure p of planted states (red circles) and the liquid EOS Eq. (A3) (black solid line) in d = 3.The regime above ϕ d = 1.776 is numerically inaccessible from conventional slow quenching procedures (see, for example, the green line obtained from the Lubachevsky-Stillinger algorithm with a growth rate γ = 3 × 10 −5[38]), because the system starts to deviate from the liquid EOS around ϕ d .(inset) Planting at ϕ = 1.78 gives the correct equilibrium pressure (red line) from t = 0.

FIG. 9 :
FIG. 9: Evolution of mean square displacement ∆ under iteration of the cavity reconstruction, at (solid lines, from bottom to top) ϕ = 2.50, 2.20, 1.95, 1.80, 1.75, 1.70, 1.60 in d = 3.The solution becomes completely unstable above ∆ d = 0.267 (dotted black line), as predicted by the replica method.If hopping is not removed, the solution diverges quickly when ϕ approaches ϕ d .See, for instance, the unfiltered results for ϕ = 1.95 (pink dashed line).

5 FIG. 10 :
FIG.10:(a)The d = 3 mean square displacement ∆(ϕ) obtained from replica method (black line), cavity reconstruction (blue squares), and MD simulations (red crosses).The replica result for ϕ d correspond to the point where the theoretical replica line has a square root singularity.(b) The theoretical results are consistent with the scaling form in Eq. (B15), but deviations are observed in the MD data close to ϕ d , due to the ambiguity in determining cage sizes when hopping is significant.

4 . 5 .
Test the dynamical behavior of ∆(t) (Eq.(B16)) below ϕ d , and determine the constants B and C. Note that here we have fixed all the other parameters, ϕ d , ∆ d , a and b from previous steps.Figure11(a) and (b)show that when ϕ < ϕ SER , Eq. (B16) is satisfied in the entire time regime; when ϕ < ϕ SER , it is only satisfied in the early β-relaxation regime.Using the same constant B, check the MCT dynamics above ϕ d using Eq.(B21).When ϕ is not too far away from ϕ d , ∆(t) does not strictly saturate to a plateau as predicted by MCT, and the scaling behavior of the intermediate time regime is modified.b.Breakdown of the Stokes-Einstein relationThe above analysis shows that the MCT scalings start to break down close to the dynamical transition ϕ d .To further investigate this property, we look at the scaling relation between the diffusion time τ D and the relaxation time

44 ϕ 47 FIG. 11 :
FIG. 11: Testing MCT scalings for the MK model in d = 6.Below ϕ d , the MD data for ϕ = 0.30, 0.35, 0.38, 0.40, 0.41, 0.42, 0.43, 0.435, 0.44 are fitted to Eq. (B16) (red lines) for (a) the early β-relaxation, and (b) the late β-relaxation together with diffusion, using fitting parameters B = 0.073 and C = 1.3.(c) Above ϕ d , the MD data ϕ = 0.45, 0.455, 0.46, 0.47 are compared to the early β-relaxation scaling in Eq. (B21) with the same value of B. A good agreement is found for the entire time regime when ϕ < ϕSER (black solid lines).When ϕ > ϕSER (green dashed lines), we only observe a good agreement for the early β-relaxation regime, which suggests that at later times hopping mixes with the MCT dynamics.

1 FIG. 13 :
FIG. 13: An example of hopping detection in d = 2 at ϕ = 2.40.(a) The particle trajectory clearly reveals two well formed cages.(b) The hopping between cages are visualized in the time series, with the two detected hopping times (dotted lines) at t = 3948.0and t = 7863.6.
− ϕ d .Therefore, the theoretical prediction of ∆ and ϕ d is fairly insensitive to the caging form and the second (or higher) moments of the cage size distribution, as well as to the method we choose (see SI Sec.IIA).

TABLE I :
Numerical values of characteristic densities and MCT exponents for the MK and the HS models.