Crystallization Instability in Glass-Forming Mixtures

Crucial to gaining control over crystallisation in multicomponent materials or accurately modelling rheological behaviour of magma flows is to understand the mechanisms by which crystal nuclei form. The microscopic nature of such nuclei, however, makes this extremely hard in experiments, while computer simulations have hitherto been hampered by their short timescales and small system sizes due to limited computational power. Here we use highly-efficient GPU simulation techniques to access system sizes around 100 times larger than previous studies. This makes it possible to elucidate the nucleation mechanism in a well-studied binary glassformer. We discover that the supercooled liquid is inherently unstable for system sizes of 10,000 particles and larger. This effect is due to compositional fluctuations leading to regions comprised of large particles only which rapidly nucleate. We argue that this mechanism provides a minimum rate of crystallisation in mixtures in general, and use our results to stabilise a model binary mixture and predict glassforming ability for the CuZr metallic glassformer.


I. INTRODUCTION
Crystallization in supercooled liquids has profound implications in fields as diverse as the development of amorphous materials [1], magma flows in volcanos [2], and aqueous solutions of ions [3]. Materials in question include metallic, inorganic, and chalcogenide glass formers, where mixtures of a number of different constituents have the effect of suppressing or controlling crystallization [4]. Alas, this tendency to crystallize places stringent limits on the size of the pieces of amorphous material that can be formed: large pieces are more likely to undergo crystal nucleation [5,6]. This "Achilles heel" of glass formation thus limits the exploitation of metallic glasses, for example, whose superior mechanical properties otherwise hold great promise [7].
It is clear that any liquid cooled below its freezing point must, for a sufficiently large system, nucleate [8,9].
However, the practical limits of cooling rate versus system size required for vitrification are not known in general. In additional to these practical considerations, crystallization is one solution to the Kauzmann paradox of vanishing configurational entropy upon which a number of theories of the glass transition rest [4,10]: crystallization avoids the need to invoke any particular theoretical description of divergent viscosity in amorphous materials [11][12][13].
It is known empirically that increasing the number of constituent species and introducing a size disparity among these components, together with a negative heat of mixing, tends to suppress nucleation-this has been the guiding principle in the development of bulk metallic glasses [5]. However, despite recent innovative approaches using model systems [14,15] and novel sampling techniques [16,17], there is still a lack of fundamental understanding of the mechanisms by which glass-forming mixtures crystallize.
Here we consider a crystallization mechanism that is always present when a glass former is produced by mixing constituents which by themselves are poor glass formers, as is often the case. We therefore expect this mechanism to be remarkably widespread. In particular, compositional fluctuations in the supercooled liquid lead to regions containing just a single constituent, and eventually such a region will occur that is large enough-and long-lived enough-that it will nucleate a crystal of that one species. Of course, depending on the specific mixture, there may be other, faster, nucleation mechanisms. Nucleation by compositional fluctuations nevertheless provides a lower bound for the nucleation rate in mixtures. We emphasize that compositional fluctuations occur even in the absence of any underlying demixing behavior driven by a thermodynamic transition. In fact, the first mixture we investigate is specifically designed not to demix, using a nonadditive attractive cross interaction between the two species. Thus the compositional fluctuations we consider are distinct from enhanced crystal nucleation rates due, for example, to density fluctuations related to a nearby critical point [18,19]. Clearly, our analysis falls within the concept of the Ostwald rule of stages, suitably generalized to mixtures [20]. In this context, we emphasize that the crystals formed through such compositional fluctuations will in general not be thermodynamically stable.
Having argued for compositional fluctuations as a relevant mechanism for crystal nucleation, we turn our attention to situations in which this mechanism may dominate. We begin with the Kob-Andersen (KA) binary Lennard-Jones mixture. Since its inception in 1994, this model, based on the metallic glass former nickel phosphorous, has been a mainstay of model systems with which to tackle the glass transition [21]. Prized for its simplicity, speed of computation, and its stability against crystallization, the KA model is among the most widely used atomistic glass formers in computer simulations. It is only recently, with the advance of high-performance graphics processing unit (GPU) computing, that the KA model has been crystallized by direct simulation [15], where an estimate was made of the nucleation rate at a single temperature and system size. Here, instead, we carry out large-scale simulations and focus on the mechanism for crystallization.
Our results reveal that nucleation in the KA model is induced by composition fluctuations as discussed above. The KA model is thus representative of systems crystallizing via this mechanism, and we expect our results to have profound consequences for the glass-forming ability of mixtures, such as metallic glasses [7] and oxides [22]. We illustrate the generality of our results by presenting results for the model metallic glass former Cu x Zr 1−x using a range of compositions from x ¼ 0.15 to 0.645.

II. FREEZING IN THE KOB-ANDERSEN MODEL GLASSFORMER
We begin the presentation of our results by studying crystallization in the KA mixture using a global structural analysis for ρ ¼ 1.204 and T ¼ 0.40 in Figs. 1(b) and 1(c). We use the NVT ensemble with a Nose-Hoover thermostat [23]. Figures 1(b) and 1(c) show, respectively, the time evolution of the population of liquid local structures (bicapped square antiprisms) and crystalline structures for system sizes of N ¼ 10 000 and N ¼ 100 000. Here and henceforth we scale time by the structural relaxation time τ α . For T ¼ 0.40, we have that the structural relaxation time τ α ¼ 2.91 × 10 5 simulation time units. A snapshot of a crystal nucleus, composed predominantly of the majority A species, is shown in Fig. 1(a). We identify particles in liquid locally favored structures (LFS) and fcc, hcp crystalline regions with the topological cluster classification (TCC) algorithm [24] and bcc crystalline regions with bond-orientational order (BOO) parameters [25,26]. Our choice of order parameter is motivated by the ability of the TCC to identify the liquid local structure (and hcp and fcc), and we have in any case confirmed that our results for identification of the crystal structures are We observe a fcc-dominated crystallite of A particles. Light green and yellow particles are fcc for A and B particles, respectively, light blue and orange are hcp, purple and dark pink are bicapped antiprism liquid locally favored structure (LFS), and dark blue and light pink particles are liquid. (b) Population of local structures as a function of time reveals rapid crystallization of fcc and hcp for N ¼ 10 000. bcc is found in very small quantities, and the bicapped square antiprism liquid LFS (which is incompatible with fcc and hcp) is also shown. The inset shows that at short times (≲5τ α ), we cannot infer any crystal growth within the fluctuations. (c) Even more rapid crystallization occurs when N ¼ 100 000; here irreversible growth in the fcc population is found within one relaxation time. Inset: Liquid LFS and hcp populations are compared with fcc population, as shown in the main figure. very similar between the two methods (see the Appendix for simulation details and order parameters).
We see from Fig. 1(b) that the liquid begins to freeze on a timescale of a few structural relaxation times τ α . Thus, for these parameters of T ¼ 0.40 and N ¼ 10 000, it is hard to regard the KA mixture as anything but a remarkably poor glass former. We further see that the population of the LFS in the liquid, the bicapped square antiprism [shown in Fig. 1(b)], reduces upon crystal growth in much the same way as in one-component hard spheres where the liquid LFS competes with the crystal symmetry [27].
Here, of course, we have a binary system, but the predominant crystal structures we find are fcc and hcp of the large A species only, and very little mixed AB bcc. The lack of bcc is consistent with predictions that the crystal nucleation barrier is much higher relative to fcc [28] and with the equilibrium KA phase diagram [29]. For the KA model, we therefore neglect the bcc structure and focus on the hcp and fcc crystals in the following.
In Fig. 1(b), we see that there seems to be very little incubation time. However, close inspection [ Fig. 1(b), inset] reveals that for timescales of a few τ α , the fluctuations in crystal population are larger than the increase, so the liquid may in fact be regarded as metastable on short timescales. In Fig. 1(c), we show that upon a further increase of system size, to N ¼ 100 000, this short time metastability vanishes, and the crystal nuclei grow immediately.
We now consider the formation of critical crystal nuclei and estimate their size. In Figs. 2(a) and 2(b), we show the number of particles N xtal in the largest connected region of crystal particles (hcp or fcc) for different system sizes. Here we select a run with a relatively long incubation period [ Fig. 2(a)]. We see that the crystal regions are smaller than 100 particles for around 40τ α before growing. These data enable us to infer a critical nucleus size of approximately 50-100 particles for T ¼ 0.40. Figure 2(b) shows the run at N ¼ 100 000, where crystal growth is immediate and thus it is difficult to infer a critical nucleus size in this case.
Next, we consider the statistics of nucleation in the KA glass former. From the ten runs we performed for N ¼ 10 000 and T ¼ 0.40 (all of which crystallized), we determine the mean nucleation time from τ X ¼ P n i¼1 t XðiÞ =n, where n is the total number of simulations, to be τ X ¼ 38.4 AE 26.8τ α (the error is the standard deviation). Here t XðiÞ is the time when the size of the largest crystal region reaches, and does not subsequently drop below, 100 particles. At higher temperatures, the driving force for crystallization is of course reduced, but the dynamics is much faster. We find that the system does crystallize at higher temperatures (we probed up to T ¼ 0.45) but that not all the runs do so. In this case, we determine the mean nucleation time for each state point following the method of Ref. [30]. In particular, we presume that nucleation is exponentially distributed in time, such that the probability of a nucleation event happening at time t is pðtÞ ¼ 1=τ X expð−t=τ X Þ. The probability that a given run of length t run crystallizes is then The fraction of runs that crystallized then gives us τ X . Errors are estimated by considering the case that one more, or one fewer, simulation runs underwent crystallization. While more sophisticated analyses have been developed, which enable accurate determination of the critical nucleus size [31], even with the considerable computational resources we have used, it has only been possible to carry out ten runs per state point. This limits the extent to which we can implement such methods.
We see from Fig. 3(a) that, when scaled by the relaxation time, the time to nucleate drops rapidly with temperature at N ¼ 10 000, and that well before the dynamical divergence temperature predicted from a Vogel-Fulcher-Tamman fit to the temperature dependence of the relaxation time (T 0 ≈ 0.30, see Supplemental Material [32]), the nucleation time τ X is expected to fall below τ α at T ≈ 0.38. Moreover in the range T ≲ 0.43, we find an exponential scaling with temperature, τ X =τ α ∼ e AT with A ≈ 97. Of course, this observation rests on only the four data points which we fit, but given the significant magnitude of the fall in τ X =τ α with temperature, we are confident that, were this trend to continue, the observation that for some T > T 0 , τ X < τ α would hold. When we simply plot the nucleation time in simulation time units [ Fig. 3(a), inset], we make two observations. Firstly, the absolute nucleation time does not change hugely (around 1 order of magnitude) throughout the temperature range in question, while the relaxation time changes by 3 orders of magnitude. Secondly, there is an upturn at the lowest temperature that we consider, T ¼ 0.395. The reason for the minimum in τ X ðTÞ presented in the inset of Fig. 3(a) is then competition between the decrease in the average nucleation barrier (for a given system size) and the increase in relaxation time upon cooling, though we emphasize that this is only one data point and more statistics would be helpful to confirm this observation. In any case, this is dwarfed by the increase in relaxation time, so the scaled quantity τ X =τ α continues to drop. Turning to the system size dependence of nucleation in Fig. 3(b), we find, as expected, a system size scaling consistent with τ X ∼ 1=N. Note that in Fig. 3(b) we consider a single temperature, T ¼ 0.40, so that τ α does not enter into the scaling.

III. COMPOSITION FLUCTUATIONS
Next, we proceed to investigate the role of compositional fluctuations in crystallization. To quantify these, we use the order parameter illustrated in Fig. 4(b). We seek to find the largest region of liquid A particles which is devoid of any B particles. We presume that such a large compositional fluctuation would be most likely to drive crystallization.
Therefore, we use the following procedure for a given snapshot.
(1) We find the A particle which is furthest away from the nearest B particle.
(2) We define a sphere, centered on the A particle, whose radius is its distance to the nearest B particle. (3) The number of particles in the sphere n s is taken as the current largest compositional fluctuation. (4) We iterate to smaller AB separations and hence smaller spheres, avoiding particles already contained in a previous sphere, and updating n s if a larger region is encountered. Since there will be small fluctuations of crystal particles in the liquid, and we are looking for the largest region of liquid A particles, we seek to avoid the effects of such A particles in crystalline environments. Therefore, we accept a maximum of 10% of the particles in the sphere to be in a crystalline environment. The time evolution of the largest compositional fluctuation in each snapshot n s is shown in Fig. 4(c). We only sample where the system has yet to crystallize, under our criterion of a nucleus size of less than 100 particles. Simply because the system has not yet crystallized does not mean that its properties are stationary, as shown in Figs. 1 and 2. However, it is still instructive to apply the same metric for the larger systems as for the smaller systems (whose properties are stationary for timescales beyond the structural relaxation time), and this we do, with the caveat that the distributions are sampled from a nonstationary system.
In Fig. 4(a), we see that for the KA system at T ¼ 0.40 and N ¼ 5000, 10 000, and 100 000, the distribution of largest composition fluctuations n s of liquid A particles has a significant dependence on the system size N. Two effects are apparent. Firstly, the typical size of compositional fluctuations increases with N. Secondly, the distribution has a "fat tail" indicating more fluctuations of larger n s than a symmetric distribution such as a Gaussian would predict. We note 50-100 particles was a rough estimate of the critical nucleus size and that fluctuations comparable to this are seen in the tails of the distributions. FIG. 3. Nucleation times τ X with respect to temperature and system size in the KA mixture. (a) Nucleation time scaled by the relaxation time τ α as a function of temperature T at N ¼ 10 000. Line is a fit to τ X =τ α ∼ e AT , with A ≈ 97. Inset: Nonscaled nucleation times τ X with a curve to guide the eye (dashed line). (b) Nucleation times as a function of system size at T ¼ 0.40. Line represents expected 1=N scaling for τ X in the case of a constant nucleation rate. We determine τ X as described in the text.

IV. STATISTICS OF COMPOSITIONAL FLUCTUATIONS
What can we say about the origin of the distribution of the largest compositional fluctuations in Fig. 4(a)? Let us suppose that the distribution of all fluctuations of A particles is exponential, Pðn A Þ ∝ expð−n A λÞ, where n A is the number of A particles around a given A particle that are closer than the nearest B particle calculated for every A particle. Here, λ is the decay constant. The extreme values of such a distribution, i.e., those fluctuations large enough to initiate nucleation, should then follow a Gumbel distribution given by P½zðn s Þ ∝ e −ðzþe −z Þ in which z ≡ ðn s − μÞ=β, where μ is the mode of the probability distribution, i.e., the highest probability point, and β is the scale of the function. Furthermore, the median of the extremes of an exponentially distributed process follows m ¼ 1=λfln n − ln½lnð2Þg in which n is the number of samples [33].
In Fig. 4(e), for several system sizes we plot the distribution of all compositional A-particle fluctuations Pðn A Þ. The dashed line indicates that, for large n A , Pðn A Þ exhibits an exponential decay virtually independent of the system size, as expected. This motivates us to fit the Gumbel distribution to Pðn s Þ in Fig. 4(a) (full lines). For N ≳ 10 000, the agreement is remarkable. Moreover, the median of the fitted Gumbel distributions hn s i exhibits a logarithmic dependence on the system size N as shown in Fig. 4(d). We conclude that because the scaling and distribution follows the Gumbel distribution, the largest compositional fluctuation is consistent with exponentially distributed fluctuations.
In fact, from Fig. 4(e) we find that the decay constant of the exponential is λ ≈ 0.22. This value of λ corresponds to a completely random distribution of A and B particles, indicating that the large regions of one species are mainly entropic and thus present irrespective of the particular system; the probability to find n A particles in a cluster of n particles, assuming indistinguishable A and B particles, Merely demonstrating the existence and size of these fluctuations is, of course, not sufficient. We need to show also that they are sufficiently long-lived to initiate the crystallization as well. In order to address this question, we now consider dynamics. At T ¼ 0.40, as noted above, the structural relaxation time τ α ¼ 2.91 × 10 5 . This is wildly in excess of the nucleation time in the one-component system at these temperatures, which is τ X ≈ 41 for a system size of N ¼ 13500 [34]. Thus, since the lifetime of the compositional fluctuations must be on the order of τ α at least, and we do not see any signs of phase separation, i.e., other mechanisms of crystallization [see Fig. 4(c) and composition-composition correlation functions in the SM [32] ], we conclude that the compositional fluctuations we identify lead to crystallization.
Before we explore the compositional fluctuations for other systems, we provide some considerations as to the crystallization mechanism. One alternative possibility is enhancement of nucleation related to density fluctuations. Now the liquid-gas binodal has been measured as lying at a temperature not much less than T ∼ 0.40 to which we simulate [35]. It is conceivable that some density fluctuations related to the proximity of liquid-gas phase separation might act to enhance nucleation, as is known for proteinlike systems [18,19]. However, the system is not in or near the two-phase region: the density of ρ ¼ 1.204 we consider is much higher than the critical isochore (around 0.3). In the  SM [32] we investigate but see little evidence for density fluctuations [35]. In any case, any such nucleation enhancement would still need to invoke a mechanism for A-B demixing, which is absent. Indeed, to observe demixing in similar binary systems, one needs to weaken the interaction between the species so that it is again nonadditive but weaker than the additive case, i.e., a positive enthalpy of mixing [36]. In fact, we see very little evidence for A-B demixing (see SM [32]), (the compositional fluctuations we have discussed notwithstanding). In short, we provide evidence that the compositional fluctuations we identify here are unrelated to the density fluctuations known to enhance nucleation in (effective onecomponent) proteinlike systems [18,19]. We also consider the consequences of our choice of an instantaneous quench protocol (see the Appendix). In Fig. S3 of the SM [32], we see that the median of the largest region of liquid A particles hn s ðTÞi shows very little dependence on temperature. Thus, as the system samples from a nearly temperature-independent distribution and due to the long mean nucleation time for T > 0.40 (more than 100τ α ), we argue that our quenching protocol does not affect our conclusions to any significant extent. The independence with respect to temperature is intriguing: we interpret this in the context that the structure of the liquid is dominated by the hard core [37], in which case a weak temperature dependence is expected.

V. DEPENDENCE OF FLUCTUATIONS ON SYSTEM COMPOSITION
We now consider other compositions of the KA mixture. The temperature independence of the compositional fluctuations suggests that the scaling leading to large compositional fluctuations may be identified at high temperature where timescales are amenable to computer simulation, without recourse to simulations of the deeply supercooled liquid. This suggests that it may be possible to use our approach to predict the glass-forming ability of mixtures in the liquid state.
Usually, as above, the 4∶1 KA mixture is simulated, but upon changing the composition to be more equimolar, we expect smaller regions of pure A particles. We focus on the 2∶1 KA mixture at zero pressure and at the higher temperature of T ¼ 0.80, where the relaxation times for 4∶1 and 2∶1 KA are comparable [38]. We also considered the 3∶1 composition, which turns out to lie close to the 2∶1 system. In Fig. 4(d), we see that at zero pressure and T ¼ 0.80 the 4∶1 mixture has a value of hn s i very similar to that at which we see crystallization (T ¼ 0.40). We infer that the change in pressure also has little effect on the compositional fluctuations, which is reasonable as they are largely random, according to the exponential distribution [ Fig. 4(e)].
As expected, the 2∶1 KA mixture in Fig. 4(d) has very much smaller values of hn s i, as its composition is closer to equimolar. To predict where crystallization might occur, we fit each composition to a logarithmic increase as indicated by the dashed lines in Fig. 4(d). From this we find that the 2∶1 KA system reaches the value of hn s ðNÞi ¼ 31 (corresponding to the 4∶1 system with N ≈ 10 000) at a system size of N ¼ 1.2 × 10 9 . Thus we expect that, for the mechanism of crystallization we consider here, the 2∶1 composition should be a very much better glass former than the usual 4∶1 system. We confirm this by very lengthy simulations of the 2∶1 (and 3∶1) KA systems at comparable supercoolings (i.e., T ¼ 0.40 for KA 4∶1) and N ¼ 10 000, 100 000, and 1 000 000, where no crystallization was observed. We simulated around 9 × 10 9 time steps for N ¼ 10 000 and 100 000 and 2 × 10 9 time steps of N ¼ 1 000 000.
Note that we are only considering the crystallization mechanism based on compositional fluctuations. While we expect the mechanism to be present in all mixtures, crystallization may be dominated by other, faster, mechanisms. For example, the 1∶1 KA mixture forms a mixed bcc crystal quite rapidly [39].

VI. CRYSTALLIZATION IN COPPER ZIRCONIUM
To address whether the mechanism described above pertains to other systems, we consider the metallic glass former copper zirconium. Here we use embedded atom model simulations (see the Appendix for more details). In Fig. 5(a), we show that, like the KA mixture, the extreme values of the composition fluctuations in CuZr also follow a Gumbel distribution. To determine the magnitude of the composition fluctuations in the liquid, we use the higher temperatures of T ¼ 1500 K or T ¼ 1270 K, respectively, so that we can run the simulations quickly. Here we use the NPT ensemble with a Nose-Hoover thermostat [23]; see the Appendix for further details.
The system size dependence of hn s ðNÞi, where we consider fluctuations of the majority species, is shown in Fig. 4(d). Again we see the logarithmic scaling; moreover, compositions such as Cu 64.5 Zr 35.5 exhibit weaker fluctuations compared to the KA model. Following our analysis of the 2∶1 KA mixture, here we estimate that the system may be susceptible to crystallization at a system size of N ¼ 3.5 × 10 16 (when n s ≈ 31). However, upon changing to the more asymmetric compositions Cu 25 Zr 75 and Cu 15 Zr 85 , we see a marked increase in the fluctuations.
Given these larger fluctuations, the logarithmic scaling in Fig. 4(d) would indicate that the metallic glass former should crystallize for those more asymmetric compositions on simulation timescales already around N ¼ 10 000, assuming similar behavior to KA. To investigate crystallization, we run simulations at lower temperatures for two compositions, Cu 25 Zr 75 and Cu 15 Zr 85 . Here the system was first equilibrated at T ¼ 2000 or 1500 K depending on the composition and then rapidly cooled to the temperature of interest. For Cu 25 Zr 75 and Cu 15 Zr 85 , and the temperatures at which we see crystallization, T ¼ 900 and 1100 K, the cooling rates are ΔT=Δt ¼ 3.0 × 10 5 and 2.0 × 10 5 K=ps, respectively.
We find that CuZr indeed crystallizes with representative runs freezing after 930τ α and 95τ α for Cu 25 Zr 75 and Cu 15 Zr 85 , respectively. Here we consider the bcc crystal, as the fcc and hcp are found only in trace quantities. In the snapshot in Fig. 5(b), we find that the nucleus for Cu 25 Zr 75 is dominated by the majority species Zr, in a manner similar to that in Fig. 1(a), although the growth is more rapid in the case of this binary metallic glass fomer (see SM [32]). Note that this higher rate of growth contrasts with slow growth previously observed in CuZr with respect to other metallic glass formers, for a binary crystal [40].
We thus infer that the mechanism for nucleation, at least for these compositions, is the same as that for crystallization in KA. Again, like KA, other mechanisms are also possible in which the crystal may be mixed [41,42]. However, we argue that we have presented a general crystallization mechanism in mixtures, which occurs in the absence of faster, specific, crystallization pathways.

VII. OUTLOOK
Before concluding, we consider the consequences for the long-term stability of supercooled mixtures. We have shown that crystal nuclei are expected in mixtures in general. But by how much should they grow? By considering the KA mixture, the growth of fcc nuclei of A particles will deplete the remaining liquid of A particles. This depletion will tend to slow and may even arrest the growth of the one-component A crystals. In the case of the KA system, we note that if the liquid approaches a 1∶1 composition, then crystallization, not of the one-component fcc, but of the 1∶1 composition bcc crystal, may be expected. We noted in the Introduction that the Ostwald rule of stages, generalized to mixtures, would provide pathways by which the nuclei may grow [20].
Given the small dimensions of the nuclei we find, and despite the developments we present here, our simulations are still small compared to experimental system sizes, and thus it seems reasonable to suppose that the final material may be composed of nanocrystals. Nanocrystals are known to have important consequences for the mechanical properties of glass-forming materials [43]. While this behavior has been seen in experiments [44], our work suggests that such nanocrystals may be rather prevalent in metallic glasses. Because identifying tiny crystalline regions is hard with x-ray scattering, requiring techniques such as 3D atom probe tomography [44], nanobeam electron diffraction [45] or fluctuation TEM [46], it is possible that such nanocrystals may go undetected. The detection of such nanocrystals is an exciting avenue for future research.

VIII. DISCUSSION AND CONCLUSIONS
We have demonstrated a general mechanism of crystallization in multicomponent systems. Our large-scale simulations of the widely used Kob-Andersen model supercooled liquid reveal that it has a fatal flaw as a glass former which is general to mixtures. Local compositional fluctuations lead to regions populated only by one species. These regions can be larger than the critical crystal nucleus size of the one-component system under similar conditions. Nucleation in these regions is fast on the timescale of this deeply supercooled liquid, apparently requiring little rearrangement of the particles, as is known to be the case for hard spheres at deep supercooling [47,48]. Our findings are important, as the results we reveal here pose a fundamental challenge for the development of glass-forming materials: mixtures whose components crystallize easily are themselves inherently unstable to crystallization and thus ultimately compromised as glass formers. Our findings rationalize the empirical rule of thumb that increasing the number of components tends to increase glass-forming  [25]. We see that the nucleus is dominated by Zr.
ablility, as the chances that a critical nucleus of one particular species is formed are reduced in that case. We find a scaling with system size which, once parametrized, may be used to predict the largest system which is stable against crystallization, and therefore the largest pieces of amorphous material which can be prepared from a given mixture. That the compositional fluctuations are rather random and insensitive to temperature suggests that simulations in the liquid at higher temperature where the dynamics are much faster may be used to predict the system size at which crystallization may be expected. We have demonstrated this principle using the 2∶1 (and 3∶1) KA mixtures and have predicted that both can reach system sizes, for comparable simulation times and supercoolings, very much larger than the usual 4∶1 mixture before crystallization occurs. These compositions may thus be used when a better glass former is needed in simulations than the standard 4∶1 model.
The binary model we use demonstrates the use of a mixture to suppress crystallization, as is typically employed in metallic and inorganic glass formers and is encountered in vitreous magmas. Although prevalent and accessible to computer simulation for the model systems we consider, we expect the same mechanism will operate for more general binary mixtures, and indeed for multicomponent systems frequently employed in the quest for ever-better glassforming alloys [7]. We demonstrate this by considering the well-studied CuZr metallic glass former, which exhibits the same behavior. Experimental evidence in support of the mechanism we find has been seen in some metallic glasses [44] and we suggest that the presence of such nanocrystals as we identify here would be worth investigating further in metallic glasses.
Crystallization via compositional fluctuations thus forms a lower bound to nucleation: other mechanisms involving more complex crystal structures may prove faster, as indeed seems to be the case for some models [41,49] and for certain compositions of the Kob-Andersen [39] and CuZr [41] models considered here. Nevertheless, we have shown that liquids which rely on mixing for their stability against crystallization are fundamentally compromised and provide a principle by which their glass-forming ability may be optimized.
In addition to the number of components, crystallization may be suppressed in alloys by the use of systems with a negative heat of mixing. Here, the Kob-Andersen mixture is engineered in that way, precisely to inhibit crystallization. However, our analysis in Sec. IV suggests that such negative heat of mixing has little impact: the Gumbel distribution assumes that the two species are randomly distributed in space, so given its success in describing the statistics, we infer that our analysis is robust to the case where there is a negative heat of mixing. Noting that small size disparities will permit rapid crystallization, and that for certain size ratios binary crystals form an additional route to crystallization as noted in Sec. VII [14,42], we expect that increasing the size ratio will inhibit crystallization. We leave the prospect of a detailed analysis of the role of size disparity for the future.
In addition to the metallic glasses we consider here, an intriguing case is aqueous ionic solutions. Here, crystallization of water occurs through segregation to ion-rich and ion-poor regions, the latter being where the ice nucleates, which appears similar to what we observe here, for the A particles [3]. However, the various anomalies in the thermodynamic behavior of water, not least increasing fluctuations, which may be related to an (avoided) liquid-liquid transition [50,51], mean that further study of that system would be needed to ascertain whether the mechanism we have identified here dominates water crystallization in some aqueous solutions.

Relaxation time determination
For the KA model, we determine the relaxation time of the liquid τ α from the self-part of the intermediate scattering function F s ðq; tÞ ≡ hexp½iqΔri of the A particles using the criterion F sA ðq; τ α Þ ¼ 0.2; the length of the wave vector is q ¼ 7.25. A system size of N ¼ 1000 is used for these simulations to suppress nucleation but has a minor effect on τ α . In the case where we cannot measure τ α directly in simulations due to extremely long simulation timescales, we extrapolate using a Vogel-Fulcher-Tamman fit (see SM for more details [32]). For Cu x Zr 1−x we obtained the intermediate scattering function from the Zr atoms and used a wave vector with q ≈ 26 nm −1 .

Identifying local structure
To detect the fcc and hcp crystals, and the bicapped square antiprism liquid locally favored structure, we use the topological cluster classification (TCC), employed previously to identify local structures in the KA mixture [24]. That is to say, we carry out a standard Voronoi decomposition and seek structures topologically identical to geometric motifs of particular interest.
For the bcc crystal, we employ a bond-orientational order (BOO) parameter analysis [25]. For each particle i we define complex order parameters q i lm ≡ 1=n b P n b j¼1 Y lm ðθ ij ; ϕ ij Þ, where Y lm is the spherical harmonic function with degree l and order m, θ and ϕ are the spherical coordinates for the vector r ij ≡ r j − r i , and n b is the number of neighbors defined from the 12 nearest neighbors. We use the complex order parameters to differentiate between solid and liquid particles using the criteria that for at least 7 nearest-neighbor bonds the scalar product q i 6 · q j 6 =jq i 6 jjq j 6 j should be greater than 0.70 to be classified as a solid particle. q i 6 is a (2l þ 1)-dimensional complex vector. The identity of each solid particle is then determined [55] using the third-order invariant order parameters, where the term in the parentheses is the Wigner 3 − j symbol and Q i lm ≡ 1=ðn b þ 1Þ P n b þ1 i¼1 q i lm is the average bond-orientational order parameter [26]. bcc particles are identified as all solid particles having W i 6 > 0. We checked that the TCC and BOO methods for the detection of fcc and hcp gave similar results.