Ludovic Berthier

Patrick Charbonneau

Daniele Coslovich

Andrea Saverio Ninarello

Misaki Ozawa

Sho Yaida

Configurational entropy measurements in extremely supercooled liquids that break the glass ceiling

Introduction-In his landmark 1948 paper, Kauzmann gathered experimental data for several glass-forming liquids and found that they all showed a steep decrease of their equilibrium configurational entropy upon lowering temperature towards their glass transition [1]. Theoretically, the nature of a thermodynamic glass transition associated with a vanishing configurational entropy is well-understood at the mean-field level [2][3][4], suggesting that glass formation is accompanied by a rarefaction of available disordered states [5]. Its pertinence beyond the mean-field framework, however, remains controversial [5][6][7][8]. In particular, it is still not known whether such entropy reduction is the core explanation for glass formation. Experimental measurements are carried out over too limited a temperature range, within boundaries that have remained essentially unchanged since Kauzmann's work and thus form a solid glass ceiling. In addition, experimental determinations of the configurational entropy are marred by approximations that influence their physical interpretation [9][10][11]. Computer simulations can potentially provide more precise estimates [12,13], but have so far been restricted to a temperature range that is not experimentally relevant.

Can the debate over the role of configurational entropy ever be settled? At first sight, closure appears unlikely for two main reasons. (i) Measuring the configurational entropy below the experimental glass transition seems logically impossible, because experiments are constrained by their own duration, which fixes an upper limit to the accessible thermalization timescale, τ . Specifically, τ /τ 0 ∼ 10 13 for molecules [14] (where the relaxation time at the onset temperature is τ 0 ≈ 10 -10 s) and τ /τ 0 ∼ 10 5 for colloids [15] (where τ 0 ≈ 10 -1 s). The situation for computer simulations is even worse. Current approaches access at most τ /τ 0 ∼ 10 5 , which is eight orders of magnitude behind molecular liquid experiments, and numerical progress has been slow. The two to three decades gained over the past 35 years [15][16][17] are mostly thanks to hardware improvements. At this pace, another century would be needed before simulations attain experimentally-relevant conditions. The glass ceiling thus appears unbreakable. (ii) There is a fundamental methodological ambiguity as to which configurational entropy should be measured in order to match theoretical calculations. Qualitatively, the configurational entropy is defined by subtracting vibrational contributions from the total entropy [1,12,18]. What is specifically meant by "vibrations" in amorphous solids, however, is ill-defined in general [5] and difficult to measure in practice [1,13]. Hence consistently determining the configurational entropy is in itself a difficult challenge, that may be underestimated in the literature.

Here, we solve both of these major problems at once. First, we take advantage of the flexibility offered by computer simulations to dramatically accelerate the equilibrium sampling of configuration space [19][20][21]. Namely, we use a system optimized for the nonlocal swap Monte Carlo (MC) algorithm, which enables its extremely fast thermalization. We establish that this approach surpasses any current alternative, and even experimental protocols. Second, we measure four proxies for the configurational entropy by deploying state-of-the-art computational tools to characterize in-silico configurations that are more deeply equilibrated than their experimental analogs [13,20,22,23] and obtain consistent results that have a clear physical interpretation. By combining these developments for a realistic model glass-former, we shift computer simulations from lagging eight decades behind experiments to exploring novel territory in glass physics. In particular, our measurements validate Kauzmann's observations that the configurational entropy decreases steeply towards the glass temperature, and extend these observations to a regime previously inaccessible.

Results-We simulate a three-dimensional polydisperse mixture of hard spheres, as in [20], which is a good model for colloids used in experiments [15,24]. We show in Appendix I that our methods and conclusions also apply to particles with soft and more complex interactions. We control the volume fraction φ, and measure pressure P to report the (unitless) reduced pressure, Z = P/(ρk B T ), where ρ is the number density, and k B T the thermal energy. This natural control variable for hard spheres plays a role akin to the inverse temperature in thermal liquids [25]. Detailed information about the simulations is provided in Appendices A, B, and C. Swap MC complements standard translational MC moves with nonlocal moves that exchange randomly-chosen pairs of particles, ensuring equilibrium sampling. Detailed tests of thermalization of all glassy degrees of freedom are reported in Appendix D, see also Ref. [21]. We demonstrate the extreme speedup actually achieved by swap MC for this model in Fig. 1, in which the structural relaxation time τ for both MC sampling methods is reported as the system approaches its glass transition. Note that the rapid increase of τ in standard MC simulations resembles the fragile super-Arrhenius behavior of standard glass-formers [5]. We can only indirectly assess fragility beyond the reported numerical regime, which we estimated to be m ≈ 50. We have fitted several empirical forms to our measurements, which thermalize up to Z ≈ 27, to estimate the experimental glass transition at τ /τ 0 = 10 13 (see Appendix E). Use of various fits reflects the well-known uncertainties associated with the empirical description of data measured over a large dynamical range [26]. The fits give consistent locations for the glass ceiling, Z g ≈ 32-34, as highlighted in Fig. 1. Remarkably, this dramatic slowdown is completely bypassed by swap MC sampling, which thermalizes the system up to Z ≈ 38 > Z g . Even most conservative extrapolation indicates that we access a dynamical range that is broader than in experiments. Meanwhile, the two-point structure barely budges (see Figs. 1B-C), which is a telltale sign of glassiness [5] and a confirmation that both crystallization and more subtle fractionation effects are absent. Visual inspection of particle configurations further confirm these conclusions (see, e.g., Figs. 1D-E). We are therefore in the unique position of studying at equilibrium a homogeneous supercooled liquid beyond the experimental glass ceiling.

We then turn to measuring the configurational entropy, s conf , in these extremely supercooled configurations. The numerical procedures leading to the four estimates of s conf are shown in Fig. 2. Further details are provided in Appendices F, G, and H. In Method 1, we determine the configurational entropy from its most conventional definition, s conf = s tot -s vib , as used in many experimental and simulation studies [1,12,18,22]. The total entropy of the equilibrium fluid, s tot , is measured by thermodynamic integration from the dilute ideal gas limit to the target volume fraction, while the vibrational contribution, s vib , is measured by Frenkel-Ladd thermodynamic integration [22,27]. The latter integration is over the amplitude of the Hookean constant, α, of a spring that constrains each particle to reside close to the position of a quenched reference equilibrium configuration. This requires estimating the mean-squared distance δr 2 between the reference and constrained systems over a broad range of α values, as illustrated in Fig. 2(A). In continuously polydisperse systems, special The Frenkel-Ladd method to obtain the vibrational entropy s vib performs a thermodynamic integration of the mean-squared distance δr 2 between a reference equilibrium configuration and a copy of the system constrained by a harmonic potential of strength α. The integration is carried out from αmax, for which the system behaves as an Einstein solid (indicated by the dashed line δr 2 = 3/(2α)) to αmin, for which particles are trapped by their own cages on the vibrational time scale. (B) Method 2: The numerically-determined Franz-Parisi potential V (Q) is used to measure the configurational entropy as

s conf = V (Q high = 0.8) -V (Q low ≈ 0.05). (C) Method 3:
The evolution of the overlap Q with the biasing field ε reveals a first-order jump at a value ε for which Q = 1/2 (dashed line). Then,

s conf = ε (Q high -Q low ). (D) Method 4:
The decay of the cavity overlap correlation function QPTS(R) with cavity radius, R, defines the point-to-set correlation length ξPTS ∝ s

-1/(d-θ) conf .
care is also needed to account for the mixing contribution to the total entropy, because this contribution formally diverges [28,29]. The mixing entropy is thus determined from an independent, additional set of simulations [29] (see Appendix F). Method 1 is equivalent to partitioning configuration space into basins of attraction of inherent structures [7]. The resulting estimate of the configurational entropy thus counts the number of energy minima [12,30], which presumably overestimates the number of relevant basins in the free energy landscape [31].

Methods 2 and 3 are both based on the Franz-Parisi theoretical construction [32], which expresses the equilibrium free energy of the liquid, V (Q), in terms of a global order parameter, the overlap Q. The overlap between two configurations is defined as Q = N -1 i,j θ(a -|r 1,i -r 2,j |), where θ(x) is the Heaviside function, r 1,i and r 2,j are the positions of particle i and j within configuration 1 and 2, and a is a fraction of the average particle diameter. By definition, Q quantifies the similarity between the coarse-grained density profiles of two configurations. To compute V (Q), we introduce a coupling between a quenched reference equilibrium configuration and a copy of the system through a field ε conjugate to Q [13,32]; ε constrains the collective density profile, whereas α in Method 1 constrains single-particle displacements. We define V (Q) = -lim ε→0 T N ln P (Q) , where P (Q) is the equilibrium probability distribution of the overlap for a given reference configuration, and brackets denote averaging over these configurations. In Method 2, we follow [13] and use the free-energy difference s conf = V (Q high ) -V (Q low ) between the global minimum at Q low ≈ 0.05 and its value at Q high = 0.8 to obtain an estimate of s conf that is closest to its theoretical definition, see Fig. 2(B). Importantly, this estimate only exists for sufficiently supercooled states, for which Q high can be defined [13]. For the present system, this happens close to the mode-coupling crossover, Z c . In Method 3, we determine the value of the biasing ε needed to 'tilt' the potential V (Q), so that a first-order phase transition, at which Q jumps from Q low to Q high , takes place as illustrated in Fig. 2(C). We use the maximum variance of the overlap fluctuations to measure ε for each volume fraction studied. In practice, this is equivalent to determining the biasing field at which the overlap reaches Q = 1/2, see Fig. 2(C).

Method 4 builds on the physical idea that the decrease of the configurational entropy is directly responsible for the growth of spatial correlations quantified by the pointto-set correlation length, ξ PTS [33][34][35]. Following what is becoming common practice [23,35], we measure ξ PTS by pinning the position of particles outside a spherical cavity of radius R, equilibrating the liquid within it, and measuring the evolution of the overlap between interior configurations, Q PTS (R), with the cavity radius R, as shown in Fig. 2(D). The decay of Q PTS (R) is controlled by ξ PTS , and the variance of the overlap fluctuations also presents a maximum [23] very close to ξ PTS . Physically, ξ PTS thus represents the cavity size above which the system starts to explore a significant number of distinct states. With minimal hypothesis [34], it can be connected to the configurational entropy through

s conf ∝ ξ -(d-θ) PTS
, with an unknown exponent θ ≤ (d -1). Various values of θ have been proposed, including θ = 2 from saturating the inequality and θ = 3/2 from a wetting argument [6,36]. However, our measurements are consistent with both of these values and thus cannot unambiguously distinguish one proposal from the other.

We gather the four estimates of the configurational entropy in Fig. 3 to produce a plot akin to the original 1948 Kauzmann representation of s conf (T ) [1]. Although in the high-temperature liquid the configurational entropy is not sensibly defined [13], three of the four measures can still be estimated. Note that only this regime was accessible in earlier simulations [12,13,22]. In the more relevant low-temperature regime, our main finding is that the important conceptual and technical differences between the four methods nevertheless result in qualitatively consistent results. In particular, the three estimates (Methods 2-4) that closely follow the theoretical definition of the configurational entropy provide numerically indistinguishable results at low temperatures. The conventional estimate of the entropy (Method 1) is larger, as expected [31], but its temperature evolution remains qualitatively consistent with the other methods. All our estimates of s conf thus exhibit a steep decrease as Z increases towards the glass phase, which is consistent with the seemingly fragile behavior of the model in Fig. 1. Although a quantitative extrapolation is hard to control, our measurements robustly suggest that s conf may vanish near Z ≈ 1/0.022 ≈ 45. We thus conclude that even for a simple glass-forming system equilibrated deeper in the landscape than any previously studied material, the trend discovered 70 years ago by Kauzmann is confirmed when more precise estimates of s conf are adopted, and persists even below the experimental glass temperature. We further show in Appendix I that similar observations can be performed for a model with a continuous pair potential, suggesting our methodological progress and physical conclusions are not restricted to hard spheres, and likely apply more generally. Note that while continuous polydisperse distributions are commonplace in colloidal suspensions, a molecular liquid with a sufficiently large number of components to approximate a continuous size distribution has yet to be considered.

Discussion-Our point-to-set measurements go beyond Kauzmann's observation by establishing that the decrease in s conf is accompanied by an increase of static spatial correlations as the glass ceiling is crossed. This result reinforces a recent experimental report based on nonlinear dielectric measurements [37]. In absolute value, the measured static length scale at the experimental glass transition appears somewhat smaller than previous estimates based on dynamical correlations [38,39], but remains compatible with the modest growth expected from general arguments based on thermally activated scaling [6,34,36] and decorrelation between static and dynamical length scales [40]. Our particle-based resolution of such correlations further provides a direct visualization of the spatial profile of the overlap within a spherical cavity (see insets in Fig. 3). In particular, within a cavity comprising about 200 particles, the positions of particles freely fluctuate near the onset pressure, but become strongly correlated over the entire cavity for the largest pressure shown. The spatial extent of static correlations is thus directly revealed.

The important methodological advances achieved here regarding the thermalization of supercooled liquids and the measurement of configurational entropy therefore support a thermodynamic view of the glass formation based on the rarefaction of metastable state accompanied by growing static correlations that is devoid of the experimental ambiguities and that extends to a temperature regime that has never been explored before.

Appendix A: Model

We study a three-dimensional hard-sphere model, for which the pair interaction is zero for non-overlapping particles and infinite otherwise. Systems have a continuous size polydispersity, with particle diameters σ randomly drawn from the distribution f (σ) = Aσ -3 , with σ ∈ [σ min , σ max ] with normalization constant A. Our model is the same as that studied in Ref. [20], with a measure of size polydispersity ∆ = σ 2 -σ 2 /σ, where

• • • = dσf (σ)(• • • ), of ∆ = 23%
, and σ min /σ max = 0.4492. The average diameter, σ, defines the unit of length. We simulate systems composed of N particles in a cubic cell of volume V under periodic boundary conditions [START_REF] Allen | Computer simulation of liquids[END_REF]. Depending on the chosen method to estimate the configurational entropy (see Appendix F), we simulate systems with either N = 1000, 8000 (Method 1) or N = 300 (Method 2 and 3). Cavities for Method 4 are carved from bulk configurations with N = 8000. The relaxation times shown in Figure 1A are obtained from samples with N = 1000. Given these parameters, the system is then uniquely characterized by its volume fraction φ = πN σ 3 /(6V ), and we frequently report the data using the reduced pressure Z = P/(ρk B T ), where ρ, k B , and T are the number density, Boltzmann constant and temperature, respectively. Without loss of generality, we set k B and T = 1/β to unity. The pressure P is calculated from the contact value of the pair correlation function properly scaled for a polydisperse system [START_REF] Santos | [END_REF].

Appendix B: Methods

To obtain equilibrium fluid configurations deep in the glassy regime, we perform Monte-Carlo (MC) simulations with both translational displacements and non-local particle swaps [19,20,[43][44][45][46][START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Brumer | [END_REF][START_REF] Gutiérrez | Procaccia[END_REF]. The two types of moves are selected randomly: with probability 0.8 we attempt a translational displacement, and with probability 0.2 we attempt a swap. Translational displacements are uniformly drawn over a cube of side 0.115. For swaps, two randomly selected particles exchange diameter. In both cases, proposed moves are accepted if no overlap is created. Following Ref. [21], we also immediately reject swaps between particles whose diameters differ by more than 0.2.

We measure the equilibrium relaxation time τ both with and without the swap moves from the time-decay of the self-intermediate scattering function, F s (k, τ ) = 1/e, where k = 5.25 is the wavenumber chosen slightly below the first maximum of the static structure factor. Note that the particles' diameters can change during a swap MC simulation, but their trajectories are continuous. Relaxation times are measured in units of MC sweeps, comprising N MC moves, irrespective of their type.

Thermalized systems at each state point are obtained in the same way for both standard and swap MC dynamics. We measure the relaxation time τ and ensure that for each state point simulations of a total duration of at least 100τ can be performed. We also check for the presence of aging effects in time correlation functions, and we measure the static structure factor, the pair correlation function and the equation of state over long simulations, paying attention to any temporal drift that could signal either improper thermalization, incipient crystallization or demixing of particles with distinct sizes. Selected results for the evolution of the structure factor with volume fraction are presented in Figs. 1B-C. Over the extreme range of densities shown here, the static structure evolves very little. Similarly, a very modest evolution is seen when the partial structure factor of the smallest particles is measured. A large increase of the low-k value of these quantities, or the emergence of discrete peaks would signal that demixing or crystallization is taking place. In fact, we have found that measuring the relaxation time for the swap simulation is the most sensitive test of thermalization, because purely static observables may appear thermalized over long simulation times, whereas the system is in fact nearly arrested in a glass state within which sampling is inefficient.

We introduce two dynamical reference states: (i) the onset of slow dynamics at φ 0 ≈ 0.56 (Z 0 ≈ 18), above which the time decay of correlation functions is non-exponential [50], and (ii) the mode-coupling crossover φ c ≈ 0.598 (Z c ≈ 23.5), at which a power-law fit extrapolates a divergence of the relaxation times [20]. Note that these particular definitions are not unique [51], but are sufficiently accurate for our purposes, where these values are simply used for qualitative reference. for the hard-sphere system studied in this work. The solid line is the empirical expression for the equation of state from Refs. [52,53], BMCSL, which describes our data well, except at very large volume fractions, where it slightly overestimates Z. A modified version of BMCSL [54] is given by the dashed line. In this section, we present evidence that strongly supports the equilibration of our samples. In addition to the tests mentioned in the Methods section of the main text, we also verified the equivalence of derivative and fluctuation expressions for the compressibility, which is a well-known equilibration check [55]. The validity of this sum rule is, however, not necessarily a strong test of equilibration, because it may appear to hold within the noise of the data, even in systems that are not fully equilibrated. A more stringent test of equilibration is to show that all the slow degrees of freedom relax within the observation time. In a structural glass former, the relevant slow degrees of freedom are those associated with collective density fluctuations and with cage breaking. To measure the decorrelation of the former, we introduce the time-dependent collective overlap function

Z
Q c (t) = N -1 i j Θ(a -| r i (t) -r j (0)|), (D1)
where Θ(x) is the Heaviside function and a = 0.346. The decay of Q c (t) quantifies the relaxation of density fluctuations up to a coarse-graining length a. In order to analyze cage breaking, here defined as the process of a particle losing some of its neighbors within a given time span [56], the following procedure is employed. First, we identify the neighbors of a particle by a radical Voronoi tessellation obtained using the Voro++ package [57]. Note that the average number of neighbors measured for our densest states is about 14. Second, we measure the probability, P

b (t), that a particle has changed all but k of its initial neighbors after a time t. For small k, the decay of P (k) b (t) provides evidence for the complete restructuring of the cages and ensures that a pseudo-molecular structure is not found in the fluid.

Representative results of Q c (t) and P

b (t), measured in swap MC simulations of the N = 1000 system, are shown in Fig. S2. All these functions decay to trivial plateaus within our simulations, with the possible, albeit mild, exception of the data point at the highest packing fraction. We also note that although the cage breaking probability decays more slowly than Q c (t), the associated relaxation times grow in a similar way. Our simulations thus pass this very stringent equilibration test, in addition to those already mentioned in the main text and in Ref. [20]. For the soft polydisperse model studied in Sec. I of the SM, a similar analysis was performed in Ref. [58]. relaxation times at most of order τ /τ 0 = 10 5 , where τ 0 = 10 4 MC sweeps is the value of τ at the onset of slow dynamics, our dynamical data have to be extrapolated to locate this ceiling. We fit measured relaxation times to various functional forms up to τ /τ 0 ≤ 10 5 , and then extrapolate up to the vicinity of the glass ceiling. In an effort to obtain an estimate as unbiased as possible, we consider a range of possible functional forms [5,25], as detailed below. We present in Fig. S3 the result of this exercise. In particular, note that whereas extrapolating the location of a putative divergence of τ is extremely delicate, extrapolating the location of the glass ceiling is in fact rather well-controlled.

The first functional form is the Vogel-Fulcher-Tammann (VFT) expression

τ = τ ∞ exp A (Z vft -Z) δ , (E1)
where τ ∞ , A, the exponent δ and the critical pressure Z vft are free parameters. Whereas δ = 1 is traditionally used, more recent experimental and numerical studies [15,25] favor δ = 2. In Fig. S3 we consider both. As noted before [20], the fit with δ = 1 yields a critical pressure Z vft ≈ 38 that falls within the range in which swap MC sampling equilibrates. Although the resulting fit is good over the range covered by standard MC dynamics, it overestimates the growth of the relaxation time beyond that range, which undermines the very logic behind the proposed scaling form. When δ = 2 is imposed the fit is still very good, but we now obtain Z vft ≈ 45, which is at least beyond the equilibrium range accessible with the swap method. Although this fit provides a more consistent description of the dynamical data, resolving one form from the other is beyond the scope of the current study.

The second functional form is the parabolic law proposed by Elmatad et al. [59,60] in the context of facilitated models,

τ = τ ∞ exp A(Z -Z 0 ) 2 , ( E2
)
where Z 0 = 17 is around the onset of slow dynamics and the fit is made over the range Z > Z 0 . In contrast to the VFT law, this parabolic expression does not invoke a divergence of the relaxation time at any finite pressure; τ only diverges when Z also diverges. Note that because Z ∝ 1/T , this form is equivalent to fitting the relaxation time of a supercooled liquid without any finite-temperature divergence [25]. As a result this expression necessarily provides a less divergent extrapolation of the relaxation time at large pressures, but still accounts well for the curvature, i.e., the fragility, of the measured dynamical data. The leading order at large Z, τ ∝ exp[AZ 2 ], indeed grows faster than the Arrhenius law. Like the VFT form, this fit is quite good over the measured range of relaxation times, see Fig. S3.

For completeness, we also include a simple Arrhenius fit to the high-Z portion of the equilibrium data

τ = τ ∞ exp[AZ]. ( E3
)
We find this fit not to be very good as it does not capture the fragility of the system. It is also inconsistent with the thermodynamic behavior of the system reported in the main text, because the configurational entropy typically does not vary much in glass-formers with an Arrhenius-like behavior. This fit is thus most probably incorrect in the sense that it underestimates considerably the evolution of the relaxation time of the system. We conclude that the first two families of expressions account well for the non-Arrhenius dependence of the relaxation time data observed in Fig. S3 and yield reasonable descriptions of the available equilibrium data over several orders of magnitude. Whereas the VFT law with δ = 1 can be logically ruled out by the swap MC measurements, the other two fits cannot be excluded on the basis of any further measurement we could perform. As is common in glass simulations, it is thus difficult to discriminate between fitting forms with or without a finite temperature singularity [5]. These forms nonetheless allow us to locate the glass ceiling Z g , for which the relaxation time is 10 13 larger than its value at φ 0 . Fits to the relevant expressions yield estimates ranging from Z g ≈ 32 (VFT law with δ = 1) to Z g = 34 (parabolic law). These values are used as boundaries of the glass ceiling box in Fig. 1 of the main text.

As discussed above, there are several reasons for which we do not expect the Arrhenius expression to describe accurately the relaxation data at high Z. If we nonetheless considered it, the equilibrium swap MC data would still be able to equilibrate the liquid at pressures higher than the glass transition it predicts. Hence, this very conservative extrapolation allows us to confidently state that we have successfully broken the glass ceiling.

Appendix F: Method 1: Conventional definition

In order to compute the configurational entropy, s conf , we define

s conf = s tot -s vib (Method 1), ( F1
)
where s tot and s vib are the total and vibrational entropy, respectively. Note that this definition is common in experimental and computational studies [1,12,18,22,61,62].

Total entropy

The total entropy is obtained by thermodynamic integration from the ideal gas limit (φ → 0) up to the target volume fraction φ,

s tot (φ) = 5 2 -ln 6φ πM 3 -ln Λ 3 - φ 0 dφ (Z(φ ) -1) φ + s mix , (F2)
where

M k = σ k (k = 1, 2, 3, • • •
) are the k-th moments and Λ = 2πβh 2 /m is the thermal de Broglie wavelength.

Without loss of generality, we here set Λ = 1. For continuous polydisperse systems one needs to pay special attention to the mixing entropy, s mix [28,63]. We get back to this point in Sec. F 3. For now, we report s tot -s mix for N = 1000 and 8000 in Fig. S4(A).

In order to validate the numerical thermodynamic integration, we also consider an analytical approximation of the equation of state (EOS). The polydisperse version of the Carnahan-Starling EOS, i.e., the so-called BMCSL EOS [52,53],

Z BMCSL (φ) = 1 1 -φ + 3M 1 M 2 M 3 φ (1 -φ) 2 + M 3 2 M 2 3 (3 -φ)φ 2 (1 -φ) 3 , ( F3
)
is known to describe experiments and simulations of polydisperse hard-sphere systems rather well in the liquid regime.

We also consider a modified version of Z BMCSL ,

Z modBMCSL (φ) = 1 + h(φ)(Z BMCSL (φ) -1), ( F4
)
where h(φ) = 0.005 -tanh(14(φ -0.79)) [54]. Both EOSs trace our simulation data very well, as can be seen in Figs. S1 and S4(A). These EOSs can thus also be used to extrapolate the configurational entropy toward very high volume fraction (Sec. F 4).

Vibrational entropy

The vibrational entropy, s vib , is obtained by Frenkel-Ladd (FL) thermodynamic integration [22,[START_REF] Frenkel | Understanding Molecular Simulation[END_REF] by performing MC simulations of a constrained system with Hamiltonian

βH(α) = βH(0) + α N i=1 (r i -r 0i ) 2 , ( F5
)
for a template configuration, {r 0i }, obtained from an equilibrium target system under H(0). In short, the FL method integrates from a large α max 1, at which particles experience a nearly pure harmonic oscillator, down to a very weak α min 1, at which particles vibrate within cages. The vibrational entropy is then obtained by

s vib = 3 2 -ln Λ 3 - 3 2 ln α max π + αmax αmin dαδr 2 (α) + α min δr 2 (α min ), ( F6
)
δr 2 (α) = 1 N N i=1 (r i -r 0i ) 2 α , (F7)
where (• • • ) α and [(• • • )] denote the thermal average with Hamiltonian H(α) and averaging over template configurations, {r 0i }, respectively, and δr 2 in Eq. (F7) is the mean-squared displacement shown in Fig. 2(A) of the main text. Numerical integration of δr 2 is performed from α max to α min . We set α max = 7.5 × 10 6 , which is well into the harmonic-oscillator scaling regime [see Fig. 2(A) of the main text],

δr 2 (α max ) = 3 2α max . ( F8
)
The choice of α min is such that δr 2 of the constrained system is comparable to the mean-squared displacement of the target system without constraint. The last term in Eq. (F6) corresponds to assuming a constant δr 2 for α ∈ [0, α min ] in the thermodynamic integration. In Fig. S4(B), we show s vib for N = 1000 and 8000, with α min = 1.5 × 10 0 , 7.5 × 10 0 , and 3.0 × 10 1 . At high volume fraction φ, the results are insensitive to the choice of α min , as expected. We can thus confidently set α min = 7.5 × 10 0 . Also, we empirically observe a linear relation between s vib and φ, which allows us to linearly extrapolate the fit and thus the configurational entropy (see Sec. F 4).

Mixing entropy

The mixing entropy of the continuous polydisperse system formally diverges in the thermodynamic limit because each particle in the system then belongs to one of an infinite number of different components [28,63]. One, however, can subtract from this quantity a physically relevant contribution, which we call the effective mixing entropy, s * mix . The main idea is that a continuous polydisperse system can be regarded as an effective M * -component system (see Ref. [64] and references therein). We then assume that the effective M * -component system shares physical properties, in particular a same (free-)energy basin, with the original continuous polydisperse system. Here, a practical description of how to estimate s * mix from this scheme in our system is provided; a full explanation is provided in Ref. [29]. To estimate M * , we first decompose the distribution f (σ) into M species, as shown in the inset of Fig. S5(B). We define M species by dividing f (σ) into equal intervals ∆σ = (σ max -σ min )/M , such that each species occupies the same fraction of the total volume, Aρπ∆σ/6 = C, where C is a constant. Note that M is an integer and that M → ∞ corresponds to the original continuous polydisperse system. We then perform a quench of the discretized system from the original configuration, and determine whether or not it remains in the same glassy basin as the original system by measuring the mean-squared displacement before and after the quench,

∆ M = 1 N N i=1 |r IS M i -r 0i | 2 , ( F9
)
where {r IS M i } is the inherent structure configuration of a M -discretized system. If M is large, the discretized system is almost identical with the original continuous polydisperse system and stays within the same basin, hence ∆ M ∆ M →∞ . If M is small, however, discretization destroys the glassy basin of the original system. The system thus structurally rearranges into a different glassy basin, and ∆ M ∆ M →∞ . This determination is done here by considering inherent structures, which for hard-sphere systems correspond to an out-of-equilibrium compression up to jamming [30]. We determine M * as the crossover between these two limit cases.

More precisely, we follow the following algorithm.

1) Obtain an equilibrium configuration of the original continuous polydisperse system, {r 0i }, for the initial configuration.

2) Discretize the diameters σ of the original system into M species, keeping φ fixed.

3) Quench the system to its inherent structure, {r IS M i }, using the algorithm described in Refs. [65,66]. 4) Repeat 1) -3) for a range of M and over different initial configurations. 5) Determine the crossover M * from ∆ M as a function of M .

In Fig. S5(A), we show how ∆ M evolves with M . At large M , ∆ M is flat and near ∆ M →∞ . The discretized system stays in its original basin. Upon decreasing M , however, ∆ M deviates from ∆ M →∞ , indicating that the discretized system escapes its basin. We estimate M * as an onset of this deviation by fitting the two linear lines from large and small M regions. We find a weak φ dependence with M * ≈ 9 -10 for all φ considered. We use linear fits to more precisely locate the crossover at M * = 9.13, see Fig. S5(A).

We map M * onto the effective mixing entropy s * mix by computing s mix for different M ,

s mix (M ) = - M m=1 X m ln X m , (F10)
X m = σmin+m∆σ σmin+(m-1)∆σ dσf (σ), (F11)
where X m is the concentration of species m. Fig. S5(B) shows the resulting s mix (M ). From the mapping indicated by the arrows, we obtain s * mix = 1.98. Note that s * mix in Fig. S5(B) varies relatively weakly with M * and so a more precise estimate of the value of M * is not needed for our purposes. The inset shows the continuous polydispersity distribution, f (σ), and the weights of delta functions for the discretized distribution with M = 5, {Xm} m=1,...,5 , as an example.

Configurational entropy

Using the above results, we obtain the configurational entropy from s tot -s vib as a function of φ. As can be seen in Fig. S6, no significant difference is observed between N = 1000 and 8000. We also find that the combination of s tot from the (modified) BMCSL EOS in Fig. S4(A) and s vib from the linear fit in Fig. S4(B) gives a reasonable extrapolation toward s tot -s vib = 0, and thus an estimate for the Kauzmann transition, φ K ≈ 0.68 (Z K ≈ 45). In Fig. 3 from the main text we use the curve extrapolated by the modified BMCSL EOS. considerably. It then becomes difficult to perform an accurate sampling of the overlap fluctuations imposed by the Hamiltonian (G2).

To accurately sample the density regime explored in the present work, we have implemented parallel tempering MC moves [69]. We conduct the n simulations needed for the umbrella sampling at volume fraction φ in parallel, and propose MC exchange moves between neighboring simulations, as characterized by nearby sets of parameters, say (k , Q ) and (k +1 , Q +1 ). An exchange between simulations and + 1 is proposed with a low frequency (typically every 10 3 MC sweeps) and they are accepted with a Metropolis acceptance rate given by the Hamiltonians H and H +1 , ensuring detailed balance.

Because each simulation now performs a random walk in the parameter space defined by {(k , Q ), = 1, • • • , n}, the sampling of overlap fluctuations is greatly enhanced. For the method to be efficient, however, we need to adjust the biasing potentials W (Q) such that the distributions P (Q) of neighboring simulations overlap significantly. We have used up to n = 20 to gather data. Thermalization was carefully checked by running long simulations (up to 10 10 MC sweeps per simulation), in order to make sure that each state point was visited several times by all simulations via the replica exchange. This represents a significant numerical effort.

Histogram reweighting

Having obtained thermalized results from n biased simulations run in parallel, we process the simulation outcome using multi-histogram reweighting methods to reconstruct the unbiased probability P (Q) from the independently measured P (Q),

P (Q, ε, φ) = n =1 P (Q, ε, φ) n =1 e -βW /Z , ( G5
)
where Z are defined self-consistently as

Z = 1 0 dQ n m=1 P j (Q , ε, φ) n m=1 e β(W -Wm) /Z m . ( G6
)
Note that the actual value of ε used in simulations plays no conceptual role because reweighting directly provides P (Q, ε , φ) from P (Q, ε, φ) for distinct field values ε and ε :

P (Q, ε , φ) = P (Q, ε, φ)e -βQ(ε -ε) 1 0 dQ P (Q , ε, φ)e -βQ (ε -ε) . ( G7
)
Two values of the field ε are nonetheless particularly relevant to determining the configurational entropy, as described in the following two subsections.

Method 2

First, the Franz-Parisi free energy V (Q) is obtained as

V (Q) = -lim ε→0 T N ln P (Q, ε, φ) , ( G8
)
where [• • • ] denotes averaging over the quenched reference configuration 2. We used 60 independent reference configurations for each value of the volume fraction. Because V (Q) is only defined up to an additive constant, we arbitrarily adjust it, such that V (Q low ) = 0, where Q low ≈ 0.05 is defined as the location of the global minimum in V (Q). This additive constant is irrelevant because only the free energy difference,

s conf = V (Q high ) -V (Q low ) (Method 2) (G9)
is needed to determine the configurational entropy. In this work, we set Q high = 0.8 as motivated in the following subsection.

Method 3

Second, we obtain the critical field value ε = ε needed to induce a phase coexistence between low and high overlap states [32,68]. In practice, we use the strength of the reweighting method to finely explore a range of ε values, and define ε as the field strength for which the overlap shows a value intermediate between Q low and Q high , where the distribution P (Q, ε , φ) shows two peaks of equal amplitude, and where the variance of the overlap fluctuations (the susceptibility) is maximal. The position of the second peak of P (Q, ε , φ) at coexistence sets Q high . Because its volume fraction dependence is very weak, we fix it close to its average value and use Q high = 0.8 for all φ. Note that a secondary minimum cannot exist in the large system size limit of finite-dimensional simulations [13], which is why we resort to the above definition of Q high .

The field ε has a simple graphical interpretation. It represents the amplitude of the field needed to 'tilt' the potential

V (Q) towards coexistence [see Fig. 2(B) of the main text]. Because the relation V (Q high ) ≈ ε (Q high -Q low )
, holds to a good approximation, ε provides an estimate of the free energy difference V (Q high ) -V (Q low ). We can thus define the configurational entropy as

s conf = ε (Q high -Q low ) (Method 3). ( G10
)
Appendix H: Method 4: The point-to-set correlation

Definition of cavity core overlap and point-to-set correlations

The similarity between two configurations, X = {x i } and Y = {y i }, is characterized by using an overlap field, q X,Y (r), computed as in Ref. [23]. For each particle x i , we find the nearest particle y inn and assign an overlap value q X;Y (x i ) ≡ w x i -y inn , where

w(z) ≡ exp - z b 2 , (H1)
with b = 0.23. This function defines overlap values q X;Y (x i ) at scattered points {x i }. We then define a continuous function passing through these points, q X;Y (r). Specifically, we first perform a Delaunay tessellation of space and, to a point r within a simplex spanned by four points {x i } i=i1,i2,i3,i4 , associate a linearly interpolated value

q X;Y (r) = i=i1,i2,i3,i4 c i q X (x i ) , (H2)
where {c i } i=i1,i2,i3,i4 satisfies r = i=i1,i2,i3,i4 c i x i with the constraint i=i1,i2,i3,i4 c i = 1. We can similarly obtain q Y;X (r), allowing us to define the overlap field

q X,Y (r) ≡ 1 2 {q X;Y (r) + q Y;X (r)} . ( H3
)
In order to capture the similarity between two configurations near the center of the cavity, we also define the cavity core overlap

q c ≡ 3 4πr 3 c |r |<rc dr q X,Y (r ) , ( H4
)
where r c = 0.576 and r = 0 is the cavity center. This integral is numerically evaluated by Monte Carlo integration using 10 4 points.

For each volume fraction φ and cavity radius R, the point-to-set correlation function, Q PTS (R; φ), is evaluated by disorder-averaging over 40 cavity centers and, within each cavity, thermal-averaging over s prod pairs of equilibrated configurations (see Tables IV-X) [23]. We extract the point-to-set correlation length through the compressed exponential fit,

Q PTS (R; φ) = A exp[-{R/ξ PTS (φ)} γ ] + Q bulk PTS (φ), (H5)
with the compressed exponent γ = 4. Here the bulk value, Q bulk PTS , corresponds to the value at R = ∞ and is evaluated by taking 10 5 pairs of independent configurations in bulk samples.

At the point-to-set length scale, we expect the probability distribution function of core overlaps to display broad fluctuations. Figure S7 bears out this expectation. In particular, as anticipated in Ref. [23], we observe that the full disorder-averaged distribution becomes bimodal at high volume fractions, showing that we have indeed entered a deeply glassy regime that remains currently inaccessible for a standard binary Lennard-Jones liquid (see also [70]).

Parallel-tempering sampling with fuzzy-ensemble replicas

Denote positions of mobile particles inside the cavity x M i with i = 1, . . . , N cav and their associated diameters σ M i . Similarly, denote positions and diameters of pinned particles outside the cavity x P j and σ P j , respectively. For cavity point-to-set measurements, one samples cavity configurations x M i with probability

P x M i σ M i , x P j , σ P j (H6) exp -βU x M i , σ M i , x P j , σ P j ,
where denotes equality up to a normalization constant and a hard wall constraint x M i < R at the edge of the cavity of size R is imposed. Note that for the standard hard-sphere potential, U , this expression does not depend on the inverse temperature β ∈ (0, ∞). Properly sampling from this weight can be computationally demanding, because the equilibration time grows rapidly both as the cavity size decreases and as the volume fraction increases (see Fig. S8 and Sec. H 4). In these systems, phase space can actually become fully disconnected, which makes proper sampling impossible through (semi-)local Monte Carlo moves.

We here overcome this difficulty through a parallel-tempering algorithm [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]71]. The simplest scheme, analogous to that employed in Ref. [23] for a thermal binary Lennard-Jones liquid, would be to couple the above ensemble with replicas governed by a probability

P λ x M i σ M i , x P j , σ P j (H7) exp -βU x M i , λσ M i , x P j , σ P j ,
where the diameter of the particles inside the cavity is uniformly shrunk by λ. This approach, however, is highly inefficient for hard interactions. Replica-swap attempts are indeed very likely to be rejected, because even small λ increments can result in hard overlaps. More precisely, in order to have appreciable replica-swapping acceptance rates, the spacing between consecutive λ's must scale as O(1/N cav ), and thus O(N cav ) replicas are needed, as opposed to O( √ N cav ) for thermal systems [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. We therefore introduce a fuzzy ensemble governed by probability weights controlled by two parameters, α and λ,

P α, λ x M i , λ σ M i , x P j , σ P j (H8) exp - α 2 N cav 2 λ -λ 2 ×exp -βU x M i , λσ M i , x P j , σ P j ,
hence the shrinking factor λ is allowed to fluctuate. With this modification, a O( √ N cav ) scaling for the number of replicas needed is recovered. One caveat is that one of the replicas must obey the original sharp distribution defined in Eq. (H6)-formally corresponding to the limit α → ∞ with λ = 1-and the aforementioned challenge imposed by hard interactions still persists for that case. This problem is here surmounted by using a large number of replica-swap attempts near the original replica-O(1000N cav ) more frequent than others-as detailed in Sec. H 3 b.

Details of Monte Carlo moves

This section details the implementation of the Monte Carlo scheme for the cavity simulations.

a. Basic moves within a replica

Within each replica, we perform Monte Carlo moves consisting of local displacements, particle identity swaps, and, for fuzzy ensembles, diameter fluctuations.

• Local displacements consist of: (i) randomly choosing a particle i from the N cav mobile particles within the cavity; (ii) displacing particle i by ∆x = ln with uniform l ∈ [0, 0.1] and uniform n on the unit sphere S 2 ; and (iii) accepting displacement only if no hard overlaps are created. Note that an additional hard spherical wall at the edge of the cavity guarantees that mobile particles cannot leave the cavity-a rare instance, even at the lowest volume fraction studied.

• Particle identity swaps consist of: (i) choosing two distinct particles i and j within a cavity; and (ii) swapping their diameters if it results in no hard overlaps. (For the highest volume fraction φ = 0.640 with the cavity size R ≥ 3.11, in order to accelerate runs, we attempt these moves only for pairs with diameter difference λ|σ i -σ j | < 0.086.)

• Diameter fluctuations consist of: (i) uniformly drawing a shift ∆λ ∈ -0.2

α √ Ncav , + 0.2 α √
Ncav , and (ii) accepting the shift with probability p = min{p Λ , 1}, where

p Λ = exp - α 2 N cav 2 λ + ∆λ -λ 2 -λ -λ 2 , (H9)
if no hard overlaps are created.

One MC sweep consists of N cav MC trial moves of the above kinds for each replica. For large cavities with R > ξ PTS +1, each MC sweep consists of N cav MC trial moves with 70% (80%) local displacements and 30% (20%) particles identity swaps (parenthesis are for φ = 0.640) alone. For small cavities, however, proper sampling requires replica exchange and diameter fluctuations (see Fig. S8).

b. Parallel tempering

For parallel tempering, we label replicas a = 1, . . . , n, where a = 1 corresponds to the original ensemble with λ1 = 1 and α 1 = ∞. For the other replicas, we choose α a≥2 = 20, and { λa } are tuned to enable appreciable replica-swap rate (see below). For a ≥ 3, each MC sweep consists of N cav MC trial moves with 60% (70%) local displacements, 30% (20%) particles identity swaps, and 10% diameter fluctuations (parenthesis are for φ = 0.640 with R ≥ 3.11). For a = 1 and a = 2, a different scheme is employed. The two cavities are run as a pair with a large number of replica-swap attempts. Each MC sweep consists of 4N cav attempts, 50% replica-swap between x M i 1 and x M i 2 -accepted only when swapping λ 1 = 1 and λ 2 does not result in hard overlaps-, 15% (17.5%) local displacements respectively for a = 1 and a = 2, 7.5% (5%) particles identity swaps respectively for a = 1 and a = 2, and 5% diameter fluctuations for a = 2.

For each a ≥ 2, a replica-identity swap between x M i a , λ a and x M i a+1 , λ a+1 is attempted every 1000 MC sweeps on average, with acceptance probability p = min{p RS , 1}, where

p RS = exp -α 2 N cav λa+1 -λa (λ a+1 -λ a ) . ( H10
)
The replica parameters, { λa } a≥2 , are tuned to ensure sufficient replica-swap rates. In order to achieve this sampling, we first define the average of the fluctuating λ

λ a = dλdx M i λ P αa, λa x M i , λ σ M i , x P j , σ P j , (H11)
For the replica a = 2, we ensure that

λ 2 = 1 + O λ 2 2 -λ 2 2
, which typically requires λ2 > 1, because of the relatively high system pressure. Replicas are then added one by one, with λ2 > λ3 > . . . > λn , each time targeting a replica-swap acceptance rate of ∼ 20%. This process is stopped upon reaching λn , such that λ n < λ dec (see Tables IV-X). Although this linear approach does not attain a globally uniform replica-swap acceptance rate (see [72,73] for more systematic approaches), the resulting scheme suffices to ensure equilibration and convergence, as defined in the next subsection.

The concerted use of replica exchange, fuzzy ensembles, and specialized sampling scheme around the original ensemble provides a sufficient number of independent cavity configurations from the desired ensemble defined in Eq. (H6).

c. Convergence criterion

The quality of equilibration within each cavity is evaluated by monitoring the convergence of two initialization schemes [23,74]: (i) from the original configuration, and (ii) from a randomized configuration prepared by first running 10 6 MC sweeps with shrunk cavity particles at ( λ, α) = (0.5, 20) and then slowly regrowing particles back to λ = 1. We then record configurations every t rec MC sweeps (see Tables IV-X), and monitor the core overlap between new configurations and the original configuration, q on c (t), as a function of the number of MC sweeps, t. The first s eq configurations are discarded, and the average overlap for the following s prod configurations,

q on c ≡ 1 s prod seq+s prod s=seq+1 q on c (t rec s) , (H12)
is computed. Convergence is deemed achieved when the results of both approaches lie within ±0.1 of each other for each cavity. This criterion also allows us to estimate the convergence time for a given value of R and φ. Replica parameters as well as s eq and s prod (see Tables IV-X) are chosen, such that at least 95% out of 40 cavities pass this convergence test, except for the most challenging data point, (φ, R) = (0.640, 3.80), where a 92.5% rate was tolerated. Because the difference between the two approaches is not systematic, however, averaging over 40 cavities results in a rather close agreement between the two schemes, i.e., they converge within ±0.01.

Swap is not enough

It has been suggested that particle identity swaps by themselves suffice to thermalize cavities, especially small ones with R < ξ PTS [35,74]. We here test this hypothesis for our polydisperse system, for which particle swaps are extremely effective at sampling bulk configurations. Figure S8 contrasts the computational time needed to equilibrate cavity configurations at φ = 0.608, with and without parallel tempering. Without replica exchange, the equilibration time rapidly grows as the cavity size decreases, which is in line with the suggestion that confinement enhances the breaking of ergodicity [23]. We observe that when R approaches the point-to-set length from above, however, the convergence time becomes too long to be measured. Tests performed with even smaller R values confirm that trend. Convergence is never achieved for cavities of the order of the point-to-set length or smaller. By contrast, using replica exchange with swap dynamics keeps convergence time within computational reach even for small cavities.

This result seems to contrast with the findings of Ref. [74] that the dynamics actually speeds up inside small cavities. We have no explanation for this discrepancy as we never observed such speed-up in any of our simulations, but we note that the timescale studied in Ref. [74] is different from the one shown in Fig. S8. Specifically, whereas we report a convergence time for the overlap inside the cavity, Cavagna et al. study the time decay of the overlap-fluctuation auto-correlation function in equilibrium. We have not tried to systematically measure this latter auto-correlation because it is already very clear that, according to the former convergence criterion, we cannot properly thermalize small cavities in simulations that use identity swaps alone. Therefore, we would not access an equilibrium correlation timescale but a nonequilibrium one. In addition, no quantitative analysis of the convergence time is provided in Refs. [35,74], which makes a direct comparison with our results impossible.

We conclude that one must generally employ a parallel-tempering scheme, as developed in Ref. [23], in order to properly sample small cavities and measure point-to-set correlations. entropy s * mix independently determined in simulations [29].

Model

We use a continuous size polydispersity, where the particle diameter σ of each particle is randomly drawn from the following particle size distribution: f (σ) = Aσ -3 , for σ ∈ [σ min , σ max ], where A is a normalization constant. We define the size polydispersity as ∆ = σ 2 -2 /σ, where • • • = dσf (σ)(• • • ). We use ∆ = 0.23, choosing σ min /σ max = 0.4492, and set σ as the unit length. We simulate systems composed of N particles in a cubic cell of volume V with periodic boundary conditions in three dimensions, d = 3. We use the following pairwise potential for the soft sphere model,

v ij (r) = v 0 σ ij r n + c 0 + c 1 r σ ij 2 + c 2 r σ ij 4 , (I2)
σ ij = (σ i + σ j ) 2 (1 -|σ i -σ j |), (I3)
where v 0 is the unit of energy, and quantifies the degree of non-additivity of the particle diameters. Non-additivity conveniently prevents crystallization and thus enhances the glass-forming ability of the model. The constants, c 0 , c 1 and c 2 , are chosen so that the first and second derivatives of v ij (r) vanish at the interaction cut-off r cut = 1.25σ ij . We employ the non-additive soft sphere model with parameters n = 12 and = 0.2, and set the number density to ρ = N/V = 1.0 for N = 1500 particles. Equilibrium configurations are produced by swap MC simulations after high temperature configurations are instantaneously quenched to the target temperature. Equilibration is ensured by the fact that particles lose memory of their initial positions and dynamical observables do not present aging. The absence of crystalline nuclei is also verified. We additionally perform standard MC simulations to obtain relaxation times down to the mode coupling crossover temperature T c . Following the procedure described in the Materials and methods section of the main text, we carry out three kinds of fits on these dynamical data. First we employ a power-law fit to extrapolate T c = 0.104. We then perform a VFT fit with δ = 1 and a parabolic fit to estimate the glass-ceiling temperature T g = 0.0720 -0.0817.

Total entropy

We perform the thermodynamic integration over the inverse temperature β from the ideal gas limit β → 0 to the target temperature β,

s tot (β) = 5 2 -ln ρ -3 ln Λ + βe(β) - β 0 dβ e(β ) + s mix , (I4)
where e(β) is the potential energy of the system and for this system we set m = 1 and h = 1 such that the thermal de Broglie wavelength Λ = √ 2πβ. Special care is needed to compute the integral in Eq. (I4) because the potential energy e(β) diverges in the high temperature limit [75,76]. To accurately calculate the integral, we thus decompose the integration range into a very high temperature regime, β ∈ [0, β 0 ], and an intermediate regime, β ∈ (β 0 , β], with β 0 the boundary between the two regimes. Therefore, the integral in Eq. (I4) can be decomposed as where I F and I N are integrals over the very high and the intermediate temperature regimes, respectively. We set β 0 = 1.68 × 10 -3 in this work. The integral I N can be performed by usual numerical integration. To obtain I F we fit the potential energy data to a polynomial function, then analytically integrate the function, which enables us to avoid the numerical integration of the diverging e(β → 0) [75,76]. In a three-dimensional system of particles interacting via v(r) ∝ r -12 , the high temperature expansion of the potential energy reads:

I =
e(β) = Aβ -3/4 + Bβ -1/2 + Cβ -1/4 + O β 0 , (I7)
where A, B, and C are constants. Using Eqs. (I6) and (I7), we get Figure S9(B) shows the results for s tot -s mix . The dashed line corresponds to the Rosenfeld-Tarazona (RT) expression [77] e(T ) ∝ T 3/5 or s tot +3 ln Λ ∝ T -2/5 , which allows us to extrapolate to lower temperatures. Figure S9(B) confirms that the RT expression works very well over our simulation range.

I F = 4Aβ

Vibrational entropy

We compute the vibrational entropy s vib by

s vib = s harm + s anh , (I10)
where s harm and s anh are the harmonic vibrational entropy and the anharmonic correction, respectively [78]. We compute s harm [12,62] using

s harm (β) = 1 N 3N a=1 {1 -ln(βhω a )} IS , (I11)
where • • • IS denotes an average over the inherent structure obtained by the conjugate gradient method, and ω a = λ a /m, and λ a is the eigenvalue of the Hessian. In Fig. S10(A), we show s harm as a function of T . In order to extrapolate s harm to lower temperatures, we follow the common scheme of fitting (1/N ) 3N a=1 {1 -ln(ω a )} IS to a second-degree polynomial in T [12,79]. We also evaluate the anharmonic contribution s anh [78] to the potential energy, where a k is a T independent coefficient. We also assume the linear term in Eq. (I14) is zero, a 1 = 0, which means that the anharmonic contribution to the specific heat vanishes at T = 0. Substituting Eq. (I14) into Eq. (I13), we obtain

s anh (T ) = k=2 k k -1 a k T k-1 . (I15)
In Fig. S10(B), we show a fit of e anh using a 2 and a 3 , and in Fig. S10(A) we show s harm + s anh . We find that the anharmonic contribution is small, |s anh | < 0.08, over our range of interest, and thus that the details of the above procedure have little influence on the final estimate of s vib .

Mixing entropy

We include the effective mixing entropy s * mix obtained by an independent set of simulations in Ref. [29]. In contrast to the hard sphere system, s * mix of the soft sphere system has a weakly linear temperature dependence on the control parameter as shown in Fig. S11. A linear fit, s * mix = b 0 + b 1 T , gives b 0 = 1.3601 and b 1 = 7.6565.

Configurational entropy

The resulting configurational entropy s conf is shown in Fig. S12. The data demonstrates that s conf decreases further below the experimental glass transition, and does not show any sign of a crossover, bending, or saturation. Also, combining the extrapolations for s tot , s harm , s anh , and s * mix , we extrapolate s conf down to zero, and estimate the Kauzmann transition to be around T K 0.04. We conclude from this analysis that the results reported in the main text for hard spheres are not specific to this interaction potential. Our methods thus apply equally well to models of supercooled liquids characterized by continuous pair potentials.

FIG. 1 .

 1 FIG. 1. Breaking the glass ceiling: thermalization beyond the experimental glass transition. (A) Structural relaxation time, τ , for both standard and accelerated swap MC dynamics as a function of the reduced pressure of a polydisperse hard sphere model, where τ is obtained from the decay of the selfintermediate scattering function at a wavenumber k = 5.25. The onset of slow dynamics occurs at Z0 ≈ 18, and the mode-coupling crossover at Zc ≈ 23.5. Times are rescaled by τ0 = τ (Z0) for standard MC. The current limit of colloidal (τ /τ0 = 10 5 ) and molecular (τ /τ0 = 10 13 ) experiments are indicated by vertical bands (this uncertainty stems from the extrapolation scheme), showing that swap MC breaks the glass ceiling. Static structure factor for Z = 18.8 ≈ Z0 (φ = 0.568) and Z = 33.2 (φ = 0.640) for (B) all particles (C) the 40% particles with the smallest diameter. (D, E) show typical snapshots for these two state points, where the smallest particles are highlighted in blue.

FIG. 2 .

 2 FIG.2. Numerical procedures leading to the four estimates of the configurational entropy. (A) Method 1: The Frenkel-Ladd method to obtain the vibrational entropy s vib performs a thermodynamic integration of the mean-squared distance δr 2 between a reference equilibrium configuration and a copy of the system constrained by a harmonic potential of strength α. The integration is carried out from αmax, for which the system behaves as an Einstein solid (indicated by the dashed line δr 2 = 3/(2α)) to αmin, for which particles are trapped by their own cages on the vibrational time scale. (B) Method 2: The numerically-determined Franz-Parisi potential V (Q) is used to measure the configurational entropy as s conf = V (Q high = 0.8) -V (Q low ≈ 0.05). (C) Method 3: The evolution of the overlap Q with the biasing field ε reveals a first-order jump at a value ε for which Q = 1/2 (dashed line). Then, s conf = ε (Q high -Q low ). (D) Method 4: The decay of the cavity overlap correlation function QPTS(R) with cavity radius, R, defines the point-to-set correlation length ξPTS ∝ s

  FIG. 3. Convergent measurements of the four estimates of the configurational entropy beyond the glass ceiling. s conf is plotted as a function of 1/Z ∝ T /P , which is equivalent to the classic Kauzmann plot. All measurements indicate a steep decrease of s conf that continues as the experimental glass ceiling is crossed. The point-to-set estimates are normalized with ξ0 = 2.0 for θ = 2 and ξ 0 = 2.1 for θ = 3/2, respectively, to match the Franz-Parisi estimates at the start of the lowtemperature regime, 1/Z = 0.04 ≈ 1/Zc. The dashed line is an extrapolation based on stot -s vib (see Appendix F). Inset: typical overlap profiles measured in a finite cavity of radius R = 3.46, with colors coding for the overlap value from low (white) to large (black). Overlap fluctuations are uncorrelated around the onset but become strongly correlated over the entire cavity at the largest pressure shown.

  FIG. S1. Equation of state Z = Z(φ) for the hard-sphere system studied in this work. The solid line is the empirical expression for the equation of state from Refs.[52,53], BMCSL, which describes our data well, except at very large volume fractions, where it slightly overestimates Z. A modified version of BMCSL[54] is given by the dashed line.

( 1 )

 1 b (t) as a function of time for the same packing fractions as in panel A.

  FIG.S3. The measured relaxation time, τ , (symbols) is fitted to the VFT form, Eq. (E1), with exponents δ = 1 and δ = 2, the parabolic law, Eq. (E2), and the Arrhenius law, Eq. (E3). The first three fits provide a good description of the data, and are used to locate the glass ceiling at τ /τ0 = 10 13 that is used in the main text.

  FIG. S4. (A)The total entropy, stot, minus the mixing contribution, smix, from simulations and from the (modified) BMCSL EOS. (B) The vibrational entropy, s vib , as a function of φ for several αmin. The solid straight line is fitted to the N = 1000 results with αmin = 7.5 × 10 0 in the high φ region.

  FIG. S5. (A)Mean-squared displacement ∆M between configurations before and after the quench. The crossover, M * , is determined from the intersection of the two linear regimes at small and large M . (B) The mixing entropy smix as a function of M from Eqs. (F10) and (F11). The dashed arrows map M * onto s * mix . The inset shows the continuous polydispersity distribution, f (σ), and the weights of delta functions for the discretized distribution with M = 5, {Xm} m=1,...,5 , as an example.

  FIG. S6.Configurational entropy stot -s vib as a function of volume fraction φ, including the effective mixing entropy smix = s * mix = 1.98. The solid (dashed) curve is a combination of stot from the (modified) BMCSL EOS and s vib from the linear fit.

  FIG. S7. Disorder-averaged probability distribution function of core overlap P (qc), at volume fractions φ = 0.568 (A), 0.597 (B), 0.618 (C), and 0.640 (D) with varying radius R = R1 (black), R2 (green), R3 (blue), and R4 (cyan). We chose (R1, R2, R3, R4)| φ=0.568 = (1.38, 1.73, 2.07, 3.11); (R1, R2, R3, R4)| φ=0.597 = (1.90, 2.25, 2.59, 3.46); (R1, R2, R3, R4)| φ=0.618 = (2.25, 2.59, 3.11, 4.15); and (R1, R2, R3, R4)| φ=0.640 = (2.77, 3.46, 3.80, 5.18).

  FIG. S8. Characteristic convergence time (in unit of MC sweeps) multiplied by the number of replicas with (blue-square) and without (red-circle) replica exchange, as a function of cavity size R. Note that the result without replica exchange at R = 2.42 is an underestimate, because particle swaps alone did not allow the results to converge even after 5 • 10 9 MC sweeps. Convergence is also not achieved for R = 1.73.

  FIG. S9. (A) Plot of β 3/4 e(β) vs. β 1/4 to extract the fit parameters, A, B, and C, in Eq. (I9). The solid line is the resulting fitting curve. The vertical arrow indicates the position of β0. (B) stot -smix as a function of T . The dashed curve is the Rosenfeld-Tarazona expression, stot -smix + 3 ln Λ ∝ T -2/5 , which is well confirmed by our data.

  can compute I F from the fitting parameters, A, B, and C. In order to obtain A, B, and C by fitting, we rewrite Eq. (I7) asβ 3/4 e(β) = A + B β 1/4 + C β 1/4 2 + O β 3/4 ,(I9)and use a quadratic fit as shown in Fig.S9(A).

  FIG. S10. (A) The harmonic vibrational entropy s harm and anharmonic correction s anh . Dashed curves are obtained by quadratic fit to (1/N ) 3N a=1 {1 -ln(ωa)}. (B) The anharmonic contribution of the potential energy e anh (T ). The solid line is obtained by fitting over a range T ∈ [0.0555, 0.125]

e

  anh (T ) = e(T ) -e IS (T )e IS is the inherent structure energy. The last term is the energy of the harmonic vibration. Using e anh (T ), s anh (T ) is given by s anh (T ) = there is no anharmonic contribution at T = 0, i.e., s anh (T = 0) = 0. Expanding e anh (T ) around zero temperature then gives e anh (T ) = k=2 a k T k , (I14)

  FIG. S12. The configurational entropy of soft spheres, s conf = stot -s vib . The glass ceiling is indicated in the blue region and the dashed curve is an extrapolation based on separately fitting the individual terms. The mode-coupling crossover temperature Tc is denoted by a vertical arrow.

TABLE I .

 I Data points for N = 300 used in the Franz-Parisi construction (Methods 2 and 3).

	17.6	19.8	22.3	24.8	27.0	28.1	30.4	32.0	34.9
	φ 0.559 0.573 0.588 0.603 0.612 0.619 0.628 0.635 0.645
	Z 23.7	24.6	25.8	28.3	30.7	32.0	34.0	37.6
	φ 0.598 0.602 0.609 0.619 0.630 0.635 0.643 0.655

TABLE II .

 II Data points for N = 1000 used in the thermodynamic integrations (Methods 1).

	Z 18.8	21.7	23.5	25.6	27.7	30.4	33.2
	φ 0.568 0.587 0.597 0.608 0.618 0.629 0.640

TABLE III .

 III Data points for N = 8000 used in the thermodynamic integration computations (Methods 1) and the point-to-set correlation (Method 4). FIG. S2. Panel A: Collective overlap time correlation function Qc(t) for several packing fractions (from left to right, φ = 0.568, 0.609, 0.640, 0.655). Panel B: Probability of losing all but one neighbor P

		Appendix D: Equilibration
		1		1
	Collective overlap correlation fct., Q c (t)	0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0 10 A 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 0.655 = φ 0.640 0.609 0.568	(t) b (1) Prob. of losing all but 1 neigbor, P	10 0 10 0 0.2 0.4 0.6 0.8 B 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 0.655 = φ 0.640 0.609 0.568
		Time, t [MC sweeps]		Time, t [MC sweeps]

TABLE IV .

 IV Cavity PTS measurement parameters φ = 0.568. Runs without parallel tempering have nave = 1.

	R	1.38 1.73 2.07 2.42 2.77 3.11 3.46 3.80 4.15
	nave	10	8	7	7	6	1	1	1	1
	λ dec 0.850 0.920 0.960 0.970 0.980 NA NA NA NA
	trec	4000 4000 4000 4000 4000 1000 400 200 200
	seq	500 500 500 500 500 500 500 500 500
	s prod 2000 2000 2000 2000 2000 2000 2000 2000 2000
	R	1.73 2.07 2.42 2.77 3.11 3.46 3.80 4.15 4.49
	nave	10	10	10	8	7	1	1	1	1
	λ dec 0.900 0.930 0.950 0.970 0.980 NA NA NA NA
	trec	2 • 10 4 2 • 10 4 2 • 10 4 10 4 4000 2000 400 400
	seq	500	500	500 500 500 500 500 500 500
	s prod 2000 2000 2000 2000 2000 2000 2000 2000 2000

TABLE V .

 V Cavity PTS measurement parameters φ = 0.587. Runs without parallel tempering have nave = 1. • 10 4 2 • 10 4 2 • 10 4 10 4 10 4 10 4 4000 2000 1000

	R	1.90 2.25 2.59 2.94 3.11 3.46 3.80 4.15 4.49
	nave	12	13	12	11	9	1	1	1	1
	λ dec 0.900 0.920 0.940 0.960 0.970 NA NA NA NA
	trec 2 seq	500	500	500 500 500 500 500 500 500
	s prod 2000 2000 2000 2000 2000 2000 2000 2000 2000

TABLE VI .

 VI Cavity PTS measurement parameters φ = 0.597. Runs without parallel tempering have nave = 1.

	R	2.07 2.42 2.77 3.11 3.46 3.80 4.15 4.49 4.84
	nave	13	14	12	12	1	1	1	1	1
	λ dec 0.910 0.920 0.950 0.960 NA NA NA NA NA
	trec	2 • 10 4 2 • 10 4 2 • 10 4 10 4 2 • 10 4 10 4 4000 2000 2000
	seq	500	500	500 500 500 500 500 500 500
	s prod 2000 2000 2000 2000 2000 2000 2000 2000 2000

TABLE VII .

 VII Cavity PTS measurement parameters φ = 0.608. Runs without parallel tempering have nave = 1.

	R	2.25 2.42 2.59 3.11 3.80 4.15 4.49 4.84
	nave	15	16	18	17	1	1	1	1
	λ dec 0.900 0.910 0.910 0.940 NA	NA NA NA
	trec	10 4 10 4 10 4 10 4 5 • 10 4 3 • 10 4 10 4 5000
	seq	2000 2000 2000 1000 500	500 500 500
	s prod 8000 8000 8000 4000 2000 2000 2000 2000

TABLE VIII .

 VIII Cavity PTS measurement parameters φ = 0.618. Runs without parallel tempering have nave = 1.

	R	2.77 2.94 3.11 3.46 3.80 4.49 4.84
	nave	17	19	18	17	16	1	1
	λ dec 0.930 0.930 0.940 0.955 0.965 NA NA
	trec	10 4 10 4	10 4 10 4 10 4 10 5 4 • 10 4
	seq	1500 3000 4000 6000 8000 500 500
	s prod 6000 12000 16000 9000 12000 2000 2000

TABLE IX .

 IX Cavity PTS measurement parameters φ = 0.629. Runs without parallel tempering have nave = 1.

	R	2.77 3.11 3.46 3.63 3.80	4.84 5.18
	nave	21	21	22	23	22	1	1
	λ dec 0.910 0.930 0.940 0.940 0.950 NA	NA
	trec	10 4 10 4 10 4	10 4	10 4 2 × 10 5 10 5
	seq	2000 6000 8000 12000 12000 1000 1000
	s prod 8000 9000 12000 18000 18000 4000 4000

TABLE X .

 X Cavity PTS measurement parameters φ = 0.640. Runs without parallel tempering have nave = 1.

  The effective mixing entropy s * mix . The dashed line is an empirical linear fit.

		2.4
		2.3
		2.2
	s* mix	2 2.1
		1.9
		1.8
		0.04 1.7	0.06	0.08	0.1	0.12	0.14
				T
	FIG. S11.

ACKNOWLEDGMENTS

The research in Montpellier was supported by funding from the European Research Council under the European Unions Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement No 306845 and by a grant from the Simons Foundation (#454933, Ludovic Berthier). The research at Duke was supported by a grant from the Simons Foundation (#454937, Patrick Charbonneau) and associated computations were carried out through the Duke Compute Cluster. Data relevant to this work have been archived and can be accessed at http://dx.doi.org/10.7924/G8ZG6Q9T.

Appendix C: Equation of state

We report in Fig. S1 the equilibrium equation of state, Z = Z(φ) (from Ref. [20]), of the system under study, thus enabling a translation of the reported results from Z to φ. The specific data points used in the four measurements of the configurational entropy are also presented in Tab. I, II, and III.

Appendix E: Locating the glass ceiling

In the context of this work, we refer to the typical relaxation time measured at the laboratory glass transition (at which, conventionally, τ /τ 0 = 10 13 [14]) as the glass ceiling, Z g . Because standard MC dynamics can only access Appendix G: Methods 2 and 3: The Franz-Parisi free energy

In this section we detail the Franz-Parisi construction [32]. We describe the order parameter of the glass transition, the umbrella sampling, as well as the parallel tempering and histogram reweighting techniques used to measure numerically the configurational entropy defined by Methods 2 and 3. Because the presentation closely follows earlier reports in which the Franz-Parisi free energy was determined numerically [13,67,68], we only provide a brief summary of the computational methodology used.

Order parameter

The order parameter of the Franz-Parisi construction is the overlap Q (or similarity) between two disordered configurations, defined as

where θ(x) is the Heaviside function, and r 1,i and r 2,j are the positions of particle i and j within configurations 1 and 2, respectively. By definition, Q is large when the density fields in configurations 1 and 2 are very close, and Q is close to zero when the density profiles are uncorrelated. The length, a = 0.23, accounts for small-amplitude thermal vibrations, so that density profiles that only differ by local vibrations caused by thermal fluctuations have Q close unity.

Umbrella sampling

First, a reference configuration, {r 2 }, is chosen at random from the equilibrium states of the system at a fixed volume fraction. In order to access the probability distribution function of the overlap Q 12 between configurations {r 1 } and {r 2 }, we conduct n distinct simulations. In each simulation, = 1, • • • , n, system 1 evolves according to the Hamiltonian

where H({r}) is the original hard-sphere Hamiltonian and ε is the thermodynamic field conjugate to the overlap Q 12 .

The biasing potential W (Q) is taken to have the form

with parameters (k , Q ) chosen to constrain the overlap Q 12 to explore values away from its average equilibrium value. We then perform local displacements and swap MC moves, which are accepted using a Metropolis acceptance rate given by the Hamiltonian (G2). Provided that the system is properly thermalized, our measurements yield the equilibrium probability distribution function of the overlap,

where • • • denotes the thermal average with Hamiltonian H in Eq. (G2) at fixed reference configuration 2. Because we perform a quenched average to compute the overlap distribution, however, a second averaging, over the reference configuration 2, is needed to determine the Franz-Parisi free energy. The idea behind the biasing potentials W (Q) in Eqs. (G2), and (G3) is that the fluctuations of the overlap in each simulation can be tailored to explore a narrow region centered around Q . Each simulation thus explores but a small range of overlap values, and it becomes unnecessary to wait for very rare overlap fluctuations to take place. Umbrella sampling enables the efficient measure of the atypical overlap fluctuations needed to reconstruct the Franz-Parisi free energy.

Parallel tempering

While the umbrella sampling technique described above considerably accelerates the measurement of overlap fluctuations, we have observed that when φ is high, N is large, and/or Q is large, the particle dynamics still slows down We compute the configurational entropy of a continuous polydisperse soft sphere system by the standard definition,

where s tot and s vib are the total and vibrational entropies [12,62]. The mixing entropy is also treated by the strategy used in Ref. [29]; we compute s tot and s vib as if the system were monodisperse, and then add the effective mixing