# Ultrastable metallic glasses in silico

Anshul D. S. Parmar, Misaki Ozawa, Ludovic Berthier

► To cite this version:

Anshul D. S. Parmar, Misaki Ozawa, Ludovic Berthier. Ultrastable metallic glasses in silico. Physical Review Letters, 2020, pp.085505. 10.1103/PhysRevLett.125.085505. hal-02986302

HAL Id: hal-02986302

https://hal.science/hal-02986302

Submitted on 2 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Ultrastable metallic glasses in silico

Anshul D. S. Parmar,¹ Misaki Ozawa,² and Ludovic Berthier¹,³

¹Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France.

²Laboratoire de Physique Statistique, École Normale Supérieure, CNRS, PSL Research University, Sorbonne Université, 75005 Paris, France

³Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom.

We develop a generic strategy and simple numerical models for multi-component metallic glasses for which the swap Monte Carlo algorithm can produce highly stable equilibrium configurations equivalent to experimental systems cooled more than 10⁷ times slower than in conventional simulations. This paves the way for a deeper understanding of thermodynamic, dynamic, and mechanical properties of metallic glasses. As first applications, we considerably extend configurational entropy measurements down to the experimental glass temperature, and demonstrate a qualitative change of the mechanical response of metallic glasses of increasing stability towards brittleness.

Glasses are obtained by cooling liquids into amorphous solids [1]. This process involves a rapidly growing relaxation time, making it difficult to investigate the nature of the glass transition in equilibrium [2, 3]. Many types of materials can form glassy states, such as molecular, oxide, and colloidal glasses, having various practical applications [4]. Among them, metallic glasses are a promising class known for higher strength and toughness [5], which is vital for applications. Computer simulations represent a valuable tool to investigate glass properties with atomistic resolution [3]. Model metallic glasses are widely used because they are simpler than molecular liquids to understand the basic mechanisms of the glass transition and accompany practical applications. However, typical cooling rates in silico are faster than in the laboratory by 6-8 orders of magnitude. Therefore, computer studies of metallic glasses may produce materials that behave differently from experimental systems. Our goal is to fill this wide gap for metallic glasses in order to access thermodynamic, dynamic, and mechanical properties that can be directly compared to experiments.

Recently, the swap Monte Carlo algorithm has enabled the production of highly stable configurations for models of continuously polydisperse soft and hard spheres [6, 7]. This was achieved by optimising the size distribution and pair interactions to produce good glassformers (preventing crystallisation) with a massive thermalisation speedup [7]. It was however found that previous popular models for metallic glasses, such as the Kob-Andersen [8] and Wahnström mixtures [9], are either not well suited for the swap algorithm [10], or crystallise too easily [7, 11–14]. Further developments are clearly needed.

Here, we show how to develop multi-component metallic glass-formers to benefit from the dramatic speedup offered by swap Monte Carlo, and thus bridge the gap between metallic glass simulations and experiments [15–18]. Our strategy differs from earlier work [7] since it is inspired by the microalloying technique used in metallic glass experiments [19, 20]. We introduce additional

species to the original binary Kob-Andersen mixture to simultaneously improve its glass-forming ability [21, 22] and swap efficiency [7]. This echoes the doping technique widely used in molecular liquids [23–25] to prevent crystallization [19, 20, 23, 24]. The speedup provided by the swap Monte Carlo algorithm depends on the concentration of the doped species. For some models, we can produce for the first time equilibrium configurations of metallic glasses at the experimental glass transition temperature in silico. Our results pave the way for the next generation of thermodynamic and mechanical studies of metallic glasses using computer simulations.

Models—The original Kob-Andersen (KA) model [8] is a 80:20 binary mixture of $N_A$ Lennard-Jones particles of type A, and $N_B$ particles of type B, mimicking the mixture Ni-P. We add a new family of particles, of type C, which can be a single type (ternary mixture) or several types (multi-component). The pair interaction is

$$v_{\alpha_i\beta_j}(r) = 4\epsilon_{\alpha_i\beta_j} \left[ \left( \frac{\sigma_{\alpha_i\beta_j}}{r} \right)^{12} - \left( \frac{\sigma_{\alpha_i\beta_j}}{r} \right)^{6} \right], \quad (1)$$

where $\epsilon$ and $\sigma$ are the energy scale and interaction range, respectively. We specify the particles index by Roman indices and the family type by Greek indices. The potential is truncated and shifted at the cutoff distance $r_{cut,ij} = 2.5\sigma_{\alpha_i\beta_j}$. For particles A and B, we use the interaction parameters of the original KA model: $\epsilon_{AB}/\epsilon_{AA} = 1.5$, $\epsilon_{BB}/\epsilon_{AA} = 0.5$, and $\sigma_{AB}/\sigma_{AA} = 0.8$, $\sigma_{BB}/\sigma_{AA} = 0.88$. Energy and length are in units of $\epsilon_{AA}$ and $\sigma_{AA}$, respectively. Given the large size and energy disparities, performing particles swaps between A and B particles is prohibited [10] which leaves the standard KA model out of the recent swap developments.

We introduce $N_C$ particles of type C. Each C particle is characterized by a continuous variable $\omega_i \in [0, 1]$ so that its interactions with A and B particles are given by

$$X_{AC} = \omega_i X_{AA} + (1 - \omega_i) X_{AB}, \\ X_{BC} = \omega_i X_{AB} + (1 - \omega_i) X_{BB}, \quad (2)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a,b) Two families of Lennard-Jones models, composed of A and B particles interacting as in the KA model, and C particles intermediate between A and B types depending on the variable $\omega$. The C particles improve the glass-forming ability and the efficiency of the swap Monte Carlo algorithm. Relaxation times for the (c) KA$_1$ and (d) KA$_2$ models using full ($\tau_{\alpha}^{(N)}$) and empty ($\tau_{\alpha}^{(S)}$) symbols. The blue box indicates the extrapolated location of the experimental glass transition temperature $T_g^*$, and $T_{mct}^*$ is the mode-coupling crossover.</figcaption></figure>

where $X$ stands for both $\epsilon$ and $\sigma$, so that C particles are identical to A (B) particles when $\omega_i = 1$ (0) and smoothly interpolate between both species for $0 < \omega_i < 1$. Two C particles $i$ and $j$ interact between each other additively:

$$X_{C_i C_j} = \omega_{ij} X_{AA} + (1 - \omega_{ij}) X_{BB}, \quad (3)$$

where $\omega_{ij} = (\omega_i + \omega_j)/2$.

This generic framework offers multiple choices for the distribution of C particles, depending on the parameters $N_C$ and on the chosen distribution $P(\omega)$ of the variable $\omega$. We have explored two simple families, illustrated in Figs. 1(a,b). The first family, KA$_1$, is obtained using a flat distribution $P(\omega)$ on the interval $[0, 1]$, see Fig. 1(a). This corresponds to a multi-component system where C particles continuously interpolate between A and B components. The second family, KA$_2$, is obtained by taking the opposite extreme where $P(\omega) = \delta(\omega - 1/2)$, see Fig. 1(b). In that case, we simulate a discrete ternary mixture. In both cases, we define $\delta = N_C/(N_A + N_B)$ and consider a range of $\delta$ values from $\delta = 0\%$ (original KA mixture), up to $\delta = 20\%$. Contrary to previous work [7], the size dispersity quantified by the variance of the diameter distribution is nearly constant across KA, KA$_1$ and KA$_2$ models. We perform simulations in a periodic cubic cell of volume $V$ in three dimensions. All models are simulated at the number density $\rho = 1.2$, denoting the number of particles in unit volume $\sigma_{AA}^3$.

**Swap Monte Carlo algorithm**-To achieve equilibration

at very low temperatures, we perform Monte Carlo (MC) simulations possessing both translational displacements and particle swaps [26, 27]. For the normal MC moves, a particle is randomly chosen and displaced by a vector randomly drawn within a cube of linear size $\delta r_{max} = 0.15$. The move is accepted according to the Metropolis acceptance rule, enforcing detailed balance. Such MC simulations show quantitative agreement with molecular dynamics simulations in terms of glassy slow dynamics [28].

When using swap MC, we also perform particle swaps. We randomly choose a C particle, say particle $i$, characterised by $\omega_i$. We then randomly choose a value $\Delta\omega$ in the interval $\Delta\omega = \pm 0.8$ and choose a second particle within this interval, say particle $j$. We estimate the energy cost to exchange the type of the two particles, $\omega_i \leftrightarrow \omega_j$, and accept the swap according to the Metropolis rule. In the swap MC scheme, we perform swap moves with probability $p = 0.2$, and translational moves with probability $1-p = 0.8$. All parameters, $(\delta r_{max}, p, \Delta\omega)$ have been carefully optimised to maximise the swap efficiency, see supplementary material (SM [29]). In particular, swaps with larger $\Delta\omega$ are essentially all rejected, confirming that direct A $\leftrightarrow$ B swaps are impossible. In essence, the C particles thus allow two-step exchanges, such as A $\leftrightarrow$ C $\leftrightarrow$ B. Although we only apply this strategy to the KA model, we expect that it should generically apply to high-entropy alloys which have more than five components [30]. In both normal and swap MC schemes, one Monte Carlo time step represents $N$ attempts to make an elementary move. Timescales are reported in this unit.

**Glass-forming ability**-Thanks to modern computer resources, the original KA model is now found to be prone to crystallisation [13, 14, 31]. We have repeated the detailed common neighbor analysis of Ref. [14]. We detected no sign of crystalline environments in our extended models, KA$_1$ and KA$_2$, across the wide temperature regime where thermalisation can be achieved using the swap MC algorithm, see SM. Thus, the extended KA models developed here are much better glass-formers than the original KA model. Similarly to experiments, the doping C particles considerably frustrate the system against crystallisation [19, 20, 23, 24].

**Equilibration speed-up**-The relaxation time $\tau_{\alpha}$ of the system is quantified from the time decay of the self-intermediate scattering function for all particles, $F_s(q, t = \tau_{\alpha}) = 1/e$. We use $q = 7.34$, close to the first diffraction peak of the static structure factor. We respectively denote $\tau_{\alpha}^{(N)}$ and $\tau_{\alpha}^{(S)}$ the relaxation times for the normal (N) and swap (S) dynamics. We finally rescale the relaxation times using its value $\tau_o = \tau_{\alpha}^{(N)}(T = T_o)$ at the onset temperature $T_o$ at which the relaxation time starts to deviate from the Arrhenius law.

We first concentrate on the physical dynamics using normal MC simulations for both models, KA$_1$ and KA$_2$, and various values of $\delta$. We find that the temperature
<table><caption>TABLE I. Characteristics of the KA₁ models. Scaling factor for the temperature ε, the lowest simulated temperature T<sup>*</sup><sub>low</sub> relative to the mode-coupling crossover T<sup>*</sup><sub>mct</sub>, thermalisation speedup, thermodynamic fragility K<sub>T</sub>, and extrapolated Kauzmann temperatures T<sup>*</sup><sub>K</sub>.</caption><thead><tr><th>KA<sub>1</sub></th><th>δ = 0%</th><th>1%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td>ε</td><td>0</td><td>0.01</td><td>0.08</td><td>0.13</td><td>0.25</td></tr><tr><td>T<sup>*</sup><sub>low</sub>/T<sup>*</sup><sub>mct</sub></td><td>0.954</td><td>0.824</td><td>0.787</td><td>0.753</td><td>0.704</td></tr><tr><td>Speedup</td><td>1</td><td>10<sup>2</sup></td><td>6 × 10<sup>3</sup></td><td>8 × 10<sup>4</sup></td><td>2 × 10<sup>7</sup></td></tr><tr><td>K<sub>T</sub></td><td>0.335</td><td>0.334</td><td>0.312</td><td>0.293</td><td>0.274</td></tr><tr><td>T<sup>*</sup><sub>K</sub></td><td>0.252</td><td>0.246</td><td>0.236</td><td>0.224</td><td>0.210</td></tr></tbody></table>

dependence of τ<sup>(N)</sup><sub>α</sub> for all models is very similar, and is only weakly affected by the C particles (see SM). The presence of the C particles changes the energy or temperature scale from the original model. To account for this perturbation and ease the comparison between models, we introduce a rescaled temperature, T* = T(1 + ε(δ)), such that the data τ<sup>(N)</sup><sub>α</sub> versus 1/T* for all models coincide, see Fig. 1(c,d) (see also the τ<sup>(N)</sup><sub>α</sub> versus 1/T plots in SM). The measured ε values reported in Table I are small and compatible with a linear growth, ε(δ) ≈ δ, suggesting that C particles simply act as a linear thermodynamic perturbation (see SM for KA₂). We confirm in SM that the pair structure is also weakly affected. From now on, we use the temperature scale T* and thus, by definition, all models studied in this paper display the same physical (normal MC) dynamics as a function of T*. They have the same reference temperatures as the original KA model: their onset temperature is T<sup>*</sup><sub>o</sub> ≈ 0.7, and the mode-coupling crossover is at T<sup>*</sup><sub>mct</sub> ≈ 0.435. These conventional MC simulations can access τ<sup>(N)</sup><sub>α</sub>/τ<sub>o</sub> ~ 10<sup>4</sup> for N = 10<sup>3</sup> particles, corresponding to the lowest simulated temperature T<sup>*</sup><sub>low</sub> = 0.415 and about 10 days of CPU time. Following earlier work [7], we locate the experimental glass transition temperature T<sup>*</sup><sub>g</sub> by extrapolating the measured dynamical date towards τ<sup>(N)</sup><sub>α</sub>/τ<sub>o</sub> = 10<sup>12</sup> using various functional forms which provide a finite range for its location. We find T<sup>*</sup><sub>g</sub> ∈ [0.3 – 0.345], see Fig. 1, and suggest T<sup>*</sup><sub>g</sub> ≈ 0.3 as our favored estimate obtained using the parabolic law [32].

Our first important achievement follows from the temperature evolution of the relaxation times when using swap MC in Fig. 1. Whereas the original KA model with δ = 0 can be thermalised down to T<sup>*</sup><sub>low</sub> ≈ 0.415, we find that thermalisation is achieved at much lower temperatures as soon as δ > 0, with a speedup that increases continuously and exponentially fast with δ. For an equivalent numerical effort, we find for δ = 1% – 20%, T<sup>*</sup><sub>low</sub> = 0.306 – 0.358 (for KA₁) and T<sup>*</sup><sub>low</sub> = 0.326 – 0.371 (for KA₂). The lowest temperature corresponds to T<sup>*</sup><sub>low</sub> ≈ 0.7T<sup>*</sup><sub>mct</sub> ≈ T<sup>*</sup><sub>g</sub>. Converting these temperatures into timescales, we estimate that the numerical speedup varies from a factor 10<sup>2</sup> for δ = 1%, up to more than 10<sup>7</sup>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Temperature evolution of configurational entropy using variables rescaled at the mode-coupling crossover temperature T<sup>*</sup><sub>mct</sub>. The horizontal lines correspond to the lowest accessible temperature for the original KA model and the δ = 20% KA₁ model.</figcaption></figure>

for δ = 20%. Thus, even a small amount of doping has a massive impact on the swap efficiency. The proposed metallic glass models considerably widen the accessible temperature regime available to computer simulations, without suffering from crystallizations.

**Configurational entropy**–We now characterise the configurational entropy, s<sub>c</sub>(T*), of the very low temperature states produced with swap MC. We determine the configurational entropy from its conventional definition, s<sub>c</sub>(T*) = s<sub>tot</sub>(T*) – s<sub>vib</sub>(T*) [33–36]. The equilibrium entropy, s<sub>tot</sub>, is straightforwardly measured by thermodynamic integration from the ideal gas to the studied state point [35]. The vibrational entropy, s<sub>vib</sub>, is obtained by a constrained Frenkel-Ladd [37] thermodynamic integration, generalised to properly quantify the mixing entropy contribution to the vibrational entropy [38]. This is a crucial point for the present models where polydispersity changes continuously with δ, and alternate approaches, for instance, using inherent structures, would be inadequate [39]. Figure 2 shows the temperature evolution of the configurational entropy of KA₁ models. We use T<sup>*</sup><sub>mct</sub> as a useful temperature scale to normalise both temperatures, T*/T<sup>*</sup><sub>mct</sub>, and entropies, s<sub>c</sub>(T*)/s<sub>c</sub>(T<sup>*</sup><sub>mct</sub>), in the spirit of Kauzmann [33]. Our data for the KA model are consistent with previous work [40], and stop at T<sup>*</sup><sub>low</sub>/T<sup>*</sup><sub>mct</sub> ≈ 0.954 where s<sub>c</sub>(T<sup>*</sup><sub>low</sub>)/s<sub>c</sub>(T<sup>*</sup><sub>mct</sub>) ≈ 0.93, corresponding to the deepest states accessible with a conventional MC algorithm. Figure 2 shows that the thermalisation speedup obtained by increasing δ is accompanied by a strong reduction of the configurational entropy. Hence, the deeply supercooled states obtained using swap MC correspond to state points where much fewer amorphous packings are available to the system, which should trans-
late into a larger point-to-set correlation length [41–43]. In earlier studies of the KA model, the putative Kauzmann transition was determined by fitting the decrease of $s_c$ with the empirical form $s_c = A(1-T_K^*/T^*)$, which also allows the determination of the thermodynamic fragility: $K_T \equiv AT_K^*$ [44]. We extend this analysis to KA₁ models and report $T_K$ and $K_T$ in Table I. Both quantities show a modest, but systematic decrease with $\delta$. Remarkably, increasing the studied time window from 4 to 11 orders of magnitude (from $\delta = 0\%$ to 20%), the steep temperature dependence of $s_c$ remains consistent with an entropy crisis taking place at $T_K \approx 0.5T_{mct}^* > 0$. In particular, we detect no sign of a new mechanism to ‘avoid’ it [45]. These data also contradict the arguments that models where swap MC works well are qualitatively distinct from those where it does not [46–48], and constitute our second important result.

**Brittle yielding**–Turning to rheology, we demonstrate that accessing highly stable glassy configurations qualitatively affects how simulated metallic glasses yield. It was recently suggested that glass stability induces a ductile-to-brittle transition, confirmed numerically in a model for soft glasses [49–51]. Here, we establish that a similar transition exists also in metallic glasses. To this end, we consider a larger system size, $N = 5 \times 10^4$, and apply the following preparation protocol for the original KA model, and the $\delta = 1\%$ KA₁ and KA₂ models. First, we thermalise the system at high temperature, $T^* = 2.0$. Second, we instantaneously quench to the temperature $T^* = 0.373$ and 0.319 for KA and KA₁,₂, respectively, where $\tau_\alpha^{(S)} \simeq 10^{10}$ MC steps, and let them age during $10^6$ swap MC steps. We expect to produce an ordinary computer glass of modest stability for the KA model, but very stable configurations for KA₁ and KA₂ models. These aged glasses are quenched to $T = 0$, and sheared using a strain-controlled athermal quasi-static protocol [52]. We apply a uniform shear along the xy-plane, with strain increments $\Delta\gamma = 10^{-4}$. We measure the xy-component of the shear stress, $\sigma_{xy}$, to obtain the stress-strain curves shown in Fig. 3(a). We visualise non-affine particle displacements [53] in Fig. 3(b-d).

For the KA model, the stress shows an initial quasi-linear increase with small plastic events, a stress overshoot punctuated by many larger plastic events, and a gradual approach to steady-state. Near yielding, plasticity is spatially heterogeneous but spreads over the entire system, see Fig. 3(b), in agreement with previous findings [54–56]. For stable initial configurations, the stress overshoot transforms into a unique, sharp, macroscopic stress discontinuity, see Fig. 3(a). This brittle behavior is accompanied by a clear system-spanning shear-band, see Fig. 3(c,d). The tendency to shear localisation upon increasing stability is well-documented [57, 58], but a genuine non-equilibrium discontinuous yielding transition only occurs for highly stable glassy systems [49], see

<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Stress-strain curves for three models, KA, KA₁, and KA₂ with $\delta = 1\%$. For each model, three individual realizations with $N = 50000$ are shown. (b,c,d) Snapshots of non-affine displacement between $\gamma = 0$ and $\gamma = 0.11$ (vertical arrow in (a)) for KA (b), KA₁ (c), and KA₂ (d).</figcaption></figure>

also Refs. [59–61]. These results extend to an experimentally relevant class of materials the observation that brittle yielding and macroscopic shear-band formation can be studied in atomistic simulations. They also suggest the universality of the random critical point controlling the brittle to ductile transition [49]. This constitutes our third important result.

**Perspectives**–The multi-component models for metallic glasses developed here can be efficiently thermalised via swap Monte Carlo simulations down to temperatures that are not currently accessible to conventional simulation techniques and are comparable to the experimental glass transition. These models fill the gap between experimental and numerical works. Considering the extensive use made of the KA model [8], the improved glass-forming ability and thermalisation efficiency will stimulate many future studies. Immediate applications concern further analysis of thermodynamic, dynamical, and mechanical properties of the stable configurations obtained here, to address questions regarding the Kauzmann temperature, the validity of the Adam-Gibbs relation (see SM for an initial attempt) and a microscopic description of shear band formation and failure in metallic glasses. More broadly, the strategy proposed here is simple and versatile, and can certainly be improved further. For example, increasing the number of components and their
concentration allows the system to reach below $T_g$ very efficiently (see SM for a four-component model). It could also be used to model some specific multi-component materials (for instance of the Ni-Pd-P type) and high-entropy alloys, to deepen our theoretical understanding of metallic glasses and help the design of novel materials with specific properties.

We thank D. Coslovich and A. Ninarello for discussions and sharing data and codes. We also thank C. Cammarota, A. Liu and S. Sastry for useful discussions. This work was supported by a grant from the Simons Foundation (#454933, L. Berthier).

[1] C. A. Angell, Science **267**, 1924 (1995).

[2] A. Cavagna, Physics Reports **476**, 51 (2009).

[3] L. Berthier and G. Biroli, Rev. Mod. Phys **83**, 587 (2011).

[4] L. Berthier and M. D. Ediger, Physics Today **69**, 40 (2016).

[5] A. Greer, Y. Cheng, and E. Ma, Mater. Sci. Eng. R Rep **74**, 71 (2013).

[6] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Phys. Rev. Lett **116**, 238002 (2016).

[7] A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X **7**, 021039 (2017).

[8] W. Kob and H. C. Andersen, Phys. Rev. E **51**, 4626 (1995).

[9] G. Wahnström, Phys. Rev. A **44**, 3752 (1991).

[10] E. Flenner and G. Szamel, Phys. Rev. E **73**, 061505 (2006).

[11] Y. Brumer and D. R. Reichman, J. Phys. Chem. B **108**, 6832 (2004).

[12] R. Gutiérrez, S. Karmakar, Y. G. Pollack, and I. Procaccia, EPL **111**, 56009 (2015).

[13] T. S. Ingebrigtsen, J. C. Dyre, T. B. Schrøder, and C. P. Royall, Phys. Rev. X **9**, 031016 (2019).

[14] D. Coslovich, M. Ozawa, and W. Kob, Eur. Phys. J. E **41**, 62 (2018).

[15] H.-B. Yu, Y. Luo, and K. Samwer, Advanced Materials **25**, 5904 (2013).

[16] D. P. Aji, A. Hirata, F. Zhu, L. Pan, K. M. Reddy, S. Song, Y. Liu, T. Fujita, S. Kohara, and M. Chen, arXiv preprint arXiv:1306.1575 (2013).

[17] P. Luo, C. Cao, F. Zhu, Y. Lv, Y. Liu, P. Wen, H. Bai, G. Vaughan, M. Di Michiel, B. Ruta, et al., Nature communications **9**, 1 (2018).

[18] T. Dziuba, Y. Luo, and K. Samwer, Journal of Physics: Condensed Matter (2020).

[19] W. H. Wang, Prog. Mater. Sci. **52**, 540 (2007).

[20] S. González, J. Mater. Res. **31**, 76 (2016).

[21] C. Tang and P. Harrowell, Nature materials **12**, 507 (2013).

[22] K. Zhang, M. Wang, S. Papanikolaou, Y. Liu, J. Schroers, M. D. Shattuck, and C. S. O'Hern, The Journal of chemical physics **139**, 124503 (2013).

[23] C. Angell and D. Smith, J. Phys. Chem. **86**, 3845 (1982).

[24] K. Takeda, O. Yamamuro, I. Tsukushi, T. Matsuo, and H. Suga, J. Mol. Struct **479**, 227 (1999).

[25] S. Tatsumi, S. Aso, and O. Yamamuro, Phys. Rev. Lett **109**, 045701 (2012).

[26] D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Vol. 1 (Elsevier, 2001).

[27] T. S. Grigera and G. Parisi, Phys. Rev. E **63**, 045102 (2001).

[28] L. Berthier and W. Kob, J. Phys. Condens. Matter **19**, 205130 (2007).

[29] See Supplemental Material at [url] for details about simulation methods, temperature scaling, determination of the experimental glass transition temperature, thermodynamics and its relation with dynamics, glass-forming ability, and an additional model system, which includes Refs. [62–70].

[30] Y. Zhang, T. T. Zuo, Z. Tang, M. C. Gao, K. A. Dahmen, P. K. Liaw, and Z. P. Lu, Prog. Mater. Sci. **61**, 1 (2014).

[31] S. Toxvaerd, U. R. Pedersen, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. **130**, 224501 (2009).

[32] Y. S. Elmatad, D. Chandler, and J. P. Garrahan, J. Phys. Chem. B **114**, 17113 (2010).

[33] W. Kauzmann, Chem. Rev. **43**, 219 (1948).

[34] F. Sciortino, W. Kob, and P. Tartaglia, Phys. Rev. Lett **83**, 3214 (1999).

[35] S. Sastry, J. Phys. Condens. Matter **12**, 6515 (2000).

[36] L. Berthier, M. Ozawa, and C. Scalliet, J. Chem. Phys. **150**, 160902 (2019).

[37] D. Frenkel and A. J. Ladd, J. Chem. Phys. **81**, 3188 (1984).

[38] M. Ozawa, G. Parisi, and L. Berthier, J. Chem. Phys. **149**, 154501 (2018).

[39] M. Ozawa and L. Berthier, J. Chem. Phys. **146**, 014502 (2017).

[40] A. Banerjee, M. K. Nandi, S. Sastry, and S. M. Bhattacharyya, J. Chem. Phys. **145**, 034502 (2016).

[41] J.-P. Bouchaud and G. Biroli, J. Chem. Phys. **121**, 7347 (2004).

[42] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, PNAS **114**, 11356 (2017).

[43] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Nat. Commun. **10**, 1 (2019).

[44] S. Sastry, Nature **409**, 164 (2001).

[45] C. P. Royall, F. Turci, S. Tatsumi, J. Russo, and J. Robinson, Journal of Physics: Condensed Matter **30**, 363001 (2018).

[46] M. Wyart and M. E. Cates, Phys. Rev. Lett **119**, 195501 (2017).

[47] L. Berthier, G. Biroli, J.-P. Bouchaud, and G. Tarjus, J. Chem. Phys. **150**, 094501 (2019).

[48] V. Lubchenko and P. G. Wolynes, The Journal of Physical Chemistry B **122**, 3280 (2017).

[49] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, PNAS **115**, 6656 (2018).

[50] M. Ozawa, L. Berthier, G. Biroli, and G. Tarjus, Physical Review Research **2**, 023203 (2020).

[51] W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Physical Review Letters **124**, 225502 (2020).

[52] C. E. Maloney and A. Lemaître, Phys. Rev. E **74**, 016118 (2006).

[53] M. L. Falk and J. S. Langer, Phys. Rev. E **57**, 7192 (1998).

[54] M. Utz, P. G. Debenedetti, and F. H. Stillinger, Phys. Rev. Lett **84**, 1471 (2000).

[55] M. Fan, M. Wang, K. Zhang, Y. Liu, J. Schroers, M. D. Shattuck, and C. S. O'Hern, Phys. Rev. E **95**, 022611 (2017).
[56] D. Rodney, A. Tanguy, and D. Vandembroucq, Model. Simul. Mater. Sci. Eng **19**, 083001 (2011).

[57] F. Varnik, L. Bocquet, J.-L. Barrat, and L. Berthier, Phys. Rev. Lett **90**, 095702 (2003).

[58] Y. Shi and M. L. Falk, Phys. Rev. Lett **95**, 095502 (2005).

[59] J. Ketkaew, W. Chen, H. Wang, A. Datye, M. Fan, G. Pereira, U. D. Schwarz, Z. Liu, R. Yamada, W. Dmowski, et al., Nature communications **9**, 1 (2018).

[60] G. Kapteijns, W. Ji, C. Brito, M. Wyart, and E. Lerner, Phys. Rev. E **99**, 012106 (2019).

[61] H. Bhaumik, G. Foffi, and S. Sastry, arXiv preprint arXiv:1911.12957 (2019).

[62] S. Sastry, Phys. Rev. Lett **85**, 590 (2000).

[63] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Physical Review E **56**, 5533 (1997).

[64] J. D. Ferry, *Viscoelastic properties of polymers* (John Wiley & Sons, 1980).

[65] M. Ozawa, C. Scalliet, A. Ninarello, and L. Berthier, J. Chem. Phys. **151**, 084504 (2019).

[66] M. P. Allen and D. J. Tildesley, *Computer simulation of liquids* (Oxford university press, 2017).

[67] G. Adam and J. H. Gibbs, J. Chem. Phys. **43**, 139 (1965).

[68] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev. A **40**, 1045 (1989).

[69] V. Lubchenko and P. G. Wolynes, Annu. Rev. Phys. Chem. **58**, 235 (2007).

[70] J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. **91**, 4950 (1987).

## SUPPLEMENTARY MATERIAL

In the supplementary material, we provide additional information regarding the following aspects: (i) details of the simulation methods, (ii) comparison of normal dynamics and static correlation function for various $\delta$ values, (iii) estimating the experimental glass transition temperature, (iv) details of the computation of entropies and $KA_2$ model, (v) test of the Adam-Gibbs relation, (vi) robustness against crystallization.

## SIMULATION METHODS

Translational Monte Carlo moves are conceptually straightforward and require only optimisation of the time decay of the structural relaxation time [26]. We find that picking a random displacement within a cube of linear length $\delta r_{max} = 0.15$ with the Metropolis acceptance criteria makes the dynamics optimal. These dynamics can be viewed as equivalent to a discretised Brownian dynamics. One Monte-Carlo sweep consists of $N$ such elementary trial moves, and timescales are reported in this unit [28].

In addition to the conventional Monte-Carlo moves, the swap Monte Carlo approach consists of adding particle swaps, where the value of $\omega$ of a randomly drawn pair of particles is exchanged. We mix translational and swap move with probability $p$ and $1-p$, respectively. As shown in Fig. 4(a), the optimal choice is near $p \approx 0.2$, for which structural relaxation decays the fastest [7].

To further optimize the swap Monte-Carlo moves, we measure the acceptance probability $P(\Delta\omega)$ as a function of the disparity of particle pairs quantified by $\Delta\omega = |\omega_i - \omega_j|$. We consider the temperature $T = 0.50$ and $\delta$ from 1% to 20%. In Fig. 4(b), we see that acceptance is larger when $\Delta\omega$ is small and particles thus strongly resemble each other. The acceptance fall to $10^{-5}$ when $\Delta\omega \approx 0.8$, showing that direct exchanges between A and B particles are strongly suppressed. The trade-off is that larger $\Delta\omega$ more efficiently thermalise the system, but are less frequently accepted. We fix a threshold $\Delta\omega = 0.8$ and do not attempt swapping particles for larger values. For a given C particle, we pick a direction for $\pm 0.80$ with equal chance, and randomly pick another particle can lie beyond $\omega = 1(A)$ or $\omega = 0(B)$ range as well, in that case, we randomly pick particle of type A or B.

Dynamics is characterized by the alpha-relaxation times obtained by the self-intermediate scattering function. Since the particle's type changes throughout the swap simulations, the relaxation time is calculated for the complete system, as the time at which the self-intermediate scattering function $F_s(q, t)$, for the first peak of static structure factor, decays to a value of $1/e$.
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) The relaxation time $\tau_{\alpha}^{(S)}$ of the KA$_1$ model as a function of the probability of swap moves $p$. (b) The swap acceptance rate, $P(\Delta\omega)$, as a function of particle disparity $\Delta\omega = |\omega_i - \omega_j|$ between the particle pairs for which the swap is attempted.</figcaption></figure>

$F_s(q, t)$ is calculated using

$$F_s(q, t) = \frac{1}{N} \left\langle \sum_{i=1}^{N} \exp \left[ -i \mathbf{q} \cdot (\mathbf{r}_i(t) - \mathbf{r}_i(0)) \right] \right\rangle, \quad (4)$$

where $\mathbf{r}_i(t)$ are the positions of particle $i$ at time $t$ and the averaging is performed over many times origins.

## TEMPERATURE SCALING FOR VARIOUS MODELS

The introduction of additional C particles represents a small thermodynamic perturbation to the original model, even if we fix the number density to $\rho = 1.2$ for all $\delta$ values. This can be seen in Fig. 5 where we show the evolution of $\tau_{\alpha}^{(N)}$ for the normal dynamics as a function of $1/T$. A systematic shift of the data with $\delta$ is clearly visible.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Relaxation times for the KA$_1$ model (a) and KA$_2$ model (b) using either the temperature, $T$ (left panels) or rescaled temperature $T^*$ (right panels).</figcaption></figure>

However, by scaling the temperature with a single constant, $T^* = T(1 + \epsilon(\delta))$, the dynamics for all models can be rescaled on a single master curve, as shown in Fig. 5. The scaling works for both KA$_1$ and KA$_2$ models. The temperature $T^*$ is thus useful to compare different systems between each other. By definition $T^* = T$ for the original KA model with $\delta = 0\%$, so that the relevant temperature scales for the KA model directly apply to KA$_1$ and KA$_2$ models.

The temperature shift quantified by $\epsilon(\delta)$ is largely due to the thermodynamic perturbation induced by the doping C particles. To show this, we present in Fig. 6(a) the evolution of the radial distribution function for the majority species, $g_{AA}(r)$ at a given temperature, $T = 0.44$ and increasing $\delta$ values. Clearly, the pair correlation function is weakly affected by the addition of C particles, in a way that increases with $\delta$. Rescaling the temperature to make the dynamics coincide removes a large part of this change, but not all of it. This is shown in Fig. 6(b) where the data for $g_{AA}(r)$ are now collected for a given value of $T^* = 0.44$. The fact that a small change of the pair correlation survives the temperature rescaling shows that despite the models have a very similar dynamics, they can still be distinguished even at the level of the pair structure. They are, therefore, not strictly identical. The improved glass-forming ability (see below) would stem from this slight difference of structure.

It has been reported that the standard KA mixture shows a gas-liquid (or gas-glass) phase separation at a lower density and lower temperature due to the nature of the attractive potential, which is characterized
<figure><img src="image_7.png" /><figcaption>FIG. 6. Pair correlation function for the majority species A for the KA₁ models at (a) T = 0.44, (b) T* = 0.44.</figcaption></figure>

by the negative slope in the isothermal pressure-density curve [62, 63]. This instability appears below the so-called Sastry density that locates the minimum of the pressure-density curve. Since the Sastry density increases with decreasing temperature, we have to care about the instability for our KA₁ and KA₁ models. We have confirmed that our models do not show such instability in both equilibrium and inherent structures, having the Sastry density below ρ = 1.2. This issue could easily be avoided by simulating these models at higher densities. We choose ρ = 1.2 to compare our results with previous studies of the original KA model.

## THE EXPERIMENTAL GLASS TRANSITION TEMPERATURE

We define the laboratory glass transition as $\tau_{\alpha}^{(N)}(T_g^*)/\tau_o = 10^{12}$. Since the standard MC can only provide relaxation times at most of the order of $\tau_{\alpha}/\tau_o \sim 10^4$, where $\tau_o(T = 0.70) \sim 3 \times 10^3$ MC sweeps is the value of structural relaxation time at the onset of slow dynamics. Considering this computational limitation, we need to extrapolate the dynamic behaviour to locate the experimental glass transition temperature $T_g^*$.

We estimate the onset temperature as the deviation from the Arrhenius form and perform a detailed study in the supercooled regime using various functional fit

forms to locate the range of the glass transition temperature [32, 42, 64]. The first functional form is the Vogel-Fulcher-Tammann (VFT) expression [64],

$$ \tau_{\alpha}(T) = \tau_{\infty} \exp \left[ \frac{1}{K_{VFT} [T/T_{VFT} - 1]} \right], \quad (5) $$

where $\tau_{\infty}$, $K_{VFT}$ and $T_{VFT}$ represent high-temperature relaxation time, the dynamic fragility and the finite temperature divergence, respectively. A second functional form is the parabolic fit [32],

$$ \tau_{\alpha}(T) = \tau_o \exp \left[ J \left( \frac{1}{T} - \frac{1}{T_o} \right)^2 \right], \quad (6) $$

where $T_o$ represents the onset temperature. In contrast to the VFT law, this parabolic fit form does not predict a divergence of the relaxation time at any finite temperature. The parabolic law can be seen as the simplest correction to the Arrhenius behaviour that includes dynamic fragility.

As a result, the VFT presumably overestimates the extrapolated relaxation time, and the parabolic law is presumably a safer extrapolation. Indeed a recent study suggests that the parabolic fit form is more reliable and consistent with the experimental data when extrapolated from a small range in the computational domain [65]. In the main text, we mark the laboratory glass transition range from the VFT and parabolic fits obtained below $T_o$. These values are used as boundaries for the possible location of $T_g$.

## THERMODYNAMICS

We provide the details of the calculation of the configurational entropy, $s_c$, which is estimated by the computation of the total, vibrational, and mixing entropies [38].

### Total entropy

The total free energy $A(\rho, T)$ of the system at a density $\rho$ and temperature $T$ can be written as the sum of the ideal gas part $A_{id}(\rho, T)$ and the excess free energy $A_{ex}(\rho, T)$. The ideal gas part $A_{id}(\rho, T)$ can be expressed as

$$ \beta A_{id}(\rho, T) = N(3 \ln \Lambda + \ln \rho - 1). \quad (7) $$

where $\Lambda$ is the thermal wavelength. The excess free energy $A_{ex}(\rho, T)$ is estimated by thermodynamic integration from a known limit [35, 66]. We consider the ideal gas as this reference state point. The computation of the excess free energy can be performed in two steps.

**Step I:** We integrate the excess pressure $P_{ex}$ from the dilute limit to the target density $\rho = 1.20$ at constant
temperature $T_r = 5.00$:

$$ \beta_r A_{ex}(\rho, T_r) = \beta_r A_{ex}(0, T_r) + N \int_0^\rho \frac{d\rho'}{\rho'} \left( \frac{\beta_r P}{\rho'} - 1 \right). \quad (8) $$

The excess free energy of the reference state $A_{ex}(\rho, T)$ contains a combinatorial term resulting from the distinct particle types. This term is the same as the mixing part in the vibrational entropy and cancels out from the configurational entropy.

**Step II:** The excess free energy at the desired temperature $T$, $A_{ex}(\rho, T)$, is calculated by integrating the average potential energy $U$ from temperature $T_r$ to $T$ [35]:

$$ \beta A_{ex}(\rho, T) = \beta_r A_{ex}(\rho, T_r) + \int_{\beta_r}^{\beta} d\beta' U(\rho, \beta'). \quad (9) $$

The total entropy $s_{tot}$ at the target state point is obtained using the relation $s_{tot} = \frac{1}{N} \frac{\delta A}{\delta T}$.

## Vibrational entropy

The vibrational entropy is estimated by employing the modified Frenkel-Ladd method [38], where the Einstein solid is used as a reference state to perform a constrained thermodynamic integration. The constrained Hamiltonian is

$$ \beta U_{\alpha}(r^{N}, r_{0}^{N}) = \beta U_{0}(r^{N}) + \alpha_{s} \sum_{i=1}^{N} |r_{i} - r_{i,0}|^{2}, $$

where $r_0^N$ is the reference equilibrium configuration. In the limit of large stiffness (say, $\alpha_{s,max}$) the system behaves as a classical non-interacting ensemble of harmonic oscillators. For sufficiently small values (say, $\alpha_{s,min} \to 0$) the liquid state is restored. The resulting vibrational entropy can be expressed as [38]

$$ s_v = (3/2) - 3 \ln \Lambda - (3/2) \ln(\alpha_{s,max}/\pi) + \lim_{\alpha_{s,min} \to 0} \int_{\alpha_{s,min}}^{\alpha_{s,max}} d\alpha_s \Delta_{\alpha_s}^{T,S} + s_{mix}^{(M)} - s_{mix}(r_0^N, \beta). \quad (10) $$

Here, $s_{mix}^{(M)}$ is the mixing entropy stemming from the combinatorial factor of the different particle type which cancels out from the ideal gas contribution. The term $\Delta_{\alpha}^{T,S}$ is the mean-squared displacement defined by

$$ \Delta_{\alpha_s}^{T,S} = \frac{1}{N} \left\langle \sum_{i=1}^{N} |r_i - r_{i,0}|^2 \right\rangle_{\alpha_s}^{T,S}, \quad (11) $$

where the upper scripts $T, S$ indicate both translational and swap MC displacements should be performed. Discretisation of the integral over $\alpha$ and boundary values are chosen as in Ref. [38].

<table><caption>TABLE II. Characteristics of the KA₂ models. Scaling factor for the temperature $\epsilon$, the lowest simulated temperature $T_{low}^*$ relative to the mode-coupling crossover $T_{mct}^*$, the factor of thermalisation speedup, the thermodynamic fragility $K_T$, and the extrapolated Kauzmann temperatures $T_K^*$.</caption><thead><tr><th>KA<sub>2</sub></th><th>δ = 0%</th><th>1%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td>ε</td><td>0</td><td>0.01</td><td>0.08</td><td>0.15</td><td>0.27</td></tr><tr><td>T<sup>*</sup><sub>Low</sub>/T<sup>*</sup><sub>mct</sub></td><td>0.954</td><td>0.852</td><td>0.807</td><td>0.777</td><td>0.750</td></tr><tr><td>Speedup</td><td>1</td><td>10<sup>2</sup></td><td>2 × 10<sup>3</sup></td><td>3 × 10<sup>4</sup></td><td>2 × 10<sup>5</sup></td></tr><tr><td>K<sub>T</sub></td><td>0.335</td><td>0.329</td><td>0.316</td><td>0.293</td><td>0.271</td></tr><tr><td>T<sup>*</sup><sub>K</sub></td><td>0.252</td><td>0.249</td><td>0.235</td><td>0.219</td><td>0.206</td></tr></tbody></table>

## Mixing entropy

The mixing entropy $s_{mix}(r_0^N, \beta)$ in Eq. (10) needs a separate estimate [38]. We perform a thermodynamic integration over a temperature range $\beta'$ from the target temperature $\beta = 1/T$ with a given reference configuration to the high temperature limit, $\beta' \to 0$. For a given configuration $(r_0^N)$, the mixing entropy can be expressed as

$$ s_{mix}(r_0^N, \beta) = \frac{1}{N} \int_{0}^{\beta} d\beta' \Delta U_{mix}(r_0^N, \beta'), \quad (12) $$

where $\Delta U_{mix}(r_0^N, \beta) = \langle U_0(r_0^N) \rangle_{\beta'}^S - \langle U_0(r_0^N) \rangle$ corresponds to the potential energy difference between the reference sample and the heated sample. In Eq. (12), the system is heated from the initial state point $\beta$ to $\beta' \to 0$, but only particle permutations are performed while the particle positions $r_0^N$ are unchanged.

## Configurational entropy

Collecting all terms, the configurational entropy is finally obtained as

$$ s_c(T^*) = s_{tot}(T^*) - s_{vib}(T^*) \quad (13) $$

so that $s_{vib} \equiv s_v - s_{mix}$. We show in Fig. 7 the temperature dependence of the various contributions to the configurational entropy as a function of $T^*$ and different values of $\delta$.

## Characteristics of KA₂ models

The table II provides additional details for the characterization of the KA₂ models.

## ADAM-GIBBS RELATION

We attempt to test thermodynamic theories of the glass transition over the experimentally relevant temper-
<figure><img src="image_8.png" /><figcaption>FIG. 7. Temperature evolution of different contributions to the configurational entropy. (a) Total entropy $s_{tot}$. (b) Vibrational entropy $s_v$. (c) Mixing entropy. The horizontal dashed line corresponds to the combinatorial mixing entropy for the $\delta = 0\%$ KA binary mixture. (d) Configurational entropy for the KA$_1$ and KA$_2$ models, lines are fits to the KA$_1$ models to estimate various thermodynamics properties. The blue area marks the laboratory glass transition $T_g^*$.</figcaption></figure>

ature regime. The Adam-Gibbs relation is a simple connection between thermodynamics and dynamics [67]. In the same vein, the random first-order transition (RFOT) theory emphasizes the role of configurational entropy and provides a generalised connection between thermodynamics and dynamics [41, 68, 69].

The RFOT theory version of the Adam-Gibbs relation can be expressed as

$$ \log(\tau_{\alpha}/\tau_{0}) \propto \frac{1}{Ts_{c}^{\alpha}}, \qquad (14) $$

where $\alpha = 1$ restores the original AG relation.

Recently, an extensive study of a range of simulation and experimental measurements suggested a possible modification of the initially proposed Adam-Gibbs relation [65]. We follow this recent analysis for the KA$_1$ models. Our results are in Fig. 8.

In Fig. 8(a), we show the relaxation times $\tau_{\alpha}^{(N)}$ combining direct measurements and extrapolations using the parabolic law, against the measured value of the configurational entropy $s_c$.

The representation shows $\log_{10} \tau_{\alpha}$ versus $1/(T^*s_c)$ such that the original Adam-Gibbs relation would predict a straight line.

The inefficient sampling and possible crystallization of the standard KA model do not allow us to perform a detailed study of the Adam Gibbs relation at very low temperatures, and the relation seems to hold over a relatively narrow dynamic range.

Instead, the KA$_1$ models can be equilibrated to much lower temperatures because the swap MC algorithm is more efficient for them. The data in Fig. 8(a) suggest that when considered over a broader dynamic range, systematic deviations from the Adam Gibbs relation become clearly visible.

As found before for other models [65], we find these small deviations from the Adam Gibbs relation can be accounted for using an exponent $\alpha < 1$ in Eq. (14). In Fig. 8(b), we show that an exponent $\alpha \approx 0.64$ actually describes the behaviour of all models from $\delta = 0\%$ (the original KA model) up to $\delta = 20\%$. Values of $\alpha$ sys-
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a) Adam-Gibbs plot for various $\delta$. The dynamics combines data measured and extrapolated using the parabolic law. The original Adam-Gibbs relation corresponds to a straight line. (b) Generalized Adam Gibbs plot using $\alpha = 0.64$ for all $\delta$ values.</figcaption></figure>

tematically smaller than 1 were also reported in many materials [65].

## GLASS FORMING ABILITY

Although considered a good glass-former for many years, progress in computer resources has led to the conclusion that the model can be relatively easily crystallised at low enough temperatures [13, 14, 31]. The KA mixture first demixes, and the large A particles partially crystallize into a face-centered cubic (FCC) structure [13, 14].

As an indicator of crystallization events in the simulations, we focus on the common neighbor analysis (CNA).

<figure><img src="image_10.png" /><figcaption>FIG. 9. The fraction of CNA-142 bonds, for the inherent structures, as a function of PT MD steps at $T = 0.373$ and $N = 1200$ from Ref. [14]. The distribution is bimodal and suggests that 10% is the boundary between crystal (>10%) and liquid. (b) The same distribution for $\delta = 1\%$ and 20% for $T^* = 0.359$ and $T^* = 0.306$, respectively. For each $\delta$, two system sizes, $N = 1000$ and 10000 are shown. The KA$_1$ model shows no sign of crystallisation.</figcaption></figure>

In this strategy, the bonds formed by neighboring particles are quantified according to the number of shared neighbors, and therefore, the local crystalline structure provides a specific signature [70]. We perform CNA analysis for the inherent structures of the KA model generated at $T = 0.373$ using parallel tempering (PT) simulation scheme (trajectories taken from Ref. [14]). To define neighbors, we consider a cutoff of $r_{cut}^{AA} = 1.4\sigma_{AA}$, which is the minimum of the first coordination shell of the pair correlation $g_{AA}(r)$. Since the majority of the population is of type A, this is a reasonable choice. At low temperatures, the standard KA model is known to form local FCC crystals, which is well characterized as “142” in the CNA analysis.

In Fig. 9(a), we report the fraction of FCC population as a function of PT steps and their distribution for $T = 0.373$, which provides an upper bound for the fraction %CNA-142 in the liquid near 0.10.

In Fig. 9(b), we report the probability distribution of the fraction of FCC particles at the lowest simulated temperatures for $\delta = 1\%$ and 20% in the KA$_1$ model. The distribution remains well below the 0.10 cutoff, suggesting
ing that 1% of impurity is enough to make the system
very robust against crystallization. Moreover, we con-
firm that a larger system size, N = 10000, does not
show the trend of crystallization. To strengthen this
point, we emphasize that 30 independent swap MC runs
have been performed for the 1% model for a duration

of about 15 $\tau_{\alpha}^{(S)}$ at the lowest temperature ($T = 0.355$) and $N = 1000$ particles. For $N = 10000$ particle case for the 1% model, we consider 10 independent swap MC runs of duration $\approx 1.5 \tau_{\alpha}^{(S)}$ at the lowest temperature ($T = 0.355$). None of them show any signatures of crystallization or demixing.