Atomistic Simulations to Predict Favored Glass-Formation Composition and Ion-Beam-Mixing of Nano-Multiple-Metal-Layers to Produce Ternary Amorphous Films

Based on the framework of long-range empirical formulas, the interatomic potentials were constructed for the Ni-Nb-Mo (fcc-bcc-bcc) and Ni-Zr-Mo (fcc-hcp-bcc) ternary metal systems. Applying the constructed potentials, atomistic simulations were performed to predict the energetically favored glass formation regions (GFRs) in the respective composition triangles of the systems. In addition, the amorphization driving forces (ADFs), i.e., the energy differences between the solid solutions and disordered phases, were computed and appeared to correlate with the so-called glass forming abilities. To verify the atomistic prediction, ion beam mixing with nano-multiple-metal-layers was carried out to produce ternary amorphous films. The results showed that the composition of ternary amorphous films obtained by ion beam mixing all locate inside the GFRs, supporting the predictions of atomistic simulations. Interestingly, the minimum ion dosage required for amorphization showed a negative correlation with the calculated ADF, implying that the predicted amorphization driving force could be an indicator of the glass formation ability.


Introduction
Amorphous alloys (i.e., metallic glasses, MGs) have attracted significant interest due to their eminent physical, chemistry and mechanical properties [1][2][3][4][5]. Several powerful non-equilibrium processing techniques, e.g., liquid melt quenching [6], ion beam mixing (IBM) [7,8], mechanical alloying [9,10] and laser-melting techniques [11][12][13], have been developed to produce amorphous alloys. Among these techniques, the ion beam mixing [14] technique is capable of fabricating MGs in both miscible and immiscible metal systems [15,16]. Since the effective cooling speed is as high as 10 12 -10 13 K/s, IBM is a very powerful and effective method to produce a number of amorphous alloys in equilibrium immiscible systems. In the field of MGs, one fundamental scientific issue is estimating glass formation ability (GFA), which describes the level of difficulty or easiness of metallic glass formation [17][18][19]. In the past decades, researchers have proposed several empirical criteria or rules to indicate the GFA of an alloy. For instance, with regard to the equilibrium phase diagram, Turnbull et al. [20] have suggested a deep eutectic criterion to predict metallic glass formation upon liquid melt quenching. By employing extensive IBM studies, Liu et al. [7] have put forward a structural difference rule to predict amorphous alloy formation by IBM, and further proposed the total width of two-phase regions derived from the equilibrium phase diagram could characterize the GFR of a binary

Construction of Ternary Interatomic Potential
The interatomic potential of a ternary system describes the interatomic interactions among the atoms involved. Therefore, if the potential is determined, most of the physical properties, including the GFA as well as the atomic configuration, can be computed or derived from relevant simulations. Great efforts have been devoted to develop realistic potentials, e.g., the embedded atom method (EAM) [45,46], the Finnis-Sinclair potential [47], the second-moment approximation of tight-binding (TB-SMA) potential [48],and their various modifications [49][50][51].
Dai et al. [52][53][54] have proposed a long-range empirical potential for bcc and fcc metals, and also extended this to hcp metals. Based on the long-range empirical potential [55,56], the potential energy E i of atom i can be calculated as follows: where r ij is the distance between atoms i and j of the alloy system at equilibrium. The pair term V(r ij ) and electron-density term ϕ(r ij ) can be expressed, respectively, by: where r c1 and r c2 are the cutoff radii for the pair and electron-density terms, respectively, α and c i are the potential parameters to be determined by the fitting procedure. The exponents m and n are integers that can be adjusted according to a specific system. For an A-B-C ternary metal system constituted by any combination of the fcc, hcp and bcc metals, there should be six sets of potential parameters, i.e., three sets for interactions of pure metals A-A, B-B and C-C, and three sets for the cross interactions of A-B, B-C and C-A. The potential parameters of pure metals can be determined by fitting to the basic physical properties, i.e., the cohesive energies, lattice constants, bulk moduli, and elastic constants of the pure metals [57]. The cross potential parameters of A-B, B-C and C-A can be determined by fitting to the basic physical properties of the intermetallic compounds with various compositions and structures. To obtain the physical properties of the compounds, the first-principle calculations were performed by employing the Cambridge Serial Total Energy Package (CASTEP) [58,59] in Material Studio, which uses the projector augmented wave (PAW) method [60] to yield higher computation efficiency. During the calculation, the exchange and correlation functions are chosen as the generalized-gradient approximation (GGA) of Perdew and Wang (PW91) [61]. Firstly, the geometry optimization is performed to calculate the lattice constants and total energies of compounds at equilibrium, and then the elastic constants and bulk moduli are calculated. Besides, the cohesive energy of the compounds is also derived from the total energy obtained by first-principle calculation.
The fitted potential parameters of the Ni-Nb-Mo (fcc-bcc-bcc) system [55] are presented in Table 1. Furthermore, the potential parameters of the Ni-Zr-Mo (fcc-hcp-bcc) system [56] constituted by three major different crystalline structures are presented in Table 2. Tables 3 and 4 present the lattice constants, cohesive energies, and bulk moduli of B2 and L1 2 compounds in the Ni-Nb-Mo and Ni-Zr-Mo systems computed from the constructed interatomic potentials. The physical properties derived from potentials match well with those via the first-principle calculations and experiments [62][63][64]. This proves the reliability of the constructed interatomic potentials. In addition, to examine whether constructed potentials can describe interactions at non-equilibrium states or not, the equation of states (EOS) derived from potentials was computed and compared with the Rose Equation [65]. From Figures 1 and 2, it can be seen that the energy curves calculated from potentials remain smooth and continuous with the Rose Equation over the whole range, for both the Ni-Nb-Mo and Ni-Zr-Mo systems. The constructed long-range empirical potentials turn out to effectively determine the interatomic       Table 3. Lattice constants (a, Å), cohesive energies (E c , eV), and bulk moduli (B 0 , Mbar) of B2 and L1 2 compounds in the Ni-Nb-Mo system derived from the potential (first line) and first-principle calculation (second line) [55].

Atomistic Simulation
Based on the constructed ternary interatomic potential, a series of molecular dynamics (MD) simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [66]. To study the relative stabilities of the Ni-Zr-Mo (fcc-hcp-bcc) ternary system, three types of solid solution models, i.e., fcc, hcp, and bcc models, were constructed according to the dominate component of the alloy composition. For example, the fcc solid solution model consists of 6912 (12 × 12 × 12 × 4) atoms, among which the main component is Ni, then the solute atoms were randomly substituted by a certain number of solvent atoms to achieve a desired concentration. Periodic boundary conditions were applied in the three directions. With a time-step of 5 × 10 −15 s, MD simulations were performed in the framework of an isothermal-isobaric ensemble. To reach a stable state, the simulation was performed at 300 K and 0 Pa for millions of time steps until all the related dynamic variables exhibited no secular variation. In addition to MD simulations, a series of MC simulations were performed to calculate the formation energy of the solid solution [25,26]. The MC simulations were conducted at 300 K and 0 Pa under an isothermal-isobaric ensemble. They were performed in the box deformation model, where the fractional coordinates of atoms in the box are fixed but the box could either expand or shrink. To reach the minimum energy state, the solid solution model is adjusted by optimizing the lattice constants with the same symmetry. According to recent studies [67], the simulated solid solution model can be compared to the liquid melt quenching process via the atomistic simulations. The inherent hierarchical structure and its effect on the mechanical property of MGs are clarified by the solid solution model and liquid melt quenching method via atomistic simulations. It was revealed that both producing techniques exhibit no pronounced differences in the local atomic structure and mechanical behavior. Therefore, the atomistic simulations in the present work reveal the physics of the thin film metallic glasses.
The structural changes for the solid solution models are monitored by the structure factor S(q). The S(q) [68] can be computed by where N is the number of atoms, q is the scattering vector, b k and r k are the scattering length and position vector of atom k, respectively. Besides, the S(q) calculated by Equation (4) is not normalized.

Ion Beam Mixing
To design the nano-multiple-metal-layers, we need to consider and determine three parameters, i.e., the total thickness of the layers, the relative thickness of the two constituent metals and the individual layer thickness in the samples. During the ion beam mixing experiment, uniform mixing was achieved when the total thickness of the layers was measured to be in the projected range (R p ), plus the projected range straggling (∆R p ) of the irradiation ions [69]. According to the transport of ions in matter (TRIM) program [70], the total thickness of the A-B-C ternary metal layers was computed to be around 40 nm, matching the irradiation ion range. The relative thickness of two constituent metals can be calculated by separating the total thickness of the layers depending on the desired alloy composition. The individual layer thickness should be thin enough to be mixed efficiently during ion beam mixing, and is 5 nm. Here, four sets of Ni-Nb-Mo nano-multiple-metal-layers [71], i.e., Ni 50  To promote the initial energetic state of the nano metal layers to a high level, 9 or 10 layers were employed to increase the interfaces free energy for each sample.
Afterwards, the A-B-C nano-multiple-metal-layers were prepared by alternatively depositing pure metal A (99.99%), B (99.99%) and C (99.99%) onto the NaCl single crystal substrate at a rate of about 0.3 Å/s. The e-gun evaporation system was under the vacuum level of 10 −6 Pa. The as-deposited nano-multiple-metal-layers were irradiated by xenon ions within a range from 8 × 10 14 to 7 × 10 15 Xe + /cm 2 in an implanter under a vacuum level lower than 5 × 10 −4 Pa. In order to avoid overheating, the sample holder was cooled by liquid nitrogen (77 K), and the current density was limited to be 15 µA/cm 2 during the irradiation process. To characterize the nano metal layer structures, the A-B-C nano-multiple-metal-layers were examined and investigated at room temperature (300 K) by a transmission electron microscopy. Also, the bright field image of samples was analyzed by high resolution transmission electron microscopy. To measure the real composition of the deposited and irradiated nano-multiple-metal-layers, X-ray fluorescence (XRF) and energy dispersive spectroscopy (EDS) were performed in the structure analysis.

Glass Formation Region for the Ni-Zr-Mo System
By performing MD simulations, two different states of Ni-Zr-Mo alloys were obtained: a crystalline state (CS) and an amorphous state (AS). We take two alloys, i.e., Ni 90 Zr 5 Mo 5 and Ni 64 Zr 36 , in the fcc-solid solution model as examples. Figure 3 presents the structure factor S(q) and projection of atomic positions for both states. For the Ni 90 Zr 5 Mo 5 alloy, the S(q) curve in Figure 3a exhibits a typical crystalline feature, which is associated with the atomic position projection in Figure 3b. For the Ni 64 Zr 36 alloy, all the crystalline peaks beyond the second peak have either flattened or disappeared as shown in Figure 3c, matching well with the S(q) curve measured by X-ray diffraction (XRD) [73]. The atomic position projection of Ni 64 Zr 36 in Figure 3d indicates that the original crystalline lattice has spontaneously collapsed and transformed into an amorphous state. It follows that by increasing the solute concentration, the crystalline lattice of the original solid solution model would be seriously distorted, and finally turn into a disordered state. Consequently, the underlying physical mechanism of the crystal-to-amorphous transformation is the collapsing of crystalline lattice in a solid solution model when the solute concentration exceeds the critical solid solubility.
According to the results of MD simulations, we have grouped all the Ni x Zr y Mo 100-x−y alloys into two structural states, i.e., crystalline state (CS) and amorphous states (AS). The constructed glass formation region at 300 K is shown in Figure 4a. The Ni-Zr-Mo composition triangle can be divided into four regions by three critical solubility lines, i.e., AB, CD and EF. When an alloy composition situates beyond the critical lines towards one of the three corners, the Ni-Zr-Mo alloy keeps its original crystalline structure, and all of three corner regions are considered as crystalline regions. When an alloy composition falls into the central hexagonal region enclosed by the critical lines, the crystalline structure becomes unstable and then completely collapses, transforming into a disordered state. This shaded area is therefore identified as the amorphous area, i.e., the GFR of the Ni-Zr-Mo system. In addition, we compare the simulated GFR with the thermodynamic predictions and experimental data. According to the other studies [72], the GFR predicted by the thermodynamic calculation within the framework of Miedema's model and Alonso's method mostly overlaps the simulated GFR. Besides, various experimental data have been collected and marked by different symbols, mostly located within the central hexagonal region. For the Ni-Mo side, the glass formation composition range of 25-68 at.% Mo predicted in Figure 4a presents a minor deviation from that of 25-75 at.% Mo, obtained from a solid-state reaction (SSR) [74]. The deviation may be caused by the factors such as the ideal simulation model, impurities, and chemical and structural fluctuations [75]. For the Ni-Zr side, the mechanical alloying (MA) [76] experiments reveal a glass formation range of 24-83 at.% Zr. Besides, XRD diffractograms of Zr x Ni 100−x thin film MGs produced by magnetron sputtering method [39] indicate that the compositions are lying within the present simulated amorphization range (10-80 at.% Zr), namely Zr 42 Ni 58 , Zr 65 Ni 35 and Zr 75 Ni 25 . Furthermore, ternary Ni-Zr-Mo MGs, e.g., (NiZr) 100−x Mo x , can be fabricated by IBM (marked as red triangles), and completely located within the central shaded region. Therefore, MD simulation can design and predict the glass formation region of the Ni-Zr-Mo system well. located within the central shaded region. Therefore, MD simulation can design and predict the glass formation region of the Ni-Zr-Mo system well.  [56]. The XRD data [73] for Ni64Zr36 MGs are also presented and marked as olive dots. Blue solid circles are for Ni, orange solid circles for Zr, and yellow solid circles for Mo.  , respectively [56]. The XRD data [73] for Ni 64 Zr 36 MGs are also presented and marked as olive dots. Blue solid circles are for Ni, orange solid circles for Zr, and yellow solid circles for Mo. located within the central shaded region. Therefore, MD simulation can design and predict the glass formation region of the Ni-Zr-Mo system well.  [56]. The XRD data [73] for Ni64Zr36 MGs are also presented and marked as olive dots. Blue solid circles are for Ni, orange solid circles for Zr, and yellow solid circles for Mo.

Glass Formation Region for the Ni-Nb-Mo System
Similarly, the whole Ni-Nb-Mo composition region as shown in Figure 4b can be divided into three regions by two CS-AS boundary polylines, i.e., GHIJ and KLM. The central shaded area surrounded by the two polylines GHIJ and KLM is considered as the amorphous region of the Ni-Nb-Mo system, within which metallic glasses are energetically favored to form. By comparison, the GFR of the Ni-Nb-Mo system derived from thermodynamics [71] covers most of the simulated region. Besides, the IBM data [71] for the Ni-Nb system indicate a glass formation range of 15-80 at.% Nb, whereas the RS data [43] present a range of 30-60 at.% Nb. Both glass formation composition regions in the experiment are close to that of 15-75 at.% Nb in Figure 4b. In addition, ternary Ni-Nb-Mo metallic glasses produced by IBM are all located within the GFR, indicating the conformity between the MD simulations and the experimental data.
To further confirm the reliability of metallic glass formation derived from the MD simulations, the structure factor of the simulated Ni 52 Nb 35 Mo 13 alloy was compared with the IBM experiments as shown in Figure 5. The first, second and third peaks are located at the q 1 = 2.75 Å −1 , q 2 = 5.5 Å −1 and q 3 = 8.05 Å −1 , respectively. By analyzing the selected area diffraction (SAD) patterns of the transmission electron microscope (TEM) image, the average diameters of the three halos appearing in the SAD patterns were calculated to be 2.24, 4.21 and 6.09 centimeters, respectively. Afterwards, the wave vector, q can be transformed from the average diameters of the diffused halos, and measured as q a = 2.804 Å −1 , q b = 5.269 Å −1 and q c = 7.622 Å −1 , respectively. By comparing q 1 to q a , q 2 to q b , and q 3 to q c , the maximum errors of the deviations are all lower than 5%, revealing that the S(q) curve corresponds to the SAD patterns. As a result, the atomic structure of the simulated alloy is similar to that of the real alloy produced by the experiments, indicating that the simulation model under the proposed long-range empirical formulas is reasonable.

Glass Formation Region for the Ni-Nb-Mo System
Similarly, the whole Ni-Nb-Mo composition region as shown in Figure 4b can be divided into three regions by two CS-AS boundary polylines, i.e., GHIJ and KLM. The central shaded area surrounded by the two polylines GHIJ and KLM is considered as the amorphous region of the Ni-Nb-Mo system, within which metallic glasses are energetically favored to form. By comparison, the GFR of the Ni-Nb-Mo system derived from thermodynamics [71] covers most of the simulated region. Besides, the IBM data [71] for the Ni-Nb system indicate a glass formation range of 15-80 at.% Nb, whereas the RS data [43] present a range of 30-60 at.% Nb. Both glass formation composition regions in the experiment are close to that of 15-75 at.% Nb in Figure 4b. In addition, ternary Ni-Nb-Mo metallic glasses produced by IBM are all located within the GFR, indicating the conformity between the MD simulations and the experimental data.
To further confirm the reliability of metallic glass formation derived from the MD simulations, the structure factor of the simulated Ni52Nb35Mo13 alloy was compared with the IBM experiments as shown in Figure 5. The first, second and third peaks are located at the q1 = 2.75 Å −1 , q2 = 5.5 Å −1 and q3 = 8.05 Å −1 , respectively. By analyzing the selected area diffraction (SAD) patterns of the transmission electron microscope (TEM) image, the average diameters of the three halos appearing in the SAD patterns were calculated to be 2.24, 4.21 and 6.09 centimeters, respectively. Afterwards, the wave vector, q can be transformed from the average diameters of the diffused halos, and measured as qa = 2.804 Å −1 , qb = 5.269 Å −1 and qc = 7.622 Å −1 , respectively. By comparing q1 to qa, q2 to qb, and q3 to qc, the maximum errors of the deviations are all lower than 5%, revealing that the S(q) curve corresponds to the SAD patterns. As a result, the atomic structure of the simulated alloy is similar to that of the real alloy produced by the experiments, indicating that the simulation model under the proposed longrange empirical formulas is reasonable. Figure 5. The structure factor S(q) of the Ni52Nb35Mo13 alloy derived from MD simulation [55] and the corresponding selected area diffraction (SAD) patterns [71]. The q1, q2 and q3 represent the first, second and third peak positions of the structure factor S(q).

Glass Formation Ability of the Ternary Systems
Theoretically, the energy difference between the solid solution and the amorphous phase serves as the driving force for the crystal-to-amorphous transition. The higher the driving force, the larger the GFA and the easier the metallic glasses can be synthesized, i.e., the issue related to evaluating the GFA of the ternary system can be transformed into calculating the driving force at a specific composition.
Let the Eam, Es.s stand for the average energy per atom of the amorphous phase and solid solution in the ternary A-B-C system. The amorphization driving force (ADF), i.e., energy difference between the solid solution and amorphous phase, can be computed by  [55] and the corresponding selected area diffraction (SAD) patterns [71]. The q 1 , q 2 and q 3 represent the first, second and third peak positions of the structure factor S(q).

Glass Formation Ability of the Ternary Systems
Theoretically, the energy difference between the solid solution and the amorphous phase serves as the driving force for the crystal-to-amorphous transition. The higher the driving force, the larger the GFA and the easier the metallic glasses can be synthesized, i.e., the issue related to evaluating the GFA of the ternary system can be transformed into calculating the driving force at a specific composition.
Let the E am , E s.s stand for the average energy per atom of the amorphous phase and solid solution in the ternary A-B-C system. The amorphization driving force (ADF), i.e., energy difference between the solid solution and amorphous phase, can be computed by The ADFs of the Ni-Zr-Mo and Ni-Nb-Mo systems were calculated and are presented in Figure 6, in addition to the ADF of the (NiZr) 100−x Mo x alloys along the vertical line XY. For the Ni-Zr-Mo system in Figure 6a, the ADF is negative within the whole GFR, suggesting that the formation energy of an amorphous phase is lower than that of the solid solution and the metallic glass formation is energetically favored. Intuitively, the larger the energy difference, the stronger the ADF, and the larger the GFA of the alloy. It was found that the alloys with compositions marked as red dots have a lower ADF than those with compositions in the other regions, suggesting a stronger ADF for the alloys. Within the red dot region, the Ni 45 The ADFs of the Ni-Zr-Mo and Ni-Nb-Mo systems were calculated and are presented in Figure  6, in addition to the ADF of the (NiZr)100−xMox alloys along the vertical line XY. For the Ni-Zr-Mo system in Figure 6a, the ADF is negative within the whole GFR, suggesting that the formation energy of an amorphous phase is lower than that of the solid solution and the metallic glass formation is energetically favored. Intuitively, the larger the energy difference, the stronger the ADF, and the larger the GFA of the alloy. It was found that the alloys with compositions marked as red dots have a lower ADF than those with compositions in the other regions, suggesting a stronger ADF for the alloys. Within the red dot region, the Ni45Zr40Mo15 alloy symbolized by the black pentagram is shown as the composition with the maximum ADF. Additionally, the composition of Ni45Zr40Mo15 obtained from the atomistic simulation is close to that of Ni48Zr48Mo4 obtained from thermodynamic calculation [56]. Therefore, the Ni45Zr40Mo15 alloy and the nearby compositions can be considered energetically favored for metallic glass formation. Figure 6. The formation energy difference between the solid solution and amorphous phase for (a) Ni-Zr-Mo system [56], and (b) Ni-Nb-Mo system [55]. The quantitative energy scale for dots with different colors are presented for the respective systems.
We also computed the ADF for the Ni-Nb-Mo system and the results are presented in Figure 6b. The maximum ADF for the Ni-Nb-Mo alloy is determined to be the composition of Ni55Nb30Mo15. In addition, the composition of Ni56Nb36Mo8 obtained from thermodynamics [71] is also close to that of Ni55Nb30Mo15 predicted by the atomistic simulations. As a result, the Ni55Nb30Mo15 alloy as well as the nearby compositions can be more thermally stable or easily produced than other Ni-Nb-Mo alloys. In particular, the ADF for the Ni-Nb binary system firstly increases by increasing the Nb concentration, and decreases after reaching the maximum ADF of the Ni-Nb alloy. The Ni45Nb55 alloy and the nearby compositions are considered to be most likely to form MGs for the Ni-Nb system. Similarly, the Ni55Mo45 alloy and the nearby compositions are most likely to form MGs for the Ni-Mo system.

Metallic Glass Formation by Ion Beam Mixing
We now present the structural transformation that took place in the Ni-Nb-Mo nano-multiplemetal-layers, which are irradiated with various doses ranging from 8 × 10 14 Xe + /cm 2 to 5 × 10 15 Xe + /cm 2 (shown in Table 5). All compositions of the four Ni-Nb-Mo nano-multiple-metal-layers locate in the glass formation region predicted by the atomistic simulation. For instance, the composition of sample Figure 6. The formation energy difference between the solid solution and amorphous phase for (a) Ni-Zr-Mo system [56], and (b) Ni-Nb-Mo system [55]. The quantitative energy scale for dots with different colors are presented for the respective systems.
We also computed the ADF for the Ni-Nb-Mo system and the results are presented in Figure 6b. The maximum ADF for the Ni-Nb-Mo alloy is determined to be the composition of Ni 55 Nb 30 Mo 15 . In addition, the composition of Ni 56 Nb 36 Mo 8 obtained from thermodynamics [71] is also close to that of Ni 55 Nb 30 Mo 15 predicted by the atomistic simulations. As a result, the Ni 55 Nb 30 Mo 15 alloy as well as the nearby compositions can be more thermally stable or easily produced than other Ni-Nb-Mo alloys. In particular, the ADF for the Ni-Nb binary system firstly increases by increasing the Nb concentration, and decreases after reaching the maximum ADF of the Ni-Nb alloy. The Ni 45 Nb 55 alloy and the nearby compositions are considered to be most likely to form MGs for the Ni-Nb system. Similarly, the Ni 55 Mo 45 alloy and the nearby compositions are most likely to form MGs for the Ni-Mo system.

Metallic Glass Formation by Ion Beam Mixing
We now present the structural transformation that took place in the Ni-Nb-Mo nano-multiple-metal-layers, which are irradiated with various doses ranging from 8 × 10 14 Xe + /cm 2 to 5 × 10 15 Xe + /cm 2 (shown in Table 5). All compositions of the four Ni-Nb-Mo nano-multiple-metal-layers locate in the glass formation region predicted by the atomistic simulation. For instance, the composition of sample 4 in Table 5 is Ni 51 Nb 19 Mo 30 , with high concentrations of Ni and Mo while the concentration of Nb is low. The corresponding SAD patterns of the phase evolution occurring in the Ni 51 Nb 19 Mo 30 nano-multiple-metal-layers upon irradiation are presented in Figure 7. In Figure 7a, the lattice constants of the bcc and fcc structure phase are computed as 3.16 Å and 3.54 Å, respectively, revealing the lattice constant of the constituent Mo and Ni in the system. After irradiation doses of 1 × 10 15 Xe + /cm 2 in Figure 7b, the Ni-rich fcc crystalline phase diffraction ring transforms into a diffuse halo, whereas the Mo-rich bcc crystalline phase still remains. As the irradiation doses increase to 5 × 10 15 Xe + /cm 2 in Figure 7c, the diffused halos without diffraction rings appear, suggesting that the nano-multiple-metal-layers turn into the uniform amorphous phase without the crystalline phase. From Table 5, as the irradiation doses increase, a similar phase transformation also take place in Ni 61 Nb 15 Mo 24 and Ni 72 Nb 20 Mo 8 nano-multiple-metal-layers upon ion beam mixing. The uniform amorphous phase is obtained by ion beam mixing experiments upon irradiation to an appropriate dosage. At lower ion dosage, the amorphous and crystalline phases coexist, while the uniform amorphous phase is obtained by increasing the ion dosage to a higher value.  Table 5 is Ni51Nb19Mo30, with high concentrations of Ni and Mo while the concentration of Nb is low. The corresponding SAD patterns of the phase evolution occurring in the Ni51Nb19Mo30 nanomultiple-metal-layers upon irradiation are presented in Figure 7. In Figure 7a, the lattice constants of the bcc and fcc structure phase are computed as 3.16 Å and 3.54 Å, respectively, revealing the lattice constant of the constituent Mo and Ni in the system. After irradiation doses of 1 × 10 15 Xe + /cm 2 in Figure 7b, the Ni-rich fcc crystalline phase diffraction ring transforms into a diffuse halo, whereas the Mo-rich bcc crystalline phase still remains. As the irradiation doses increase to 5 × 10 15 Xe + /cm 2 in Figure 7c, the diffused halos without diffraction rings appear, suggesting that the nano-multiplemetal-layers turn into the uniform amorphous phase without the crystalline phase. From Table 5, as the irradiation doses increase, a similar phase transformation also take place in Ni61Nb15Mo24 and Ni72Nb20Mo8 nano-multiple-metal-layers upon ion beam mixing. The uniform amorphous phase is obtained by ion beam mixing experiments upon irradiation to an appropriate dosage. At lower ion dosage, the amorphous and crystalline phases coexist, while the uniform amorphous phase is obtained by increasing the ion dosage to a higher value.    [71] (a) at as-deposited state, after being irradiated to a dose of (b) 1 × 10 15 Xe + /cm 2 and (c) 5 × 10 15 Xe + /cm 2 . The experimental observation of amorphous phase formation can be explained by destroying the structure of crystalline lattices by ion beams. The stable structure of crystalline lattice makes it difficult for the amorphous phase to form. For example, if high ion dosage is needed to destroy the crystalline lattice and form the uniform amorphous phase, the GFA of the corresponding alloy can be considered as poor under the same circumstance (e.g., beam flux, accelerating voltage, etc.). If a low ion dosage is required to produce amorphous alloys, the GFA of an alloy is regarded as good. No matter whether the GFA of an alloy is good or poor, there exists a minimum required ion dosage D m , below which no uniform amorphous phase can be formed [77]. Furthermore, there exists a negative correlation between the GFA of an alloy and the D m , i.e., the lower the minimum required ion dosage, the better the GFA of an alloy. To examine this negative correlation, it can be seen from Table 5 that Ni 52 Nb 35 Mo 13 amorphous alloy could be easily formed, and possesses the lowest D m among all of the Ni-Nb-Mo nano-multiple-metal-layers. Coincidentally, the composition of Ni 52 Nb 35 Mo 13 also presents the largest ADF calculated by atomistic simulations, of all the samples, indicating the highest GFA. Therefore, the negative correlation is not only supported by IBM experiments, but also proven by atomistic simulation results. Table 6 summarizes the phase evolution that takes place in the four sets of the Ni-Zr-Mo nano-multiple-metal-layers. For the (NiZr) 60 Mo 40 nano-multiple-metal-layers, the SAD patterns and corresponding bright field image are presented in Figure 8a,b, respectively, after the ion dosage was increased to 7 × 10 15 Xe + /cm 2 . As can be seen from Figure 8a, there exists a diffused halo reflecting the amorphous phase, which is also verified by the gray matrix of uniform amorphous phase in Figure 8b. Nevertheless, after the Mo concentration increased to 60 at.%, no uniform amorphous phase could be obtained in the (NiZr) 40 Mo 60 nano-multiple-metal-layers until irradiated to the dose of 7 × 10 15 Xe + /cm 2 . The SAD patterns in Figure 8c show an amorphous halo as well as several diffraction rings derived from the Mo-rich bcc crystalline phase with the lattice constant indexed as 3.12 Å. The corresponding bright field image confirms the dual-phase nature, i.e., an amorphous-crystalline coexisting compound. 8b. Nevertheless, after the Mo concentration increased to 60 at.%, no uniform amorphous phase could be obtained in the (NiZr)40Mo60 nano-multiple-metal-layers until irradiated to the dose of 7 × 10 15 Xe + /cm 2 . The SAD patterns in Figure 8c show an amorphous halo as well as several diffraction rings derived from the Mo-rich bcc crystalline phase with the lattice constant indexed as 3.12 Å. The corresponding bright field image confirms the dual-phase nature, i.e., an amorphous-crystalline coexisting compound.

Connections between ADF and GFA
The ADF predicted from the atomistic simulations is also connected with the experimental data. From Figure 9a, it can be seen that the ADF of (NiZr) 100 − x Mo x alloy starts to increase gradually by adding an appropriate amount of Mo, but decreases after reaching the peak value with further addition of the Mo content. When the Mo content is higher than 60 at.%, the crystalline lattice remains stable and no uniform amorphous phase can be produced, presenting the GFR of (NiZr) 100−x Mo x as 0-60 at.% Mo. By comparing the results of ion beam mixing experiments in Table 6, it can be seen that as the Mo concentration increases, the uniform amorphous phase could be formed in (NiZr) 100−x Mo x (x = 10, 20, 40), with its compositions all locating within the predicted GFR of the Ni-Zr-Mo system. This indicates that the simulation prediction is relevant to the experimental data. Nevertheless, when the Mo concentration increases to 60 at.%, the (NiZr) 40 Mo 60 alloy becomes a mixture of the amorphous phase and Mo-rich bcc crystalline phase. The above experimental observation can be explained by the predicted GFR. Since the composition of (NiZr) 40 Mo 60 is near the critical solid solubility line AB, no uniform amorphous phase could be formed. In other words, excessive adding of the Mo concentration would depreciate the GFA of an alloy.
Compared to the (NiZr) 100−x Mo x alloys, a ADF changing rule similar to that of the Mo addition to the Ni (35−x/2) Zr (65−x/2) Mo x and Ni (65−x/2) Zr (35−x/2) Mo x alloys can be observed in Figure 9b.Therefore, it is of great importance to explain the effect of Mo addition on the glass formation of Ni-Zr alloy from a theoretical perspective. There are two dominant explanations for the observed phenomenon. On the one hand, given that the atomic radius of Mo (r Mo = 1.39 Å) is between those of Ni (r Ni = 1.25 Å) and Zr (r Zr = 1.58 Å), the Mo addition promotes three elements with different crystalline structures to regulate the coordination polyhedron. Therefore, it improves the formation of short-range compositional order [78]. In addition, the mixing heats ∆H f of the Zr-Mo and Ni-Mo systems are -9 and -11 kJ/mol, respectively. As the negative mixing heats ∆H f of Mo with Zr and Ni force the Mo atoms to fully spread in the Ni-Zr matrix, it is able to promote interaction among the constituent elements and stabilize the amorphous phase [79], i.e., improve the GFA. On the other hand, a negative mixing heat ∆H f could also favor compound formation, leading to restricted diffusion and counteracting glass formation [80]. Besides, as the melting temperature of Mo is 2850 K, such a high melting point would stabilize the crystalline structure and make the uniform amorphous phase formation difficult, i.e., it reduces the GFA. Once the influence of a high melting temperature is beyond that of size difference and the heat of formation, no uniform amorphous phase can be energetically fabricated, downgrading the GFA of an alloy by further increasing the Mo content.
The ADF predicted from the atomistic simulations is also connected with the experimental data. From Figure 9a, it can be seen that the ADF of (NiZr)100 − xMox alloy starts to increase gradually by adding an appropriate amount of Mo, but decreases after reaching the peak value with further addition of the Mo content. When the Mo content is higher than 60 at.%, the crystalline lattice remains stable and no uniform amorphous phase can be produced, presenting the GFR of (NiZr)100−xMox as 0-60 at.% Mo. By comparing the results of ion beam mixing experiments in Table 6, it can be seen that as the Mo concentration increases, the uniform amorphous phase could be formed in (NiZr)100−xMox (x = 10, 20, 40), with its compositions all locating within the predicted GFR of the Ni-Zr-Mo system. This indicates that the simulation prediction is relevant to the experimental data. Nevertheless, when the Mo concentration increases to 60 at.%, the (NiZr)40Mo60 alloy becomes a mixture of the amorphous phase and Mo-rich bcc crystalline phase. The above experimental observation can be explained by the predicted GFR. Since the composition of (NiZr)40Mo60 is near the critical solid solubility line AB, no uniform amorphous phase could be formed. In other words, excessive adding of the Mo concentration would depreciate the GFA of an alloy. Compared to the (NiZr)100−xMox alloys, a ADF changing rule similar to that of the Mo addition to the Ni(35−x/2)Zr(65−x/2)Mox and Ni(65−x/2)Zr(35−x/2)Mox alloys can be observed in Figure 9b.Therefore, it is of great importance to explain the effect of Mo addition on the glass formation of Ni-Zr alloy from a theoretical perspective. There are two dominant explanations for the observed phenomenon. On the one hand, given that the atomic radius of Mo (rMo = 1.39 Å) is between those of Ni (rNi = 1.25 Å) and Zr (rZr = 1.58 Å), the Mo addition promotes three elements with different crystalline structures to regulate the coordination polyhedron. Therefore, it improves the formation of short-range compositional order [78]. In addition, the mixing heats ΔHf of the Zr-Mo and Ni-Mo systems are -9 and -11 kJ/mol, respectively. As the negative mixing heats ΔHf of Mo with Zr and Ni force the Mo atoms to fully spread in the Ni-Zr matrix, it is able to promote interaction among the constituent

Conclusions
The long-range empirical potential was constructed for the Ni-Nb-Mo (fcc-bcc-bcc) and Ni-Zr-Mo (fcc-hcp-bcc) systems, especially for the Ni-Zr-Mo system consisting of three different crystalline structures. Based on the constructed n-body potential, atomistic simulations not only revealed the underlying process of metallic glass formation as the spontaneous collapse of the crystalline lattice while the solute content exceeds the critical value, but also predicted the composition range energetically favored for metallic glass formation. Besides, the amorphization driving force, defined as the energy difference between the solid solution, and the amorphous phase for each alloy was computed and found to be positively correlated with glass forming ability, which is frequently defined in terms of critical size or cooling rate, etc.
To prove the atomistic simulation prediction for metallic glass formation, ion beam mixing of nano-multiple-metal-layers was performed to produce Ni-Nb-Mo and Ni-Zr-Mo metallic glass films. The composition of these ternary metallic films produced by ion beam mixing all located within the predicted glass formation region. Ion-beam mixing experiments indicated that there exists a minimum ion dosage for amorphization and the minimum ion dosage showed somewhat negative correlation with the glass forming ability. It was also found that the amorphization driving force of (NiZr) 100−x Mo x alloy gradually increases with the proper addition of Mo concentration, whereas it decreases with excessive addition of the Mo content. Therefore, the amorphization driving force of each alloy can be defined as an indicator of its glass forming ability.