# Numerical investigation of the equilibrium Kauzmann transition in a two-dimensional atomistic glass

Gerhard Jung, Misaki Ozawa, Giulio Biroli, Ludovic Berthier

► To cite this version:

Gerhard Jung, Misaki Ozawa, Giulio Biroli, Ludovic Berthier. Numerical investigation of the equilibrium Kauzmann transition in a two-dimensional atomistic glass. 2026. [https://hal.science/hal-05469873v1](https://hal.science/hal-05469873v1)

HAL Id: hal-05469873

https://hal.science/hal-05469873v1

Preprint submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Numerical investigation of the equilibrium Kauzmann transition in a two-dimensional atomistic glass

Gerhard Jung,¹ Misaki Ozawa,¹ Giulio Biroli,² and Ludovic Berthier³

¹Laboratoire Interdisciplinaire de Physique (LIPhy),
Université Grenoble Alpes, 38402 Saint-Martin-d'Hères, France

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

³Gulliver, CNRS UMR 7083, ESPCI Paris, PSL Research University, 75005 Paris, France

(Dated: July 8, 2025)

Dense liquids gradually transform into non-equilibrium amorphous solids as they pass through the experimental glass transition. Experimentally, ergodicity is lost because measurements are conducted within a finite time window. More than seventy years ago, Kauzmann posed a fundamental question: If experiments could run indefinitely, would there exist a critical temperature at which an ergodicity-breaking phase transition occurs? Random first-order transitions represent the modern theoretical framework for this idea, rigorously established in the mean-field limit of high-dimensional atomistic systems and several idealized physical models. However, achieving theoretical understanding in finite dimensions is challenging, while experimental and numerical limitations on accessible timescales hinder direct observation of the putative Kauzmann transition. Here, we overcome this longstanding barrier by developing a computational strategy that combines three advanced Monte Carlo methods to access the equilibrium thermodynamic properties of a two-dimensional atomistic glass-former down to zero temperature across a range of system sizes. This enables us to directly measure thermodynamic and structural observables that provide unambiguous evidence that the system undergoes a Kauzmann transition at a temperature that vanishes in the thermodynamic limit. This transition is towards an ideal glass state characterized by a complex energy landscape with a hierarchical organization of low-lying states. Our results are the first demonstration that computer simulations can fully probe the statistical mechanics of the bulk transition to a non-ergodic glass state. We anticipate that our study will serve as a foundation for future simulation work on larger systems, three-dimensional materials, and more complex glass-forming models to fully elucidate the nature of the glass state of matter.

## INTRODUCTION

Conventional glasses are obtained when liquids fall out of equilibrium below the experimental glass transition temperature $T_g$ [1], see Fig. 1a. Below $T_g$, relaxational dynamics becomes slower than the observation time, and the liquid appears frozen in a non-ergodic amorphous state. In his landmark paper [2], Kauzmann famously articulated fundamental questions regarding the fate of the liquid state if thermodynamic equilibrium and ergodicity could be maintained far below $T_g$. While Kauzmann debated the relative stability of liquid and crystalline states [3], his data collection also hinted at a possible thermodynamic phase transition at the (now called [4, 5]) Kauzmann temperature, $T_K < T_g$, between a liquid and an ideal glass state, see Fig. 1a.

The hypothesis of an experimentally inaccessible Kauzmann transition associated with an entropy crisis has been a central theme in glass research over the last decades [1, 4–9]. From a theoretical viewpoint, Kauzmann’s key contribution was the proposal of a novel type of thermodynamic phase transition between an ergodic liquid and a non-ergodic, non-periodic solid. Developing accurate theories of the glass transition remains a fundamental open problem in classical statistical physics, with a growing body of work impacting various scientific fields (including protein folding [10], biological processes [11],

computer science [12] and machine learning [13]) where analogs of a Kauzmann transition were identified.

The transition envisioned by Kauzmann has been described in refined theoretical analysis, from the early work of Gibbs-diMarzio [6] to its modern incarnation as a random first-order transition [14, 15]. The patient elaboration of a consistent theoretical framework of the liquid-glass transition culminated in the last decade with the unification [16] of the statistical physics of liquids in large dimensions with the Parisi theory of disordered systems, first developed for spin glasses. A random first-order phase transition with a vanishing equilibrium free energy difference between liquid and glass phases (interpreted as a vanishing configurational entropy $TS_c$, see Fig. 1a) is mathematically demonstrated in large dimensions, with known features [16]. In particular, the crystal no longer play a role in this analysis.

There is a profound disconnect between the detailed understanding of the equilibrium phase transition between liquid and glass states at mean-field level, and the dramatic lack of solid results in two and three dimensions where this should apply. Several reasons explain this unsatisfactory situation. First, even modern experimental techniques fail to produce equilibrated samples down to $T_K$ [17]. Second, analytic progress beyond mean-field is difficult, as this requires generalizing Parisi’s theory to finite dimensions [18]. Third, despite notable recent progress [19–23], numerical simulations have so far failed
to approach $T_K$ in bulk systems because equilibrium sampling of the complex energy landscape of glassy liquids at low enough temperatures is a computationally hard problem [24].

We overcome for the first time the longstanding sampling challenge in a realistic off-lattice glass model. We developed an original computational strategy that integrates the complementary strengths of several advanced Monte Carlo techniques. This approach enables ergodic sampling of the configuration space of a two-dimensional glass model [25] down to $T = 0$, across a range of system sizes. This methodological breakthrough therefore allows us to fully characterize the thermodynamics of the system, including its ground state and organization of low-lying excited states. Tailored free energy measurements [26] demonstrate the existence of a random first-order transition towards a low-temperature equilibrium glass phase characterized by a hierarchical organization of low-temperature states. By increasing the system size, the transition temperature approaches zero. Our work offers the first glimpse of the equilibrium liquid-glass phase transition in finite dimensions. It should serve as a foundation for future simulation work using larger systems, three-dimensional and molecular models [27] of glass-formers to fully illuminate the nature of the glass state of matter.

## NUMERICAL STRATEGY TO ACHIEVE EQUILIBRATION AND ERGODIC SAMPLING

We developed a computational strategy to successfully sample the configuration space of a realistic glass-former across the bulk Kauzmann transition, see Fig. 1b. We study a ternary mixture of Lennard-Jones particles in two dimensions [28] (see Methods Sec. I). We recently developed this computational model to fulfill multiple constraints relevant to the present study. (i) It is microscopically realistic (it resembles metallic glasses and previously studied glass models). (ii) It is robust against ordering (demixing, crystallisation). (iii) The swap Monte Carlo (MC) algorithm [19] alone provides a significant equilibration speedup compared to conventional molecular dynamics. (iv) This model was recently used to benchmark sampling algorithms [25] and to devise the appropriate tools (based on specific heat measurements) that validate equilibration and ergodic sampling.

For reference, the onset temperature for slow dynamics is $T_o \approx 0.5$, the mode-coupling crossover (where molecular dynamics becomes inefficient) is $T_{mct} \approx 0.3$, and the experimental glass transition (where typical experiments fall out of equilibrium) is $T_g \approx 0.15$. Using swap MC, the system falls out of equilibrium near $T_{min} \approx 0.12 < T_g$. Although it performs below $T_g$, swap MC is insufficient to approach closely, or even cross, a putative thermodynamic Kauzmann singularity. We show in Fig. 1c the average inherent state energy $\langle e_{inh} \rangle$ obtained after quenching to $T = 0$ configurations obtained at temperature

$T$ [29]. These data confirm that swap MC equilibrates much better than conventional molecular dynamics, but falls out of equilibrium near $T_{min}$.

We discovered that a carefully designed combination of three Monte Carlo techniques is exceptionally efficient at sampling the rugged configuration space at low $T$. Our approach, sketched in Fig. 1b, is composed of two main steps. First, we combine swap MC to parallel tempering [30, 31] to produce large sets of $N_s = 10^6$ independent equilibrium configurations. For system sizes $N = 44, 55, 66, \text{ and } 77$, we can construct such equilibrium sets down to $T = 0.105 < T_{min}$, thus already improving over the performances of swap MC alone. At this temperature, however, the system is still in the supercooled liquid phase. The second step uses the population of $N_s$ equilibrium configurations as starting points for the population annealing algorithm [32, 33]. In this algorithm, the population is annealed to lower temperatures in a stepwise manner, with population resampling and reweighting steps performed at each $T$. Figure 1c demonstrates that population annealing yields much lower inherent state energies, suggesting that it pushes the equilibration limit to much lower temperatures.

To decide whether ergodic sampling is achieved, we follow the time evolution of the specific heat, $c_V$. A correct determination of $c_V$ requires proper ergodic exploration of the configuration space, which represents a demanding and sensitive test [25]. In Fig. 1d we show, for selected temperatures, that the specific heat measured between in the interval $0 \le T \le 0.105$ always converges to its equilibrium value. This convergence test was successful for all investigated sizes down to $T = 0$, but failed at a finite temperature $T > 0$ for larger $N$ values. This test sets the limit to the $N$ values analyzed. See also Methods Secs. III, IV and Supplemental Figs. S2, S3) for additional equilibration tests confirming these conclusions.

In short, our strategy combining swap, parallel tempering and population annealing Monte Carlo algorithms allows us to determine the equilibrium thermodynamic behavior of the system at any $T \ge 0$ for $N \le 77$. In addition, the $T \to 0$ limit offers a direct determination of the ground state and opens a window into low-lying excitations, thus allowing a detailed analysis of the structure of the potential energy landscape.

## CONFIGURATIONAL ENTROPY

In his original work [2], Kauzmann defined a configurational (or ‘excess’) entropy as the difference between the supercooled liquid and crystal entropies. Its extrapolation to temperatures below $T_g$ revealed what is now called the Kauzmann temperature, as the moment where this excess entropy may vanish. Our current understanding, however, has shifted: the configurational entropy has to be defined independently of the crystal. What is physically more relevant is the entropy difference with respect to amorphous glassy states [26, 34, 35]. Several
<figure><img src="image_2.png" /></figure>

FIG. 2. Numerical measurements of the configurational entropy. **a.** Configurational entropy for $N = 77$ calculated using various techniques: expanding around inherent states including anharmonic corrections ($s_c^{\text{harm}}$), expanding around equilibrium states using the Frenkel-Ladd method ($s_c^{\text{FL}}$), and Shannon entropy obtained from the distribution of inherent states ($s_c^{\text{inh}}$). When only using swap Monte Carlo, $s_c^{\text{harm,SMC}}$ displays a non-equilibrium saturation at low $T$ when the system falls out of equilibrium. The agreement between the three methods is remarkable, and all equilibrium estimates indicate that the configurational entropy vanishes near $T = 0$. **b.** Inherent state configurational entropy $s_c^{\text{harm}}$ for different system sizes $N$. Equilibrium is achieved down to $T = 0$ for $N = 44$ and 77. For $N = 176$ and 704, data is shown down to the lowest temperature for which equilibration is reached. The dotted line is an empirical extrapolation of the $N = 704$ data using a quadratic polynomial fit.

dependence on $N$ over a wide range. Our measurements from $N = 44$ to $N = 77$ are all consistent with a smooth evolution and an entropy crisis taking place at $T = 0$ exactly. For larger systems, $N = 176$ and $N = 704$ (where the bulk limit is nearly reached) we can no longer equilibrate the system at very low temperatures. The dashed line through the $N = 704$ data is a hypothetical extrapolation used earlier [41] to estimate the entropy at vanishing temperature, again consistent with $T = 0$.

## FREE ENERGY MEASUREMENTS AND THE KAUZMANN TRANSITION

Having located an entropy crisis at $T = 0$ by direct equilibrium measurements, we are left with two hypotheses. This behavior could suggest a zero-temperature critical point associated with a new kind of phase transition, such as a random first-order transition. Or, this behavior

<figure><img src="image_3.png" /></figure>

FIG. 3. First-order equilibrium Kauzmann transition revealed by the evolution of the Franz-Parisi free energy. **a.** The Franz-Parisi potential $\beta V(Q)$ as a function of the overlap order parameter $Q$, between independent pairs of equilibrium configurations for $N = 55$. The potential is vertically shifted arbitrarily for visibility. For $T$ far above $T_K$ the potential has a single minium at small $Q$ corresponding to the liquid phase. A secondary minimum appears at large $Q$ above $T_K$, corresponding to a metastable glass. At $T = T_K$, both minima have the same free energy. Below $T < T_K$, the large $Q$ glass minimum dominates. The evolution of the free energy is reminiscent of conventional first-order phase transitions. **b.** Temperature-dependent average overlap $\langle Q \rangle$ for different system sizes $N$ reveals a sharp change from small to large at a size-dependent transition temperature. The peculiar case of $N = 66$ is discussed in Methods Sec. VIII. **c.** Evolution of the critical Kauzmann temperature $T_K$ with system size, suggesting that $T_K \to 0$ for large $N$. This is consistent with the entropy extrapolation for $N = 704$ in Fig. 2b.

simply reveals that, at very low temperatures, the system is confined near its ground state, possibly decorated with localized excitations, as proposed by Stillinger [4]. The configurational entropy is unable to distinguish between these two situations [36], and a different analysis is needed to investigate whether finite-size systems exhibit a Kauzmann transition into an ideal glass phase, and how this putative transition evolves with $N$.

We employ the formalism first introduced by Franz and Parisi in the context of mean-field glass models [42], and extensively used in recent numerical studies of finite dimensional supercooled liquids [26, 43, 44]. We study an effective Landau free energy, $V(Q)$, expressed as a function of the overlap $Q$. The overlap quantifies the similarity of the density fields in two independent configurations (see Methods Sec. VI), and it is the correct order parameter of the liquid-glass Kauzmann transition. The function $V(Q)$, known as the Franz-Parisi potential, is the free energy cost to observe the value $Q$ of the overlap with respect to a fixed reference equilibrium configuration.
In this framework, the Kauzmann transition resembles a conventional first-order phase transition associated to a jump in the order parameter $Q$ [45]. In the liquid above $T_K$, the global minimum of $V(Q)$ is at low overlap, $Q \approx Q_{\text{liq}}$, indicative of a liquid phase in which the system accesses many different configurations with small mutual overlap. In the glass below $T_K$, $V(Q)$ displays a different global minimum at a much larger value, $Q \approx Q_{\text{glass}}$, revealing that the equilibrium glass only samples a limited number of amorphous states. In this regime, there is a finite probability to remain close to a reference configuration, thus leading to a large value of $Q$. The overlap jumps discontinuously from $Q_{\text{liq}}$ to $Q_{\text{glass}}$ at $T_K$, where the free energy difference $\Delta V = V(Q_{\text{glass}}) - V(Q_{\text{liq}})$ vanishes. In contrast to the estimates of the configurational entropy shown in Fig 2, the definition $Ts_c = \Delta V$ can lead to a genuine entropy crisis [26].

To determine $V(Q)$, we measure the mutual overlap $Q_{mn}$ between two independent equilibrium configurations $m$ and $n$. Given two independent sets of $N_s$ equilibrium structures at temperature $T$ we evaluate the quenched Franz-Parisi potential as:

$$V(Q) = -\frac{T}{N} \sum_{m=1}^{N_s} \ln \left( \sum_{n=1}^{N_s} \delta(Q - Q_{mn}) \right). \quad (1)$$

Unlike previous approaches, which relied on biased Monte Carlo techniques to probe rare fluctuations of $Q$ much above $T_K$ [43], we directly construct $V(Q)$ from unbiased equilibrium sampling at any $T$. This is possible because the system is equilibrated near $T_K$, where the relevant fluctuations of $Q$ become typical.

Measuring the Franz-Parisi potential down to $T = 0$ reveals our main physical result, see Fig. 3. A Kauzmann transition from the liquid (with a small average overlap) to the ideal glass phase (with a large average overlap) takes place when the glass free energy becomes lower than that of the liquid, as directly observed in Fig. 3a. At temperatures much larger than $T_K$, the system exhibits low-overlap configurations with $Q \approx Q_{\text{liq}}$, characteristic of the liquid state. As the temperature decreases, a second local minimum appears at high overlap $Q_{\text{glass}}$, becoming degenerate with the liquid minimum at a well-defined temperature that we identify as the Kauzmann temperature: $T_K = 0.057$ for $N = 55$. Below $T_K$, the system enters in the low-temperature equilibrium glass, characterized by the dominance of the high-overlap minimum. Correspondingly, the average overlap jumps abruptly from low to large value in a narrow temperature regime near $T_K$, as shown in Fig. 3b. We insist that these low-temperature free energy and overlap measurements at $T < T_K$ are all performed in equilibrium conditions.

Phase transitions are rounded in finite systems and the free energy $V(Q)$ not necessarily convex [44]. To understand finite size corrections, we studied four system sizes. The evolution of both the average overlap and of $V(Q)$ indicates a systematic shift of $T_K$ to lower temperatures with increasing $N$, see Figs. 3b,c. These data show that

$T_K$ rapidly shifts towards zero as $N$ increases. Although the accessible system sizes do not permit a quantitative finite-size scaling analysis, they indicate that $T_K < 0.01$ in the large-$N$ limit, which is substantially lower than the glass transition temperature $T_g = 0.15$ and conspicuously close to $T_K = 0$. Overall, these findings provide compelling evidence for the existence of a zero-temperature critical point at $T_K = 0$ in the thermodynamic limit.

## AMORPHOUS GROUND STATES AND HIERARCHY OF LOW-LYING EXCITATIONS

We now study the lowest energy configurations to investigate whether the Kauzmann transition is associated with a rich and complex structure of the free energy landscape, as theoretically expected. For all investigated $N$, we obtain global ground states that are completely amorphous, as illustrated in Fig. 4f for $N = 44$ and 77. This suggests that crystallization is efficiently suppressed. We cannot exclude that a periodic structure with a very large unit cell eventually emerges for much larger $N$.

To explore the organization of low-energy minima, we classify all identified inherent states into families based on structural similarity between them. The bimodality of the Franz-Parisi potential indicates that pairs of configurations are either very similar or markedly distinct. We define a family as a set of inherent states with mutual overlap $Q > 0.8$. This threshold value was chosen based on the results shown in Fig. 3.

Each inherent state is characterized by two excitation energies: the excitation energy relative to the global ground state, $\Delta E_{\text{inh}} = E_{\text{inh}} - E_{\text{gs}}$, and an intra-family excitation energy, $\Delta E_{\text{intra}} = E_{\text{inh}} - E_{\text{fam}}$, where $E_{\text{fam}}$ denotes the minimum energy within the same family. We can also define the position of a given family relative to the ground state, $\Delta E_{\text{fam}} = E_{\text{fam}} - E_{\text{gs}}$. By construction, one has $\Delta E_{\text{inh}} = \Delta E_{\text{intra}} + \Delta E_{\text{fam}}$. This decomposition is useful, as we can disentangle intra-state (localised) excitations from inter-states (collective) excitations.

We focused our analysis on low-lying minima, specifically those with excitation energies $\Delta E_{\text{inh}} \ll |E_{\text{gs}}|/N \approx 3.0$, as illustrated in Fig. 4c. The probability distributions of inherent states as a function of the excitation energies defined above are shown in Fig. 4c, and their decompositions into localized and collective excitations in Figs. 4d,e for three system sizes. We observe a proliferation of low-lying states as $N$ increases, with a nearly equivalent contribution from both types of excitation energies. Thus, the bottom of the potential energy landscape is composed of several distinct families, themselves decorated by localized excitations.

To visualize the organization of states and localized excitations, we compute the matrix of mutual overlap between low-lying states [46] alongside the corresponding clustering trees in Figs. 4a,b. In the tree representation, each circle represents an inherent state, with the vertical position reflecting its inherent energy, while connections
<figure><img src="image_4.png" /><figcaption>FIG. 4. Hierarchical organization of low energy states. a. Tree and matrix representation of the bottom of the potential energy landscape for $N = 44$. States with mutual overlaps $Q > 0.8$ are grouped into families. The family minimum is shown with a green circle and the ground state with a yellow circle, excited states as blue circles. In the tree, the x-axis is arbitrary, the y-axis is the inherent energy. Using the same ordering of states, the matrix represents the value of the mutual overlap and reveals the hierarchical organisation of the bottom of the energy landscape. b. Same as in a for $N = 77$. c. Probability of finding states with inherent energy $\Delta E_{\text{inh}}$ for various system sizes $N = 44, 55$, and 77. d. Probability of inter-family excitation energies $\Delta E_{\text{fam}}$ for the same $N$ values. e. Probability of intra-family excitation energies $\Delta E_{\text{intra}}$ for the same $N$ values. The proliferation of both distinct minima and localized excitations for increasing $N$ is obvious. f. Snapshot of the ground states for $N = 44$ and $N = 77$. The color encodes the particle type. g. Excerpt of the $N = 77$ tree, and corresponding snapshots of the selected excited states showing various degree of overlap $Q$ with the ground state. The color encodes the local overlap for each particle using the same color code as in a.</figcaption></figure>

in the tree are based on mutual overlaps and constructed via a clustering algorithm (see Methods Sec. VII). The resulting structure reveals a clear hierarchy: inherent states are grouped into families, and families in turn cluster into higher-order families of families. This organization reveals a complex, hierarchical structure at the bottom of the energy landscape, even after localized excitations have been coarse-grained.

Finally, a real-space view of the hierarchical phase-space organization is provided in Fig. 4g. We provide one snapshot of a localized excitation within the ground-state family, two snapshots of collective excitations from distinct families with a finite overlap with the ground state, and one snapshot from a structurally distinct family. These snapshots reveal that intra-family excitations correspond to localized rearrangements, while transitions between different families involve increasingly collective up to system-spanning reorganizations. These observations are in agreement with the bimodal overlap distribution obtained from the Franz-Parisi free energy potential. Once again all these structures are fully amorphous.

DISCUSSION

Because we can directly access all thermodynamic quantities down to arbitrarily low temperatures, our results clarify the behavior of glass-forming materials when equilibrated at temperatures that have remained so far been inaccessible to both experiments and simulations.

First, our analysis demonstrates that the existence of the crystal phase can be, and should be when possible, separated from the discussion of the Kauzmann transition. Second, we find that localized excitations, invoked by Stillinger [4, 47], indeed proliferate with system size. This is expected since their energy cost is not extensive. However, the Franz-Parisi free energy naturally coarse-grains these excitations and allows us to probe the properties of collective states. The finite-size Kauzmann transition described in Fig. 3 characterizes the temperature evolution of those states, and is unaffected by local defects.

Our findings provide compelling evidence that the landscape of low-lying states is intrinsically complex and organized in a hierarchical manner, as envisioned in Parisi's theory of disordered systems [16, 46]. We uncover a very similar picture with clusters of low-energy configurations which are nested within one another, forming a multi-level hierarchy. The bottom of the energy land-
scape is populated by minima whose energy differences are much smaller than the energy per particle. Yet, these minima correspond to completely different configurations in real space, with very small mutual overlap, thus reflecting global, system-wide rearrangements.

The natural interpretation is that a non-trivial Kauzmann transition occurs at $T_K = 0$ in our system. While this has been suggested using extrapolations before [40, 41], this conclusion is now supported by direct equilibrium measurements down to $T = 0$. Testing whether the same picture holds for more complex

models and three-dimensional systems is an exciting and important challenge, for which our work lays the conceptual and methodological foundation. This conclusion should also motivate theoretical work to incorporate finite-dimensional fluctuations into the framework of random first order transitions, as current attempts and analogies to random field Ising models [48, 49] do not point towards the possibility of a zero-temperature critical point in two-dimensional glasses. We hope that our solid numerical results will provide useful guidance and motivation for future work.

[1] M. D. Ediger, C. A. Angell, and S. R. Nagel, Supercooled liquids and glasses, *The Journal of Physical Chemistry* **100**, 13200 (1996).

[2] W. Kauzmann, The nature of the glassy state and the behavior of liquids at low temperatures, *Chemical Reviews* **43**, 219 (1948).

[3] H. Tanaka, Possible resolution of the kauzmann paradox in supercooled liquids, *Phys. Rev. E* **68**, 011505 (2003).

[4] F. H. Stillinger, Supercooled liquids, glass transitions, and the kauzmann paradox, *The Journal of Chemical Physics* **88**, 7818 (1988).

[5] R. J. Speedy, Kauzmann's paradox and the glass transition, *Biophysical chemistry* **105**, 411 (2003).

[6] J. H. Gibbs and E. A. DiMarzio, Nature of the glass transition and the glassy state, *The Journal of Chemical Physics* **28**, 373 (1958).

[7] P. G. Debenedetti and F. H. Stillinger, Supercooled liquids and the glass transition, *Nature* **410**, 259 (2001).

[8] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, *Rev. Mod. Phys.* **83**, 587 (2011).

[9] R. S. Welch, E. D. Zanotto, C. J. Wilkinson, D. R. Casar, M. Montazerian, and J. C. Mauro, Cracking the kauzmann paradox, *Acta Materialia* **254**, 118994 (2023).

[10] J. D. Bryngelson and P. G. Wolynes, Spin glasses and the statistical mechanics of protein folding, *Proceedings of the National Academy of Sciences* **84**, 7524 (1987).

[11] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, *Phys. Rev. X* **6**, 021011 (2016).

[12] R. Monasson, R. Zecchina, S. Kirkpatrick, B. Selman, and L. Troyansky, Determining computational complexity from characteristic 'phase transitions', *Nature* **400**, 133 (1999).

[13] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to machine learning for physicists, *Physics Reports* **810**, 1 (2019).

[14] V. Lubchenko and P. G. Wolynes, Theory of structural glasses and supercooled liquids, *Annu. Rev. Phys. Chem.* **58**, 235 (2007).

[15] G. Biroli and J.-P. Bouchaud, The random first-order transition theory of glasses: A critical assessment, in *Structural Glasses and Supercooled Liquids* (John Wiley and Sons, Ltd, 2012) Chap. 2, pp. 31–113.

[16] G. Parisi, P. Urbani, and F. Zamponi, *Theory of simple glasses: exact solutions in infinite dimensions* (Cambridge University Press, 2020).

[17] M. D. Ediger, Perspective: Highly stable vapor-deposited glasses, *The Journal of Chemical Physics* **147**, 210901 (2017).

[18] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, Glass and jamming transitions: From exact results to finite-dimensional descriptions, *Annual Review of Condensed Matter Physics* **8**, 265 (2017).

[19] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation of glass transition studies, *Phys. Rev. X* **7**, 021039 (2017).

[20] V. Bolton-Lum, R. C. Dennis, P. Morse, and E. Corwin, The ideal glass and the ideal disk packing in two dimensions, *arXiv preprint arXiv:2404.07492* (2024).

[21] X. Fan, D. Xu, J. Zhang, H. Hu, P. Tan, N. Xu, H. Tanaka, and H. Tong, Ideal noncrystals: A possible new class of ordered matter without apparent broken symmetry, *arXiv preprint arXiv:2404.17675* (2024).

[22] F. Leoni, J. Russo, F. Sciortino, and T. Yanagishima, Generating ultrastable glasses by homogenizing the local virial stress, *Phys. Rev. Lett.* **134**, 128201 (2025).

[23] F. Ghimenti, L. Berthier, and F. van Wijland, Irreversible monte carlo algorithms for hard disk glasses: From event-chain to collective swaps, *Phys. Rev. Lett.* **133**, 028202 (2024).

[24] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, *Nature Reviews Physics* **5**, 102 (2023).

[25] G. Jung, G. Biroli, and L. Berthier, Normalizing flows as an enhanced sampling method for atomistic supercooled liquids, *Machine Learning: Science and Technology* **5**, 035053 (2024).

[26] L. Berthier and D. Coslovich, Novel approach to numerical measurements of the configurational entropy in supercooled liquids, *Proceedings of the National Academy of Sciences* **111**, 11668 (2014).

[27] R. Simon, J.-L. Barrat, and L. Berthier, Molecular motion at the experimental glass transition, *arXiv preprint arXiv:2506.21290* (2025).

[28] G. Jung, G. Biroli, and L. Berthier, Predicting dynamic heterogeneity in glass-forming liquids by physics-inspired machine learning, *Phys. Rev. Lett.* **130**, 238202 (2023).

[29] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid, *Nature* **393**, 554 (1998).

[30] R. H. Swendsen and J.-S. Wang, Replica monte carlo simulation of spin-glasses, *Phys. Rev. Lett.* **57**, 2607 (1986).
[31] K. Hukushima and K. Nemoto, Exchange monte carlo method and application to spin glass simulations, Journal of the Physical Society of Japan 65, 1604 (1996).

[32] K. Hukushima and Y. Iba, Population annealing and its application to a spin glass, in *AIP Conference Proceedings*, Vol. 690 (American Institute of Physics, 2003) pp. 200–206.

[33] J. Machta, Population annealing with weighted averages: A monte carlo method for rough free-energy landscapes, Phys. Rev. E 82, 026704 (2010).

[34] F. Sciortino, W. Kob, and P. Tartaglia, Inherent structure entropy of supercooled liquids, Phys. Rev. Lett. 83, 3214 (1999).

[35] S. Sastry, The relationship between fragility, configurational entropy and the potential energy landscape of glass-forming liquids, Nature 409, 164 (2001).

[36] L. Berthier, M. Ozawa, and C. Scalliet, Configurational entropy of glass-forming liquids, The Journal of Chemical Physics 150, 160902 (2019).

[37] F. Sciortino, Potential energy landscape description of supercooled liquids and glasses, Journal of Statistical Mechanics: Theory and Experiment 2005, P05015 (2005).

[38] M. Ozawa, G. Parisi, and L. Berthier, Configurational entropy of polydisperse supercooled liquids, The Journal of Chemical Physics 149, 154501 (2018).

[39] D. Frenkel and A. J. C. Ladd, New monte carlo method to compute the free energy of arbitrary solids. application to the fcc and hcp phases of hard spheres, The Journal of Chemical Physics 81, 3188 (1984).

[40] A. D. Parmar and A. Heuer, Depleting states dictate the ideal glass and physics of glass transition, arXiv preprint arXiv:2307.10143 (2023).

[41] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Zero-temperature glass transition in two dimensions, Nature Communications 10, 1508 (2019).

[42] S. Franz and G. Parisi, Phase diagram of coupled glassy systems: A mean-field study, Phys. Rev. Lett. 79, 2486 (1997).

[43] L. Berthier, Overlap fluctuations in glass-forming liquids, Phys. Rev. E 88, 022313 (2013).

[44] B. Guiselin, L. Berthier, and G. Tarjus, Statistical mechanics of coupled supercooled liquids in finite dimensions, SciPost Phys. 12, 091 (2022).

[45] B. Guiselin, G. Tarjus, and L. Berthier, Is glass a state of matter?, Physics and Chemistry of Glasses-European Journal of Glass Science and Technology Part B 63, 136 (2022).

[46] L. Berthier, G. Biroli, P. Charbonneau, E. I. Corwin, S. Franz, and F. Zamponi, Gardner physics in amorphous solids and beyond, The Journal of Chemical Physics 151, 010901 (2019).

[47] S. Sastry, Numerical test of stillinger's argument concerning the kauzmann paradox, The Journal of Physical Chemistry B 108, 19698 (2004).

[48] J. D. Stevenson, A. M. Walczak, R. W. Hall, and P. G. Wolynes, Constructing explicit magnetic analogies for the dynamics of glass forming liquids, The Journal of Chemical Physics 129, 194505 (2008).

[49] G. Biroli, C. Cammarota, G. Tarjus, and M. Tarzia, Random field ising-like effective theory of the glass transition. ii. finite-dimensional models, Physical Review B 98, 174206 (2018).

[50] A. D. S. Parmar, M. Ozawa, and L. Berthier, Ultrastable metallic glasses in silico, Phys. Rev. Lett. 125, 085505 (2020).

[51] H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, Unveiling dimensionality dependence of glassy dynamics: 2d infinite fluctuation eclipses inherent structural relaxation, Phys. Rev. Lett. 117, 245701 (2016).

[52] S. Büchner and A. Heuer, Potential energy landscape of a model glass former: Thermodynamics, anharmonicities, and finite size effects, Phys. Rev. E 60, 6507 (1999).

[53] K. Kim and R. Yamamoto, Apparent finite-size effects in the dynamics of supercooled liquids, Phys. Rev. E 61, R41 (2000).

[54] L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli, Finite-size effects in the dynamics of glass-forming liquids, Phys. Rev. E 86, 031502 (2012).

[55] D. Gazzillo and G. Pastore, Equation of state for symmetric non-additive hard-sphere fluids: An approximate analytic expression and new monte carlo results, Chemical Physics Letters 159, 388 (1989).

[56] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019).

[57] D. J. Evans and B. L. Holian, The nose-hoover thermostat, The Journal of Chemical Physics 83, 4069 (1985).

[58] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271, 108171 (2022).

[59] R. Yamamoto and W. Kob, Replica-exchange molecular dynamics simulation for supercooled liquids, Phys. Rev. E 61, 5473 (2000).

[60] S. T. Tokdar and R. E. Kass, Importance sampling: a review, Wiley Interdisciplinary Reviews: Computational Statistics 2, 54 (2010).

[61] D. Gessert, W. Janke, and M. Weigel, Resampling schemes in population annealing: Numerical and theoretical results, Phys. Rev. E 108, 065309 (2023).

[62] B. Coluzzi, G. Parisi, and P. Verrocchio, Lennard-jones binary mixture: A thermodynamical approach to glass transition, The Journal of Chemical Physics 112, 2933 (2000).

[63] J. M. Polson, E. Trizac, S. Pronk, and D. Frenkel, Finite-size corrections to the free energies of crystalline solids, The Journal of Chemical Physics 112, 5339 (2000).

[64] B. Coluzzi and G. Parisi, On the approach to the equilibrium and the equilibrium properties of a glass-forming model, Journal of Physics A: Mathematical and General 31, 4349 (1998).

[65] H. W. Kuhn, The hungarian method for the assignment problem, Naval research logistics quarterly 2, 83 (1955).

[66] W. Kob and L. Berthier, Probing a liquid to glass transition in equilibrium, Phys. Rev. Lett. 110, 245702 (2013).

[67] E. Schubert and L. Lenssen, Fast k-medoids clustering in rust and python, Journal of Open Source Software 7, 4183 (2022).
# METHODS

## I. TWO-DIMENSIONAL ATOMISTIC GLASS-FORMER

The atomistic glass-former simulated in this work is very similar to the model introduced in Ref. [28]. The model corresponds to a modified Kob-Andersen mixture in two-dimensions (KA2D) in which particles interact via a Lennard-Jones potential,

$$V_{\alpha\beta}(r_{ij}) = \begin{cases} 4\epsilon_{\alpha\beta} \left[ \left( \frac{\sigma_{\alpha\beta}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{\alpha\beta}}{r_{ij}} \right)^{6} + C_0 \right. \\ + C_2 \left( \frac{r_{ij}}{\sigma_{\alpha\beta}} \right)^2 + C_4 \left( \frac{r_{ij}}{\sigma_{\alpha\beta}} \right)^4 \\ \left. 0 \right] & r_{ij} < r_{\alpha\beta}^{\text{cut}} \\ 0 & \text{otherwise,} \end{cases}$$

for particle distance $r_{ij} = |r_i - r_j|$, where $r_i$ is the position of the i-th particle. The ternary KA2D model consists of three types $\alpha, \beta = \{1, 2, 3\}$ where types 1 and 2 interact via the usual Kob-Andersen non-additive interactions, $\epsilon_{11} = 1.0$, $\epsilon_{12} = 1.5$, $\epsilon_{22} = 0.5$ and $\sigma_{11} = 1.0$, $\sigma_{12} = 0.8$, $\sigma_{22} = 0.88$. Additionally, we introduce a third species [50] with interaction parameters $\epsilon_{13} = 0.75$, $\epsilon_{23} = 1.5$, $\epsilon_{33} = 0.75$ and $\sigma_{13} = 0.9$, $\sigma_{23} = 0.8$, $\sigma_{33} = 0.94$. The cutoff depends on the particle type $r_{\alpha\beta}^{\text{cut}} = 2.5\sigma_{\alpha\beta}$ and we define the constants $C_0 = 0.04049023795$, $C_2 = -0.00970155098$ and $C_4 = 0.00062012616$ to make the potential continuous up to the second derivative. Differently from Ref. [28] we slightly change the type composition to a ratio 5:3:3. The system parameters were empirically adjusted to avoid any signatures of local crystallization, and indeed no crystalline structure was found in any of the simulation runs performed. Results are reported in reduced Lennard-Jones units defined by $\epsilon_{11}$ (energy scale), $\sigma_{11}$ (length scale) and $\sigma_{11}/\sqrt{m/\epsilon_{11}}$ (time scale), with mass $m = 1$ for all types.

The particles evolve in a square box of dimensions $L_x = L_y$ with constant particles density $\rho = N/(L_x L_y) = 1.19$, where N is the number of particles. We use N = 44, 55, 66, 77, 176, and 704 in this study. The linear box sizes are therefore $L_x(N = 44) = 6.08$, $L_x(N = 55) = 6.8$, $L_x(N = 66) = 7.44$, $L_x(N = 77) = 8.036$, $L_x(N = 176) = 12.16$ and $L_x(N = 704) = 24.32$.

The relaxation dynamics of the KA2D model using MD dynamics is shown in Supplemental Fig. S1. We find the expected relaxation dynamics as featured by the incoherent scattering function (ISF) for both N = 704 (Fig. S1a) and N = 77 (Fig. S1b). We note that the effects of Mermin-Wagner fluctuations [51] are small and negligible at these system sizes [28]. In particular, the ISF shows a decay towards a plateau, highlighting vibrational dynamics in cages followed by a long-time structural relaxation. Importantly, the N = 77 system behaves dynamically as a typical bulk supercooled liquid with negligible finite size effects [52]. As expected, the relaxation times for N = 77 are slightly higher than for N = 704 due to the small box size [53, 54], which includes additional constraints (compare Figs. S1c and d) but both systems show

an apparent Arrhenius behavior, $\tau_{\alpha}^{\text{ISF}}(T) \propto \exp(6.1/T)$, for small temperature. The glass transition temperature $T_g$ can be identified as, $\tau_{\alpha}^{\text{ISF}}(T_g) = 10^{12}$, yielding $T_g = T_g(N = 704) = 0.147$ and $T_g(N = 77) = 0.156$. We will report the former value throughout this work since it is the value we expect in the thermodynamic limit, but the difference is not very important, considering that we have mainly focused on temperatures $T < 0.1$ in this manuscript.

## II. INHERENT STATE ENERGY

Given a set of $N_c$ equilibrium configurations, we can extract ensemble averages. In particular, we calculate the ensemble-averaged inherent energy density, $\langle e_{\text{inh}} \rangle = \langle E_{\text{inh}} \rangle/N$. For each configuration $m$ with particle coordinates $\mathbf{R}_m = (\mathbf{r}_{m,1}, \mathbf{r}_{m,2}, \dots, \mathbf{r}_{m,N})$ we determine the nearest inherent state, $\mathbf{R}_{\text{inh},m}$. This is achieved by performing steepest decent with a maximal step size per particle $\Delta r = 0.01$ to minimize the potential energy $E_{\text{pot}}$. Finally, we define $E_{\text{inh}} = \sum'_{i,j} V(r_{\text{inh},ij})/2$, as the potential energy of the energy minimum.

## III. CRITERIA FOR EQUILIBRATION AND SAMPLING

Based on a time-dependent specific heat $c_V(t)$ we have recently established a reliable and independent measure to test whether an ensemble of supercooled liquid configurations represents an equilibrium ensemble or not [25]. Here, the 'time' $t = M\tau$ corresponds to the number $M$ of configurations used for averaging, multiplied by the relaxation time $\tau$ used between the production of two independent configurations. To define this measure we introduce the average, $\langle X \rangle_M = M^{-1} \sum_{m=1}^M X_m$ of an observable $X_m$ evaluated for configuration $m$. Using this definition, we obtain the potential energy, $\langle E_{\text{pot}} \rangle_M$ and specific heat at temperature $T$, $c_V(t = M\tau) = (\langle E_{\text{pot}} \rangle_M - \langle E_{\text{pot}} \rangle_M^2)/(NT^2)$.

If a set of $N_c$ configurations represents an equilibrium ensemble, we expect $c_V(t)$ to increase from $c_V(t) = 1$, corresponding, by equipartition in a two-dimensional system, to thermal vibrations close to an arrested structure, and reach a long-time plateau $c_V = c_V(t \to \infty)$ that corresponds to the equilibrium value of the specific heat. In the case that $c_V(t = N_c\tau)$ has not yet reached a plateau, this clearly indicates that the set is not sufficiently equilibrated [25]. For N = 77, which was by far the most difficult system size to equilibrate (by comparison with smaller N values) we observe from Fig. 1d in the main text, that $c_V(t)$ attains a plateau at all investigated temperatures. This provides compelling evidence that we can create equilibrium ensembles down to very low temperatures. This test fails whenever we do not use the combination of Monte Carlo techniques described in the main text.
# IV. NUMERICAL TECHNIQUES ENABLING EQUILIBRATION DOWN TO ZERO TEMPERATURE

Equilibration is achieved by the unique combination of three different Monte-Carlo techniques: swap Monte-Carlo (SMC), parallel tempering (PT) and population annealing (PA).

## A. Swap Monte Carlo algorithm

SMC is based on the idea of accelerating equilibration by swapping the positions of pairs of particles of different types [19, 55, 56]. The model introduced in the previous section was constructed to enable the best usage of SMC via the insertion of the intermediate type, which massively increases the swap probability [50].

First, we create configurations with the composition defined above. Then, we use SMC in combination with molecular dynamics (MD) simulations for short thermalization runs of approximately $N_{t,th} = 10^8$ time steps at $T = 0.148$. The molecular dynamics time step is set to $\Delta t = 0.005$ and the system is thermalized using a Nosé-Hoover thermostat [57]. Every 10 MD steps we include a phase of SMC moves. In each phase we perform $2N$ swap attempts, by randomly choosing two particles of different type, swapping their position, and accepting the attempt via the Metropolis acceptance criterion [56]. All simulations are performed using LAMMPS [58].

## B. Parallel tempering

The thermalized configurations at $T = 0.148$ by SMC are then used as input for PT simulations. The main idea behind PT is to perform parallelized SMC simulations of several replicas of the system [31]. Each replica evolves at a different temperature. Regularly, exchange events are attempted in which two replicas exchange their current configuration, which is accepted via the Metropolis acceptance criterion. If successful, configurations will travel in temperature space by visiting different replicas [30, 59]. Configurations originally simulated at a low temperature can thus randomly be exchanged between replicas, relax at high temperature, and possibly travel back to low temperature. This procedure has been shown to accelerate sampling in many instances.

We perform simulations of 8 replicas in parallel (9 for $N = 77$) at temperatures $T = 0.256, 0.229, 0.205, 0.183, 0.164, 0.148, 0.134$, and $0.122$ (and $T = 0.111$ for $N = 77$). Exchange events are attempted every $10^3$ MD steps. The microscopic dynamics is equivalent to the one described above in the SMC section. We first perform equilibration runs for $N_{t,eq} = 2 \times 10^9$ steps, and afterwards sample the equilibrium ensemble for $N_{t,prod} = 10^{10}$ steps, by extracting $10^4$ configurations equally distributed over time. This procedure amounts

to around $9 \times 180$ CPU hours for $N = 77$ (nine cores for roughly one week). We have performed 2 (for $N = 44, 55$), 24 (for $N = 66$), and 78 (for $N = 77$) independent PT runs.

In Fig. S2 we show an analysis of the equilibrium ensembles created by this sampling procedure. For any of the sampled temperatures, the time-dependent specific heat attains a well-defined plateau, clearly indicating equilibration (see Fig. S2a,b). We additionally observe that two completely independent sets of configurations produce the same ground state and low energy excitations, similar to the observation in Ref. [40]. We additionally establish the more demanding criterion that the states are sampled with the same probability (see Fig. S2c-f). These ergodic sets will be the starting point of the population annealing procedure.

## C. Population annealing

Population annealing is based on using an initial set of configurations, representing an equilibrium ensemble at temperature $T_0$ and perform successive reweighting and population resampling to create ensembles at lower temperatures [32, 33, 60, 61].

We use the $N_c$ configurations sampled using the PT simulations described above as input for PA in an attempt to reach even lower temperatures. In population annealing step $n$, we choose a temperature $T_n < T_{n-1}$ slightly smaller than $T_{n-1}$ and assign the Boltzmann weight $W_m = \exp(-(\beta_n - \beta_{n-1})E_{pot,m})$ to each configuration $m$ in the population. Here, $E_{pot,m}$ is the total potential energy of configuration $m$ and $\beta_n = (k_B T_n)^{-1}$ the inverse temperature with $k_B = 1$. Afterwards, we create $\tau_m = N_c W_m / \sum_{s=1}^{N_c} W_s$ copies of configuration $m$. To maintain an ensemble of constant size $N_c$, we use the “systematic resampling” scheme of Ref. [61]. We then use these annealed configurations and perform short thermalization runs of $5 \times 10^3$ MD steps at temperature $T_n$, using the SMC dynamics.

The PA procedure starts at $T_0 = 0.122$ ($T_0 = 0.111$ for $N = 77$) and is applied down to $T_{nmax} = 0.001$ where $nmax = 25$. In other words, we perform 25 PA annealing steps to reach the lowest temperature investigated in this work. Temperature steps are slightly non-uniformly distributed with larger steps in the beginning, but the results were insensitive to this choice, as long as a sufficient number of small steps was used.

In Fig. 1d of the main manuscript we have shown that the time-dependent specific heat $c_V(t = N_c\tau)$ attains a plateau, indicating equilibration. Here, $\tau = 5 \times 10^3$ corresponds to the relaxation time between the output of two separate configurations in the PT run, which were then used as input for the PA procedure. In Supplementary Fig. S3a we observe the same for the time-dependent potential energy $\langle E_{pot} \rangle_{N_c}$. We found that the specific heat is a superior measure for equilibration compared to the potential energy since its equilibration time is much
longer, consistent with the findings in Ref. [25].

For $T = 0.057$ we perform a reweighting (RW) procedure, which corresponds to a single population annealing step, starting from $T_0 = 0.111$ quenching directly to $T_1 = 0.057$. We observe that using RW only, the equilibration time is significantly increased, and a plateau is not reached when studying $\langle E_{\text{pot}} \rangle_{N_c}$ within our simulation time window. The reason for the improved performance of PA compared to RW can be inferred from the histogram of inherent states, as shown in Supplementary Fig. S3b,c. PA systematically leads to an increased probability of finding configurations with low inherent energies. As described above, this is ensured by the stepwise selection of configurations with low thermal energy and resampling their thermal energies. In contrast, when solely performing a single reweighting step, thermal fluctuations are blurring the selection of the proper low-temperature states, and reweighting is unable to correctly sample low temperatures. The above analysis highlights the importance of using a stepwise PA procedure to obtain the results presented in this manuscript.

More generally, we have only managed to create equilibrium ensembles at $T = 0.001$ when using the combination of the three described Monte Carlo techniques. All three of them play a key role in this success, and each one of them is unable to achieve equilibration by itself. In Ref. [25] we have intensively benchmarked enhanced sampling techniques. While PA alone performs well, we still found that it fell out of equilibrium at around $T_{\text{PA}} = 0.19$, and thus at temperature much higher than the glass transition temperature $T_g$. Similarly, without the usage of PT, we were not able to equilibrate the system with $N = 77$. As a conclusion, combining SMC/PT/PA in the order described above is vital to ensure equilibration down to $T \to 0$.

## V. CONFIGURATIONAL ENTROPY

We have used three different techniques to evaluate the configurational entropy.

### A. Expansion near inherent states

The idea of this method is to evaluate the vibrational entropy $s_{\text{glass}}$ by first performing a harmonic approximation around each inherent state. Afterwards, an harmonic correction is applied to improve the results. Combined with the measurement of the total entropy, this procedure then allows to calculate the configurational entropy,

$$s_c = s_{\text{tot}} - s_{\text{glass}}. \quad (S1)$$

Our methodology mainly follows the lines of Ref. [36], but we had to derive several finite size corrections which are relevant due to the small system sizes investigated in the present work, compared to earlier efforts.

We start by calculating the total entropy $s_{\text{tot}}$. The total Helmholtz free energy can be obtained by thermodynamic integration starting from infinite temperature $\beta = 0$:

$$\beta F_{\text{tot}} = \beta F_{\text{id}} + \int_{0}^{\beta} d\beta' \langle E_{\text{pot}}(\beta') \rangle, \quad (S2)$$

with the ideal gas free energy,

$$\beta F_{\text{id}} = -\ln Z_{\text{id}} - Ns_{\text{mix}}, \quad (S3)$$

derived from the partition sum,

$$Z_{\text{id}} = \frac{1}{N!} \left( \frac{L^2 m}{2\pi \beta \hbar^2} \right)^{dN/2}, \quad (S4)$$

the mixing entropy,

$$Ns_{\text{mix}} = \ln \left( \frac{N!}{\prod_{\alpha} N_{\alpha}!} \right), \quad (S5)$$

and the average potential energy $\langle E_{\text{pot}}(\beta) \rangle$. Here, $N$ is the total number of particles, $N_\alpha$ denotes the number of particles of type $\alpha$, $d=2$ is the space dimension, and $L = L_x = L_y$ the linear box length. Using the relation $s_{\text{tot}} = \frac{d}{2} + N^{-1}\beta \langle E_{\text{pot}}(\beta) \rangle - N^{-1}\beta F_{\text{tot}}$ we find the entropy per particle,

$$s_{\text{tot}} = \frac{d}{2} + \frac{\ln Z_{\text{id}} + \beta \langle E_{\text{pot}}(\beta) \rangle}{N} + s_{\text{mix}} - \int_{0}^{\beta} d\beta' \langle E_{\text{pot}}(\beta') \rangle. \quad (S6)$$

Note that using the standard terms for $s_{\text{mix}}$ based on Stirling's approximation leads to a difference $\Delta s_{\text{tot}}(N = 44) \approx 0.15$ compared to the above expressions. Therefore, using the full formula in Eq. (S5) is essential, in particular for systems with $N \lesssim 100$.

In practice, we separate the integral,

$$\int_{0}^{\beta} d\beta' \langle E_{\text{pot}} \rangle = \int_{0}^{\beta_0} d\beta' \langle E_{\text{pot}} \rangle + \int_{\beta_0}^{\beta} d\beta' \langle E_{\text{pot}} \rangle \\ = I_{\text{F}} + I_{\text{N}}. \quad (S7)$$

For $I_{\text{N}}$, we use standard numerical integration schemes, for example Simpson integration, from $\beta_0$ to $\beta$. For $I_{\text{F}}$, we fit the temperature-dependence of the potential energy based on the high-temperature expansion [62]:

$$\beta \langle E_{\text{pot}}(\beta) \rangle = A\beta^x + B(\beta^x)^2 + C(\beta^x)^3 + \dots, \quad (S8)$$

where $A$, $B$, and $C$ are coefficients, and $x = d/n$ with $n$ being the exponent of the repulsive power-law potential ($x = 1/6$ in our case). By using this expansion, we can calculate $I_{\text{F}}$ via

$$I_{\text{F}} = \frac{1}{x} A \beta_0^x + \frac{1}{2x} B (\beta_0^x)^2 + \frac{1}{3x} C (\beta_0^x)^3 + \dots. \quad (S9)$$

The different contributions to the total entropy are shown in Supplemental Fig. S4a. We observe that even
for $N = 77$ the finite-size corrections by not using Stirling's approximation account to $\Delta s_{tot}(N = 77) > 0.03$ which is a very relevant contribution, especially at low temperatures where we seek to study the vanishing of the configurational entropy. Supplemental Fig. S4a also highlights that the non-trivial temperature dependence of $s_{tot}$ emerges from the energetic contributions.

Having calculated the total entropy, we now evaluate the glass entropy, $s_{glass}$, which denotes the entropy contribution from thermal vibrations. The harmonic approximation is based on expanding the potential energy around its inherent state, $E_{pot} \approx E_{inh} + \frac{1}{2}\Delta R^T H \Delta R$, where $\Delta R = R - R_{inh}$ and $H$ is the Hessian matrix. This enables the calculation of the partition sum,

$$Z_{\text{harm}} = \Lambda^{-Nd} e^{-\beta E_{\text{inh}}} \int d\mathbf{R} e^{-\frac{1}{2} \beta \Delta \mathbf{R}^T H \Delta \mathbf{R}}, \quad (\text{S10})$$

where $\Lambda = \sqrt{2\pi\hbar^2\beta}$ is the thermal de Broglie wavelength. We do not have the mixing term because each inherent state exists $\prod_{\alpha} N_{\alpha}$! times. $H$ is diagonalized by an orthogonal matrix $U$. We then obtain the eigenvalues of $H$ (sorted by ascending order), $0 = \lambda_1 = \lambda_2 = \cdots = \lambda_d < \lambda_{d+1} \le \lambda_{d+2} \le \cdots \le \lambda_{Nd}$. The first $d$ eigenvalues are zero due to translational invariance. Using a variable transformation, $\mathbf{X} = U\Delta\mathbf{R}$, one gets

$$ 
\begin{aligned}
Z_{\text{harm}} &= \Lambda^{-Nd} e^{-\beta E_{\text{inh}}} \prod_{s=1}^{Nd} \int X_s e^{-\frac{\beta}{2} \lambda_s X_s^2} \\
&= \Lambda^{-Nd} e^{-\beta E_{\text{inh}}} \left( \prod_{s=1}^{d} \int X_s \right) \prod_{s=d+1}^{Nd} \int X_s e^{-\frac{\beta}{2} \lambda_s X_s^2} \\
&= \Lambda^{-Nd} e^{-\beta E_{\text{inh}}} N^{\frac{d}{2}} L^d \prod_{s=d+1}^{Nd} \sqrt{\frac{2\pi}{\beta\lambda_s}}.
\end{aligned} 
\quad (\text{S11})
$$

In the last equation we used $X_s = \sqrt{N}r_s^{\text{CM}}$, where $r_s^{\text{CM}}$ is the $s$-th component of the center of mass, due to the orthonormality of the eigenvectors associated with the zero modes. This gives rise to $\int dX_s = \sqrt{N}L$.

Finally, we calculate the glass entropy, $s_{glass}^{\text{harm}} = \frac{d}{2} + N^{-1}\beta(\langle E_{\text{pot}} \rangle - \overline{E_{\text{harm}}})$ by subtracting the free energy $\beta F_{\text{harm}} = -\ln Z_{\text{harm}}$ from the kinetic and potential energies. Here, $(\cdots)$ denotes the average over independent inherent states. Consistent with previous work, we assume within the harmonic approximation $\langle E_{\text{pot}} \rangle \approx \overline{E_{\text{inh}}} + \frac{(N-1)d}{2}T$. Collecting all the terms, we find

$$s_{\text{glass}}^{\text{harm}} = \frac{1}{N} \left[ \sum_{s=d+1}^{Nd} \left( 1 - \overline{\ln(\beta\hbar\omega_s)} \right) + \frac{d}{2} \left( 1 - \ln\left(\frac{2\pi\hbar^2\beta}{NL^2}\right) \right) \right] \quad (\text{S12})$$

where $\omega_s = \sqrt{\lambda_s}$. This term is identical to the typical expression (as used for instance in Ref. [36]), including finite-size corrections of order $O(N^{-1} \ln N)$.

The above calculations were performed within the harmonic approximation and we have neglected anharmonic

contributions, $\langle E_{\text{anh}} \rangle = \langle E_{\text{pot}} \rangle - \left( \overline{E_{\text{inh}}} + \frac{(N-1)d}{2}T \right)$. To improve upon the harmonic approximation, we fit $\langle E_{\text{anh}} \rangle$ using polynomials of the form $\langle E_{\text{anh}}(T) \rangle / N = \sum_{k=2} a_k T^k$. This allows us to evaluate the anharmonic contribution,

$$s^{\text{anh}} = N^{-1} \int_{0}^{T} dT' \frac{1}{T'} \frac{\partial \langle E_{\text{anh}}(T') \rangle}{\partial T'}, \quad (\text{S13})$$

which yields, $s_c^{\text{anh}} = \sum_{k=2} \frac{k}{k-1} a_k T^{k-1}$. Finally, we collect all terms to obtain the configurational entropy:

$$s_c^{\text{harm}} = s_{\text{tot}} - s_{\text{glass}}^{\text{harm}} - s^{\text{anh}}. \quad (\text{S14})$$

The different contributions to $s_c^{\text{harm}}$ are shown in Supplemental Fig. S4b. Importantly, we observe that the vanishing of the configurational entropy at $T = 0$ is far from trivial as it results from a cancellation of three large terms in Eq. (S14). In particular, $s_{tot}$ is obtained from a thermodynamic integration starting at $T = \infty$, while $s_{glass}^{\text{harm}} + s^{\text{anh}}$ is based on the harmonic approximation around the $T \to 0$ inherent state. Both contributions thus emerge from opposite and independent limits. The observation that $s_c^{\text{harm}}(T \to 0) \to 0$ is thus a strong indication that we have correctly sampled the equilibrium configurations at low temperatures, and that the proposed estimate of the configurational entropy is physically meaningful.

## B. Expansion near equilibrium configurations: Frenkel-Ladd method

A major disadvantages of the harmonic approach described above is that it relies on the concept of inherent states. While well-established, it is desirable to have an independent measure of configurational entropy. This can be achieved by applying a method close to the one proposed by Frenkel-Ladd (FL) [39] to compute the free energy of solids. While being physically similar to the harmonic approach, the FL method calculates vibrational contributions via a thermodynamic integration from an Einstein solid without the need to determine inherent states. In this approach, the configurational entropy is

$$s_c^{\text{FL}} = s_{\text{tot}} - s_{\text{glass}}^{\text{FL}}, \quad (\text{S15})$$

where $s_{tot}$ is the total entropy derived above. The FL glass entropy is determined from the Helmholtz free energy via thermodynamic integration from a system in which each particle in a reference configuration is harmonically trapped. The interaction energy becomes $\beta U_{\alpha}(\mathbf{R}, \mathbf{R}_0) = \beta E_{\text{pot}}(\mathbf{R}) + \alpha \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{r}_{0,i})^2$. Here, $\mathbf{R}_0 = (\mathbf{r}_{0,1}, \mathbf{r}_{0,2}, \dots, \mathbf{r}_{0,N})$ denotes the positions of the particles in a given reference equilibrium configuration. The corresponding free energy and partition function are
given by

$$ \beta F_{\alpha} = -\ln Z_{\alpha}, \qquad (S16) $$

$$ Z_{\alpha} = \Lambda^{-Nd} \int \mathbf{R}e^{-\beta U_{\alpha}(\mathbf{R}, \mathbf{R}_{0})}. \qquad (S17) $$

Special care is required when computing the free energy for small system sizes using the Frenkel-Ladd method. The thermodynamic integration requires accurate evaluation of the mean-squared displacement from the reference configuration under a harmonic potential, as explained below. The Monte Carlo estimation of this quantity must be performed with the center of mass (CM) motion removed. Otherwise, the integration diverges when the spring constants become small. However, our goal is to estimate the glass entropy including the degrees of freedom associated with the center of mass, since this entropy must be compared with the total entropy, which also accounts for the CM contribution. For sufficiently large systems, the distinction between CM-fixed and CM-free treatments is negligible. In our case, however, where the system size is small, this distinction has a non-negligible impact on the precise determination of the configurational entropy. Therefore, it is essential to first relate the numerically accessible free energy of the CM-fixed system to that of the CM-free system, which is what we wish to obtain [63].

In practice, we rewrite $\beta F_{\alpha}$ as

$$ \beta F_{\alpha} = -\ln Z_{\alpha} = -\ln Z_{\alpha}^{\text{CM}} + \ln \frac{Z_{\alpha}^{\text{CM}}}{Z_{\alpha}}, \qquad (S18) $$

where

$$ Z_{\alpha}^{\text{CM}} = \Lambda^{-Nd} \int \mathbf{R}e^{-\beta U_{\alpha}(\mathbf{R}, \mathbf{R}_{0})} \delta(\mathbf{r}^{\text{CM}}). \qquad (S19) $$

Here $\delta(\mathbf{r}^{\text{CM}})$ constrains the position of the center of mass, $\mathbf{r}^{\text{CM}} = N^{-1} \sum_{i=1}^{N} \mathbf{r}_i$, at the origin $\mathbf{0}$ which is set by the center of mass in the reference configuration: $N^{-1} \sum_{i=1}^{N} \mathbf{r}_{0,i} = \mathbf{0}$. The last term in Eq. (S18) corresponds to the probability to find $\mathbf{r}^{\text{CM}}$ at the origin in the CM-free ensemble, which is nothing but $L^{-d}$:

$$ \frac{Z_{\alpha}^{\text{CM}}}{Z_{\alpha}} = \langle \delta(\mathbf{r}^{\text{CM}}) \rangle_{\alpha} = \frac{1}{L^{d}}. \qquad (S20) $$

This term, along with the delta function in Eq. (S19), represents a finite-size correction to the standard approach.

When $\alpha_{\text{max}}$ is sufficiently large, the system corresponds to an Einstein solid and we have

$$ \beta F_{\text{glass}}^{\text{FL}} = \beta \overline{E_{\text{pot}}(\mathbf{R}_0)} - \ln Z_{\text{Ein}} - \lim_{\alpha_{\text{min}} \to 0} N \int_{\alpha_{\text{min}}}^{\alpha_{\text{max}}} d\alpha \Delta(\alpha) + \ln \langle \delta(\mathbf{r}_{\text{CM}}) \rangle, \qquad (S21) $$

with the partition sum,

$$ Z_{\text{Ein}} = \Lambda^{-Nd} \int d\mathbf{R} e^{-\alpha_{\text{max}} \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{r}_{0,i})^2} \delta(\mathbf{r}_{\text{CM}}), \qquad (S22) $$

and the average displacement from the reference configuration, $\Delta(\alpha) = N^{-1} \langle \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{r}_{0,i})^2 \rangle_{\alpha}^{\text{CM}}$, where $\langle \cdots \rangle$ denotes averaging over reference configurations. We fix the CM motion during the Monte-Carlo evaluation. $Z_{\text{Ein}}$ can be solved analytically using the inverse Fourier transform of the delta function.

$$ Z_{\text{Ein}} = \Lambda^{-Nd} \int \frac{d\mathbf{k}}{(2\pi)^d} \prod_{i=1}^{N} \int d\Delta \mathbf{r}_i e^{-\alpha_{\text{max}} \Delta \mathbf{r}_i^2 - \frac{i}{N} \mathbf{k}^T \Delta \mathbf{r}_i}, \qquad (S23) $$

where $\Delta \mathbf{r}_i = \mathbf{r}_i - \mathbf{r}_{0,i}$. Subsequently, we perform the Gaussian integrals to get

$$ Z_{\text{Ein}} = \Lambda^{-Nd} \left( \frac{\pi}{\alpha_{\text{max}}} \right)^{Nd/2} \left( \frac{\alpha_{\text{max}} N}{\pi} \right)^{d/2} . \qquad (S24) $$

This term is identical to the finite N-correction derived in Ref. [63]. We finally get the glass entropy:

$$ s_{\text{glass}}^{\text{FL}} = \frac{d}{2} - \frac{d}{2} \left(1 - \frac{1}{N}\right) \ln\left(\frac{\alpha_{\text{max}}}{\pi}\right) - \ln \Lambda^d + \frac{d}{2} \frac{\ln N}{N} + \frac{\ln L^d}{N} + \lim_{\alpha_{\text{min}} \to 0} \int_{\alpha_{\text{min}}}^{\alpha_{\text{max}}} d\alpha \Delta(\alpha). \qquad (S25) $$

We evaluate $\Delta(\alpha)$ using Monte Carlo simulations and integrate the results to obtain the final term in Eq. (S25). We also take the size dispersity into account following the methodology derived in Ref. [38]. In practice, we compute $\Delta(\alpha)$ using swap Monte Carlo and estimate the additional mixing entropy contribution with Monte Carlo sampling that employs only particle-swap moves with no translational displacements. In our system, the influence of polydispersity is negligible, in particular at lower temperatures, so it does not affect much our results.

In Supplemental Fig. S4 we compare the glass entropy calculated from the Frenkel-Ladd method, $s_{\text{glass}}^{\text{FL}}$, the harmonic approximation evaluated above and find a very good agreement. This shows that inherent states are not required to define a configurational entropy.

### C. Enumeration of inherent states

At very low temperatures, the equilibrium ensemble of thermalized configurations populates a small number of inherent states. This enables us to directly enumerate them [40]. Since we have direct access to the probability $p_m(T)$ to observe inherent state $m$ at temperature $T$, we can estimate the configurational entropy via the Shannon entropy of the probability distribution,

$$ s_c^{\text{inh}}(T) = - \sum_m p_m(T) \ln p_m(T). \qquad (S26) $$

At higher temperatures we expect this measure to deviate from the true configurational entropy, since we will not be able to enumerate all inherent states anymore. Therefore, we have not shown data for $s_c^{\text{inh}}(T)$ above $T = 0.075$ in Fig. 2a.
## VI. OVERLAP AND FRANZ-PARISI POTENTIAL

The overlap $Q$ between two independently sampled configurations $\mathbf{R}_1$ and $\mathbf{R}_2$ with $N$ particles is calculated as [64],

$$Q(\mathbf{D}) = \frac{1}{N} \sum_{i=1}^{N} F(D_i) \qquad (S27)$$

with the minimal distance between configurations $\mathbf{D} = (D_1, D_2, ..., D_N)$. Here, $\mathbf{D}$ is the distance vector between two configurations, $\mathbf{D} = (|\mathbf{r}_{1,1} - \mathbf{r}_{2,1}|, |\mathbf{r}_{1,2} - \mathbf{r}_{2,2}|, ..., |\mathbf{r}_{1,N} - \mathbf{r}_{2,N}|)$, where $\mathbf{r}_{m,i}$ is the position of particle $i$ in configuration $m$. The function $F(x) \in [0, 1]$ should be unity for very small distances and quickly decays to zero for larger distances. We maximize the overlap $Q$ over possible translations $\mathcal{T}$, rotations $\mathcal{R}$, and mirror transformations $\mathcal{M}$ for the configuration $\mathbf{R}_2$. Importantly, labeling of the particles in configuration $\mathbf{R}_2$ is defined such that it minimizes the total distance, $\mathbf{D}^T \mathbf{D}$, via the Hungarian algorithm [65], which gives rise to $\mathbf{D}_{\min}$. We therefore find

$$\max_{\mathcal{T}, \mathcal{R}, \mathcal{M}} Q(\mathbf{D}_{\min}). \qquad (S28)$$

This expression defines the quantity analyzed throughout this manuscript. We have chosen the specific form,

$$F(x) = \frac{1}{2} \left( e^{-\left(\frac{x}{a_0(T)}\right)^2} + e^{-\left(\frac{x}{a_1}\right)^2} \right), \qquad (S29)$$

where $a_1 = 0.1$ and $a_0(T)$ corresponds to the square root of the plateau height of the mean-squared displacement at temperature $T$. None of the results shown in this manuscript depend qualitatively on this specific choice of $F(x)$, but both terms together ensure a stable position of the potential minima observed in Fig. 3a for the various temperatures by accounting for decreasing thermal fluctuations as $T \to 0$. Without the first term, at very low temperatures, all overlaps would be $Q = 1$, and without the second term the small-$Q$ peak shifted strongly towards $Q = 0$.

We maximize the overlap by iterating over various transformations. In two dimensions, we need to use three rotations $(\pi/2, \pi$ and $3\pi/2)$ and four mirror symmetries $(x = 0, y = 0, x = y$ and $x = -y)$. We find the optimal translation by making steps of $\Delta r = 0.5\sigma_{11}$. This implies that we need to verify $(2L/\sigma_{11})^2$ copies of configuration $\mathbf{R}_2$. For each of these copies we apply the Hungarian algorithm to find the minimal distance vector $\mathbf{D}$ between the copy and configuration $\mathbf{R}_1$. Subsequently, we subtract the average distance, i.e., $D_i \to D_i - N^{-1} \sum_j D_j$. This step explains in hindsight why it is sufficient to validate translations in coarse steps $\Delta r$: if two configurations strongly overlap, all particles will have the identical distance, which will then be removed in the subsequent correction step. Using the corrected distance vector, we

calculate $Q$ and finally choose the maximum after iterating over all copies.

The optimization procedure is illustrated in Supplemental Fig. S5a. It shows how two seemingly very different configurations are actually the same after mirroring the $x$-coordinates and shifting. The above algorithm iteratively validates various possibilities and determines the optimal transformations and shift.

We extract the overlap $Q_{mn}$ by selecting $N_s = 10^3$ independent configurations $\mathbf{R}_m$ and, similarly, $N_s$ independent configurations $\mathbf{R}_n$. It has been ensured that $\mathbf{R}_m$ and $\mathbf{R}_n$ originate from two independent SMC+PT+PA runs. We can thus exclude the possibility that large overlaps stem from insufficient sampling or slow decorrelation within a single run.

Finally, we calculate the quenched Franz-Parisi potential,

$$\beta V(Q) = -\frac{1}{N} \sum_{m=1}^{N_s} \ln \left( \sum_{n=1}^{N_s} \delta(Q - Q_{mn}) \right), \qquad (S30)$$

where the $\delta$-function is discretized with a bin size $\Delta Q = 0.01$. The term inside the logarithm is the histogram of overlap values sampled for a given reference configuration $m$, and the first sum averages over independent $m$ configurations.

In addition to the Franz-Parisi potential shown for $N = 55$ in Fig. 3a, we show the ones extracted for $N = 44$ and $N = 77$ in Supplemental Figs. S5b,c, respectively. While all systems show a clear transition at a critical temperature $T_K$ (which depends on $N$) the potentials look quantitatively slightly different. This is mainly caused by the local minima observed in Supplemental Fig. S5b,c. These local minima correspond to configurations coming from related families with overlaps $Q$ which are between $Q_{liq}$ and $Q_{glass}$ (see Fig. 4 and discussion). For $N = 55$ only one family contributed at lower temperatures, giving the Franz-Parisi potential a simpler shape without local minima.

Transitions such as shown in Fig. 3 have only been observed for structural glass-forming liquids using biased Monte-Carlo simulations [44] or systems with pinned particles [66]. Here, we find it using a large set of independent equilibrium structures in the bulk.

## VII. CONSTRUCTING FAMILIES OF LOW-ENERGY STATES

Using Eq. (S28) we calculate the overlap between all the low-energy states analyzed in Fig. 4 of the main manuscript. The ground state defines the first family (yellow sphere in Fig. 4a). Subsequently, we iterate over all low-energy states, sorted by energy and starting from the ground state. If the state has an overlap $Q < 0.8$ with all previously identified family members, a new family is created and the current state is identified as the family
minimum (purple spheres in Fig. 4a). If it has an overlap $Q > 0.8$ with one of the family minima it is instead identified as an excited state and included into an existing family, chosen from the state with which it has the largest overlap.

Each state is characterized by its inherent energy difference, $\Delta E_{\text{inh}} = E_{\text{inh}} - E_{\text{gs}}$, where $E_{\text{gs}}$ is the ground state energy (shown in Fig. 4b). Each family has an excitation energy $\Delta E_{\text{fam}} = E_{\text{fam}} - E_{\text{gs}}$, where $E_{\text{fam}}$ is the energy of the family minimum (shown in Fig. 4c). Finally, each state also has an excitation energy $\Delta E_{\text{intra}} = E_{\text{inh}} - E_{\text{fam}}$, denoting its energy relative to the family minimum (shown in Fig. 4d). As a consequence, the energy difference to the ground state is formed by adding the two individual excitation energies, $\Delta E_{\text{inh}} = \Delta E_{\text{fam}} + \Delta E_{\text{intra}}$.

To create the tree-based and hierarchical representation in Figs. 4a,b we use the k-Medoids clustering algorithm provided by the python library FasterPAM [67]. This clustering algorithm is usable for non-Euclidean metrics such as the overlap $Q$. Different from typical clustering algorithms, the non-Euclidean nature of $Q$ requires the most central point within the cluster to be an actual object, also called the medoid. The total number of clusters (and thus medoids) is set to equal the number of families found using the procedure described above. In the k-Medoids algorithm the medoids are iteratively adapted to minimize the total distance of each object to the nearest medoid. We then order the states within one family according to their inherent energies. Subsequently, families themselves are clustered if the states belonging to each of the families have an average overlap $Q > 0.25$. The final order of states is thus iteratively determined by starting with the state with lowest inherent energy (i.e., the ground state in the first iteration), followed by the states belonging to the same family. Then we include all those states that belong to families which are similar to the selected family. This procedure is iterated until all states have been selected.

## VIII. THE SPECIAL CASE OF 66 PARTICLES

In Fig. 3b it can be seen that the average overlap $\langle Q \rangle$ for $N = 66$ behaves differently from that of the other system sizes. This can be explained by investigating the distribution of inherent energies and the corresponding ground state and excited states (see Supplemental Fig. S6). For this system size, we observe that at $T = 0.035$ a majority of the configurations have an inherent energy $-199.95 < E_{\text{inh}} < -199.90$, and these states have a high mutual overlap as well as a significant overlap with the ground state ($Q = 0.7$). As a consequence, the average overlap grows $\langle Q \rangle \approx 0.5$. However, the final transition towards $\langle Q \rangle > 0.7$ when all configurations correspond to the ground state only happens at lower temperatures, $T < 0.02$, thus explaining the peculiar two-step increase in Fig. 3b for $N = 66$.

This behavior can also be rationalized by analyzing the

configurations themselves (see Supplemental Fig. 3a,c). The ground state features a pronounced symmetry along the $x = y$ and $x = -y$ axes. The configurations remain amorphous but they already remind of a quasicrystal. The excited state featured in Fig. 3c has the same symmetries. One can thus conclude that $N = 66$ features some kind of order or symmetry unlike the completely disordered structures seen for the other system sizes, and is therefore an outlier.

## ACKNOWLEDGMENTS

MO thanks the support by MIAI@Grenoble Alpes and the Agence Nationale de la Recherche under France 2030 with the reference ANR-23-IACL-0006). This work was also supported by a grant from the Simons Foundation (LB: #454933, GB: #454935). GB acknowledges funding from the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). LB acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grants ANR-20-CE30-0031 (project THEMA) and ANR-24-CE30-0442 (project GLASSGO).
<figure><img src="image_5.png" /><figcaption>FIG. S1. Dynamics of the KA2D model. a. Incoherent scattering function (ISF), $F_s(t) = \langle N^{-1} \sum_i e^{-iq \cdot (r_i(t) - r_i(0))} \rangle$, for $|q| = 2\pi\sigma_{11}^{-1}$ and $N = 704$. b. Temperature-dependence of the structural relaxation time $\tau_{\alpha}^{\text{ISF}}$ defined as $F_s(\tau_{\alpha}^{\text{ISF}}) = e^{-1}$. The purple line shows the extrapolation to $\tau_{\alpha}^{\text{ISF}} = 10^{12}$ which is a common definition of the glass transition temperature, $T_g(N = 704) = 0.147$. $\tau_{\alpha}^{\text{BB}}$ is the relaxation time of the bond-breaking correlation [28], showing very similar scaling behavior. c/d. Same as panels a/b for $N = 77$. The extrapolation in d yields, $T_g(N = 77) = 0.156$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. S2. Sampling configurations using parallel tempering. a. Time-dependent $c_V(t)$ for different temperatures $T$ for $N = 44$. All curves reach a clear plateau, thus indicating equilibration. b. Probability distribution of inherent energies for two independent sets for $N = 44$. Both sets yield basically the same probabilities. c. Same as b but zoomed to emphasize the lowest inherent energies. d/e/f. Same as a/b/c for $N = 77$.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

FIG. S3. Sampling configurations using population annealing (PA) and reweighting (RW) for $N = 77$. **a.** Time-dependent $\langle E_{\text{pot}} \rangle_{N_c}$ for different temperatures for $N = 77$. Here, 'time' $t = N_c\tau$ corresponds to the number $N_c$ of initial configurations used for PA and RW. All curves clearly reach a plateau, thus indicating equilibration. Reweighting shows a much slower decay for $T = 0.057$ (dotted line) thus indicating insufficient equilibration. **b.** Probability distribution of inherent energies for PA at different temperatures. Low energies are increasingly populated, until only the ground state is populated at $T = 0.006$. **c.** Same as **b** for RW. This procedure produces incorrect results at low temperatures $T < 0.07$, finally selecting a single state away from the ground state.

<figure><img src="image_8.png" /></figure>

FIG. S4. Different contributions to the configurational entropy for $N = 77$. **a.** Temperature-dependence of different contributions to the total entropy $s_{\text{tot}}$. The term $s_{\text{tot}}^{\text{old}} - s_{\text{tot}}$ describes the difference between the total entropy using finite size corrections, $s_{\text{tot}}$, and Eq. (10) in Ref. [36] without finite-size corrections, $s_{\text{tot}}^{\text{old}}$. All other terms are defined in Methods Section V. **b.** Temperature-dependence of the different contributions to the configurational entropy $s_c$.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. S5. Calculating the overlap $Q$. a. Illustration of how the maximal overlap $Q$ is calculated. Configurations 1 and 2 are independently sampled and appear distinct. However, mirroring the $x$-coordinate and subsequently apply the shift $x \to x-2.5$, $y \to y+1.5$ reveals that the two configurations are actually very similar with overlap $Q \approx 1$. b. Same as Fig. 3a for $N = 44$. c. Same as Fig. 3a for $N = 77$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. S6. The special case $N = 66$. a. Snapshot of the ground state for $N = 66$. Color code represents the particle type. b. Probability distribution of inherent energies for $N = 66$ at $T = 0.035$. c. Snapshot of the fourth excited state with energy $E_{inh} = -199.935$. Color code for the overlap with ground state is as in Fig. 4g.</figcaption></figure>