Connecting glass-forming ability of binary mixtures of soft particles to equilibrium melting temperatures

The glass-forming ability is an important material property for manufacturing glasses and understanding the long-standing glass transition problem. Because of the nonequilibrium nature, it is difficult to develop the theory for it. Here we report that the glass-forming ability of binary mixtures of soft particles is related to the equilibrium melting temperatures. Due to the distinction in particle size or stiffness, the two components in a mixture effectively feel different melting temperatures, leading to a melting temperature gap. By varying the particle size, stiffness, and composition over a wide range of pressures, we establish a comprehensive picture for the glass-forming ability, based on our finding of the direct link between the glass-forming ability and the melting temperature gap. Our study reveals and explains the pressure and interaction dependence of the glass-forming ability of model glass-formers, and suggests strategies to optimize the glass-forming ability via the manipulation of particle interactions.

I n studies of glass transition and jamming transition [1][2][3][4][5] , binary mixtures of particles are widely employed to avoid crystallization. If the two components mix up randomly, the particle size mismatch can frustrate the global structural order 6 . However, under certain circumstance, some binary mixtures may undergo phase separation or demixing during the solid formation, i.e., the same types of particles aggregate. Although undesirable in disordered systems, phase separation has attracted much attention in many fields [7][8][9][10][11][12][13][14] .
The glass-forming ability (GFA), i.e., the capacity of a material to resist crystallization and maintain glassy, is fundamental in studies of glasses. A binary mixture prone to phase separation tends to form crystallites of the same type of particles, and hence has a poor GFA. Because glasses are diverse and out of equilibrium, it is difficult to establish a common understanding of the GFA. Moreover, phase separation is one of the multiple forms of crystallization. Under certain conditions, binary mixtures can also form complex crystalline or quasicrystalline structures 15,16 , which complicates further the understanding of the GFA to fight against them. Among various interpretations, geometric and energetic frustrations are thought to be essential to the determination of the GFA [17][18][19][20][21][22][23] .
For binary mixtures, varying the size ratio γ of the large to small particles is an effective way to cause geometric frustration. However, it has been shown that a large γ often promotes phase separation and thus results in poor GFAs for binary mixtures of hard particles, mainly arising from entropy effects 14,[24][25][26][27][28][29] . To suppress phase separation, a moderate size ratio, 1 < γ < 2, is usually adopted. However, it is still unclear if such ratios can prevent phase separation and maintain glassy states after extremely long equilibration. Even with these size ratios, the GFA is sensitive to the particle composition. It has been suggested that the GFA of binary mixtures is the best near the eutectic point or triple point when particle composition is varied 23 .
In this work, we pay more attention to the energetic frustration on the GFA by studying binary mixtures of soft particles interacting via finite-range repulsions. In the zero temperature (T = 0) and zero pressure (p = 0) limit, such soft particles behave like hard ones 30,31 . Taking a mixture in the hard particle limit as the reference, which has a decent GFA with certain values of particle size ratio and composition, we vary the pressure (density) and track the evolution of the GFA. When pressure increases, the potential energy plays a more important role and offers soft particles extra opportunities of phase separation. With the intervention of potential energy, the particle size ratio γ, the cause of geometric frustration, may excite and affect energetic frustration as well and hence affect the GFA. We aim at looking for some simple energetic criteria of the GFA and proposing some energy strategies to manipulate it, via the study of the effects of pressure and particle interaction.
We study two types of widely employed model glass-formers, binary mixtures of soft particles interacting via harmonic or repulsive Lennard-Jones (RLJ) repulsion. Both types of systems with a diameter ratio γ = 1.4 and a 50:50 particle composition have been used in many studies of glasses. At low pressures, both systems exhibit a pressure-independent GFA, equal to that of the hard particle counterparts. At high pressures, the GFA of RLJ systems still remains almost constant, but remarkable phase separation occurs in harmonic systems. Therefore, γ = 1.4 does cause not only geometric but also energetic frustrations in soft particle systems. By performing analytic calculations, we find that the behaviors of the GFA can be explained by the pressure dependence of the melting temperatures of the two components, confirmed by our simulation results. This thus builds up a bridge between non-equilibrium and equilibrium quantities and suggests that the pressure dependence of the melting temperatures of constituent components can be an energetic precursor to evaluate the GFA of glass-formers. We also show that the GFA of the hard particle limit can be achieved at high pressures by a proper modulation of the particle stiffness. Combining our work and previous studies on the composition effects 23 , we expect to achieve a much more comprehensive understanding of the GFA of binary mixtures with effects of composition, pressure, and particle interaction being all included.

Results
Control variables of binary mixtures. We denote the two types of particles in binary mixtures as A and B particles. The particle composition is quantified by the concentration of B particles, c B = N B /(N A + N B ), with N A and N B being the numbers of A and B particles, respectively. There are two quantities to make the constituent A and B particles different. One is the particle size. When the diameter ratio γ = σ A /σ B > 1, A and B particles differ in size, where σ A and σ B are diameters of A and B particles, respectively. The other is the particle stiffness, characterized by the interaction energy scale ϵ AA and ϵ BB as defined in Eqs. (17) and (18) of the Methods section. A particles are softer than B particles if ϵ AA < ϵ BB . Here we manipulate the particle stiffness by letting ϵ AA = (1 + Δ)ϵ AB and ϵ BB = (1 − Δ)ϵ AB with Δ2[−1, 1]. One of our primary goals is to study and understand the pressure dependence of the GFAs of widely employed model glass-formers with γ = 1.4, c B = 0.5, and Δ = 0 (ϵ AA = ϵ BB = ϵ AB ). This manipulation of particle stiffness can facilitate the analytic calculations with Δ being an independent variable, from which the solutions of Δ = 0 can be straightforwardly obtained.
In a rather different perspective from most of the previous studies, we are mainly concerned about the effects of pressure and particle interaction on the GFA. Here we mainly show results of N = N A + N B = 4096 systems in two dimensions. Because the simulations are rather expensive, we have only repeated part of the simulations for larger systems and three-dimensional systems and find consistent results. In the following, we will first present results for harmonic and RLJ systems with γ = 1.4, c B = 0.5, and Δ = 0. Then we will attack the special case with γ = 1 and c B = 0.5. By varying Δ, we are able to sort out potential energy effects without the interference of geometric frustration caused by the particle size distinction. Based on observations of γ = 1, we will derive a generalized picture for γ > 1, from which the results of systems with γ = 1.4, c B = 0.5, and Δ = 0 can be understood. Finally, we will show that the same picture applies to different values of c B .
Characterization of glass-forming ability. For given values of γ, c B , Δ, and p, we start with a liquid equilibrated at about four times of the melting temperature of B particles as defined later, and then decrease the temperature by a small step δT and relax the system for a duration Δt by performing molecular dynamics simulations under constant temperature and pressure. The same procedure is repeated until a solid-like state at a temperature about one-tenth of the melting temperature is achieved. This leads to a quench rate κ = δT/Δt. We compare the GFAs over a wide range of pressures, which have been rarely studied before. It is tricky how to choose a reasonable quench rate to include the pressure effects. Here, we use a dimensionless quench rateκ as defined and discussed in the Methods section, in purpose of giving supercooled liquids at different pressures comparable opportunities to relax their structures.
In the parameter space considered in this work, we have not observed the formation of complex crystalline structures within our simulation time window, so the crystallization mainly takes the form of phase separation 23 . Therefore, the GFA can be well characterized by the degree of mixing of A and B particles against phase separation in the resultant solid-like states, which is quantified here by the parameters where the sum is over all A or B particles, δ n si ;n i is the Kronecker delta, n i is the number of nearest neighbors of particle i, and n si is the number of nearest neighbors which are the same type as particle i. We calculate χ A and χ B for A and B particles separately. Their values are close to 0 when two components mix up well. With more particles being separated, the values grow up and approaches 1 for well-separated states. Therefore, smaller values of χ A and χ B mean better GFAs.
Although χ A and χ B can characterize the GFA well in our work, we should realize that they may not work well to distinguish glasses from complex crystals with A and B particles being mixed 15,16 . This is the limitation of the parameters when studying the GFA against the formation of complex crystals. In that case, one needs to select appropriate parameters to characterize the structural order to distinguish complex crystals from amorphous solids, which is out of the scope of current work.
Glass-forming ability for γ = 1.4 and c B = 0.5 with Δ = 0. Figure 1a, b compare the GFAs of binary mixtures of harmonic and RLJ particles in two dimensions with γ = 1.4, c B = 0.5, and Δ = 0, so there is no cause of strong energetic frustration. These binary mixtures have been widely employed as good glassformers in previous studies.
In the T → 0 and p → 0 limit when particle overlap is tiny, both harmonic and RLJ systems are equivalent to hard particle systems 30,31 and exhibit similar GFAs to that of the hard particle counterpart with γ = 1.4 and c B = 0.5. As shown in Fig. 1a, b, both χ A (p) and χ B (p) tend to approach a constant at low pressures for both systems.
Because of aging, χ A (p) and χ B (p) evolve with the quench ratẽ κ. Fig. 1 compares three values ofκ. At low pressures, χ B (p) (for small particles) remains small whenκ decreases, while χ A (p) (for large particles) grows up. As illustrated by snapshots in Fig. 1c, A particles form clusters whenκ is small, even at low pressures. It is expected that with even smallerκ the clustering or phase separation will be stronger. Here A and B particles have the same concentration. The separation can be weakened by increasing c B to around 0.6 23 in order to leave smaller rooms for A particles to aggregate, as will be shown later, but it remains elusive whether the separation is inevitable as long as the waiting time is sufficiently long 32,33 .
The decrease ofκ at fixed pressure is analogous to the decrease of the compression rate at fixed temperature. The emergence of phase separation at low pressures thus suggests that with sufficiently slow compression rates binary mixtures of hard particles with γ = 1.4 and c B = 0.5 will undergo phase separation. Therefore, we show direct evidence suggesting that such binary mixtures widely employed as good glass-formers cannot prevent phase separation, so that thermodynamically stable states may be phase-separated. Figure 1 also shows that, although behaving similarly at low pressures, harmonic and RLJ systems have significantly different GFAs at high pressures. The GFA of RLJ systems still remains almost constant in pressure. In contrast, harmonic systems undergo apparent phase separation, with χ B (p) increasing quickly when pressure increases. Apparently, at high pressures, some energetic frustration is excited strongly in harmonic systems, destabilizing the mixing of particles. The particle size distinction triggers such energy effects. Moreover, note that harmonic and RLJ systems have almost identical GFAs up to p ≈ 10 −2 , which is already far beyond the hard particle limit where particle interactions are not negligible. Then it is interesting to know why harmonic and RLJ systems can have similar GFAs beyond the hard particle limit, why RLJ systems can maintain an almost constant GFA to high pressures, and why the GFA of harmonic systems quickly drops at high pressures.
Harmonic and RLJ repulsions are distinct in their respectively bounded (soft-core) and unbounded natures. The soft-core nature of harmonic repulsion induces rich phase behaviors at high densities [34][35][36][37][38] . Unlike RLJ particles whose melting (crystallization) or glass transition temperature always increases with the  Fig. 1 Pressure and quench rate dependence of the glass-forming ability. a, b are for harmonic and RLJ systems, respectively, with γ = 1.4, c B = 0.5, and Δ = 0. The solid and empty symbols are χ A and χ B , respectively. Circles, squares, and diamonds are for dimensionless quench rateκ ¼ 8:16 10 À11 , 8.16 × 10 −10 , and 8.16 × 10 −9 , respectively. The lines are guides for the eye. The vertical dot-dashed line in a shows the crossover pressure p n ≈ 0.14 (also shown in Fig. 2a for comparison), at which the melting temperatures of A and B particles intersect and χ A is approximately equal to χ B , as discussed in the text. c Snapshots of solid-like states for harmonic systems. From top to bottom,κ decreases. From left to right, pressure p increases with the values being shown at the bottom. The yellow and blue disks are A (large) and B (small) particles. To distinguish particles, we have moderately decreased the particle diameters by 10-25%.
Melting temperature gap between components. In this subsection, we will introduce two effective melting temperatures, T m,A (p) and T m,B (p), for A and B particles, respectively. We will see that they are not the actual equilibrium melting temperatures, which should be obtained from complicated calculations of the equilibrium phase diagram and would vary with c B 23,50 . They are instead derived from the equilibrium melting temperature of monocomponent systems and the simple conversion of units. However, they turn out to work well to characterize the GFA.
An interesting feature shown in Fig. 1a is that χ A (p) and χ B (p) of harmonic systems intersect at roughly the same pressure p n at different quench rates. When p < p n , χ A > χ B , so A particles are easier than B particles to form clusters. When p > p n , χ A < χ B . The clustering or nucleation starts below the melting temperature. If χ B > χ A , B particles may experience a longer time than A particles to nucleate, suggesting that the melting temperature of B particles is higher, and vice versa. Therefore, we doubt whether Fig. 1a implies that B particles have a lower (higher) melting temperature than A particles when p < p n (p > p n ). If this is the case, melting temperature would play an important role in the determination of the GFA, but the question is how A and B particles can feel different melting temperatures in the same mixture.
Note that temperature is the energy, and we are concerned about pressure effects. As defined in the Methods section, for two-dimensional mixtures, the temperature and pressure are in units of ϵ AB k À1 B and ϵ AB σ À2 B , respectively. However, A and B particles also have their own temperature and pressure units, which are ϵ AA k À1 B and ϵ AA σ À2 A for A particles and ϵ BB k À1 B and ϵ BB σ À2 B for B particles, respectively. Therefore, in the mixture at a pressure p in units of ϵ AB σ À2 B , A and B particles effectively feel different pressure values, P A and P B , in their own units: Let us denote T m (p) as the melting temperature of monocomponent systems, i.e., when γ = 1 and Δ = 0 so that A and B particles are identical. In the mixture, since A and B particles feel different pressures in their own units, the corresponding melting temperatures are T m (P A ) and T m (P B ), respectively. In the common units of ϵ AB k À1 B , the melting temperatures that A and B particles feel are thus This leads to a melting temperature gap Apparently, both γ and Δ can affect ΔT m (p). Note that there is a special case: When T m (p) is linear, i.e., T m (p) = Cp with C being the system-dependent coefficient, so the melting temperature gap is no longer a function of Δ. We will see later that such a linear behavior is crucial to the understanding of the pressure and interaction dependence of the GFA. Given T m (p), Eq. (4) indicates that T m,A (p) can be obtained by multiplying the pressure and temperature of T m (p) curve by γ −2 (1 + Δ) and 1 + Δ, respectively. Correspondingly, Eq. (5) shows that T m,B (p) can be obtained by multiplying the pressure and temperature of T m (p) curve by 1 − Δ and 1 − Δ, respectively. For the cases shown in Fig. 1 with γ = 1.4 and Δ = 0, we have T m,A (p) = T m (γ 2 p) and T m,B (p) = T m (p). Fig. 2 shows the melting temperatures against pressure for both harmonic and RLJ systems with γ = 1.4 and Δ = 0. In the log-log scale, T m,A (p) is obtained by just shifting the T m (p) curve horizontally by an amount of log 10 (γ −2 ) = −2log 10 1.4, while T m,B (p) is the same as T m (p). If T m monotonically increases with p, T m,A (p) is always larger than T m,B (p), as shown in Fig. 2b for RLJ systems. If T m (p) is nonmonotonic, T m,A (p) and T m,B (p) may intersect. As shown in Fig. 2a, T m (p) of harmonic systems is non-monotonic in p, so T m,A (p) = T m,B (p) at p ≈ 0.14. In Fig. 1a, we display this pressure as the vertical dot-dashed line. Interestingly, it roughly agrees with p n at which This agreement indicates that our hypothesis that the difference between χ A and χ B , as shown in Fig. 1, is related to the difference between melting temperatures felt by A and B particles is valid. This is further supported by RLJ systems. For RLJ systems, T m,A (p) is always larger than T m,B (p) because T m (p) monotonically increases, consistent with the fact that χ A (p) is always larger than χ B (p), as shown in Fig. 1b.
Although T m,A (p) and T m,B (p) are not the actual equilibrium melting temperatures, the concurrence of T m,A = T m,B and χ A = χ B at p = p n for harmonic systems suggests that ΔT m (p) defined here at least qualitatively reflects the correct pressure evolution of the gap between equilibrium melting temperatures. This can also be seen from the resultant solid-like states visualized in Fig. 1c. For harmonic systems at p > p n , T m,A < T m,B and strong phase separation occurs. From the crystallization mechanism, the liquidus line is in coexistence with the B-solid. When T m,A < T < T m,B , B particles crystallize, while A particles can crystallize until T < T m,A . For phase-separated states obtained from finite  Fig. 1a, it is roughly p n at which χ A = χ B . The inset of a shows how the equilibrium melting temperature T m (p) is determined from simulation, which is also T m,B (p) with Δ = 0 by definition. At a fixed pressure p, the melting temperature T m (marked as the vertical dashed line) is the temperature at which the density ρ undergoes an abrupt change with the decrease of temperature.
quench rates, this sequence of the two-step crystallization would result in purer B-solids than A-solids, i.e., isolated A (B) particles are poor (rich) in B-solids (A-solids), as shown by the snapshots at p = 0.2 in Fig. 1c. The snapshots at other pressures show the opposite behavior when p < p n and T m,A > T m,B .
Seen from Fig. 2, another interesting and important feature is that T m (p)~p at low pressures for both harmonic and RLJ systems. The melting temperature curves deviate from linear at high pressures, where harmonic and RLJ systems exhibit significantly different GFAs. The comparison between Figs. 1 and 2 shows that the GFAs maintain almost constant at fixed dimensionless quench rateκ roughly in the pressure regime where T m (p) is linear. Then the question is whether this is just a coincidence or implies some underlying correlations.
Mixing and demixing for γ = 1 and c B = 0.5. In this work, we focus on the potential energy effects on the GFA. Therefore, before understanding the results of systems with γ = 1.4, c B = 0.5, and Δ = 0, we would like to discuss first a simpler case with γ = 1 and c B = 0.5, so that A and B particles have the same size and geometric frustration induced by particle size difference is absent. Then, it would be clearer to sort out the underlying physics of the evolution from mixing to demixing of two types of particles with only energetic frustration being involved, which provides us with crucial clues to understand the γ = 1.4 case. In the next subsection, we will show that the findings for γ = 1 can be generalized to γ > 1.
With the variation of Δ, A and B particles have different stiffness, leading to energetic frustration affecting the mixing of particles. Because of the trivial symmetry of γ = 1, now we only need to vary Δ from 0 to 1. When Δ = 0, A and B particles are trivially the same and should statistically mix up well. When Δ increases from 0, the distinction in particle stiffness leads to the variation in particle overlap, in analogy to the evolution with the growth of γ 51 . It is thus expected that phase separation will occur when Δ is large. Fig. 3a shows examples of χ B (Δ) of the resultant solid-like states at different pressures and a fixed dimensionless quench rate for harmonic systems. χ A (Δ) behaves similarly to χ B (Δ), so we do not show it here. In this and the next subsections, because we extend the pressure to much lower values, to use the same dimensionless quench ratesκ as in Fig. 1 is far beyond our computational capacity. Therefore, we use fasterκ. We have verified (not shown here) that the results to be presented are reproducible with different values ofκ. At fixed pressure, Fig. 3a shows that there is a crossover Δ c below which χ B remains small and constant. When Δ > Δ c , χ B grows up, so Δ = Δ c signals the onset of phase separation. Fig. 3a also shows that Δ c increases when pressure decreases.
The appearing consistency between the constant GFA at low pressures and the linear T m (p) as discussed in the previous subsection stimulates us to investigate whether the emergence of phase separation at Δ c is also correlated with the linearity of T m (p). respectively, in the temperature-pressure plane with log-log scale. In Fig. 3b, we denote p c as the crossover pressure above which T m (p) becomes nonlinear. Here, we set p c ≈ 0.004 at which T m (p) deviates from linear by 5%. Correspondingly, Eqs. (4) and (5) indicate that T m,A (p) and T m,B (p) become nonlinear at p c,A = p c (1 + Δ) and p c,B = p c (1 − Δ), respectively, when γ = 1. Because Δ > 0, T m,A (p) and T m,B (p) collapse and are both linear when p < p c,B . Seen from Eq. (7), the melting temperature gap ΔT m (p) = 0 when p < p c,B , and ΔT m (p) > 0 otherwise. Therefore, in the pressure regime where p < p c , p = p c,B = p c (1 − Δ † ) sets a crossover below (above) which ΔT m = 0 (ΔT m > 0). In the pressure regime where p > p c , any nonzero Δ leads to ΔT m ≠ 0, so Δ † = 0. Note that both Δ † and Δ c increase when pressure decreases. It is then natural to ask whether they are related. Figure 3c shows Δ † (p) together with the iso-χ B contours for harmonic systems at fixedκ. With the decrease of χ B , the contours approach Δ † (p). The contours may shift upward gradually with the decrease ofκ. With current computational capacity, we expect from Fig. 3c that Δ † agrees with Δ c . Next, we will show that this expectation is valid.
At fixed pressure, the melting temperature gap ΔT m is the function of Δ, as shown by Eq. (6). Therefore, χ B (Δ) shown in Fig. 3a can be converted to χ B (ΔT m ). This functional relation establishes the connection between the GFA characterized by χ B and the melting temperature gap ΔT m for γ = 1. Fig. 3d shows χ B (ΔT m ) curves at different pressures for both harmonic and RLJ systems. At all pressures, χ B has a minimum when ΔT m = 0. With the increase of ΔT m , χ B increases and phase separation emerges. Interestingly, in the pressure regime where p < p c and T m (p) is linear, all χ B (ΔT m ) curves can collapse nicely onto the same master curve, when χ B is plotted against ΔT m /p ν with ν ≈ 1.02 for both harmonic and RLJ repulsions. Figure 3d also shows that the scaling collapse stops working when p > p c , highlighting the important role of the linear T m (p).
The scaling collapse indicates that at different pressures lower than p c , in order to reach the same degree of particle mixing or demixing, the variation of particle stiffness (Δ) needs to cause a melting temperature gap ΔT m proportional to pressure. In this pressure regime, T m,A (p) is always linear, i.e., T m,A (p) = Cp. From Eq. (6), we have when Δ > Δ † and ΔT m > 0. The right hand side of Eq. (9) is constant in pressure for a given χ B . Because T m,B (p) is nonlinear when ΔT m > 0 (so is T m ð p 1ÀΔ Þ on the right hand side of Eq. (9)), p 1ÀΔ must be a constant. This leads to where a is a χ B -dependent prefactor. Eq. (10) sets the pressuredependent values of Δ to reach the same χ B at p < p c . Note that Δ † (p) in Eq. (8) shows exactly the same form of pressure dependence. Therefore, Eq. (8) is actually a direct consequence of the scaling collapse in the ΔT m → 0 limit. On the other hand, Δ c defined above follows Eq. (10) as well, corresponding to the minimum χ B where ΔT m = 0. Therefore, the scaling collapse naturally suggests that Δ † = Δ c . Consequently, for γ = 1, phase separation can occur at a given pressure p only when Δ > Δ † and ΔT m > 0. Equation (8) also indicates that Δ † → 1 when p → 0, so ϵ BB → 0 and B particles cannot feel the existence of each other. In the p → 0 limit, as long as Δ < 1 and there are still repulsions among B particles, A and B particles become identical hard particles when γ = 1 and there will be no phase separation. This is consistent with our ΔT m argument, because it is impossible to cause ΔT m > 0 in the p → 0 limit.
Glass-forming ability for γ > 1 and c B = 0.5. Inspired by the results of γ = 1, in this subsection, we will generalize the connection between χ A (χ B ) and ΔT m to γ > 1. The results of γ = 1 are thus naturally incorporated into the generalized picture. Our arguments will be examined by simulation results of γ = 1.4 and c B = 0.5 without loss of generality. The observations of systems with γ = 1.4, c B = 0.5, and Δ = 0 presented above will then be explained. Here we also vary Δ as done for γ = 1. Because γ > 1 and A and B particles have different sizes, Δ and −Δ apparently correspond to two different systems. We thus need to restore the range of Δ to [ −1, 1].
According to Eq. (6), at fixed pressure p, the melting temperature gap ΔT m varies with Δ, so we are able to establish the functional relation between the GFA and ΔT m . Figure 4a Analogous to γ = 1, Fig. 4c, d show that for both harmonic and RLJ systems χ A (ΔT m ) and χ B (ΔT m ) curves at different pressures lower than p c can also collapse onto the same master curves, when χ A and χ B are plotted against ðΔT m À ΔT Ã m Þ=p ν with ν ≈ 1.02. Therefore, the scaling collapse shown in Fig. 3d for γ = 1 is just a special case with ΔT Ã m ¼ 0. When p > p c , the scaling collapse breaks down.
For γ = 1, we have also shown that there is a range of Δ where Then it is natural to ask whether there also exists a range of Δ within which ΔT m ¼ ΔT Ã m for γ > 1. Another interesting question is what determines ΔT Ã m . Because T m (p) is linear when p < p c , according to Eqs. (6) and (7), at a given p < p c , there indeed exist a range of Δ within which ΔT m is constant, as long as both T m,A (p) and T m,B (p) are still linear. As shown in Fig. 2, γ = 1.4 leads to the separation of T m, A (p) and T m,B (p) when Δ = 0. When Δ varies from 0, T m,A (p) and T m,B (p) curves in Fig. 2 will shift in both horizontal and vertical directions simultaneously by the same amount of log 10 (1 + Δ) and log 10 (1 − Δ), respectively, in the log-log scale. According to Eqs. (4) and (5), at a given p < p c , T m,A (p) or T m,B (p) become nonlinear when Δ is below or above In the p → 0 limit, Δ y 1 and Δ y 2 approach −1 and 1, respectively, consistent with our previous discussions for γ = 1. At low pressures, Δ y 2 is apparently larger than Δ y 1 . When Δ y 1 ≤ Δ ≤ Δ y 2 at a given p, ΔT m remains constant. If such a constant ΔT m is ΔT Ã m , the scaling collapse shown in Fig. 4c, d also validates that Δ y 1 and Δ y 2 are the onsets of the growth of χ A and χ B and the weakening of the GFA, as observed for γ = 1.
With the increase of pressure, Δ y 1 and Δ y 2 approach each other, until arriving at a critical pressure When p > p e , any variation of Δ normally causes the change of ΔT m , so there is no interval of Δ within which ΔT m remains constant. Therefore, the condition p < p c used above should be generally modified to p < p e . When γ = 1, Eqs. (11)- (14) lead to exactly the same results as shown in the previous subsection, with Δ y 2 ¼ ÀΔ y 1 ¼ Δ y , p e = p c , and Δ* = 0. For γ = 1, Δ* = 0 corresponds to the best mixing of particles with ΔT Ã m ¼ 0 at all pressures. Is it possible that Δ* defined in Eq. (14) sets ΔT Ã m for γ > 1 as well?
When plugging in Δ = Δ*, Eqs. (4) and (5) turn to Interestingly, T m,A (p) = γ 2 T m,B (p), so in the log-log scale T m,A (p) is obtained by simply shifting T m,B (p) upward by an amount of log 10 (γ 2 ). Then both T m,A (p) and T m,B (p) become nonlinear when p > p e . In the inset of Fig. 4e, we show T m,A (p) and T m,B (p) for harmonic systems with γ = 1.4 and Δ = Δ*. In the log-log scale, they are parallel along the vertical direction. We calculate the melting temperature gap ΔT m;Δ Ã ðpÞ at Δ = Δ*. Surprisingly, it agrees well with ΔT Ã m ðpÞ over the whole range of pressures studied and for both harmonic and RLJ systems, as shown in Fig. 4e. This agreement confirms our conjecture that there exists a pressure-independent Δ* at which the GFA or the degree of particle mixing is the strongest, equal to that of the hard particle counterpart. When p > p e , any deviation of Δ from Δ* causes ΔT m ≠ΔT Ã m and the weakening of the GFA. When p < p e , Δ* lies in the interval of Δ 2 ½Δ y 1 ; Δ y 2 , within which ΔT m ¼ ΔT Ã m and the GFA can maintain the best. Now we are able to understand the results of systems with γ = 1.4, c B = 0.5, and Δ = 0 discussed earlier. From Eqs. (11) and (12), we can see that Δ = 0 lies in ½Δ y 1 ; Δ y 2 when p < γ −2 p c , so the GFA can maintain the best. When p > γ −2 p c , Δ = 0 lies outside of ½Δ y 1 ; Δ y 2 and is away from Δ * when p > p e . Therefore, ΔT m,0 , the melting temperature gap with Δ = 0, deviates from ΔT Ã m and the GFA becomes worse. How much ΔT m,0 deviates from ΔT Ã m determines how worse the GFA becomes. Fig. 4e also compares ΔT m,0 (p) with ΔT Ã m ðpÞ. At low pressures, they agree very well, as expected. At high pressures, ΔT m,0 deviates from ΔT Ã m . For RLJ systems, the deviation is small, so the GFA is still close to the best. In contrast, for harmonic systems, the deviation significantly increases with the increase of pressure, so the GFA becomes much worse and remarkable phase separation occurs. The more essential reason of such a distinction between harmonic and RLJ systems is that harmonic systems have a much more nonlinear T m (p).
Seen from Fig. 4a, b, the minimum values of χ A and χ B with Δ = Δ* are almost constant in pressure. This explains why the GFA remains constant at low pressures for harmonic and RLJ systems with γ = 1.4, c B = 0.5, and Δ = 0. More interestingly, this suggests that the GFA of hard particle systems can be maintained to high pressures far beyond the hard particle limit, by properly modulating the particle interactions and causing the interplay between geometric and energetic frustrations. To apply Δ* in Eq. (14) is one solution for binary mixtures, which should work generally for other types of particle interactions.
Dependence on c B . In previous subsections, we have focused on c B = 0.5. Fig. 5 verifies that our major findings of the pressure and particle interaction effects on the GFA hold for other values of c B . For given c B , γ, and p, we tune Δ and vary the melting temperature gap ΔT m , exactly as done for c B = 0.5.
In this subsection, we take harmonic systems with γ = 1.4 as the example. As shown in Fig. 5a, b, for c B = 0.3 and 0.7, χ A and χ B reach the minimum at a pressure-dependent ΔT Ã m . Curves of χ A (ΔT m ) and χ B (ΔT m ) at different pressures lower than p e defined in Eq. (13) can collapse well when plotted against ðΔT m À ΔT Ã m Þ=p ν with ν ≈ 1.02. Therefore, as for c B = 0.5, the GFA of the hard particle limit can be sustained with the increase of pressure, as long as ΔT m (p) remains linear. For all values of c B studied, Δ = 0 leads to a ΔT m equal to ΔT Ã m when p < p e , as already discussed for c B = 0.5 in the previous subsection. Therefore, Δ = 0 can maintain the GFA of the hard particle limit until p = p e , above which the GFA becomes weaker. However, when Δ = Δ*, Figs. 4e and 5c show that ΔT m;Δ Ã % ΔT Ã m , so the GFA of the hard particle limit can be maintained over a wide range of pressures and beyond p = p e . Seen from Fig. 5a, b, a remarkable change with c B is the variation of the minimum values of χ A and χ B , reflecting the variation of the GFA of the hard particle limit with c B . The GFA can be approximately quantified by χ = χ A (1 − c B ) + χ B c B . Fig. 5a, b suggest that the GFA of c B = 0.3 is apparently weaker than those of c B = 0.5 and 0.7, while the GFAs of c B = 0.5 and 0.7 are similar with that of c B = 0.5 seeming a little bit stronger. The GFAs in the c B = 0 and 1 limits are trivially the weakest. We also find (not shown here) that the best GFA associated with the minimum χ occurs between c B = 0.5 and 0.7.
The evolution of the GFA with the increase of c B characterized in our work is consistent with previous results of binary mixtures of hard disks with γ = 1.4 from the calculations of the interface energy 23 . Based on our results, we can now qualitatively extend the c B dependence of the GFA to various pressures. In Fig. 5d, we schematically plot the GFA against c B for different pressures and for Δ = 0 and Δ*, respectively. For Δ = 0, it can be expected that the GFA(c B ) curve of the hard particle limit can be maintained until p = p e . When p > p e and the GFA becomes weaker, the whole GFA curve will shift down with the increase of pressure. For Δ = Δ*, we would expect that the GFA(c B ) curve of the hard particle limit can be maintained beyond p = p e .
In the perspective of the empirical argument of the connection between the GFA and eutectic point, the GFA can be evaluated from the depression of the equilibrium melting temperature with the variation of c B 6,23 . However, the argument may not be easily applied to compare the GFAs at different pressures or densities, purely from the comparison of melting temperatures. One may then suggest to perform the calculations of the equilibrium phase diagram and the energy barriers (and the rate) of the nucleation of corresponding crystals to evaluate the GFA, as done for binary mixtures of hard disks 23,50 . Because we have multiple parameters, c B , p, and Δ at fixed γ, it would be a rather difficult task.
Here we realize the comparison of the GFAs over pressures, interactions, and particle compositions in a much simpler way.
The key is the finding of the underlying connection between the GFA and the c B -independent melting temperature gap between species derived from the simple conversion of units. The melting temperatures adopted here are not the actual ones from direct simulations or calculations of equilibrium systems, which should depend on not only p but also c B for given γ and Δ, but they turn out to work well to reveal the pressure and interaction dependence. Another important finding is Δ*, which suggests an energy strategy to fight against the weakening of the GFA caused by the softness of particles.
Although we do not show exact results of the equilibrium phase diagram and melting temperatures in the complicated parameter space, the schematic plots in Fig. 5d already to some extend indirectly reflect the pressure and interaction evolution of the equilibrium phase diagram in the T − c B plane. It would be expected that when p < p e the phase diagrams at various pressures look similar for Δ = 0 and Δ*, only that the melting temperatures grow linearly with the pressure. When p > p e , the phase diagram starts to deviate for Δ = 0, while the similarity may be still maintained for Δ = Δ*.

Discussion
By investigating the pressure and interaction dependence of the GFA, we find similarities and distinctions between two types of widely studied model glass-formers, binary mixtures of harmonic and RLJ particles with γ = 1.4, c B = 0.5, and Δ = 0. We focus on the energetic frustration, taking the hard particle limit as the reference. For a given γ, the GFA is the best in the hard particle limit, purely determined by geometric frustration. The involvement of the potential energy normally weakens the GFA of hard particles. We find that the GFAs of both harmonic and RLJ systems remain constant and identical at low pressures, but bifurcate at high pressures. In contrast to RLJ systems which still maintain an almost constant GFA at high pressures, significant phase separation occurs in harmonic systems.
With the variation of both γ and Δ, we come up with a generalized picture, which explains the pressure and interaction dependence of the GFA. Our major findings include (i) the nonequilibrium GFA of binary mixtures of soft particles is connected to the equilibrium melting temperature, (ii) the GFA of the hard particles can be maintained to high pressures with a proper modulation of the particle interaction such as using Δ* in Eq. (14), and (iii) melting temperatures more linear in pressure are better to suppress phase separation and maintain a good GFA. We find that these results are valid for different values of c B . In combination of the c B dependence studied as well in previous approaches and the pressure and interaction dependence studied in this work, we are able to have a much more comprehensive picture of the GFA of binary mixtures of soft particles. Harmonic and RLJ potentials studied here are rather different, but their behaviors can both be understood by the same picture. Moreover, the derivations in this work do not aim at any specific interaction. Therefore, we expect that our findings are valid to other types of interaction. Note that in this work ϵ AA , ϵ BB , and ϵ AB are connected by the variable Δ and are thus not independent. Therefore, the major conclusions made here may not be directly generalized to other conditions, for example, when all three ϵ's are independence of each other. However, the robust and consistent evidence shown in this work suggests that the underlying connections between the GFA and the melting temperature should not be a coincidence. Follow-up studies are required to find out whether a more general picture can be established.
Since the melting temperature affects the GFA, it is then straightforward to expect that it also plays some role in dynamics of supercooled liquids. It is interesting to know whether A and B particles exhibit different dynamics related to their melting temperature difference. We are also curious about how the dynamics change when Δ varies from 0 to Δ* with the best GFA, from which we may reveal some underlying connections between the GFA and dynamical properties of glass-formers, such as kinetic fragility and dynamic heterogeneity.

Methods
System information. Our systems contain N/2 A and N/2 B particles with the same mass m and a size ratio γ defined earlier. Periodic boundary conditions are applied in all directions. We consider two types of particle interactions, harmonic: and repulsive Lennard-Jones (RLJ): where r ij and σ ij are the separation between particles i and j and sum of their radii, ϵ ij is the characteristic energy scale, and Θ(x) is the Heaviside step function. We set the units of mass, energy, and length to be m, ϵ AB , and σ B , so the units of time and temperature are σ B m 0:5 ϵ À0:5 AB and ϵ AB k À1 B with k B being the Boltzmann constant.
Dimensionless quench rate. In this work, we compare the GFA over a wide range of pressures from 10 −5 to above 0.1 and thus have to confront the challenge to choose a reasonable quench rate. To our knowledge, this issue has not been seriously considered, because people rarely compare the GFAs at different pressures. It has been shown that, when glass-forming liquids are quickly quenched to T = 0, the properties of the resultant inherent structures, e.g., potential energy and stability, depend on the parent temperature T p prior to the quench [52][53][54][55] . For glassforming liquids, there are two characteristic temperatures, the onset temperature T onset and the glass transition temperature T g 40,55 . When T p > T onset and the liquids still exhibit the Arrhenius relaxation behavior, the potential energy of the inherent structures does not vary much with T p . When T g < T p < T onset and the liquids are supercooled and exhibit super-Arrhenius behavior, the potential energy of the inherent structures decreases with the decrease of T p . Previous results have suggested that T onset is around 2T g 40 .
In our study, we start with an equilibrium liquid above T onset and apply a quench rate κ = δT/Δt. Apparently, how long the system stays in the temperature window, (T g , T onset ), is crucial to the structures and properties of the final solid-like states. In contrast, above T onset or below T g , the system is either an equilibrium liquid or a glass with extremely slow structural relaxation, so how long the system stays at T > T onset and T < T g in the accessible time scales should not significantly affect the final solid-like state.
Previous studies have shown that T g~p at low pressures for systems studied in this work with γ = 1.4 and Δ = 0 30 . If we use the same quench rate κ for all pressures, the time for the system to stay in the supercooled regime, (T g , T onset ), is (T onset − T g )/κ~p/κ. Compared with high-pressure systems, the low-pressure systems almost undergo no time in the supercooled regime, which is unfair for them to explore lower-energy inherent structures. Therefore, to compensate this loss at low pressures and give systems at quite different pressures comparable chances to explore lower-energy inherent structures, we need to let κ~p. This leads to an updated quench rate κ Ã ¼ κ=ðpσ d B =k B Þ, where d is the dimension of space. This is actually to nondimensionlize the temperature step δT in the expression of κ by pσ d B =k B . Moreover, it has also been shown that the structural relaxation time τ of the supercooled liquids studied in this work satisfies the scaling at low pressures 30 : is the scaling function. Therefore, in order for systems at different pressures to undergo comparable structural relaxations in the supercooled liquid regime, the quench rate is required to be further divided by p 1/2 . This is actually to nondimensionlize the time duration Δt in the expression of κ by ðpσ dÀ2 B =mÞ À1=2 .
Then we finally obtain a dimensionless quench ratẽ For two-dimensional systems mainly studied in this work,κ ¼ k B m 1=2 σ 2 B κ p 3=2 . By using the sameκ at different pressures, we are able to compare the GFAs of systems with the pressure varying over several orders of magnitude. This is particularly crucial to compare the low-pressure regimes where both the glass transition and melting temperatures are linear in pressure. The robust scaling collapse of χ(ΔT m ) curves strongly validates the use ofκ. Whenκ is fixed, κ~p −3/2 , so the computational cost dramatically increases with the decrease of pressure. As a compromise, we use larger values ofκ when focusing on the low-pressure systems in Figs. 3-5 than those for higher pressures in Fig. 1.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The computer codes of this study are available from the corresponding author upon reasonable request.