Efficient measurement of point-to-set correlations and overlap fluctuations in glass-forming liquids

A well-known difficulty in understanding the increasing sluggishness of glass-forming liquids upon lowering temperature T is that no obvious change to the static structure of the liquid accompanies it.¹ In stark contrast to the critical slowing down observed near a standard critical point, static correlation functions of glass formers barely budge while the structural relaxation timescale grows by more than 15 orders of magnitude. The glassy slowdown has thus instead been attributed to a different type of criticality, one at which the free-energy landscape becomes very rugged, leading to the emergence of many metastable states separated by growing free-energy barriers.^2–4 

In order to capture this ruggedness, the concept of a point-to-set (PTS) correlation was introduced about a decade ago.^5,6 PTS correlations generalize multi-point correlations to their infinite-point limit, which makes them sensitive to the structure of the full configuration. Various PTS geometries have since been considered, including cavity,^7–11 random pinning,^12–18 and walls.^19,20 Despite subtle differences between them,^21,22 the overarching observation remains. For glass formers the PTS correlation length grows more rapidly than lengths extracted from two-point correlation functions^{8,10,12,20,21} and is thus an important aspect of glassy physics. However, extracting PTS correlations from numerical simulations has thus far been limited by the computational challenge of measuring them.^21,22

In order to better understand the nature of this challenge, let us consider the case of PTS correlation within the cavity geometry. First consider a liquid at equilibrium, specify a finite cavity, say a sphere of radius R, and pin everything outside of it — thus fixing a set of particles. Then allow the particles inside the cavity to explore phase space under the constraints exerted by the set of pinned particles, now acting as an effective quenched disorder. Particles inside the cavity thus explore the local free-energy landscape. A standard way of characterizing these local landscapes, adapted from the study of spin glasses, is to consider the statistics of the overlap, $q r$⁠, between two equilibrium configurations at a point r inside the cavity. The materials we consider are characterized by a very slow dynamics in the bulk, and confining them with amorphous boundaries is found to make the natural dynamics of these systems even slower.^20–23 The major computational problem is thus that of properly sampling the equilibrium fluctuations of the confined fluid.

The origin of the dynamical slowdown may be attributed to the confining amorphous boundaries provoking an entropy crisis qualitatively similar to the one which could be happening at the Kauzmann temperature, T_K, in the bulk liquid.²⁴ Increasing the confinement would reduce the number of accessible states from being exponentially large in the number of particles for large cavities, to being sub-exponential for small cavities.⁵ This shift would then be evidenced in the evolution of overlap statistics in the local landscapes, which is exactly what PTS correlations purport to capture, and its location defines the PTS correlation length, ξ_PTS(T). From this perspective, the computational difficulty encountered in previous studies for R ≲ ξ_PTS should not come as a surprise. Confinement should be inducing the analog of an equilibrium (or “ideal”) glass transition in a system containing a finite number of particles.⁵ Proper sampling then requires very large computer resources, as recent studies of constrained systems have evidenced.^17,18,25 As a result, the entropy crisis picture has only been partially validated. In particular, although solid evidence for a growing PTS correlation length has been obtained in various systems,^8–12 several questions remain unexplored regarding the nature of the crossover between high and low overlap regimes, the associated fluctuations, and the connection to dynamical relaxation in the bulk.

In this paper, we develop a generic computational method to more quickly and reliably sample PTS correlations. As a proof of principle, we study the Kob-Andersen binary Lennard-Jones model,^26,27 which is a classical glass-forming liquid. Yet our approach is sufficiently generic to be applied to other glassy systems, including hard spheres. We use the efficiency of our approach to record novel physical observables, beyond the PTS correlation function, that rely on efficiently sampling thermodynamic fluctuations inside a cavity. In particular, we follow the evolution with cavity size and temperature of the complete probability distribution function of the overlap. Its variance defines an overlap susceptibility that allows us to directly locate the PTS crossover scale without any empirical fitting. In addition, we are able to spatially resolve overlap fluctuations inside the cavity, giving us access both to their radial profile and their orthoradial fluctuations, and thus to make connection with PTS measurements in other geometries, such as flat walls.

The plan of this paper is as follows. In Sections II and III, we detail the simulation model and the parallel-tempering methodology, respectively. In Section IV, we present the approach for computing PTS correlations as well as the results, and in Sec. V, we evaluate results for the wall and wandering lengths from the cavity PTS setup. We conclude in Section VI.

We study the behavior of a glass-forming liquid first proposed by Kob and Andersen.^26,27 The Kob-Andersen binary Lennard-Jones (KABLJ) model contains two particle species, denoted A and B, with equal mass m, and interacting via the pair potential

where α, β ∈ {A, B} with parameters ε_AB/ε_AA = 1.5, ε_BB/ε_AA = 0.5, σ_AB/σ_AA = 0.8, and σ_BB/σ_AA = 0.88. The interaction potential is cut off at $r α β cut =2.5 σ α β$ and shifted, such that it vanishes at the cutoff. The relative number of particles is N_A : N_B = 4 : 1 and the overall number density is $ρ=N/V=1.2 σ A A − 3$⁠. Note that we report below lengths and temperatures in standard dimensionless Lennard-Jones units set by σ_AA and ε_AA/k_B.

We use bulk samples with N = 135 000 particles in a periodic cubic box, generated as described in Ref. 28, at temperatures T = 1.00, 0.80, 0.60, 0.51, 0.45. At each temperature, we take 50 equilibrated snapshots, separated by more than 2τ_α, where τ_α is the bulk structural relaxation time. Then, within each snapshot, we randomly select a position as cavity center. For each center, we fix the particles outside the cavities of radius R and continue sampling configurations for the N_cav particles inside the cavity, as described in Sec. III.

Note that in order to directly evaluate how the behavior of the cavity changes as function of R, we preserve each of the cavity centers as R increases. Note also that the system size used here is far larger than necessary, because the length scales we extract in Sections IV and V are rather small. Sampling 50 cavity centers within a single snapshot would have sufficed to ensure statistical independence of the resulting cavities and would have thus provided similar results.

A key hurdle to measuring PTS observables in simulations has been the large computational effort needed to sample different equilibrium configurations inside a cavity. Here, we design a Monte Carlo parallel-tempering²⁹ scheme (a type of Hamiltonian replica exchange³⁰) that sidesteps this difficulty by dramatically lowering barriers within the rugged free-energy landscape. The approach goes as follows. We prepare a = 1, …, n replicas of a given configuration inside the cavity. The replicas evolve at different temperatures T_a and “shrinkage” parameters λ_a ≤ 1 (see the  Appendix for actual values) with a deformed Hamiltonian

where $λ ̃ a = λ a$ for a pair of mobile particles within a cavity and $λ ̃ a = 1 + λ a 2$ for a pair containing one mobile (inside the cavity) and one pinned particle (outside the cavity). We draw cavity configurations from the bottom replica with T₁ = T and λ₁ = 1, whereas the other replicas evolve at higher temperatures T_a > T and with smaller particles, λ_a < 1, in order to speed up their dynamics.

Within each replica, we perform simple Monte Carlo (MC) moves that consist of: (i) choosing a particle i from N_cav mobile particles inside the cavity; (ii) displacing particle i by $Δx=l n ˆ$⁠, where l is uniformly drawn from [0, 0.3] and $n ˆ$ uniformly drawn on the sphere S²; and (iii) accepting/rejecting the displacement according to the Metropolis criterion. Note that we put a hard spherical wall at the edge of the cavity, so that all moves that take a mobile particle outside the cavity radius are rejected. Each MC sweep consists of N_cav MC trial moves, hence on average each particle attempts to move once. More crucially, in order to release the disorder constraint (frustration) induced by the pinned particles, identity-exchange of a pair of adjacent replicas [in (T_a, λ_a) space] is attempted every 1000 MC sweeps, on average. This replica swapping is again accepted or rejected according to the Metropolis criterion, so that our MC algorithm ensures proper equilibrium sampling. Compared to earlier schemes for cavity studies, the proposed method is distinct from the simpler annealing procedure used before for the same model¹⁰ and is more generally applicable than the local particle swap Monte Carlo moves that are only efficient for specific glass-forming models.⁸ 

The quality of the equilibration within each cavity is evaluated by employing two schemes, following the approach developed by Cavagna et al.:²² we start the system from (i) the original configuration and (ii) a randomized configuration prepared by putting the cavity at T = 1.00 and λ = 0.6 for 10⁴ MC sweeps. We then record the re-equilibrated configurations every t_rec = 10⁴ MC sweeps and monitor the core overlap (as defined in Sec. IV) between the new configuration and the original cavity configuration, $q c on t$⁠, as a function of the number of MC sweeps, t. Their running average,

decreases for the first scheme and increases for the second, converging upon equilibration (Fig. 1). The first s_eq configurations are discarded, and the overlap for the following s_prod,

is computed. Convergence is deemed obtained when the results of both approaches lie within ±q_tol of each other for each cavity. Replica parameters as well as s_eq and s_prod are chosen, such that for q_tol = 0.1, at least 98% of cavities pass this convergence test (see the  Appendix for actual values). However, because the difference between the two approaches is not systematic, averaging over 50 cavities results in a very close agreement between the two schemes, namely in a convergence of $q ̄ ( t )$ within ±0.005.

This algorithm behaves differently in different cavity size regimes. The most impressive results are obtained for small cavities, even at low temperatures. Only a few replicas are indeed needed for a small system to jump over the relatively high free-energy barriers and thus quickly equilibrate and properly sample configurations. This efficiency is illustrated in Fig. 1, which shows that the algorithm time to convergence is smaller for smaller cavities. This outcome markedly contrasts with simpler MC schemes, the sampling efficiency of which is akin to that of molecular dynamics. As an example, results from Ref. 10 suggest that for low temperatures T ≤ 0.51 and small radii R ≤ 3.0, at least 10¹⁰ MC sweeps would be needed to equilibrate the system. Here, by contrast, we can equilibrate a similar system within ∼10⁷ MC sweeps using ∼10 replicas. (Note, however, that the computational efficiency of our algorithm depends sensitively on the details of the parallel-tempering parameters, which require more fine-tuning than the parameters of simpler schemes.) We thus conservatively achieve at least a 100-fold speedup, and that speedup would be even stronger for state points where simpler MC schemes have not been found to converge. By contrast, for R ≫ ξ_PTS, confinement effects are negligible and even simple Monte Carlo moves suffice to sample configurations. The most problematic regime is R ∼ ξ_PTS at low temperatures, where even relatively large cavities display fine-tuning bottlenecks. That this regime sets the lower limit on the convergence criterion is unsurprising, as the optimal spacing between adjacent replica parameters generically scales as $1/ N cav$⁠. We thus need $∼ ξ PTS 3 / 2$ replicas to properly equilibrate a sample at this most resilient regime.

In earlier work,²² an acceleration of the overlap dynamics was reported in a binary soft sphere mixture, where traditional Monte Carlo moves are supplemented by binary particle exchanges (or “swap” moves), a method which was used in a series of numerical works.^7,8,22 We cannot directly compare our work to these dynamic measurements, as our parallel-tempering scheme does not allow us to measure time correlation functions at fixed temperature, and no evolution of the convergence time was reported using particle swaps only.²² We believe, however, that the speedup of the convergence time that we report in Fig. 1 is of a different nature and should uniquely be attributed to the merits of the parallel-tempering scheme and not to any natural dynamical process taking place in the glass-forming liquid at various degrees of confinement. Methods developed in Ref. 31 could be used in the future to quantify more precisely the convergence time of the parallel-tempering scheme using time-correlation functions associated to the random walk of the various copies in the (T_a, λ_a) space.

Once the glass-forming liquid is confined by amorphous boundaries in a finite cavity, we need to analyse the equilibrium fluctuations of the overlap field q(r) inside the cavity, which is the principal observable used in this paper. To this end, we proceed as follows. Denote a pair of configurations by $X= x i$ and $Y= y i$⁠. For each particle x_i, find the nearest particle y_inn of the same species and assign an overlap value $q X;Y x i ≡w | x i − y i nn |$⁠, where

with b = 0.2. This function defines overlap values $q X;Y x i$ at scattered points $x i$⁠. We then define $q X;Y r$ to be a continuous function precisely passing through these points. Specifically, we first perform a Delaunay tessellation of space and, to a point r within a simplex spanned by four points $x i i = i 1 , i 2 , i 3 , i 4$⁠, we associate a linearly interpolated value

where $c i i = i 1 , i 2 , i 3 , i 4$ satisfies $r= ∑ i = i 1 , i 2 , i 3 , i 4 c i x i$ with $∑ i = i 1 , i 2 , i 3 , i 4 c i =1$⁠. Similarly, we obtain $q Y ;X r$ and define

which provides the overlap near the core of the cavity

where r_c = 0.5 and r = 0 is the cavity center. This integral is numerically evaluated by Monte Carlo integration with 10⁴ points. Note that our results are insensitive to both the choice of overlap function w(z) as long as its range is kept of the order of cage size and to distinguishing or not between species. It is, however, important for the overlap field, $q X,Y r$⁠, to be continuous in order to attain good spatial resolutions within a cavity (see Sec. V). We thus avoid using a step function.

While we monitored the overlap between the original configuration and a re-equilibrated configuration in checking for a good equilibration, for the rest of the paper, we evaluate the overlap between two statistically independent re-equilibrated configurations. Overlaps are computed by comparing two configurations obtained from two different schemes, which guarantees their statistical independence. In the following, we denote 〈…〉_J(R), the thermal average inside the cavity with the effective quenched disorder J(R) set by a pinned external configuration, evaluated by averaging over s_prod pairs of configurations, and […] the average over disorder, evaluated by averaging over 50 disorder realizations.

As alluded to in the introduction, we expect a crossover from high to low core overlap around the cavity PTS length scale. The PTS correlation function,

indeed decays on a characteristic length scale, ξ_PTS, that clearly grows as temperature decreases, as shown in Fig. 2(a). This behavior is qualitatively consistent with earlier work^8,10 and reveals a rarefaction of metastable states leading to the growing PTS length ξ_PTS, as originally envisioned by Bouchaud and Biroli.⁵ 

The standard way to quantify ξ_PTS has thus far been to fit the decay of $G PTS R$ to a compressed exponential form $G PTS ∼exp − ( R / ξ ) η$⁠, where η is an adjustable exponent. When reliable information about small cavity is available, however, this ad hoc fitting form is not computationally optimal. In order to obtain a more robust determination of the PTS length scale, we can instead appeal to another important aspect of the PTS physics: the appearance of an increasing number of accessible configurations as the cavity size increases. Existence of these competing states results in configurational fluctuations, in addition to vibrational fluctuations around these configurations. Hence, the overlap function can take a larger number of possible values, and their relative contribution depends on the probability with which alternate configurations can be accommodated by a given cavity. In the large cavity limit, most configurations are different from one another and the overlap is very small with very little fluctuations. In the interesting crossover regime, however, we expect the probability distribution function (PDF) of core overlaps to display a very broad range of fluctuations, while for small cavities configurations are all very similar, hence the overlap is high and does not fluctuate much.

These expectations are validated by our numerical results, as illustrated in Fig. 3. For each cavity and temperature, we find the PDF to be narrow at large and small R, with a maximum in the width of the fluctuations for an intermediate cavity size. For some (though not all) cavities, we find that the PTS crossover size corresponds to a bimodal distribution of overlap values, e.g., Fig. 3(b). As temperature is lowered, we find that the fraction of bimodal distributions increases, which results in a disorder-average PDF that is itself nearly bimodal, e.g., Fig. 3(c). We expect that at even lower temperatures, most PDF should become bimodal, so that the full disorder-averaged distribution then also becomes bimodal. Bimodal PDF of the overlap have indeed been found in other constrained systems.^{17,18,25,33,34} We note in passing that a bimodal distribution of overlaps in a cavity is an interesting qualitative signature of the correspondence between the PTS crossover length and a crossover analogous to an equilibrium glass transition in a finite size system.

The shape of the PDF encodes all the information we need about the PTS correlation length. In particular, the breadth of the distribution can be easily detected from the PTS susceptibility

The numerical results in Fig. 2(b) show that this function has a clear non-monotonic evolution with R, with a peak that gives an unambiguous, fitting-form-free definition of ξ_PTS. In addition, this approach to determining ξ_PTS is computationally advantageous because it does not require a careful fit of the PTS correlation tails, where equilibrium sampling is slower within our Monte Carlo scheme, and where the overlap value is small (and hence difficult to precisely measure). Using the susceptibility χ_T instead allows us to estimate the PTS length scale from only a handful of very precise measurements around the susceptibility maximum. Given the limited size of the PTS length for most of the computationally accessible bulk temperatures and the efficiency of parallel-tempering for small cavities, this approach offers an especially effective scheme for determining the PTS length in generic model glass-forming liquids.

It is comforting to note that the peak location of the PTS susceptibility, $ξ PTS peak$⁠, roughly agrees with that obtained through the compressed exponential fit, $ξ PTS fit$⁠. The correspondence is also very good with $ξ PTS th$⁠, defined as $G PTS ξ PTS th ≡ q th$ with q_th = e⁻¹. The quantitative agreement between these three definitions is demonstrated in Fig. 2(c), which shows that all three estimates give comparable results, except at the (uninteresting) highest temperature studied, T = 1.0, where the very definition of a PTS length becomes somewhat problematic.²⁵ The similarity between the different approaches allows us to provide a reasonable estimate of ξ_PTS at the lowest temperature considered, T = 0.45, which is fairly close to estimates of the mode-coupling crossover temperature for this model, T_MCT ≈ 0.435.²⁶ At that temperature, the susceptibility peak occurs at cavity sizes that are slightly beyond the computational reach of our algorithm. The correlation decay nonetheless provides a reasonable estimate of ξ_PTS.

Although there are systematic quantitative differences in ξ_PTS obtained by the current approach compared with earlier results, these can be explained by the choice of overlap function and other algorithmic details. Qualitatively, nothing much differs. The representation of the temperature dependence of ξ_PTS in Fig. 2(c) emphasizes the analogy between cavity PTS measurements and other ways of constraining the available phase space of the system, which are being actively investigated, such as random pinning^13–18 or coupling to a reference configuration.^33–38 In all cases, an external constraint with quenched disorder allows one to induce a sharp transition (pinning, coupling) or crossover (cavity) between an unconstrained phase with low overlap and a localized phase with large overlap. The data in Fig. 2(c) thus resemble qualitatively published phase diagrams for the physics of constrained glass-forming liquids.^17,18,33,34 Note that, in all of these approaches, one is interested in understanding the emergence of metastable states in the physically relevant temperature regime much above the putative Kauzmann transition, whose existence/absence is therefore completely immaterial.

The availability of well-averaged configurations for cavities at R ≈ ξ_PTS enables us to have a closer look into the PTS correlation physics. We specifically study two length scales associated with the spatial structure of overlap profiles under pinning. To define them, let us first consider the PTS physics in the limit R/ξ_PTS → ∞, where configurations are frozen beyond a flat wall. The first length scale of interest is the wall (or penetration) length, $ξ ∞ wall$⁠, which corresponds to the characteristic distance from the wall for the decay of the overlap.^19,39,40 The spatial fluctuations of this decay length also result in a wandering of the iso-overlap surface. The second length scale is thus the wandering length, $ξ ∞ ⊥$⁠, which characterizes the size of these fluctuations in the direction orthogonal to the wall.⁴¹ A priori, there is no relation between the growth of these two scales. The wandering length could even decrease as temperature decreases. Recent theoretical work,⁴¹ however, suggests that these two length scales grow at the same rate, albeit remain subdominant to the cavity PTS correlation length, ξ_PTS. Here, we analyze avatars of these two length scales, $ξ ⋆ wall$ and $ξ ⋆ ⊥$⁠, in cavities of size R ≈ R_⋆ ≡ ξ_PTS. We emphasize that these avatars are distinct from their counterparts in the flat wall limit, although they are smoothly connected to $ξ ∞ wall$ and $ξ ∞ ⊥$ as functions of R/ξ_PTS.

To extract the wall length, we first consider the radial dependence of the overlap profile within a cavity, $[ 〈 q r 〉 J ( R ) ]$⁠, which can be obtained by binning the overlap profile away from the cavity center (Fig. 4(a)). When the cavity size R ≈ ξ_PTS, nonlocal constraints imposed by disorder is felt throughout the full cavity, via the rarefaction of metastable states in the local free-energy landscapes. The near collapse observed in the inset of Fig. 4(a) suggests that the extent of the pinning effect in this regime, $ξ ⋆ wall$⁠, is proportional to ξ_PTS. When confined on the scale of ξ_PTS the system thus appears fully correlated, as if it were at a sharp transition point. By contrast, the data in Fig. 4(b), obtained for R ≈ 2ξ_PTS, show that the penetration length measured in a geometry that interpolates between the finite cavity and the flat wall seems to have a milder temperature dependence than ξ_PTS. This behavior is as expected, given that $ξ ∞ wall$ is observed to grow subdominantly to ξ_PTS.

To analyze the avatar of wandering length, we consider the orthoradial fluctuations of the overlap field inside the cavity. Recent field-theoretic calculations suggest that the overlap profile $〈q r 〉 J ( R )$ for each individual cavity is aspherical,^41,42 even though it appears spherically symmetric when averaged over the disorder of the cavity pinning. A hint of this effect can be seen in Fig. 5 by examining the iso-overlap surfaces,

which are generically rather bumpy, even after smoothing particle-scale fluctuations. In order to quantify asphericity, we calculate $〈q r 〉 J ( R )$ on 10⁶ points uniformly distributed within a cavity. We then consider a point to be part of the iso-overlap surface at q_iso if it falls within $[ q iso − Δ q 2 , q iso + Δ q 2 ]$ with Δq = 0.01 and denote the resulting set of points as $r i i = 1 , … , n$⁠. The variance, $χ ⊥ q iso J ( R )$⁠, of distances to the edge of the cavity, $R − | r i | i = 1 , … , n$⁠, directly quantifies the sample-to-sample deviations from spherical iso-overlap surfaces. The average of this quantity gives an estimate of the distance between the iso-overlap surface and the cavity wall, $R−r q iso J ( R )$⁠. After averaging over the pinning disorder, we obtain

The corresponding numerical results for this quantity at R ≈ ξ_PTS are shown in Fig. 6. Near the edge of the cavity, χ_⊥ is small. As roughness accumulates toward the core, χ_⊥ grows until the iso-overlap surface becomes so rough that, as it is defined, χ_⊥ does not approximate (ξ^⊥)² well anymore. Note indeed that in Fig. 5(c) the surface even becomes topologically distinct from a sphere. The collapse in the inset of Fig. 6 suggests again the extent of orthoradial fluctuations of the overlap in this regime, $ξ ⋆ ⊥$⁠, may be collapsed by the PTS length scale.

On the whole, these results suggest that the spatial structure of overlap profiles in this regime is governed by a single length scale, ξ_PTS, although a conclusive interpretation of this effect has yet to be reached. It would be interesting to examine how this structure morphs into in the flat wall geometry limit, where $ξ ∞ wall$ and $ξ ∞ ⊥$ are expected to grow subdominantly to ξ_PTS. The largest cavities we equilibrated in this paper are of the order of 2ξ_PTS, and our preliminary analysis did not conclusively determine whether they belong to the flat wall regime or not, in particular for the wandering length. A study of the evolution of the overlap profiles for a broader range of cavity sizes, including the careful analysis of the wandering length in the wall geometry, should be performed in future work.

The results of this work suggest that parallel tempering is a method of choice for studying PTS correlations in generic model glass-forming liquids. This scheme has indeed allowed us to quantify PTS correlations in a canonical glass former over the temperature regime that is generically studied in the bulk, bypassing the slowdown due to confinement by amorphous boundaries. This scheme thus provides an efficient solution to a major computational obstacle that has prevented progress in studies of static correlations in glass formers for nearly a decade.

In addition, by defining a novel overlap susceptibility for finite cavities, we have obtained PTS correlation information well beyond the regime typically attained in simulations, and from a reduced computational effort. Our results are broadly consistent with a growing static cavity PTS length scale as the temperature of a glass-forming liquid is lowered, and thus with glassiness being associated with a roughening of the free energy landscape.

As expected from extrapolating earlier PTS results, however, the measured static length scale does not grow by much more than a factor of two upon approaching the mode-coupling crossover temperature. In Fig. 7, we replot our data for ξ_PTS from Fig. 2, and compare its temperature evolution with several dynamic correlation lengthscales. Given apparent subtleties in extracting dynamical length scales, we record three independent sets of data taken from Refs. 20, 43, and 44, taken from either measurements of four-point dynamic correlations^43,44 or from dynamic profiles near a wall²⁰ for the same model. This comparison shows that all three dynamic lengths grow more rapidly than the static one toward low temperatures, but the difference is relatively modest. A stronger decoupling between static and dynamic quantities was observed in hard-sphere systems^12,21 than what we find in the KABLJ model. The distinction between the two models is in harmony with previous findings that signatures of the mode-coupling crossover (associated to a rapid growth of a dynamic correlation length) appear weaker in the KABLJ model than in hard or quasi-hard spheres.^20,45

While there have been recent attempts at experimentally measuring PTS correlations in colloidal glass-formers,^46,47 corresponding studies for molecular glass-formers, roughly 10 orders of magnitude more sluggish, still remain inaccessible. Analysis of such extremely sluggish systems would be particularly helpful in properly assessing the theoretical framework that surrounds PTS correlations. In future work, we will thus couple the present approach with other enhanced sampling techniques in an attempt to push back the computationally accessible boundary for the study of glassiness. Such studies will also be useful for systematically investigating the behavior of various length scales ξ_PTS, $ξ ⋆ , ∞ wall$⁠, and $ξ ⋆ , ∞ ⊥$ and how they relate to one another deep in the dynamically sluggish regime.

We acknowledge many fruitful exchanges and correspondence with G. Biroli, C. Cammarota, G. Hocky, T. Kawasaki, V. Martin-Mayor, D. Reichman, G. Szamel, and G. Tarjus. We also acknowledge the Duke Compute Cluster for computational resources. The research leading to these results has received funding from the European Research Council under the European Union Seventh Framework Programme (No. FP7/20072013)/ERC Grant Agreement No. 306845 (L.B.). P.C. and S.Y. acknowledge support from the National Science Foundation Grant No. NSF DMR-1055586 and the Sloan Foundation.

The lists of replica parameters used for each temperature and radius explored in this paper, along with s_eq, s_prod, and equilibration success rate are provided in Tables I–V. In general, we choose $( T a , λ a ) a = 1 , … , n$⁠, such that they satisfy (when n > 1) the linear relation

where (T_dec, λ_dec) is chosen such that the particles can decorrelate the overlap effectively. We always choose the last replica parameter λ_n ≥ λ_dec.

The recording time is fixed as t_rec = 10⁴ MC sweeps. As described in the main text, we discard the first s_eq configurations, keep the following s_prod configurations, and compare the thermal averages [computed with Eq. (4)] obtained from two schemes: one starting from the original configuration and the other from a randomized configuration. The convergence criterion requires that results from both approaches lie within ±q_tol of each other. For each data point, we record as success rate how many of 50 cavities satisfy this criterion for a particular q_tol. Globally, for q_tol = 0.1, 98% of cavities are deemed to have converged, for all temperatures and radii.

Note that the parameters chosen are neither unique nor optimal. Some of the chosen parameters result in near bottlenecks in the exchanges of replicas. Tuning the replica parameters for each cavity by hand could certainly help sampling low-temperature configurations. A more promising and general way would be to algorithmically tune the replica parameters by monitoring upward and downward flows of replicas,⁴⁸ in order to reduce the human time investment.