Structural order as a genuine control parameter of dynamics in simple glass formers

Glass transition is characterised by drastic dynamical slowing down upon cooling, accompanied by growing spatial heterogeneity. Its rationalisation by subtle changes in the liquid structure has been long debated but remains elusive, due to intrinsic difficulty in detecting the underlying complex structural ordering. Here we report that structural order parameter characterising local packing capability can well describe the glassy dynamics not only macroscopically but also microscopically, no matter whether it is driven by temperature or density. A Vogel-Fulcher-Tammann (VFT)-like relation is universally identified between the structural relaxation time and the order parameter for supercooled liquids with isotropic interactions. More importantly, we find such an intriguing VFT-like relation to be statistically valid even at a particle level, between spatially coarse-grained structural order and microscopic particle-level dynamics. Such a unified description of glassy dynamics based solely on structural order is expected to contribute to the ultimate understanding of the long-standing glass-transition problem.


SUPPLEMENTARY NOTE 1. DETAILS OF SIMULATION METHODS
We have performed molecular dynamics simulations in square boxes in 2D and cubic boxes in 3D, with periodic boundary conditions. Four types of interaction potentials are considered, including harmonic, Lennard-Jones (LJ), Weeks-Chandler-Andersen (WCA), and purely hard interactions. For the former three, the simulations are carried out using velocity Verlet algorithms in the N V T ensemble with a Berendsen thermostat [1,2], whereas the event-driven molecular dynamics is used for the fourth [3]. We study both binary and continuously polydisperse systems for the harmonic systems and only polydisperse systems for the rest. In total, sixteen different glass formers are studied, among which, for the nine listed in Table 1 of the main text, we characterize both macroscopic and microscopic properties, and for the other seven listed in Supplementary Table 1 we characterize only the macroscopic properties. Therefore, we have accessed the degrees of freedom in terms of spatial dimensions, interactions, compositions, and also regimes in the phase space controlled by temperature or density, and confirmed consistent results. Although we have not studied all combinations of these degrees of freedom exhaustively, generality is expected for what we call "hard-sphere-like" systems. More discussion is given accordingly in the text.
Three sets of simulations are carried out to obtain results shown in the paper: (1) continuous cooling from a hightemperature liquid state to glass state, to mimics the experimental glass transition process; (2) long-time simulations at constant temperatures, to access structural and dynamic properties in thermal equilibrium; (3) isoconfigurational simulations [4,5], to access microscopic dynamical information. In the following, we first give the simulation details for the harmonic systems. Then, we describe the systems with LJ, WCA, and purely hard interactions, with slightly different settings in simulation details.
Harmonic systems. For all cases, the systems are first equilibrated at high temperature (where τ α < 10) for t = 1000. For the continuous cooling simulations, the systems are then cooled down to T = 10 −4 in 2D and T = 2 × 10 −5 in 3D at a constant cooling rate γ = dT /dt. For the equilibrium simulations with which we calculate the self-intermediate scattering functions on the fly, the systems are cooled down to target temperatures within a time scale of several tens of the structural relaxation time τ α and then equilibrated for at least several tens of τ α before the production sampling. For the isoconfigurational ensemble [4,5] with which we calculate the microscopic structure relaxation time, the systems are cooled down to target temperatures within around 30τ α and then equilibrated for another 30τ α . After the equilibration run, we quench the system to the nearest inherent structure by using the FIRE algorithm [6] to remove the thermally-excited random displacements. Using this inherent structure as the initial configuration, 100 trajectories are simulated with different momenta assigned randomly from the appropriate Maxwell-Boltzmann distribution. Results of selected state points using 200 trajectories confirm good convergence in a statistical sense. For the 2D polydisperse system with ∆ = 11% and 3D polydisperse system with ∆ = 8%, 20 initial configurations are used in the calculations to ensure good statistics; for the rest, 30 initial configurations are used. Because of the limitation of the computational power, for the state point of the lowest temperatures under study, 10 initial configurations are used. We set the target temperatures covering a wide range of τ α (approximately, τ α ∈ [30, 10 5 ] in 2D and [100, 10 5 ] in 3D).
Lennard-Jones and Weeks-Chandler-Andersen systems. We have investigated the Weeks-Chandler-Andersen (WCA) and Lennard-Jones (LJ) systems and confirmed the generality of our findings. The interaction potential between particles i and j is when r ij /σ ij < R c and zero otherwise, where r ij is the particle separation, σ ij is the sum of the particle radii, and f (r ij ) guarantees that the potential and its first derivative are zero at r ij = R c σ ij . Here R c = 2 1/6 and 2.5 correspond to WCA and LJ systems, respectively. For 2D polydisperse systems (N = 4096 particles and particle size polydispersity ∆ = 13%), both macroscopic and microscopic properties are characterized. We set the number density ρ = 0.89, so that both WCA and LJ systems have a positive pressure at zero temperature. The simulation methods are the same as soft repulsive systems, except that a smaller time step is used to ensure a good convergence of numerical integration. Without loss of generality, we simulate three initial states in the isoconfigurational ensemble. For 3D systems as listed in Supplementary Table 1, only the macroscopic properties are characterized. In Methods of the main text, we have defined the structural order parameters for hard-sphere-like systems, for which particles have a well-defined surface, as given by the range of the repulsive interaction. In such systems, it is clear that the main driving force for structural ordering comes from packing and excluded volume effects, which are Supplementary (3)- (5) in the main text). The situation can be complex, considering the recent studies which suggest a nonperturbative effect of the attractive forces [7,8]. However, it has been shown that an LJ liquid at typical condensed-phase state points can be categorized as a simple liquid in the sense that the intermolecular interactions may be ignored beyond the first coordination shell [9]. This means that its structure is dominated by what may be termed "packing effects". Therefore, it is reasonable to expect that our structural order parameters characterizing the local packing capability also work in LJ systems. Motivated by the fact that, over a wide range of density in the typical condensed phase, the pair correlation functions g(r) of WCA and LJ systems resemble each other [7], we define the effective radius of an LJ particle at the minimum of the interaction potential. So σ i,eff = 2 1/6 σ i for both WCA and LJ systems. It is also possible to define the effective radius on the repulsive side of the potential at which the potential energy with respect to the minimum equals the thermal energy or to pick up the position of the first peak of g(r/ σ ) at R 1 and set σ i,eff = R 1 σ i . We find that the values of Θ are quite insensitive to these choices because Θ measures orientational rather than translational correlations. Therefore, we use σ i,eff = 2 1/6 σ i in this study.
Hard disk systems. We have also investigated 2D polydisperse hard disks systems and confirmed the generality of our findings. The mixture of N = 10000 disks is simulated with event-driven molecular dynamics using the DynamO package [3]. Similar to the other polydisperse systems studied in this work, the particle size is extracted from a Gaussian distribution with polydispersity ∆ = 11%, for which no transition to a hexatic phase is observed. All particles have the same mass. The unit of length is set by the average disk diameter σ and the time scales are reported using the event-driven unit. All simulations are run at fixed densities ρ = N/L 2 in square boxes. After the equilibration runs, we generate the isoconfigurational ensemble from the equilibrated configurations. Without loss of generality, we have simulated three independent initial states at each density, and for each initial state, 100 trajectories are simulated. SUPPLEMENTARY NOTE 2. ADDITIONAL RESULTS IN HARMONIC SYSTEMS 1. Analysis of dynamics. Here we characterize the structure relaxation in both 2D and 3D glass-forming liquids [2]. In 2D case, we use the position of particle j relative to its n j neighbouring particles l, r j (t) = r j (t) − l r l (t)/n j , to characterize the dynamics, which helps to remove the long-wavelength Mermin-Wagner fluctuations [15][16][17]. In Supplementary Figure 1, we shown the temperature dependence of self-intermediate scattering functions F s (k, t) in both 2D and 3D, where typical features of glassy dynamics, i.e. the two-step relaxation, are generally observed.
The structure relaxation time τ α is measured from the time decay of the self-intermediate scattering function: F s (k, τ α ) = e −1 . Supplementary Figure 2 shows the temperature dependence of τ α for the six major systems under study. The drastic dynamical slowing down is well described by the Vogel-Fulcher-Tammann (VFT) relation τ α = τ 0 exp[DT 0 /(T − T 0 )], which suggests a divergence of τ α at the ideal glass transition temperature T 0 . The fitting parameters are listed in Table 1 of the main text.
We may also plot τ α as a function of the inverse of temperature 1/T , as shown in Supplementary Figure 3, from which a crossover from the Arrhenius to non-Arrhenius behaviours is identified as the onset temperature T on of sluggish glassy dynamics. The appearance of the non-Arrhenius behaviour suggests a growth of the activation energy, which is a clear indication of cooperative dynamics typical for the so-called "fragile" glass-forming liquids. From the potential energy landscape (PEL) formalism, T on also signals the onset of a significant influence from the underlying inherent structures and a change in the manner of exploration of the PEL [18].
Supplementary Figure 4 further shows power-law fittings of the temperature dependence of τ α , as predicted by the      mode-coupling theory (MCT). MCT captures the initial stage of slowing down upon cooling and predicts a divergence of τ α toward T mct , which however does not take place in reality. As seen from Supplementary Figure 4, the growth of τ α eventually deviates from the MCT prediction at low temperatures. Nevertheless, T mct provides the reference information on the temperature range, or the degree of supercooling, at which one is looking. This analysis suggests that our simulations cover the range of temperature down to around T mct . We note that the determination of T mct from power-law fittings is purely empirical and not rigorous. Within numerical precision, our result that T mct = 5.4 × 10 −4 for the 3D binary system is consistent with Ref. [19], where T mct = 5.2 × 10 −4 was found. The slight difference may originate from the fact that big and small particles are analyzed separately in Ref. [19] but all together in current work.
2. Glassy structure formation towards glass transition. Corresponding to Fig. 1a in the main text, in Supplementary Figs. 5 and 6, we show the glassy structure formation for the rest five systems in 2D and 3D, respectively. Universal behaviours are observed, including linear temperature dependence of structural order of instantaneous states in the supercooled regime, and a constant structural order of inherent states in both the hightemperature simple liquid and low-temperature glass regimes. This confirms the general validity of our results. Fig. 1a of the main text, here we further unveil the role of structural order in glassy dynamics by a close comparison of instantaneous and inherent states. In Supplementary Figure 7a, we show the average number of neighbour change N c between the corresponding configurations in instantaneous and inherent states. Here the neighbouring particles are defined by the radical Voronoi tessellation [20]. We can see that N c turns to increase rapidly above T g , which is a clear indication that local structural order is destructed under thermal noise. Supplementary Figure 7b further shows that the more disordered structures are less resistive to destruction due to thermal noise. Since Θ can be either bigger or smaller than Θ IS , it is essential to note that the thermal fluctuation may either increase or decrease local structural order in instantaneous configurations, and hence its role is highly nontrivial in supercooled liquids. The substantial structural difference between instantaneous and inherent states in the supercooled regime (see Fig. 1a in the main text) indicates that the inherent structures do not properly reflect the real liquid structures, which are under intrinsic influence of thermal fluctuation (or, entropy). Thus, we argue that it is the instantaneous state rather than the inherent one that we should look at for a quantitative relation between structure and dynamics. We note that this fact has not been taken seriously in many previous works.

A close comparison of instantaneous and inherent states. Corresponding to
We further illustrate this point via a comparison of structure-dynamics correlations based on instantaneous (thermalized) and inherent states. As shown in Supplementary Figure 8 for a 2D polydisperse system (∆ = 13%) at T = 0.0018, the peak structure-dynamics correlation is increased by more than 20% and reaches C r = 87% by starting from instantaneous (thermalized) states. Such a high structure-dynamics correlation is unprecedented and crucial to the determination of the quantitative VFT-like relation. Therefore, even though it might appear as a minor change that we switch from inherent structure to instantaneous (thermalized) structure to define the structural order parameter, it has a very fundamental physical significance. The important point is that the inherent structure is a state that is never visited by a system in the liquid state. So the structural order parameter defined for instantaneous structures should be regarded as a true measure of the liquid structure. Because of this feature, this order parameter can have the role of a genuine control parameter, or an effective thermodynamic intensive variable. and inherent states. For 2D polydisperse system (∆ = 13%) at T = 0.0018, we characterize the structural order in thermalized states and the corresponding inherent states, and then quantify the structure-dynamics correlation respectively. We find that the peak correlation is increased by 21.2% and reaches Cr = 87% by starting from thermalized states, which is unprecedented in previous studies. We calculate the total entropy through thermodynamic integration starting from the low-density ideal-gas limit (T = 0.01, ρ → 0) along the isotherm T = 0.01, up to the density studied ρ = 0.675. a and b, Density dependence of pressure P and the excess pressure over the ideal-gas value P id , respectively. The solid lines indicate the ideal-gas law P ∼ ρ and the predicted virial correction P − P id ∼ ρ 2 in a and b, respectively. The good consistence between our data and the theoretical prediction suggests that we have properly accessed the ideal-gas limit. c, Temperature dependence of the potential energy U at fixed density ρ = 0.675. This is used for the calculation of the total entropy at target temperatures [10][11][12][13].

Relation between configurational entropy and our structural order parameter.
Here we explore the relation between configurational entropy S conf and our structural order parameters. In particular, we focus on binary glass formers to avoid possible issues related to the polydisperse distribution of particle sizes [21][22][23], the resolution of which is beyond the scope of the current work. S conf under certain conditions is determined by the multiplicity of local potential energy minima sampled by the liquid. Following the well-established method in the literature [10][11][12][13], we evaluation S conf by calculating the difference of the total entropy and the vibrational entropy of the basins We calculate S tot via thermodynamic integration starting from the ideal-gas reference point. For 3D harmonic systems under study, we first integrate from the ideal-gas limit (ρ → 0, T = 0.01) to the reference state (ρ = 0.675, T = 0.01) along the isotherm. The density dependence of pressure P shown in Supplementary Figs. 9a and 9b confirms that we have properly accessed the ideal-gas limit. The total entropy at the target temperature along the isochoric path can . This observation is consistent with the previous studies in 2D [14]. d, Following the generalized Adam-Gibbs relation, the structural order parameter 1/(Θ − Θ0) is plotted as a function of (1/T S conf ) α , with the same α as c. The solid line suggests a linear relation between them.
then be calculated using the temperature dependence of the potential energy U shown in Supplementary Figure 9c. For all results presented in the following, the vibrational entropy is calculated by approximating each basin as a harmonic well [10][11][12][13]24]. This approach has been confirmed to be valid at low temperatures [10][11][12][13]24], which agrees reasonably well with other methods [12,23]. We refer to Refs. 10-13 for the detailed formulism. Supplementary Figure 10a shows the temperature dependence of S tot and S vib in 3D binary systems, from which we calculate S conf according to Eq. (2). Supplementary Figure 10b shows the Adam-Gibbs plot of structure relaxation time τ α as a function of 1/T S conf , which confirms the Adam-Gibbs relation τ α = τ 0 exp(A/T S conf ) in this system [14,24,25]. In comparison with Eq. (2) in the main text, this result motivates a direct relation between our structural order parameter and configurational entropy. This point is confirmed in Supplementary Figure 10c.
We further explore the relation between S conf and our structural order parameter Θ in 2D binary systems. Previously, it has been shown that τ α does not following the standard Adam-Gibbs relation in 2D, but rather a generalized Adam-Gibbs relation τ α = τ 0 exp[(A/T S conf ) α ] with α = 1 [14]. As shown in Supplementary Figs. 11a-c, our results are also described by the generalized form of Adam-Gibbs relation. Accordingly, we find a linear relation between our structural order parameter 1/(Θ − Θ 0 ) and (1/T S conf ) α . Although a fundamental understanding of such a relation lacks at present, which is interesting to be explored further [26], this result nevertheless point to an intriguing connection between structural order as measured by our order parameter and configurational entropy.

5.
Relation between Θ and local structure entropy. Here we study the local structure entropy s 2 at a particle level, which is a local version of the two-body translational correlation contribution to the excess entropy [27]: where k B is the Boltzmann constant, ρ is the number density, and g i (r) is the pair correlation function between particle i and the other particles. For simplicity, we focus on a 2D polydisperse system (∆ = 13%), but generality is also expected for other systems. Supplementary Figure 12a shows the spatial distribution of s 2 , which is highly correlated with that of bare Θ (see Fig. 3d in the main text). Such a direct correlation can be more easily visualized after spatial coarse-graining of both quantities, as shown in Supplementary Figure 12b for s 2,CG and Fig. 3e in the main text for Θ CG . This result indicates that similar to hard-sphere crystallization, sterically favoured structures sacrifice configurational entropy to gain more vibrational entropy and thus low free energy locally. This supports free-energy-driven formation of locally favoured structures [28]. The temperature dependence of P (s 2 ) is shown in Supplementary Figure 12c. A shift to a lower value of s 2 is seen with decreasing temperature, indicating a growth of structural correlation. From the inset, we can see that ordered particles contribute to the development of the tail of s 2 towards a negative direction. At the same time, we note that local structure entropy is not local in a strict sense Supplementary Figure 12: Local structure entropy. a,b, Spatial distributions of bare two-body excess entropy s2 and that coarse-grained at L = 3.8, s2,CG, respectively, for a 2D PM (∆ = 13%) at T = 0.0018. See the correspondence between these two panels and Figs. 3c-e in the main text. c, Probability distribution of s2 for a range of temperatures (decreasing temperature from right to left, with T = 0.0037, 0.0028, 0.0023, 0.002, and 0.0018). Inset: For T = 0.0018, probability distribution of s2, P (s2), for the whole system (blue solid), half particles with smallest Θ hence most ordered (black dashed), and the other half particles with largest Θ hence most disordered (black dotted).
(by taking the upper integration limit of Eq. (3) to long-distance). This makes a quantitative comparison between s 2 and our order parameter illusive. 6. Nonlocal scenario for structure relaxation. Corresponding to Fig. 3b in the main text, in Supplementary Figs. 13 and 14, we show the correlation between the microscopic relaxation time (τ α for each particle) and the structural order as functions of the coarse-graining length L for 2D and 3D systems, respectively. We generally observe that C r significantly increases initially with increasing the coarse-graining length and maximizes around a temperature-dependent length scale, which is identified as the characteristic static correlation length ξ of the underlying structure. We notice from a close inspection of the correlation curves C r (L) that they are quite smooth in 2D while more fluctuated in 3D. We speculate that this is due to the large number of particles in each neighbour shell in 3D, and hence a significant influence from the density fluctuations when the coarse-graining length covers an additional neighbour shell. An improved version of coarse-graining might be necessary to smooth out such fluctuations. 7. Relationship between structure relaxation time τ α and static correlation length ξ. An important question concerning the nature of glass transition is whether the drastic dynamical slowing down originates from the growth of the underlying static order. In Supplementary Figure 15, we show the relation between structure relaxation time τ α and the static correlation length ξ for both 2D and 3D systems. In general, the relation τ α = τ 0 exp[D(ξ/ξ 0 ) d/2 ] with d being the spatial dimension gives a reasonable description of the data, suggesting a common structural origin of slow glassy dynamics. This observation is consistent with previous studies [27,28] and can be rationalized in an Ising-type critical scenario for glass transition [27][28][29]. We note that this relation is also consistent with the RFOT scenario [30][31][32], although such a relation is usually expected for below T mct in that scenario.
8. Microscopic relation between structural order and relaxation in 3D. Corresponding to Figs. 4c and 4d in the main text, in Supplementary Figure 16, we show how the microscopic structure relaxation is related to the structural order in 3D binary and polydisperse (∆ = 13%) systems. The same behaviours as the 3D polydisperse systems with ∆ = 8% are observed that, after coarse-graining, the microscopic relation between structural order and relaxation tends to follow the macroscopic one in the case of deep supercooling.
9. Intrinsic fluctuations of microscopic relaxation dynamics. In Supplementary Figure 16 and also Fig. 4 in the main text, we have shown the correspondence between single-particle structural order and microscopic relaxation time. Good correlations are reached at lowest temperatures (see also Supplementary Figs. 13 and 14), but certain degrees of fluctuations, meaning a certain width of scattering of the microscopic data points around the macroscopic relation, are seen. To understand the origin of such fluctuations, here we study the microscopic relaxation dynamics in more detail. In addition to the self-intermediate scattering function for the whole system and each particle defined through the isoconfigurational average (see Methods in the main text), we calculate the self-intermediate scattering function for each particle and in each run in the isoconfigurational ensemble. This additional analysis would provide further information about the intrinsic fluctuations of single-particle relaxation dynamics. For particle j in the a th run, F j,a s (k, t) = cos ik · [r a j (t) − r a j (0)] , with the superscript a indicating the a th trajectory. Supplementary Figure  17 shows F s (k, t) averaged over the whole system (dark blue circles), for an immobile (green squares, left) or mobile (red squares, right) particle averaged within the isoconfigurational ensemble, and for that particle in different runs of the isoconfigurational ensemble (fluctuating curves in the background). Several interesting observations can be made from Supplementary Figure 17: (1) Not only that each particle may relax differently from the global average, but Fs(k, t) for 2D polydisperse system (∆ = 13%) at T = 0.0018. a, Fs(k, t) averaged over the whole system, for an immobile particle averaged within the isoconfigurational ensemble, and for that particle in different runs of the isoconfigurational ensemble (fluctuating curves in the background, 10 different realizations are shown). This immobile particle has a structural order ΘBare = 0.07344 and ΘCG = 0.09228. b, The same plot as in a but with the selected particle being mobile. This mobile particle has a structural order ΘBare = 0.1224 and ΘCG = 0.1166. Note that larger Θ means less ordered. We see that not only each particle may relax differently from the global average, but also the same particle with exactly the same initial static structural order may relax differently due to different initial velocity fields in the isoconfigurational ensemble. See text for more discussions.
also the same particle with exactly the same initial static structural order may relax very differently due to different initial velocity fields in the isoconfigurational ensemble. (2) The immobile particle has a high structural order, and its relaxation dynamics is more determinate. The structure relaxation is found to take place with a narrow time window in different runs. (3) In comparison, the mobile particle has a low structural order, and the relaxations significantly differ in different runs. This observation can be rationalized from the competition between structural order (determinate effect) and thermal fluctuation (stochastic effect). Thermal fluctuation plays an intrinsic role in the microscopic dynamics; therefore, in the case of weak local structural order, the particle tends to relax randomly. Consequently, the microscopic relaxation time that we measure for each particle should be understood in a statistical sense, which is intrinsically affected by the thermal fluctuations. Meanwhile, we should mention that our structural order parameters are adequate, but by no meaning, perfect and coarse-graining is a working but approximate procedure. These also contribute to the scattering of data points in Fig. 4 in the main text and Supplementary Figure 16, in addition to the intrinsic thermal fluctuations.