The structural origin of the hard-sphere glass transition in granular packing

Glass transition is accompanied by a rapid growth of the structural relaxation time and a concomitant decrease of configurational entropy. It remains unclear whether the transition has a thermodynamic origin, and whether the dynamic arrest is associated with the growth of a certain static order. Using granular packing as a model hard-sphere glass, we show the glass transition as a thermodynamic phase transition with a ‘hidden' polytetrahedral order. This polytetrahedral order is spatially correlated with the slow dynamics. It is geometrically frustrated and has a peculiar fractal dimension. Additionally, as the packing fraction increases, its growth follows an entropy-driven nucleation process, similar to that of the random first-order transition theory. Our study essentially identifies a long-sought-after structural glass order in hard-sphere glasses.

W hen a liquid is cooled towards the glass transition, the dynamics slows down dramatically. The mechanism of this phenomenon has been extensively investigated for decades, but there is no consensus on it yet 1,2 . Recently, the discovery of dynamic heterogeneity and the rapid increase of its correlation length near the glass transition suggest the collective nature of the dynamics, which has inspired hope that a corresponding static correlation length associated with some critical behaviour similar to an ordinary phase transition can be identified 2,3 . This length scale is also supposed to have a configurational entropy origin, as originally suggested by the Adam-Gibbs theory 4 . Recent searches for this static correlation length have identified crystalline, icosahedron-like, or 'point-toset' type of orders which show significant increases of the correlation lengths as dynamic arrest is approached 5 . However, for all existing approaches, either the relationship with the configurational entropy or the structural nature of the order remains to be understood.
The hard-sphere system is a popular model glass former because it can simulate systems such as metallic glasses, granular systems and colloidal suspensions 6 . The close relationship between static granular packing and hard-sphere glasses dates back to the pioneering work of Bernal 7 , who experimentally investigated granular packing to model liquid structures. Additionally, agitated granular systems exhibit slow dynamics that resemble those of thermal glassy systems [8][9][10][11] . The analogy between granular packing and a thermal glass can be formally understood by Edwards' ensemble, on which a statistical framework similar to the equilibrium statistical mechanics can be established for static granular packing 12 . Additionally, the hard-sphere system has a zero-temperature geometric phase transition, the jamming transition, which demonstrates interesting properties like marginal stability and critical scaling behaviours 13,14 . Jamming transition and its distinction with the glass transition have recently been characterized from both thermodynamic 15 and rheological 16 points of view. A recent theoretical approach has tried to incorporate the jamming transition into the framework of glass theory through a Gardner transition to fractal sub-basins in the free energy landscape 17 . However, the possible structural changes associated with these two transitions and the structural nature of the fractal sub-basins remain to be explored.
Using granular packing as a model hard-sphere glass system, we identify a geometrically frustrated polytetrahedral structural order which fulfils the requirements of a static glass order. By using synchrotron X-ray imaging techniques 18,19 , we carry out a systematic study of the dynamics, thermodynamics, and structures of tapped granular packing. We demonstrate that the system exhibits all key phenomena of a thermal glassy system. Particularly, a polytetrahedral structural order grows rapidly as the packing fraction increases and it is spatially correlated with the slow relaxation dynamics. The non-trivial fractal dimension and length scale of this polytetrahedral order are consistent with an entropy-driven nucleation model similar to the random first-order transition (RFOT) theory 20 .

Results
Thermodynamic variables based on Edwards' ensemble. To obtain the thermodynamic variables and detailed local structures of the packing, we carry out X-ray microtomography scans on packing of different average packing fractions F prepared by three different protocols (Methods). Based on the statistical framework originally proposed by Edwards and coworkers 12  thermodynamic variables, such as entropy S and compactivity w (similar role as temperature), are evaluated from the fluctuations of the reduced Voronoi cell volumes v voro ¼ V voro /V g , where V g and V voro are the particle volumes and their Voronoi cell volumes respectively, according to 21,22 where The constant, which is similar to the Boltzmann constant, is set to unity. V g is also set to unity for convenience. We further assume that the compactivity of random loose packing is infinitely large, that is, 1 wðF¼0:55Þ ¼ 0, to complete the integral of equation (1). In equation (2), the entropy S(F) can only be identified up to a constant and we assume S(F ¼ 0.64)E1.1 by using the Shannon entropy calculation results 22 . The integral of equation (1) is calculated numerically starting from a second-order polynomial fitting of s 2 v as a function of F: s 2 v ¼ 2:387F 2 À 3:063F þ 0:997 (Fig. 1a). Our results are in quantitative agreements with previous study 22 (Fig. 1b). It is worth noting that the configurational entropy defined above corresponds to the complexity in hard-sphere glass terminology and should not be confused with the definition as the entropy difference between liquid and crystal phases 6 .
Like a thermal hard-sphere glass, S decreases with decreasing w or increasing F. We extrapolate S(F) and S(w) curves to S ¼ 0 to obtain the corresponding glass-close-packing (GCP) packing fraction F GCP ¼ 0.671 or GCP compactivity w GCP ¼ 0.0554 (refs 6,23,24). F GCP is close to the jammed ideal glass transition density 0.68 from replica theory calculations 6 . In the following, structural relaxation time and structural correlation length will be expressed as functions of these thermodynamic variables to draw analogies between the tapped granular system and a thermal hard-sphere glass.
Relaxation time. To study the slow relaxation dynamics in granular packing, we use X-ray absorption imaging to measure the time evolution of the average packing fraction F under tapping (Methods). The structural relaxation times t are calculated from the compaction curves at different tapping intensities G. The packing is first tapped at G ¼ 15 for 1,000 times to reach F 0 ¼ 0.615. Then, the compaction curves are measured by tapping the packing at different G to reach the corresponding reversible-branch packing fractions F N ¼ F(G) (Fig. 1c). The packing that has reached reversible-branch is at steady state and memoryless. Each compaction curve can then be fitted using the empirical Kohlrausch-Williams-Watts law 8 where t is the number of taps and b is the stretching exponent (Fig. 1d). The fitted values of b lie in the range of 0.5-0.9 if both t and b are allowed to vary. To be consistent, we fix b ¼ 0.7 to obtain t for different compaction curves. Nonetheless, the t values remain essentially unchanged as compared with the case when b is allowed to vary. As shown in Fig. 2, t increases with decreasing w or increasing F, similar to a thermal hard-sphere glass 25 .
Polytetrahedral order. As shown above, the great analogies in slow dynamics and thermodynamics between tapped granular packing and a thermal hard-sphere glass suggest a common origin. In the following, we propose that the formation of local geometrically frustrated quasi-regular tetrahedra is the microscopic mechanism for the dynamic arrest in both systems 26,27 .
We demonstrate that the polytetrahedral order associated with these quasi-regular tetrahedra corresponds to the long-soughtafter glass order in hard-sphere glasses 26,27 , by showing: the polytetrahedral order is spatially correlated with the slow relaxation dynamics; the static correlation length of this polytetrahedral order increases rapidly as F increases; and the size and shape of the polytetrahedral order are consistent with a configurational entropy-driven nucleation model similar to the RFOT theory 20 .
Similar to previous studies 28 , a quasi-regular tetrahedron is defined as a Delaunay simplex whose shape is close to a regular tetrahedron, with the shape deviation less than some threshold value of a polytetrahedral order parameter d ¼ e max À 1. In this expression, e max , in units of mean particle diameter s, is the length of the longest edge of the tetrahedron. The d values lie in a range between zero and an upper limit around one. There exist other measures to define a quasi-regular tetrahedron, such as tetrahedricity D (ref. 23). However, it turns out that our general results are not sensitive to any particular definition.
Correlation between order and dynamics. We first show that the quasi-regular tetrahedra are spatially correlated with the slow relaxation dynamics. In the tapping experiment with G ¼ 4, we track the trajectories of all particles for a number of taps by conducting a tomographic scan after each tap (Methods). We define Dr i ¼ dr i À v(r i ) as the diffusing displacement after one tap of particle i located at r i , where the absolute displacement dr i is subtracted by an averaged steady-state convection displacement v(r i ) (Methods). To characterize the mobility of particles belonging to tetrahedra of different d values, we define the d-mobility where the average is taken over all particles composing tetrahedra with d. As shown in Fig. 3a, a positive correlation between d-mobility and d value is clearly visible. The average value of ffiffiffiffiffiffiffi Dr 2 p after one tap is B0.04s, which is close to the typical cage size in the packing 30 . This is owing to the fact that one tap duration at G ¼ 4 is on similar timescale as the structural relaxation time t, that is, the particles have undergone many collisions upon one tap. It is worth noting that similar correlation between the particle mobility and the shape of tetrahedron it belongs to has also been observed in colloidal systems 31 , which suggests that the inherent friction of granular system is not the cause of the structure-dynamics correlation. ARTICLE As an alternative evidence, we define a tetrahedron correlation function p d (t) as the probability that one tetrahedron (Delaunay simplex) with d at t ¼ 0 is composed of the same four particles from t ¼ 0 to the t th tapping. This function captures the dependency of the local structural relaxation time upon d. As shown in Fig. 3b, tetrahedra with smaller d relax much more slowly than those with larger d.
Spatial correlation of polytetrahedral order. The fraction of quasi-regular tetrahedra grows as F increases, which is accompanied by increasing spatial correlations among them, that is, they tend to aggregate with each other. This spatial correlation can be demonstrated explicitly in terms of a percolation analysis on the Delaunay networks 32 . Tetrahedra are coloured according to their d values (tetrahedra with dod c are coloured where d c is a threshold) or coloured randomly but with the same number of tetrahedra. In both cases, face-adjacent coloured tetrahedra are joined together to form clusters. We cut a cubic region out of the packing and define L cluster as the longest spanning range of each cluster in directions parallel to the three axes of the cube. Max(L cluster ) grows as more cells can be coloured when the threshold value of d c is gradually relaxed, and it can ultimately reach the size of the cubic region L box at the percolation limit. As shown in Fig. 4, it turns out that tetrahedra chosen based on their d (or D) values show higher likelihood to percolate compared with the random-colouring case, which suggests their spatial correlations. To emphasize the unique relevancy of this polytetrahedral order, similar percolation analyses are carried out over various other structural order parameters, which show almost identical percolation behaviours with their random counterparts, suggesting that there exist negligible correlations 29 (Methods). Next, we use d c ¼ d* ¼ 0.245 to select out quasiregular tetrahedra at different F. The threshold value adopted is similar to previous studies 28,33 . The polytetrahedral order associated with these quasi-regular tetrahedra has a polytetrahedral structure 28 . The specific d* value chosen corresponds to a percolation transition of the polytetrahedral order at random close packing (RCP) (F ¼ 0.64) associated with the jamming transition of frictionless particles 13,28,34 .
We quantify how the spatial correlations among these quasiregular tetrahedra vary with F by calculating the correlation length x at different F. The average size of the polytetrahedral clusters x c is evaluated using the radius of gyration R g of all where N c is the number of particles belonging to a cluster, and R 2 g ¼ 1 2 h r i À r j 2 i ij is half the average square distance between all pairs of particles in a cluster 35 . We note that x c includes contributions from both the intrinsic correlation length x of the polytetrahedral order and a trivial correlation length x r associated with the random percolation process, which also increases with F (Fig. 5a). Therefore we define x ¼ x c À x r which equals to zero when the tetrahedra are uncorrelated. As shown in Fig. 5a, x, in units of s, also increases with decreasing w or increasing F, similar to t.
To further characterize the structural change associated with the growth of x, we study the evolution of P, which is the fraction of particles belonging to at least one quasi-regular tetrahedron. P increases from 80 to 98% as F increases from 0.572 to 0.634 (Fig. 5b), whereas x increases from 0.44 to 8.14 accordingly. This almost 20-fold growth of x can therefore only be induced by the merging of smaller polytetrahedral clusters into bigger ones instead of the simple inclusion of more particles. Meanwhile, the fraction of quasi-regular tetrahedra increases from 13 to 27%, as shown in Fig. 6a-d.
Dependencies of s and n on thermodynamic variables. In the following, we analyse the dependencies of t, x on thermodynamic variables w and S, similar to those in a thermal glassy system.
The relation between t and w can be well described by the Vogel-Fulcher-Tammann form (Fig. 2) The fitted w t ¼ 0.049 ± 0.003 agrees nicely with the 0.045 value from previous hard-sphere simulation 25 . The corresponding F t ¼ 0.678 is consistent with F GCP from above entropy extrapolation 1,6 , suggesting a similar relationship between dynamics and thermodynamics as a thermal hard-sphere glass 1 . t 0 ¼ 0.024±0.014 tap turns out to be the microscopic timescale in our system which is much smaller than either the tap duration or the structural relaxation time t.
Additionally, t and S follows an Adam-Gibbs type of relation where s c is the configurational entropy density, defined as S/hV voro i, with hV voro i being the average Voronoi cell volume (Fig. 2). The exponent yc d À y ¼ 1:0 AE 0:2 is surprisingly close to the Adam-Gibbs relation 4 . Here, the exponent yc d À y follows the convention of the RFOT theory 20,36 .   x also shows a diverging behaviour with decreasing w that can be fitted using a power-law function 37 which yields w x ¼ 0.051 ± 0.016 (F x ¼ 0.675 ± 0.023) and n ¼ 1.4±0.6 (Fig. 5a). The value of n is different from the critical exponent (E2/3) of the three-dimensional (3D) Ising universality class as suggested in recent experiments 38 . We further demonstrate that the growth of the polytetrahedral order is consistent with a configurational entropy-driven nucleation model similar to the RFOT theory 20 . In RFOT theory, the static correlation length or the mosaic size x mosaic , is determined by a competition between the gain in configurational entropy S and a free energy cost proportional to x y mosaic because of surface mismatch between different mosaics 20 . Specifically, the competition results in a typical length scale x mosaic / ð 1 Tsc Þ 1 d À y , where T is the temperature and d is the spatial dimension 37 . Motivated by this nucleation model, we plot x versus 1 wsc and fit the data according to (Fig. 5c). The fitting yields 1 d À y ¼ 1:78 AE 0:34 (y ¼ 2.44±0.11). At first sight, the y value obtained was incompatible with the original RFOT theory (y ¼ d/2) 37 or Adam-Gibbs relation (y ¼ 0) 4 . However, as shown in Fig. 6e, this discrepancy can be naturally reconciled with the fractal nature of our polytetrahedral order 39 , since it has a fractal surface dimension y s ¼ 2.57 (for F ¼ 0.634 packing) which is compatible with the extracted y from above nucleation model (Fig. 6e). Notably, this unusual y value has been observed before in both experiments and simulations 36 . Interestingly, y s decreases slightly with increasing F, suggesting that the polytetrahedral order will have less rough surfaces as F increases (Fig. 6f). These fractal polytetrahedral clusters fail to tile space because they are frustrated geometrically 27,40 . Furthermore, cluster size N c shows a power-law distribution: p N c ð Þ / N À m c . It turns out that the hyperscaling relationship d f (m À 1) ¼ d roughly holds, where d f is the cluster fractal dimension (Fig. 6e). Similar fractal clusters have been observed for immobile particle clusters in other glass-forming liquids 41 .
In addition to the proceeding cluster-analysis, we also attempt to extract the correlation length using a spatial correlation function of d (Fig. 5d). This function behaves rather similarly to the standard pair correlation function g(r) (ref. 19), which displays a finitelength decaying behaviour but no discernable differences for packing with different F. This indicates that protocols based on pair correlation analysis are incapable of capturing correlations in our system 34 . Same problem could also exist in the 'point-to-set' type of analysis, because pinning a smooth boundary might not be the best way to capture a fractal phase inside 42 .
where the average is taken for each pair of tetrahedra separated by distance r, and d i and d j are their d values. The location of a tetrahedron is defined as the average location of its four particles. The lines are shifted for clarity.

Discussion
The main conclusion of the current study is that quasi-regular tetrahedra are the structural elements of glass order in weakly polydisperse hard-sphere glass-particle systems. The order grows by following an entropy-driven nucleation process which is reminiscent of the diffusion limited cluster aggregation (DLCA) in kinetic gelation process, where independent clusters following DLCA growth processes touch in forming a global percolating fractal structure and acquire mechanical rigidity suddenly 43 . The fact that the growth is correlated can therefore induce cooperativity and non-Arrhenius behaviour in the system.
Similar to the gelation process in systems with attractive interactions 18,43 , we suggest that the jamming transition corresponds to a rigidity percolation transition of glass order for systems with repulsive interactions, that is, at RCP, the percolated polytetrahedral clusters acquire an infinite mechanical correlation length abruptly. The fractal shape of the percolated polytetrahedral order could therefore be related to the marginal solid behaviour 44 and unique scaling behaviours of the jamming transition 13 . Our jamming transition picture therefore suggests a unified scenario for rigidity transitions in systems with attractive or purely repulsive interactions: both are driven by the rigidity percolation of an underlying glass order.
The above percolation mechanism of jamming transition can also provide a simple geometric explanation of the continuous range of jamming density (the J-line) observed in numerical simulations 6,15,45 . In our system, the highest packing density we can theoretically achieve is B0.64, which corresponds to the RCP state. Interestingly, the RCP state has a non-zero entropy and consists of many polytetrahedral clusters, and therefore is not the ideal glass state 46 . At RCP, the average internal packing fraction of each cluster is approximately 0.67, which suggests that the rather low global packing fraction originates from the existence of cluster boundaries. As a result, if we can extrapolate the configurational entropy towards zero in obtaining a single large polytetrahedron spanning the whole system, that is, the jammed ideal glass state or GCP, then in principle, we can obtain a jammed packing density B0.67. Interestingly, this is exactly the upper limit of the J-line as has been predicted by the simulations 6,45 .
Additionally, since the packing in the current work is close to the lower-density-limit of the J-line, their equilibrium counterparts correspond to supercooled liquid states which are not very deep in the free energy landscape. This supercooled liquid picture is consistent with both the mean-field study by Mari et al. 45 , in which they found that the J-point correspond to the system just entering the landscape regime, and the fractal cooperatively rearranging regions (CRRs) found in glass-forming liquids near the dynamical crossover temperature 39 , where the CRRs bear great similarity to the fractal polytetrahedra in our system. The location in the free energy landscape also naturally explains our anomalous scaling exponent y as compared with that of the original RFOT theory, which mainly deals with mosaic states very deep in the free energy landscape.
In the current study, we implicitly assume the validity of Edwards' ensemble for packing prepared by both tapping and flow pulse protocols 21,47 . Despite the fact that packing prepared under these protocols has previously been established as ergodic, history-independent, and is therefore prone to a valid statistical analysis 21,47 , there still exist some ongoing debates regarding the basic assumptions of the Edwards' framework, especially the flat measure assumption [48][49][50][51][52][53] . As of today, a direct correspondence between the Edwards' entropy and the configurational entropy of a thermal hard-sphere glass is still waiting to be established theoretically. Therefore, despite the self-consistency of a thermodynamic analysis of our granular hard-sphere glass and its great analogy with a thermal hard-sphere glass, a direct one-to-one correspondence should not simply be taken for granted. This should also be born in mind when comparing our entropy-driven nucleation picture with RFOT theory.
Overall, our current study illustrates the origin of fragile glass behaviour in one type of model glass former. It suggests that other fragile glass systems can potentially be categorized by the growth of different types of structural orders, similar to the studies of crystalline orders. Methods X-ray projection imaging. We use X-ray projection absorption imaging to measure the average packing fraction F. The granular particles used for all experiments in this study are glass particles (Duke Scientific, USA) with 200 ± 15 mm particle diameter and a slight polydispersity of around 3%. The experiment is conducted with a 27-keV monochromatic X-ray beam at the BL13W1 beamline of the Shanghai Synchrotron Radiation Facility (SSRF). F are measured at ten different tapping numbers (evenly spaced on the logarithmic scale) for each compaction curve.
To obtain F, we also take flat-field images when the packing is outside the X-ray field-of-view. The projection images with and without the packing have intensity distributions I(x,z) and I 0 (x,z), respectively, where (x,z) denotes the coordinates of a pixel. According to Beer's light absorption law, F of the packing can be calculated as where l 0 is the attenuation length of glass at an X-ray energy of 27 keV, R and H are the radius and height of the container in the field-of-view. The integration is taken over the area covered by glass particles. The influence of the acrylic container has also been corrected. The F values obtained are consistent with independent tomography measurements as shown in Fig. 1c.
X-ray microtomography. We use microtomography to obtain the 3D packing structures prepared by tapping, hopper deposition and flow pulse protocols. These three protocols can cover a wide range of F from 0.572 to 0.634. Using micro-tomography, we also investigate the correlation between dynamics and structure by tracking the displacements of all particles and the corresponding structural evolution of the packing for a consecutive number of tapping steps in the tapping experiment.
In the tapping protocol, we fill a 9-mm ID acrylic cylindrical container with particles to B1 cm in height, and we use an electromagnetic exciter to tap the container. Packing with different F ranging from 0.618 to 0.634 is obtained by varying the tapping intensity G, which is measured by an accelerometer as the ratio between the peak-to-peak acceleration and the gravitational acceleration. The tapping consists of a single cycle of 60-Hz sine wave spaced with 0.5 s intervals to allow the system to relax completely. A total of 1,000 taps are applied on each packing with different G to reach steady state. In the hopper deposition protocol, we first place a hopper with its outlet touching the base of the cylindrical container, and fill the hopper with particles. The hopper is then slowly lifted up with a step motor to let the particles drain gradually from the outlet. The packing formed is cylindrical at the bottom, with a conical top. The F prepared by this protocol is 0.598. In the flow pulse protocol, the particles are placed in an acrylic cylindrical tube filled with water. The tube is sealed with a fine copper mesh at the bottom from a water inlet. The packing is prepared by subjecting the particles to a sequence of flow pulses generated from a syringe pump. The particles are allowed to fully settle between pulses. F ranging from 0.572 to 0.588 are obtained by varying the flow velocity.
The X-ray micro-tomography experiment is carried out at both the 2BM beamline of the Advanced Photon Source (APS) at Argonne National Laboratory, and the BL13W1 beamline of SSRF. At Advanced Photon Source, the 'pink' X-ray beam from a bending magnet source with a median energy of 27 keV is used for the high-speed tomography image acquisitions, and the single exposure time is 3 ms. Each tomographic scan consists of 1,500 projection images. At SSRF, the monochromatic 27 keV X-ray beam is used and the single exposure time is 40 ms. Each tomographic scan consists of 720 projection images. The 3D structures are first reconstructed using the conventional filtered back-projection algorithm. A marker-based watershed image segmentation technique is then implemented to obtain all particles' positions and sizes.
For the packing prepared using the tapping and flow pulse protocols, each reconstructed 3D structure consists of approximately 17,000 particles after excluding particles within four particle diameters from the container boundary. For the packing prepared using the hopper deposition protocol, B2,600 particles are used after excluding those in the conical top and boundaries. For each reconstructed packing structure, we conduct a Voronoi tessellation, and the average packing fraction F is obtained by averaging the local packing fractions F loc , which is defined as the ratio between the volume of each particle and its Voronoi cell. NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9409 ARTICLE Calculation of the convection displacement v(r i ). We coarse-grain the whole packing into sub-volumes of cubic shape with the size of each cube about three particle diameters. Since the convection is steady after extensive tapping, we calculate the convection displacement by simply averaging spatially and temporarily of the displacements of all particles inside each cube during the full tapping sequence. We also prove that the results are not sensitive to the coarsegraining size by varying the cube size from two to five particle diameters, and it turns out that the value of Dr 2 only slightly depends on the cube size, and its correlation with d remains approximately unchanged.
Percolation of various structural orders. Various structural mechanisms of the glass transition have suggested the existence of different local structural orders, such as icosahedral order, crystalline order as the driving mechanism of the glass transition. To test all these different mechanisms, we define their corresponding structural order parameters to determine whether there exist significant increases in their spatial correlation lengths as F increases. The spatial correlations are analysed based on percolation analyses on the Delaunay or Voronoi networks of the packing. In both networks, each cell is identified as a site, and the common surface between two sites as a bond. Sites are coloured according to the values of various structural order parameters, or they are coloured randomly. Coloured sites that are connected to each other through bonds can form clusters, and they finally percolate the whole network when enough sites are coloured.
We define two structural order parameters, d and tetrahedricity D based on the Delaunay cells 23 .
d ¼ e max À 1, where e max is the longest tetrahedron edge in units of average particle diameter. Tetrahedricity is the average tetrahedron edge length.
We define the local packing fraction F loc and local contact number Z as two possible order parameters for the free volume theory. We also calculate standard bond orientational order (BOO) parameters that have been commonly used to identify local crystalline or icosahedral orders 55 . Additionally, in our previous study, we found that there exists correlation between the local anisotropy index b 0;2 1 of the shape of the Voronoi cell and certain locally favoured structures with fivefold symmetry 19 . We therefore also include it as a possible order parameter.
Local packing fraction F loc ¼ Vg Vvoro . Local contact number Z. Neighbouring particles whose centre-to-centre distance to one particle is shorter than a threshold r c ¼ 1.011s are defined as its contacting neighbours. s is the average particle diameter. r c is determined from a previously established analysis protocol.
Bond orientational order. We use a modified definition of BOO in which each bond is weighed by the area of its corresponding Voronoi facet 56 : Q lm ¼ P k i¼1 Ai A Y lm ðy i ; j i Þ where (y i ,j i ) is the angular position of the i th Voronoi neighbour in the spherical system of the central particle, A i is the area of the Voronoi facet shared by the central particle and its i th neighbour, and A is the total surface of the Voronoi cell. Y lm are spherical harmonics. The local BOO parameters are defined as: q l ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4p 2l þ 1 P l m¼ À l Q lm j j 2 q and where ð l m1 l m2 l m3 Þ are Wigner 3j symbols. Specifically, we calculated the BOO parameters q 4 , q 6 ,ŵ 4 andŵ 6 , whereŵ l ¼ w l =ð P l m¼ À l Q lm j j 2 Þ 3=2 . The local Voronoi cell anisotropy index b 0;2 1 . It is calculated using a Minkowski tensor W 0;2 1 ¼ R n ndA, defined as the surface integral of the tensor-valued self-product of the bounding surface normal n (ref. 57). b 0;2 1 is defined as the ratio of the smallest and largest eigenvalues of W 0;2 1 . The value of b 0;2 1 ranges from one (isotropic shape) to zero (a line or a plane).
The purpose of the percolation analyses is to identify the spatial correlation properties of these local structural order parameters, that is, whether they form clusters or distribute randomly. To achieve this goal, we colour cells based on the values of above structural order parameters. Although the exact nature of the amorphous structural order remains unknown, it is reasonable to assume they are locally compact. Therefore, they have smaller d and D, and larger F loc and Z. For other structural order parameters, we calculated their Pearson correlation coefficients with F loc and found that b 0;2 1 and q 6 show positive correlations, whereas q 4 ,ŵ 4 andŵ 6 show negative correlations. Therefore, larger b 0;2 1 and q 6 , and smaller q 4 ,ŵ 4 andŵ 6 correspond to more compact local structures.
The colouring process works as follows: cells or simplexes with xo(4)x c (o for d, D, q 4 ,ŵ 4 andŵ 6 while 4 for F loc , Z, b 0;2 1 and q 6 ) are coloured, where x ¼ d; D; F loc ; b 0;2 1 ; q 4 ; q 6 ;ŵ 4 ;ŵ 6 is one of the above structural order parameters, and x c is the corresponding threshold. To compare with a random-colouring process, for each x c value, we also color the same number of cells or simplexes as those with xo(4)x c , except that these cells are chosen randomly.