Beyond packing of hard spheres: The effects of core softness, non-additivity, intermediate-range repulsion, and many-body interactions on the glass-forming ability of bulk metallic glasses

Introduction

When metallic liquids are cooled at rates $R$ exceeding the critical cooling rate $R_{c}$, crystallization can be bypassed and amorphous alloys are formed Telford 2004. Pure metals and most alloys are extremely poor glass formers with $R_{c}>10^{10}$ K/s. In contrast, a number of bulk metallic glasses (BMGs) have been identified with $R_{c}<1$ K/s and critical casting thickness $d_{c}>1$ cm, which enables them to be employed in commercial applications Schroers 2013; Zhong et al. 2014. The discovery of novel BMGs with optimized casting thickness and mechanical properties has largely been a trial-and-error process Greer 1995; Suryanarayana and Inoue 2011. Although combinatorial deposition and characterization techniques Tsai and Flores 2014; Ding et al. 2014 now allow efficient exploration of parameter space, there are an exponentially large number of possible BMG-forming atomic compositions Zhang et al. 2015a. Thus, a quantitative and predictive understanding of the GFA of BMG-forming alloys is necessary to narrow down the vast parameter space.

Silica and polymers possess critical cooling rates that are more than $15$ and $10$ orders of magnitude lower, respectively, than those for pure metals (Fig. 1). Network bonding in silica and chain entanglement in polymers provide the physical mechanisms to inhibit crystallization Zachariasen 1932; Johnson 1999; Hoy and O’Hern 2012. In contrast, the main source of geometric frustration in alloys is the mismatch between atomic sizes Egami and Waseda 1984; Miracle et al. 2003; Miracle 2004; Miracle 2013; Jalali and Li 2004; Jalali and Li 2005. Molecular dynamics simulations of binary hard spheres have shown that tuning the atomic size ratio can decrease $R_{c}$ by more than $13$ orders of magnitude Zhang et al. 2014. Packing of hard spheres can also rationalize the correlation between the number of components, their atomic size ratios, and the GFA of BMGs Zhang et al. 2015a.

Although the packing of hard spheres plays an important role in determining the GFA of alloys, it is obvious that metals possess additional features that are not represented by hard-sphere interactions. Other features of metallic interactions, such as metallic bonding Zhang et al. 2015b, the form of the interatomic pair potential, and many-body interactions Daw and Baskes 1983, can change the crystalline structure that competes with glass formation and change the prediction of $R_{c}$ by several orders of magnitude from the hard-sphere value. Compared to the $\sim 13$ orders of magnitude variation in $R_{c}$ that results from the packing of hard-spheres, changes to $R_{c}$ are small, but not negligible and may explain the crucial differences between an amorphous film and a bulk metallic glass. Since the critical casting thickness $d_{c}$ is negatively correlated with $R_{c}$ and increasing $R_{c}$ by two orders of magnitude can reduce $d_{c}$ by one order of magnitude Inoue 2000, more accurate models of intermetallic potentials are needed to identify BMGs with $d_{c}>1~{\rm cm}$ (Fig. 1).

The interatomic potential in the embedded atom model (EAM) is frequently implemented in computational studies of the structural and mechanical properties, as well as the dynamics, of metallic systems Daw and Baskes 1983. The EAM potential energy includes a pairwise-additive term, which is in general different from the hard-sphere and Lennard-Jones pair potentials (Fig. 2 (a)), and a many-body contribution from the electron charge density, which is fitted to ab initio calculations of lattice parameters, elastic constants, and other thermodynamic properties Mendelev and Ackland 2007; Sheng et al. 2011.

In this manuscript, we seek to identify the key features of the pairwise and many-body interactions that strongly influence the GFA of alloys. For example, we investigate the effects of the softness of the pairwise repulsive core, pairwise non-additivity, and the form of the pairwise intermediate-range repulsion on the GFA. We then measure the GFA for the full embedded atom models of several pure metals and BMGs to determine the contribution of the many-body interactions on the GFA. We find that the changes in the GFA arising from variations in the pair and many-body contributions of the embedded atom model are small compared to the $13$ orders of magnitude change in GFA between monoatomic and binary and ternary hard-sphere systems. However, these peturbations to the GFA may still be important for discovering new bulk metallic glass formers.

The manuscript includes three additional sections after the introduction. First, in Sec. 2, we describe the hard-sphere, repulsive Lennard-Jones, Lennard-Jones, and Dzugutov-Shi potentials used to vary the form and non-additivity of the pairwise interactions. We also introduce the embedded atom model for pure metals and alloys. For each interatomic potential, we discuss the methods employed to measure the critical cooling rate $R_{c}$. We then report the results for the GFA for all interaction potentials in Sec. 3. We conclude the manuscript in Sec. 4.

models and methods

As described above, the embedded atom model for metallic systems includes pairwise and many-body interactions. In this section, we define three metrics (core softness, non-additivity, and intermediate-range repulsion) to characterize the form of the pairwise interactions. We describe molecular dynamics simulations of monodipserse and binary systems interacting via generalized Lennard-Jones or Dzugutov-Shi Dzugutov 1992; Shi et al. 2014 potentials to quantify the effects of the softness of the repulsive core and strength of the intermediate-range repulsion on the GFA. We also introduce molecular dynamics simulations of binary hard spheres to study variations in the GFA from non-additive pairwise interactions. We estimate values for the pairwise core softness, non-additivity, and form of the intermediate-range repulsive interactions from fits to the pairwise contributions of the EAM for pure metals and binary BMGs. We also introduce the Lennard-Jones and full EAM potentials that we employ to study the effects of many-body interactions on the GFA.

To tune the softness of the pairwise repulsive core Zhen and Davies 1983, we employ the generalized $m$-$n$ Lennard-Jones (LJ) potential (Fig. 2 (b)),

where $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, $\sigma_{i}$ is the diameter of atom $i$, and $\epsilon$ is the energy scale of the interaction. The interaction potential has a minimum $u_{m}=-\epsilon$ at $r_{m}=2^{1/6}\sigma_{ij}$. The exponent $m$ (or equivalently the curvature $\kappa$ of the pair potential at the minimum) controls the softness of the repulsive core, where smaller $m$ corresponds to softer interactions. Note that the generalized Lennard-Jones potential is fixed at $u_{\rm LJ}(r_{ij})\equiv u_{\rm LJ}^{12-6}(r_{ij})$ for $r_{ij}>r_{m}$. To separate the effects of the attractive interactions from the repulsive core, we also studied the generalized $m$-$n$ repulsive Lennard-Jones (RLJ) potential Weeks et al. 1971 as shown in Fig. 2 (b):

To obtain physical values for the softness exponent $m$, we fit the repulsive part of the EAM pair potential of typical BMG-forming elements to $u_{\rm RLJ}^{m-6}(r)$. As shown in Table 1, we find that $m$ varies from approximately $3$ to $14$. The repulsive cores for most metals appear softer than Lennard-Jones interactions with $m=12$.

To investigate the effects of softness of the pairwise repulsive core on the GFA of metallic systems, we performed molecular dynamics (MD) simulations of $N=1372$ spherical atoms with mass $m_{0}$ that interact via the generalized Lennard-Jones and repulsive Lennard-Jones potentials with $n=6$ and a range of $m$ values. We studied three binary LJ systems with softness exponents $m_{A}=m_{B}=m_{AB}=12$ (LJ12-6), $m_{A}=m_{B}=m_{AB}=5$ (LJ5-6), and $m_{A}=12$, $m_{B}=5$, and $m_{AB}=8$ (LJ12-6/LJ5-6). We set the atomic diameter ratio to be $\alpha=\sigma_{B}/\sigma_{A}=0.95$ and varied the number fraction of small atoms $x_{B}=N_{B}/N$ from $0$ to $1$. Temperatures and times are given in units of $\epsilon/k_{B}$ and $\sigma_{A}\sqrt{m_{0}/\epsilon}$, respectively. After equilibrating the systems at high temperature $T_{0}=2$, the liquids were cooled exponentially $T(t)=T_{0}\exp(-Rt)$ with rate $R$ to low temperature, $T_{f}=0.01$, using the Gaussian constraint thermostat Allen and Tildesley 1987 with time step $\Delta t=0.001$. Constant volume $V$ simulations at number density $\rho\sigma_{A}^{3}=N\sigma_{A}^{3}/V=1$ were performed for both the LJ and RLJ models. For LJ systems, we also cooled systems with the constraint that the pressure $p$ (in units of $\epsilon/\sigma_{A}^{3}$) decreased exponentially in time from an initial pressure $p_{0}=1$ to final pressure $p_{f}=0.001$ according to

using a Gaussian constraint barostat Allen and Tildesley 1987. A cooling rate of $R=1$ in the units used in the MD simulations corresponds to a cooling rate of $10^{15}$ K/s using $\sigma_{A}\sim 3\times 10^{-10}$ m, $\epsilon/k_{B}\sim 10^{3}$ K, and molar mass $M\sim 10^{-1}$ kg/mol, which are typical values for BMGs Zhen and Davies 1983.

2 Non-additive binary hard spheres

The sizes of metallic atoms are often estimated from the first peak of the radial distribution function $g(r)$ of crystalline and disordered solids Hume-Rothery 1950. In binary alloys with species $A$ and $B$, the repulsive core $\sigma_{AB}$ between atoms $A$ and $B$ can differ from the average diameter $\overline{\sigma}_{AB}=(\sigma_{A}+\sigma_{B})/2$. We quantify the non-additivity of the pairwise repulsive core using the parameter

Many binary alloys possess $\Sigma<0$, which indicates that the repulsive core $\sigma_{AB}$ between $A$ and $B$ atoms is smaller than the average diameter. We list $\sigma_{A}$, $\sigma_{B}$, $\sigma_{AB}$, and $\Sigma$ for several binary alloys obtained from EAM calculations of $g(r)$ in Table 2. Non-additive binary hard spheres have been shown to form exotic crystalline structures, in particular intermetallic compounds Punnathanam and Monson 2006; Woodcock 2011. In addition, non-additivity due to bond shortening with $\Sigma<0$ can lead to unusual intermediate-range order in BMGs Cheng et al. 2009; Senkov et al. 2012; Sheng et al. 2006. The well-studied Kob-Andersen model for $\rm{Ni}_{80}{\rm P}_{20}$ glasses also has $\Sigma=-0.149$ Kob and Andersen 1995.

To study the effects of nonadditivity on the GFA, we compressed $N=500$ binary hard spheres with mass $m_{0}$ that interact pairwise via

over a range of diameter ratios $\alpha$ and number fractions of the small sphere $x_{B}$ using event-driven MD simulations. We first equilibrated liquid states at packing fraction $\phi=0.25$. To compress the system, we ran the MD simulations at constant volume for a time interval $\tau$, and then compressed the system instantaneously until the closest pair of spheres came into contact Jalali and Li 2004; Zhang et al. 2014. We performed successive compressions until the pressure increased to $10^{3}$, which corresponds to $(\phi_{J}-\phi)/\phi_{J}<10^{-3}$, where $\phi_{J}$ is the packing fraction at the onset of jamming. We varied the compression rate $R\equiv 1/\tau$ over $5$ orders of magnitude Zhang et al. 2014. We report $R$ in units of $\sqrt{k_{B}T/m_{0}\sigma_{A}^{2}}$. Note that in these units $R=1$ corresponds to a cooling rate of $10^{12}$ K/s for alloys Truskett et al. 2000.

3 Dzugutov-Shi (DZ) potential

The pair potential of many metallic systems includes intermediate-range repulsive interactions Friedel 1958 in addition to short-range attractive interactions, which can give rise to intermediate-range positional order Fujima et al. 2007; Wu et al. 2013. Intermediate-range pairwise repulsive interactions are often modeled using the Dzugutov potential Dzugutov 1992; Dzugutov 1993; Roth and Denton 2000; Doye and Wales 2001. Shi et. al. introduced a modified version of the original Dzugutov potential that allows one to continuously tune the interaction potential between the LJ potential to one that includes intermediate-range repulsion Shi et al. 2014. The Dzugutov-Shi (DZ) potential is given by

where the “bump” potential $u_{\rm bump}(r_{ij})$ models the intermediate-range repulsive interactions using a sinusoidal pulse,

of the strength $\xi$ within the range $\lambda\sigma_{ij}\leq r_{ij}\leq\delta\sigma_{ij}$. The location of the peak and width of $u_{\rm bump}$ are given by $(\lambda+\delta)/2$ and $\delta-\lambda$. To obtain physical values for $\xi$, $\lambda$, and $\delta$, we fit the DZ potential to the EAM pair potential for several elements. We show values of $\xi$, $\lambda$, and $\delta$ for elements commonly found in BMGs in Table 3. Pb, Pd, Pt, Mg, Fe, Ta, Au, Ti, Mo, W, and Nb do not have significant intermediate-range repulsive interactions.

To study the effects of intermediate-range repulsive interactions on the GFA, we performed MD simulations of $N=1372$ spherical atoms that interact pairwise via the DZ potential. We followed the same cooling protocol as used for the simulations of Lennard-Jones systems with pressure that decreases exponentially in time as discussed in Sec. 2.1. We fixed the strength of the intermediate-range repulsive interactions at $\xi=0.35\epsilon$ and varied $\lambda$ and $\delta$ to tune the location of the peak $(\lambda+\delta)/2$ and range $\delta-\lambda$ of $u_{\rm bump}$. We also studied binary mixtures composed of $A$ atoms that interact via the DZ potential with $\xi=0.35\epsilon$, $\lambda=1.2$, and $\delta=2.15$, and $B$ atoms that interact via the LJ potential with diameter ratio $\alpha=0.95$. The number fraction of small atoms $x_{B}$ is varied from $0$ to $1$ in steps of $0.2$.

4 LJ-EAM and EAM potential

The total potential energy $U$ employed in the embedded-atom model for metals includes pairwise and many-body contributions:

where the many-body embedding function $F_{i}$ depends on the electron density associated with each atom $i$ (normalized by $e/\sigma_{A}^{3}$) and $\overline{\rho}^{e}_{i}=\sum\limits_{j\neq i}\rho^{e}(r_{ij})$ Daw and Baskes 1983; Mendelev and Ackland 2007; Sheng et al. 2011. To quantify the effects of the many-body interactions on the GFA, we focused on the LJ-EAM potential, where $u(r_{ij})=u_{LJ}(r_{ij})$, $F_{i}(\overline{\rho}^{e}_{i})=A\overline{\rho}^{e}_{i}(\ln\overline{\rho}^{e}_{i}-r_{m}/\sigma_{A})/2$ and $\rho^{e}(r_{ij})=C\exp[-\beta(r_{ij}-r_{m})]$, where $C$ and $r_{m}$ are calibrated to experimental data on alloys Baskes 1999; Nam et al. 2007. We set the atomic diameter $\sigma_{A}=2.8$ ${\AA}$ and attraction depth $\epsilon=0.2$ eV for the LJ potential to match the pair potential of typical metals such as Zr. The parameters $A$ and $\beta$ control the many-body interaction strength and inverse decay length of the electron density, respectively.

We performed MD simulations of the LJ-EAM for several pure metals and of the full EAM for several binary alloys using the LAMMPS simulation software Plimpton 1995. We cooled systems in the liquid state to low temperature at constant zero pressure at different rates $R$. The initial and final temperatures for several systems (specified by $A$ and $\beta$) are summarized in Table 4. For our studies of the full EAM potential, we set $N=4000$ and fixed the initial and final temperatures at $T_{i}=2000K$ and $T_{f}=300K$.

5 Critical cooling rate

To calculate the critical cooling rate $R_{c}$ for each metallic system, we initialized the liquid state at high temperature, cooled the system exponentially to low temperature at a given rate $R$ at either fixed volume or exponentially decaying pressure as in Eq. 3, and measured the global bond orientational order parameter $Q_{6}$ Steinhardt et al. 1983. For hard-sphere interactions, we compressed the systems so that the packing fraction approached that at jamming onset exponentially, which is thermodynamically equivalent to cooling systems exponentially Parisi and Zamponi 2010. For all systems studied, the average global bond orientational order parameter $Q_{6}$ versus $\log R$ possesses a sigmoidal shape with a midpoint defined by $R_{c}$. Below, we show results for $R_{c}$ for the pair potentials described in Secs. 2.1-2.3 and the full and LJ-EAM potential in Sec. 2.4.

Results

To investigate the effects of softness of the repulsive core on the GFA, we first measured the critical cooling rate $R_{c}$ for monodisperse systems that interact via the generalized LJ (Eq. 1) and RLJ (Eq. 2) pairwise potentials as a function of the softness exponent for $m=1$, $3$, $5$, $8$, $10$, and $12$. As shown in Fig. 3, when cooling at constant number density $\rho\sigma_{A}^{3}=1$, the GFA increases weakly ($R_{c}$ decreases by less than an order of magnitude) as the repulsive core becomes softer ($m$ decreases). When cooling a LJ system with a pressure that decays exponentially in time as in Eq. 3, the dependence of $R_{c}$ on the softness exponent $m$ is even weaker, except for systems with extremely soft core repulsions with $m=1$. In contrast, most atomic species that are found in BMGs possess $m>4$ (Table 1).

As shown in Fig. 3, the crystalline structures that compete with glass formation in systems with core-softened RLJ interactions at $\rho\sigma_{A}^{3}=1$ are face-centered cubic (FCC) for all exponents $m$ studied. In addition, FCC crystals compete with glass formation in LJ systems, but as the repulsive core softens, body-centered cubic (BCC) crystals become more stable Hoover et al. 1972. We find that BCC is the crystal type that competes with glass formation for $m=3$ LJ systems cooled at constant density $\rho\sigma_{A}^{3}=1$ and for $m=3$ and $5$ LJ systems cooled such that the pressure obeys Eq. 3.

Structural characterizations of atomic systems that interact via the generalized LJ potential are shown in Fig. 4 for cooling rates $R>R_{c}$. As the repulsive core of the potential becomes softer (i.e. $m$ decreases), the attractive well of the potential widens to include second-neighbor attractive interactions, which can compensate repulsive first-neighbor interactions. Indeed, LJ systems with $m=1$ and $3$ exhibit phase separation into dilute and compressed regions when cooled at fixed density $\rho\sigma_{A}^{3}=1$ and volume contraction, where the first neighbor separations are smaller than the location of the potential minimum, when cooled such that the pressure obeys Eq. 3. In fact, the $m=3$ LJ system displays two isostructural glassy states, contracted and expanded, with different densities as shown in the inset to Fig. 4. Similar isostructural transitions have been found in equilibrium systems with narrow-ranged attractive interactions Frenkel 2006. Large density differences between polymorphs in metallic glasses such as those found in Ce55Al45 are often attributed to electronic many-body interactions Sheng et al. 2007. However, here we show that softening the pairwise repulsive core (which increases the range of the attractive well) can also give rise to polymorphs with different densities.

We also investigated the effects of core softness on the glass-forming ability in binary mixtures that interact via the generalized $m$-$6$ LJ potential. We focused on three mixtures with diameter ratio $\alpha=\sigma_{B}/\sigma_{A}=0.95$: (1) conventional LJ systems with $m=12$, (2) core softened LJ systems with $m=5$, and (3) mixtures of LJ systems with $m=12$ ($A$ species) and $m=5$ ($B$ species). While FCC is the crystalline structure that competes with glass formation for binary LJ systems with $m=12$, BCC is the competing crystalline structure for binary mixtures with $m=5$ for all number fractions $x_{B}$ as shown in Fig. 5. For both $m=12$ and $m=5$ systems, the variation in $R_{c}(x_{B})$, which is less than an order of magnitude, is controlled by the diameter ratio $\alpha=0.95$. In binary mixtures of LJ systems with $m=12$ and $m=5$ interactions, FCC remains the crystalline structure that competes with glass formation, except when $x_{B}\approx 1$. However, because of the incompatibility between FCC and BCC crystalline structures, the GFA for the $m=12$ and $m=5$ LJ mixtures is significantly enhanced compared to glasses with $m=12$ or $m=5$ interactions alone. For example, Ni-Ta is a good glass former despite the fact that it possesses a diameter ratio near unity ($\alpha\approx 0.9$) Wang et al. 2010. Incompatibility between competing BCC and FCC crystal structures is a possible cause of the enhanced GFA. As shown in Table 1, Ni has a relatively large pairwise repulsive exponent ($6<m<10$) with equilibrium FCC structure, while Ta has a relatively small exponent ($3<m<5$) with equilibrium BCC structure Hume-Rothery 1950. Since the softness exponents of the pairwise interactions vary significantly from one element to another (Table 1), softness-induced competing crystal incompatibility can enhance the GFA of binary and multi-component BMG-forming alloys.

2 Non-additivity

We performed event-driven molecular dynamics simulations of binary non-additive hard spheres (Sec. 2.2) to investigate the effects of non-additivity of the pairwise repulsive interactions on the GFA of alloys. We measured the critical cooling rate $R_{c}$ of non-additive binary hard spheres with diameter ratios $\alpha=\sigma_{B}/\sigma_{A}=1.0$, $0.97$, $0.95$, $0.93$, $0.9$, and $0.5$ and number fractions of the small spheres $x_{B}=0.5$ and $2/3$ over a range of non-additivity parameters $\Sigma$. Since $\Sigma>0$ is rare among binary alloys (Table 2), we expect that hard-sphere systems with positive non-additivity are poor glass-formers. For example, we find that systems with $\alpha=1$ and $\Sigma=0.05$ display strong demixing between $A$ and $B$ particles and are not good glass formers.

Our previous studies of additive binary hard spheres ($\Sigma=0$) have shown that well-mixed FCC solid solutions are the crystal structures that compete with glass formation when $\alpha\gtrsim 0.8$, while the systems tend to demix when $\alpha\lesssim 0.8$ Zhang et al. 2014. For $\Sigma<0$ and $\alpha=1.0$, $0.97$, $0.95$, $0.93$, and $0.9$, the GFA improves as $\Sigma$ becomes more negative, and the competing crystal structure remains the FCC solid solution (Fig. 6). The change in $R_{c}$ with decreasing $\Sigma$ also increases as $\alpha$ decreases with roughly an order of magnitude difference in $R_{c}$ between systems with $\Sigma=0$ and $\Sigma=-0.05$ at $\alpha=0.9$. Enhancement of the GFA arising from non-additivity of the repulsive cores ($\Sigma<0$) has also been observed in LJ systems Zhang et al. 2013.

For binary systems with large atomic size differences (i.e. $\alpha\ll 0.8$), the variation of $R_{c}$ with $\Sigma$ is opposite to that obtained for binary systems with small atomic size differences. As shown in Fig. 6, we find that $R_{c}$ grows with increasing $\Sigma$ at $\alpha=0.5$. For $\alpha=0.5$ and $\Sigma<0$, compound crystals are the ordered structures that compete with glass formation since negative non-additivity promotes mixing. As an example, although the $AB_{2}$ compound is the densest crystal for binary hard spheres with $\alpha=0.5$ and $\Sigma=0$, it is not kinetically accessible during compression due to the strong drive for demixing Filion and Dijkstra 2009; Hopkins et al. 2011; Zhang et al. 2014. However, when $\Sigma$ becomes negative (e.g. $\Sigma=-0.05$), we find that the $AB_{2}$ compound forms easily for the compression rates that we studied, as shown in the inset to Fig. 6. Thus, the formation of intermetallic compounds in alloys can be enhanced by pairwise negative non-additivity among different atomic species.

3 Intermediate-range Repulsive Interactions

We also investigated crystallization and glass formation as a function of the form of intermediate-range repulsive pairwise interactions (Sec. 2.3). We first performed molecular dynamics simulations of monodisperse spheres interacting via the DZ potential (Eq. 6) at fixed strength $\xi=0.35\epsilon$ and varying peak location $(\lambda+\delta)/2$ and width $\delta-\lambda$. In Fig. 7, we plot the critical cooling rate $R_{c}$ as a contour plot versus $(\lambda+\delta)/2$ and $\delta-\lambda$ over ranges that are relevant to BMGs (Table 3). We find several regions of good glass-forming ability (small $R_{c}$) and different crystal structures that compete with glass formation. For a large region of parameter space, FCC is the competing crystal structure. BCC is the competing crystal structure when the location of the peak in $u_{\rm bump}$ approaches third-neighbor separations at $r_{ij}\approx\sqrt{3}r_{m}$. We also find an “8-4” crystal structure that competes with glass formation, with atom positions located on embedded octagons and squares when they are projected into two dimensions. (See the inset of Fig. 7). In three dimensions, one can see that the atoms forming the octagons and squares are located in alternating stacked layers. (See Fig. 8 for a comparison of the radial distribution functions for FCC, BCC, and $8$-$4$ crystals.) When the intermediate-range repulsion becomes too strong (i.e. large $\delta$), microphase separation becomes energetically favorable compared to macroscale phase separation Brazovskii 1975; Seul and Andelman 1995.

We also studied the critical cooling rate $R_{c}$ for binary mixtures (e.g. Zr-Cu alloys), in which one component possesses intermediate-range repulsive interactions and the other component does not. We focused on binary systems with atoms that interact via the DZ ($A$ species) and LJ potential ($B$ species) with diameter ratio $\sigma_{B}/\sigma_{A}=0.95$. For the DZ potential, we set the parameters $\xi/\epsilon\approx 0.4$, $\lambda\approx 1.2$, and $\delta\approx 2.2$ to mimic those of Zr atoms (Table 3). As shown in Fig. 5, $R_{c}$ for this binary mixture is suppressed by more than two orders of magnitude compared to the pure system with LJ or DZ interactions alone because the two species possess incompatible equilibrium crystal structures (i.e. FCC and BCC). This mechanism of incompatible equilibrium crystal structures may explain the exceptionally good glass-forming ability of the Zr-Cu system, even though it is a binary, rather than, multi-component alloy.

4 LJ EAM for Monoatomic Systems

To determine the relative contributions of the pairwise and many-body interactions to the GFA of alloys, we performed molecular dynamics simulations of the LJ-EAM potential (Sec. 2.4) as a function of the many-body interaction strength $A$ and electron density inverse decay length $\beta$ for monoatomic systems. In Fig. 9, we show the critical cooling rate $R_{c}$ for monodisperse LJ-EAM systems as a function of $A$ for $\beta=2$, $4$, and $6$ ${\AA}^{-1}$. We find that $R_{c}\approx 10^{13}~{\rm K/s}$. $R_{c}$ changes by less than one order of magnitude as $A$ and $\beta$ are varied over the range that is relevant for elements found in BMGs even though the total potential energy per atom $U/N$ varies linearly with $A$. We also find that FCC crystals are the ordered structures that compete with glass formation in monoatomic LJ-EAM systems over the full parameter range for $A$ and $\beta$. Thus, we argue that many-body interactions have a weak influence on the GFA compared to the pairwise interactions for monoatomic systems.

5 Full EAM for Binary Alloys

We also measured the critical cooling rate $R_{c}$ for several binary alloys as a function of the number fraction $x_{B}$ of the small atomic species using the full EAM potential. We focused on Zr-Cu, Mg-Al, and Cu-Ni alloys with atomic diameter ratios that range from $\alpha=0.79$ to $0.98$. In Fig. 10, we compare $R_{c}$ versus $x_{B}$ from simulations of the full EAM potential for these alloys to $R_{c}$ obtained from simulations of additive hard spheres with comparable values of $\alpha$ Zhang et al. 2014.

As expected, $R_{c}$ for binary alloys with $\alpha\sim 1$ (i.e. Cu-Ni) is nearly independent of $x_{B}$. In addition, when the hard-sphere simulations with $\alpha=1$ are calibrated to Ni, $R_{c}$ from simulations of the hard-sphere and EAM potentials agree semi-quantitatively. From our previous simulations of hard spheres Zhang et al. 2014, we know that $R_{c}(x_{B})$ develops a deep minimum that shifts to larger $x_{B}$ as $\alpha$ decreases from unity. For example, when $\alpha=0.9$, $R_{c}$ for hard-sphere systems at $x_{B}\approx 0.6$ is two orders of magnitude less than the value when $\alpha=1$. Although we are not able to simulate sufficiently slow rates, it appears that $R_{c}$ at the minimum in $x_{B}$ for Mg-Al with $\alpha=0.94$ will decrease by at least two orders of magnitude and the minimum in $R_{c}(x_{B})$ will occur at $x_{B}>0.5$. We also find similar results for $R_{c}$ for hard spheres with $\alpha=0.79$ and for EAM of Zr-Cu with a deep minimum in the range $0.2<x_{B}<0.8$.

We also determined the crystal structures that compete with glass formation in the full EAM simulations of binary alloys. We find that FCC (or HCP) is most often the competing crystal structure, as in simulations of additive binary hard spheres, but we also find exceptions. In particular, we show that on the Zr-rich side of Zr-Cu, BCC crystal structures compete with glass formation. The BCC equilibrium structure for the Zr-Cu alloys can likely be attributed to the pairwise part of the EAM potential. For example, the pair potential for Zr possesses intermediate-range repulsive interactions with the location of peak $(\lambda+\delta)/2=1.70$ and width $\delta-\lambda=1.08$ (Table 3) in a region of parameter space that has been shown to display BCC crystal structure (Fig. 7).

Conclusion

The hard-sphere model has provided a predictive description of crystallization and glass formation in simple liquids Weeks et al. 1971. In addition, we have shown in recent studies that the additive hard-sphere model can explain more than $13$ orders of magnitude variation in the critical cooling rate $R_{c}$, which nearly spans the full range of GFA from that for pure metals to that for the best BMGs Zhang et al. 2014. We also showed that the best binary and ternary BMGs occur in the region of parameter space (i.e. diameter ratio and number fraction) with the smallest values of $R_{c}$ for hard spheres.

However, in metallic systems, there are a number of additional features of the interatomic potential beyond hard-core repulsions, including softness, non-additivity, and range of the pairwise interactions. For example, metallic atoms typically appear softer (with smaller values of the exponent of the repulsive core) than the commonly used LJ pair potential and possess several per cent negative non-additivity due to shortening of metallic bonds Cheng et al. 2009. In addition, Friedel oscillations in metals give rise to intermediate-range repulsion at separations beyond the short-range attractive well Friedel 1958. The interatomic potential for metals also includes many-body interactions from the electronic degrees of freedom. In this manuscript, we investigated how these additional features affect the GFA of pure and binary metallic systems.

We performed molecular dynamics simulations of several model systems to study the effects on the GFA for each of the key features of the interatomic potential separately. For example, we performed simulations of monodisperse and binary spheres that interact via the generalized LJ and DZ pair potentials to quantify the effect of the softness of the repulsive core and form of the intermediate-range repulsive interactions on the GFA. We also performed MD simulations of non-additive binary hard spheres to quantify the effects of non-additivity on the GFA. We found that softness, non-additivity, and form of the intermediate-range repulsions cause deviations in $R_{c}$ that are only $1\sim 2$ orders of magnitude from the additive hard-sphere predictions.

While FCC is the most stable crystal structure for LJ and hard-sphere systems, softening of the repulsive core gives rise to novel contracted disordered structures, as well as the formation of BCC crystals. We also showed that negative non-additivity of the repulsive core in binary alloys improves the GFA when the competing crystal structures are solid solutions. However, when the atomic size ratio is in the demixing regime ($\alpha<0.8$), negative non-additivity can favor the formation of compound crystals and decrease the GFA. The crystal structure that competes with glass formation, and thus the GFA, also depends sensitively on the form of the intermediate-range repulsive interactions. We find that when the competing crystal structures of each component in an alloy are incompatible (e.g. FCC and BCC), the GFA can be enhanced compared to hard-sphere predictions.

We also investigated the relative contributions of the pairwise and many-body interactions to the GFA by performing molecular dynamics simulations of the LJ-EAM potential. We found that including the many-body interactions only changes $R_{c}$ by less than one order of magnitude compared to that when the many-body interactions are not included. We also calculated $R_{c}$ for several binary alloys using the full EAM potential and found qualitatively the same results as for binary hard spheres. Thus, we argue that hard-sphere interactions provide a qualitatively accurate model for predicting the GFA of alloys. Other features of the interatomic potential (beyond additive hard-core repulsion) give rise to only $1$-$2$ orders of magnitude variation of $R_{c}$, which is small compared to the more than $13$ orders of magnitude variation predicted by hard-sphere systems. Despite this, including additional features to the interatomic potential beyond hard-sphere interactions is important for the design of new BMGs since precise quantification of the critical casting thickness can determine whether a new BMG is commercially viable.

References