Hierarchical landscape of hard disk glasses

We numerically analyse the landscape governing the evolution of the vibrational dynamics of hard disk glasses as the density increases towards jamming. We find that the dynamics becomes slow, spatially correlated, and starts to display aging dynamics across an avoided Gardner transition, with a phenomenology that resembles three dimensional observations. We carefully analyse the behaviour of single glass samples, and find that the emergence of aging dynamics is controlled by the apparition of a complex organisation of the landscape that splits into a remarkable hierarchy of minima as jamming is approached. Our results show that the mean-field prediction of a Gardner phase characterized by an ultrametric structure of the landscape provides a useful description of finite dimensional systems, even when the Gardner transition is avoided.


I. INTRODUCTION
The Gardner transition has a long history in the spin glass literature. It was first discovered in the context of mean-field spin glass models characterized by a random first order transition, where the glassy phase was found to destabilise at low temperature into a more conventional spin glass phase [1][2][3], characterized by a complex hierarchical organisation of pure states. The existence of a Gardner phase was more recently also proven for particle models in the absence of quenched disorder. In the mean-field limit, dense fluids first undergo a random first order transition to a glassy phase that may also transform in a marginal spin glass phase [4][5][6][7]. The existence of a Gardner transition in finite dimensional glass and spin glass models is still a debated issue, that neither theory nor simulations have been able to resolve [8][9][10][11][12][13][14].
This debate is particularly relevant in the context of hard sphere glasses, where the mean-field theory has been more precisely developed [5,7]. When the density is increased, the hard sphere fluid first undergoes a glass transition to an arrested glassy state [6]. Upon further compression, the hard sphere pressure increases very fast and diverges at a jamming transition. Jamming itself is a critical phenomenon, whose properties can be obtained in terms of marginal stability [15,16]. In the context of mean-field theory, a marginally stable Gardner phase is encountered between the glass and jamming transitions [4,7,17]. The marginal properties of the Gardner phase are thus key to properly describe the jamming criticality within mean-field theory. This is a very puzzling observation, as the critical exponents of jamming seem to remain roughly unchanged between d = ∞ and d = 2 [16,18], whereas instead the Gardner transition appears fragile against finite dimensional fluctuations [10-12, 19, 20].
These two sets of observations are hard to reconcile. On the one hand, the precise coincidence between meanfield predictions for jamming and numerical measurements in finite d would suggest that there exists a large pressure regime where mean-field predictions must become accurate [4,7]. On the other hand, it is unlikely that a Gardner phase can exist in dimensions as low as d = 2, where jamming criticality remains similar to the one in d = ∞ [15,16,18,21], suggesting that the Gardner transition, in itself, is not needed for a theory to capture the jamming criticality. Resolving this paradox is important, as mean-field theory can be used to tackle a large number of physical questions from a fully microscopic perspective, such as thermodynamic properties [4,19,[22][23][24], the structure of phase space, the evolution of vibrational dynamics [25,26], or the rheology of hard sphere glasses [27][28][29]. In addition, it is also important to understand under which conditions a complex free energy landscape may become physically relevant, in order to make novel predictions for experimental work dealing with the glassy dynamics of granular, colloidal, and molecular systems.
To tackle these questions we study the properties of hard disk glasses in the regime between glass and jamming transitions. For this system, it is known that jamming has the same critical properties as in any larger dimension up to d = ∞. It is also certain that a sharp Gardner transition cannot exist [10][11][12], although numerical simulations in d = 2 spin glasses indicate that the spin glass phenomenology can indeed be observed [30,31]. This suggests that a relatively sharp crossover should exist between two types of glass states, a simple and a marginal glass (see Refs. [22,24] for earlier hints for 2d hard disks systems), associated to the standard spin glass phenomenology and an emerging complexity of phase space, but that this phenomenology is not the result of a phase transition.
We find that hard disk glasses indeed display the same phenomenology as three dimensional hard spheres, associated with growing (but non-diverging) timescales and lengthscales, as they are compressed towards jamming [19,24]. By investigating the properties of single glass samples, we find that the emergence of collective vibrational dynamics corresponds to the development of a complex organisation of the landscape and a hierarchy of barriers and lengthscales that are strongly reminiscent of the mean-field description of spin glass phases [1, 3,17]. From our results, it is not possible to decide whether the complexity of the free-energy landscape near jamming stems from the marginality of jammed states, or, rather, whether the existence of an avoided Gardner phase transition is responsible for the original properties of the jamming transition. Our study suggests that it is not possible to observe jamming criticality without the Gardner physics. To our knowledge there is also no evidence that a finite dimensional particle system exhibits marginal stability far from a jamming transition [32], contrary to predictions obtained within mean-field theory [33].
This paper is organized as follows. In Sec. II, we introduce the model, simulation methods, and define the main observables used in this work. We present the numerical evidence for an avoided Gardner transition in Sec. III. In Sec. IV, we present a detailed analysis of the free energy landscape and the non-equilibrium dynamics observed in individual glass samples. In Sec. V, we demonstrate the emergence of a hierarchical landscape and show that the location of the Gardner avoided transition is strongly dependent on the chosen nonequilibrium protocol. In Sec. VI, we present results for larger systems and a careful analysis of finite-size effects. We conclude the paper in Sec. VII.

A. Model and methods
We simulate dense assemblies of N hard disks in a d = 2 square box with periodic boundary conditions. To avoid crystallization, polydisperse mixtures are adopted, where the diameters of the disks are drawn from the continuous distribution P (σ) ∝ σ −3 , with σ min /σ max = 0.45, σ min and σ max are the minimum and the maximum diameters, respectively. The simulation is carried out using Monte Carlo (MC) simulations, using two types of algorithms. Throughout this work, the time unit is defined as a MC sweep, i.e. a series of N attempts to perform a translational move for a single particle chosen at random. The units of length and energy are set by the average diameterσ = 1 N i σ i and temperature k B T , respectively. The pressure is reduced by the factor k B T ρ, where ρ is the number density. The dimensionless pressure Z = P/(k B T ρ) is extracted from the rescaled pair distribution, where r ij is the separation between particles i and j, D ij = (σ i + σ j )/2, with σ i the diameter of particle i. We perform simulations using different system sizes, N = 1024, 4096 and 16384.
Waiting time evolution of the volume fraction ϕ after the pressure is suddenly increased from Zg = 31 to the targeted value for N = 1024. For each pressure shown, a steady state value of the density is reached after tw ≈ 10 4 . For Z = 400 we show the compressions obtained for larger systems with N = 4096 and 16384, which exhibit no system size dependence. (b) Equation of state for the equilibrium fluid, and glass equation of state obtained by compression of an initially equilibrated system at ϕg = 0.820. For reference, we show the location of the mode-coupling crossover, and of the avoided Gardner transition ϕG.
Our first task is to prepare equilibrium configurations of the system at large packing fraction, defined as ϕ = π i σ 2 i /(4L 2 ). Following previous work [34,35], we use a SWAP MC technique which combines translation MC moves for single particles and SWAP moves for randomly chosen pairs of particles. To optimize the thermalization efficiency, the ratio of SWAP versus translational moves is set to 0.2 [35].
At thermal equilibrium, the equation of state of the polydisperse hard disk fluid is well described by the following empirical equation of state [36], where M 1 and M 2 are the first and second moments of the diameter distribution. The measured equation of state in shown in Fig. 1. By measuring the translational dynamics of the equilibrium supercooled fluid, we observe that the dynamics slows down rapidly as ϕ increases. By fitting the relaxation time to a power law divergence, we estimate the location of the mode-coupling crossover to be ϕ M CT = 0.795, see Fig. 1. Above this density, particle diffusion becomes very slow when standard MC simulations are employed, but the SWAP MC method allows us to prepare much denser systems in equilibrium conditions. In this work, we prepare N s independent equilibrium configurations at ϕ g = 0.820, where the pressure is Z g ≈ 31. At this packing fraction, the standard MC dynamics is nearly frozen and diffusive dynamics is not observed over the duration of our simulations. Therefore, each of the N s configurations corresponds to a different glass sample.
Having prepared equilibrium glass configurations, we then employ N P T standard MC simulations to rapidly compress the system. In our simplest protocol, the pressure is instantaneously changed from its equilibrium value Z g = 31 to a much larger pressure Z > Z g to rapidly compress the system. This is analogous to an instantaneous temperature quench. We denote t w the 'waiting time' spent since the pressure has been suddenly increased to its targeted value.
It is important to perform rapid compressions to avoid the formation of localised "defects" [20,32,37], which correspond to small groups of particles which may undergo a small rearrangement in the initial steps of the compression when the pressure is not yet very large. Our simulations indicate that this happens with a low but finite probability, and this is distinct from the more collective Gardner physics we wish to analyse at larger pressures.
In order to obtain very fast compressions, we adjust the amplitude of individual compression moves in the N P T MC protocol to obtain an acceptance rate of about 0.1. We also frequently attempt compression moves, once every 100 MC steps. The amplitude of the translational moves is carefully chosen to have an acceptance rate of about 0.6. Empirically, we find that these choices optimise both the measured compression rate and the CPU time. Examples of fast compressions from Z g = 31 to Z = 100, 200, 400 and 10000 are shown in Fig1. For each final pressure, we observe that after a short time of about 10 4 MC sweeps, the density reaches a steady state value that increases with Z. The relation Z(ϕ) in that steady state defines the glass equation of state, which we report in Fig. 1.
As noted in a recent study [19], the efficiency of MC compression moves decreases rapidly as the system size increases. If one keeps the compression parameters fixed, this leads to an unwanted slowdown of the compression time, which scales proportionally to N . To keep the compression rate roughly constant, we introduce a slightly different compression move for larger systems in which we randomly select a subset of N = 1024 particles (corresponding to the smallest system size studied in this work) at each compression move, and attempt to change their sizes. As can be seen in Fig. 1, this approach removes the N -dependence of the compression protocol. It is therefore much easier to compare results obtained for different system sizes.
In order to understand the structure of the landscape of hard disk glasses, we need to perform an average over the N s distinct glass samples produced in equilibrium, and for each glass sample, we also need to perform a thermal average over a large number N c of "clones". Each clone is prepared by replicating the particle positions of the glass at density ϕ g = 0.820, but the different configurations are simulated using different random numbers for the MC moves, thus leading to independent trajectories within the same glass metastable basin. Hence, for each simulation, we need to choose a pair (N s , N c ) for the number of glass samples and clones. To study macroscopic quantities, we have used two sets of values, (N s , N c ) = (100, 20) for general aspects of the physics related to the Gardner transition, while we use (N s , N c ) = (1, 100 − 1000) when we study individual glass samples.

B. Physical observables
In this section, we introduce the physical observables used throughout this work. Because we want to perform various levels of average, we need to be precise with the notations. We consider N s independent glass samples, that we design using Greek letters: α = 1, · · · , N s . For each glass sample α, we introduce N c clones, which we design using latin letters: a = 1, · · · , N c . And we use i = 1, · · · , N to design the particle index.
We then introduce two types of averages. We use an overline O to define a disorder average over the N s distinct glass samples, and use brackets O to denote averages over the clones within a given glass sample. To avoid complex notations, we use brackets for both averaging single clone quantities (over N c clones) and two-clone quantities (over N c (N c − 1)/2 pairs of clones).
The basic observable we consider is the particle coordinate. We define R α a,i (t) the coordinate of particle i within clone a and sample α at time t. Tracking the motion at single particle level, we can measure the full mean-squared displacement (MSD), where t w is the waiting time since the beginning of the compression. It is important to keep the two-time notations, since the emergence of aging dynamics is a powerful tool to probe the Gardner transition [24]. There is an important complication which emerges in two-dimensional solids, which takes the form of longrange cooperative motion [38]. These so-called Mermin-Wagner fluctuations happen at finite temperature and pressure, and may dramatically affect the single particle dynamics we wish to study in the glass phase [39,40]. These fluctuations happen in addition to the vibrational motion within the cage, but they do not change the nature of the glass phase itself. To disentangle the true emergence of marginal glasses (which involve spatially correlated motion) from the more mundane Mermin-Wagner fluctuations, we introduce cage-relative coordi- nate as follows [22,40] where N i is the number of nearest neighbors of particle i. We define neighbors as particles separated by a distance smaller than 2σ, which essentially includes the first shell of neighbors. Using the cage relative coordinates, we can then define the cage-relative mean-squared displacement as We also introduce the instantaneous MSD ∆ α a (t w , t) for a given clone a within sample α, ∆ α (t w , t) = ∆ α a (t w , t) for a given sample α, such that ∆(t w , t) = ∆ α (t w , t).
In Fig. 2, we show some representative data for both the full MSD D(t w , t) and the cage relative MSD ∆(t w , t) after a quench from Z g = 31 to Z = 100. Not only is the MSD much larger than the cage relative one, but its system size dependence is very different, as D(t w , t) seems to grow without bound as N increases [40]. By contrast, the cage relative MSD is insensitive to system size, showing that Mermin-Wagner fluctuations have properly been removed from our measurements. In the rest of this work we quantify all the dynamic observables using cage relative definitions only. Whereas this distinction is crucial to quantify displacements (and, thus, dynamics), it is immaterial as far as particle positions are concerned.
The study of complex landscape makes heavy use of two-clone quantities, in the spirit of spin glass physics where the order parameter is an overlap between cloned configurations. In the present case, we introduce the dis-tance between clones a and b as Similarly to the MSD, we introduce the instantaneous value ∆ α ab of the distance between clones a and b within the sample α, and the typical distance between ∆ α AB = ∆ α ab within a given sample α. From now on, we omit the time label for those one-time observables.
We are interested in the spatial fluctuations of the distance between clones. It is convenient to measure them in the Fourier domain and access the following structure factor, where the distance field is defined in the Fourier domain as Here both the indexes a and b range from 1 to N c , to the keep the symmetry. We have implicitly performed an angular average in Eq. (7). For visualisation purposes, we represent snapshots of the real space version of the displacement field, A simpler quantity measuring the global fluctuations of the distance field is to obtain the q → 0 limit of the structure factor S AB (q), which corresponds to the following susceptibility [32] where the normalisation is given by to ensure χ AB ∼ O(1) for a fully uncorrelated field. It is also useful to decompose the susceptibility χ AB into two distinct components, which, respectively, quantify the fluctuations within a sampleχ AB and the sampleto-sample fluctuations χ AB , defined as follows: It is easy to verify that χ AB =χ AB + χ AB . Finally, the susceptibility within a single sample α is measured as, As before, we can see that The final set of variables that we study in this work corresponds to performing a time average over the coordinates. To this end, we define time averaged cage relative positions for each particle after a quench: Here, we average the coordinates between time t w and t w + τ . Computing the difference between clones, we have The physical idea is that the long-time average will provide an estimate of the average position around which particle i performs vibrations in clone a within sample α.
We expect that all average positions would be the same from one clone to the next in a simple glass, whereas in a marginal phase characterized by a complex landscape, the average positions can be distinct for different clones. Finally, the time dependence of ∆ AB (t w , τ ) will inform us about how slow is the exploration of the landscape by the different clones within a given glass sample.

III. EVIDENCE OF AN (AVOIDED) GARDNER TRANSITION
In this section, we use the system size N = 1024, and perform averages over both samples and clones. We use (N s , N c ) = (100, 20) to obtain a reasonable statistical accuracy (requiring these numbers to be as large as possible) within the available computer resources [14,19,41]. We follow earlier work and compress the system instantaneously using the protocol shown in Fig. 1, starting from the equilibrium state at (ϕ g , Z g ) = (0.820, 31) along the glass equation of state to pressures Z = 100, · · · , 10000.

A. Signatures of a Gardner crossover
In Fig. 3 we reproduce several key signatures of the Gardner transition for the 2d hard disk glass, reported earlier for the 3d hard sphere system, which signal the emergence of marginal stability and a rough free energy landscape at large enough pressure [24,32]. For our 2d system, this can only signal a sharp Gardner crossover, not a real phase transition. Our first observation is the evolution of the MSD ∆(t w , t) and distance ∆ AB (t w ) following a quench, as shown in Fig. 3 (a). For small pressure, both quantities rapidly converge at long time to the same limit, indicating that the typical distance traveled by one clone is similar to the typical distance between clones. This shows that the structure of the glass remains simple. For pressures larger than Z G ≈ 200, by contrast, the MSD ∆(t w , t) is strongly dependent on the waiting time, indicating slow dynamics. And the MSD no longer converges to the typical distance between clones, which indicates that the dynamics within the glass is no longer ergodic.
From these data, we compute the long time limit (in practice, we use 10 8 MC sweeps) of both ∆ and ∆ AB obtained after a quench. The results are shown in Fig. 3 (b), as a function of the pressure where the quench is performed. As expected, we find that for Z < Z G both quantities converge to the same long-time limit, but depart gradually from one another as Z increases beyond Z G .
In Fig. 3 (c), we show the probability distribution of distances P (∆ α ab ) measured at long waiting time t w = 10 8 and various pressures. The distribution has a simple Gaussian form for Z < Z G , but becomes strongly non-Gaussian for Z > Z G , which can be taken as another sign that a complex structure of the glass basin emerges beyond Z G characterized by a broad distribution of dis-tances between the clones, as found before [24]. We observe that the distribution P (∆ α ab ) first broadens as Z → 1000, but then it narrows down again, as the pressure increases further, Z > 1000, which corresponds to a non-monotonic evolution of the susceptibility χ AB (Z), as shown in Fig. 3 (d). In addition, for intermediate pressures Z = 200 · · · 1000, the distribution P (∆ α ab ) develops a clear shoulder. This indicates that different pairs of clones can be either very close, or, instead, very distant from one another, and suggests already that some kind of self-organisation of the glass basin emerges near Z G .
Finally, we also employ a procedure used before [24] to locate the Gardner transition, and compute the skewness Γ AB (Z) of the distributions P (∆ α ab ), see Fig. 3 (d). The skewness presents a clear peak, but we find that the position of the peak depends on the chosen timescale, shifting from Z ≈ 200 at t w = 10 7 to Z ≈ 300 at t w = 10 8 with no sign of a saturation.
To summarise, the results reported in Fig. 3 indicate that all key signatures of the Gardner transition seem to survive in our 2d hard disk system, which appears to behave similarly to the 3d and mean-field systems [4,19,24]. Notice that some of these important signatures are absent for the nearly 1d system studied in Ref. [20] and 3d and 2d systems of soft spheres [32,37], where collective aging effects are not observed. At this stage, it is unclear that the Gardner transition is indeed avoided in our system, but the absence of a sharp phase transition will become more obvious in Sec. V.

B. Growing time scales and length scales
Next, we characterize the time scales and length scales associated to the emergence of a rough landscape inside hard disk glasses.
First, we follow the time evolution of the global fluctuations, χ AB (t w ) in Fig. 4 (a). For all pressures, the initial growth of χ AB up to t w = 10 4 is identical, which corresponds to the initial compression of the system. For Z = 100 · · · 200, the susceptibility rapidly reaches a constant limiting value, which increases by a factor of about 2. For Z > 400, no steady state value can be observed, and the susceptibility remains time dependent. It grows to reach large values of about 20 for Z = 1000 with a slow time dependence at long times compatible with a nearly logarithmic increase with time (the data are shown in a log-lin representation in Fig. 4 (a)).
Interestingly, for Z = 3000 and 10000 the susceptibility becomes somewhat smaller, indicating a nonmonotonic behaviour with Z, which was already apparent in Fig. 3 (d). This overall behaviour is qualitatively similar to observations performed in 3d spin glass models [19]. Roughly speaking, the non-monotonicity corresponds to the competition between the emergence of an increasingly complex landscape that tends to increase the susceptibility, and the slowing down of the dynamics which makes the exploration of that landscape more difficult at large

pressures.
A similar non-monotonic behaviour can be observed in real space, as shown in Fig. 4 (b) which represents the spatial fluctuations of the distance field S AB (q) for a time t w = 10 8 and various pressures. For each pressure, the structure factor has a peak near q ≈ 6, which corresponds to the interparticle distance. More interesting is the development of a peak at q → 0, which corresponds in real space to spatial fluctuations happening at length scales much larger than the interparticle scale. As Z increases, this peak becomes higher, and a saturation to a low q plateau, which is clearly visible for Z < 200 be-comes difficult to observe. This directly indicates that the fluctuations of the distance field between clones inside glass basins become spatially correlated over a large correlation length that grows as the pressure increases. Just as for the susceptibility, we observe that the height of a low-q peak has a non-monotonic behaviour with Z, reaching its maximum near Z = 1000. For N = 1024, the linear size of the system is L ≈ 30, and so the observations in Fig. 4 (b) indicate that a collective correlation length of at least L/2 ≈ 15 develops in the system near Z G . Given the relatively modest range of wavevector at hand, it is difficult to quantitatively extract the value of a correlation length from these data, a task which is known to be difficult in glass-formers [42]. We will nevertheless confirm these observations using larger systems in Sec. VI.
We now turn to the growing time scales associated to the signatures of a Gardner transition. In Fig. 4 (c) we report the evolution of the distance ∆ AB (t w , τ ) between time-averaged positions, see Eq. (15). We choose t w = 4×10 5 , which is large enough for the density to converge to its steady state value. For Z = 100 · · · 200, the distance decays to zero at long times, indicating that the structure of the glass basin remains simple at those pressure, and all clones can freely explore the entire basin. However, the time scale for this exploration increases rapidly with Z, from τ ∼ 10 6 at Z = 100 to τ > 10 8 at Z = 200. For even larger pressures, ∆ AB (t w , τ ) does not decay to zero, indicating that different clones are now confined in different restricted part of the glass basin that are not dynamically accessible at these pressures. In Sec. IV we revisit the physical meaning of this increasing time scale in terms of emerging barriers that appear near Z G .

C. Thermal average versus disorder average
In the quantities presented above, the fluctuations have two distinct origins: there are clone-to-clone fluctuations within each glass sample, and there are sample-to-sample fluctuations because each glass corresponds to a distinct particle packing.
We have performed a more detailed study of a small number (N s = 10) of individual glass samples, using for each glass a large number of clones, N c = 100. We find that each glass behaves dynamically similarly, and variations between samples are mostly quantitative. In particular, each individual sample exhibits a clear aging behaviour, indicating that for all samples, a non-trivial collective dynamics emerges at large enough pressures. Strikingly, all the key signatures of a Gardner transition shown in Fig. 3 can be reproduced for a single glass.
However, the precise location of the crossover pressure differs from one sample to another. For N = 1024, the range of these variations is estimated to vary from Z G ≈ 100 to Z G ≈ 600, for a set of glass configurations drawn from our pool of equilibrium systems at (ϕ g , Z g ) = (0.820, 31).
In order to estimate the relative influence of the two sources of the fluctuations, we decompose the susceptibility χ AB into its two contributions defined in Eq. (12). For a quench to Z = 1000, we show the result in Fig. 5. The data in Fig. 5 indicate that the two sources of fluctuations are of the same order, and we find in addition that the ratio χ AB /χ AB increases only weakly with the pressure.
This result suggests that the analysis of single glass samples provides an important part of the fluctuations associated to the avoided Gardner transition. In addition, averaging over distinct glass basins washes out the detailed information related to the structure of the landscape within each glass. In the next two sections, we thus perform a more detailed analysis of the evolution of the landscape in a single glass basin.

IV. DETAILED ANALYSIS OF A SINGLE GLASS SAMPLE
In this section, we analyse the physics associated to the Gardner transition for a single glass sample selected randomly. The protocol and parameters are the same as in Sec. III. We generate N c = 1000 clones to understand better the structure of this particular glass basin. For this sample, the key signatures shown in Fig. 3 are reproduced with Z G ≈ 200. We first analyse the physics in the vicinity of the Gardner crossover near Z G and then move to larger pressures.

A. Emergence of sub-basins at Gardner crossover
To understand the emerging structure of the glass basin, we use a large number N c of clones issued from the same equilibrium configurations at Z g = 31 and we independently compress them instantaneously at some pressures Z > Z G . After a given waiting time, we then measure the relative distance ∆ α ab between clones a and b for this sample α. It is useful to cluster the clones according to their relative distances. With the simple procedure explained in Appendix A, we can relabel the clones and construct heat maps such as the ones shown in Fig. 6 (a), where the color is coding for the distance between the clones. For sample α in the condition of Z = 600 and t w = 2 18 × 100, the map reveals a clustering of the clones into two subpopulations, which we refer to as group m and group n . We present maps for different times and pressures below. This indicates that, after a fast compression, the different clones now explore two distinct parts of the glass basin that are not easily accessed dynamically. The underlying physical picture is that of a glass basin that is now fractured into two sub-basins separated by a barrier between them.
We now measure the corresponding probability distribution functions P (∆ α ab ) shown in Fig. 6 (b). As ex-pected from the structure of the heat map in Fig. 6 (a), the full distribution of distances between clones is bimodal, which reflects the clustering of the clones into two families. Clones within a given group are close to each other, while they are relatively far from all clones belonging to the other group, which leads to two peaks in the distributions. Decomposing the clones into group m and group n , we also measure the distribution of distances between clones belonging to the same group in Fig. 6 (b). These distributions are nearly Gaussian, centered around a small value roughly equal to the cage size at the pressure Z = 600. Finally, in Fig. 6(c1-c3), we visualise the displacement fields between clones, ∆ α ab (R), generated by choosing a pair of clones belonging to group m and group n as in Fig.6  (c1), or to the same group m (Fig.6 (c2)) or group n (Fig.6  (c3)). The method to construct these snapshots is detailed in Appendix B. We observe long-range correlations in Fig.6 (c1) which are instead much shorter-ranged in panels Fig.6 (c2-c3). These snapshots suggest that the emergence of large length scale across the Gardner crossover Z G is a direct consequence of the fracturation of the glass basin into multiple sub-basins.

B. Physical interpretation of aging dynamics
We now extend the analysis of the heat map to a larger range of timescales and pressures. Our goal is to explore the connection between the emergence of barriers within the glass basin and the aging dynamics reported in Fig. 3 (a) and the growing time scale reported in Fig. 4 (c).
We first display the evolution of the heat maps for several pressures Z and various waiting times t w in Fig. 7. The first row of heat maps shows the results obtained for a relatively low pressure, Z = 100 < Z G , where aging is absent and the vibrational dynamics is featureless. Accordingly, we observe that all clones within the glass sample resemble each other and can freely explore the basin. As a result, the heat map is both featureless and time independent. For Z = 300, 400 and 600, one can see the emergence of the two subpopulations group m and group n discussed in Sec. IV A. Fixing the pressure, the time dependence observed shows that as t w increases the size of the group m enlarges, whereas group n shrinks. In addition, this dynamics slows down as Z increases so that for a fixed t w , the number of clones surviving within group n increases with Z.
The emergence of a bimodal structure that slowly fades away with time is confirmed in Fig. 8 (a-c), where we report the probability distributions function corresponding to the maps shown in Fig. 7. In agreement with the maps, the distributions develop a bimodal structure at intermediate waiting times, where the peak at small distance corresponds to the amplitude of thermal vibrations, whereas the peak (or shoulder for Z = 300) at large distance corresponds to the typical distance between the two sub-basins identified by group m and group n . For a fixed Heat maps of distances between clones ∆ α ab (tw) for sample α with N = 1024, organized by Nc = 100 clones. The waiting time increases from left to right, the pressure from top to bottom, as indicated. Whereas the maps at Z = 100 < ZG are featureless and time-independent, they exhibit an organisation into two sub-basins for Z ≥ 300, which slowly disappears as tw increases. t w , increasing the pressure tends to broaden the distributions, which is related to the growth of the susceptibility χ AB near Z G shown in Fig. 3 (d).
The time evolution of these bimodal distributions can be easily described by computing the evolution of the fraction of clones P n (t w ) that survive in the less stable group n after a time t w . This can be done when the pressure is in the vicinity of the Gardner crossover. We report the time evolution of P n (t) in Fig. 8 (d). In agreement with the time evolution of the maps and of the distributions, we find that group n becomes empty after t w ≈ 10 7 for Z = 300, while the corresponding population survives much longer and remains finite even after t w = 2 × 10 8 for Z = 600, suggesting that the dynamic evolution from one basin to the other becomes slower as the pressure increases. The time decay of P n (t w ) is quantitatively similar to the evolution of the distance between average positions ∆ AB (t w , τ ) discussed in Fig. 4 (c).
We now close the loop between the structure of the glass landscape, the emerging barriers near Z G , and the physical interpretation of the aging dynamics arising for Z > Z G . A careful analysis of the MSD at the single clone level reveals that the existence of two subpopulations, group m and group n , implies also two types of MSD. Clones that belong at time t w to the more stable group m are characterized by simple MSD smoothly increasing to reach the cage size, and never leave that long-time plateau. These clones do not contribute to the aging dynamics. By contrast, clones belonging at time t w to the less stable group n may undergo some activated FIG. 9. Connection between landscape structure and aging dynamics for sample α at Z = 600. (a) Blue lines: aging dynamics for a single clone a undergoing activated event from groupn to groupm via MSD ∆ α a (tw, t). Red square: the distance ∆ α ab (tw) between clone a and a randomly chosen clone belonging to groupm jumps discontinuously after the rearrangement. (b) The cage-relative displacements by the configurations of clone a measured before and after the activated event detected in (a). This displacement field is very similar to the snapshot of ∆ α ab in Fig. 6 (c1). event at time t + t w later to join group m . The example of such a clone, clone a, is shown in Fig. 9 (a) (blue lines). The MSD first grows to the plateau corresponding to the cage size, but then a rearrangement takes place that increases the MSD further to a second plateau value. Notably, this second plateau corresponds to the typical distance ∆ α ab between group m and group n (red squares). After the rearrangement, the clone a belongs to group m .
Obviously, the time when the sudden rearrangement takes place in clone a fluctuates from one clone to another. Therefore, after averaging over clones and over glass samples, the MSD eventually exhibits the smooth aging behaviour reported in Fig. 3 (a).
We also confirm that the spatial correlations that develop during the aging are also directly controlled by the underlying structure of the phase space. We show in Fig. 9 (b) a snapshot representing the real space displacement field corresponding to the activated event detected in clone a in Fig. 9 (a), by measuring ∆ α a (t w , t) for well-chosen times. This displacement field is very similar to the distance field shown in the snapshot of Fig. 6 (c1), corresponding to ∆ α ab between clones belonging to group m and group n . The remarkable similarity between Fig. 9 (b) and Fig. 6 (c1) indicates that the long-time aging dynamics in the vicinity of the Gardner crossover is fully controlled by activated relaxations between the emerging sub-basins.
Remarkably, this aging dynamics is spatially correlated over a relatively large length scale, comparable to the system size L/2 ≈ 15, and involves a large number of particles. Thus, we conclude that the sub-basins structure that we detect for 2d hard disks, is qualitatively distinct from the highly localised defects found in both nearly 1d hard disks and 3d and 2d soft spheres [20,32,37].
In summary, the detailed analysis of a single glass sample shows that as Z increases across Z G , the glass basin fractures into two sub-basins separated by a collective barrier. After a sudden quench, clones are initially randomly trapped into one of the sub-basins, and as t w increases, the clones belonging to the less stable sub-basin undergo a collective, spatially correlated rearrangement to join the more stable sub-basin. This physical picture explains the emergence of growing timescales and lengthscales in the vicinity of the Gardner crossover. In particular, if fully accounts for the growing timescale detected directly in Fig. 4 (c).

C. Proliferation of states at large pressures
We now present results of the same sample α, obtained by quenching to much larger pressure, Z = 10000 Z G ≈ 200. We again use N c = 1000 clones and record the time evolution of the probability distribution function P (∆ α ab ); see Fig. 10 (a). The distribution shows a peak at a large value (much larger than the cage size), and it starts to develop a tail at smaller values as t w increases. However, the simple bimodal structure found at smaller pressures is now totally absent. Physically, this means that after a quench to a large pressure, all clones are trapped within a large number of distinct structures, so that two clones chosen at random are typically separated by a large distance. Thus, instead of finding two subbasins separated by a single large barrier, the N c clones in sample α are now trapped in a number of O(N c ) distinct states. Correspondingly, when we perform a cluster analysis of the clones to construct a heat map, we find no interesting structure emerging from the map, see Fig. 10 (b). We conclude that the number of basins in- creases dramatically upon approaching jamming. The crossover between the simple bimodal structure observed for Z = 600 in Fig. 8 and the very complex structure for Z = 10000 in Fig. 10 is not sudden: the change of pressure between these two state points is quite large (a factor of 16), and quenches performed at pressure values intermediate between those reveal a gradual change from one limit to the other.
Tracking the time evolution of the clones, we observe that there is still some activated dynamics taking place, which is responsible for the aging dynamics observed at the level of the averaged MSD. If one follows the time evolution of the packing fraction in individual clones, as shown in Fig. 10 (c), one can observe that the density performs small jumps corresponding to activated events taking place within each clone. Similarly, the time evolution of the MSD for individual clones shown in Fig. 10 (d) shows that the MSD exhibit small jumps, which presumably correspond to crossing small barriers between neighboring states. The amplitude of these jumps is smaller than at lower pressure, suggesting that aging at large pressures corresponds to crossing smaller barriers.
In summary, a direct quench to large pressures reveals a proliferation of accessible states, indicating that the glass basin breaks at large Z into many sub-basins. Dynamically, the clones at such large pressure can only cross small barriers and explore a very small part of the entire glass basin. Altogether, these results therefore point towards a hierarchical organisation of the glass landscape where the large sub-basins that appeared at smaller pressure are themselves broken into many smaller sub-basins, involving therefore a hierarchy of barriers.
The analysis of this complex structure emerging in a single glass sample is not easy. To get more insight into this structure, we turn to more complicated compression protocols to more evidently reveal the hierarchical structure of the landscape of hard disk glasses.

V. EVIDENCE OF HIERARCHICAL LANDSCAPE APPROACHING JAMMING
In the previous section, we observed a fracturation of the glass basin of sample α in two basins near the Gardner crossover, and a proliferation of dynamically inaccessible states at larger pressure. In this section, we demonstrate that the glass basin is actually organised in a hierarchical way where, as the pressure increases, the main basin breaks into sub-basins that themselves break into several sub-basins, etc. This description suggests the existence of a broad distribution of barriers separating these different states.

A. Compressing the system in several steps
To illustrate the finer structure of the glass landscape, we continue our exploration of a single representative sample α, but change our protocol. As an example, we use the following history: Z = 31 → 600 → 1000 → 10000, where the applied pressure is kept constant for a duration dt = 2 18 × 100 at each step. In Fig. 11 (a) we present the evolution of the susceptibility χ α ab measured in sample α during this multi-step protocol, and the corresponding heat maps of distance in Fig. 11 (b-e).
During the first stage Z = 31 → 600, we recover the physics obtained during a direct quench across Z G , namely a slow growth of the susceptibility associated to a bimodal organisation of the states, as discussed in Sec. IV. If we waited a larger time at this pressure, group n would eventually disappear and all clones would be located inside group m , which is the most stable sub-basin. Such a protocol is used below in Sec. V B. Here instead, we again increase the pressure to Z = 1000 after a waiting time dt when both group m and group n are still populated, followed by another quench to Z = 10000 at time t w = 2dt. In this multi-step protocol, we find that the susceptibility increases rapidly after each of the steps, indicating that the distribution of distances between clones broadens at each step.
The most striking observation comes from the corresponding heat maps in Fig. 11 (c-e), which show that the bimodal structure observed at Z = 600 fractures into a finer structure hierarchically. Specifically, we observe that each blue square from Fig. 11 (b) at Z = 600 becomes decomposed into sub-squares at Z = 1000 in Fig. 11 (c), which themselves break into smaller clusters at Z = 10000 in Fig. 11 (d). Whereas the clones were scattered across many different states with little organisation between them after a direct quench to Z = 10000 (recall the heat map in Fig. 10 (b)), the states reached at the same pressure in a multi-step compression protocol are now clearly clustered and organised in a hierarchical manner in Fig. 11 (d). Physically the difference is that the subpopulation of clones in different sub-regions of the landscape are better thermalised in a multi-step protocol than by a direct quench to a large pressure, where the clones are instead trapped in a myriad of states and can barely evolve dynamically across the landscape. The difference between these two protocols is also clear if one compares the susceptibilities measured after the same time by direct quench or in a multi-step protocol, see Fig. 11 (a).
In Fig. 12 we provide a real space view of the hierar-FIG. 12. Distance fields ∆ α ab (R) for sample α at Z = 10000, using the four pairs of clones shown in Fig. 11 (d). From (i) to (iv), the various levels of the hierarchy of states are associated to larger displacements and longer-ranged spatial correlations.
chical landscape shown in Fig. 11 (d) at Z = 10000. To this end, we select pairs of clones belonging to the different layers of states shown in the maps, see the notations (i-iv) in Fig. 11 (d), and we build a snapshot showing the corresponding distance fields ∆ α ab (R). In the snapshot we observe that the overall shade increases as the main distance between clones increases, and also that the spatial displacements become correlated over increasing lengthscales. This suggests that the hierarchical organisation of the states corresponds, in real space, to a hierarchy of length scales, presumably associated to a hierarchy of barriers between these states. This hierarchical structure should give rise to rejuvenation and memory effects [43], as discussed extensively in the spin glass literature. The various compression steps discussed above are accompanied by an aging dynamics corresponding to barriers of smaller and smaller amplitude, akin to rejuvenation effects. Given the distinct time scales mentioned above one should find the memory of the previous state by decompressing the system. This is shown in Fig. 11 (a) (red squares). When the pressure is decreased from Z = 1000 to Z = 600, the susceptibility goes back to the value it had at the end of the first compression step at Z = 600, which indicates memory. Accordingly, the heat map in Fig. 11 (e) is similar to the one in Fig. 11 (b), showing that the clones have retained the memory of the organisation they had at that time.
It is remarkable that such a complex organisation of the glass basin emerges in our system of 2d hard disks, given that no Gardner transition is expected to take place at finite pressure in this system. In the regime of pressure considered here Z ≤ 10000, the system is still pretty far from jamming, in the sense that the jamming criti- cality itself is not yet well-developed (for instance, the power laws characterizing the pair correlation function at contact cannot be observed at all), but our results are quite consistent with the organisation of the landscape obtained in mean-field descriptions of spin glass phases [1]. In particular, it is striking that our heat maps resemble the pictures drawn to describe the ultrametric organisation of phase space in spin glasses [44,45].

B. Direct evidence for absence of Gardner transition
The description of the glass landscape as a hierarchical organization of sub-basins suggests a strategy to maintain equilibrium up to the pressures that are above the Gardner crossover shown in Fig. 3. The simplest way to detect the Gardner crossover is by monitoring the difference between the long-time limits of the MSD and the averaged distance between clones. For the single sample α considered extensively in this work, we report those data in Fig. 13 (a), which confirms that Z G ≈ 200 for this sample.
Next, we again use the idea of a multi-step protocol, but now choose carefully the intermediate pressures and duration spent at each step as follows. We first perform a quench from the initial liquid state at Z g = 31 up to Z = 300. At this pressure, for the single glass sample α, we have shown in Fig. 8 (d) above that it takes about t w ≈ 10 7 for all clones to cluster inside the most stable sub-basin, which we called group m . To make sure that all the clones end up in the same part of the landscape, we spend a time t w = 2 × 10 8 at this pressure.
Using the same set of clones, we then repeat the protocol of quenching the system to larger pressures, and we again measure the long-time limits of both ∆ and ∆ α ab in this two-step protocol. Notably, we observe that they now start to differ from one another only when Z becomes larger than Z ≈ 900. In the language used in the previous section, it corresponds to the pressure where the underlying sub-basin of group m splits into distinct sub-basins that become dynamically inaccessible above Z ≈ 900. This resembles a second avoided Gardner transition.
We then perform a three-step protocol Z g = 31 → 300 → 1000 and wait a long enough time at Z = 1000 (we use t w = 2 × 10 8 ) so that all clones again gather into the same sub-basin at this pressure at the end of the second step. Then, we again quench further the clones to larger Z. The data in Fig. 13 suggest that a third Gardner crossover now occurs at Z ≈ 2000 in this protocol, which is about 10 times larger than the Gardner crossover detected using a direct quench.
In those three protocols, we measure the susceptibility χ α ab at a fixed time t w = 2 × 10 8 , as reported in Fig. 13  (b). The non-monotonic behaviour observed before for a direct quench is again present in the multi-step protocols, but the location of the maximum is shifted to larger pressures. More interestingly, the amplitude of the maximum decreases as well, which agrees well with the above findings that the hierarchy of states that appear at larger pressures correspond to smaller and smaller lengthscales as well (recall the snapshots in Fig. 12).
In Sec. III, we obtained key signatures for the existence of a Gardner transition in our 2d system, with a transition possibly rounded by finite time effects. Here, we observe that by crossing the transition more slowly, these signs are shifted to another location and become less prominent. If the crossover was underlaid by a genuine phase transition, a modest shift of its location would be expected, together with sharper dynamic signatures. The results in Fig. 13 are thus direct evidence for a transition that is avoided. In the following section, we additionally show that increasing the system size does not provide diverging correlation lengthscales either.

VI. RESULTS FOR LARGER SYSTEMS
In this final section, we expand the above analysis to larger system sizes to understand how robust the results for N = 1024 can be. We study how the measurements related to the Gardner physics change with increasing the system size to N = 4096 and 16384, quenching the system from (ϕ g , Z g ) = (0.820, 31) to different pressures, using (N s , N c ) = (100, 20) and maximal waiting time t w ∼ 3 × 10 7 .

A. Large but finite correlation length
First, we present the N dependence of several quantities after a quench to both Z = 100 and Z = 400 starting from the equilibrium system at Z g = 31 in Fig. 14, which implies that the results presented earlier for N = 1024 do not strongly depend on the system size. The susceptibility in Fig. 14 together with the MSD (not shown) confirm that in large systems the dynamics is fast (no aging) and trivial (the susceptibility remains O(1)) at Z = 100. Looking at the structure factor S AB (q) the absence of finite size effects is also obvious (data for larger N nearly superpose) and a clear plateau at low-q is observed suggesting a small correlation length of a few particle diameters at most.
The situation is different at Z = 400, where the MSD shows aging (not shown) and the susceptibility is large (of magnitude 15 at t w = 10 7 ), and the low-q part of the structure factor saturates at a much lower wavevector, indicating a much larger correlation length of the distance field. Therefore, the signatures of the Gardner crossover found in N = 1024 are also observed in larger systems. A key observation from Fig. 14 is the presence of a finite size effect between N = 1024 and larger sizes, compared with the superposition of curves for N = 4096 and 16384. This suggests that the maximal correlation length observed in our system is of the order of the linear size of the system with N = 4096, that is, L/2 ≈ 30, which is consistent with the length extracted from the low-q part of S AB (q). The presence of such a large correlation length provides the explanation why the physics related to a phase transition is obvious in our data, as it requires much larger system sizes to realise that the correlation length is actually not divergent. It would be interesting to revisit 3d hard sphere glasses to understand how large the correlation length can become in d = 3.
Regarding timescales, we revisit the measurements of the distance between averaged positions, ∆ AB (t w , τ ) from (Eq. (15)), for different system sizes in Fig. 14 (c). As for other quantities, we again find that the timescales extracted from these data change weakly with N , indicating in particular that the measured relaxation times do not appear to grow with the system size.

B. Self-averaging and landscape in large systems
The data in Fig. 14 indicate that the maximal correlation length in the 2d hard disk glass after an instantaneous quench is about ξ ≈ 30, the linear size of the system N = 4096. Hence, a larger system with N = 16384, corresponding to L ≈ 120 should be understood as the superposition of several independent subsystems. Correspondingly, the numerical measurements with N = 16384 should be self-averaging, with the fluctuations within each glass sample, and sample-to-sample fluctuations should be much smaller. A second consequence is that it should become less relevant to discuss the landscape of single glass samples in such large systems, because the total landscape is in fact the convolution of independent landscapes. We therefore expect that the hierarchical structure revealed when studying systems with N = 1024 should of course still be physically relevant (for instance to account for memory and rejuvenation effects), but should also be more difficult to visualise.
To illustrate these ideas, we randomly select one glass sample with N = 16384 and repeat the measurements presented in Fig. 6 above. Specifically, we quench this sample from Z g = 31 to Z = 600 using N c = 100 clones and analyse the results after a waiting time t w = 2 18 × 100. The heat map of distances between clones in Fig. 15 (a) appears featureless, which is different from the heat map constructed for the same parameters with N = 1024 in Fig. 6 (a). Accordingly, the distribution of distances between clones shown in Fig. 15 (b) present no remarkable structure, in particular, it centers around a value between the two peaks of the distribution of the N = 1024 sample in Fig. 6 (b). When we randomly select two clones to visualise the distance field in real space, as shown in Fig. 15 (c), one can observe large corre- lated domains, but the largest domains are clearly smaller than the linear size of the system, L = 120. Therefore, the large system behaves as the superposition of many smaller parts, which explains the simple forms of the heat map and the distribution.
The self-averaging suggested by these measurements can be directly confirmed by comparing the MSD measured after a quench for a single clone (symbols), to the results obtained after averaging over N c clones (lines) in Fig. 15 (d). Remarkably, the behaviour of one clone is similar to the averaged behaviour, and the individual events that could be clearly observed for N = 1024 are now washed out by the self-averaging taking place inside large systems.

VII. CONCLUSION
We show that 2d hard disks exhibit all the phenomenology expected for systems undergoing a Gardner transition inside the glass phase [7], even though no such phase transition is expected to occur in the thermodynamic limit [10][11][12]. We observe, in particular, that hard disks display slow and aging dynamics, spatially correlated motion, non-Gaussian global fluctuations, and large sample-to-sample fluctuations, which are all reminiscent of generic observations performed in systems possessing a spin glass phase at low temperatures [8,13,14,30,43]. In that sense, 2d hard disks are more similar to the d = 3 and d = ∞ hard sphere models [4,19,24] than, for in-stance, nearly 1d hard disks and soft spheres for which several important indicators of a Gardner transition are absent [20,32,37]. We show that these observations are natural, given that a large correlation is measured as the pressure increases, estimated to be ξ ≈ 30. Thus, we conclude that there is a sharp Gardner crossover in 2d hard disk glasses.
By carefully analysing the behaviour of individual glass samples, we demonstrate that the Gardner crossover is associated with the emergence of a complex free energy landscape, where the metastable glass basin breaks into smaller sub-basins with increasing pressure. We also discovered a hierarchical organisation of the landscape, associated to a hierarchy of time scales and length scales. Our study suggests that the analysis of a single glass sample with a modest system size does yield strong evidence of a phase transition [22], but our analysis of sample-tosample fluctuations, of a range of different system sizes and of different thermal histories has shown that no phase transition actually exists in 2d hard disks.
This complex organisation of the landscape contrasts strongly with the behaviour of stable glasses formed using pair potentials with no cutoff, where no jamming transition is present [32,37]. Therefore our findings suggest that systems approaching the jamming transition have a landscape that is more complex than ordinary amorphous solids. At the theoretical level, this complexity is signaled in the mean-field limit by the concept of a Gardner phase that is entered before reaching the jamming critical point [17]. In this approach the marginality of the Gardner phase is central to account for the properties of the jamming transition itself [4,33]. In physical dimensions, there is no strong evidence that a Gardner phase exists in d = 3, and no transition is expected for d < 3 [10][11][12]. However, our simulations in d = 2 suggest that the phenomenology associated to a Gardner phase can be present even when a real phase transition does not exist. It would be interesting to understand better why the man-field physics is so strongly relevant in d = 2. One possibility is that long-ranged elastic interactions between very dense hard sphere packings make the system closer to its mean-field limit.
Our results suggest, therefore, that the concepts and tools associated to the Gardner transition are not only useful to account for the jamming point, but also describe qualitatively new physical features characterizing the physical properties of a large variety of amorphous materials even away from jamming. We expect, therefore, that aging, slow dynamics, correlated particle motion, and the existence of a hierarchy of timescales, length scales and barriers to be relevant for future numerical and experimental studies of broad range of amorphous solids.
We suggest, in particular, that dense amorphous materials made of simple colloidal and granular particles should display very complex dynamic features not only when approaching the glass transition but also very deep in the arrested phase.
Where Θ is step function, (X α a,i , Y α a,i ) = R α a,i is the coordinate of particle i in clone a and sample α, and (X, Y ) = R is the center of each site.
[1] Marc Mézard, Giorgio Parisi, and Miguel Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, Vol. 9 (World Scientific