Growing timescales and lengthscales characterizing vibrations of amorphous solids

Introduction -

Understanding the nature of the glass transition, which describes the gradual transformation of a viscous liquid into an amorphous solid, remains an open challenge in condensed matter physics . As a result, the glass phase itself is not well understood either. The main challenge is to connect the localised, or ‘caged’, dynamics that characterizes the glass transition to the low-temperature anomalies that distinguish amorphous solids from their crystalline counterparts . Recent theoretical advances, building on the random first-order transition approach , have led to an exact mathematical description of both the glass transition and the amorphous phases of hard spheres in the mean-field limit of infinite-dimensional space . A surprising outcome has been the discovery of a novel phase transition inside the amorphous phase, separating the localised states produced at the glass transition from their inherent structures. This Gardner transition , which marks the emergence of a fractal hierarchy of marginally stable glass states, can be viewed as a glass transition deep within a glass, at which vibrational motion dramatically slows down and becomes spatially correlated . Although these theoretical findings promise to explain and unify the emergence of low-temperature anomalies in amorphous solids, the gap remains wide between mean-field calculations and experimental work. Here, we provide direct numerical evidence that vibrational motion in a simple three-dimensional glass-former becomes anomalous at a well-defined location inside the glass phase. In particular, we report the rapid growth of a relaxation time related to cooperative vibrations, a non-trivial change in the probability distribution function of a global order parameter, and the rapid growth of a correlation length. We also relate these findings to observed anomalies in low-temperature laboratory glasses. These results provide key support for a universal understanding of the anomalies of glassy materials, as resulting from the diverging length and time scales associated with the criticality of the Gardner transition.

Preparation of glass states –

Experimentally, glasses are obtained by a slow thermal or compression annealing, the rate of which determines the location of the glass transition . We find that a detailed numerical analysis of the Gardner transition requires the preparation of extremely well-relaxed glasses (corresponding to structural relaxation timescales challenging to simulate) in order to study vibrational motion inside the glass without interference from particle diffusion. We thus combine a very simple glass-forming model – a polydisperse mixture of hard spheres – to an efficient Monte-Carlo scheme to obtain equilibrium configurations at unprecedentedly high densities, i.e., deep in the supercooled regime. The optimized swap Monte-Carlo algorithm , which combines standard local Monte-Carlo moves with attempts at exchanging pairs of particle diameters, indeed enhances thermalization by several orders of magnitude. Configurations contain either $N=1000$ or $N=8000$ (results in Figs. 1-3 are for $N=1000$, and for $N=8000$ in Fig. 4) hard spheres with equal unit mass $m$ and diameters independently drawn from a probability distribution $P_{\sigma}(\sigma)\sim\sigma^{-3}$, for $\sigma_{\rm min}\leq\sigma\leq\sigma_{\rm min}/0.45$. We similarly study a two-dimensional bidisperse model glass former and report the main results in the Appendix.

We mimic slow annealing in two steps (Fig. 1). First, we produce equilibrated liquid configurations at various densities $\varphi_{\rm g}$ using our efficient simulation scheme, concurrently obtaining the liquid equation of state (EOS). The liquid EOS for the reduced pressure $p=\beta P/\rho$, where $\rho$ is the number density, $\beta$ is the inverse temperature, and $P$ is the system pressure, is described by

where $s_{k}$ is the $k$-th moment of $P_{\sigma}(\sigma)$, and $f(\varphi)=0.005-\tanh[14(\varphi-0.79)]$ are fitted quantities. The structure of the equilibrium configurations generated by the swap algorithm has been carefully analyzed. Unlike for other glass formers , no signs of orientational or crystalline order were observed . Following the strategy of Ref. , we also obtain the dynamical crossover $\varphi_{\rm d}=0.594(1)$ (see the Appendix). We have not analyzed the compression of equilibrium configurations with $\varphi_{\rm g}<\varphi_{\rm d}$, as done in earlier studies , because structural relaxation is not well decoupled from vibrational dynamics, although the obtained jammed states should have equivalent properties.

where the constant $C$ weakly depends on $\varphi_{\rm g}$.

Our numerical protocol is analogous to varying the cooling rate – and thus the glass transition temperature – of thermal glasses, and then further annealing the resulting amorphous solid. Each value of $\varphi_{\rm g}$ indeed selects a different glass, ranging from the onset of sluggish liquid dynamics around the mode-coupling theory dynamical crossover , $\varphi_{\rm d}$, to the very dense liquid regime where diffusion and vibrations ($\beta$-relaxation processes) are fully separated . For sufficiently large $\varphi_{\rm g}$, we thus obtain unimpeded access to the only remaining glass dynamics, i.e., $\beta$-relaxation processes .

Growing timescales –

These effects suggest a complex vibrational dynamics. Aging, in particular, provides a striking signature of a growing timescale associated with vibrations, revealing the existence of a “glass transition” deep within the glass phase.

To determine the timescale associated with this slowdown, we estimate the distance between independent pairs of configurations by first compressing two independent copies, $A$ and $B$, from the same initial state at $\varphi_{\rm g}$ to the target $\varphi$, and then measuring their relative distance

Global fluctuations of the order parameter –

The evolution of the probability distribution functions, $P(\Delta_{AB})$ and $P(\Delta)$, as well as their first moments, $\langle\Delta_{AB}\rangle$ and $\langle\Delta\rangle$, are presented in Figs. 3a,b for a range of densities across $\varphi_{\rm G}$. For $\varphi<\varphi_{\rm G}$, dynamics is fast, $\langle\Delta_{AB}\rangle$ and $\langle\Delta\rangle$ coincide, and $P(\Delta_{AB})$ and $P(\Delta)$ are narrow and Gaussian-like. For $\varphi>\varphi_{\rm G}$, however, the MSD does not converge to its long-time limit, $\langle\Delta\rangle<\langle\Delta_{AB}\rangle$, which indicates that configuration space explored by vibrational motion is now broken into mutually inaccessible regions. Interestingly, the slight increase of $\langle\Delta_{AB}\rangle$ with $\varphi$ in this regime (Fig. 3b) suggests that states are then pushed further apart in phase space, which is consistent with theoretical predictions . When compressing a system across $\varphi_{\rm G}$, its dynamics explores only a restricted part of phase space. As a result, $\Delta_{AB}$ displays pronounced, non-Gaussian fluctuations (Fig. 3a). Repeated compressions from a same initial state at $\varphi_{\rm g}$ may end up in distinct states, which explains why $\Delta_{AB}$ is typically much larger and more broadly fluctuating than $\Delta$ (Fig. 3a). These results are essentially consistent with theoretical predictions , which suggest that for $\varphi>\varphi_{\rm G}$, $P(\Delta_{AB})$ should separate into two peaks connected by a wide continuous band with the left-hand peak continuing the single peak of $P(\Delta)$. The very broad distribution of $\Delta_{AB}$ further suggests that spatial correlations develop as $\varphi\rightarrow\varphi_{\rm G}$, yielding strongly correlated states at larger densities.

Growing correlation length –

The rapid growth of $\chi_{AB}$ in the vicinity of $\varphi_{\rm G}$ suggests the concomitant growth of a spatial correlation length, $\xi$. Its measurement requires spatial resolution of the fluctuations of $\Delta_{AB}$, hence for each particle $i$ we define $u_{i}=\frac{|\boldsymbol{r}_{i}^{A}-\boldsymbol{r}_{i}^{B}|^{2}}{{\left\langle{{\Delta}_{AB}}\right\rangle}}-1$ to capture its contribution to deviations around the average $\langle\Delta_{AB}\rangle$. A first glimpse of these spatial fluctuations is offered by snapshots of the $u_{i}$ field (Fig. 1), which appear featureless for $\varphi<\varphi_{\rm G}$, but highly structured and spatially correlated for $\varphi\gtrsim\varphi_{\rm G}$. More quantitatively, we define the spatial correlator,

where $\boldsymbol{r}_{i,\mu}$ is the projection of the particle position along direction $\mu$. Even for the larger system size considered, measuring $G_{\rm L}(r)$ is challenging because spatial correlations quickly become long ranged as $\varphi\rightarrow\varphi_{\rm G}$ (see Fig. 4a). Fitting the results to an empirical form that takes into account the periodic boundary conditions in a system of linear size $L$,

where $a$ and $b$ are fitting parameters, nonetheless confirms that $\xi$ grows rapidly with $\varphi$ and becomes of the order of the simulation box at $\varphi>\varphi_{\rm G}$ (Fig. 4b). Note that although probed using a dynamical observable, the spatial correlations captured by $G_{\rm L}(r)$ are conceptually distinct from the dynamical heterogeneity observed in supercooled liquids , which is transient and disappears once the diffusive regime is reached.

Experimental consequences –

The system analysed in this work is a canonical model for colloidal suspensions and granular media. Hence, experiments along the lines presented here could be performed to investigate more closely vibrational dynamics in colloidal and granular glasses, using a series of compressions to extract $\Delta$ and $\Delta_{AB}$. Experiments are also possible in molecular and polymeric glasses, for which the natural control parameter is temperature $T$ instead of density. Let us therefore rephrase our findings from this viewpoint. As the system is cooled, the supercooled liquid dynamics is arrested at the laboratory glass transition temperature $T_{\rm g}$. As the resulting glass is further cooled its phase space transforms, around a well-defined Gardner temperature $T_{\rm G}<T_{\rm g}$, from a simple state (akin to that of a crystal) into a more complex phase composed of a large number of glassy states (see the Appendix for a discussion of the phase diagram as a function of $T$).

Around $T_{\rm G}$, vibrational dynamics becomes increasingly heterogeneous (Fig. 1), slow (Fig. 2), fluctuating from realization to realization (Fig. 3), and spatially correlated (Fig. 4). The $\beta$-relaxation dynamics inside the glass thus becomes highly cooperative and ages . The fragmentation of phase space below $T_{\rm G}$ also gives rise to a complex response to mechanical perturbations in the form of plastic irreversible events, in which the system jumps from one configuration to another . This expectation stems from the theoretical prediction that the complex phase at $T<T_{\rm G}$ is marginally stable , which implies that glass states are connected by very low energy barriers, resulting in strong responses to weak perturbations .

A key prediction is that the aforementioned anomalies appear simultaneously around a $T_{\rm G}$ that is strongly dependent on the scale $T_{\rm g}$ selected by the glass preparation protocol. Annealed glasses with lower $T_{\rm g}$ are expected to present a sharper Gardner-like crossover, at an increasingly lower temperature. Numerically, we produced a substantial variation of $\varphi_{\rm g}$ by using an efficient Monte-Carlo algorithm to bypass the need for a broad range of compression rates. In experiments a similar or even larger range of $T_{\rm g}$ can be explored , using poorly annealed glasses from hyperquenching and ultrastable glasses from vapor deposition . We expect ultrastable glasses, in particular, to display strongly enhanced glass anomalies, consistent with recent experimental reports . Interestingly, a Gardner-like regime may also underlie the anomalous aging recently observed in individual proteins .

Conclusion –

Since its prediction in the mean-field limit, the Gardner transition has been regarded as a key ingredient to understand the physical properties of amorphous solids. Understanding the role of finite dimensional fluctuations is a difficult theoretical problem . Our work shows that clear signs of an apparent critical behaviour can be observed in three dimensions, at least in a finite-size system, which shows that the correlation length becomes at least comparable to the system size as $\varphi$ approaches $\varphi_{\rm G}$. Although the fate of these findings in the thermodynamic limit remains an open question, the remarkably large signature of the effect strongly suggests that the Gardner phase transition paradigm is a promising theoretical framework for a universal understanding of the anomalies of solid amorphous materials, from granular materials to glasses, foams and proteins.

All the results discussed in this Appendix have been obtained using molecular dynamics (MD) simulations, starting from the initial states produced using the swap algorithm as explained in the main text.

Appendix A Dynamical crossover density

We follow the strategy developed in Ref. to determine the location of the dynamical (mode-coupling theory – MCT) crossover $\varphi_{\rm d}$. (i) We obtain the diffusion time $\tau_{D}=\bar{\sigma}^{2}/D$, where $D$ is the long-time diffusivity and the average particle diameter, $\bar{\sigma}$, is also the unity of length. At long times, the mean-squared displacement (MSD) $\Delta(t)=\frac{1}{N}\sum_{i=1}^{N}{\left\langle{|{\boldsymbol{r}}_{i}(t)-{\boldsymbol{r}}_{i}(0)|^{2}}\right\rangle}$ is dominated by the diffusive behavior $\Delta(t)=2dDt=2d\bar{\sigma}^{2}(t/\tau_{D})$ (Fig. 5a). Note that we here ignore the dependence of $\Delta(t)$ on $t_{\rm w}$ (compared with Eq. (1) in the main text), because we are interested in equilibrium liquid states below $\varphi_{\rm d}$, where no aging is observed. (ii) We determine the structural relaxation time $\tau_{\alpha}$ by collapsing the mean-squared typical displacement (MSTD) $r^{2}_{\rm typ}(t/\tau_{\alpha})$ in the caging regime (Fig. 5b), where the typical displacement $r_{\rm typ}(t)$ is defined as $r_{\rm typ}(t)=\lim_{z\rightarrow 0}\frac{1}{N}\sum_{i=1}^{N}{\left\langle{|{\boldsymbol{r}}_{i}(t)-{\boldsymbol{r}}_{i}(0)|^{z}}\right\rangle}^{1/z}$. (iii) We find the density threshold $\varphi_{\rm SER}=0.56(1)$ for the breakdown of Stokes-Einstein relation (SER), $D\propto\eta^{-1}$, where $\eta$ is the shear viscosity. Because $\tau_{D}\propto 1/D$ and $\tau_{\alpha}\propto\eta$ in this regime, the SER can be rewritten as $\tau_{D}\sim\tau_{\alpha}$ (Fig. 5c). (iv) We fit the time $\tau_{D}$ in the SER regime ($\varphi<\varphi_{\rm SER}$) to the MCT scaling $\tau_{D}\propto|\varphi-\varphi_{\rm d}|^{-\gamma}$ (or equivalently, $D\propto|\varphi-\varphi_{\rm d}|^{\gamma}$) to extract $\varphi_{\rm d}=0.594(1)$ (Fig. 5d).

The equilibrium liquid configurations obtained from the Monte-Carlo swap algorithm are in the deeply supercooled regime $\varphi_{\rm g}>\varphi_{\rm d}$, where the structural $\alpha$-relaxation and thus diffusion are both strongly suppressed. As long as $\varphi_{\rm g}$ is sufficiently far beyond $\varphi_{\rm d}$, the MSD for $\varphi\geq\varphi_{\rm g}$ exhibits a well-defined plateau, and the diffusive regime is not observed in the MD simulation window (see Figs. 2a and 2b in the main text). To further reveal the separation between the $\alpha$- and $\beta$-relaxations, we decompress the equilibrium configuration and show that the resulting equation of state (EOS) follows the free-volume glass EOS ( Eq. (3) in the main text) up to a threshold density, at which the system melts into a liquid (Fig. 6). This behavior suggests that our compression/decompression is slower than the $\beta$-relaxation and much faster than the $\alpha$-relaxation, such that the system is kept within a glass state. If the $\alpha$-relaxation were faster than the decompression, the state would follow the liquid EOS instead of the glass EOS under decompression. Note that a similar phenomenon has been reported in simulations of ultrastable glasses .

Appendix C Particle size effects

Suppressing the $\alpha$-relaxation and diffusion is crucial to our analysis. Besides pushing $\varphi_{\rm g}$ to higher densities, we find that it is useful to filter out the contribution of smaller particles, which are usually more mobile, from the calculation of the observables. For example, the MSD of the smaller half of the particle size distribution grows faster and diffuses sooner than that of the larger half (Fig. 7). The diffusion of smaller particles,

however, vanishes as $\varphi_{\rm g}$ increases, which suggests that the effect is not essential to the underlying physics but an artifact of our choice of system. For example, Fig. 8 shows that the smaller particles have very similar aging behavior as the larger particles (compared with Fig. 2a). For this reason, $\Delta(t,t_{\rm w})$ and $\Delta_{AB}(t,t_{\rm w})$ in this work are always calculated using only the larger half of the particle distribution.

Appendix D Distribution of single particle cage sizes

It is well known that, at the jamming point in finite dimensions, not all particles are part of the mechanically rigid network. Particles that are excluded from this network rattle relatively freely within their empty pores, hence the name “rattlers”. Because these localized excitations are not included in the infinite-dimensional theory, one could wonder whether the particles destined to become rattlers at jamming might play a role in our determination of the Gardner transition at finite dimensions. We argue here that it is not the case. Indeed, as shown in previous studies (see e.g. ), the effect of rattlers becomes important only for reduced pressure $p\gtrsim 10^{4}$. The Gardner line detected in this work covers much lower reduced pressures, $30\lesssim p_{\rm G}\lesssim 500$, which allows us to ignore rattlers.

Appendix E Caging timescale

We define a timescale $\tau_{\rm cage}$ to characterize the onset of caging. The ballistic regime of MSD at different $\varphi$ is described by a master function $\Delta(t,t_{\rm w})/\Delta_{\rm m}\sim\Delta_{\rm ballistic}(t/\tau_{\rm m})$, independent of waiting time $t_{\rm w}$, where the microscopic parameters $\tau_{\rm m}$ and $\Delta_{\rm m}$ correspond to the peak of the MSD (see Fig. 10a). To remove the oscillatory peak induced by the finite system size, we introduce $\tau_{\rm cage}$ slightly larger than, but proportional to $\tau_{\rm m}$, so that $\tau_{\rm cage}$ corresponds to the beginning of the plateau. The same collapse is obtained for any $t_{\rm w}$; as a result, $\tau_{\rm cage}$ is independent of $t_{\rm w}$. Above $\varphi_{\rm G}$, $\tau_{\rm cage}$ is the time needed for relaxing the fastest vibrations. The dependence of $\tau_{\rm cage}$ on $\varphi_{\rm g}$ is summarized in Fig. 10b. Note that $\tau_{\rm cage}\sim{\cal O}(1)$ with weak variation for all considered state points.

By contrast, the mean-squared distance between two copies, $\Delta_{AB}(t)$, depends only weakly on $t$ (see Fig. 2b in the main text). Our choice of $t=\tau_{\rm cage}$ therefore does not affect substantially the value of $\Delta_{AB}\equiv\Delta_{AB}(\tau_{\rm cage})$ (Fig. 3 in the main text). Note that, $\Delta_{AB}$ basically describes the asymptotic long time behavior of $\Delta(t)$, i.e, $\Delta_{AB}\approx\Delta_{AB}(t\rightarrow\infty)\simeq\Delta(t\rightarrow\infty,t_{\rm w}\rightarrow\infty)$.

Appendix F Absence of crystallization and of thermodynamic anomalies at the Gardner density

It is quite obvious that, because the system is not diffusing away from the original liquid configuration at $\varphi_{\rm g}$ during the simulation time window, no crystallization can happen in the glass regime. Indeed, when crossing $\varphi_{\rm G}$ no sign of incipient crystallization or formation of comparable anomaly appears in the pair correlation function. Also, $\frac{d(1/p)}{d\varphi}$ is essentially constant in the glass regime; nothing special happens to this quantity at $\varphi_{\rm G}$. The crossover would thus remain invisible if we only considered the compressibility and not more sophisticated observables.

has a stronger theoretical motivation, but at the cost of requiring more fitting constants. Data is nonetheless also well described by this functional form, as we show in Fig. 11a, where the solid lines are obtained from fits to the functional form Eq. (8) and the dashed lines are the logarithmic fits from Fig. 2c in the main text.

MCT further suggests that there should be a relation $\gamma_{\rm G}=1/a_{\rm G}$ between the exponent $\gamma_{\rm G}$ obtained from the fit of Eq. (9) and the exponent $a_{\rm G}$ from the fit of Eq. (8). At this point, we could not fit the data using a constant value of $a_{\rm G}$ for all densities. Instead, the exponent decreases as $\varphi$ increases (see Fig. 11c), and becomes increasingly incompatible with $\gamma_{\rm G}$ extracted from the fit Eq. (9). Although this fact seems to be inconsistent with the mean-field theory, the same behavior of $a_{\rm G}(\varphi)$ (also quantitatively) was recently reported in a similar study in a mean-field model of hard spheres (HS) over comparable timescales . This inconsistency might thus be related to the inner technical difficulty of fitting the exponent $a$ using Eq. (8).

Appendix H Time dependence of the skewness and determination of the Gardner density

Appendix I System-size dependence of the susceptibility at the Gardner transition

From a theoretical viewpoint, whether the mean-field Gardner transition persists in finite dimensions is still under debate . In this work, we have shown the existence of a crossover (reminiscent of the mean-field Gardner transition) at two system sizes, $N=1000$ and $N=8000$. However, the proof of the existence of the Gardner transition in the thermodynamic limit would require a systematic use of finite-size scaling techniques , which is beyond the scope of this paper. Previous studies further suggest that this kind of analysis might be extremely challenging. For example, symmetry arguments suggest that the Gardner transition should be in the same universality class as the de Almeida-Thouless line in mean-field spin-glasses in a field , whose finite-dimensional persistence is still the object of active debate even after intensive numerical scrutiny . One way to test whether our data are compatible with a true phase transition is by checking that the caging susceptibility at the transition point, $\chi_{\rm G}\equiv\chi_{AB}(\varphi=\varphi_{\rm G})$, appears to divergence at $N\rightarrow\infty$. Considering that $\chi_{\rm G}$ must be finite in a finite system (as shown it Fig. 3c for $N=1000$), the divergence requires that $\chi_{\rm G}$ increases with $N$. We can see that this requirement is fulfilled in Fig. 13, which compares susceptibilities for $N=1000$ and 8000.

Appendix J Compression-Rate Dependence

In this section, we consider how our overall analysis depends on the compression rate $\gamma_{\rm g}$ used for preparing samples. In principle, a proper $\gamma_{\rm g}$ should be such that particles have sufficient time to equilibrate their vibrations but not to diffuse. In other words, the timescale associated with compression, $\tau_{\rm g}\sim 1/\gamma_{\rm g}$, should lie between the $\alpha-$ and $\beta-$relaxation times, $\tau_{\beta}<\tau_{\rm g}<\tau_{\alpha}$. For our system, we observe that when $10^{-3}\leq\gamma_{\rm g}\leq 10^{-4}$ and $\varphi<\varphi_{\rm G}$, both $\Delta(t,t_{\rm w})$ and $\Delta_{AB}(t)$ reach flat plateaus that are essentially independent of $\gamma_{\rm g}$ (see Figs. 14a and b). Thus in this range of compression rates, restricted equilibrium within a glass state is reached, while keeping the $\alpha-$relaxation sufficiently suppressed. When $\varphi>\varphi_{\rm G}$, however, $\Delta(t,t_{\rm w})$ and $\Delta_{AB}(t)$ display $\gamma_{\rm g}$-dependent aging effects consistent with a growing timescale in the Gardner phase. As a result, the order parameters $\Delta$ and $\Delta_{AB}$, which are defined at the time scale of $\tau_{\rm cage}\sim{\cal O}(1)$, slightly depend on $\gamma_{\rm g}$ when $\varphi\gtrsim\varphi_{\rm G}$ (see Fig. 14c). This mild $\gamma_{\rm g}$-dependence has nonetheless relatively little impact on our analysis of $\varphi_{\rm G}$. In particular, the location of the peak position of the caging skewness, based on which we determine the value of $\varphi_{\rm G}$, is independent of $\gamma_{\rm g}$ within the numerical accuracy (see Fig. 14d). Note that $\gamma_{\rm g}^{-1}$ plays a role akin to the waiting time $t_{\rm w}$. Varying $\gamma_{\rm g}$ is thus equivalent to varying $t_{\rm w}$ (see Fig. 2a and Fig. 14a).

Appendix K Spatial correlation functions and lengths

The susceptibility $\chi_{AB}=N\frac{{\left\langle{\Delta^{2}_{AB}}\right\rangle}-{\left\langle{\Delta_{AB}}\right\rangle}^{2}}{{\left\langle{\Delta_{AB}}\right\rangle}^{2}}$ discussed in the main text is directly associated with the unnormalized point-to-point spatial correlation function computed between two copies, $A$ and $B$,

where $r=|\boldsymbol{r}|$ and the denominator is essentially the pair-correlation function between two clones,

In a similar way, we define the normalized line-to-line spatial correlation function

where $\boldsymbol{r}_{i,\mu}$ is the projection of the particle position along the direction $\mu$. This last definition is also Eq. (6) in the main text.

The results for this second correlation length are shown in Fig. 17b. Both estimators are expected to measure the same object, that is, both $\xi_{\rm L}$ and $\xi_{\rm P}$ should be proportional to the true correlation length, $\xi$. The actual values obtained from both fits must, however, be regarded merely as indicators of the correlation growth, because the extraction of $\xi$ is rather inaccurate. The linear size of simulation box should be several times larger than the correlation length $\xi$ to obtain accurate estimations.

Appendix L Phase diagram for thermal glasses

In order to connect our $1/p-\varphi$ phase diagram for HS (Fig. 1 in the main text) with traditional presentations of thermal glass results (see, e.g, Fig. 1 in Ref. ), we present an alternate version of that phase diagram (Fig. 19). In this different representation, we plot the specific volume $1/\varphi$ as the $y$ axis and the ratio between temperature and pressure $T/P=1/(\rho p)$ as the $x$ axis. It essentially describes how the specific volume changes with the temperature, at constant pressure, in different phases. We expect this phase diagram to be qualitatively reproducible in thermal glass experiments.

Appendix M Bidisperse hard disk results and analysis

We also study a two-dimensional bidisperse model glass former , using the same approach as for HS described the main text. The system consists of an equimolar binary mixture of $N=1000$ hard disks (HD) with diameter ratio $\sigma_{1}:\sigma_{2}=1.4:1$. In this case, we do not use the swap algorithm. Equilibrium configurations are obtained by slow relaxations during MD runs, so that particles all diffuse, i.e., $\Delta(t)\geq 10\sigma_{1}^{2}$. For each $\varphi_{\rm g}$, $N_{\rm s}=100$ samples are obtained. The liquid EOS is fitted to

is the 2$d$ Carnahan-Starling (CS) form, and

is a fitted function with parameters $c_{1}=0.52$, $c_{2}=1.0$, $c_{3}=2.7$, and $c_{4}=14$. The estimated dynamical crossover is $\varphi_{\rm d}=0.790(1)$. The $d=2$ version of $\varphi_{\rm G}$ is obtained from the peak of caging skewness of big particles. The results are summarized in Table 2 and the phase diagram is reported in Fig. 20.

Appendix N Summary of numerical results

We summarize numerical values of our main results for HS in Table 1 and for HD in Table 2.

References