Zero-temperature glass transition in two dimensions

Liquids cooled towards the glass transition temperature transform into amorphous solids that have a wide range of applications. While the nature of this transformation is understood rigorously in the mean-field limit of infinite spatial dimensions, the problem remains wide open in physical dimensions. Nontrivial finite-dimensional fluctuations are hard to control analytically, and experiments fail to provide conclusive evidence regarding the nature of the glass transition. Here, we develop Monte Carlo methods for two-dimensional glass-forming liquids that allow us to access equilibrium states at sufficiently low temperatures to directly probe the glass transition in a regime inaccessible to experiments. We find that the liquid state terminates at a thermodynamic glass transition which occurs at zero temperature and is associated with an entropy crisis and a diverging static correlation length. Our results thus demonstrate that a thermodynamic glass transition can occur in finite dimensional glass-formers.

which is denoted as a blue box.

Supplementary Note 1. STRUCTURAL CORRELATIONS
In the main text, we show that ξ PTS increases as temperature decreases. Ref. [1], however, showed that for some computational models made of polydisperse particles, correlation lengths related to the degree of order present increase faster than ξ PTS . In particular, Ref. [1] analyzed the two-points positional and bond-orientational correlations, paying particular attention to the radial decay of the functions g(r) − 1 and g 6 (r)/g(r), respectively Results for these two quantities are reported in Supplementary Figure 2. Here, following Ref. [1], Delaunay neighbors are obtained from a radical Voronoi tessellation. Both functions exhibit clear peaks at distances corresponding to the correlation shells, but their temperature evolution is relatively mild. We fit the peak points with an exponential function of the form C s,6 exp(−r/ξ s, 6 ) in order to extract a correlation function both for positional ξ s and bond- orientational ξ 6 correlations. The temperature evolution of the resulting static correlation lengths is shown in Supplementary Figure 3 together with that of ξ PTS . Over the whole temperature range, we observe an increase by a factor ≈ 2.2 and ≈ 2.7 for ξ s and ξ 6 with saturation at low temperature, which is considerably smaller than the factor of ≈ 5.1 increase observed for ξ PTS . Coupled with the additional verifications for potential crystallization and fractionation, this result rules out the presence of significant structural order in our system, even at extremely low temperatures. Our observations are also remarkably different from those of Ref. [1]; they show that good glass formers are not affected by increases in positional and bond-orientational order.

Supplementary Note 2. CONFIGURATIONAL ENTROPY
The configurational entropy, s conf , is defined as where s tot and s glass are the total entropy and the entropy of a typical glass state, respectively.
We separately measure s tot and s glass by thermodynamic integration based on the scheme developed in Ref. [2]. positional order length, ξ s , bond-orientational order length, ξ 6 , and point-to-set length, ξ PTS (see Supplementary Note 3). Errorbars for ξ s and ξ 6 correspond to the 95% confidence interval of the associated fitting. The increase of the first two lengths is mild compared to that of ξ PTS .

A. Setting
Consider a M -component polydisperse system. (A system with M = N is said to have a continuous polydispersity.) If N m is the number of particles of the m-th species, then the fraction of the m-th species is X m = N m /N , and hence M m=1 N m = N and M m=1 X m = 1. For simplicity, we set all particles masses to unity. We denote particle positions as r N = (r 1 , r 2 , · · · , r N ), and the set of their diameter as Σ N = {σ 1 , σ 2 , · · · , σ N }. In order to consider permutations of particle diameters as additional degrees of freedom, we introduce a permutation π to the set Σ N . A specific sequence of particle diameters is denoted Σ N π , e.g., Σ N π = (σ 3 , σ 8 , σ 5 , · · · ). A total of N ! possible such permutations exists, and for a system with continuous polydispersity, all such permutations are distinguishable.
The system potential energy, U , depends both on particle positions r N and on the permutation π, and is thus formally denoted U (Σ N π , r N ). For notational simplicity, we write U (r N ) = U (Σ N π * , r N ) for the reference system with Σ N π * . The resulting canonical partition function at inverse temperature β = 1/T is where Λ = 2πβ 2 is the thermal de Broglie wavelength with the unit mass. Without loss of generality, we set the Planck constant = 1. Note that Eq. (2) should be distinguished from the conventional partition function, Z, in which only particle positions r N are degrees of freedom, The following subsections describe how Eq. (2) can be used to compute both the total and the glass entropies.

B. Total entropy
The partition function Z in Eq. (2) for the target system βU (Σ N π , r N ) reduces to the conventional partition function Z without permutations in Eq. (3), because diameter permutations are always compensated by position permutations in absence of constraint, i.e., The total entropy computation is therefore equivalent to what has been observed in previous studies [3,4].
Using a high-temperature β → 0 ideal gas as an exactly solvable reference system, we perform a thermodynamic integration over (inverse) temperature up to the target temperature β, where s N m ! is the ideal gas mixing entropy per particle and e pot (β) is the average potential energy per particle. The integration in Eq. (5) requires special care, because e pot (β) diverges in the high-temperature limit [3,4]. We sidestep the difficulty by introducing an intermediate temperature β 0 that separates the very high temperature regime, β ∈ [0, β 0 ], from the rest, β ∈ (β 0 , β]. We thus write where I N is obtained by usual thermodynamic integration, and I F is obtained by fitting the e pot (β) to a polynomial, and then analytically integrating the resulting function [3,4].
The specific polynomial form we use for the high-temperature expansion of a system of soft where the constants A, B, C, and D are determined by fitting, as in Supplementary Figure 4(a). Using Eqs. (6) and (7), we then get which only depends on the fit parameters, A, B, C, and D. Supplementary Figure 4

C. Glass entropy
We evaluate the entropy of glass states by Frenkel-Ladd (FL) thermodynamic integration [5][6][7][8], which requires imposing a harmonic potential with spring constant α on particle positions. The process then entails integrating the long-time limit of the mean-squared displacement starting from a strong α max , at which the system behaves as an Einstein solid, and reaching a weak α min , at which the system is self caged. More specifically, we set where r N 0 is the template configuration from the equilibrium configuration of the target system.
As for the total entropy, we start from the partition function in Eq. (2) for the glass state, Note that the numerator of Eq. (10) is now multiplied by N !, because a given template configuration, r N 0 , selects a single glass basin from the position phase space, while there exists N ! identical such choices, generated by permuting r N 0 . Note also that the presence of the template configuration r N 0 prevents diameter permutations from being compensated by position permutation. The identity in Eq. (4) therefore does not hold in the glass state. The integration limit, lim α min →0 , also requires special conceptual and practical considerations.
Although for FL integration of a crystal α min is chosen to be infinitesimally small, here an additional constraint is that the system should remain within a glass basin and should thus not melt. The practical implementation of this constraint is detailed below.
We compute the entropy s α = βe tot,α − βf α , where e tot,α is the total energy and f α = −(βN ) −1 ln Z α is the free energy. The glass entropy of the target system is then where · · · here denotes averaging over template configurations r N 0 .
For convenience, we also define the following statistical averages, where the superscripts denote statistical averages over positions (T) and permutations (S), evaluated by Monte Carlo (MC) simulations with standard translations and diameter swaps, respectively. Note that any diameter permutation can be expressed as the product of the swaps of two diameters, hence permutation-phase space is properly sampled by swap MC simulations.
Following the conventional Frenkel-Ladd prescription [5] for Eq. (10), we obtain where ∆ T,S α are constrained mean-squared displacements and s mix (r N 0 , β) is a mixing entropy contribution defined by This generalization of the standard FL integration method to systems with continuous polydispersity includes two novel physical features. First, the mean-squared displacement ∆ T,S α has to be evaluated by MC simulations of both translational and swap displacements, and is thus generally distinct from the standard mean-squared displacement, ∆ T α . Because ∆ T,S α accounts for the non-vibrational contributions due to diameter permutations as well as for the purely vibrational contribution, ∆ T,S α ≥ ∆ T α . Including the non-vibrational contribution also markedly improves the estimation of the glass entropy [2,9], as we will see below. Second, the expression contains terms related to the mixing entropy, s Eq. (5). The remaining mixing entropy contribution, s mix in s conf , is finite even for systems with continuous polydispersity. Therefore, with this scheme continuous polydispersity does not present any conceptual or technical difficulty [2].
The key remaining tasks in order to compute s glass involve measuring the mixing entropy contribution s mix and integrating ∆ T,S α . Both are detailed below.

Mixing entropy s mix
The mixing entropy contribution, s mix , is determined by thermodynamic integration, where ∆U mix is a potential energy difference defined by In practice, to get ∆U mix (r N 0 , β ) the system is gradually heated from the target temperature β to an infinite temperature β → 0 using MC simulations with a fraction p swap = 1 of the diameter swaps. Particles are thus kept at the same position as in the template configuration with α. At large α, the system is an Einstein solid with ∆ T,S α = 1/α, but upon decreasing α, ∆ T,S α plateaus. The system is then self caged. Further decreasing α, however, makes the harmonic constraint too weak to prevent the glass state from melting, thus implicitly defining α min . The ensuing particle diffusion explains the upturn of ∆ T,S α . In order to perform the integration in Eq. (15), a practical manipulation of the limit must be used for α < α min . We consider While α max should straightforwardly be chosen in the Einstein solid regime, e.g., we use α max 1 × 10 7 , the choice of α min is not unambiguous. Based on the above discussion, we understand that α min should be within the plateau regime of ∆ T,S α , where ∆ T,S   (shaded region).

D. Potential energy landscape approach
We also consider an alternate approach for estimating s conf based on the potential energy landscape (PEL) [12]. In this approach, the glass entropy is obtained from information about the inherent structures (IS) of the glass state. In order to evaluate the impact of polydispersity on s conf , we employ an effective M * -component approximation as in Ref. [13]. This approach provides an effective mixing entropy s * mix = s (M * ) mix . (The numerical determination of M * is explained below.) We then compute the glass entropy s glass by s glass = s harm + s anh , where s harm and s anh are the harmonic vibrational entropy and its anharmonic correction, respectively [12]. The harmonic term is computed as where · · · IS is an average over IS configurations obtained by the conjugate gradient method and ω a = λ a /m is the square root of eigenvalue λ a of the Hessian of this IS. Supplementary where we used the fact that the system is perfectly harmonic at low T , i.e., s anh (T = 0) = 0.
A low-temperature expansion, e anh (T ) = k=2 c k T k , has T -independent coefficients, c k . Substituting this expansion into Eq. (22) gives The fit of e anh with parameters c 2 and c 3 is shown in Supplementary Figure 8(b), and the resulting s harm + s anh is shown in Supplementary Figure 8(a). The resulting anharmonic contribution is |s anh | < 0.1 in the temperature range of interest.

MW fluctuations effects
The glass entropy measured using the PEL approach also is not affected by MW fluctuations. Consider first the mean-squared displacement of standard solids, |u| 2 . For a monodisperse crystalline solid, one can write where g(ω) and c are the vibrational density of states and the velocity of sound, respectively.
Because one expects a Debye scaling g(ω) ∝ ω d−1 at low ω, in d = 2, |u| 2 ∼ ln L → ∞ diverges in the thermodynamic limit. Writing Eq. (21) using the density of state formalism, by contrast, in d = 2 gives the L-dependent term, ln L L 2 , that vanishes as the system size increases. Additionally, s anh does not depend on system size because e and e IS (and thus e anh ) display no system-size dependence at large enough L. The glass entropy, s harm + s anh , is therefore system-size independent and hence unaffected by MW fluctuations.

Determination of M *
The effective component, M * , is determined based on the potential energy landscape.
As explained in Ref. [13], M * should be such that (i) particle diameter swaps within a single effective species leave the potential energy basin unaffected, and (ii) particle diameter swaps between different species drive the system out of the original basin. To determine M * in practice, we prepare equilibrium configurations of the original continuously polydisperse system characterized by the distribution f (σ). We then decompose f (σ) into M species . This observation indicates that at smaller M , the impact of particle swaps is so strong that the original basin is destroyed, and the system moves to another basin.
From the e IS vs. x = log 10 M plot, the clear crossover between large and small M behaviors determines M * as the intersection of two linear fits as shown in Supplementary Figure 9(a).
We show the resulting s * mix = s (M * ) mix as a function of the temperature in Supplementary  Figure 9(b). We confirm the absence of size dependence by comparing results for systems with N = 1000 and N = 20000.
We also employ an exponential fit as an alternative way to extract M * from the crossover.

A. PTS observables
Similarity between two configurations within the cavity is characterized by the cavity core overlap, q c , computed as in Refs. [14][15][16][17]. (i) We assign a local overlap value to each particle through the overlap estimator function w(z) ≡ exp − z b 2 with b = 0.2; (ii) we perform a linear interpolation through a Delaunay tessellation to define a continuous overlap field; and (iii) we measure the cavity core overlap by taking the average of the field values within the radius r c = 1.0 from the cavity center, evaluated by MC integration with 10 3 points.
For each temperature T and cavity radius R, the PTS correlation function One way to extract the PTS correlation length is through the compressed exponential fit, with the bulk value, Q bulk PTS , evaluated by taking 10 5 pairs of independent configurations in bulk samples. Note that differently from Ref. [14][15][16][17], the compression exponent γ [see Supplementary Figure 11(a)] is here not fixed but treated as an additional fit parameter.
Its value ranges roughly from 2 to 5 from high to low temperatures. Another definition of the PTS length, ξ th PTS , is given by the relation Q fit PTS (ξ th PTS ; T ) − Q bulk PTS ≡ e −1 . A third estimate comes from the peak location of the PTS susceptibility [14] [see Supplementary   Figure 11 Specifically the peak location, ξ peak PTS , is estimated through polynomial extrapolation of five maximal values. All three estimates qualitatively support the conclusion that the PTS

B. PTS equilibration
In order to properly and efficiently sample the cavity configurations, we employ a paralleltempering scheme [18,19] adapted to the cavity sampling as in Refs. [14] with varying tem- As in Ref. [14], we impose the linear relation between replica temperatures and shrinking factors as Ta−T 1 T dec −T 1 = λa−λ 1 λ dec −λ 1 with T dec and λ dec chosen appropriately (see Supplementary  Tables 1-11). The chosen shrinking factors, {λ a } a≥2 , ensure sufficient replica-swap rates. In order to achieve this sampling, replicas are added one by one, with λ 1 = 1 > λ 2 > . . . > λ n , each time targeting a replica-swap acceptance rate of ∼ 20% [17]. This process is stopped upon reaching λ n < λ dec . In Supplementary Tables 1-11, the average number of replicas used, n ave = [n], is recorded for each given temperature and radius.
The quality of the equilibration within each cavity is assessed from monitoring the convergence of two preparation schemes [14,20]: one starting from the original configuration and the other starting from a randomized configuration prepared by running 10 4 MC sweeps with shrunk and heated cavity particles, with (λ, T ) = (0.6, 0.5). Convergence is deemed achieved when obtained through both approaches lie within ±0.1 of each other for each cavity. Here q on c (t) is the cavity core overlap between the original configuration and the equilibrated configuration after t MC sweeps, recorded each t rec = 10 4 MC sweeps. The first s eq configurations are discarded, and thermal averages are taken over the following s prod configurations. With our choice of parallel-tempering parameters (see Supplementary Tables 1-11), for all temperatures and radii, at least 96% of all cavities pass the convergence test. Averaging over cavities results in an even closer agreement between the two schemes, i.e., overlap estimates converge to within ±0.01.
In obtaining PTS correlation functions and PTS susceptibility in Eqs. (26) and (28), respectively, we evaluate core cavity overlaps for s prod pairs of configurations obtained through the two different schemes.

C. Glassiness
As detailed in Ref. [16], PTS observables and equilibration diagnose glassiness by accessing information about the underlying free-energy landscape. On the static side, the probability distribution function of cavity core overlaps exhibits broad fluctuations at the PTS length scale, with bimodal distribution in the deeply glassy regime (see Supplementary   Figure 12). This nontrivial signature of confinement in turn leads to a peak in the PTS susceptibility and to a nonconvex dependence of the PTS correlation, as functions of the cavity radius R (see Supplementary Figure 11) [14]. In nonglassy systems, by contrast, these nontrivial behaviors are absent [16].
Proper sampling within cavity confinement grows increasingly challenging as R decreases.
Without parallel-tempering, the relaxation time explodes for decreasing R (see Supplemen- tary Figure 13). This dynamical observation also bears out that the slowdown in our polydisperse soft-disk system is triggered by the rugged free-energy landscape characteristic of glassiness.  Supplementary Table 1 Cavity PTS measurement parameters T = 0.400, with 100 cavities.
(b) Without parallel tempering. The equilibration time rapidly grows as the cavity radius shrinks, which is interpreted as a finite-size echo of a glass transition [14]. For instance, for R ≤ 3.0 equilibration is not attained even after 10 9 MC steps.

E. Results for hard disks
For the point-to-set length measurement, we also study a two-dimensional hard-disk model, for which the pair interaction is zero for non-overlapping particles and infinite otherwise. The system has the same size distribution f (σ) and size polydispersity δ as the soft-disks described in the main text. Given these parameters, the system is then uniquely characterized by its area fraction ϕ = πN σ 2 /(4V ), 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 to unity for the hard-disks. The pressure P is calculated from the contact value of the pair correlation function properly scaled for a polydisperse system [21]. We use N = 1000 for this model.  Results for hard disks are presented in Supplementary Figure 15. Most technical details are the same as for the soft-disk case. The most notable difference concerns the paralleltempering algorithm, which is adapted from that for d = 3 hard spheres [17], treating the two replicas a = 1 and 2 differently from the rest. Randomized configurations are here prepared by 10 6 MC sweeps with shrunk particles at λ = 0.5, and λ dec s are chosen appropriately (see         It is now clear that, in contrast to d = 3, dynamics in d = 2 is indeed influenced by the presence of the long-wavelength density fluctuations that enhance the mean-squared displacement of particles and thus seemingly breaks the standard cage picture of glassiness.
It has further been established, however, that such dynamical differences can be eliminated by studying bond-orientational relaxation or by introducing the cage-relative mean-squared displacement, which disentangles the Mermin-Wagner fluctuations from the underlying development of glassiness. Upon such disentanglement, the cage picture can be recovered in d = 2 as well.
Here, we disentangle Mermin-Wagner fluctuations from the measurements of the configurational entropy using approaches in the same spirit as those used in previous dynamical studies. Even after disentangling effects of these fluctuations, our determination of the configurational entropy and its comparison with d = 3 results suggest that the glass transition in d = 2 and 3 are fundamentally different. In particular, the latter occurs at a finite