Glass Forming Ability in Systems with Competing Orderings

Some liquids, if cooled rapidly enough to avoid crystallization, can be frozen into a nonergodic glassy state. The tendency for a material to form a glass when quenched is called “glass-forming ability,” and it is of key significance both fundamentally and for materials science applications. Here, we consider liquids with competing orderings, where an increase in the glass-forming ability is signaled by a depression of the melting temperature towards its minimum at triple or eutectic points. With simulations of two model systems where glass-forming ability can be tuned by an external parameter, we are able to interpolate between crystal-forming and glass-forming behavior. We find that the enhancement of the glass-forming ability is caused by an increase in the structural difference between liquid and crystal: stronger competition in orderings towards the melting point minimum makes a liquid structure more disordered (more complex). This increase in the liquid-crystal structure difference can be described by a single adimensional parameter, i.e., the interface energy cost scaled by the thermal energy, which we call the “thermodynamic interface penalty.” Our finding may provide a general physical principle for not only controlling the glass-forming ability but also the emergence of glassy behavior of various systems with competing orderings, including orderings of structural, magnetic, electronic, charge, and dipolar origin.


I. INTRODUCTION
Glasses are formed from liquids when the transition to the equilibrium crystalline phase is avoided and the relaxation time of the disordered state drastically increases upon cooling [1]. Revealing the physical origin of the glassforming ability is of fundamental significance for the understanding of important materials such as organic materials [2], silicon, water [3], oxides, metals [1,4], semiconductors, and phase-change materials [5], which are used in both crystalline and amorphous states. Previous studies showed that glass formation is often linked to a suppression of the local order that promotes crystal formation [6][7][8][9][10][11][12][13][14][15][16][17][18]. Examples include hard-sphere-like systems, where polydispersity increases the liquid-crystal interface tension [19] by reducing the crystal-like structural order and increasing icosahedral order in the supercooled liquid state [7,20], enhancing the glass-forming ability and decreasing the fragility of the liquid [7]; polymeric glass formers, where the crystallization is controlled by stereoregularity [7,21]; water-salt mixtures, where the salt suppresses the tetrahedral ordering [22,23]; and bulk metallic glasses, such as rare-earth-based multicomponent alloys [4,24], which have competing components and interactions.
From a thermodynamic point of view, destabilization of the crystal state is accompanied by a depression of the melting temperature [7,9,25], as in the eutectic point of metal alloys (which is due to the entropy of mixing contribution). Perhaps the most interesting property of eutectic points is that they usually denote regions of the phase diagram where the system can be vitrified with ease, providing a simple criterion to aid in the search for new glassy states. The enhancement of the glass-forming ability near the eutectic point is demonstrated, for example, by metallic glasses, like the Cu-Zr binary alloy [26] and the Zr-Ti-Cu ternary alloy [27]. Such enhancement of the glass-forming ability can also be observed near triple points, in which two different crystalline phases coexist with the liquid [25], an idea that led Bhat et al. to successfully vitrify germanium as a metallic amorphous state, the first example of glass obtained from a monoatomic liquid [28]. The enhanced glass-forming ability near eutectic points is consistent with Turnbull's criterion [29], which classifies materials in terms of the reduced glass transition temperature T rg ¼ T g =T m , where T g is the glass transition temperature and T m the melting temperature. A high value of T rg (around 2=3) is commonly found in glass formers. At the eutectic point, T m goes through a minimum, and thus, T rg is expected to have a maximum signaling a region of high glass-forming ability. At the same time, glass-forming ability is often displaced from the eutectic composition in binary mixtures [30], such as Ni-Nb [31] and Cu-Zr [32] alloys, meaning that a melting point depression alone cannot be used as a general criterion for glass-forming ability. Such behavior may be observed when there is additional competition from local ordering (such as icosahedral ordering) taking place besides the competition between the two crystal orderings at the eutectic point [7,33].
Turnbull's criterion is a rule of thumb that encompasses several distinct thermodynamic and dynamical effects into a single quantity. The aim of this article is to study the glass-forming ability of two simple systems, where we can control glass-forming ability via an external control parameter. By studying these systems, we can exploit the fact that the glass-forming ability has a nonmonotonic behavior as a function of this control parameter and look for those factors that also display analogous nonmonotonic behavior. In order to encompass a broad family of glass formers, we consider two very different model systems, which are simple and yet general, targeting some important family of glass formers: tetrahedral materials and multicomponent mixtures. We reveal that the key physical factor controlling the glass-forming ability is the structural difference between a crystal and its melt, which can be quantified by a single adimensional quantity, the interface energy scaled by the thermal energy, which we call the "thermodynamic interface penalty".

II. RESULTS AND DISCUSSION
A. Systems studied The first system we study is the Stillinger-Weber (SW) model [34], which is a monoatomic model with two-and three-body interactions (see Appendix A 1 for details). The spherically symmetric two-body interaction favors dense local arrangements and the three-body interactions favor open local arrangements. This potential is characterized by the parameter λ controlling the relative importance of the three-body term [U 3 ðrÞ] over the two-body term [U 2 ðrÞ] (see Appendix A 1), and it is a model for materials with locally tetrahedral coordination such as silicon, germanium, carbon, or water, which are characterized by different values of λ. Molinero et al. [10] introduced the important concept of tuning the relative strength of two-body and three-body interactions in the SW potential in order to control the glass-forming ability of the model. Importantly, in Ref. [10], it was shown that the melting temperature T m reaches a minimum for λ ≃ 18.5, where the stable thermodynamic phase changes from body-centred cubic (bcc) (for λ < 18) to diamond cubic (dc) or diamond hexagonal (dh) (for λ > 19). In this system, the glassforming ability is high in a specific region of the phase diagram where two different crystalline structures compete, and the fluid-to-crystal transition is driven by energy.
The second system we study is a two-dimensional binary disk (BD) system, which is a mixture of large (L) and small (S) hard disks with their diameters σ L and σ S chosen such that σ L =σ S ¼ 1.4 (see Appendix A 1 for details). This system has a eutectic phase diagram, with two competing crystals that are hexagonal phases of large and small disks, respectively. The parameter controlling the glass-forming ability is the composition of the mixture c S ¼N S =ðN L þN S Þ [35], along which the system displays a eutectic point. Here, glass-forming ability originates from atomic strains that are due to size differences between the components, and the fluid-to-solid transition is entirely driven by entropy.
In the first system, crystallization is induced by lowering temperature T, whereas in the second system by increasing pressure at a fixed T. Thus, we may regard the former and latter as thermal and athermal systems, respectively. These different characters of the two systems provide a clue to elucidate the physical factor controlling the glass-forming ability, as will be discussed later.

B. Theoretical basis
To understand the physical origin of the glass-forming ability of these systems, we need to isolate all the factors that contribute to the homogeneous nucleation rate (I).
Here, we note that the temperature is the control parameter for the thermal SW system but not for the athermal BD system. A general definition of a glass former is a system for which I as a function of the control parameters is always low enough such that the system can be quenched to low temperatures without the appearance of a crystal nucleus of critical size, that would otherwise trigger irreversible crystallization. The inability to transform into the stable phase leaves the system in a disordered state. According to classical nucleation theory, the crystal nucleation rate I is obtained as where f is a numerical factor (which depends on terms like the Zel'dovich factor [1]), D=D 0 is the adimensional selfdiffusion constant, ΔF Ã is the free-energy barrier that the system has to overcome in order to crystallize, and β ¼ 1=k B T (where k B is the Boltzmann's constant). Here, we note that the crystallization kinetics is controlled by the self-diffusion and not by the viscosity, which are decoupled for T ≤ T m (see, e.g., Ref. [7]). We write the free-energy cost to form a crystal nucleus, of size N n and interface area S n , in adimensional form by scaling it with the thermal energy, where γ is the interface tension between crystal and liquid, and Δμ ¼ μ solid − μ liquid is the chemical potential difference between solid and liquid. It is important to note that the terms competing in the dimensionless free-energy cost scaled by the thermal energy are βΔμ and βγ, and not Δμ and γ themselves (for convenience, we set the unit of length to 1). Here, we refer to these scaled adimensional quantities as the "thermodynamic driving force" (βΔμ) and "thermodynamic interface penalty" (βγ). Thus, the nucleation rate should not be described by the absolute values of Δμ and γ, but by those scaled by the available thermal energy. This also allows us to have a unified description of the glassforming ability for both thermal (SW) and athermal (BD) systems, as will be shown later. In this picture, βΔF Ã is expressed as where d is the dimensionality and c is a numerical factor depending on the shape of the nucleus.
Here, we note that the macroscopic definition of the thermodynamic interface penalty can be directly linked to the structural difference between two phases via the general Ginzburg-Landau theory [7], where the thermodynamic interface penalty is obtained from the gradients of the order parameter fields, where ρ and Q are density and structural order parameters, x is the coordinate perpendicular to the interface, and K ρ and K Q are positive coefficients associated with the spatial gradients of ρ and Q, respectively. This also shows that the order parameter gradient, which is an intrinsic physical quantity characterizing the penalty associated with the structural difference between the two phases, is characterized by βγ and not by γ. We will see the importance of this fact later. Thus, the nucleation rate I is controlled by three dimensionless physical factors: (i) D=D 0 , the dynamics of mass transport from the melt to the nucleus; (ii) βΔμ, the thermodynamic driving force; and (iii) βγσ d−1 , the thermodynamic interface penalty of forming an interface of area σ d−1 (σ is the typical size of the particles, which we set to unity for convenience) between the melt and the crystal. In principle, all three factors can contribute to the decrease of the nucleation rate in the glass-forming region, but in the following, we will show that the glass-forming ability is mainly controlled by the last term, βγσ d−1 , in the region of competing orderings. Notice that, in the SW system, γ refers to an interface tension, while in the BD system to a line tension.

C. Phase diagrams and thermodynamic driving force
We start our discussion by presenting the phase diagrams of both the Stillinger-Weber (SW) and binary hard disks mixture (BD) systems. Details on the computational techniques are given in Appendix A 1. In Fig. 1(a), we   [36]) to obtain the lines of constant driving force, βΔμ ¼ k, for k ¼ 0.2, 0.4, 0.6. The first observation is that the gradient of the thermodynamic driving force increases towards the melting point depression region. This can be seen from the spacing of the βΔμ ¼ k lines getting more dense towards the glass-forming region. The same is true regardless of the stable crystalline phase, which changes with λ: the driving force grows faster with decreasing T when approaching the melting point minimum from either side. We note that the behavior of the thermodynamic driving force is often studied as a function of T=T m or at constant supercooling T=T m ¼ const (the homologous temperature). Thus, we also study such behavior and the results are shown in Figs. 5-8 in Appendix B 1. As shown in Fig. 7(a), we find that, along the line T=T m ¼ 0.8, the thermodynamic driving force has no anomalous behavior in the glass-forming region, as it just monotonically increases as a function of λ. This is true for any supercooling T=T m , as shown in Fig. 5. Full symbols in Fig. 1(a) represent the homogeneous nucleation lines for both the bcc and dc crystals, here approximately defined as the loci of state points where homogeneous crystallization in the corresponding crystal occurs directly in simulations. Both homogeneous nucleation lines have a stronger λ dependence than the lines of constant driving force (dashed lines), meaning that higher driving forces are necessary to overcome the free-energy barrier for nucleation as the melting point depression is approached. We can conclude that homogeneous nucleation is not hindered by a lack of driving force, but either by dynamical slowing down or by an increase of interfacial free-energy cost. Interestingly, a similar result was obtained for ice nucleation [37]. Given the observation that the crystallization rate is dramatically (by many orders of magnitude) suppressed in the glassforming region [10], we can conclude that the thermodynamic driving force is not responsible for the glass-forming ability of the SW system. In Fig. 1(b), we plot the phase diagram for the BD system as a function of the pressure P and the concentration of small disks, c S (see Appendix A 1). The BD system behaves like a typical eutectic mixture, with the hexagonal L and S solids coexisting with a fluid phase: The melting pressure increases as the eutectic composition is approached from either L-rich and S-rich sides. As in the SW phase diagram, also in Fig. 1(b), we plot lines of constant thermodynamic driving force βΔμ S ¼ βΔμ L ¼ k (dashed lines) and also observe that there is no loss of driving force as the glass-forming region is approached.

D. Transport kinetics
We next consider the dynamics of the system, starting from the SW system. Figure 2 with decreasing temperature is partly compensated by the dynamic anomaly. The net effect is that the suppression of the dynamics along a line of constant thermodynamic driving force is less than a factor of 5, far too low to account for the suppression of the nucleation rate I by many orders of magnitude, as also pointed out in Ref. [40], where the increase of 48 orders of magnitude was observed from λ ¼ 21 to 24. Figure 9 in Appendix B 2 reports on the λ dependence of the isodiffusivity lines, showing that the thermodynamic driving force along these lines becomes lower around the melting point depression. We fit our diffusivity data with the Vogel-Fulcher-Tammann (VFT) is the socalled fragility index and T 0 is the Kauzmann temperature. Note that a higher B value means a less fragile (stronger) liquid. The results are displayed in the inset of Fig. 2(a), where both B and T 0 are plotted as a function of λ, showing that the liquid always remains fragile (B ≪ 100), but around the melting point, the depression region becomes less fragile, with both a decrease of T 0 and an increase of B, as predicted [25]. It is worth noting that we experimentally observed similar behavior [23], including the unimportance of the transport kinetics in the enhancement of the glassforming ability and the decrease of the fragility near the eutectic point, in LiCl-water mixtures, where the controlling factor (equivalent to λ) is the water concentration. In Fig. 2(b), we plot the diffusion coefficient computed from event-driven molecular dynamics simulations of the BD system along different isobars. In this case, we do not distinguish between the two species L and S, as they share the same trends as the average diffusivity reported in Fig. 2

(b). Diffusion increases monotonically with c S ,
showing no anomalous slowing down around the glassforming region. The slowing down of the dynamics is, thus, only due to the increase of pressure on approaching the eutectic point. In the figure, we plot as full symbols the diffusivity computed along the liquidus line, showing that D drops around the eutectic composition by one order of magnitude. As in the case of the SW model [ Fig. 2(a)], this decrease does not account for the suppression of the nucleation rate I in the glass-forming region, which, as we will see in the next section, is almost ten orders of magnitude.

E. Thermodynamic interface penalty
We now turn to the last factor affecting glass-forming ability, i.e., the thermodynamic interface penalty βγ. For the SW system, we employ two different methods to compute the interface tension γ. In the first method, we directly compute the free-energy barriers for nucleation as a function of nucleus size n and apply classical nucleation theory to estimate βγ [41][42][43][44]. In order to disentangle the contribution of βΔμ and βγ on the nucleation barriers, we follow lines where βΔμ is constant, so that any change in the free-energy barrier can only be ascribed to a change in βγ. The results are shown in Fig. 3(a), where the free-energy barrier βΔF is plotted for state points at different λ with constant βΔμ ¼ 0.6. The different colors represent state points where the crystal being grown is bcc (red lines) or dc (blue lines). We note that other forms of crystals (like β-tin) are not nucleated in our simulations due to their high free-energy barriers. The figure shows that the barriers drastically increase as one approaches the glass-forming region from either the bcc or the dc side. for the formation of the bcc crystal (red lines) and the dc crystal (blue lines) for state points at constant thermodynamic driving force, βΔμ ¼ 0.6. The system size for these simulations is N ¼ 4000 particles. The inset shows βγσ 2 extracted from the fit of the free-energy barrier with the CNT expression: where R is the radius of the nucleus [R ¼ ð3n=4πρÞ 1=3 for a spherical nucleus], jβΔμj ¼ 0.6. (b) Order parameter probability distribution functions Πðc S Þ for the BD system at coexistence at different pressures, P ¼ 12, 13, 15 for the L crystal and P ¼ 19, 20, 21 for the S crystal. The inset shows the height of the barrier (circles) as a function of concentration c S , which is directly proportional to the line tension. Also plotted are the liquidus line (dashed red line) and the state points (diamond symbols) at which Πðc S Þ is calculated (notice that, in this case, the y-axis represents pressures). Pressure is in adimensional units with σ L ¼ 1 and β ¼ 1.
Fitting the barriers with the classical nucleation theory (CNT) expression (and with the Tolman length [44]), one obtains a height of 23.5k B T for λ ¼ 23.15 (a good crystalforming system) that rapidly grows to 3800k B T at λ ¼ 20.15 (completely suppressed crystallization). The same behavior is observed approaching the melting point depression region from low λ. The inset of Fig. 3(a) shows the value βγσ 2 extracted from a classical nucleation theory fit of the free-energy barriers at different values of λ. This result shows that the increase of the nucleation barrier is due to a steep increase in the interfacial free-energy penalty scaled by the available thermal energy βγ in the melting point minimum region.
In Appendix B 1, we do the same calculation of freeenergy barriers along a line of constant supercooling, T=T m ¼ 0.8 (Fig. 6). Figure 7 reports both the thermodynamic driving force (a) and the thermodynamic interface penalty (b), showing that, while βΔμ has a monotonous behavior across the glass-forming region, the same is not true for βγ, which is sharply peaked in that region. βγ is, thus, the control parameter for glass-forming ability. It is important to note that the same is not true for the bare interface tension γ. This is plotted in Fig. 8(b), where we see that γ behaves differently on the bcc-stable and the dcstable sides: It decreases towards the glass-forming region from the bcc side, while it increases from the dc side.
For consistency, we also check if the same behavior is observed for the interface tension at coexistence in a thin slab geometry. We compute the free-energy difference between a two-phase system with a planar interface and the bulk fluid (or crystal) phase at melting. This calculation is reported in Appendix B 3 (see Fig. 10), where we show the interfacial free-energy penalty along the melting lines (i.e., where the driving force is null, βΔμ ¼ 0). The inset of Fig. 10 indicates that βγσ 2 has a steep maximum in the melting point depression region, in agreement with the behavior of the same quantity calculated using CNT [see the inset of Fig. 3(a)]. In Fig. 11 in Appendix B 3, we also show the λ dependence of the unscaled bare values of the interface tension γ estimated by the two methods.
In Fig. 12(b) in Appendix B 4, we show the density difference Δρ between the melt and the stable crystalline phase. By comparing the λ dependence of Δρ in Fig. 12(b) and that of γ in Fig. 11, we can see that the quantity γ correlates with Δρ (see also Appendixes B 3 and B 4). As can be seen in Eq. (4), the information on the spatial gradient of the order parameters is hidden by the λ dependence of the available thermal energy. The above results show that the density (ρ) difference between the liquid and crystal is asymmetric between the bcc and dc sides [see Fig. 12 According to Eq. (4), this means that the density and structure order parameter differences should contribute asymmetrically and symmetrically to βγ, respectively. The rather symmetric shape of βγ [see the inset of Fig. 3(a)] indicates that βγ is largely dominated by the contribution from the structural order parameter gradient (j∇Qj) across the interface, and is rather weakly affected by the density gradient (∇ρ). For the BD system, we compute the full density of states of the system for different pressures. From the density of states, we thus obtain the order parameter distribution Πðc S Þ, as shown in Fig. 3(b) (see Appendixes A 3 and B 6 for the details). It is known [45] that the ratio log ½Π max =Π min is directly proportional to the free-energy cost of forming an interface between the two phases and can be regarded as an adimensional measure of the interface tension. We plot log ½Π max =Π min in the inset of Fig. 3(b). We note that it was not numerically feasible to obtain the barrier heights closer to the eutectic point than what is plotted in Fig. 3(b). But the trend of increase of the barriers can be expected to continue when approaching the eutectic point even closer. Thus, the BD system experiences a rapid increase of the bare interfacial free-energy cost (∝ γσ) in the glass-forming region. The free-energy barrier to nucleation depends exponentially on the square of the interface tension [note that d ¼ 2 in Eq. (3)], which allows us to estimate around ten orders of magnitude suppression in the range considered in Fig. 3(b).
Here, we emphasize again that, as shown in the Introduction, it is crucial to consider the behavior of Δμ and γ scaled by the thermal energy, i.e., βΔμ and βγ, respectively, to correlate them with glass-forming ability. In particular, the penalty coming from the structural difference should be measured by the thermodynamic interface penalty βγ and not by the bare interface tension γ [see Eq. (4)]. The control parameter in the SW and BD systems is different (temperature for SW and pressure for BD), and the nature of the melting point minimum is also different in the two cases (it is a triple point for SW and a eutectic point for BD), but the mechanism controlling the glass-forming ability is the liquid-crystal structural difference, macroscopically quantified by the adimensional thermodynamic interface penalty βγ, in both systems. To further strengthen this common feature and provide the underlying microscopic mechanism behind this macroscopic rule, we now consider the glass-forming ability from a microscopic perspective.

F. Microscopic origin of glass-forming ability
It is often suggested that the similarity between a crystal and its melt is an important factor controlling crystallization. Several works on a variety of simple potentials [44,[46][47][48][49] have suggested that the relevant symmetry expressing this similarity is bond-orientational order. The fluctuations in bond-orientational order that promote crystallization are referred to as "crystalline precursors." Like crystal nuclei, particles in precursors have high bondorientational order, and they are orientationally correlated with a finite correlation length. However, differently from nuclei, the local density of crystalline precursors is still close to the one of the melt, and they do not possess local translational order. Nucleation is observed to occur from these regions [44,[46][47][48][49]. We show here that the increase in the glass-forming ability towards the melting point minimum is linked with the suppression of crystal precursors. Crystalline precursors share the same local symmetry of the corresponding crystalline structures: For the SW system, we identify them based on the same Q 12 order parameter used to identify crystals (see Appendix A 4); for the BD system, we instead employ local hexatic order ψ 6 as a measure of crystallinity in the system (see Appendix A 4). Both Q 12 and ψ 6 measure bond-orientational order. Figure 4 plots the fraction of precursors (open circles) as a function of λ for state points with βΔμ ¼ 0.6 in the SW model (a), and as a function of c S along the liquidus line of the BD model (b). In both cases, we see that the glassforming region corresponds to a minimum in the number of precursors, or the degree of crystalline bond-orientational order. The top panels in Fig. 4(a) display precursor particles in snapshots at low (left), intermediate (center), and high (right) values of λ, with the glass-forming region clearly devoid of precursors for the SW model. The decrease of precursors towards the triple point is consistent with the above argument that the structural difference between the liquid and crystal increases towards the triple point, which is inferred from the steep increase of βγ and Eq. (4). For the BD system, in Fig. 4(b) we plot the average hexatic order ψ 6 for different isobars (open symbols) and for state points along the liquidus line (full symbols). Differently from that observed in the analysis of the diffusion coefficients [ Fig. 2(b)], the hexatic order shows a pronounced minimum for each isobar around the eutectic composition. Again, structural order of the fluid phase is strongly suppressed in the glass-forming region, in accordance to what found in the SW system. Figure 4(b) also shows snapshots where precursor regions whose disks i with ψ 6 ðiÞ > 0.6 are colored in red; we notice that, indeed, the length scale of the precursor regions becomes smaller in the glassforming region. This analysis highlights that the glassforming region of the liquid has the largest structural difference from both crystals, since it contains two competing local arrangements, while each crystal contains only one. In other words, our analysis confirms that structural similarity between the liquid and the crystal, which can be characterized by βγ [see Eq. (4)], is a key feature of good crystal formers.
Near the melting point minimum, the liquid structure is more disordered due to competing ordering, suggesting a larger difference in the free-energy difference between the liquid and crystal, βΔμ. Glass-forming ability can thus only originate from a concomitant increase of the interfacial term βγ.
Here, it is worth pointing out that there is an important difference in the coupling of structural ordering with density between the SW and BD model. In the above, we showed that the bare interface tension γ correlates with the density difference between the melt and the stablecrystalline phase in the SW model (see Appendixes B 3 and B 4). For liquids of directional bonding, such as the SW model including water, local structural ordering accompanies a decrease in the local density [50]. Thus, the crystal-liquid density difference linked to the degree of local structural ordering in the liquid is an important factor controlling the crystal-liquid interfacial energy. This is the case for the SW model for λ > 18, where low-density tetrahedral order is the local structural order. It should be noted, however, that the contribution of ∇Q to βγ in Eq.  Fig. 3(a)]. In other words, the dominant contribution to βγ comes from the structural order parameter difference between the liquid and crystal rather than the density difference, as discussed above. On the other hand, for liquids with isotropic interactions such as the BD model, there is little coupling between structural ordering and density [46][47][48][49]. In this case, the interface tension is controlled solely by the difference in bond-orientational order.
In Fig. 4(a) (right scale), we also plot both the configurational (filled squares) and vibrational (open diamonds) contributions to the entropy of the melt as a function of λ. As can be seen in Fig. 4(a), the glass-forming region shows an anomalous increase of s conf in the melt, suggesting a close connection between βγ and the configurational entropy difference between the liquid and crystal (with zero configurational entropy). In Appendix B 5 (see Fig. 13), we show that the configurational entropy s conf is indeed the largest contribution to the free-energy barrier to crystallization. We note that this result is plotted along a line of constant βΔμ, showing that the increase of the configurational entropy due to competing orderings (bcc vs dc) overcomes the rapid decrease of T toward the triple point (usually entropy decreases with lowering T, as seen, for example, in the behavior of the vibrational entropy). We stress that since the entropy s contributes to free energy as −Ts, s conf should be linked to βγ rather than γ. We note that this behavior of the configurational entropy is fully consistent with that of crystal precursors: Both are the measures of the structural complexity of a supercooled liquid. In the simple Adam-Gibbs scenario [51], the increase of configurational entropy is in agreement with the decrease of fragility that was observed in the inset of Fig. 2(a), where the glassy dynamics becomes stronger in the glass-forming region.

III. CONCLUSIONS
To conclude, we have investigated the origin of the glassforming ability in two model systems where glassiness is signaled by a depression of the melting point as a function of some parameter (λ for the SW model, and c S for the BD model). By isolating all the factors that control the crystal nucleation rate, we show that neither the difference in bulk free energy nor the transport kinetics can account for the enhanced glass-forming ability, while the drastic suppression of the nucleation rate is in fact due to a rapid increase of thermodynamic interface penalty, βγσ 2 , which is a macroscopic measure of the liquid-crystal structural differences [see Eq. (4)]. From a microscopic point of view, the increased interface penalty mirrors the increased structural differences between the crystal and the melt, i.e., the disappearance of crystal precursors. For the SW model, by further investigating the thermodynamics and structure of the melt and the crystals, we also find an increase in configurational entropy in the melt around the melting point depression region, reflecting the large number of local arrangements that are possible when two (or more) local orderings are in competition, which also makes a liquid less fragile (i.e., stronger) [7].
In the models studied here, we have considered both the competition between an isotropic attraction and an anisotropic repulsion (SW model), and the competition induced by disorder factors (concentration in the BD model). Furthermore, the former is a thermal system, whereas the latter is an athermal one. Despite their completely different characters, the phenomenology of glass-forming ability is essentially the same between the two models. This provides strong evidence that the enhancement of glassforming ability is caused by an increase in the structural difference between the liquid and crystal, which is signaled by the thermodynamic interface penalty βγ (not γ). This point should be checked carefully for various systems, including realistic systems, in the future.
In relation to the above, it is worth considering other types of frustration to crystallization. Since we put our main focus on the glass-forming ability in systems with competing crystalline orderings in this work, we have not considered systems where a different type of noncrystalline ordering can come into play (e.g., icosahedral ordering), which can represent another important source of glassforming ability [7,18,33,52,53]. Reference [18] is particularly important, as it shows the effect of the competition of different locally favored structures on the crystallizability of a large variety of spin liquids with a common orderdisorder transition (one crystal). In their work, they show that locally favored structures can have both an agonist or antagonist behavior towards crystallization, and also, in their case, crystal affinity correlates linearly with the interface tension. We note that, even in such cases, our theoretical consideration (see Sec. II.B) suggests that the key controlling factor of glass-forming ability may be the thermodynamic interface penalty βγ.
The general expression for the thermodynamic interface penalty in the Ginzburg-Landau theory, Eq. (4), establishes the direct connection of βγ (not γ) to the difference of the order parameters between the disordered and ordered phases. The relevant order parameters are determined by the type of ordering, such as crystalline order parameters (more precisely, the density and crystalline bondorientational parameter) for the systems considered in this work. But a similar mechanism should hold also for systems where glassy behavior stems from a competition between competing orderings and a disorder effect, including materials with first-order magnetic, electronic, charge, orbital, and dipolar ordering [54][55][56][57][58][59]. Finally, we note that the principle found here is consistent with available experiments on glass-forming ability [22,23,25,28]. We hope that our finding will provide a general physical principle for controlling the glass-forming ability in various materials, including metallic alloys and phase-change materials, and open a new avenue towards physics-based material design of various types of glassy materials.

Simulation methods a. SW model
The SW potential is given by where λ is the parameter controlling the relative importance of the three-body term [U 3 ðrÞ] over the two-body term [U 2 ðrÞ]. The two-body term comprises a steep repulsion at very short separations and a short-range attraction: The three-body term is a directional repulsion that promotes specific conformations in between triplets of particles, modeling the effect of directional interactions: All microscopic parameters are the same as in the SW parametrization of silicon [34], except for λ, which becomes the control parameter. The system size is N ¼ 1000 particles if it is not mentioned. In the following, we will always consider isobars at P ¼ 0; see Ref. [36] for full phase diagrams of the model. We use reduced units where ϵ and σ are the units of energy and length, respectively. In the original SW parametrization for silicon, the units correspond to ½T silicon ¼ 25, 157 K and ½p silicon ¼ 377, 674 bar. Simulations are conducted in the NPT ensemble, allowing anisotropic responses to the pressure [60]. Absolute entropies are calculated via standard freeenergy techniques (Widom insertion method and thermodynamic integration for the melt, and Einstein crystal for the solid) and checked against grand-canonical simulations [36]. The vibrational contribution to the entropy is computed via normal-mode analysis. Configurations are quenched to the inherent state via a conjugate-gradient algorithm, and the vibrational entropy is computed in the harmonic approximation by calculating the eigenvalues (ω 2 i ) of the Hessian matrix: We have checked that the anharmonic correction is negligible at the temperatures investigated. Integration of coexistence points and constant βΔμ lines in the ðT; λÞ plane are conducted with Hamiltonian integration [61]. Free-energy barriers are computed with the CNT-US method described in detail in the Supplemental Material of Ref. [44]. Diffusion coefficients are computed with NVE molecular dynamics simulations.

b. BD model
The interaction potential of the BD system is given by where the two species have diameters σ L and σ S , respectively. In the following, we fix σ L =σ S ¼ 1.4, which gives rise to a eutectic phase diagram, with two competing crystals that are hexagonal phases of large (L) and small (S) disks, respectively. Monte Carlo simulations are conducted in the NPT ensemble. The system size is N ¼ 1058 disks. Cluster moves are employed for the volume step [62], and event-chain moves with swap moves for the displacement step [63]. The density of states is computed with Wang-Landau simulations, which we have adapted for the semigrand-canonical ensemble [64,65]. Wang-Landau simulations are further improved with parallel tempering simulations between windows with different compositions.

Integration of coexistence points is obtained with
Gibbs-Duhem integration techniques in the semigrandcanonical ensemble [66]. We note here that, in the thermodynamic limit, two-dimensional ordering of hard disks happens via a weakly first-order phase transition to a hexatic phase, and a second-order phase transition to the crystalline phase [63], but these transitions are only observed with a very large number of disks, while, at our system sizes, the transition appears as a first-order transition between a liquid and a solid phase, as shown in Appendix B 6. For our purposes, we can thus safely ignore the difference between the hexatic and solid phase, as they appear as a single peak in the order parameter distribution functions. Diffusion coefficients are computed with event-driven simulations with the DynamO package [67]. See Appendix B 6 for further details on simulations.

Estimation of βγ for the SW model
In order to isolate the interfacial effect from that of the driving force, we compute the free-energy barriers along a path of constant driving force, βΔμ ¼ 0.6. In this way, the difference between the free-energy barriers can only be ascribed to a change in the term βγσ 2 , where again σ is the unit of length. The main difficulty in the calculation of the barriers is the choice of a proper order parameter for the identification of crystalline nuclei. In our case, this issue is made more difficult by the fact that we need an order parameter that is able to detect all four crystal phases that are relevant in the glass-forming region (cubic and hexagonal diamond, bcc, and β-tin). Very recently [44], we developed a novel order parameter, called Q 12 (see Appendix A 4), which is able to detect all the relevant crystalline structures of tetrahedrally coordinated materials. We thus employ the Q 12 order parameter to bias the formation of crystalline nuclei in umbrella sampling (US) simulations and, thus, compute the free-energy barriers for crystallization. In particular, we use the CNT-US scheme, a variation of the US scheme that allows the computation of the nucleation barrier in a very efficient way [44].

Estimation of βγ for the BD model
For the BD system, we compute the full density of states of the system for different pressures. From the density of states, we thus obtain the order parameter distribution Πðc S Þ. A detailed explanation of the method is contained in Appendix B 6. At coexistence, Πðc S Þ is a doubly peaked distribution, one peak representing the liquid phase, and the other the solid (either L or S solid), with a minimum in between the two peaks (see Fig. 14 in Appendix B 6). The probability distributions for several state points along the coexistence line are plotted in Fig. 3(b). The Πðc S Þ function is normalized such that its peak value at coexistence is unity, Πðc S S Þ ¼ Πðc L S Þ ¼ 1, where c S S and c L S are the concentration of the solid and liquid phases, respectively. In this way, it is easily observed that the height of the barrier is increasing as we increase the pressure on either side of the eutectic point. As is well known [45], the ratio log ½Π max =Π min is directly proportional to the free-energy cost of forming an interface between the two phases, βF s ¼ lim L→∞ log ½Π max =Π min =2L, where F s is the line free-energy cost (i.e., the line tension) and L is the size of the system. We thus use the ratio log ½Π max =Π min as an adimensional measure of the interface tension.

Precursors and crystallinity
To identify crystal particles in the SW model, we use local bond-order analysis, following Ref. [44]. A (2l þ 1)dimensional complex vector (q l ) is defined for each particle i as q lm ðiÞ ¼ f1=ðN b ðiÞg P N b ðiÞ j¼1 Y lm ðr ij Þ, where we set l ¼ 12, and m is an integer that runs from m ¼ −l to m ¼ l.
The functions Y lm are the spherical harmonics, andr ij is the normalized vector from particle i to particle j. The sum goes over the first N b ðiÞ ¼ 16 neighbors of particle i. We then introduce a spatial coarse-graining step Q lm ðiÞ¼ f1=½N b ðiÞg P N b ðiÞ k¼0 q lm ðkÞ [68]. The scalar product between Q 12;m of two particles is defined as Q 12 ðiÞ · Q 12 ðjÞ ¼ P m Q 12;m ðiÞQ 12;m ðjÞ. If the scalar product ½Q 12 ðiÞ=jQ 12 ðiÞj · ½Q 12 ðjÞ=jQ 12 ðjÞj between two neighbors exceeds 0.75, then the two particles are deemed connected. A connection between two particles expresses the fact that their local environment, as quantified by Q 12 , is similar. We then identify particle i as crystalline if it is connected with at least 12 neighbors. A precursor is then defined as a particle with five or more connected neighbors.
For the BD model, we use the definition of local hexatic order, where n j are the neighbors of particle j, and θ jk is the angle that the bond between particle j and k makes with a reference axis. In the above definition, neighboring particles are defined as particles sharing an edge in the radical Voronoi diagram obtained from the particles' positions and sizes. In Fig. 5, we plot the thermodynamic driving force βΔμ as a function of T=T m for different values of λ, in the bcc (red lines), β-tin (black lines), and dc (blue lines) stable regions. The figure shows that the rate of change in driving force is a monotonous function of λ, with the rate being slowest in the bcc region and fastest in the dc region.
The nucleation barriers in Fig. 3(a) were computed along lines of constant driving force βΔμ ¼ 0.6. In Fig. 6, we instead compute the barriers along a line of constant supercooling T=T m ¼ 0.8, showing a rapid increase of the activation energy on approaching the glass-forming region, either from the bcc side (a) or the dc side (b).
In Fig. 7, we plot the thermodynamic factors along the line T=T m ¼ 0.8: both the thermodynamic driving force (a), obtained via thermodynamic integration, and the thermodynamic interface penalty (b), obtained by a oneparameter fit of the barriers in Fig. 6. We observe that, while the thermodynamic driving force (βΔμ) has a monotonous behavior across the glass-forming region, the thermodynamic interface penalty (βγ) steeply increases towards the glass-forming region, indicating that βγ is the factor controlling glass-forming ability.
We highlight here the importance of using adimensional parameters. Figure 8 shows the bare Δμ (a) and γ (b) factors along the homologous line. We observe that the behavior of these thermodynamic quantities is very different between the bcc-stable side and the dc-stable side. The chemical potential difference Δμ decreases sharply from the dc side with decreasing λ, but in the glass-forming region, it is still bigger than the low values measured on the bcc side. This is consistent with the estimates in Ref. [10]. The interface free energy γ increases while approaching the glass-forming region from the dc-stable side. This is similar to the pressure effect observed for γ in simulations of water [37]. But on the bcc-stable side, we instead observe a decrease of γ when approaching the glass-forming region. These results clearly indicate that γ alone is not a good control parameter for glass-forming ability, and the thermodynamic interface penalty (βγ) should be used instead, as shown in Fig. 7(b). We stress that only the behavior of the thermodynamic interface penalty (not the bare interface tension) provides a coherent picture for the thermal SW and athermal BD systems.

Isodiffusivity lines in the SW model
We show in Fig. 9 the phase diagram of Fig. 1(a) Fig. 2(a), the dynamics becomes anomalous around the melting point depression, and less thermodynamic driving force is required to produce the same diffusivity compared to the regions away from the melting point depression. Departure from βΔμ ¼ k lines happens closer to the melting point depression as the system becomes slower.

Excess free energy for interface formation in the SW model
In the main text [ Fig. 3(a)], we computed the free-energy barriers along lines of constant driving force βΔμ ¼ 0.6, which displayed a steep increase in the triple-point region. From a CNT fit to the barriers, we then obtained numerical estimates of the adimensional interface tension contribution, or the thermodynamic interface penalty βγσ 2 , plotted in the inset of Fig. 3(a). While unambiguously showing an increase of interface free-energy cost in the triple-point region, the value of βγσ 2 depends on the particular choice of the fitting function. In this section, we compute βγσ 2 without relying on classical nucleation theory. It has been recently demonstrated that it is possible to reversibly construct a path between a bulk fluid phase and a twophase system with a planar interface separating the fluid from the crystal [69,70]. We obtain this by use of umbrella sampling simulations with the order parameter Q 12 defined in Appendix A 4. We use a thin-box configuration, where the width, height, and length are in the ratio 1∶3∶10. In this configuration, the crystal nucleus transforms from a cylindrical shape to a slab perpendicular to the length direction. Shape transitions are continuous and allow a reversible sampling in the fluid-solid coexistence region, from the fluid to the slab geometry, where the interface is flat and there is no free-energy penalty in moving the interface [71]. This is represented in Fig. 10, where the free-energy profile as a function of the number of solid particles n, and for different values of λ, is reported. The fluid basin corresponds to n ¼ 0, while the two-phase system with the crystal in the slab geometry is encountered at high values of n > n Ã , where βΔFðnÞ becomes flat. Here, the flat profile of βΔFðnÞ signals that we are correctly at the melting temperature, where both fluid and crystalline phases are equally likely. The free-energy difference represents the cost of the interface formation where S is the area of the interface (the width and height directions). The numerical results are reported in the inset of Fig. 10 as red squares, showing again a steep increase in the triple-point region.
The results for the bare interface tension γσ 2 are compared in Fig. 11. Red squares represent the interface tension computed from the free-energy cost of forming a thin slab interface at coexistence. The black circles are instead the interface tension computed with a CNT fit to the free-energy cost of forming spherical nuclei at βΔμ ¼ 0.6. The dependence on the tetrahedral parameter λ is similar for both measurements, but a relative discrepancy of 15% exists between them. Since our focus is on the free-energy cost of nucleation, we extract the thermodynamic interface penalty directly from the nucleation barriers with classical nucleation theory. In principle, the solid-liquid interface tension depends on the orientation of the interface, but most simulations point to the fact that this dependency is very weak, and nuclei are mostly spherical [72]. It is nevertheless interesting to consider the interface tension at coexistence. Above, for numerical convenience, we have only considered the case of a thin-slab geometry, but recent extensions of the capillary wave method [73] or the mold integration technique [74] allow for a full calculation of the interface tension of a flat interface.
The biggest source of error originates from the CNT fit, which is applied only to relatively small nuclei, where the capillary approximation is more likely to fail.

Pair correlators and liquid-crystal density difference as a function of λ in the SW model
Here, we consider common parameters that fail to capture the onset of glassy behavior. In particular, we consider the two-body radial distribution function gðrÞ and the density difference between the melt and the crystal.
In Fig. 12(a), we plot the radial distribution function for different liquid state points along a line of constant thermodynamic driving force to crystallization, βΔμ ¼ 0.6. Also plotted are the gðrÞ functions for the corresponding stable crystals at λ ¼ 14.5 (bcc), λ ¼ 18.1 (β-tin), and λ ¼ 23.15 (dc). The figure shows that the melt undergoes considerable restructuring as a function of λ, and that this change closely mirrors the local structure of the corresponding crystals. At low values of λ, the local order is a dense structure with eight nearest neighbors, while, at high values of λ, it is characterized by tetrahedral ordering with four nearest neighbors. As we move to intermediate values of λ (the triple-point region), the melt progressively acquires competing local orderings, as the relative strength of three-body interactions changes [75]. We note that a similar change takes place as a function of pressure P (instead of λ), which is the origin of the V-shaped T-P phase diagram of watertype liquids such as water, Si, and Ge [25]. The maximum 14  rate of change in the structure occurs around λ ≈ 18, where we observe a jump in the position of the first minimum [ Fig. 12(a)]. This region of fast structural rearrangement corresponds to the location of maximum diffusivity observed in Fig. 2. While the radial distribution functions in Fig. 12(a) show significant structural changes in the liquid, it is not evident why these would lead to a large increase in the thermodynamic interface penalty in the glass-forming region. Another simple measure of structural difference between a liquid and a solid is their difference in density. In Fig. 12(b), we plot the density of the liquid (continuous line) and that of all relevant crystalline structures (symbols) for state points at different λ but constant thermodynamic driving force, βΔμ ¼ 0.6. The difference between the melt and solid densities is represented by the dashed line (shifted for clarity). The dc crystal is an open crystalline structure, with a density lower than the density of the melt, while both the β-tin and bcc structures are denser than the melt. The behavior of the dc crystal is different from the other structures: While the dc crystal barely changes its density as a function of λ, both β-tin and bcc have density changes that mirror the one of the melt. So while the density difference increases noticeably in the dc region going towards the triple-point point, it is instead far less pronounced on the bcc/β-tin side, where it is almost constant with λ. The density difference between the crystal and the melt is very well correlated with the λ dependence of the bare interface tension γσ 2 shown in Fig. 11. This confirms that the density difference plays an important role in determining the interface free-energy cost, as is well known; however, the contribution of the density difference to the thermodynamic interface penalty βγσ 2 is much smaller than that of the structural order parameter difference.
Here, it is worth discussing a link between structural ordering and density in the liquid state of λ > 18. For models with directional bonding, such as the SW model, local structural ordering automatically accompanies a decrease in the local density [50]. So high λ liquids with higher structural order have a lower density, which is closer to the density of the crystal. The degree of structural ordering in a liquid state decreases with a decrease in λ for λ > 18 since the degree of frustration increases due to competing ordering. This explains the increase of the density difference between the liquid and the solid with a decrease in λ (>18).

Energy and entropy contribution to the free-energy barrier in the SW model
In Fig. 13, we plot separately the free-energy barrier (βΔF) and its entropic contribution (s) as a function of nucleus size n for λ ¼ 23.15. The entropic contribution is obtained by calculating the average energy of particles in nuclei of size n (the inset of Fig. 13) and subtracting it from the total free-energy barrier obtained from CNT-US simulations. The figure shows that the transition is energy driven and that most of the barrier is of entropic origin, meaning that the energy and entropy of the bulk crystal are both lower than those of the liquid. The same is true for all values of λ considered in this study.
Through normal mode analysis, we can further separate the entropy in its vibrational (s vib ) and configurational (s conf ) contributions. Loosely speaking, s conf measures the number of equilibrium configurations that are separated by some structural rearrangement; importantly, a crystal has only one such configuration and, thus, zero configurational entropy. s vib , on the other hand, measures the extent of the vibrations about each of these configurations. We find that the largest contribution to the barrier is due to s conf (s vib ≪ s conf for all n), which is entirely lost when the fluid orders in the crystalline state. In Fig. 4(a) (right scale), we plot both the configurational (filled squares) and vibrational (open diamonds) contributions to the entropy of the melt as a function of λ. The glass-forming region shows an anomalous increase of s conf in the melt, which, being the largest contributor to the nucleation barrier, disfavors the transition to the crystalline state. The increase in s conf signals an increased number of possible global arrangements of the system that arises when two local environments have an unresolved competition. This can be rationalized by thinking of the melt as a mixture of many small patches of different (although rapidly interconverting) local symmetry and size. To this mixture, one can associate a sort of "mixing entropy," which is absent in a homogeneous liquid. This mixing contribution explains the increase of configurational entropy and, in turn, the increase in the interface tension and, finally, the glassforming ability of the system. Note that, as mentioned in the main text, s conf and βγσ 2 are both adimensional quantities, which can be compared on common ground. FIG. 13. Free-energy barrier (βΔF) and its entropic contribution (s) as a function of nucleus size n for λ ¼ 23.15. The inset shows the size n dependence of the energy per particle eðnÞ ¼ EðnÞ=N.

Hard disks' phase transition
In the thermodynamic limit, one-component hard disks undergo a first-order phase transition from the liquid to the hexatic phase and a second-order phase transition from the hexatic phase to the solid phase. These transitions were observed in systems of N ¼ 1024 2 disks [63], a size that is necessary to resolve the very narrow density window in which the hexatic phase is stable. In our study, the system size (N ¼ 1058) is such that the hexatic and solid phases are not resolved, and the transition appears as a first-order transition between a fluid and a solid/hexatic phase [76]. This is seen, for example, with Wang-Landau simulations in the isothermal-isobaric ensemble, in which different density states are sampled uniformly. Starting from the density of states, the probability distribution of the order parameter (the density ρ) can be computed for different values of the pressure P. The coexistence pressure is such that the area under the two peaks is equal, and, in our case, we obtain βP coex σ 2 ≈ 9.09. The ρ distribution at coexistence is shown in Fig. 14. Only two peaks are resolved, one corresponding to a fluid phase at ρ ¼ 0.884, and the other one corresponding to a hexatic/solid phase at ρ ¼ 0.907. The units of length are such that the diameter of the disks is unity σ ¼ 1. In the case of the BD system, we are limited to small system sizes, as the calculation of the phase diagram is computationally very expensive: both simulations in the isobaric ensemble and density of state calculations are prohibitively expensive for hard interactions at large system sizes (despite our use of cluster moves and parallel tempering techniques). We thus limit our study to sizes of N ¼ 1058 and consider the transition to be effectively first order, as the presence of an intermediate hexatic phase would not impact our arguments.