Structural-dynamical transition in the Wahnstr\"om mixture

In trajectory space, dynamical heterogeneities in glass-forming liquids correspond to the emergence of a dynamical phase transition between an active phase poor in local structure and an inactive phase which is rich in local structure. We support this scenario with the study of a model additive mixture of Lennard-Jones particles, quantifying how the choice of the relevant structural and dynamical observable affects the transition in trajectory space. We find that the low mobility, structure-rich phase is dominated by icosahedral order. Applying a nonequilibrium rheological protocol, we connect local order to the emergence of mechanical rigidity.


Introduction
Supercooled liquids show emergent dynamical and structural heterogeneities when cooled towards the glass transition [1,2,3,4]. The relation between slow dynamics and some form of short-range (local) order, however, is still poorly understood. On the one hand, the efficient filling of space with atoms of different sizes requires a certain degree of topological order [5] and the dynamic slowdown can rigorously be linked to emerging static lengthscales [6]; on the other hand, computer simulations have shown that the correlation between local structural features and slow dynamics is strongly model dependent [7,8]. In experiments, colloidal [9,10,11,12] and metallic glasses [13,14,15] provide evidence for emerging local order as well as, on the contrary, support for purely dynamical scenarios where local structure has limited influence on the dynamics [16,17]. Historically, the study of local structure with complex higher order metrics has played a decisive role in understanding amorphous systems and packings since the times of Bernal and Finney [18,19,20,21] and has contributed to a geometric and thermodynamic interpretation of the emerging frustration in glasses [22,23]. However, alternative approaches which disregard structural features and focus on dynamical [24] or vibrational/elastic aspects [25] of the problem or relaxation have been proposed, in striking contrast with established thermodynamic theories of the glass transition [26,27]. It is therefore important to understand what drives strong or weak coupling between structure and dynamics in different supercooled liquids.
A major difficulty encountered in the investigation of the role of structural changes in dynamic arrest is the fact that particle-resolved studies (and in particular conventional computer simulations) can only access a limited dynamic range of slow relaxation. Typically, this encompasses 4 to 5 orders of magnitude in time, meaning that such studies mainly capture the onset of the mechanisms that characterise the deeply supercooled and glassy regimes (when the relaxation times are 10 to 20 orders of magnitude larger with respect to the liquid regime) [28]. Therefore, alternative sampling routes to explore the deeply supercooled regime from a structural and/or dynamical point of view have been developed in recent years, including pinning fields [29,30,31], particle-swap Monte-Carlo on particular models [32,33] or biased dynamical ensembles [24,34,35,36].
A potential route to study dynamical and structural heterogeneities in glassformers is provided by efficient sampling methods in trajectory space, where novel dynamical phase transitions have been uncovered and connected to the dynamical slowdown observed in supercooled liquids [24]. The study of trajectory space in glassy systems has been originally promoted in the context of the dynamical facilitation theory of slow dynamics [37,38,34]. Within this framework, on-lattice idealised models [37,39,40,41] as well as more realistic models of structural glasses [34,35,42,43,44] have been shown to undergo a first-order dynamical phase transition in trajectory space between an active phase with high mobility (fast relaxation) and an inactive phase with low mobility (slow dynamics).
However, this purely dynamical picture has been more recently complemented by a structural aspect: active/inactive phases correspond to trajectories particularly poor/rich in local structure [45,46] and can be seen as representative of the low temperature state of the supercooled liquid [36]. Dynamical transitions are therefore understood to correspond to a structural-dynamical transitions, where the slowdown of the dynamics becomes intimately related to the growth of short-range-order domains.
Still, much of the evidence for structural-dynamical phase transitions in atomistic models of glassformers up to now is restricted to only two model systems (the Kob-Andersen mixture [34,42,45,35,36], a popular Lennard-Jones mixture with weak structural dynamical correlations [47], and the moderately polydisperse hard spheres [46]). In order to understand how system-dependent this picture is, it is important to extend the scope of these studies to other model systems.
In the present numerical work, we consider the case of a popular atomistic glassformer originally introduced by Göran Wahnström as a simple model for supercooled liquids [48]. It consists in a binary mixture of Lennard-Jones particles whose parametrization has been found to provide a good model of fragile glasses, with a particularly strong coupling between its slow dynamics and the emergence of local geometrical motifs [47,49,50]. These are typically icosahedra, a very common arrangement in simple models of glass-forming liquids composed of spherically symmetric particles.
The article is structured as follows: in Section 2 we present the model studied and the importance sampling technique employed for trajectory sampling; in Section 3 we introduce the relevant observables and the phase transitions in trajectory space that can be probed through the dynamical s-ensemble and the structural-dynamical µ-ensemble; in Section 4 we show that it is possible to connect the structural-dynamical transition to the emergence of rigidity in the glass, as the icosahedra-rich phase presents distinctive rheological properties; finally, we conclude the article with an overview of the results and their implications.

Model and sampling technique 2.1 The Wahnström binary mixture
We study the Wahnström binary mixture of Lennard-Jones particles. The model is a 50:50 mixture of large (A) and small (B) particles with parameters σ A /σ B = 1.2, m A /m B = 2, A / B = 1 and cutoff r cut = 2.5σ at number density ρ = 1.296. Lengths, temperature and times are reported in units of σ A , A /k B , and (m A σ 2 A / A ) 1/2 , respectively. The mixing rule for the interaction is additive, i.e. it follows the Lorentz-Berthelot rules This atomistic supercooled liquid has been extensively studied since its original design [48]. The model reproduces to a good degree the relaxation behaviour of socalled fragile glasses, as its structural relaxation time τ α (as measured from the decay of the intermediate scattering function [49]) undergoes a non-Arrhenius (superexponential) increase when the system is cooled below the crossover or onset temperature T onset = 1.0 [51,49]. Furthermore, as the temperature is decreased, the disordered structure of the liquid changes with the formation of five-fold symmetric domains and in particular of local particle motifs with icosahedral coordination [49,51,52] which contribute to the emergence of strong frustration [53]. Equilibration of the liquid in conventional simulations around and below the so-called mode-coupling temperature T M C = 0.56 is computationally expensive, making the low temperature, activated regime (crucial for testing theoretical predictions [2]) unreachable. Divergence of the relaxation times, if modelled by the super-Arrhenius Vogel-Fulcher-Tamman law ln τ α /τ ∞ = DT 0 /(T − T 0 ), is predicted at temperature T 0 ≈ 0.46. For reference, we report in Fig. 1(a) the temperature dependence of both the structural and dynamical properties of the model.
Beyond local structural order, the model has been shown to crystallise, under suitable conditions, into a MgZn 2 Laves phase formed by icosahedral motifs and socalled Frank-Kasper bonds [54] but in the supercooled liquid regime the contribution of such a large unit cell to the increased degree of local order has been shown to be limited [49].

Replica exchange in trajectory space
As in previous work [35,36,45,46], in order to sample large fluctuations of the time-integrated observables we employ an importance sampling technique that extends equilibrium replica exchange methods to ensembles of trajectories.
We sample space and time extensive observables O x on systems of N = 512 evolving for a finite observation time t obs . A generic time-integrated observable O x is defined as a double sum over the number of particles and a discretization of time into L intervals for a total of K = L+1 points: where f t,i is a specific microscopic observable (e.g. a single particle indicator function). The goal of the importance sampling technique is to efficiently measure the probability distribution P (O x ; T s ) for a given value of the thermostat temperature T s . In particular, we are interested in the large deviations from the typical value of the probability distribution. In order to calculate such rare fluctuations in trajectory space, new trajectories are generated through shifting and shooting moves (inspired by Transition Path Sampling [55]). Hence, the algorithm performs a random walk in trajectory space with acceptance probability determined by a Metropolis rule which ensures detailed balance and where are the values of a biasing pseudopotential which is a function of the extensive observable O x computed over old and new trajectories. We choose Ψ to have a parabolic form where O j 0 is the reference (typical) value associated to replica j. Depending on the observable, we take a number of distinct replicas varying from 8 to 16, with equally spaced values of O j and values for the harmonic constant ω that ensure good mixing of neighboring replicas. Mixing is also enhanced by 2500 swap attempts among all (not-necessarily neighboring) replicas.
The Monte-Carlo algorithm in trajectory space simulation starts with an equilibrated trajectory assigned to all replicas at temperature T s . A new trajectory is then generated via Transition Path Sampling moves (1/4 shifting, 3/4 shooting [55]) independently for every replica, accepted or rejected according to Eq. 3. Swap attempts between different replicas are then performed, completing the cycle. During the sampling, we employ a velocity-Verlet integrator with timestep dt = 0.005 to resolve the equation of motion and the Andersen thermostat to keep the temperature constant.
We perform several tens of thousands of cycles and collect statistics and block averages from three to eight non-overlapping blocks of data whose size ranges between 1.2 · 10 4 and 3 · 10 4 trajectories, which deliver estimates for the averages and standard errors. Crucially, depending on the sampling temperature T s , correlations may be very long-lived and the number of Monte-Carlo cycles spent during the equilibration in trajectory space can be very large (∼ 6 · 10 4 ) Monte-Carlo sweeps), as shown in Fig. 2. We then discard trajectories produced during equilibration and collect data from the converged, late steps of the Monte-Carlo.
From the collected ensemble of trajectories, we calculate distributions and expectation values using the Multistate Bennett Acceptance ratio (MBAR) method extended to ensembles of trajectories [57]. This technique allows us to obtain the unbiased probability distribution P (O x ; T s ) and expectation values for any quantity A as where y is the conjugated field to the observable O x and indicate averages according to the unbiased distribution P (O x ; T s ). Notice that the denominator in Eq. 5 corresponds to the moment generating function of the probability distribution and is a generalization of the partition sum to trajectory spaces. In this work, we focus on two particular ensembles:  ) γ ] −1 fit for τα and n respectively, with γ = 6.6, T 1/2 = 0.47 from Ref. [56]. Vertical lines also indicate the location of relevant temperatures in the Wahnström model: the onset of slow dynamics Tonset, the mode-coupling transition temperature TMC and the VFT temperature T0. (b) Threedimensional rendering of the four local motifs considered in this work: the icosahedron, the defective icosahedron (10B) and two nine-particle motifs unrelated to five-fold symmetry.
For the icosahedron and 10B, a pentagonal ring of particles is highlighted in gold. and the observable O x is the time-integrated mobility of the particles; and the µ-ensemble, where y = µ and the relevant observable is a time-integral over the number of particles in a particular local motif (here the icosahedron).
In the presence of transitions in trajectory space, we expect to measure probability distributions for the timeintegrated observables that are not Gaussian and display long, eventually exponential tails. For variables that follow a Gaussian probability distribution, the kurtosis the ratio between the fourth central moment and the squared second moment) has value κ 4 = 3. Therefore, the excess kurtosis κ exc = κ 4 −3 is often employed as a benchmark for the deviations from a Gaussian distribution. So-called leptokurtic (fat-tailed) distributions correspond to positive κ exc while platykurtic (thin-tailed) distributions correspond to negative κ exc .  Table 1. Expectation values for the mean, the variance and the excess kurtosis κexc for several observables as measured via trajectory sampling for different values of the trajectory length at temperature Ts = 0.67. 3 Dynamical and structural phase transitions

Observables
We analyse the emergence of phase transitions in trajectory space by monitoring a variety of observables. We perform importance sampling in trajectory space according to time-integrated observables that are either dynamical (such as the mobility excitations) or structural (a selection of geometrically different structural motifs, see Fig. 1(b)). Furthermore, in order to relate the trajectory-space picture back to the thermodynamic picture, we also monitor the inherent state energy of the selected configurations, whose statistics in the trajectory ensemble has been proven to closely reproduce the equilibrium properties.
Structures are detected employing the Topological Cluster Classification algorithm and we refer to Reference [58] for a more detailed discussion of the geometries considered here.
In particular, for the time-integrated quantities we have: number of excitations: To quantify the number of mobile particles, we compute the observable where is the single particle displacement, Θ is the Heaviside function and a is a scale for cage motion, here set to a = 0.3σ large . number of particles in icosahedral motifs: Given the important role of icosahedral order in the Wahnström mixture, we track this specific local motif along the trajectories. Additionally, we perform importance sampling according to the number of icosahedra. The corresponding time-integrated extensive structural-dynamical observable is then where h ico i is an indicator function, which takes value 1 if a particle is found in an icosahedral environment or 0 if it is not. With a certain abuse of language, we will interchangeably refer to the population of icosahedra or the population of particles in icosahedral motifs when considering the intensive quantity O ico /N K.
number of particles in 9A motifs: We compute O 9A performing the summation as in Eq.7, but with a different indicator function h 9A i . In this case, we consider the 9A structure of the Topological Cluster Classification, which is composed of six particles combined to form three four-folded rings, surrounded by three further spindle particles on each quadrangular facet (forming a tricapped trigonal prism). According to previous studies [49], we do not expect this motif to be a good predictor of structural-dynamical heterogeneity for the Wahnström mixture. However, in the case of other simple liquids dominated by five-fold symmetric local order, such as moderately polydisperse hardsphere, 9A motifs have been shown to be complementary to local icosahedral order, becoming less frequent when the packing fraction (and the population of icosahedra) increase [59].
number of particles in BCC motifs: As a further test, we compute the time-integrated observable O BCC considering a nine particle structure that (weakly) correlates with body centered cubic local order and anticorrelates strongly with icosahedral and five-fold symmetric order. number of particles in five-fold symmetric motifs: Finally, to track five-fold symmetric local order that is not fully icosahedral, we consider the defective icosahedron structure 10B, composed of three interlaced pentagonal rings. This structure is characteristic of hardsphere mixtures, and has been shown both in simulations and experiments to drive a clear structuraldynamical phase transition [46].
We also measure a static observable, i.e. not timeintegrated. This is the inherent state energy (ISE) of configurations located at the centre of each trajectory, chosen in order to avoid finite-time effects on the statistics [60]. Inherent state energies are obtained minimising the potential energy of the system for a maximum of 1000 iterations of the FIRE algorithm [61].

s-ensemble
First, we consider the response of the system to a dynamical bias. This means that we collect trajectories according the observable O exc , i.e. the time-integrated number of mobility excitations. We employ the large deviation formalism and notation, and we define s as the dynamical conjugate field related to the excitations, so that positive/negative values of s correspond to atypically small/large densities of mobility excitations, hence the name of s-ensemble [37]. As we sample the mobility large deviations, we track all the other dynamical and static order parameters.
In Fig. 3 and in Table 1 we summarise our findings for a particular thermostat temperature T s = 0.67T onset ≈ 1.2T MC , where T MC and T onset indicate the transition temperature predicted by the power-law fit to the relaxation of mode-coupling theory and the onset of the two-step relaxation dynamics respectively. In Fig. 3 we compare the (scaled) logarithm of the probability distributions (i.e. the rate function) of the considered observables O x for increasing values of the trajectory length t obs = 1.9, 3.8, 5.7, 9.5τ α (K = 20, 40, 60, 100). At the considered temperature, we expect to observe deviations from Gaussian fluctuations in the tails (i.e. large deviations) of the probability distributions. With this comparative analysis, we want to stress that the choice of the observable is non-trivial and different observables present characteristic features.
First we notice that the population of excitations (which is the reaction coordinate along which we perform importance sampling) shows mostly Gaussian fluctuations around the mean value O x /KN = 0.096(2) for all the sampled trajectory lengths. However, the variance computed at different trajectory lengths appears to slowly converge to smaller values, with the tails of the probability distributions gradually narrowing. This indicates that very short trajectories of length t obs = 1.9, K = 20 are affected by finite size effects that enhance the observation of large fluctuations.
Higher order moments converge even more slowly but point to the emergence of non-Gaussian features. For example, the excess kurtosis is negative for short trajectories and becomes mildly positive for the longest trajectories K = 100. This underlines that even longer trajectories are needed to obtain more marked signatures of a dynamical phase transition in terms of population of excitations at the relatively high temperature considered here, with an non-negligible increase of the computational cost. Notice that it is only in the long time limit that a large deviation principle holds and rate functions converge [62], and therefore it is only in this limit that a formal phase transition in trajectory space is expected. measured on the same trajectories produced in the sensemble? In the following, we analyse them one by one.
For the time-integrated population of particles in icosahedral motifs, we observe that average values do not depend on K; however, higher order moments show a dependence on the trajectory length. The values of the excess kurtosis κ exc show a marked increase in non-Gaussian features of the trajectory probability distribution, as confirmed by direct inspection of the probability distribution. The excess kurtosis is positive (i.e. fat-tails) and goes approximately from 0.48 to 2.0 when the trajectory length increases from K=20 to K=100. For a comparison, notice that for a common leptokurtic distribution of positive random variables such as the Rayleigh distribution, the excess kurtosis is κ exc = −(6π 2 − 24π + 16)/(4 − π) 2 hence κ exc ≈ 0.24 showing that the distribution for the icosahedra is even more leptokurtic. Compared to the response of the mobility excitations, the time-integrated population of icosahedra provide a much stronger signature for a dynamical phase transition. In particular we observe that populations of icosahedra of order 0.2 are only two orders of magnitude less likely than the converged typical value O ico /N K = 0.09, with a strong exponential tail in the probability distribution. Non-Gaussian fluctuations are therefore stronger when tracking the time-integrated population of icosahedra than in the case of excitations.
These results are consistent with previous literature [49,53,47] where the role of icosahedral motifs as locally favoured structures (LFS) of the Wahnström mixture as been discussed and their strong correlation with dynamical heterogeneities measured. They also confirm the scenario originally suggested for another popular glassformer (the Kob-Andersen mixture), whereby trajectories sampled according to time-integrals of the LFS delivered stronger signatures for a dynamical transition that mobility excitations [35,36].
An icosahedral motif is detected in the TCC via the combination of seven five-fold symmetric rings [58], and the statistics of the number of icosahedra appears to strongly indicate the presence of non-Gaussian fluctuations related to a structural-dynamical phase transition in the system. How does such transition change if we take into account a less restrictive observable that still identifies five-fold symmetry? To answer this question, we consider the so-called defective icosahedron structure 10B (see Sec.3.1 above). We first notice that the average population of particles in 10B per trajectory is much larger than the population of icosahedra (0.59 vs 0.089) and the variance again slowly converges with increasing t obs . However, The excess kurtosis is much smaller in absolute values, changing sign from negative towards positive values (leptokurtic distributions) as the trajectory length is increased. This matches the dynamical notion of locally favoured structures: icosahedra are not only the minimum energy structure for the Wahnström interaction, they also are the individual motif (among the several options of the Topological Cluster Classification) that displays the longest persistence time [49]. The indicators for a structural-dynamical transition in terms of 10B motifs are much weaker than in the case of icosahedral order. Yet, they confirm that the inactive (low population of excitations) regime is dominated by long-lived five-fold symmetric motifs.
Is it possible to detect signs of the transition in other structural observables? We consider the two exemplary cases of the 9A and BCC9 structures. These motifs both correspond to arrangements of 9 particles with different symmetries which are not minimum energy clusters of the potential. The average populations of the two motifs are very different (∼ 0.07 for 9A and 0.74 for BCC9). The 9A probability distribution is well approximated by a Gaussian for all the trajectory lengths considered here, and the corresponding excess kurtosis are (in absolute value) the smallest among all the considered structures. The BCC9 motif, conversely, presents relatively large but negative excess kurtosis, indicating that the tails of the distributions decay more rapidly than in the case of a Gaussian distribution.
For a given trajectory length, we consider the sensemble averages as a function of the field s to highlight correlations and anticorrelations between the observables. With trajectories of length K = 60, we show in Fig. 4 that for negative s 0 we sample trajectories characterised by large densities of excitations (active phase) while for s 0 we have trajectories with low densities of excitations (inactive phase). These correspond respectively to trajectories that are poor and rich in icosahedra. The anticorrelation between mobility and five-fold symmetry is reflected also in the negative correlation between mobility and 10B structures. On the other hand, mobility positively correlates with the remaining motifs (9A and BCC9).
Finally, we consider how the active/inactive transition is translated in terms of the energy landscape of the system. To do so we also track the inherent state energy (ISE) of the central configuration of every single trajectory and plot the corresponding probability distribution. This (as expected) does not show dependence in the trajectory length and it is well reproduced by a Gaussian fit, see Fig. 3. Normal fluctuations are confirmed by the analysis of the respective excess kurtosis, which are by far the smallest measured throughout our analysis (as small as κ exc = −0.003). In Fig. 4, we do observe a transition to trajectories whose central configurations display typically much more negative energies with respect to the equilibrium typical value at s = 0. This is consistent with the finding that in a different binary mixture (Kob-Andersen) low mobility is a good predictor of low inherent state energies [60].

µ-ensemble
The direct route to access structural-dynamical phase transitions is to sample trajectories according to a relevant time-integrated structural observable. From the previous discussion, and in particular from the magnitude of the non-Gaussian fluctuations as measured by the excess kurtosis, it is evident that icosahedral motifs are well suited to this purpose.
Therefore we perform additional trajectory sampling according to the time-integrated number of icosahedral motifs. As in the case of the s-ensemble, we sample trajectories following the replica exchange scheme, with quadratic pseudo-potentials for the replicas with suitable spring constant ω.
In the new ensemble of trajectories, the conjugate field related to the number of particles in icosahedral motifs is termed µ. Consistently with previous works in the literature [35,36], averages of any arbitrary quantity in the µ-ensemble are defined as In the previous section, we have shown that in the s-ensemble an emergent active/inactive transition is mirrored by a rapid increase of the population of particles in icosahedral motifs. In the µ-ensemble we sample such structural transition directly. In Fig. 5(a,b) we plot the µdependence of the average mobility O exc µ /N K and the average population of icosahedra O ico µ /N K for several thermostat temperatures T s , from T s = 0.72 to T s = 0.65. At different temperatures, we perform simulations of different trajectory lengths t obs . Since the relevant timescale for the dynamics is the structural relaxation time τ α , we plot the first moments as a function of the nondimensional scaled conjugate field µτ obs := µt obs /τ α . Just below the onset temperature we observe signs of a phase transition at large µτ obs between trajectories poor in icosahedra with high mobility and trajectories rich in icosahedra with low mobility. As we reduce the thermostat temperature, the transition moves to values closer to µ = 0. Through a spline fit and the estimate of the maximum in the derivative, we obtain the value µ * τ obs at which the transition takes place. The very small values of µ * τ obs at relatively high temperatures compared to T 0 obtained from the Vogel-Fulcher-Tammann fit or the mode-coupling T M C temperatures suggest that trajectories with an exceptionally high population of icosahedra should be highly likely, and signatures of bi-modality in the probability distribution of the time-integrated observables should become accessible even to conventional simulations as the temperature is reduced.
In Fig. 5 (c,d) we plot such probability distributions both for µ = 0 and the critical value µ = µτ * obs , shifting and rescaling the abscissa axis by the mean and the standard deviation. We observe that, as temperature decreases, the structure-rich tail of the probability distribu- tions raises of several orders of magnitude. Signs of bimodality are weak at low temperatures, due to the relatively short observation time t obs , but clearer at higher temperatures. Moreover, if we evaluate the probability distributions at coexistence µ = µ * , Fig.5(d), a peak at high population of structures emerges more clearly. The knowledge of µ * τ obs (T s ) allows us to draw an approximate structural-dynamical phase diagram, Fig. 6, identifying the locus of points where the transition from icosahedra-poor to icosahedra-rich trajectories occurs. For the considered temperatures, we observe that most of the data points lie on a straight line. An extrapolation of the line to µτ obs = 0 would imply that at temperature T = 0.64 coexistence between the two structuraldynamical phases would be observable at µ = 0, i.e. in conventional simulations with no need for importance sampling. Previous numerical studies of the model [63,49] managed to equilibrate the supercooled liquid down to temperature T = 0.58, with no signature of a transition while decreasing the temperature, but with a rapid increase of the population of icosahedra. This excludes the possibility of a transition at µ = 0 for at least T > 0.58. As discussed in [36], several alternative scenarios can be obtained with different extrapolations at low temperatures, including ones where the transition asymptotically reaches µ = 0 only in the T → 0 limit [64]. Here we notice that as we reduce the temperature, the critical field µ * τ obs is reduced by progressively smaller amounts for the successive temperatures. Lower temperature sampling is partly hindered by the long convergence times of the Monte-Carlo in trajectory space, see Fig. 2.
In the icosahedra-rich regime, approximately 50% of the particles can be found in a local icosahedral environment. However, a complex unit cell formed by several icosahedra and Frank-Kasper bonds has been shown to drive the system towards crystallisation [54]. We check this possibility by monitoring the concentration of Frank-Kasper bonds, here defined as pairs of large A particles surrounded by six common B particles. In Fig. 7 we plot the average fraction of particles involved in Frank-Kasper bonds for the increasing reference concentration of icosahedra O j ico in the replica-exchange scheme at an exemplary temperature T s = 0.67. We observe a rapid increase in the number of Frank-Kasper bonds as we consider replicas with very high concentrations of icosahedra. This is consistent with the overall behaviour of the Wahnström supercooled liquid at low temperatures, where Frank-Kasper bonds are very common [54]. However, in order to form a crystalline phase, four-fold Frank-Kasper bonds between the large particle species are necessary. If we focus on the fraction of A particles in four-fold bonds, this increases very mildly across all of the replicas, and stays below 5% in the highest bias replica, excluding crystal formation in the icosahedra-rich phase.
In conclusion, both the s and the µ ensemble calculations provide evidence for an inactive and icosahedra-rich dynamical phase that becomes progressively more likely to be observed for T < T onset . We now study the icosahedral phase more in detail to understand its relation with the emergence of rigidity in the glass.  many orders of magnitude. Such a phenomenological transition is accompanied by the emergence of solid behaviour: the glass behaves like a solid, in the sense that it can be probed through rheological measurements, proving a finite elastic response and shear modulus.

Rheological response of the inactive/icosahedra-rich phase
We have shown that as the temperature is decreased, the Wahnström mixture explores more and more frequently trajectories that are exceptionally rich in structure. Moreover, the icosahedra-rich trajectories not only are characterised by low mobility (inactive trajectories) but they also tend to have configurations with low inherent state energies. Is it possible to connect these structural and dynamical changes to the emergence of solidity, i.e. to the rheological response of the system?
We test this idea realising an ensemble of configurations extracted from the trajectories produced in the µ ensemble at the thermostat temperature T s = 0.65. From every umbrella i of the replica-exchange algorithm we extract a population of configurations that are representative of the fluctuations, in trajectory space, around a specific value of the population of icosahedra In particular, we produce a discrete group of 8 sets with 75 initial configurations each at the following typical population of icosahedra n 0 To understand the purely mechanical response of the different sets of configurations we study the linear shearing of the system in the Athermal Quasi-static limit (AQS) [65,66]. Under this protocol, the system is slowly deformed in a chosen direction at a fixed shear rateγ = 0.005 for a small time interval ∆t = 0.005; subsequently, the FIRE energy minimisation algorithm [61] is employed to lead the particles to the closest inherent state. The two steps are repeated until the system reaches a maximum total strain of γ = 0.5 In Fig. 8 and Fig. 9 we plot the response of the system in terms of shear stress σ xy and fraction of particles in icosahedral domains for different typical values of the initial population of icosahedra n 0 i . The first striking result is that the yield stress σ yield = max γ σ xy (γ) strongly depends on n 0 i , and it approximately doubles as the typical population of icosahedra quadruples. The yield strain γ yield (the value of strain at which the maximum stress is reached) is not sensitive to the different starting conditions and is located at approximately γ yield ≈ 0.12 for the chosen strain rate. At the same time, we notice that the shear protocol induces a sudden increase of the population of icosahedra at very early times (very small strains) and a progressive decay of the population which accelerates as the yield strain is reached. The overall, instantaneous increase of the population of icosahedra can be  understood as a consequence of the minimisation procedure, which destroys thermal fluctuations present in the initial configurations and promotes the formation of local minimum energy motifs, such as the icosahedron. This implies that the overall population of icosahedra n ico (γ) can be distinguished into two families: the first refers to the subset of particles that are located in icosahedral domains in the original starting configurations produced in the µ ensemble, and it is identified by the boolean vector η µ of length N ; the second refers to all the remaining particles in icosahedral domains, resulting from the AQS protocol, identified by the vector ημ.
As the system is sheared, the number of icosahedra changes very mildly for strains below the yield strain, and only later declines, supporting the idea that the popula-b a Fig. 10. Auto-correlation of the probability for a particle to be in icosahedral domain for the ηµ and the ημ populations (see main text for definition) for different values of the typical initial population of icosahedra.The ημ family is not defined at γ = 0 so, we take the smallest γ as the reference state. The vertical dashed line corresponds to the yield strain.
tion of icosahedra is related to the rigid, elastic response of the system. Having defined two subpopulations of icosahedra, we now quantify their respective differences in the mechanical response.
To do so, we compute separate auto-correlation functions η x (γ)η x (0) for the η µ and the ημ populations, Fig. 10(a,b). We notice that only for large initial populations of icosahedra the autocorrelation functions start close to unity. This shows that the reorganisation induced by the AQS protocol not only forms new icosahedral motifs, but it also initially destroys a fraction of them. The two families of autocorrelation functions show distinctively different behaviours: the icosahedra present in the initial µ-ensemble, Fig. 10(a), show a long plateau that terminates only when the yield strain is attained; the icosahedra generated via AQS, Fig. 10(b), continuously decorrelate at earlier times (smaller strains).
A further confirmation of the different responses between the η µ and the ημ icosahedra is provided by the distribution of the potential energy of the individual particles constituting the two families. In Fig. 11 we plot the overall energy distributions, for the population, collected all along the shearing protocol. We clearly observe that not all icosahedral motifs are energetically the same: particles located in icosahedral motifs purely emerging from energy minimisation have energies that are typically higher than particles identified in icosahedral motifs in the original µ ensemble configurations. The energy gap between the two families widens as we consider initial configurations with larger concentrations of icosahedra: at very high initial concentrations (55%), the η µ subpopulation has energies that are 8% lower than the ημ subpopulation.

Conclusions
Through numerical simulations, we have discussed a third example of structural-dynamical phase transition in a model of atomistic glassformer, after the previously considered cases of the Kob-Andersen mixture [34,42,45,35,36] and the moderately polydisperse hard-spheres [46].
A quantitative analysis of the probability distributions of time-integrated observables demonstrates that wellchosen time-integrated structural motifs can be used to perform efficient importance sampling. In particular, it makes possible to explore structure-rich trajectories (representative of colder temperature states [36]) that are otherwise hard to reach. At the same time, we find confirmation for a sharp (first-order) transition in trajectory space, that becomes measurable below the onset temperature, between a structure-rich and a structure poor dynamical phase, the former becoming more and more likely as the temperature is reduced, similarly to what has been previously observed in the Kob-Andersen mixture [36]. Within the range of the explored temperatures, it is impossible to assess what the low temperature fate of the transition may be: the reduction of the critical conjugate field value µ * suggests that, as the temperature decreases, the structure-rich phase would prevail. However, it is unclear whether the previously reported crystallisation into complex Laves phases [54] would interfere with the emergence of the icosahedra-rich phase. In our simulations, we monitored the evolution Frank-Kasper bonds (an essential element of the complex crystalline phase) and do not find a significative increase in the icosahedra-rich phase compared to the icosahedra-poor phase.
In Reference [36], the study of the alternative Kob-Andersen mixture at low temperature indicated possible scenarios for the temperature dependence of the structural-dynamical transition. In the present case of the Wahnström mixture, we observe that at 0.72T onset the structure-rich phase is highly metastable, while a relatively modest decrease of the temperature to 0.67T onset makes the exploration of the structure-rich basin 5 to 7 times more likely, see Fig. 5 (d). The temperature dependence of the critical value µ * τ obs shows a decrease towards µ = 0 which becomes less pronounced as the temperature is decreased. This is accompanied by strong correlations between successive steps in the trajectory-space Monte-Carlo that slow down equilibration and make lower temperature sampling particularly challenging. The present data support the narrowing of the free-energy gap (in trajectory space) between the structure-rich and structurepoor states when the temperature is decreased, and do not exclude the possibility that the transition terminates at a lower critical point at finite temperature, as in kinetically constrained models with additional softness [67].
In order to better understand the importance of the structure-rich phase, configurations obtained through trajectory sampling have been probed through an outof-equilibrium rheological protocol, effectively regarding these configurations as samples of an amorphous material at T = 0. Consistently with previous studies of the Wahnström mixture based on conventional simulations [68], we find that icosahedra play a major role in the emergence of rigidity: icosahedra-rich configurations display much larger yield stresses than icosahedra-poor ones. However, we nuance this statement, as we are able to split the overall family of icosahedral motifs according to the preparation protocol: well-thermalised configurations from trajectory sampling have icosahedral regions that are more robust to shear and with lower energies than the icosahedral domains obtained via energy minimisation. This highlights that the requirement of sampling long-lived structural motifs (implicit in trajectory sampling) allows us to explore metabasins that are not just richer in structure, but more stable as well.
CPR acknowledges the Royal Society for funding. FT and CPR acknowledge the European Research Council (ERC consolidator grant NANOPRS, project number 617266). This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol. FT contributed to the generation and analysis of the data and the writing of the manuscript. FT, TS and CPR contributed to the editing of the manuscript.