Emergent solidity of amorphous materials as a consequence of mechanical self-organisation

Amorphous solids have peculiar properties distinct from crystals. One of the most fundamental mysteries is the emergence of solidity in such nonequilibrium, disordered state without the protection by long-range translational order. A jammed system at zero temperature, although marginally stable, has solidity stemming from the space-spanning force network, which gives rise to the long-range stress correlation. Here, we show that such nonlocal correlation already appears at the nonequilibrium glass transition upon cooling. This is surprising since we also find that the system suffers from giant anharmonic fluctuations originated from the fractal-like potential energy landscape. We reveal that it is the percolation of the force-bearing network that allows long-range stress transmission even under such circumstance. Thus, the emergent solidity of amorphous materials is a consequence of nontrivial self-organisation of the disordered mechanical architecture. Our findings point to the significance of understanding amorphous solids and nonequilibrium glass transition from a mechanical perspective.

1. Basic analysis of glassy dynamics. Here we characterize the structure relaxation in glass-forming liquids [1]. In particular, we show results of the 2D harmonic systems with N = 4096 particles focused in this work. We use the position of particle j relative to its n j neighboring particles l, r j (t) = r j (t) − l r l (t)/n j , to characterize the dynamics in 2D [2,3], which helps to remove the long-wavelength Mermin-Wagner fluctuations [3][4][5][6]. Then the self part of the intermediate scattering function is given as Here r j (t) is the relative position of particle j at time t, k is the scattering vector with k = |k| close to the first peak of the static structure factor, and brackets indicate an ensemble average. The mean square displacement (MSD) is another probe of the dynamics: (2) Figure 1 shows F s (k, t) and MSD for a wide range of temperature. The emergence of the intermediate-time plateau of both characterisations in the supercooled regime signals the transient localisation, i.e., caging, a universal dynamic feature of supercooled liquids.
The structure relaxation time τ α is measured from the time decay of the self-intermediate scattering function: F s (k, τ α ) = e −1 . Figure 2 shows the temperature dependence of τ α . From Fig. 2a, we find that 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 hypothetical divergence of τ α at the ideal glass transition temperature T 0 ≈ 0.63. We may also plot τ α as a function of the inverse of temperature 1/T , as shown in Fig. 2b, from which a crossover from the Arrhenius to non-Arrhenius behaviours is identified as the onset temperature T on ≈ 2.1 of sluggish glassy dynamics. The appearance of  Here v(t) = {vi(t)} is the velocity field of the whole system. From the relaxation curves, the period of oscillation τo ≈ 5 can be identified, which has a weak temperature dependence.
the non-Arrhenius behaviour suggests a growth of the activation energy, which is an indication of cooperative dynamics typical for the so-called "fragile" glass-forming liquids. Figure 2c further shows a power-law fitting of the temperature dependence of τ α , as predicted by the mode-coupling theory (MCT). The MCT temperature is empirically estimated to be T mct = 1.21.
In addition, we characterise the velocity correlation in a temperature range covering both simple liquids and deep in the solid phase, as shown in Fig. 3. From the oscillation of velocity correlation, we estimate the period of oscillation of the short-time vibrational motions τ o ≈ 5, which shows only a weak temperature dependence.
2. Examination of finite-size effects. Here we examine the finite-size effect on the long-range stress correlation in thermal amorphous solids. In principle, the long-range character cannot be proven but rather suggested by numerical simulations with finite-size systems. In this context, the method of finite-size scaling, used extensively in numerical studies of critical phenomena, is particularly suited for investigations of the emergence of the long-range correlations. In the main text ( Fig. 4), we show through finite-size scaling that a well-defined percolation transition of the forcebearing particles exists in the large system size limit. The result also suggests that the disordered solid states below the non-equilibrium glass transition temperature are not affected by system sizes, so that a general relation between the long-range stress correlation and the percolating force-bearing network is observed. Here we further illustrate this point in Fig. 4. The angle-averaged shear-stress correlation C σxy (r) for a range of system sizes from N = 4096 to N = 262144 are plotted, which in turn correspond to linear sizes from L ≈ 72 to L ≈ 579. As we can see, C σxy (r) from a smaller system closely matches those from larger ones. This suggests that finite-size effects can be neglected in the characterisation of stress correlations in the solid phase. Therefore, we mainly focus on the systems with N = 4096 particles, unless otherwise specified, in order to extensively explore a wide range of cooling rates, the convergence of numerics with very large ensembles, and additional properties like the shear modulus.
3. The necessity of large ensemble. It is worth pointing out that the decay of power-law correlation functions 1/r d in d dimensions is relatively fast from the viewpoint of numerics. In order to obtain clear evidence of such longrange correlations, a large ensemble of independent configurations are necessary to suppress statistical fluctuations. We note that, particularly when entering the glassy regime, even extremely long-time average does not work so efficiently to sample the configuration space. To emphasise this point, we compare the results calculated by using a different number of configurations in Fig. 5. At low temperature (T = 0.1) deep in the solid phase, we find that the data from a smaller ensemble of configurations is more fluctuated, although the power-law behaviour is still noticeable (Fig. 5a). At higher temperature (T = 1.0), slightly below the γ-dependent T g , we see that while the power-law behaviour is still clear for a large ensemble of configurations, it is basically hidden behind the noise for a small ensemble of configurations (Fig. 5b). In Fig. 5c, we further show C σxy (r) from an ensemble of N conf = 100 configurations for a range of temperatures covering both above and below T g ≈ 1.3 (see Fig. 1c in the main text). The results are essentially indistinguishable for different temperatures, so the crossover at T g cannot be resolved. However, if the power-law scaling is plotted unadvisedly together with the data as in Fig. 5c, the long-range character may be incorrectly claimed (we note that this can be found in the literature). ∆P ∼ ∆r is observed for small ∆r, which is suggestive of a harmonic regime close to IS. A strong deviation can be observed when the system is slightly away from IS, exhibiting similar behaviours as the potential energy shown in Fig. 2c in the main text. Therefore, this observation further confirms the giant anharmonic effects in thermal amorphous solids.
5. Initial development of long-range stress correlation after an instant quench. Until now, we have focused on the static properties related to the stress correlations in amorphous solids. For a better understanding of the long-range character as an emergent property of the metastable glass states, it is also of interest to study how it emerges, starting from an intrinsically random state. To this aim, we bring the system from a high-temperature liquid state to the target temperature by an instant quench and follow the evolution of stress correlation during further equilibration. The evolution of C σxy (r) deep in the solid phase after a quench is shown in Fig. 7a. We find that a power-law feature emerges already at t = 10 which is in the fast β (or, cage-rattling) time scale and corresponds to only several periods of oscillations (see Figs. 1 and 3). Nevertheless, the correlation continues to grow with time, reflecting the underlying optimisation of the force network for partial mechanical equilibrium. To better quantify this behaviour of time evolution, we follow the integrated shear-stress correlations for a range of target temperatures, as shown in Fig. 7b. In particular, the curve corresponds to  Figure 6: Giant anharmonic effects in thermal amorphous solids seen from pressure. Corresponding to Figs. 2b and 2c in the main text, here we shown the relationship between the potential pressure difference ∆P and direct distance ∆r with respect to IS during energy minimizations for different initial temperatures Ti. The solid black lines indicate the linear relation ∆P ∼ ∆r. Note that the average pressure in IS decreases from 0.023 to 0.021 when the initial temperature is lowered from 3.0 to 0.2. is the surprisingly short time scale required for the development of long-range (or quasi-long-range) stress correlations, which strongly suggests their relevance when considering the properties of realistic metastable glassy states [7]. This result indicates that the establishment of partial mechanical equilibrium after an instant quench in temperature, i.e., solidification, is a fast process facilitated by vibrational motions rather than diffusive motions. It also validates that our simulations, which cover much longer time scales, are adequate to study this process.
At temperatures close to the non-equilibrium glass transition, we find nonmonotonic evolutions of the correlation, which is a clear sign of the ageing effect. This is most clearly seen from the curve of open orange circles at T = 1. in Fig. 7b (see Fig. 1c in the main text that at this temperature, the system is expected to go back to the liquid state after long-time relaxations). This result unveils the real-space mechanical aspect of physical ageing, which was previously characterised by geometrical structural changes [8] or macroscopic quantities [9]. This is an interesting point which definitely deserves careful studies in the near future.
6. Further percolation analysis of the force-bearing network. Here we perform further percolation analysis of the force-bearing network, concerning the nature of observed percolation transition within the percolation theory [10]. First, we characterise the probability distribution of cluster sizes at the critical occupation fraction f c , i.e. at the non-equilibrium glass transition. As shown in Fig. 8, the cluster size shows a power-law distribution n s (s) ∼ s −τ with τ = 1.85, which is significantly different from the random-percolation exponent τ rand = 187/91 ≈ 2.05 [10]. As indicated by the dashed line, such a difference may not be accounted for by numerical uncertainties. It is known that the critical exponents of random percolation are very robust in changing the percolation model details [10,11]. They do not depend on the short-range structure of the lattice (e.g., square or triangular) or the type of percolation (site, bond or even continuum) [10]. In other words, only in the presence of long-range correlations that the critical exponents of percolation transition can be changed [11]. Therefore, the value of the scaling exponent τ different from τ rand suggests that the observed percolation of force-bearing particles does not belong to the universality class of random percolation but is controlled by certain long-range correlations. Here, we expect that the emergence of longrange stress correlation in thermal amorphous solids across the non-equilibrium glass transition may be the underlying mechanism of this unusual nature of percolation transition.
We further characterise the topology of the largest cluster across the percolation transition by plotting its size s 1 as a function of the system size N in Fig. 8b. At f c , as indicated by the filled diamonds and solid line, we find that s 1 ∼ N d f /d with d f ≈ 1.97. Within numerical precision, the observed d f is not different from d = 2 (although, within numerical precision, we cannot rule out the possibility that it coincides with d f,rand = 91/48 ≈ 1.9 of random percolation). This is consistent with the behaviour of the relative size of the largest cluster ψ = s 1 /N shown in Fig. 4c of the main text. There, the finite-size scaling suggests an abrupt formation of system-spanning force-bearing network at f c in the large system size limit with ψ c ≈ 0.4. The finite value of ψ c suggests that the percolating cluster is compact, meaning d f ≈ 2. Below f c , i.e., at high temperatures in the liquid phase, the growth of s 1 does not follow the system size. Therefore, stress correlation should be short-ranged. Above f c , i.e., at low temperatures in the solid phase, s 1 is proportional to the system size, as expected from a compact nature after (and away from) the percolation transition.
Finally, we mention that similar observations, namely τ < 2 and d f = 2, have been made in several percolation models, e.g., the no-enclaves percolation [12,13] and the percolation of cluster holes [14]. In addition, the abrupt formation of a percolating cluster at f c with a finite ψ = s 1 /N is also lively discussed in the so-called "explosive percolation" models [15][16][17][18][19]. Whether there is a deep link between our physical systems and these mathematical models is an interesting topic to be explored in future studies.

SUPPLEMENTARY NOTE 2. RESULTS FOR ADDITIONAL SYSTEMS
1. The emergence of long-range stress correlation in systems with different volume fractions. The emergence of long-range stress correlation driven by temperature is characterised in three more 2D harmonic systems with φ = 0.86, 0.95 and 1.0, in addition to the one with φ = 0.91 on which we have focused in the main text. Figure 9 shows the temperature dependence of integrated shear-stress correlations for all these four systems at a cooling rate of γ = 10 −5 . As expected, the non-equilibrium glass transition signalled by the sudden growth of shear stress correlation upon cooling is generally observed, with a higher T g for larger φ. This result confirms that our observation is general with respect to the volume fraction and not restricted to systems close to the jamming transition.  2. The emergence of long-range stress correlation driven by volume fraction. In addition to the nonequilibrium glass transition driven by temperature, which has a close correspondence to that in atomic systems, here we further study the glass transition driven by volume fraction (equivalently, density or pressure), which is more closely related to the colloidal systems [20,21]. Similar to the case of temperature control (see Methods in the main text), the systems are first equilibrated at φ = 0.7 and T = 1.0 for t = 10000 and then compressed to φ = 1.0 in a quasi-continuous fashion. The compression rate is defined as γ φ = ∆φ/∆t, where the compression step is fixed to ∆φ = 5 × 10 −3 and ∆t is the residence time at each volume fraction, which controls γ φ . Figure 10 shows the φ-dependence of integrated shear-stress correlations for different compression rates γ φ . It is seen that the nontrivial shear-stress correlation emerges across a γ φ -dependent glass transition, which is essentially similar to that driven by temperature (see Fig. 1c in the main text). This result confirms that the emergence of long-range stress correlation is a general consequence of the non-equilibrium glass transition, no matter it is driven by temperature or volume fraction.
3. Emergence of long-range stress correlation in 2D Lennard-Jones systems. In the main text, we have focused on the system with a finite-range repulsive interaction (see Methods of the main text), which has been widely studied in numerical simulations and shown to behave as canonical glass formers [22,23]. The simplicity of this model allows us to carry out extensive analyses to explore carefully the parameter dependence of results, the convergence of numerics and the finite-size effects. However, our major findings, that the long-range stress correlation emerges as a consequence of the non-equilibrium glass transition even under giant anharmonic effects, are not restricted to this simplified model but rather generic for glass formers with different interactions. To support this point, we consider an equimolar binary mixture system in which particles i and j interact through a Lennard-Jones (LJ) potential: 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 we set R c = 2.5, particle diameter ratio=1.4, and consider a 2D system with a number density ρ = 0.61, which ensures positive pressure at zero temperature. Temperature is in unit of /k B . This and similar models have been widely studied as models of atomic glass formers [24][25][26][27][28], e.g., the famous Kob-Anderson model to describe amorphous Ni 80 P 20 [24] and the equimolar binary version for metallic glass Cu 50 Zr 50 [25,27,28]. Figure 11 shows the temperature dependence of integrated shear-stress correlations in LJ systems. Nontrivial stress correlations are found to emerge across the cooling-rate dependent glass transition, with the same features as in the  Figs. 2b and 2c in the main text, finite-temperature configurations are relaxed to the inherent states (IS) of potential energy landscape (PEL) using energy minimisation (EM) methods. Top panels are for 2D Lennard-Jones systems. a, Relationship between contour distance ∆s and direct distance ∆r with respect to IS during EM for different initial temperatures Ti. Power-law fits are labelled by their corresponding exponent. A crossover from a linear relation ∆s ∼ ∆r close to IS to a fractal one ∆s ∼ ∆r 1.35 away from IS is observed. b, Relationship between energy ∆E with respect to IS and ∆r. The solid black lines indicate the quadratic (harmonic) relation ∆E ∼ ∆r 2 . A strong deviation from the harmonic approximation is observed at approximately the same ∆r where the crossover into a fractal-like PEL happens. Bottom: c,d, The same analyses as in a and b for 3D harmonic systems. A crossover from a linear relation ∆s ∼ ∆r close to IS to a fractal one ∆s ∼ ∆r 1.35 away from IS, and correspondingly a strong deviation from harmonic approximation is observed. Taken together, these results suggest a fractal potential energy landscape with a universal fractal dimension dPEL = 1.35 in different disordered systems. We emphasise that the thermal system (initial point of each curve in the top-right end) in physically relevant temperatures are completely out of the harmonic regime.
harmonic systems (see Fig. 1c in the main text). This result strongly suggests that the long-range stress correlation is a generic feature emerging when the system falls out of equilibrium and becomes trapped in the metastable glass state, which is independent of the interacting potential.
4. The fractal-like potential energy landscape in 2D Lennard-Jones and 3D harmonic systems. In Fig. 2 of the main text, we show that the potential energy landscape (PEL) of a 2D harmonic system has a fractal geometry, with a nontrivial fractal dimension d PEL = 1.35. Here we confirm that such a fractal-like PEL also exists in 2D Lennard-Jones and 3D harmonic systems. Most strikingly, as shown in Figs. 12a and 12c, the fractal dimension of PEL in both systems are found to be approximately d PEL = 1.35, the same as the 2D harmonic system. This result strongly suggests a universal fractal character of PEL, which represents an important aspect of the peculiar nature of amorphous solids. Moreover, a strong deviation from the harmonic energy expansion is found to coincide with the crossover into the fractal-like PEL, as shown in Figs. 12b and 12d. In conclusion, the general existence of giant anharmonicity controlled by the fractal-like PEL definitively rules out the effectiveness of harmonic approximations to understand the thermal amorphous solids from their inherent states.  Figure 13: Percolation of force-bearing particles gives rise to long-range stress correlation in 2D Lennard-Jones systems. a, Temperature dependence of percolation probability P for different γ. The γ-dependent nonequilibrium glass transition is signalled by the percolation of force-bearing particles. Inset: Collapse of P from different γ when replotted as a function of the fraction of force-bearing particles f . b, Temperature dependence of the relative size of the largest force-bearing cluster ψ = s1/N for different γ, with s1 being the number of particles involved in the largest cluster. A strong similarity can be seen in comparison with the integrated stress-stress correlation in Fig. 11. Inset: Collapse of ψ from different γ when replotted as a function of f . c, Integrated shear-stress correlations for different cooling rates γ are plotted as a function of the relative size of the largest force-bearing network ψ = s1/N . The filled circle and lines in purple have the same meaning as those in Fig. 4a of the main text. The nice collapse of the data at different γ suggests that the shear-stress correlation, more specifically, its far-field behaviour, is uniquely controlled by the largest force-bearing network. This result points to the general mechanism that the percolation of force-bearing particles gives rise to long-range stress correlation, irrespective of the interaction potential.  Figure 14: Finite-size scaling analysis for the percolation of force-bearing particles in 2D Lennard-Jones systems. a, Percolation probability P as functions of the fraction of force-bearing particles f for different system sizes N at γ = 10 −6 . Inset: Scaling collapse of P from different system sizes. Here d = 2 is the spatial dimension, ν = 4/3 is the scaling exponent, and fc = 0.48 is the critical occupation fraction. b, Corresponding to panel a, the relative size of the largest cluster ψ = s1/N as functions of f for different system sizes. The vertical and horizontal arrows indicate fc = 0.48 and the corresponding ψc, respectively. Inset: Scaling collapse of ψ from different system sizes below the percolation transition, using the same fc and ν as the inset of a. The horizontal dashed line indicates ψc ≈ 0.4 at the percolation threshold fc. Taken together, both panels a and b indicate that, when N → ∞, an abrupt formation of system-spanning force-bearing network takes place across the non-equilibrium glass transition.
5. Percolation analysis of the force-bearing network in 2D Lennard-Jones systems. Here, we explore a possible generalisation of our force network analysis in 2D Lennard-Jones systems. Due to the long-range attractive tail in the interaction, each particle may interact with many neighbours, making the force network highly complicated. However, it was shown through a study of boson peak in Lennard-Jones systems that a relatively few strong repulsive interactions (stiff contacts) play a dominate role, whereas the more numerous weaker interactions can be treated as perturbations [29]. This means that, although intrinsically complex, the identification of strong backbone may capture the main effects of the force network even in Lennard-Jones systems. Our definition of the force-bearing network is indeed based on the same physical consideration. However, we find that it is necessary to extend the criterion to F str < 1.28 to properly identify the force-bearing particles, which may be due to the renormalising effect of long-range attractions. Figure 13a shows the temperature dependence of the percolation probability for three different cooling rates. It is seen that the percolation of force-bearing particles takes place at a γ-dependent glass transition temperature, in good agreement with the emergence of the long-range stress correlation shown in Fig. 11. In addition, the temperature dependence of the relative size of the largest force-bearing cluster ψ = s 1 /N is shown in Fig. 13b, where a close similarity compared with Fig. 11 can be observed. In Fig. 13c, We further plot the integrated shear-stress correlations as a function of the relative size of the largest force-bearing network ψ and find a nice collapse of the data from different cooling rates. In conclusion, these results suggest that the long-range stress correlation in LJ systems is also uniquely controlled by the underlying percolating force-bearing network, which therefore points to a general understanding of the non-equilibrium glass transition from the mechanical perspective.
Finite-size scaling analysis for the percolation of force-bearing particles is carried out to support the plot of Fig. 13. As shown in Fig. 14a, the transition of percolation probability P from 0 to 1 as a function of the fraction of forcebearing particles f becomes steeper with increasing the system size N . Scaling collapse of P is realised in the inset of Fig. 14a, indicating a well-defined percolation transition at f c = 0.48 in the large system size limit. Correspondingly, the relative size of the largest force-bearing network ψ as a function of f is shown in Fig. 14b, with the finite-size scaling collapse of ψ below f c given in the inset. This result again indicates that a system-spanning force-bearing network with ψ ≈ 0.4 emerges abruptly at f c as N → ∞. Therefore, a well-defined percolation transition of the force-bearing network in the large system size limit is also expected in LJ systems. Taken together results from both harmonic and LJ systems, the percolation of the force-bearing network that gives rise to the long-range stress correlation is speculated to be the universal physical mechanism underlying the non-equilibrium glass transition.