Lifetimes and lengthscales of structural motifs in a model glassformer

Abstract

We use a newly-developed method to identify local structural motifs in a popular model glassformer, the Kob–Andersen binary Lennard-Jones mixture. By measuring the lifetimes of a zoo of clusters, we find that 11-membered bicapped square antiprisms, denoted as 11A, have longer lifetimes on average than other structures considered. Other long-lived clusters are similar in structure to the 11A cluster. These clusters group into ramified networks that are correlated with slow particles and act to retard the motion of neighbouring particles. The structural lengthscale associated with these networks does not grow as fast as the dynamical lengthscale ξ₄ as the system is cooled, in the range of temperatures our molecular dynamics simulations access. Thus we find a strong, but indirect, correlation between static structural ordering and slow dynamics.

---

1 Introduction

The nature of the rapid increase in viscosity as liquids are cooled toward the glass transition is the subject of many theoretical approaches, however there is no consensus on its fundamental mechanism.^1–4 One plausible scenario is the emergence of self-induced memory effects upon supercooling of liquids, which causes slow dynamics.⁵ However the recent discovery of dynamic heterogeneities, i.e., spatial heterogeneities in the relaxation dynamics that emerge on supercooling,^6–9 is suggestive of the importance of a growing dynamic length scale in the slowing down approaching the glass transition.

In addition to these dynamical phenomenona, the idea of a structural change leading to vitrification has a long history.¹⁰ Sir Charles Frank suggested that upon supercooling liquids would form polyhedral motifs such as icosahedra that do not fill space. A related approach, geometric frustration,¹¹ suggests that the glass transition can be thought of as a manifestation of a crystallisation-like transition that would occur in curved space (where the structural motifs tessellate¹²), but that growth of the “crystal nuclei” of polyhedral motifs is frustrated in Euclidean space. Another approach based on frustration against crystallization (due to random disorder effects or competing orderings)^13,14 addresses the origin of slow dynamics and the physical factors controlling glass-forming ability within the same framework.

It has become clear that a range of glass formers exhibit a change in structure upon the emergence of slow dynamics,^15,16 and it is debated as to whether there exists a static lengthscale that underlies the growing lengthscale for the dynamical correlations. Broadly speaking two types of structure have been identified: spatially extendable crystal-like ordering,^15,17–19 and non-extendable polyhedral ordering.^20–28 For the former it has been suggested that critical-like fluctuations of crystalline order are the origin of dynamic heterogeneities in certain classes of supercooled liquid.¹⁷ The latter concerns particles organised into polyhedra that cannot tile Euclidean space due to geometrical constraints.^10,11 Instead they form ramified structures with a fractal dimension that is less that the dimensionality of the system, i.e. non-extendable ordering. Some metallic glasses have been shown to exhibit this second type of ordering.^29,30 The exact relationship of the polyhedral order to the dynamic heterogeneities is unclear, however measurements have shown that the polyhedral domains are slow to relax.¹⁶ Although the two types of orderings have a different nature, we note that both are induced to lower the free energy locally in the situation where its global minimization (crystallization) is prohibited.¹⁴ In relation to this, we note that even polyhedral ordering often has some connection to crystalline order, for example icosahedra in quasicrystal approximants.²⁶ Furthermore, even in the case where extendable order dominates slow dynamics, non-extendable polyhedral order competes with extendable crystal-like order in systems such as 2D spin liquids¹⁵ and polydisperse hard spheres.¹⁹

One of the main difficulties with identifying structural correlations in supercooled liquids is that it is not known a priori which (if any) type of static order is important for the dynamic slowdown. Frequently studies have employed structure detection methods, such as Voronoi face analysis,³¹ common neighbour analysis,³² bond orientational order analysis,³³ or the topological cluster classification employed here^34–36 that search for predefined types of structural “motif” and investigate the change in population as a function of temperature.^{20,22,28,35,37} However recently a number of “order-agnostic” schemes have been devised to identify structural correlations without first having to define what structures will be searched for.^38–44 In relation to this, it has recently been pointed out that to access such hidden structural ordering in apparently random structures it is crucial to focus on many-body correlations.^45,46

To strengthen the link between structure and glassy behaviour three types of evidence have been presented. Firstly, many-body structural correlation functions have been presented that clearly show structural changes occur on supercooling towards the glass transition.^16,26,46 Secondly, dynamically slow regions have been correlated with different types of local ordering.^16,24 Thirdly, the presence of structural and dynamic length scales that grow similarly has been sought.^15,17,18 This is motivated by strong evidence that, for sufficient cooling (which leads to relaxation timescales that may or may not be accessible to computer simulations as we employ here), an increasing dynamical lengthscale necessitates an increasing lower bound to a structural lengthscale.⁴⁷ The case of growing structural and dynamic length scales is controversial. One of us has identified a direct correspondence between the growing dynamical and structural lengthscales in polydisperse systems displaying crystal-like ordering,^17,18,45 while others have claimed that structural lengthscales are decorrelated from the dynamic lengthscales and only grow weakly on cooling.^48–50

Local structural ordering has yet to be found in some glassforming systems. For these systems no one-to-one correspondence between an order-agnostic structural lengthscale and the dynamical lengthscale has been found,^41,44,48,51 but static perturbation analysis has suggested a growing structural lengthscale.⁵²

The Kob–Andersen binary Lennard-Jones system⁵³ of interest here is known to display polyhedral ordering.¹⁶ In this system, weakly growing structural lengthscales have been identified by the order-agnostic “point-to-set” analysis,⁴³ while other approaches using static perturbation of inherent structures³⁹ and finite size scaling⁵⁴ find a stronger increase in static lengthscales. The polyhedral ordering in the Kob–Andersen (KA) mixture,⁵³ take the form of the bicapped square antiprism [“11A” in the topological cluster classification (TCC) nomenclature^34–36,55 owing to its original identification as minimum energy cluster of the Morse potential^56,57]. This structure has been linked to slow dynamics¹⁶ and frustrates crystallisation, which in the KA mixture occurs spontaneously by phase separation into two face-centred cubic lattices.⁵⁸ We note that crystals based on 11A could in principle coexist with other structures [the stoichimetry of the 11A crystal with a small particle in the centre (see Section 4.3 and Fig. 5) is not compatible with the KA mixture⁵⁹]. Another possibility is four-fold symmetric crystals which have been predicted as low-lying energy minima for the KA mixture.⁶⁰ An order parameter associated with the formation of 11A clusters has been shown to control the dynamical phase transition in trajectory-space, a hallmark of dynamical facilitation theory for the glass transition,⁶¹ indicating that the transition has both structural and dynamical character.⁶² Here we study the spatial correlations between the domains of 11A in the supercooled liquid.

We consider the lifetimes of a multitude of structures in the supercooled liquid using the TCC algorithm. We detail our simulation protocol in section 2, and briefly review the KA model's dynamical behaviour in section 3. We use the TCC to identify any structural changes that occur in these mixtures on cooling towards the glass transition in section 4. We then study in section 4.2 the lifetimes of the clusters that are found at deeply supercooled state points in order to gauge which structures are likely candidates to be associated with slow domains of dynamic heterogeneities. Finally we analyse how correlation lengths for the domains of structured particles are related to the growing dynamic lengthscale in section 5 before concluding.

2 Model and simulation details

The Kob–Andersen (KA) binary mixture is composed of 80% large (A) and 20% small (B) particles of the same mass m.⁵³ The nonadditive Lennard-Jones interactions between each species, and the cross interaction, are given by σ_AA = σ, σ_AB = 0.8σ, σ_BB = 0.88σ, ε_AA = ε, ε_AB = 1.5ε, and ε_BB = 0.5ε. The results are quoted in reduced units with respect to the A particles, i.e. we measure length in units of σ, energy in units of ε, time in units of , and set Boltzmann's constant k_B to unity. The interactions are truncated and smoothed using the Stoddard–Ford method.⁶³ The truncation lengths are in proportion to the interaction lengths,⁵³i.e. r^AA_tr = 2.5, r^AB_tr = 2.0 and r^BB_tr = 2.2. The simulations consist of N = 10 976 particles in 3D with periodic boundary conditions such that N_A = 8781 and density ρ = 1.2. The α-relaxation time τ^A_α for each state point is defined by fitting the Kohlrausch–Williams–Watts stretched exponential to the decay of the intermediate scattering function (ISF) of the A-type particles.

Equilibrated samples at each temperature were prepared by simulating for 100τ^A_α in the canonical NVT-ensemble using the Nosé–Poincaré thermostat with coupling parameter 1.0.⁶⁴ The thermostat was switched off, then further equilibration was performed in the microcanonical NVE-ensemble using the velocity Verlet algorithm for 1000τ^A_α. On completion of the equilibration process, trajectories of length 300τ^A_α were sampled for analysis. Colder state points were obtained in a step-wise fashion by quenching instantaneously from an equilibrated configuration of the previous higher temperature state point. The stability of the deep quenches was checked by ensuring there was no time evolution in the ISF, the partial radial distribution functions g_AA(r), g_AB(r) and g_BB(r), or the number of clusters detected by the TCC algorithm across the trajectories. Crystallisation was not seen in any of the simulations.

3 Dynamical behaviour

In Fig. 1 the relaxation time τ^A_α as a function of inverse temperature is plotted. For the equilibrium liquid state points at high temperatures, the relaxation times are well fitted by an Arrhenius function. As the temperature is lowered a cross-over in the dynamical behaviour occurs and the relaxation times increase faster than predicted by the Arrhenius equation.^2,65,66

Fig. 1 Increase in the alpha relaxation time, τ^A_α, on cooling. The data are fitted with a hybrid Arrhenius-VFT fit (solid lines). The dotted lines indicate the relaxation times predicted by the high-T Arrhenius fit. Orange dashed line indicates the crossover from Arrhenius to super-Arrhenius (T*).

We fit the two regimes for the relaxation time delimited by an onset temperature for slow dynamics T*.¹⁶ For T > T* an Arrhenius form is used, while for lower temperatures the Vogel–Fulcher–Tammann (VFT) equation is fitted.^67–69

We set T* = 1.00 and fit the Arrhenius equation finding τ_∞ = 0.0693 and E_∞ = 2.91. For the VFT fit the fragility parameter is found to be D = 7.48, the VFT temperature is T₀ = 0.325 and τ_∞′ is set to ensure continuity of the fit at T*. While it is possible to use other fitting forms,^70,71 and although VFT may be physically reasonable,¹ there remains no clear consensus as to which form best describes data such as those plotted in Fig. 1.⁷² Nonetheless, following our previous study,²⁸ here we choose this fitting procedure to reflect the onset of slow dynamics (super-Arrhenius) for T < 1^2,66 and, over the range of temperatures we consider, find good agreement with the assumption of a crossover to a VFT regime.

4 Structural analysis

4.1 Fraction of particles participating within clusters

We analyse how the particles in the supercooled liquids are structured using the topological cluster classification algorithm. This algorithm identifies a number of local structures as shown in Fig. 2, including those which are the minimum energy clusters for m = 5 to 13KA particles in isolation. The first stage of the TCC algorithm is to identify the bonds between neighbouring particles. The bonds are detected using a modified Voronoi method with a maximum bond length cut-off of r_c = 2.0 for all types of interaction (AA, AB and BB).^34,36,55 A parameter which controls identification of four- as opposed to three-membered rings f_c is set to unity thus yielding the direct neighbours of the standard Voronoi method.^{36,55,73–75}

Fig. 2 Clusters detected by the topological cluster classification. Highlighted are minimum energy clusters for the Kob–Andersen system. The colours of the particles and the bonds are pertinent to the detection method of the clusters.^36,55

In Fig. 3 we plot the fraction of particles detected within each type of cluster, N_C/N. The onset temperature for slow dynamics is indicated by the orange dotted line on each of the plots. It is clear from these order parameters that the liquid sees continuous changes in local structure as it is cooled. The majority of clusters see an increase in their numbers upon cooling. All of the particles are identified within certain simple structures such as the 5A triangular bipyramid irrespective of temperature, while other more complicated clusters, such are 10K, 11B and HCP, are almost never seen.

Fig. 3 The fraction of particles participating in each cluster type for the KA mixture. The dotted orange lines mark the onset temperature of slow dynamics T*. (a) Clusters 5A to 8K, (b) 9A to 11A, (c) 11B to 12D, (d) 12E to 13K and the crystal clusters.

Given the variety and range of structural changes that occur, it is not clear a priori which of the structures, if any, are important for the formation of dynamic heterogeneities on cooling. It should not be assumed that just because structural changes occur within the supercooled regime that they are responsible for the formation of dynamic heterogeneities, as structural changes also occur in the Arrhenius regime which is not characterised by glassy behaviour. The numbers of particles within 11A, 11W, 12B, 12E and 12K clusters show the largest relative increases in the super-Arrhenius regime. Moreover, the rate of increase in these clusters grows upon further supercooling.

We also find clusters in which almost all particles participate, while other clusters are only seen in trace quantities. In general the trend for clusters where N_C/N changes significantly on cooling is for it to increase monotonically. However this is not always the case, as seen for the numbers of particles within 9A and 9X clusters which decrease on cooling. The question remains as to how to determine the contribution of the different structures to the glassy behaviour and dynamical heterogeneities of the supercooled liquid.

4.2 Cluster lifetime distributions

In order to identify which clusters might be relevant to the slow dynamics, we employ the dynamic topological cluster classification algorithm^28,36 to measure the lifetimes of the different TCC clusters at the lowest temperature state point. A lifetime τ_{[small script l]} is assigned to each “instance” of a cluster, where an instance is defined by the unique indices of the particles within the cluster and the type of TCC cluster. Each instance cluster occurs between two frames in the trajectory and the lifetime is the time difference between these frames. Any periods where the instance is not detected by the TCC algorithm are shorter than τ^A_α in length, and no subset of the particles becomes un-bonded from the others during the lifetime of the instance.

The measurement of lifetimes for all the instances of clusters in these N = 10976 simulations is intensive in terms of the quantity of memory required to store the instances, and the number of searches through the memory required by the algorithm each time an instance of a cluster is found to see if it existed earlier in the trajectory. Therefore we do not measure lifetimes for the clusters where N_C/N > 0.8, since the vast majority of particles are found within such clusters and it is not immediately clear how dynamic heterogeneities could be related to structures that are pervasive throughout the whole liquid.

In Fig. 4 we plot the lifetime autocorrelation function P(τ_{[small script l]} ≥ t) for T = 0.498. Fig. 4 clearly shows that the most persistent or the longest lived of the different types of clusters are the 11A bicapped square antiprisms. All other clusters display lifetime autocorrelation functions that decay more quickly than 11A. The long-time tail of the 11A autocorrelation function indicates that some of these clusters preserve their local structure on timescales far longer than τ^A_α. As we shall see below, this effect is enhanced when the 11A group into domains.

Fig. 4 Lifetime autocorrelation functions for the clusters P(τ_{[small script l]} ≥ t) for the lowest temperature state point T = 0.498.

The two next-slowest decaying clusters are 11W and 12K. 12K is a KA minimum energy cluster, formed by bonding an additional particle to 11A. There is a high degree of overlap between the 11W clusters and the 11A clusters as their bonding is similar, and only small fluctuations are required for an 11A to be reclassified as an 11W. However the faster-decaying clusters also contain KA minimum energy clusters, for example 13K and 10K. Moreover the lifetimes of all cluster types hold no simple relationship to their size and frequency of occurrence. For example the n = 11 particle 11F cluster is much more numerous than 11A, yet displays far quicker decay of P(τ_{[small script l]} ≥ t). There is also no monotonic trend in the lifetime of the ground state clusters with the cluster size, as the 10K and 13K decay faster than the 12K and 11A. These results clearly demonstrate that the average lifetime of each type of cluster is a property of the local ordering of the particles rather than the size of the cluster or its pervasiveness.

The fast initial drops of P(τ_{[small script l]} ≥ t) reflect the existence of large numbers of clusters with lifetimes τ_{[small script l]} ≪ τ^A_α. The lifetimes of these clusters are comparable to the timescale for beta-relaxation where the particles fluctuate within their cage of neighbours. It could be argued that these clusters arise spuriously due to the microscopic fluctuations within the cage, and that the short-lived clusters are not representative of the actual liquid structure. However almost no 11A are found at higher temperatures, cf.Fig. 3(b), where microscopic fluctuations in the beta-regime also occur. We have not yet found a way to distinguish between the short and long-lived 11A structurally, so we conclude that the measured distribution of 11A lifetimes, which includes short-lived clusters, is representative of the true lifetime distribution. However, given the structural similarity of 11A and 11W, small fluctuations leading to reclassification could contribute to the drop at short times. We will see below that as 11A overlap, one particle may be a member of multiple clusters and that the majority of particles found in short-lived 11A also participate in longer-lived 11A. In other words short- and long-lived 11A mainly lie in the same regions of the liquid.

Another interpretation for the initial drop in P(τ_{[small script l]} > t) of the clusters is that our intuition that there is constant local structure in the beta-regime, which then relaxes on the timescale of the alpha-regime, is incorrect. It may be the case that microscopic ballistic and “cage-rattling” motions are enough to reorder local structures without relying on the “cage-hopping” motions of the alpha-regime. It has been seen in previous studies that deeply supercooled liquids can crystallise on a timescale before the diffusive range of the mean squared displacement is reached, and occurs with most particles moving by less than one diameter.^76–79 Those results demonstrate that significant changes in local structure (i.e. liquid-like to crystalline) are possible with only small movements of the particles, which could explain the initial drops of P(τ_{[small script l]} ≥ t) that occur on a timescale ≃0.1τ^A_α.

4.3 Composition and dynamics of particles in long-lived clusters

The structural analysis above was performed by treating all particles identically with the TCC algorithm. Here we examine the composition of the 11A bicapped square antiprisms in terms of A- and B-species. The 11A cluster consists of a central particle surrounded by 10 outer (or shell) particles. We find that the central particle is a B-species in >99% of all instances of 11A clusters. This is a different composition than the crystal structures found to be low-lying energy minima for the KA mixture,⁶⁰ and may be related to the stoichiometry of the system. Thus the 11A we find are unrelated to any underlying crystal.

In Fig. 5 we plot the compositions of the shell particles of the 11A clusters. The majority of 11A clusters have m_A = 10 A-species in the shell of the cluster. We note that this arrangement maximises the number of AB bonds for the central B-particle, which is energetically favourable for the central B-particle with the KA Lennard-Jones interactions.

Fig. 5 Composition of the 11A clusters at T = 0.498. Almost all 11A have a small B-particle at the centre. m_A is the number of A-species in the shell of the cluster. The height of the bars show the relative proportions that each of the compositions occurs.

In Fig. 6 we examine how the dynamics of the 11A clusters translates into the dynamics of individual particles. We show in Fig. 6(a) that the number of particles within 11A clusters as a function of the cluster lifetime. Although there is a fast initial drop in the lifetime autocorrelation function of 11A clusters on the beta-relaxation timescale (Fig. 4, solid purple line), as the 11A overlap there remains a significant fraction of the particles these clusters with lifetimes comparable to the dynamic heterogeneities ≈τ^A_α. The difference between N_11A(τ_{[small script l]} ≥ 0.1τ^A_α)/N and N_11A(τ_{[small script l]} ≥ 0)/N indicates that only 6% of the particles are members of 11A with τ_{[small script l]} < 0.1τ^A_α and not a member of an 11A with a longer lifetime as well.

Fig. 6 Dynamics of the particles within 11A clusters in the KA mixture. (a) The fraction of particles participating in 11A clusters with lifetime τ_{[small script l]} > t. N_11A(τ_{[small script l]} ≥ 0)/N = 0.24. (b) The mean squared displacement of particles identified initially within 11A polyhedra of various lifetimes.

Fig. 6(b) shows the mean squared displacement (MSD) of the particles identified initially within 11A clusters (coloured lines) and compares this to the system-wide MSD (black line). The MSD is defined as the ensemble average 〈δr²(t)〉 = 〈|r_i(t + t₀) − r_i(t₀)|²〉 for the subset of particles of interest (indexed by i). All of the particles within 11A relax more slowly than the system-wide average (black line), and the time they take to attain diffusive motion increases as the lifetime of the 11A in which they participate in at t₀ increases. In other words the longer the lifetime of the 11A cluster, the slower the particles become. Since some 11A last for very long times [Fig. 4], it is expected that these particles may exhibit very low mobilities as they maintain some of their nearest neighbours throughout (e.g. the central particle in an 11A cluster will always have the same shell particles as its nearest neighbours). For the longest lived 11A (blue lines) there appears to be a super-diffusive regime after the initial sub-diffusive regime, indicating that the particles in these clusters may be hopping out of their cage of neighbours as the 11A structure relaxes.

4.4 Analysis of structured domains

On cooling, the number of 11A clusters in the KA mixture increases. At high temperatures the clusters are generally isolated from one another [Fig. 7(a)]. As the temperature is lowered and the number of clusters increases, domains of clustered particles form. We now analyse the character of these domains of clustered particles and determine the effect the domains have on individual particle dynamics.

Fig. 7 Analysis of the domains of 11A clusters. (a)–(c) Domains form on cooling from high to low temperature (slices through 3D simulation box). Particles in 11A clusters are shown full size in red, other particles are blue dots. (d) The radius of gyration R_G of the domains versus the number of particles in the domain n for T = 0.498. R_G is well fitted by n^0.48 indicating the domains have a fractal dimension d_f ≃ 2. (e) The mean lifetime of 11A clusters [small tau, Greek, macron]_{[small script l]}versus the domain size n. (f) 11A domains affect the motion of neighbouring particles. The MSD 〈δr²(τ_h)〉 of non-11A particles as a function of distance from 11A domains d (solid line). The dotted line is the MSD over τ_h of all particles not in 11A clusters and independent of the distance from an 11A domain (d).

The domains of 11A that form on cooling in the KA mixture are shown in Fig. 7(a)–(c). At high temperature 11A are predominantly isolated [Fig. 7(a)]. Upon cooling, the 11A overlap and join together [Fig. 7(b)] to form (transient) networks at low temperatures [Fig. 7(c)].

In order to investigate the structure of the domains, we calculate their radius of gyration,

where n is the number of particles in the domain and the double sum extends over all pairs of particles in the domain. For configurations where the domains of 11A percolate throughout the simulation box R_G cannot be defined. These configurations are rare at T = 0.498 and are excluded from the analysis. The consequences of percolating 11A domains are discussed further below.

The radius of gyration shows a power law growth in the size of the domain n with an exponent of 0.48 [Fig. 7(d)]. In other words the domains have a fractal dimension d_f ≃ 2, indicating that they are not space-filling. The individual 11A have enhanced stability as the size of the 11A domains grow [Fig. 7(e)]. The time [small tau, Greek, macron]_{[small script l]} is the mean lifetime of 11A clusters constituting a domain of size n in a configuration. For T = 0.498 the general trend is that the mean 11A lifetime increases with the size of the 11A domains. Fig. 7(e) indicates that there is particularly stable arrangement of two overlapping 11A clusters (n ≃ 16) relative to domains from similar size. The mean lifetime [small tau, Greek, macron]_{[small script l]} doubles between isolated 11A (n = 11) and extended domains (n ≈ 1000).

We now consider the effect the domains have on the remainder of the system. In Fig. 7(f), we plot the MSD of the particles not in 11A domains 〈δr²(τ_h)〉 against the distance d from the nearest 11A particle at time t = t₀. The time τ_h ≃ τ^A_α is the time of the maximum in the dynamic susceptibility χ₄(t) defined below.⁸⁰ For distances d < 0.96 the non-11A particles are mainly B-species due to the nonadditive nature of the KA Lennard-Jones interactions. These particles are more mobile than the majority A-species, independent of d, due to their smaller size. The number of particles with d < 0.96 of a 11A domain is around 10% of the system, i.e. half of all the B-species.

Moving further away from the 11A domains, there is then a region of particles with d ≃ 1 with reduced mobility compared to the average for non-11A particles [dotted line in Fig. 7(f)]. Around 37% of the particles are found in this region. Subsequent neighbours of the 11A domains for d > 1.26 have increased mobility relative to the average. Therefore, excluding the minority B-species, the first nearest neighbours of the domains have suppressed mobility compared to the average, indicating coupling between the structured domains of 11A clusters and the dynamics of the neighbouring particles. Correspondingly the second and third shells of neighbours to the 11A domains have higher mobility compared to the average, indicating a hierarchy of spatial dynamics related to the domains of 11A clusters.

5 Correlation lengths

5.1 Dynamic correlation lengths

Finally we consider whether the structured domains of particles are related to the increasing dynamic correlation lengths in supercooled liquids. In order to do this, we calculate the dynamic correlation length ξ₄, following Lačević et al.⁸⁰ Details of this procedure can be found elsewhere.^28,80,81ξ₄ has been previously calculated for the Kob–Andersen model and our values correspond closely to those in the literature.⁸¹ The dynamical correlation length ξ₄ is obtained by analogy to critical phenomena.⁸⁰ A (four-point) dynamical susceptibility is calculated aswhere

The overlap function w(|r_j(t + t₀) − r_l(t₀)|) is defined to be unity if |r_j(t + t₀) − r_l(t₀)| ≤ a, 0 otherwise, where a = 0.3. The dynamic susceptibility χ₄(t) exhibits a peak at t = τ_h, which corresponds to the timescale of maximal correlation in the dynamics of the particles. We then construct the four-point dynamic structure factor S₄(k, t):

where j, l, m, n are particle indices and k is the wavevector. For time τ_h, the angularly averaged version is S₄(k, τ_h). The dynamic correlation length ξ₄ is then calculated by fitting the Ornstein–Zernike (OZ) function to S₄(k, τ_h), as if the system were exhibiting critical-like spatio-temporal density fluctuations,to S₄(k, τ_h) for k < 2.⁸⁰ The resulting ξ₄ are plotted in Fig. 8.

Fig. 8 Comparison of static and dynamic correlation lengths. Structural ξ_RG (squares) and ξ_SC (crosses), and dynamical ξ₄ (circles) correlation lengths. ξ₄ is fitted with a power law (solid line), which diverges at a value T_C = 0.47.

We carried out an unconstrained fit to the ξ₄ data, according to ξ₄(T) = ξ⁰₄(T − T_C)^−ν. The line is plotted in Fig. 8. We find the “critical exponent” is ν = 0.588 ± 0.02, the “critical temperature” is T_C = 0.471 ± 0.002 and the prefactor is ξ⁰₄ = 0.59 ± 0.02. Under the caveat that obtaining ξ₄ from fitting S₄ in limited size simulations is notoriously problematic^81,82 and thus any numerical values should be treated with caution, we observe that the value of T_C is not hugely different to the glass transition temperature found by fitting Mode-Coupling theory to this system, around 0.435.^53,83 We also note that ν = 0.588 lies between mean field (ν = 0.5) and 3D Ising (ν = 0.63) criticality. Here it is worth noting that Onuki and his coworkers pointed out that the dynamical correlation length estimated from four-point density correlations may be affected seriously by thermal low-frequency vibration modes, which may lead to a strong system-size dependence.⁸⁴ They also showed that a bond-breakage correlation length is free from these vibrational modes, and thus a suitable measure of dynamical coherence.

5.2 Static correlation lengths

We consider two static correlation lengths for the domains of particles in 11A clusters. The first method allows for direct comparison with the dynamic lengthscale ξ₄. We define a structure factor restricted to the particles identified within 11A:where N_11A is the number of particles in 11A clusters. We then fit the Ornstein–Zernike equation (eqn (5)) to the low-k behaviour of the angularly-averaged S_11A(k) in order to extract a structural correlation length ξ_11A. This is plotted in Fig. 8. This procedure is akin to the calculation of the dynamic lengthscale ξ₄: first a structure factor is calculated from a selected fraction of the particles (either immobile or structured), and the Ornstein–Zernike expression used to extract a correlation length.

The second lengthscale we consider for the structured particles is derived from the radius of gyration of the domains of clusters. We define

ξ_RG = R^11A_G(〈n〉/m)^1/d_f, (7)

where R_G^11A is the radius of gyration of a single cluster, 〈n〉 is the ensemble average of the domain size, m is the number of particles in the cluster, and 1/d_f is the exponent of the power law fitted to R_Gversus n. This correlation length does not probe the correlations between the domains, as per ξ_11A, rather it characterises the growth in size of the domains on cooling until a percolation transition is reached.

The temperature behaviour of the different correlation lengths is shown in Fig. 8. All three correlation lengths increase on cooling, however the manner in which each of the lengths increases is quite different. The main result is that the growth in the dynamic correlation length ξ₄ is not matched by the growth in the structural correlation length ξ_11A. Thus we do not find one-to-one correspondence between the behaviour of structural and dynamic correlation lengths.

5.3 Discussion

The fact that different behaviour between the dynamic and static lengths is found is in agreement with some recent studies on 2D and 3D systems,^{43,44,48–50} however we note that a one-to-one correspondence in the growth in a lengthscale relating to static crystalline order and the dynamic correlation length for polydisperse quasi-hard sphere systems has been found.^15,17,18,85

The structural order we find is distinct from the crystalline structure and is thought to frustrate crystallisation. We note that in the studies on 3D systems that have found one-to-one correspondence between the structural and dynamical lengthscales,^17,85 the relative increase in the static lengthscales going from state points in the Arrhenius regime to into the supercooled regime is no more than a factor of 2. Increases of comparable magnitude have been found in studies using point-to-set measures,⁴³ and more indirect methods^39,52 rather than focusing on explicit structures. In Fig. 8 our static length ξ_11A shows an increase of the same order between Arrhenius and supercooled state points, however the increase in ξ₄ is relatively much greater (factor of ≈6 or more between high and low T). We note a recent 2D study that suggests one-to-one correspondence in lengthscales for crystalline order and dynamic heterogeneities on quasi-hard spheres breaks down with the addition of attractions between the particles.⁸⁶ However, we also note that 3D polydisperse hard-sphere and Lennard-Jones systems have the same link between the correlation length of crystal-like order and slow dynamics, albeit over a limited range of correlation length.¹⁷ These points need to be clarified in the future.

The structured domains form rarefied networks with d_f ≃ 2. These networks do not fill space, whereas it is thought that any critical-like nature of the dynamic heterogeneities would imply that the domains have d_f ≃ 2.5.⁸⁷ Therefore geometrically it appears that the domains of slow and structured particles may be different, at least at the state points we have accessed. Moreover the number of particles in the domains is unlikely to coincide across a range of temperatures, as the N_11A for the structured particles increase monotonically on approaching the glass transition, while the number of particles selected by Q(t) is relatively constant as a function of temperature.

Here we note that Mosayebi et al.^39,52 found behaviour consistent with critical-like static correlations in the same system as ours, which implies d_f ∼ 2.5. This seems to suggest that 11A together with other clusters or kinetically pinned particles may correspond to their static structures with solid-like nature. Now the temperature range over which we were able to equilibrate our simulations is less than that for which Mosayebi et al.^39,52 present data. Indeed over our range of temperature, the relative increase in static lengthscale they measure is consistent with ours. However, since our d_f ∼ 2.0, our findings are not consistent with Ising-like criticality, although critical-like behaviour at lower temperatures than we have been able to access cannot be ruled out. We also note that the structural motifs investigated here are of intrinsically of discrete nature, whereas the structural measures showing critical-like behaviours (d_f ∼ 2.5), such as bond orientational order^14,17 and static structures^39,52 have a continuous nature.

The MSD of the 11A particles indicates that there is not a one-to-one correspondence between the particles selected in S₄ and S_11A. The particles selected by w for S₄ all have δr²(τ_h) ≤ 0.09 strictly (by the definition of w). However, as can been seen from the red line in Fig. 6(b), the MSD of the 11A particles over the same timescale is 〈δr²(τ_h)〉 ≈ 0.1. Only the particles in the longest-lived 11A clusters [blue lines in Fig. 6(b)] have squared displacements over τ_h comparable to the immobile particles on which S₄ is measured. Thus on this basis direct correspondence between ξ₄ and ξ_11A should not necessarily be expected. We note that this might also be related to the effects of low-frequency vibrational modes.⁸⁴

Note χ₄ has been shown to exhibit dependence upon system size for N ≲ 1000.⁸⁸ We expect that such effects are reasonably small here, and have taken care to only consider temperatures where all our measured lengths are smaller than the system size. While system size effects cannot be ruled out, we do not believe these make a significant impact on the conclusions we draw.

The lengthscale ξ_RG indicates how the size of the domains grows on cooling. In our related study on the Wahnström binary Lennard-Jones glass former,²⁸ where the concentration of clusters was higher, a percolating network of icosahedra was formed. In a study on the Kob–Andersen model considered here,⁶² some of us found evidence for an increase in the population of 11A clusters at lower temperatures than we have been able to access here. Thus if the 11A clusters percolate, this would suggest a diverging structural lengthscale at a temperature higher than either the VFT or MCT temperatures (∼0.325 and ∼0.435 respectively). In any case, percolating domains of 11A clusters do not imply structural arrest since each cluster has a finite lifetime. This scenario contrasts with colloidal gels where a percolating network of local structures leads to dynamic arrest.³⁵

We also note that Fig. 7(f) strongly indicates that the effect of the domains of clustered particles on the surrounding liquid extends around one particle diameter from the domains. This shows that the dynamical effect of the structured particles is hierarchical and not solely limited to the structured domains themselves. However the correlation lengths we measured from a structure factor of the domains and their first nearest neighbours was no greater than the lengthscale ξ_11A for the domains themselves.

Thus the question remains as to the most appropriate structural (and dynamical) correlation length to explain the viscous slowing down in supercooled liquids, and how order-specific correlation lengths, such as ξ_RG and ξ_11A, are related to order-agnostic structural correlation lengths.^{38–40,44,52,89} We conclude this discussion with the following observations.

1. Static correlation lengths have been measured in a variety of systems. Most seem to grow less strongly than dynamic correlation lengths in the regime accessible to computer simulation (and real-space colloid experiments).^{28,38–40,43,44,52,89} Those that are observed to grow in a way comparable to the dynamic correlation length ξ₄ are often related to crystalline order.^15,17,90 Here we note that a bond orientational order parameter is a continuous variable, whereas structural motifs are discrete by the definition. There might also be some effects of low-frequency vibrational modes on the estimation of the dynamical correlation length in our system.⁸⁴ These points need further investigation.

2. We note that geometric frustration suggests a term in R⁵ to be included in a classical nucleation theory like equation, where R is the size of a growing domain of the preferred structure (11A here).¹¹ This R⁵ term strongly suppresses growth of such domains, consistent with ramified structures as we (and others^16,24) find. Now the arguments supporting geometric frustration tended to focus on icosahedra, not 11A bicapped square antiprisms as we find here. However, in our study of the Wahnström model (whose local structural motif is the icosahedron),²⁸ we found very similar behaviour to that reported here. In any case, crystals can be formed of 11A, though not for the 80–20 composition of the KA mixture,⁵⁹ and icosahedra (with the inclusion of Frank–Kasper bonds).²⁶ Assuming crystallisation is avoided, given the presence of such an R⁵ term, one might enquire as to why growth of polyhedral domains is expected in the first place.

3. Apart from ξ₄, dynamic correlation lengths have also been measured by other means.^41,43 In particular, evidence has been found for non-monotonic behaviour of such dynamic correlation lengths around the Mode-Coupling transition.⁴¹ While ξ₄ itself does not exhibit non-monotonic behaviour as such, there is evidence that a different scaling is followed for quenches below the Mode-Coupling transition.⁹¹ Moreover, that our fit of ξ₄ implies divergence at T = 0.47 suggests that some other scaling might be found at lower temperatures than we have accessed. Another way to approach the discrepancy in Fig. 8 is to enquire whether ξ₄ is the “right” dynamical correlation length.⁹² One could even speculate that dynamical correlations are enhanced around the Mode-Coupling transition (as our results imply, along with those of Kob et al.⁴¹), and that at least the lengthscale does not grow significantly at lower temperatures. This would rationalise the estimates of dynamic correlation lengths close to the molecular glass transition (some 10 orders of magnitude slower in relaxation time than our simulations are able to access) which suggest correlation lengths only of a few molecular diameters.^93,94 However, we should note that there is a possible deficiency of the standard dynamical correlation length ξ₄.⁸⁴ We point out that the wavenumber dependence of the transport coefficient, viscosity, is a physically appealing method for estimating the intrinsic dynamical correlation length.^95–97

6 Summary and conclusions

We have demonstrated that by studying the lifetimes of different structural orderings within the Kob–Andersen supercooled liquid, the relatively stable orderings of particles can be detected unambiguously. This method alleviates some of the difficulties in identifying structural correlations relevant to glassy behaviour from the temperature dependency of the number of particles participating in clusters. The most stable cluster found in the Kob–Andersen supercooled liquid is the 11A bicapped square antiprism. The relaxation of particles within these clusters proceeds more slowly as the lifetime of the cluster increases. This is consistent with previous work based on Voronoi polyhedra.¹⁶

The long-lived clusters form rarefied domains on cooling with a fractal dimension d_f ≃ 2, i.e. the structured domains are non-space-filling, at least in the regime we have accessed. The lifetime of the 11A clusters increases markedly with the size of the domains of these clusters. The non-11A particles neighbouring these domains have reduced mobility compared to particles further from the domains, suggesting a link between structure and dynamic heterogeneity at this level of detail. In other words, the network of 11A clusters acts to “pin” its neighbours.

We examined the relationship between the structured domains and the dynamic heterogeneities by considering static and dynamic correlation lengths of structured/slow particles. A static correlation length calculated in a like-for-like manner with the dynamic correlation length was found to grow moderately on cooling, however its increase was outmatched by the growth in the dynamic correlation length ξ₄. The difference in behaviour of the correlation lengths was rationalised by noting that the structured domains grow in a non-space filling manner, and that the correlation between the structured and slow particles is not perfect. The relationship between the our static and dynamic correlation lengths, and other lengthscales for static order, remains an open question, (see section 5.3).

Finally we consider a possible direction for future study. It has been shown recently that an order parameter associated with the population of 11A clusters can be used to drive a first order transition in an ensemble of trajectories.⁶² The susceptibility of the transition to the field coupled to the structural order parameter was found to be higher than when biasing with a field coupled to the dynamical activity, which is the usual method that the transition is accessed. Furthermore a recent study by Singh et al. has shown that “ultra-stable” KA glasses prepared by a vapour deposition technique have high numbers of clusters equivalent to 11A polyhedra.⁹⁸ Together these results provide further evidence for a connection between the atomic level structure, most easily accessed with high-order structural correlation functions, and the glass transition. The biasing of fields coupled to structural order parameters in trajectory space, and possibly in configuration space as well, thus opens up a new route for the preparation of ultra-stable glassy states⁹⁸ pertinent to temperatures well below the Mode-Coupling temperature. The relaxation times inferred for these states are many orders of magnitude higher than those which can currently be prepared with conventional simulations. Study and characterisation of the properties of these states will shed further light on the nature and role of local structure in the glass transition. In other words, while Fig. 8 (like much of the recent literature^{28,39,41,43,44,52}) indicates a decoupling between structural and dynamical lengthscales, other static and/or dynamical measures might finally lead to coupling of structure and dynamics.

Acknowledgements

We gratefully acknowledge stimulating discussions with Rob P. Jack, Thomas Speck and Stephen Williams. A. M. is funded by EPSRC grant code EP/E501214/1. C. P. R. thanks the Royal Society for funding. H. T. acknowledges support from a grant-in-aid from the Ministry of Education, Culture, Sports, Science and Technology, Japan and the Aihara Project, the FIRST program from JSPS, initiated by CSTP. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol.

References

1. A. Cavagna, Phys. Rep., 2009, 476, 51–124.
2. L. Berthier and G. Biroli, Rev. Mod. Phys., 2011, 83, 587–645.
3. F. H. Stillinger and P. G. Debenedetti, Annu. Rev. Condens. Matter Phys., 2013, 4, 263–285.
4. G. Biroli and J. P. Garrahan, J. Chem. Phys., 2013, 138, 12A301.
5. W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory, Oxford University Press, Oxford, 2008.
6. M. M. Hurley and P. Harrowell, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 1694–1698.
7. R. Yamamoto and A. Onuki, J. Phys. Soc. Jpn., 1997, 66, 2545–2548.
8. M. Ediger, Annu. Rev. Phys. Chem., 2000, 51, 99–128.
9. L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti and W. van Saarloos, Dynamical Heterogeneities in Glasses, Colloids, and Granular Media, Oxford University Press, Oxford, 1st edn, 2011.
10. F. C. Frank, Proc. R. Soc. London, Ser. A, 1952, 215, 43–46.
11. G. Tarjus, S. A. Kivelson, Z. Nussinov and P. Viot, J. Phys.: Condens. Matter, 2005, 17, R1143–R1182.
12. D. R. Nelson, Defects and Geometry in Condensed Matter Physics, Cambridge University Press, 2002, p. 392.
13. H. Tanaka, J. Phys.: Condens. Matter, 1998, 10, L207–L214.
14. H. Tanaka, Eur. Phys. J. E, 2012, 35, 113.
15. H. Shintani and H. Tanaka, Nat. Phys., 2006, 2, 200–206.
16. D. Coslovich and G. Pastore, J. Chem. Phys., 2007, 127, 124504.
17. H. Tanaka, T. Kawasaki, H. Shintani and K. Watanabe, Nat. Mater., 2010, 9, 324–31.
18. F. Sausset and G. Tarjus, Phys. Rev. Lett., 2010, 104, 065701.
19. M. Leocmach and H. Tanaka, Nat. Commun., 2012, 3, 974.
20. H. Jónsson and H. Andersen, Phys. Rev. Lett., 1988, 60, 2295–2298.
21. T. Kondo and K. Tsumuraya, J. Chem. Phys., 1991, 94, 8220.
22. T. Tomida and T. Egami, Phys. Rev. B: Condens. Matter, 1995, 52, 3290–3308.
23. R. Jullien, P. Jund, D. Caprion and D. Quitmann, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 6035–6041.
24. M. Dzugutov, S. I. Simdyankin and F. H. M. Zetterling, Phys. Rev. Lett., 2002, 89, 195701.
25. E. Lerner, I. Procaccia and J. Zylberg, Phys. Rev. Lett., 2009, 102, 125701.
26. U. R. Pedersen, T. B. Schrøder, J. C. Dyre and P. Harrowell, Phys. Rev. Lett., 2010, 104, 1–4.
27. D. Coslovich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 8.
28. A. Malins, J. Eggers, C. P. Royall, S. R. Williams and H. Tanaka, J. Chem. Phys., 2013, 138, 12A535.
29. T. Schenk, D. Holland-Moritz, V. Simonet, R. Bellissent and D. M. Herlach, Phys. Rev. Lett., 2002, 89, 75507.
30. D. B. Miracle, Nat. Mater., 2004, 3, 697–702.
31. M. Tanemura, Y. Hiwatari, H. Matsuda, T. Ogawa, N. Ogita and A. Ueda, Prog. Theor. Phys., 1977, 58, 1079–1095.
32. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem., 1987, 91, 4950–4963.
33. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B, 1983, 28, 784–805.
34. S. R. Williams, arXiv:0705.0203, 2007.
35. C. P. Royall, S. R. Williams, T. Ohtsuka and H. Tanaka, Nat. Mater., 2008, 7, 556–61.
36. A. Malins, J. Eggers and C. P. Royall, 2013, in preparation.
37. R. Della Valle, D. Gazzillo, R. Frattini and G. Pastore, Phys. Rev. B: Condens. Matter, 1994, 49, 12625–12632.
38. G. Biroli, J. P. Bouchaud, A. Cavagna, T. S. Grigera and P. Verrochio, Nat. Phys., 2008, 4, 771–775.
39. M. Mosayebi, E. Del Gado, P. Ilg and H. C. Öttinger, Phys. Rev. Lett., 2010, 104, 205704.
40. F. Sausset and D. Levine, Phys. Rev. Lett., 2011, 107, 045501.
41. W. Kob, S. Roldán-Vargas and L. Berthier, Nat. Phys., 2011, 8, 164–167.
42. C. Cammarota and G. Biroli, Europhys. Lett., 2012, 98, 36005.
43. G. M. Hocky, T. E. Markland and D. R. Reichman, Phys. Rev. Lett., 2012, 108, 225506.
44. A. J. Dunleavy, K. Wiesner and R. C. P., Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 041505.
45. M. Leocmach, R. J and H. Tanaka, J. Chem. Phys., 2013, 138, 12A536.
46. D. Coslovich, J. Chem. Phys., 2013, 138, 12A539.
47. A. Montanari and G. Semerjian, J. Stat. Phys., 2006, 125, 23–54.
48. B. Charbonneau, P. Charbonneau and G. Tarjus, Phys. Rev. Lett., 2012, 108, 4.
49. B. Charbonneau, P. Charbonneau and G. Tarjus, J. Chem. Phys., 2013, 138, 12A515.
50. P. Charbonneau and G. Tarjus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042305.
51. A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett., 2006, 96, 185701.
52. M. Mosayebi, E. Del Gado, P. Ilg and H. C. Öttinger, J. Chem. Phys., 2012, 137, 024504.
53. W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 4626–4641.
54. S. Karmakar and I. Procaccia, arXiv:1204.6634, 2012, p. 6.
55. A. Malins, PhD thesis, University of Bristol, 2013.
56. J. P. K. Doye, D. J. Wales and R. S. Berry, J. Chem. Phys., 1995, 103, 4234–4249.
57. J. P. K. Doye and D. J. Wales, J. Chem. Soc., Faraday Trans., 1997, 93, 4233–4243.
58. S. Toxvaerd, T. B. Schrøder and J. C. Dyre, arXiv:0712.0377, 2007.
59. J. R. Fernandez and P. Harrowell, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2003, 67, 011403.
60. T. Middleton, J. Hernandez-Rojas, P. N. Mortenson and D. J. Wales, Phys. Rev. B: Condens. Matter, 2001, 64, 184201.
61. D. Chandler and J. P. Garrahan, Annu. Rev. Phys. Chem., 2010, 61, 191–217.
62. T. Speck, A. Malins and C. P. Royall, Phys. Rev. Lett., 2012, 109, 195703.
63. S. Stoddard and J. Ford, Phys. Rev. A: At., Mol., Opt. Phys., 1973, 8, 1504.
64. S. Nose, J. Phys. Soc. Jpn., 2001, 70, 75.
65. M. Goldstein, J. Chem. Phys., 1969, 51, 3728–3739.
66. L. Berthier and J. P. Garrahan, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2003, 68, 041201.
67. H. Vogel, Phys. Z., 1921, 22, 645.
68. G. S. Fulcher, J. Am. Ceram. Soc., 1925, 8, 339–355.
69. G. Tammann and W. Hesse, Z. Anorg. Allg. Chem., 1926, 156, 245–257.
70. V. K. De Souza and D. J. Wales, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 134202.
71. H. Tanaka, J. Non-Cryst. Solids, 2005, 351, 3385–3395.
72. T. Hecksler, A. I. Nielsen, N. Boye Olsen and J. C. Dyre, Nat. Phys., 2008, 4, 737–741.
73. J. L. Meijering, Philips Res. Rep., 1953, 8, 270.
74. W. Brostow, J.-P. Dussault and B. L. Fox, J. Comput. Phys., 1978, 29, 81–92.
75. N. Medvedev, J. Comput. Phys., 1986, 67, 223–229.
76. E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates and P. N. Pusey, Phys. Rev. Lett., 2009, 103, 135704.
77. I. Saika-Voivod, R. Bowles and P. H. Poole, Phys. Rev. Lett., 2009, 103, 225701.
78. E. Sanz, C. Valeriani, E. Zaccarelli, W. C. K. Poon, P. N. Pusey and M. E. Cates, Phys. Rev. Lett., 2011, 106, 215701.
79. J. Taffs, S. R. Williams, H. Tanaka and C. P. Royall, Soft Matter, 2013, 9, 297–305.
80. N. Lačević, F. W. Starr, T. B. Schrųder and S. C. Glotzer, J. Chem. Phys., 2003, 119, 7372.
81. E. Flenner and G. Szamel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 051502.
82. E. Flenner and G. Szamel, J. Phys.: Condens. Matter, 2007, 19, 205125.
83. W. Kob and H. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 4134–4153.
84. H. Shiba, T. Kawasaki and A. Onuki, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 041504.
85. T. Kawasaki and H. Tanaka, J. Phys.: Condens. Matter, 2010, 22, 232102.
86. W.-S. Xu, Z.-Y. Sun and L.-J. An, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 041506.
87. A. Onuki, Phase Transition Dynamics, Cambridge University Press, Cambridge, 2002.
88. S. Karmakar, C. Dasgupta and S. Sastry, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 3675.
89. C. Cammarota and G. Biroli, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 8850–5.
90. T. Kawasaki, T. Araki and H. Tanaka, Phys. Rev. Lett., 2007, 99, 2–5.
91. E. Flenner and G. Szamel, J. Chem. Phys., 2013, 138, 12A523.
92. P. Harrowell, Dynamical Heterogeneities in Glasses, Colloids and Granular Materials, OUP, Oxford, 2010.
93. L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. El Masri, D. L'Hote, F. Ladieu and M. Pierno, Science, 2005, 310, 1797–1800.
94. D. Fragiadakis, R. Casalini and R. C. M., Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 042501.
95. A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2009, 103, 135703.
96. A. Furukawa and H. Tanaka, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 061503.
97. A. Furukawa and H. Tanaka, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 030501.
98. S. Singh, M. D. Ediger and J. J. de Pablo, Nat. Mater., 2013, 12, 139–44.

---