Growing timescales and lengthscales characterizing vibrations of amorphous solids

Preparation of Glass States

Experimentally, glasses are obtained by a slow thermal or compression annealing, the rate of which determines the location of the glass transition (1, 2). We find that a detailed numerical analysis of the Gardner transition requires the preparation of extremely well-relaxed glasses (corresponding to structural relaxation timescales challenging to simulate) to study vibrational motion inside the glass without interference from particle diffusion. We thus combine a very simple glass-forming model––a polydisperse mixture of hard spheres––to an efficient Monte Carlo scheme to obtain equilibrium configurations at unprecedentedly high densities, i.e., deep in the supercooled regime. The optimized swap Monte Carlo algorithm (12), which combines standard local Monte Carlo moves with attempts at exchanging pairs of particle diameters, indeed enhances thermalization by several orders of magnitude. Configurations contain either N = 1,000 or N = 8,000 (results in Figs. 1–3 are for N = 1,000; results in Fig. 4 are for N = 8,000) hard spheres with equal unit mass m and diameters independently drawn from a probability distribution Pσ(σ)∼σ−3, for σmin≤σ≤σmin/0.45. We similarly study a 2D bidisperse model glass former and report the main results in SI Appendix.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Two glass phases. Inverse reduced pressure–packing fraction (1/p v φ) phase diagram for polydisperse hard spheres. The equilibrium simulation results at φg (green squares) are fitted to the liquid EOS (Eq. 1, green line). The dynamical cross-over, φd, is obtained from the liquid dynamics. Compression annealing from φg up to jamming (blue triangles) follows a glass EOS (fit to Eq. 3, dashed lines). At φG (red circles and line) with a finite p, stable glass states transform into marginally stable glasses. Snapshots illustrate spatial heterogeneity above and below φG, with sphere diameters proportional to the linear cage size and colors encoding the relative cage size, ui (see the text).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Emergence of slow vibrational dynamics. (A) Time evolution of Δ(t,tw) for several tw and φ (from top to bottom, φ=0.645,0.67,0.68,0.684,0.688), following compression from φg=0.643. For φ≳φG=0.684, Δ(t,tw) displays strong aging. (B) Comparison between Δ(t,tw) (points) and ΔAB(t+tw) (lines) for the longest waiting time tw=1,024. For φ<φG both observables converge to the same value within the time window considered, but not for φ>φG. (C) The time evolution of δΔ(t,tw) at tw=2 displays a logarithmic tail, which provides a characteristic relaxation time τ. (D) As φ approaches φG, τ grows rapidly for any tw.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Global fluctuations of the order parameter. (A) Probability distribution functions for ΔAB and Δ above, at, and below the Gardner cross-over, φG=0.670(2) for φg=0.630. Vertical lines mark 〈Δ〉 (solid) and 〈ΔAB〉 (dashed), which also represent the peak positions. (B) Comparing 〈Δ〉 and 〈ΔAB〉 shows that the average values separate for φ≳φG (Data are multiplied by 5k, where k=0,1,…,5 for φ=0.655,0.643,…,0.598, respectively.) Around φG, (C) the global susceptibility χAB grows very rapidly, and (D) the skewness ΓAB peaks. Numerical estimates for φG are indicated by vertical segments.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Growing correlation length. (A) Spatial correlator GL(r) (Eq. 6) for different φ (from bottom to top, φ=0.65,0.67,0.675,0.677,0.68,0.682,0.687) annealed from φg=0.640, with N=8,000 (larger systems are used here to significantly measure the growth of ξ). (B) Fitting GL(r) to Eq. 7 (lines in A) provides the correlation length ξ, which grows with φ and becomes comparable to the linear system size upon approaching φG=0.682(2) (dashed line).

We mimic slow annealing in two steps (Fig. 1). First, we produce equilibrated liquid configurations at various densities φg using our efficient simulation scheme, concurrently obtaining the liquid equation of state (EOS). The liquid EOS for the reduced pressure p=βP/ρ, where ρ is the number density, β is the inverse temperature, and P is the system pressure, is described bypliquid(φ)=1+f(φ)[pCS(φ)−1],[1]with pCS(φ) from ref. 13pCS(φ)=11−φ+3s1s2s3φ(1−φ)2+s23s32(3−φ)φ2(1−φ)3,[2]where sk is the kth moment of Pσ(σ), and f(φ)=0.005−tanh[14(φ−0.79)] are fitted quantities. The structure of the equilibrium configurations generated by the swap algorithm has been carefully analyzed. Unlike for other glass formers (14, 15), no signs of orientational or crystalline order were observed (16, 17). Following the strategy of ref. 18, we also obtain the mode-coupling theory dynamical cross-over φd=0.594(1) (SI Appendix). We have not analyzed the compression of equilibrium configurations with φg<φd, as done in earlier studies (19, 20), because structural relaxation is not well decoupled from vibrational dynamics, although the obtained jammed states should have equivalent properties.

Second, we use these liquid configurations as starting points for standard molecular dynamics simulations during which the system is compressed out of equilibrium up to various target φ>φg (21). Annealing is achieved by growing spheres following the Lubachevsky–Stillinger algorithm (21) at a constant growth rate γg=10−3 (see SI Appendix for a discussion on the γg -dependence). The average particle diameter, σ¯, serves as unit length, and the simulation time is expressed in units of βmσ¯2. To obtain thermal and disorder averaging, this procedure is repeated over Ns samples (Ns≈150 for N=1,000 and Ns=50 for N=8,000), each with different initial equilibrium configurations at φg, and over Nth=64−19,440 independent thermal (quench) histories for each sample. Quantities reported here are averaged over Ns×Nth quench histories, unless otherwise specified. The nonequilibrium glass EOSs associated with this compression (dashed lines) terminate (at infinite pressure) at inherent structures that correspond, for hard spheres, to jammed configurations (blue triangles). To capture the glass EOSs, we use a free-volume scaling around the corresponding jamming point φJpglass(φ;φg)=CφJ(φg)−φ,[3]where the constant C weakly depends on φg.

Our numerical protocol is analogous to varying the cooling rate––and thus the glass transition temperature––of thermal glasses, and then further annealing the resulting amorphous solid. Each value of φg indeed selects a different glass, ranging from the onset of sluggish liquid dynamics around the dynamical cross-over (1, 2), φd, to the very dense liquid regime where diffusion and vibrations (β-relaxation processes) are fully separated (2). For sufficiently large φg, we thus obtain unimpeded access to the only remaining glass dynamics, i.e., β-relaxation processes (4).

Growing Timescales

A central observable to characterize glass dynamics is the mean-squared displacement (MSD) of particles from position ri(tw)Δ(t,tw)=1N∑i=1N〈|ri(t+tw)−ri(tw)|2〉,[4]averaged over both thermal fluctuations and disorder, where time t starts after waiting time tw when compression has reached the target φ. The MSD plateau height at long times quantifies the average cage size (SI Appendix). Because some of the smaller particles manage to leave their cages, the sum in Eq. 4 is here restricted to the larger half of the particle size distribution (SI Appendix). When φ is not too large, φ≳φg, the plateau emerges quickly, as suggested by the traditional view of caging in glasses (Fig. 2A). When the glass is compressed beyond a certain φG, however, Δ(t,tw) displays both a strong dependence on the waiting time tw, i.e., aging, and a slow dynamics, as captured by the emergence of two plateaus. These effects suggest a complex vibrational dynamics. Aging, in particular, provides a striking signature of a growing timescale associated with vibrations, revealing the existence of a “glass transition” deep within the glass phase.

To determine the timescale associated with this slowdown, we estimate the distance between independent pairs of configurations by first compressing two independent copies, A and B, from the same initial state at φg to the target φ, and then measuring their relative distanceΔAB(t)=1N∑i=1N〈|riA(t)−riB(t)|2〉,[5]so that ΔAB(t→∞)≃Δ(t→∞,tw→∞), as shown in Fig. 2B. The two copies share the same positions of particles at φg, but are assigned different initial velocities drawn from the Maxwell–Boltzmann distribution. The time evolution of the difference δΔ(t,tw)=ΔAB(tw+t)−Δ(t,tw) indicates that whereas the amplitude of particle motion naturally becomes smaller as φ increases, the corresponding dynamics becomes slower (Fig. 2C). In other words, as φ grows particles take longer to explore a smaller region of space. In a crystal, by contrast, δΔ(t,tw) decays faster under similar circumstances. A relaxation timescale, τ, can be extracted from the decay of δΔ(t,tw) at large t, whose logarithmic form, δΔ(t,tw)∼1−ln⁡t/ln⁡τ, is characteristic of the glassiness of vibrations. As φ→φG, we find that τ dramatically increases (Fig. 2D), which provides direct evidence of a marked cross-over characterizing the evolution of the glass upon compression.

Global Fluctuations of the Order Parameter

This sharp dynamical cross-over corresponds to a loss of ergodicity inside the glass, i.e., time and ensemble averages yield different results. To better characterize this cross-over, we define a timescale τcage for the onset of caging [τcage≈O(1); SI Appendix], and the corresponding order parameters ΔAB≡ΔAB(τcage) and Δ≡Δ(τcage,tw=0).

The evolution of the probability distribution functions, P(ΔAB) and P(Δ), as well as their first moments, 〈ΔAB〉 and 〈Δ〉, are presented in Fig. 3 A and B for a range of densities across φG. For φ<φG, dynamics is fast, 〈ΔAB〉 and 〈Δ〉 coincide, and P(ΔAB) and P(Δ) are narrow and Gaussian-like. For φ>φG, however, the MSD does not converge to its long-time limit, 〈Δ〉<〈ΔAB〉, which indicates that configuration space explored by vibrational motion is now broken into mutually inaccessible regions. Interestingly, the slight increase of 〈ΔAB〉 with φ in this regime (Fig. 3B) suggests that states are then pushed further apart in phase space, which is consistent with theoretical predictions (11). When compressing a system across φG, its dynamics explores only a restricted part of phase space. As a result, ΔAB displays pronounced, non-Gaussian fluctuations (Fig. 3A). Repeated compressions from a same initial state at φg may end up in distinct states, which explains why ΔAB is typically much larger and more broadly fluctuating than Δ (Fig. 3A). These results are essentially consistent with theoretical predictions (9, 11), which suggest that for φ>φG, P(ΔAB) should separate into two peaks connected by a wide continuous band with the left-hand peak continuing the single peak of P(Δ). The very broad distribution of ΔAB further suggests that spatial correlations develop as φ→φG, yielding strongly correlated states at larger densities.

To quantify these fluctuations we measure the variance χAB and skewness ΓAB (SI Appendix and ref. 11) of P(ΔAB) (Fig. 3 C and D). The global susceptibility χAB is very small for φ<φG and grows rapidly as φG is approached, increasing by about two decades for the largest φg considered (Fig. 3C). Whereas χAB quantifies the increasing width of the distributions, ΓAB reveals a change in their shapes. For each φg we find that ΓAB is small on both sides of φG with a pronounced maximum at φ=φG (Fig. 3D). This reflects the roughly symmetric shape of P(ΔAB) around 〈ΔAB〉 on both sides of φG and the development of an asymmetric tail for large ΔAB around the cross-over, a known signature of sample-to-sample fluctuations in spin glasses (22) and mean-field glass models (11). Note that because the skewness maximum gives the clearest numerical estimate of φG, we use it to determine the values reported in Fig. 1.

Growing Correlation Length

The rapid growth of χAB in the vicinity of φG suggests the concomitant growth of a spatial correlation length, ξ. Its measurement requires spatial resolution of the fluctuations of ΔAB, hence for each particle i we define ui=(|riA−riB|2/〈ΔAB〉)−1 to capture its contribution to deviations around the average 〈ΔAB〉. A first glimpse of these spatial fluctuations is offered by snapshots of the ui field (Fig. 1), which appear featureless for φ<φG, but highly structured and spatially correlated for φ≳φG. More quantitatively, we define the spatial correlatorGL(r)=〈∑μ=13∑i≠juiujδ(r−|ri,μA−rj,μA|)〉〈∑μ=13∑i≠jδ(r−|ri,μA−rj,μA|)〉,[6]where ri,μ is the projection of the particle position along direction μ. Even for the larger system size considered, measuring GL(r) is challenging because spatial correlations quickly become long-ranged as φ→φG (Fig. 4A). Fitting the results to an empirical form that takes into account the periodic boundary conditions in a system of linear size L,GL(r)∼1rae−(r/ξ)b+1(L−r)ae−[(L−r)/ξ]b,[7]where a and b are fitting parameters, nonetheless confirms that ξ grows rapidly with φ and becomes of the order of the simulation box at φ>φG (Fig. 4B). Note that although probed using a dynamical observable, the spatial correlations captured by GL(r) are conceptually distinct from the dynamical heterogeneity observed in supercooled liquids (23), which is transient and disappears once the diffusive regime is reached.