Common mechanism of thermodynamic and mechanical origin for ageing and crystallization of glasses

The glassy state is known to undergo slow structural relaxation, where the system progressively explores lower free-energy minima which are either amorphous (ageing) or crystalline (devitrification). Recently, there is growing interest in the unusual intermittent collective displacements of a large number of particles known as ‘avalanches’. However, their structural origin and dynamics are yet to be fully addressed. Here, we study hard-sphere glasses which either crystallize or age depending on the degree of size polydispersity, and show that a small number of particles are thermodynamically driven to rearrange in regions of low density and bond orientational order. This causes a transient loss of mechanical equilibrium which facilitates a large cascade of motion. Combined with previously identified phenomenology, we have a complete kinetic pathway for structural change which is common to both ageing and crystallization. Furthermore, this suggests that transient force balance is what distinguishes glasses from supercooled liquids.

A geing and devitrification are very slow dynamical processes taking place in glasses. Ageing leads to a gradual change in the physical properties of glasses 1,2 ; a prominent example is a change in sample size, potentially very harmful for applications to high precision devices. Devitrification, the transformation of a glass to a crystal, has also been a subject of intensive study for many years 3 . For material scientists, the tendency for a supercooled system to devitrify during ageing or on heating can be detrimental to the stability of organic 4,5 and metallic glasses 6 , silicates 7 , macromolecules 8,9 and aqueous systems 10 . The latter is particularly relevant to cryogenics, where ice formation can have adverse effects on biological samples 11 . For many industrial products, including pharmaceuticals 4 , devitrification during storage is also a very serious problem. Despite the technological importance of these phenomena, their fundamental mechanisms remain elusive, because of their slow dynamics and apparently stochastic nature.
To access the microscopic mechanisms behind ageing and devitrification, hard-sphere (colloidal) glasses provide an ideal system due to the experimental accessibility to particle-level information [12][13][14][15][16] and the simple nature of the hard-core interaction. Different studies have considered ageing and devitrification separately in hard-sphere glasses. In the case of ageing, an increase in structural order over time was found 17 . For hard spheres (both monodisperse and weakly polydisperse), structural order is characterized not only by density but also by bond orientational order (BOO), a measure of angular order between neighbouring particles. This ordering originates from a thermodynamic driving force to lower the free energy of the system. Note that higher BOO means larger vibrational (or correlational) entropy, that is, lower free energy for hard spheres (see, for example, ref. 18). Regarding devitrification, on the other hand, Pusey et al. 19 showed that crystallization can take place beyond the glass transition point in over-compressed hard-sphere glasses. This may apply more generally, as jammed or glassy metastable states are also seen to undergo ordering.
On deeper overcompression, Sanz et al. 20 found that crystallization occurs through discrete collective events, where groups of particles suddenly undergo large displacements, accompanied by an equally sudden increase in the proportion of crystalline particles in the system: these events were termed 'avalanches' for their intermittent, collective nature. They first found that the randomization of particle velocities in their molecular dynamics simulation averted the incidence of the events, and second, that polydispersity suppressed the growth of crystallinity while maintaining the same intermittent dynamics. Thus, they concluded that avalanches 'mediate' crystallization, or more precisely, that 'chance' collective motion in the system seemed to trigger the displacements that led to crystallization. At the same time, they noted that there was a spatial heterogeneity to which such events are linked, and suggested 'soft spots' 21 as a potential candidate for what might distinguish them from other regions. So while ageing is proposed to be thermodynamically driven 17 , it is argued that devitrification and avalanche events are of a strongly stochastic nature with an unidentified structural signature for initiation; crystallization is their by-product. Thus, the connection between ageing and devitrification, if any, remains elusive: despite being the two principal dynamic phenomena taking place in glasses, their physical mechanisms also remain unknown.
For example, what makes the kinetics of the phenomena in glasses special remains unclear, since the relation of ageing and devitrification to corresponding processes in the supercooled liquid state, that is, structural relaxation and crystallization, respectively, are yet to be ascertained.
Here, we address these fundamental problems by studying structural evolution during avalanches in both ageing and devitrifying hard-sphere systems. We show that ageing and devitrification are both characterized by the same kinetic processes, and that the special feature of dynamics in glasses originates from a temporal mechanical balance characteristic of solids, which has not been clearly recognized so far. We distinguish the very first particles to take part in avalanche events (AIs (avalanche initiator)) and go on to show that they are characterized by a structural signature, with a lower local density and low BOO. We then show that this motion leads to a transient loss in the force balance of the system, and it is this that gives rise to the cascade of particle motion which we perceive as an 'avalanche'. Finally, we find that this motion triggers displacements throughout the system in the vicinity of pre-ordered regions. We believe this to correspond to the avalanche devitrification identified by Sanz et al., though we go on to show that this happens in ageing systems as well.

Results
Progression of an avalanche. We simulate supercooled monodisperse and weakly polydisperse hard-sphere-like systems in the glassy regime to study devitrification and ageing, respectively. Details regarding the system and the simulation are given in the Brownian dynamics simulation section in the Methods. For the monodisperse case, we see a clear increase in the crystallinity, but not for the polydisperse case, as shown in Supplementary Fig. 1. We proceed to study how an avalanche event is triggered in a glass state, how this triggering event spreads over the system, or induces an avalanche, and how structural ordering is finally enhanced by the avalanche. We first observe the incidence and progression of avalanches, looking at the mean squared displacement (MSD) of particles in the system from some arbitrary time t 0 . An example is shown in Fig. 1a for an ageing system. We can see that the sudden collective displacement of particles takes place intermittently. We proceed to look at the circled event in more detail in Fig. 1b. First, we identify four different times during the event: a time just preceding the event, t i ; the time at which the first particles are displaced and the MSD suddenly rises, t a ; an intermediate time t m , purely for illustrative purposes; the time at which the MSD finally reaches its next plateau, t f . Snapshots of the system at each of these times are given in Fig. 1c-f, with blue shaded regions corresponding to regions with high BOO parameter Q 6 40.25 (see the Structural analysis section in the Methods for the definition of BOO) at time t i . Only particles which are displaced more than s l /3 (s l being the size parameter of particle l (see the Analysis of avalanche events section in the Methods) between t i and the time of the snapshot are shown.
Starting from a quiescent, metastable state in Fig. 1c, we see that a small, localized cluster of these particles appear in panel d. When the number of nearby particles displaced more than the s l /3 threshold becomes 45 for the first time after time t i , we call these particles AIs. These will be discussed at length below. Note that they appear in disordered regions with low local crystalline BOO (low Q 6 ). This is followed by particles in the vicinity being subsequently displaced, as shown in panel e, before the system reaches the next plateau at panel f. The avalanche particles identified at the end of the event are called APs, and, with few exceptions, are significantly more numerous than AIs.
Structural origin and evolution of AIs and APs. Note that we now have three populations of particles, AIs, APs and the whole set. These can be analysed separately with a large enough sample of events. With lists of AIs and APs over all the events, we go on to look at the typical local structure of AIs and APs at t i . We start with local volume density f l , defined as 1 6 pd 3 l =V l vor , where d l is the effective particle diameter (see the Brownian dynamics simulation section in the Methods) and V l vor is the local volume of particle l as found from a radical Voronoi tesselation. A distribution of f l values for AIs, APs and all particles at t i , is shown in Fig. 2a,d for devitrifying and ageing systems, respectively. It is clear that AIs are found in regions with a markedly lower density than the global average. APs exhibit this to a much lesser extent.
A lower density can arise from either larger local Voronoi volumes or simply smaller particle size. In the monodisperse (devitrifying) case, only the former is possible, whereas both are possible for the polydisperse case. We look at the distribution of local Voronoi volumes for the polydisperse (ageing) case, only to find that the volume for AIs is in fact skewed towards lower values, that is, towards higher volume fraction, albeit only slightly. This is shown in Supplementary Fig. 2a, where the distributions of local Voronoi volumes are given for AIs, APs and all particles. The lower local volume fraction in the polydisperse case must thus come primarily from AIs having smaller particle size, as confirmed in the particle size distributions given in Supplementary Fig. 2b.
Similar intermittent dynamics have been reported for ageing in glasses in simulations and experiments [22][23][24][25] . These works identified avalanche events as barrier-crossing events in a complex energy landscape under the influence of internal stresses. In this energy-landscape picture, AIs and APs can be related to the activation stage (from the initial states to the nearby saddle states) and the relaxation stage (from the saddle states to the final states), respectively. Egami et al. 25 also found cascade events spreading over the whole system in unstable glasses formed by an instant quench, which may be caused by largeamplitude density inhomogeneities frozen in by such a quench. A link between such nonlocal events and self-organized criticality was also suggested 23,25 , but it is beyond the scope of our paper. Here we focus on localized events in a rather stable glass, which is crucial for identifying structural characteristics that are independent of the box size. This allows us to reveal a specific mechanism for event initiation, propagation and consequence.
Going on to look at BOO parameters Q 6 and w 6 , we superimpose the values of these parameters for AIs and APs over all trajectories for devitrifying and aging cases over a typical distribution over the whole system at some arbitrary time point (Fig. 2b,e). Q 6 is associated with crystalline BOO in hard-sphere systems, while a lower (more negative) w 6 is associated with dodecahedral and icosahedral BOO 26 . It is clear that AIs can be found concentrated in regions where neither Q 6 nor À w 6 are notably large, that is, in disordered regions. The average value over the AIs is noted as a white square on each, both below the thresholds associated with crystals and locally favoured structures. The distribution of APs, on the other hand, are indistinguishable from the typical distribution.
From these two findings, we can conclude that avalanches are more likely to initiate from a low degree of structural order and a low local Voronoi density (see Fig. 3 on the correlation between the structural order and local density). This region may be regarded as a kind of 'defect' in a glass.
We look further, and see whether certain changes in translational order or BOO are more likely for AIs and APs compared to the population as a whole. This can be shown by calculating p(DO) subset À p(DO) whole , where O is the order parameter of interest. Looking at the devitrifying case first in Fig. 2c, we find several things. First, it is clear that both AIs and APs, with their large displacements, generally entail a greater change in the local volume fraction f l than average, giving the curves their distinctive 'cursive V' shape. AIs and APs are less likely to have a local volume fraction change |Df l |o0.04 than average, and more likely to have a f l change 40.04, with a maximum around Df l j j$0:06. This implies that particles with large amplitude displacements are most likely accompanied by a Second, we note that the AIs are very clearly skewed towards an increase in local volume fraction f l . What is particularly interesting is that this feature is not observed for APs, that is, they are equally likely to find a positive and negative local volume change. We will show later that this is a direct consequence of the fact that the motion of APs is of mechanical (not thermodynamic) origin. Looking to the ageing case in Fig. 2f, we once again see this asymmetry for AIs and APs. The statistical robustness of these results is discussed in the Event statistics section in the Methods. It is worth noting that local-structural order and local density are not entirely independent of each other at such a high packing density 18 . We are clearly in a regime where translational ordering and bond orientational ordering are closely related to each other, particularly for hard-sphere-like systems, where the local symmetry is selected by packing. Instead of looking at the structure at t i , we now look at the change in local volume fraction f l and BOO over the time interval t i otot f for the three subsets, starting with AIs. To sharpen the sensitivity of the BOO parameter to localized change, we choose q 6 instead of its coarse-grained counterpart Q 6 . As seen in Fig. 3a,b, we see a positive correlation between an increase in q 6 and f l . Note that the colour map simply reflects the number of particles in the plot for particular Df and Dq 6 ranges. It should be noted that the trend is much stronger for the devitrifying case.
As seen in Fig. 2, AIs originate from regions of lower local volume fraction (or, density) and low BOO. Thus, the increase in f l is not so much a densification beyond the system average, but rather a relaxation of a density inhomogeneity present in the heterogeneous, over-compressed structure. The causal relationship between them is clearly highlighted by the difference in the degree of asymmetry seen between AIs and APs. Since AIs have a structural origin, their states before and after are distinct. There is thus a clear structural distinction between AIs and APs: the asymmetry in structural evolution of AIs indicates the presence of a local, thermodynamic driving force for the initiation of avalanche events. On the other hand, it would appear that APs are rather a kinetic by-product of the AIs (discussion continued below). Here, it is worth noting that it is known 27 that there is no discernible difference between the distributions of Voronoi volumes and mobility for a supercooled liquid: it is our ability to identify AIs and APs that uncovers this structural distinction.
These conclusions are further confirmed for the monodisperse case by using an iso-configurational ensemble. We produced an iso-configurational ensemble (B2,000t B length), initiated from a configuration taken from one of the many MSD plateaux in the monodisperse case. A total of 100 trajectories were found, and 38 events were recorded. Further to reproducing the spatial heterogeneity identified in ref. 20, we apply the same treatment as above to produce Fig. 4a-c: the trends are exactly the same as what was found for 50 independent trajectories. As seen in Fig. 4d, there is good agreement between the unique AIs (yellow) and regions which obey all of the following criteria: Q 6 o0.25, w 6 4 À 0.01, f l o0.62, where the AI f l probability distribution becomes larger than the distribution over the whole system. Low density and low BOO thus provide us with the regions where AIs are more likely to appear, as they are characterized by lower activation energies for particle motion, where the caging effect is locally weaker. As in any stochastic process, the study of regions show the threshold values separating particles with highly crystal-like order (Q 6 ZQ X 6 ) and those with high dodecagonal or icosahedral-like order (w 6 rw dod 6 ) 26 . Note that crystals have Q 6 larger than 0.4 (see the Structural analysis section in the Methods for details). The white thick open square represents the average values of w 6 and Q 6 for AIs. On average, thus, AIs possess neither significant crystalline order nor locally favoured structures. This indicates that AIs are particles localized in disordered regions. (c,f) Difference between the probabilities of Df l for AI or AP subsets and the distribution for the system as a whole. There is a bias towards larger change for both AIs and APs, but there is a skew to positive change for AIs which is not seen for APs. For both devitrification and ageing, APs almost equally gain or lose density, but AIs have a strong tendency to increase local density.
of low density and low BOO order only provides a probability map of where AIs will occur. The actual location cannot be predicted in a deterministic way, as the underlying motion is Brownian, without memory effects (inertia).
Function of APs in structural evolution. Does this mean that APs do not have any role to play in crystallization or ageing? If the only significant structural evolution in the system was derived from the AIs, the total change would be minuscule, as AIs associated with an event make up only around 0.2-0.3% of the system (B10 particles in 4,000). Consider the events shown in This last point is worth considering carefully. Though there is some overlap, the particles which in fact experience the largest change are those which move less than the avalanche threshold, but more than average. Figure 5b,d show the fraction of particles with a particular displacement size for particles undergoing a large change in Q 6 , for monodisperse and polydisperse cases, respectively, using the same thresholds for Q 6 changes as in panels a and c. The distribution of displacements for particles which undergo large structural change are peaked below the avalanche threshold, Dx4s l /3. This is not simply a reflection of the overall distribution of displacements in the system (see the red dashed curve), which has no peak and monotonically decreases with an increase in the displacement. Also note that there is a clear bias towards positive changes (orange curve) in Q 6 compared to negative changes (blue curve) for the devitrifying case.
In the above we have shown that APs are not the particles which crystallize, as also reported in ref. 20. In the final stage of the avalanche identified in our work, APs cause movements which are smaller than the avalanche threshold and yet induce significant BOO development. This raises the fundamental question of what is the physical mechanism behind the avalanche mediated devitrification 20 in monodisperse systems, or why such small particle displacements (less than the particle size) can facilitate crystallization. The natural answer to this question arises from a crystallization pathway where the first step is development of spatial coherency in BOO and the fact that small displacements are enough for this enhancement 28 . On noting that enhancement of BOO spatial coherency cannot induce translational ordering in systems with a high enough polydispersity 18 , we can now also claim the same mechanism for ageing systems.
In supercooled liquid and glassy states with weak frustration against crystallization, like the one studied here, the system always seeks a way to lower the free energy even at a local scale by using every possible kinetic path. In a glassy state, a system tries to increase structural order locally, which results in density inhomogeneities accumulated in the disordered region overcoming the energy barrier to motion (AIs), causing a cascade of particle motion (APs). important question arises: what are the conditions underlying the intermittency of the dynamics? In an ordinary supercooled liquid, particle motion proceeds in a continuous manner. It is only in a highly supercooled or glassy state that we see intermittent, avalanche-like particle motion. We argue that the latter is a characteristic of a glassy state, where mechanical equilibrium is satisfied both locally and globally, at least transiently. Suggestions have been made 29,30 and evidence given 31,32 that mechanical balance is attained in a glassy state or in the inherent structure, resulting in the emergence of a long-range elastic correlation like in jammed granular matter 33 . Before an avalanche, a system is in a mechanically stable state even with thermal noise: the local ordering events which trigger the motion of AIs lead to the breakdown of this transient mechanical equilibrium, or the loss of force balance, until the system attains a new mechanically stable state.
Evidence for this in our system is given in Fig. 6, where we look at how particles change their neighbours over time. We can define two types of neighbour, simply using proximity (nearest neighbours, NNs) or those with which a repulsive force is acting (force neighbours, FNs). The average number of neighbours kept compared to some initial time is shown in Fig. 6a, both for FNs and NNs, superimposed on the progression of an avalanche event in a polydisperse trajectory. We can immediately see that the system is in a state of mechanical equilibrium before, with an extended plateau over which particles keep its neighbours, and maintains the force chain structure of the system. It is only when the avalanche initiates that this structure changes. The existence of a clear plateau afterwards shows that mechanical equilibrium has been restored, with a different connectivity between particles as before.
The best evidence for the role of mechanical stability is spatial correlation between AIs/APs and new FN connections. This is shown in Fig. 6c, where new force connections formed over the course of an avalanche have been shown as bold red bonds, superimposed on transparent blue bonds representing the original force network before the event. The AIs (yellow) and APs (red) for the event are also shown. The abrupt formation of these new FNs is also shown dynamically in the Supplementary Movie.
It is worth noting that the average number of FNs surrounding AIs and APs is around 6, the Maxwell rigidity criterion 34 in three dimensions. Distributions for the number of FNs for AIs, APs and all particles for both monodisperse and polydisperse cases are given in Supplementary Fig. 3. The role of this in the dynamics becomes clearer when different volume fractions are considered. Short trajectories are produced for volume fractions f ranging from 0.58 to 0.66 using the same number of particles n, and the average number of FNs recorded after tE300t B . It can be seen that the average number of FNs, n FN , first passes this threshold when the volume fraction reaches f ¼ 0.65 (Fig. 6b). This coincides with the volume fraction above which the avalanches can be observed, that is, sudden jumps are seen in the MSD preceded and followed by clear plateaus (see Supplementary  Fig. 4). Also note that the increase in n FN starts from the nominal glass transition volume fraction, f g E0.58-0.59 (Fig. 6b), corroborating the role that mechanical contacts might have to play in glass formation.
Note that the frequency of these jumps also changes with volume fraction. Avalanche events take place in an intermittent manner when a system is able to overcome the barrier between one metastable basin and another with the help of thermal noise. Since the height of this barrier steeply increases with an increase in the degree of supercooling, the frequency of avalanches should decrease for deeper supercooling; we confirm this in our simulations, as shown in Supplementary Fig. 4, where avalanche incidence is increasingly suppressed at higher volume fractions. This trend is also seen qualitatively in ref. 35 using a much sharper inter-particle potential.

Discussion
We have now seen that avalanche initiation has a structural precursor, and that the propagation of APs and the intermittency of the dynamics is related to the temporary loss of mechanical stability of the system for both devitrifying and ageing systems. Avalanches further induce small-amplitude motion of particles (os l /3) in regions which already have a high structural order. It is however worth noting what makes the devitrifying and aging systems different. In the devitrifying systems considered here, the crystal growth is interface-limited, meaning that the barrier for crystallization comes from the addition of particles to the crystalline surface. This process does not require large-scale translational motion 28,36,37 , and can thus occur under the smallamplitude motion set in place by AP particles. In most ageing systems (mixtures or highly polydisperse samples), crystal growth is instead diffusion-limited, thus requiring large-scale diffusional motion for the interface to grow, and it is thus not easily activated by avalanche motion. This is shown by the fact that particles undergoing large changes in BOO for the polydisperse case do not seem to reflect the clear bias to positive changes in q 6 seen in the monodisperse case, as shown in Fig. 5b,d. This is also evident for AIs and APs in Supplementary Fig. 5a,b, where there is no clear bias to order development like in Fig. 2c,f. These show that BOO development in ageing systems is not in the expected thermodynamic direction for every event. Of the 50 polydisperse runs, 32 cases experience an overall increase in Q 6 , the others experience a decrease. There is, however, an overall positive trend: the mean change in Q 6 per event is DQ 6 ¼ (2.24 ± 1.94) Â 10 À 4 , where the error is an unbiased standard error in D Q 6 . This corresponds to E87% confidence in a Q 6 increase. This is in contrast to a DQ 6 ¼ (3.83 ± 1.39) Â 10 À 4 for the monodisperse case, with a 99.8% confidence in a Q 6 increase. This bias to positive changes reflects the thermodynamic direction, which is hindered by the suppression of long-range diffusion in the polydisperse case.
As seen in Fig. 5c, regions where large changes in BOO occur are localized to pre-ordered regions (MRCO). Previous work on ageing in 3D hard-sphere systems at lower volume fractions 17 has shown that there is a slow development in the spatial correlation length of BOO parameter Q 6 . Such an increase is most easily achieved if re-structuring takes place near pre-ordered sites, explaining the large changes in Q 6 close to MRCO regions induced by nearby AP motion. This is analogous to the localization of Q 6 growth to the MRCO in the devitrifying case, with an added stochasticity factor.
Finally, we give evidence for the fact that APs form strings, as briefly noted in ref. 20. Supplementary Fig. 6 shows the degree of overlap with the former position of neighbours at an arbitrarily chosen initial point over time, a decrease in which denotes stringlike displacement. The exact definition is given in the String-like collective motion section in the Methods. The event is the same as the one used in Fig. 6a; the coincidence with avalanche incidence is clear. It is particularly interesting that the collective, string-like motion seen in the Adam-Gibbs scenario for quasi-equilibrium supercooled liquids 38 could be involved in the structural evolution of the glassy state of unstable glasses like the ones shown here. String-like motion may be the only possible way for particles to move in a highly packed state 38 . This similarity in dynamics lends further credence to a thermodynamic picture of the glass transition and the intrinsic link between crystallization and vitrification under the influence of frustration 18 .
The emerging picture of an avalanche in glasses is thus the following. Avalanches occur between mechanically stable states, where the force network between particles does not change with time. They are initiated by small clusters of particles (AIs) that have lower than average local density, within environments with a low crystalline BOO, a lack of locally favoured structures (icosahedra), and where the number of constraints is close to the threshold for mechanical stability. In polydisperse systems, AIs are more often small particles. Due to the loss of mechanical stability, the displacement of the AIs is propagated through the system in a cascade-like fashion; these particles are called APs and, on average, do not directly experience distinct structural evolution, that is, locally, they will almost equally likely decrease or increase the amount of crystalline order in the system, although slightly biased towards the increase. This makes them unique from AIs. However, the APs induce significant change in their surroundings which move less, but cause a development of Q 6 in the vicinity of nearby highly ordered clusters. This triggering of ordering is analogous to the crystallization of supercooled liquids under some external perturbation. In short, we find that the ageing and devitrification of glasses proceed via a common pathway involving the interplay of thermodynamics and mechanics: thermodynamically-initiated mechanical facilitation of particle motion (that is, avalanche) and the enhancement of structural order by the cascade of small-scale displacements induced by this avalanche event.
This work highlights a kinetic pathway via which strongly supercooled glassy systems can age or devitrify, highlighting the particular role that inhomogeneities in density and structural order have to play in their stability-a stronger frustration (c) Re-arrangement of the force network during the avalanche event shown in a. Yellow and red particles are AIs and APs respectively, as in Fig. 5a,c. The stable force network that is not involved in the re-arrangement is shown by the blue network, whereas newly formed force bonds resulting from a re-arrangement since a time preceding the avalanche are shown as thick red bonds. We can see some correlation between APs and AIs and a significant local re-arrangement of FNs. This is a snapshot after an avalanche event; please see the Supplementary Movie to see the correlation over the course of the whole event.
towards this kind of relaxation would make a glassy state more stable against ageing or devitrification. This is a useful guiding principle for the stabilization of the amorphous state in various materials. For a multi-component glass, for example, local demixing is required for devitrification. The mechanism of devitrification in such a mixture remains for future investigation.

Methods
Brownian dynamics simulation. We simulated supercooled monodisperse and weakly polydisperse glasses in the glassy regime using Brownian dynamics simulation for particles interacting through a Weeks-Chandler-Andersen (WCA) potential 39 in a three dimensional box with periodic boundary conditions. We employed the monodisperse system as an example of a devitrifying system and a weakly polydisperse system (6% Gaussian) as an example of an ageing system 40 .
The difference between them is clearly shown in Supplementary Fig. 1, where the proportion of crystalline particles p X (see below) is shown over time. The devitrifying system experiences growth in crystallinity, whereas the ageing system does not. Initial over-compressed states were attained using rapid compression, expanding the size of the particles from 1% their intended size to their final dimension in less than a Brownian time t B , including noise to allow frustration limited structures to relax. This is similar to the Lubachevsky-Stillinger algorithm 41 . Although this volume fraction exceeds random close packing fraction (B64% and B66% for monodisperse and polydisperse cases, respectively 18 ), local-structural ordering such as BOO allows a system to retain mobility. We note that the volume fraction is well above the glass transition volume fraction (f g B0.58-0.59) and structural relaxation in the system, including the time scale of avalanches, would be significantly longer than the Brownian time t B . As our primary interest is in persistent changes in structure as opposed to transient fluctuations, we average all positions and BOO parameters over 3t B intervals.
A Brownian Dynamics simulation code was developed in FORTRAN 95 using the Ermak-McCammon algorithm 42 , with no hydrodynamic interactions and no inertial term. Particles interact via a WCA potential, U WCA ðrÞ¼ 4E s=r ð Þ 12 À s=r ð Þ 6 þ 1=4 È É for ro2 1 6 s, where s is the diameter of the particle, and the temperature kept such that k B T=E ¼ 0.025 (k B being the Boltzmann constant), as in refs 17,43. Volume fractions are calculated using an effective diameter d ¼ 1.0953s where the WCA potential decays to k B T. This is in agreement with a mapping of the WCA potential phase diagram to the hard sphere one to within o0.5% (ref. 44). The diffusivity D is chosen to mimic experimental colloidal systems, assuming a 1 mm particle radius in aqueous viscous surroundings with a viscosity of 1 mPas.
For all simulations shown in the main text, the number of particles is n ¼ 4,000 unless otherwise stated. Here we mention the dependence of the observed dynamics, particularly the MSD, on system size. We produced eight trajectories with a system size n ¼ 16,000, and compared their MSDs with eight randomly chosen trajectories from the n ¼ 4,000 set in the monodisperse case. This is given in Supplementary Fig. 7. It can be seen that events are smaller and less well-defined for larger system size. This is to be expected, since the properties of individual events should be independent of system size, that is, since the MSD is normalized by particle number, and the events give smaller jumps in the MSD. Also, a larger system size makes it more likely for multiple events to be occurring at once, hence the lack of clear, flat plateaux, as events initiate before APs from another event have relaxed.
Structural analysis. The metrics we employ include local density as characterized by a radical Voronoi construction, and the q 6 order parameter which has found wide acceptance as a suitable order parameter to detect crystalline order in hard spheres 45,46 . Five-fold symmetric structures, such as icosahedral and dodecahedral structures, are targeted with the w 6 order parameter 26 . BOO parameters such as q 6 and w 6 are calculated as described in ref. 47. To reduce fluctuations in the distribution of q 6 , a coarse-grained version, commonly written as Q 6 , is also employed 48 . All the definitions of the order parameters can be found in ref. 28. Here we note that crystals have Q 6 larger than 0.4 (refs 45,49). The proportion of crystalline particles p X is defined as the proportion of particles in the system which have seven or more 'solid' bonds, that is, bonds between NNs for which the bond coherence d 6ij ¼ q 6 ðiÞ Á q 6 ðjÞ= q 6 ðiÞ j jq 6 ðjÞ j j(ref. 28) exceeds 0.7. Radical Voronoi constructions were carried out using the VORO þ þ package available at http://math.lbl.gov/voro þ þ / (ref. 50).
Analysis of avalanche events. We analyse three populations of particles, AIs, APs and the whole set, with a large enough sample of events as follows. In 32 independent simulations of devitrifying and ageing systems, each B7,000t B in length, we identified 50 avalanche events each for both cases. Care was taken to identify events which had a clear beginning and end, that is, with a clear plateau in the mean-squared displacement at both ends. Each event was labelled manually with a t i and t f , and APs were found by comparing the state of the system at these two times; AIs for each event were found by comparing particle positions for progressively later times from t i until the condition above was met.
Event statistics. We briefly comment on the statistics of this result. Histograms were generated by taking probability distributions over AIs, APs and all particles for individual events and then averaging these, avoiding bias of data to larger events among the 50 that were identified. The marker size (in both Fig. 2c,f and Supplementary Fig. 5a,b) is an indication of the maximum contribution possible from each event. In total, 315 AIs and 6,121 APs were counted for the monodisperse case, 320 AIs and 3,718 APs for the polydisperse case. Given the relatively small number of AIs, the width of bins could not be made significantly narrower than what is shown. To be sure that these trends are not due to some particularly favourable events mixed into the set, we focused on the two peaks in Fig. 2c,f, took a random selection of 25 events amongst the sets of 50, and measured the asymmetry between the AI and AP curves, marked D AI and D AP in the figure. This was repeated 100 times. In the monodisperse case, D AI ¼ 0.0976 ± 0.0235 and D AP ¼ 0.0384±0.0148; in the polydisperse case, D AI ¼ 0.1044±0.025 and D AP ¼ 0.0099 ± 0.0108. In both cases, D AI 4D AP .
String-like collective motion. Following the formalism of Donati et al. 38 , we calculate the value of the collective motion parameter d¼s À 1 ij min X i ðtÞ À X j ðt 0 Þ ; À X j ðtÞ À X i ðt 0 Þ Þ, averaged over NNs i and j, and then over all particles in the system over time. s ij is the sum of the radii of particles i and j. Here X i is the position vector of particle i. The original work set a threshold for d below which strings were considered to move collectively. Thus, we can interpret a decrease in the global average value of d to be a collective displacement. As shown in Supplementary Fig. 6 (same polydisperse trajectory as in Fig. 6a), the global value of d decreases during the course of an avalanche, with well-defined metastable states before and after.s Data availability. The raw data analysed to derive the findings of this paper is available from the corresponding author upon reasonable request.