Efficient measurement of point-to-set correlations and overlap fluctuations in glass-forming liquids
EEfficient measurement of point-to-set correlations and overlap fluctuations inglass-forming liquids
Ludovic Berthier, Patrick Charbonneau,
2, 3 and Sho Yaida a) Laboratoire Charles Coulomb, UMR 5221 CNRS and Universit´e de Montpellier, Montpellier,France Department of Chemistry, Duke University, Durham, North Carolina 27708,USA Department of Physics, Duke University, Durham, North Carolina 27708, USA
Cavity point-to-set correlations are real-space tools to detect the roughening of the free-energy landscape thataccompanies the dynamical slowdown of glass-forming liquids. Measuring these correlations in model glassformers remains, however, a major computational challenge. Here, we develop a general parallel-temperingmethod that provides orders-of-magnitude improvement for sampling and equilibrating configurations withincavities. We apply this improved scheme to the canonical Kob-Andersen binary Lennard-Jones model fortemperatures down to the mode-coupling theory crossover. Most significant improvements are noted for smallcavities, which have thus far been the most difficult to study. This methodological advance also enables us tostudy a broader range of physical observables associated with thermodynamic fluctuations. We measure theprobability distribution of overlap fluctuations in cavities, which displays a non-trivial temperature evolution.The corresponding overlap susceptibility is found to provide a robust quantitative estimate of the point-to-setlength scale requiring no fitting. By resolving spatial fluctuations of the overlap in the cavity, we also obtainquantitative information about the geometry of overlap fluctuations. We can thus examine in detail how thepenetration length as well as its fluctuations evolve with temperature and cavity size.PACS numbers: 64.70.Q-, 05.10.-a, 05.20.Jj
I. INTRODUCTION
A well-known difficulty in understanding the increas-ing sluggishness of glass-forming liquids upon loweringtemperature T is that no obvious change to the staticstructure of the liquid accompanies it . In stark contrastto the critical slowing down observed near a standardcritical point, static correlation functions of glass form-ers barely budge while the structural relaxation timescalegrows by more than 15 orders of magnitude. The glassyslowdown has thus instead been attributed to a differenttype of criticality, one at which the free-energy landscapebecomes very rugged, leading to the emergence of manymetastable states separated by growing free-energy bar-riers .In order to capture this ruggedness, the concept of apoint-to-set (PTS) correlation was introduced about adecade ago . PTS correlations generalize multi-pointcorrelations to their infinite-point limit, which makesthem sensitive to the structure of the full configuration.Various PTS geometries have since been considered, in-cluding cavity , random pinning , and walls .Despite subtle differences between them , the overar-ching observation remains. For glass formers the PTScorrelation length grows more rapidly than lengths ex-tracted from two-point correlation functions and is thus an important aspect of glassy physics. How-ever, extracting PTS correlations from numerical simu-lations has thus far been limited by the computational a) Electronic mail: [email protected] challenge of measuring them .In order to better understand the nature of this chal-lenge, let us consider the case of PTS correlation withinthe cavity geometry. First consider a liquid at equilib-rium, and specify a finite cavity, say a sphere of radius R , and pin everything outside of it – thus fixing a set ofparticles. Then allow the particles inside the cavity to ex-plore phase space under the constraints exerted by the setof pinned particles, now acting as an effective quencheddisorder. Particles inside the cavity thus explore the localfree-energy landscape. A standard way of characterizingthese local landscapes, adapted from the study of spinglasses, is to consider the statistics of the overlap, q ( r ),between two equilibrium configurations at a point r insidethe cavity. The materials we consider are characterizedby a very slow dynamics in the bulk, and confining themwith amorphous boundaries is found to make the naturaldynamics of these systems even slower . The majorcomputational problem is thus that of properly samplingthe equilibrium fluctuations of the confined fluid.The origin of the dynamical slowdown may be at-tributed to the confining amorphous boundaries provok-ing an entropy crisis qualitatively similar to the one whichcould be happening at the Kauzmann temperature, T K ,in the bulk liquid . Increasing the confinement wouldreduce the number of accessible states from being expo-nentially large in the number of particles for large cav-ities, to being sub-exponential for small cavities . Thisshift would then be evidenced in the evolution of overlapstatistics in the local landscapes, which is exactly whatPTS correlations purport to capture, and its location de-fines the PTS correlation length, ξ PTS ( T ). From thisperspective, the computational difficulty encountered in a r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n previous studies for R < ∼ ξ PTS should not come as a sur-prise. Confinement should be inducing the analog of anequilibrium (or ‘ideal’) glass transition in a system con-taining a finite number of particles . Proper samplingthen requires very large computer resources, as recentstudies of constrained systems have evidenced . Asa result the entropy crisis picture has only been partiallyvalidated. In particular, although solid evidence for agrowing PTS correlation length has been obtained in var-ious systems , several questions remain unexplored re-garding the nature of the crossover between high and lowoverlap regimes, the associated fluctuations, and the con-nection to dynamical relaxation in the bulk.In this paper we develop a generic computationalmethod to more quickly and reliably sample PTS cor-relations. As a proof of principle, we study the Kob-Andersen binary Lennard-Jones model , which is aclassical glass-forming liquid. Yet our approach is suf-ficiently generic to be applied to other glassy systems,including hard spheres. We use the efficiency of our ap-proach to record novel physical observables, beyond thePTS correlation function, that rely on efficiently sam-pling thermodynamic fluctuations inside a cavity. In par-ticular, we follow the evolution with cavity size and tem-perature of the complete probability distribution func-tion of the overlap. Its variance defines an overlap suscep-tibility that allows us to directly locate the PTS crossoverscale without any empirical fitting. In addition, we areable to spatially resolve overlap fluctuations inside thecavity, giving us access both to their radial profile andtheir orthoradial fluctuations, and thus to make connec-tion with PTS measurements in other geometries, suchas flat walls.The plan of this paper is as follows. In Sections IIand III we detail the simulation model and the parallel-tempering methodology, respectively. In Section IV wepresent the approach for computing PTS correlations aswell as the results, and in Sec. V we evaluate resultsfor the wall and wandering lengths from the cavity PTSsetup. We conclude in Section VI. II. SIMULATION MODEL
We study the behavior of a glass-forming liquid firstproposed by Kob and Andersen . The Kob-Andersenbinary Lennard-Jones (KABLJ) model contains two par-ticle species, denoted A and B , with equal mass m , andinteracting via the pair potential V αβ ( r ) = 4 ε αβ (cid:20)(cid:16) σ αβ r (cid:17) − (cid:16) σ αβ r (cid:17) (cid:21) , (1)where, α, β ∈ { A, B } with parameters ε AB /ε AA = 1 . ε BB /ε AA = 0 . σ AB /σ AA = 0 .
8, and σ BB /σ AA = 0 . r cut αβ = 2 . σ αβ andshifted, such that it vanishes at the cutoff. The relativenumber of particles is N A : N B = 4 : 1 and the overall number density is ρ = N/V = 1 . σ − AA . Note that we re-port below lengths and temperatures in standard dimen-sionless Lennard-Jones units set by σ AA and ε AA /k B .We use bulk samples with N = 135 ,
000 particles ina periodic cubic box, generated as described in Ref. 28,at temperatures T = 1 . , . , . , . , .
45. At eachtemperature, we take 50 equilibrated snapshots, sepa-rated by more than 2 τ α , where τ α is the bulk structuralrelaxation time. Then, within each snapshot, we ran-domly select a position as cavity center. For each center,we fix the particles outside the cavities of radius R , andcontinue sampling configurations for the N cav particlesinside the cavity, as described in Sect. III.Note that in order to directly evaluate how the behav-ior of the cavity changes as function of R , we preserveeach of the cavity centers as R increases. Note also thatthe system size used here is far larger than necessary, be-cause the length scales we extract in Sections IV and Vare rather small. Sampling 50 cavity centers within asingle snapshot would have sufficed to ensure statisticalindependence of the resulting cavities, and would havethus provided similar results. III. PARALLEL-TEMPERING METHODOLOGY ANDCONVERGENCE
A key hurdle to measuring PTS observables in simula-tions has been the large computational effort needed tosample different equilibrium configurations inside a cav-ity. Here, we design a Monte Carlo parallel-tempering scheme (a type of Hamiltonian replica exchange ) thatsidesteps this difficulty by dramatically lowering barrierswithin the rugged free-energy landscape. The approachgoes as follows. We prepare a = 1 , ..., n replicas of agiven configuration inside the cavity. The replicas evolveat different temperatures T a and ‘shrinkage’ parameters λ a ≤ V αβ ( r ; ˜ λ a ) = 4 ε αβ (cid:32) ˜ λ a σ αβ r (cid:33) − (cid:32) ˜ λ a σ αβ r (cid:33) , (2)where ˜ λ a = λ a for a pair of mobile particles within a cav-ity and ˜ λ a = λ a for a pair containing one mobile (insidethe cavity) and one pinned particle (outside the cavity).We draw cavity configurations from the bottom replicawith T = T and λ = 1, whereas the other replicasevolve at higher temperatures T a > T and with smallerparticles, λ 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) displacingparticle i by ∆ x = l ˆ n , where l is uniformly drawn from[0 , .
3] and ˆ n uniformly drawn on the sphere S ; and(iii) accepting/rejecting the displacement according tothe Metropolis criterion. Note that we put a hard spher-ical wall at the edge of the cavity, so that all moves that (a) t ¯ q R = 1 . R = 1 . R = 2 . R = 2 . R = 2 . R = 2 . R = 3 . R = 3 . (b) t ¯ q R = 1 . R = 1 . R = 2 . R = 2 . R = 2 . R = 2 . R = 3 . R = 3 . FIG. 1. (a) Running average of core overlaps, ¯ q , after t MC sweeps from both the original (solid lines) and a randomized(dashed lines) configurations at T = 0 .
51 for a cavity of radius R . (b) The same quantities averaged over 50 cavities. Note thatconvergence is systematically faster in small cavities, where parallel-tempering is more efficient. take a mobile particle outside the cavity radius are re-jected. 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. Thisreplica swapping is again accepted or rejected accordingto the Metropolis criterion, so that our MC algorithm en-sures proper equilibrium sampling. Compared to earlierschemes for cavity studies, the proposed method is dis-tinct from the simpler annealing procedure used beforefor the same model , and is more generally applicablethan the local particle swap Monte Carlo moves that areonly efficient for specific glass-forming models .The quality of the equilibration within each cavity isevaluated by employing two schemes, following the ap-proach developed by Cavagna et al. : we start the sys-tem from (i) the original configuration and (ii) a ran-domized configuration prepared by putting the cavity at T = 1 .
00 and λ = 0 . MC sweeps. We thenrecord the re-equilibrated configurations every t rec = 10 MC sweeps and monitor the core overlap (as defined inSect. IV) between the new configuration and the originalcavity configuration, q onc ( t ), as a function of the numberof MC sweeps, t . Their running average¯ q ( t ) ≡ t/t rec ) ( t/t rec ) (cid:88) s =1 q onc ( t rec s ) (3)decreases for the first scheme and increases for the sec-ond, converging upon equilibration (Fig. 1). The first s eq configurations are discarded, and the overlap for thefollowing s prod , (cid:104) q onc (cid:105) ≡ s prod s eq + s prod (cid:88) s = s eq +1 q onc ( t rec s ) (4) is computed. Convergence is deemed obtained when theresults of both approaches lie within ± q tol of each otherfor 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 Appendix foractual values). However, because the difference betweenthe two approaches is not systematic, averaging over 50cavities results in a very close agreement between the twoschemes, namely in a convergence of ¯ q ( t ) within ± . T ≤ .
51 andsmall radii R ≤ .
0, at least 10 MC sweeps would beneeded to equilibrate the system. Here, by contrast, wecan equilibrate a similar system within ∼ MC sweepsusing ∼
10 replicas. (Note, however, that the computa-tional efficiency of our algorithm depends sensitively onthe details of the parallel-tempering parameters, whichrequire more fine-tuning than the parameters of simplerschemes.) We thus conservatively achieve at least a 100-fold speedup, and that speedup would be even strongerfor state points where simpler MC schemes have not beenfound to converge. By contrast, for R (cid:29) ξ PTS , confine-ment effects are negligible and even simple Monte Carlomoves suffice to sample configurations. The most prob-lematic regime is R ∼ ξ PTS at low temperatures, whereeven relatively large cavities display fine-tuning bottle-necks. That this regime sets the lower limit on the con-vergence criterion is unsurprising, as the optimal spacingbetween adjacent replica parameters generically scales as1 / √ N cav . We thus need ∼ ξ / replicas to properly equi-librate a sample at this most resilient regime.In earlier work , an acceleration of the overlap dy-namics was reported in a binary soft sphere mixture,where traditional Monte Carlo moves are supplementedby binary particle exchanges (or ‘swap’ moves), a methodwhich was used in a series of numerical works . Wecannot directly compare our work to these dynamic mea-surements, as our parallel-tempering scheme does not al-low us to measure time correlation functions at fixed tem-perature, and no evolution of the convergence time wasreported using particle swaps only . We believe, how-ever, that the speedup of the convergence time that we re-port in Fig. 1 is of a different nature, and should uniquelybe attributed to the merits of the parallel-temperingscheme and not to any natural dynamical process tak-ing place in the glass-forming liquid at various degreesof confinement. Methods developed in Ref. 31 could beused in the future to quantify more precisely the conver-gence time of the parallel-tempering scheme using time-correlation functions associated to the random walk ofthe various copies in the ( T a , λ a ) space. IV. POINT-TO-SET OBSERVABLES
Once the glass-forming liquid is confined by amorphousboundaries in a finite cavity we need to analyse the equi-librium fluctuations of the overlap field q ( r ) inside thecavity, which is the principal observable used in this pa-per. To this end, we proceed as follows. Denote a pair ofconfigurations by X = { x i } and Y = { y i } . For each par-ticle x i , find the nearest particle y i nn of the same species,and assign an overlap value q X ; Y ( x i ) ≡ w (cid:0)(cid:12)(cid:12) x i − y i nn (cid:12)(cid:12)(cid:1) ,where w ( z ) ≡ exp (cid:20) − (cid:16) zb (cid:17) (cid:21) , (5)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 passingthrough these points. Specifically, we first perform a De-launay tessellation of space and, to a point r within asimplex spanned by four points { x i } i = i ,i ,i ,i , we asso-ciate a linearly interpolated value q X ; Y ( r ) = (cid:88) i = i ,i ,i ,i c i q X ( x i ) , (6)where { c i } i = i ,i ,i ,i satisfies r = (cid:80) i = i ,i ,i ,i c i x i with (cid:80) i = i ,i ,i ,i c i = 1. Similarly, we obtain q Y ; X ( r ), anddefine q X , Y ( r ) ≡ { q X ; Y ( r ) + q Y ; X ( r ) } , (7) which provides the overlap near the core of the cavity q c ≡ πr (cid:90) | r | 00 (red-cross), 0 . 80 (green-circle), 0 . 60 (cyan-square), 0 . . 45 (black-plus). Solid lines are fits to a compressed exponential, G PTS = A exp[ − ( R/ξ fitPTS ) η ] + G bulkPTS ,where the bulk overlap value, G bulkPTS , is obtained from 4000 pairs of independent configurations in bulk samples, and η = 3is found to work reasonably well for all T . (b) PTS susceptibilities with cavity radius R . Solid lines are guides for theeyes. The dashed line for T = 0 . 45 is extrapolated from R = 3 . R ≈ . 8. See text fordetails. (c) Temperature evolution of inverse ξ fitPTS (red-square), ξ thPTS (green-circle), and ξ peakPTS (blue-diamond), rescaled tounity at T = 0 . 8, along with an extrapolation to prior estimates of bulk T K ≈ . 30 for this model (dashed line) . Thisrepresentation emphasizes the similarity with other constrained phase diagrams, such as random pinning and coupling to areference configuration. (PDF) of core overlaps to display a very broad range offluctuations, while for small cavities configurations areall very similar, hence the overlap is high and does notfluctuate much.These expectations are validated by our numerical re-sults, as illustrated in Fig. 3. For each cavity and tem-perature, we find the PDF to be narrow at large andsmall R , with a maximum in the width of the fluctua-tions for an intermediate cavity size. For some (thoughnot all) cavities, we find that the PTS crossover size cor-responds to a bimodal distribution of overlap values, e.g. Fig. 3b. As temperature is lowered, we find that the frac-tion of bimodal distributions increases, which results ina disorder-average PDF that is itself nearly bimodal, e.g. Fig. 3c. We expect that at even lower temperatures, mostPDF should become bimodal, so that the full disorder-averaged distribution then also becomes bimodal. Bi-modal PDF of the overlap have indeed been found inother constrained systems . We note in pass-ing that a bimodal distribution of overlaps in a cavityis an interesting qualitative signature of the correspon-dence between the PTS crossover length and a crossoveranalogous to an equilibrium glass transition in a finitesize system.The shape of the PDF encodes all the information weneed about the PTS correlation length. In particular, thebreadth of the distribution can be easily detected fromthe PTS susceptibility χ T ( R ) ≡ (cid:104) (cid:104) q (cid:105) J ( R ) − (cid:104) q c (cid:105) J ( R ) (cid:105) . (10)The numerical results in Fig. 2b show that this func-tion has a clear non-monotonic evolution with R , with apeak that gives an unambiguous, fitting-form-free def-inition of ξ PTS . In addition, this approach to deter-mining ξ PTS is computationally advantageous because it does not require a careful fit of the PTS correlationtails, where equilibrium sampling is slower within ourMonte Carlo scheme, and where the overlap value issmall (and hence difficult to precisely measure). Us-ing the susceptibility χ T instead allows us to estimatethe PTS length scale from only a handful of very pre-cise measurements around the susceptibility maximum.Given the limited size of the PTS length for most of thecomputationally-accessible bulk temperatures and the ef-ficiency of parallel-tempering for small cavities, this ap-proach offers an especially effective scheme for determin-ing the PTS length in generic model glass-forming liq-uids.It is comforting to note that the peak location ofthe PTS susceptibility, ξ peakPTS , roughly agrees with thatobtained through the compressed exponential fit, ξ fitPTS .The correspondence is also very good with ξ thPTS , definedas G PTS (cid:0) ξ thPTS (cid:1) ≡ q th with q th = e − . The quantita-tive agreement between these three definitions is demon-strated in Fig. 2c, which shows that all three estimatesgive comparable results, except at the (uninteresting)highest temperature studied, T = 1 . 0, where the verydefinition of a PTS length becomes somewhat problem-atic . The similarity between the different approachesallows us to provide a reasonable estimate of ξ PTS at thelowest temperature considered, T = 0 . 45, which is fairlyclose to estimates of the mode-coupling crossover tem-perature for this model, T MCT ≈ . . At that tem-perature, the susceptibility peak occurs at cavity sizesthat are slightly beyond the computational reach of ouralgorithm. The correlation decay nonetheless provides areasonable estimate of ξ PTS .Although there are systematic quantitative differencesin ξ PTS obtained by the current approach compared withearlier results, these can be explained by the choice of (a) q c P (b) q c P (c) q c P FIG. 3. (a,b) PDF of core overlaps, P ( q c ), at T = 0 . 51 for R = 1 . . . . R ∼ ξ PTS . overlap function and other algorithmic details. Qual-itatively, nothing much differs. The representation ofthe temperature dependence of ξ PTS in Fig. 2c empha-sizes the analogy between cavity PTS measurements andother ways of constraining the available phase space ofthe system that are being actively investigated, such asrandom pinning or coupling to a reference config-uration . In all cases, an external constraint withquenched disorder allows one to induce a sharp transition(pinning, coupling) or crossover (cavity) between an un-constrained phase with low overlap and a localized phasewith large overlap. The data in Fig. 2c thus resemblequalitatively published phase diagrams for the physics ofconstrained glass-forming liquids . Note that, inall of these approaches, one is interested in understand-ing the emergence of metastable states in the physically-relevant temperature regime much above the putativeKauzmann transition, whose existence/absence is there-fore completely immaterial. V. PENETRATION AND WANDERING LENGTHS The availability of well-averaged configurations for cav-ities at R ≈ ξ PTS enables us to have a closer look into thePTS correlation physics. We specifically study two lengthscales associated with the spatial structure of overlap pro-files under pinning. To define them, let us first considerthe PTS physics in the limit R/ξ PTS → ∞ , where con-figurations are frozen beyond a flat wall. The first lengthscale of interest is the wall (or penetration) length, ξ wall ∞ ,which corresponds to the characteristic distance from thewall for the decay of the overlap . The spatial fluc-tuations of this decay length also result in a wandering ofthe iso-overlap surface. The second length scale is thusthe wandering length, ξ ⊥∞ , which characterizes the sizeof these fluctuations in the direction orthogonal to thewall . A priori, there is no relation between the growthof these two scales. The wandering length could evendecrease as temperature decreases. Recent theoretical work , however, suggests that these two length scalesgrow at the same rate, albeit remain subdominant to thecavity PTS correlation length, ξ PTS . Here, we analyzeavatars of these two length scales, ξ wall (cid:63) and ξ ⊥ (cid:63) , in cav-ities of size R ≈ R (cid:63) ≡ ξ PTS . We emphasize that theseavatars are distinct from their counterparts in the flatwall limit, although they are smoothly connected to ξ wall ∞ and ξ ⊥∞ as functions of R/ξ PTS .To extract the wall length, we first consider the ra-dial dependence of the overlap profile within a cavity,[ (cid:104) q ( r ) (cid:105) J ( R ) ], which can be obtained by binning the over-lap profile away from the cavity center (Fig. 4a). Whenthe cavity size R ≈ ξ PTS , nonlocal constraints imposedby disorder is felt throughout the full cavity, via therarefaction of metastable states in the local free-energylandscapes. The near collapse observed in the inset ofFig. 4a suggests that the extent of the pinning effect inthis regime, ξ wall (cid:63) , is proportional to ξ PTS . When con-fined on the scale of ξ PTS the system thus appears fullycorrelated, as if it were at a sharp transition point. Bycontrast, the data in Fig. 4b, obtained for R ≈ ξ PTS ,show that the penetration length measured in a geom-etry that interpolates between the finite cavity and theflat wall seems to have a milder temperature dependencethan ξ PTS . This behavior is as expected, given that ξ wall ∞ is observed to grow subdominantly to ξ PTS .To analyze the avatar of wandering length, we considerthe orthoradial fluctuations of the overlap field inside thecavity. Recent field-theoretic calculations suggest thatthe overlap profile (cid:104) q ( r ) (cid:105) J ( R ) for each individual cav-ity is aspherical , even though it appears sphericallysymmetric when averaged over the disorder of the cavitypinning. A hint of this effect can be seen in Fig. 5 byexamining the iso-overlap surfaces (cid:104) q ( r ) (cid:105) J ( R ) = q iso , (11)which are generically rather bumpy, even after smooth-ing particle-scale fluctuations. In order to quantify as-phericity, we calculate (cid:104) q ( r ) (cid:105) J ( R ) on 10 points uniformlydistributed within a cavity. We then consider a point (a) r [ h q i ] X (b) r [ h q i ] X FIG. 4. Spherical (or s -wave) radial profile of [ (cid:104) q ( r ) (cid:105) J ( R ) ] for(a) R ≈ ξ PTS and (b) R ≈ ξ PTS . Colors as in Fig. 2. Inboth panels, insets represent the overlap profile as a functionof rescaled variable X ≡ ( R − r ) /ξ PTS , i.e. , the distance fromthe cavity edge rescaled by PTS length. Rescaling nearlycollapses the results in (a) but not in (b), which is hinting ata decoupling between ξ wall ∞ and ξ PTS . to be part of the iso-overlap surface at q iso if it fallswithin [ q iso − ∆ q , q iso + ∆ q ] with ∆ q = 0 . 01, and de-note the resulting set of points as { r i } i =1 ,...,n . Thevariance, χ ⊥ ( q iso ) J ( R ) , of distances to the edge of thecavity, (cid:8) R − (cid:12)(cid:12) r i (cid:12)(cid:12)(cid:9) i =1 ,...,n , directly quantifies the sample-to-sample deviations from spherical iso-overlap surfaces.The average of this quantity gives an estimate of thedistance between the iso-overlap surface and the cavitywall, R − r ( q iso ) J ( R ) . After averaging over the pinningdisorder, we obtain χ ⊥ ( q iso ) ≡ (cid:104) χ ⊥ ( q iso ) J ( R ) (cid:105) . (12) (a) (b)(c) FIG. 5. Iso-overlap surfaces (cid:104) q ( r ) (cid:105) J ( R ) = q iso for R = 3 . T = 0 . 51. Roughness propagates from (a) high overlap q iso = 0 . q iso = 0 . q iso = 0 . q is o χ ⊥ X √ χ ⊥ / ξ P T S FIG. 6. Asphericity, χ ⊥ ( q iso ), of the iso-overlap surface atvarious values q iso of the average overlap. Colors as in Fig. 2.For each disorder realization, each increment ∆ q of q iso (seetext for details) contains at least 100 sampling points overthe range of q iso shown above. Inset: Rescaled radial extentof the fluctuating iso-overlap surface (cid:112) χ ⊥ ( q iso ) as a functionof rescaled average distance from the cavity edge (see Fig. 4)for cavities with R ≈ ξ PTS . Within the thermal noise, areasonably good data collapse is observed. 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 roughthat, as it is defined, χ ⊥ does not approximate ( ξ ⊥ ) well anymore. Note indeed that in Fig. 5c the surfaceeven becomes topologically distinct from a sphere. Thecollapse in the inset of Fig. 6 suggests again the extentof orthoradial fluctuations of the overlap in this regime, ξ ⊥ (cid:63) , may be collapsed by the PTS length scale.On the whole, these results suggest that the spatialstructure of overlap profiles in this regime is governed bya single length scale, ξ PTS , although a conclusive inter-pretation of this effect has yet to be reached. It wouldbe interesting to examine how this structure morphs intoin the flat wall geometry limit, where ξ wall ∞ and ξ ⊥∞ areexpected to grow subdominantly to ξ PTS . The largestcavities we equilibrated in this paper are of the order of2 ξ PTS , and our preliminary analysis did not conclusivelydetermine whether they belong to the flat wall regime ornot, in particular for the wandering length. A study ofthe evolution of the overlap profiles for a broader rangeof cavity sizes, including the careful analysis of the wan-dering length in the wall geometry, should be performedin future work. VI. CONCLUSIONS The results of this work suggest that parallel tem-pering is a method of choice for studying PTS correla-tions in generic model glass-forming liquids. This schemehas indeed allowed us to quantify PTS correlations in acanonical glass former over the temperature regime thatis generically studied in the bulk, bypassing the slow-down due to confinement by amorphous boundaries. Thisscheme thus provides an efficient solution to a major com-putational obstacle that has prevented progress in studiesof static correlations in glass formers for nearly a decade.In addition, by defining a novel overlap susceptibilityfor finite cavities, we have obtained PTS correlation in-formation well beyond the regime typically attained insimulations, and from a reduced computational effort.Our results are broadly consistent with a growing staticcavity PTS length scale as the temperature of a glass-forming liquid is lowered, and thus with glassiness be-ing associated with a roughening of the free energy land-scape.As expected from extrapolating earlier PTS results,however, the measured static length scale does not growby much more than a factor of two upon approachingthe mode-coupling crossover temperature. In Fig. 7, wereplot our data for ξ PTS from Fig. 2, and compare itstemperature evolution with several dynamic correlationlengthscales. Given apparent subtleties in extracting dy-namical length scales, we record three independent setsof data taken from Refs. 20, 43, and 44, taken from eithermeasurements of four-point dynamic correlations orfrom dynamic profiles near a wall for the same model.This comparison shows that all three dynamic lengthsgrow more rapidly than the static one toward low tem-peratures, but the difference is relatively modest. Astronger decoupling between static and dynamic quanti-ties was observed in hard-sphere systems than whatwe find in the KABLJ model. The distinction between T ˜ ξ FIG. 7. Temperature evolution of various normalized lengthscales, ˜ ξ ≡ ξ/ξ ( T ), in the KABLJ model, rescaled to unityat T = 0 . ξ as determined from standard fitting of afour-point dynamic correlation function based on single par-ticle displacements, data taken from Ref. 43 (cyan-circle) andRef. 44 (green-cross); ξ dyn (blue-diamond) obtained from fit-ting dynamic profiles near an amorphous wall ; and ξ PTS (red-square) is the PTS length scale measured in this work.A modest decoupling between static and dynamic lengths isobserved. the two models is in harmony with previous findings thatsignatures of the mode-coupling crossover (associated toa rapid growth of a dynamic correaltion length) appearweaker in the KABLJ model than in hard or quasi-hardspheres .While there have been recent attempts at experi-mentally measuring PTS correlations in colloidal glass-formers , corresponding studies for molecular glass-formers, roughly 10 orders of magnitude more sluggish,still remain inaccessible. Analysis of such extremely slug-gish systems would be particularly helpful in properlyassessing the theoretical framework that surrounds PTScorrelations. In future work, we will thus couple thepresent approach with other enhanced sampling tech-niques in an attempt to push back the computationallyaccessible boundary for the study of glassiness. Suchstudies will also be useful for systematically investigat-ing the behavior of various length scales ξ PTS , ξ wall (cid:63), ∞ , and ξ ⊥ (cid:63), ∞ and how they relate to one another deep in the dy-namically sluggish regime. ACKNOWLEDGMENTS We acknowledge many fruitful exchanges and corre-spondence with G. Biroli, C. Cammarota, G. Hocky,T. Kawasaki, V. Martin-Mayor, D. Reichman, G. Sza-mel, and G. Tarjus. We also acknowledge the Duke Com-pute Cluster for computational resources. The researchleading to these results has received funding from theEuropean Research Council under the European UnionSeventh Framework Programme (FP7/20072013)/ERCGrant agreement No. 306845 (LB). PC and SY ac-knowledge support from the National Science FoundationGrant No. NSF DMR-1055586 and the Sloan Foundation. L. Berthier and G. Biroli, Rev. Mod. Phys. , 587 (2011). T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys.Rev. A , 1045 (1989). M. M´ezard and G. Parisi, J. Chem. Phys. , 1076 (1999). V. Lubchenko and P. G. Wolynes, Ann. Rev. Phys. Chem. ,235 (2007). J.-P. Bouchaud and G. Biroli, J. Chem. Phys. , 7347 (2004). A. Montanari and G. Semerjian, J.Stat.Phys. , 23 (2006). A. Cavagna, T. S. Grigera, and P. Verrocchio, Phys. Rev. Lett. , 187801 (2007). G. Biroli, J.-P. Bouchaud, A. Cavagna, T. Grigera, and P. Ver-rocchio, Nat. Phy. , 771 (2008). F. Sausset and D. Levine, Phys. Rev. Lett. , 045501 (2011). G. M. Hocky, T. E. Markland, and D. R. Reichman, Phys. Rev.Lett. , 4626 (2012). G. Biroli, S. Karmakar, and I. Procaccia, Phys. Rev. Lett. ,165701 (2013). B. Charbonneau, P. Charbonneau, and G. Tarjus, Phys. Rev.Lett. , 035701 (2012). C. Cammarota and G. Biroli, Proc. Nat. Acad. Sci. U. S. A. ,8850 (2012). R. L. Jack and L. Berthier, Phys. Rev. E , 021120 (2012). B. Charbonneau, P. Charbonneau, and G. Tarjus, J. Chem.Phys. , 12A515 (2013). P. Charbonneau and G. Tarjus, Phys. Rev. E , 042305 (2013). W. Kob and L. Berthier, Phys. Rev. Lett. , 245702 (2013). M. Ozawa, W. Kob, A. Ikeda, and K. Miyazaki, Proc. Nat. Acad.Sci. U.S.A. , 6914 (2015). W. Kob, S. Rold´an-Vargas, and L. Berthier, Nat. Phys. , 164(2012). G. M. Hocky, L. Berthier, W. Kob, and D. R. Reichman, Phys.Rev. E , 052311 (2014). L. Berthier and W. Kob, Phys. Rev. E , 011102 (2012). A. Cavagna, T. S. Grigera, and P. Verrocchio, J. Chem. Phys. , 204502 (2012). C. Cammarota, G. Gradenigo, and G. Biroli, Phys. Rev. Lett. , 107801 (2013). W. Kauzmann, Chem. Rev. , 219 (1948). L. Berthier and D. Coslovich, Proc. Nat. Acad. Sci., U.S.A. ,11668 (2014). W. Kob and H. C. Andersen, Phys. Rev. Lett. , 1376 (1994). W. Kob and H. C. Andersen, Phys. Rev. E , 4626 (1995). E. Dyer, J. Lee, and S. Yaida, (2013), arXiv:1309.5085 . D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, New York, 2nd Ed., 2001). H. Fukunishi, O. Watanabe, and S. Takada, J. Chem. Phys. ,9058 (2002). L. A. Fernandez, V. Martin-Mayor, S. Perez-Gaviro, A. Taran-con, and A. P. Young, Phys. Rev. B , 024422 (2009). F. Sciortino, W. Kob, and P. Tartaglia, Phys. Rev. Lett. ,3214 (1999). L. Berthier, Phys. Rev. E , 022313 (2013). L. Berthier and R. L. Jack, Phys. Rev. Lett. , 205701 (2015). S. Franz and G. Parisi, Phys. Rev. Lett. , 2486 (1997). C. Cammarota, A. Cavagna, I. Giardina, G. Gradenigo, T. S. Grigera, G. Parisi, and P. Verrocchio, Phys. Rev. Lett. ,055703 (2010). G. Parisi and B. Seoane, Phys. Rev. E , 022309 (2014). R. L. Jack and J. P. Garrahan, (2015), arXiv:1508.06470 . C. Cammarota, G. Biroli, M. Tarzia, and G. Tarjus, Phys. Rev.Lett. , 115705 (2011). G. Gradenigo, R. Trozzo, A. Cavagna, T. S. Grigera, and P. Ver-rocchio, J. Chem. Phys. , 12A509 (2013). G. Biroli and C. Cammarota, (2014), arXiv:1411.4566 . A. Adams, T. Anous, J. Lee, and S. Yaida, Phys. Rev. E ,032148 (2015). E. Flenner, H. Staley, and G. Szamel, Phys. Rev. Lett. ,097801 (2014). H. Shiba, T. Kawasaki, and K. Kim, (2015), arXiv:1510.02546 . L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli,Phys. Rev. E , 031502 (2012). S. Gokhale, K. H. Nagamanasa, R. Ganapathy, and A. K. Sood,Nat. Commun. , 4685 (2014). K. H. Nagamanasa, S. Gokhale, A. K. Sood, and R. Ganapathy,Nat. Phys. , 403 (2015). H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer, J. Stat.Mech. , P03018 (2006). Appendix: Parallel-tempering parameter choices The lists of replica parameters used for each temper-ature and radius explored in this paper, along with s eq , s prod , and equilibration success rate are provided in Ta-bles I-V. In general, we choose { ( T a , λ a ) } a =1 ,...,n , suchthat they satisfy (when n > 1) the linear relation T a − T T dec − T = λ a − λ λ dec − λ , (A.1)where ( T dec , λ dec ) is chosen such that the particles candecorrelate the overlap effectively. We always choose thelast 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 configura-tions, and compare the thermal averages [computed withEq. (4)] obtained from two schemes: one starting fromthe original configuration and the other from a random-ized configuration. The convergence criterion requiresthat results from both approaches lie within ± q tol of eachother. For each data point, we record as success rate howmany of 50 cavities satisfy this criterion for a particular q tol . Globally, for q tol = 0 . 1, 98% of cavities are deemedto have converged, for all temperatures and radii.Note that the parameters chosen are neither uniquenor optimal. Some of the chosen parameters result innear bottlenecks in the exchanges of replicas. Tuningthe replica parameters for each cavity by hand could cer-tainly help sampling low-temperature configurations. Amore promising and general way would be to algorithmi-cally tune the replica parameters by monitoring upwardand downward flows of replicas , in order to reduce thehuman time investment.0 R . . . . . . . s eq s prod T dec λ dec { λ a } T = 1 . 00. Success rate isdetermined for q th = 0 . R . . . . . . . s eq s prod T dec λ dec { λ a } T = 0 . 80. Success rate isdetermined for q th = 0 . R . . . . . . . . s eq s prod T dec λ dec { λ a } T = 0 . 60. Success rate isdetermined for q th = 0 . R . . . . . . . . . s eq s prod T dec λ dec { λ a } T = 0 . 51. Success rate isdetermined for q th = 0 . R . . . . . . . s eq s prod T dec λ dec { λ a } T = 0 . 45. Success rate isdetermined for q th = 0 ..