Nanometer-scale gradient atomic packing structure surrounding soft spots in metallic glasses

The hidden order of atomic packing in amorphous structures and how this may provide the origin of plastic events have long been a goal in the understanding of plastic deformation in metallic glasses. To pursue this issue, we employ here molecular dynamic simulations to create three-dimensional models for a few metallic glasses where, based on the geometrical frustration of the coordination polyhedra, we classify the atoms in the amorphous structure into six distinct species, where “gradient atomic packing structure” exists. The local structure in the amorphous state can display a gradual transition from loose stacking to dense stacking of atoms, followed by a gradient evolution of atomic performance. As such, the amorphous alloy specifically comprises three discernible regions: solid-like, transition, and liquid-like regions, each one possessing different types of atoms. We also demonstrate that the liquid-like atoms correlate most strongly with fertile sites for shear transformation, the transition atoms take second place, whereas the solid-like atoms contribute the least because of their lowest correlation level with the liquid-like atoms. Unlike the “geometrically unfavored motifs” model which fails to consider the role of medium-range order, our model gives a definite structure for the so-called “soft spots”, that is, a combination of liquid-like atoms and their neighbors, in favor of quantifying and comparing their number between different metallic glasses, which can provide a rational explanation for the unique mechanical behavior of metallic glasses. Liquid-like atoms in a gradient atomic packing structure might determine plasticity in metallic glasses. A team led by Yanqing Su at the Harbin Institute of Technology in China and Robert Ritchie at the University of California, Berkeley, in the USA, used molecular dynamics simulations to classify atoms in different metallic glasses according to their stacking, thereby accounting for the amorphous medium-range order in metallic glasses. They found that local atomic structures gradually transitioned from a loose to a dense stacking of atoms via liquid-like, transition, and solid-like regions, and that liquid-like atoms and their neighbors were equivalent to ‘soft spots’ and associated with initiation of local irreversible atomic arrangements. Modeling this gradient atomic packing structure may help us better understand and design plasticity in glassy alloys.


INTRODUCTION
Compared to their crystalline counterparts, metallic glasses (MGs) are vitrified solids in the metastable state. [1][2][3] The atomic structure of such amorphous matter and its relevance to mechanical behavior is a fundamental and intriguing problem. As reflected in the central materials science paradigm that "structure determines properties", it is important to understand any hidden order in seemingly disordered glassy alloys in order to establish a causal link between such local structure and macroscopic properties. Such an objective, which is universally viable for MGs, has yet to be successfully achieved in condensed matter physics. 4,5 Despite a lack of long-range order, MGs do possess short-range (SRO) and medium-range order (MRO) that have been extensively characterized by experiments 6-10 and atomic simulations. 5,11,12 Many models which describe the general basic principles of an amorphous structure, such as polytetrahedral packing model 13,14 and the efficient cluster packing model, 15,16 have established some insight into the local order in MGs. However, both of these models fail to give an exact pattern for MRO, and cannot give a global description of the amorphous structure. The extended order of atomic packing, on the length-scale of several nanometers and beyond, has become a topic of modeling interest, but with many aspects still uncertain. 5,17,18 The internal structure of MGs is inherently inhomogeneous, although macroscopically uniform, with any structural heterogeneity inevitably leading to dynamic heterogeneity. [19][20][21][22][23] By employing a coarse-grained model, recent computer simulations have demonstrated that there is stiffness fluctuation in MGs. 21,22 This mechanical non-uniformity, over areas much larger than in any simulation, has also been confirmed by experiment using atomic force microscopy (AFM) techniques. 23 It is believed that such nanoscale heterogeneity plays an important role in the macroscopic mechanical behavior of MGs. 5,22 The MRO, specifically the organization of the SRO motifs, is a key factor that dictates the nature of the heterogeneity in mechanical properties observed in MGs. It has been demonstrated that the MRO covers a length-scale comparable to the size of a shear transformation zone (STZ), which in MGs represents the basic event associated with plastic deformation. [24][25][26] The detected "soft spots", defined as aggregates of atoms that intensively tend to be engaged in low-frequency vibrational modes, exhibit a higher propensity to record such local rearrangements. 24,27 However, the structural feature of these "soft spots" is only based on the fraction of "geometrically unfavored motifs" (GUMs) with no clear boundary 25 ; there is no allowance for how they interconnect with one another and what the adjacent motifs are, i.e., any one-to-one correspondence has not been delineated. 17 Therefore, the ambiguous structure of these "soft spots" makes it difficult to define their precise role, i.e., their existence is not currently conducive to explaining the origin of the unique mechanical properties, such as high strength coupled with the large plastic strains (good ductility and fracture toughness), observed in certain microalloyed MGs. 28 Accordingly, the understanding of the structural origins of the deformation behavior of monolithic MGs still requires the contribution of MROs to be defined. Indeed, the limited understanding about the role of MRO in the formation of "soft spots" has compromised the definition of any link between the mechanical properties of MGs and the heterogeneity inherent in amorphous structures.
To this end, the current study focuses on the hidden order in amorphous structures, on the length-scale of a couple of nanometers (MRO) and beyond, with the rationale of defining the structural feature of "soft spots". We examine different kinds of amorphous alloys using molecular dynamics (MD) simulations with an embedded-atom method (EAM) potential. Our structural analysis shows that a powerful spatial distribution pattern within the nanometer-scale, termed as "gradient atomic packing structure", which considers all kinds of atoms rather than several specific ones, exists in the seemingly disordered atomic arrangement of the MGs. According to this particle spatial distribution pattern, the microscopic features and anelastic behavior of welldistinguished regions in MGs are revealed and discussed. Further, we describe the fundamental structural origin of the increased strength and deformability after microalloying in MGs, based on the discernible and comparable numbers of "soft spots".

Classification of total atoms in MGs
All atoms in MGs are divided into several different classes in terms of the geometrical frustration of their corresponding coordination polyhedra at the shortest length-scales. The non-directional metallic bonding and high density of MGs connote a high degree of SROs such that the basic packing unit is mostly a tetrahedron with four atoms. Thus, we can distinguish and study the coordination polyhedra based on polytetrahedral packing that incorporate disclinations. To fill three-dimensional (3D) space, rotational "defects" (termed four-fold and six-fold disclinations) are required for bulk polytetrahedral packing. 14 Additionally, for Voronoi polyhedra with polytetrahedral packing, the condition of 2n 4 + n 5 = 12 must be satisfied, 29 where n i denotes the number of i-edged polygon in Voronoi polyhedra. Meanwhile, the single disclination is prohibited so that <0, 1, 10, 0> and <0, 0, 12, 1> do not exist. 29 In our work, we term these clusters as Kasper polyhedra. According to the different disclination density of each Kasper polyhedra, the relevant ones for MGs are divided into five groups, 5,17 which are listed in Table 1 for coordination numbers (CNs) ranging from 8 to 17. A comprehensive review of polytetrahedral packing with disclinations can be found in ref. 5 For efficient packing of in a solute-centered cluster, the Z clusters with fewest disclinations are the most densely and geometrically stable configuration, 5,30 with the favored CN being determined by the effective size ratio of solute and solvent. 31,32 For example, the local populous motif for Cu is the Z12 full icosahedra (<0, 0, 12, 0>), and Z16 type (<0, 0, 12, 4>) for Zr-centered polyhedra in the Cu 64 Zr 36 alloy. 33 With increasing disclination density, the clusters become more distorted and loose. Apart from Kasper polyhedra, in MGs there are inevitably non-Kasper ones as well, such as <1, 2, 5, 4>, <1, 0, 9, 3>, <1, 2, 6, 3>, etc. Indeed, the clusters found in MGs can be grouped into six types. The non-Kasper polyhedra, whose basic packing unit is not all of tetrahedral structure, will also have a lower packing efficiency. Here we adopt Roman numerals to represent the atoms at the centers of the corresponding coordination polyhedra, as shown in Table 1.
Gradient atomic packing structure in MGs Note that some scholars found that similar clusters tend to group together into nanoscale zones by combining fluctuation electron microscopy experiments and hybrid reverse Monte Carlo simulations. 34,35 Here we gain insight into the statistical correlation between these six kinds of atoms in as-quenched and annealed samples by calculating a nearest-neighbor correlation index, 36 where p ij is the probability of atom i and j to be the nearest neighbors and p 0 ij is the probability for a model that the distributions of atoms are spatially uncorrelated. The p ij is given by p ij = m ij /p total , where m ij and p total denote, respectively, the number of nearest-neighbor pairs of type i and j and the total number of nearest-neighbor atomic pairs. The value of p 0 ij is calculated by where n i (n j ) is the number of i (j) types of atoms, respectively, and N is the total number of atoms. As such, a random distribution of atom i and j implies a zero value of C ij , and positive or negative values of C ij , respectively, indicates a preference or avoidance of atoms with i and j being the nearest neighbors. Figure 1a, b, respectively, show the correlation matrix of C ij for the six kinds of atoms in as-quenched and annealed Cu 64 Zr 36 ; the results show that the correlation intensity between different types of atoms presents a distinct gradient variation, a result that is highlighted by black arrows (from left to right in Table 1) and red arrows (from right to left in Table 1) in Fig. 1a. The similar appearance can also be seen in Fig. 2, which presents a correlation matrix of C ij for the atoms centered in different coordination polyhedra (both Cu-centered and Zr-centered ones) within the asquenched Cu 64 Zr 36 sample. The colorized map is segmented into a grid of 6 × 6 based on the classification presented in Table 1. It is of note that the VI atoms centered in non-Kasper polyhedra have a variety of Voronoi indices; only five ones are listed whose fraction exceeds 1%. The I atoms centered in Z clusters, which account for a main part in real MGs, 30,33,37 tend to gather together; these form the backbone of the amorphous alloy. According to the correlation level with I atoms, three different regions can be selected, as marked by dashed lines in Fig. 1b. Figure 1c, d, respectively, display three selected clusters in asquenched state and a thin slab in annealed state, with different (false) colors representing the six kinds of atoms. It is clear that a regular pattern of atomic arrangement can be observed directly in local nanometer-scale regions-we term this as "gradient atomic packing structure", i.e., the I atoms tend to gather into dense areas, the VI atoms collect with IV and V atoms to form loose areas, together with II and III atoms located between these two areas. In Supplementary Fig. S1 and S2, we also give the correlation matrix of C ij for the six kinds of atoms in some other glassy alloys, such as the Mg-Cu-Y system, where a parallel pattern can also be found. Thus, we believe that the "gradient atomic packing" structure is a characteristic of glassy alloys.
As we know, a tetrahedron with equilateral triangle faces has the highest packing efficiency among the Bernal's canonical holes. 38 However, in high-temperature amorphous liquids, the local atomic packing is extremely loose and the VI atoms, whose coordination shell is not all triangulated (2n 4 + n 5 ≠ 12), will dominate. Upon cooling, the density of tetrahedral structure gradually increases; furthermore, the spatial aggregation of tetrahedra leads to many five-edged faces in Voronoi polyhedra. Thus, the I atoms whose coordination polyhedra has maximum possible number of fivefold bonds n 5 at the corresponding CN should aggregate to form stable regions in the final alloys. Moreover, as the high cooling rate required for the formation of a MG inevitably will result in insufficient time for the above process to occur, the VI atoms must be inherited from the liquid state to form unstable regions. As there are large structural differences between the stable and unstable regions, there should be an appropriate transition to adjust this mismatch. Thus, based on an  Table 1. The white arrows highlight that the "gradient nanostructure" existed in MGs increasing disclination density of coordination polyhedra from II to V atoms, a gradient atomic packing local structure will emerge to help atoms fill the 3D-space naturally and stably.
Gradient evolution of atomic performance The fractions of the six kinds of atoms, as a function of temperature during continuous cooling, are plotted in Fig. 3a. It is found that the collective IV, V, and VI atoms in the final MGs are naturally inherited from the liquid phase; accordingly, we refer to these unstable atoms as liquid-like atoms as they constitute liquidlike regions. Conversely, the most stable and collective I atoms increase in number sharply in the supercooled liquid region, relative to the other two incremental ones; we term these atoms as solid-like atoms as they constitute solid-like regions. Correspondingly, the II and III atoms are termed transition atoms as they constitute transition regions.
In Fig. 3b, we present the tendency of structural evolution during the aging of a supercooled liquid for 400 ns at 800 K. Among all the atoms, only solid-like atoms (I atoms) experience a persistent and pronounced increase in population at the cost of others. The lower the correlation intensity with I atoms, the higher the extent of reduction for corresponding atoms, especially the VI ones, can be achieved. With regard to Fig. 3a, upon continuous cooling, the liquid-like atoms are not only transformed into solidlike atoms, but also turned into transition atoms. However, during the annealing of the supercooled liquid, the fraction of transitional II atoms holds constant, together with a diminishing proportion of III atoms, which means a dynamic balance must exist for this process. Figure 4 provides several selected clusters at the different annealing moments to view the dynamic transformation of atomic types, i.e., the transition atoms (II and III atoms) act as intermediates during the transformation process from liquid-like atoms to solid-like atoms. Therefore, we can deduce that the trajectory of the glass transition is from liquid-like atoms, through transition atoms, to solid-like atoms; this is the dynamic factor for the "gradient atomic packing structure" formed in MGs, as illustrated in Supplementary Fig. S4. Taking this one step further, the gradient characteristic of amorphous structure may be the true essence of the nature of glass transition. Furthermore, it can be expected that MGs formed in laboratory experiments at several orders of magnitude lower cooling rates will contain even more solid-like atoms, resulting in a more ordered and stable alloy.
Directing our attention to the potential energy (E) of atoms within different groups, we introduce differences between potential energy E j i of each atom and mean potential energy 〈E j 〉 over the same atom type, constructed as δ i ðEÞ ¼ E j i À hE j i, where i delegates the atom index and j denotes the element type. In Fig. 3c, from left to right, each solid bar represents a bin that contains a number of atoms, in an ascending order from the lowest to the highest δ(E). Interestingly, the share of each type of  Table 1 Nanometer-scale gradient atomic packing structure surrounding soft B Wang et al. atoms in these bars also follows the gradient behavior; the solidlike atoms dominate the low δ(E), the transition atoms hold an average scale in each bar, whereas the liquid-like atoms dominate the high δ(E). Among transition atoms, because of the difference in correlation level with solid-like or liquid-like atoms (Fig. 1a), the II atoms take a higher proportion in low δ(E) than the III ones. The critical issue here is to try and relate such a gradient structure to the soft spots in MGs. We examine the response of such internal structure to the external stress stimulus. A local minimum non-affine displacement 39 D 2 min À Á was adopted to identify the atoms involved in any local irreversible rearrangement, such as a shear transformation (ST). 40,41 This anelastic strain of atoms in different groups was tracked during shear deformation. The shear stress-strain curves for the as-quenched and annealed Cu 64 Zr 36 samples are plotted in Supplementary Fig. S5. In Fig. 3d (Supplementary Fig. S6), we show how D 2 min develops with increasing strain, such that the gradient evolution can be deduced. Apparently, the stable solid-like atoms deform in a more affine manner, with less atomic rearrangement, as compared to the higher incidence of STs which are promoted with the higherenergy liquid-like atoms; the transition atoms connecting these solid-like and liquid-like atoms are always in an intermediate state.
Similarly, among transition atoms, the average D 2 min for II atoms is less than that for III atoms, following the difference in δ(E) between these two kinds of atoms. Figure 5 (Supplementary Figs. S7 and S8) map the 2D spatial distribution of the solid-like, transition, and liquid-like regions inside slabs of the as-quenched (annealed) sample, each with a thickness of 2.5 Å, roughly equivalent to the average atomic spacing. Structural heterogeneity, with an approximate correlation length of~1-2 nm, can be deduced which is a little smaller than the spatial resolution mapped by AFM 20 ; such a discrepancy should be attributed to the inherent limitations of cooling rates (fast and uniform) and sample size (nanoscale) in MD simulations. For a direct comparison, the local motifs of atomic rearrangements in sheared MGs are superimposed in the contoured maps, where the white circles represent atoms which have experienced the most accumulative non-affine strains (or clear STs) after global strains of 2% (Figs. 5a) and 4% (Fig. 5b). In our simulations, we chose a critical value of D 2 min (e.g., 0.1, with reference to the initial configuration) as the threshold to judge whether atoms are taken to be the ones that have experienced clear STs D 2 min ! 0:1 Arguably, the "gradient atomic packing structure", beyond local SRO, comes into play in controlling the nucleation and propagation of STs. The visual correlation in Fig. 5 (Supplementary Figs. S7 and S8) establishes that: (i) the solid-like regions, especially the sites enriched a large proportion of I atoms, deform in a more affine manner with the least plastic events; (ii) the STs have an increased tendency to originate from the transition regions; and (iii) the liquid-like regions are most likely to record the non-affine displacement than the other regions. In addition, comparing Fig.  5b with Fig. 5a, the propagation of STs, which are trapped by solid-like regions ("backbone"), prefers to select liquid-like regions first, and then transition regions, as highlighted by the arrows in Fig. 5b. It should be noted here that not all liquid-like regions would experience STs for a given loading regime, as the stress field (tensor) and loading direction are extrinsic factors which influence the local strain response of atoms in additional to the intrinsic atomic structure. 42,43 More importantly, as different kinds of regions own unique atomic arrangements, the interface between them can be considered to be essentially equivalent to a grain boundary in crystalline solids, i.e., the atoms adjoined to the interface are likely to act as nucleation sites for STs, as presented by Fig. 5a and Supplementary Fig. S7. The solid-like or transition atoms, which also tend to display anelastic relaxation, must be located close to the liquid-like regions, where the coordination environment is characterized by a high distortion energy.

DISCUSSION
In ref. 24 the structure of "soft spots" is determined simply based on the fraction of GUMs, which are generally defined as the polyhedra that deviate most significantly from the geometrically preferable clusters. Especially, this definition of "soft spots" does not consider the contribution of MRO such that the neighboring motifs surrounding these regions are undefined. Based on the above results, the gradient atomic packing structure in MGs, covering the MRO, provides a means to redefine the composition of these "soft spots"; specifically, we refer to the regions containing liquid-like atoms and their neighbors as the "soft spots", as presented by yellow circles in Fig. 6. Therefore, we can regard such regions as locations where the local structure is ideally activated to induce relaxation and deformation. Further, as the "soft spots" must include liquid-like regions, from where there is a one-to-one correspondence between these two objects, we can adopt the well-recognized liquid-like regions to precisely represent the "soft spots". Next, by inspecting Fig. 5 (Supplementary Figs. S7 and S8), it is apparent that the dimensions of these liquid-like regions display a non-uniform distribution, with the larger bulk regions not always providing preferred sites for STs, i.e., each region has the equal chance of being a nucleation site for a plasticity event. Consequently, we can actually argue that how many shear bands form in MG during deformation only depends on the number of liquid-like regions in it; however, the number of these regions is determined, not just by the fraction of liquid-like atoms, but by the spatial distribution of these atoms. For the sake of understanding, two idealized scenarios are employed to make the comparison, as shown in Fig. 6. Each arrow hypothetically represents one direction of a shear band. We can presume that a ductile MG, whose liquid-like atoms distribute evenly, will form more liquid-like regions ("soft spots"), further displaying more fertile sites for STs, and forming multiple shear bands; in contrast, the brittle sample will be one whose liquid-like atoms gather together into one region such that it will merely form single shear band, which accounts for its limited ductility. Recently, Wang et al. 44 utilized a mesoscale shear-transformation zone dynamics modeling framework to explore the effect of spatial correlation of the nanoscale elastic heterogeneity on the mechanical behavior of MGs; a critical correlation length was identified at which the ductility of amorphous alloy is optimum. Above the critical value, as the correlation length increases, the enlarged "soft spots", although reduced in number, enhance the extent of strain localization and further decrease the number of shear bands, which is consistent with our statement above. Conversely, below the critical value, the increased correlation length facilitates the formation of shear bands. Hence, it can be deduced that for a given MG, a critical size of "soft spots" is needed for shear bands to nucleate, below which the material will display poor toughness properties, despite the fact that there may be a large number of "soft spots" present. Further work should investigate how to predict this critical size of "soft spots" for the nucleation of shear bands.
It is generally accepted that an amorphous alloy with lower atomic-packing densities (liquid-like) will show greater plasticity; conversely, higher glass-forming ability (GFA) and strength are invariably associated with higher atomic-packing densities (solidlike). However, specific microalloyed MGs have high GFA and high strength yet still display large plastic strains, as for example with the addition of minor amounts of Al to Cu-Zr binary alloys ( Table  2). In Fig. 7a, we illustrate the shear stress-strain response for the as-quenched Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 samples. The simulated mechanical behavior is in good agreement with experimental results, such as a larger yield strain, higher strength, and elastic modulus after minor additions of Al. Meanwhile, the local shear strain for each atom was monitored during the deformation; only atoms with large local strain (anelastic atoms, threshold is 0.1) are shown in Fig. 7b. As the shear strain increases, the variation in the fraction of these anelastic atoms in both samples can be seen in Fig. 7a. At the initial elastic stage, compared to Cu 50 Zr 50 , more evenly distributed plastic sites occur in Cu 47.5 Zr 47.5 Al 5 , even though the fraction of anelastic atoms is almost the same. With increase in strain, the deformed Cu 50 Zr 50 contains more anelastic atoms, which are more likely to gather together, until the yield point is reached; conversely, a relatively slower strain concentration in Cu 47.5 Zr 47.5 Al 5 should be responsible for the delayed yield (larger yield strain) phenomenon. Beyond the yield point, in turn, there are more anelastic atoms, but with a relatively uniform The critical diameter d c is taken from ref. 43 Nanometer-scale gradient atomic packing structure surrounding soft B Wang et al.
distribution, in the yielding Cu 47.5 Zr 47.5 Al 5 alloy. In general, minor additions of Al can facilitate more sites for STs to further weaken the degree of strain localization during the deformation. Figure 8 presents a comparison of the proportion of different types of atoms between the Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 glassy alloys at different cooling rates. In order to eliminate statistical error, we employ 100 atomic configurations to calculate averaged structural parameters. Compared to Cu 50 Zr 50 , Cu 47.5 Zr 47.5 Al 5 has a higher fraction of solid-like atoms (especially the atoms centered in full icosahedra, as shown in Supplementary Fig. S9), yet conversely, has a slightly lower proportion of liquid-like atoms. Accordingly, this ternary MG possesses higher atomic-packing densities and stronger stiff backbone, which contribute to its higher GFA and strength. Moreover, the transition atoms in our "gradient atomic packing structure" model act as an intermediary such that the change in number of solid-like and liquid-like atoms will not be precisely the opposite, i.e., an MG can have a higher fraction of the solid-like atoms without sacrificing its share of the liquid-like ones.
To explain the reason for the weaker strain localization in Cu 47.5 Zr 47.5 Al 5 , we firstly pay attention to the spatial correlation length 20,45 of the liquid-like atoms in these two samples. The autocorrelation function is calculated by the program OVITO, 46 which employs the fast Fourier transform algorithm to compute the convolution. In Fig. 9, the correlation coefficient c (normalized by covariance) for the liquid-like atoms in Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 samples is shown; the value c decays exponentially with respect to the distance r from a reference position, fitted by the two-phase exponential decay function, cðrÞ ¼ c 0 þ A 1 e Àr=a1 þ A 2 e Àr=a2 , where c 0 is the offset, A 1 (A 2 ) is the amplitude, and a 1 (a 2 ) is the decay constant. Clearly, at different cooling rates, the liquid-like atoms in Cu 50 Zr 50 always possess a larger correlation length and are more likely to aggregate to form larger liquid-like regions.
In addition, we also employ the inhomogeneity parameter h 47 to quantitatively assess the spatial distribution of the liquid-like atoms in Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 alloys. To achieve this parameter, a given 3D box was subdivided into grids with the In order to eliminate the dependence of the inhomogeneity parameter on the grid number generated in different systems, a weighted-sum technique was employed. The edges of the 3D box were divided into 2, 4, 8, …, 2 r segments with the same length; the corresponding numbers of grids accordingly were 2 3 , 4 3 , 8 3 , …, 2 3r . The value of r should meet the criterion that each grid contains at most one atom. Thus, the inhomogeneity parameter h can be expressed as 47 where N is the number of investigative atoms in the box, w is the weight factor whose value is about 8.89, m i is the actual number of atoms in each grid, and m is the mean value of atoms calculated for 2 3r grids. To calculate the value of h for the liquid-like atoms in each Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 sample, a value of r was either taken as 5 or 6. The value of h ascends with increasing degree of inhomogeneity for the distribution of characteristic atoms. Of the two samples, the h value calculated for Cu 47.5 Zr 47.5 Al 5 is relatively small, as displayed in Fig. 10. In addition, the h value for both samples increases with decreasing the cooling rate, but the difference in the value of this parameter for the Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 glasses becomes larger. That is to say, the liquid-like atoms in Cu 47.5 Zr 47.5 Al 5 formed in laboratory experiments with several orders of magnitude lower cooling rate should distribute more uniformly than in Cu 50 Zr 50 MG processed under the same experimental conditions. Based on the above results, a smaller correlation length of the liquid-like atoms with a more uniform distribution must lead to more liquid-like regions ("soft spots") in Cu 47.5 Zr 47.5 Al 5 MG. Taking this one step further, minor doping with proper Al atoms can form more "soft spots", weaken the strain localization, and then promote shear band formation, 48 which in turn improves the ductility of the Cu 50 Zr 50 -based alloy.
In brief, the amorphous structures that contain a higher fraction of solid-like atoms coupled with more evenly distributed liquidlike atoms will not only reinforce the "backbone", but also generate more fertile sites for plasticity events; as this is the fundamental essence of plastic deformation in glassy solids, such amorphous structure will in turn simultaneously display higher strength, higher ductility, and improved fracture toughness properties. 28,48 Finally, we should note that although previous studies 5,16,37 have shown that different MGs are comprised of different favored motifs (Z clusters) due to their respective atomic size ratios, the "gradient atomic packing structure" is a common feature associated with various glassy systems. Consequently, our findings here should be representative of all amorphous metal alloys.
In conclusion, we have revealed here a global structural characteristic of MGs, namely that the atomic packing in MGs follows a regular local pattern, that of the nanometer-scale "gradient atomic packing structure". For amorphous structures,  The correlation coefficient c as a function of distance r for the asquenched Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 samples at different cooling rates. The correlation length is selected to be the distance when the value of correlation coefficient decays to 0.02 the nature of atomic packing in local regions is to stack liquid-like atoms first, then the transition atoms, and finally the solid-like atoms. Based on this notion, we propose a model for the amorphous structure of MGs as comprising solid-like regions, transition regions, and liquid-like regions. The discernible regions, including liquid-like atoms and their neighbors, tend to be soft and fertile locations for ST bands, corresponding to so-called "soft spots". In addition, each ST should percolate through liquid-like and transition regions, and finally be frustrated by the solid-like regions ("backbone") whose intensity is mainly determined by the degree of enrichment of solid-like atoms.
This notion of a "gradient atomic packing structure" provides a clear boundary to demarcate the structural components of the "soft spots" in MGs. Specifically, quantifying the fraction and spatial distribution of the liquid-like atoms, which constitute the prime nature of "soft spots", can indicate the propensity for strain localization in different amorphous alloys. We regard these conclusions as a step forward in establishing concrete relationships between nanoscale structure and macroscale properties in MGs, i.e., to defining the role of medium-range order in establishing the origin of plastic deformation in these materials. We also trust that this understanding can aid the structural underpinning of dynamic heterogeneity mapped out in experiments.

Generation of amorphous samples by MD simulation
Simulations were performed using the LAMMPS code 49 36 sample with 10,000 atoms was generated by cooling the liquid to 800 K with the same rate of 0.1 K/ps, and then hold the supercooled liquid for structural relaxation at 800 K for 400 ns (under an NPT ensemble) before quenching to 50 K. Six smaller as-quenched Cu 50 Zr 50 and Cu 47.5 Zr 47.5 Al 5 samples, each with 10,000 atoms, were quenched at 0.01, 0.1, and 1 K/ps, following the procedure stated above. Then the ensemble was switched to isothermal-isovolumetric ensemble (NVT); each sample was relaxed for 1 ns and 100 atomic configurations were collected for structure analysis. The time step employed in all the simulations was 2 fs. The temperature was controlled by the Nose-Hoover thermostat, 52 and the external pressure was held at zero by employing a Parinello Rahman barostat. 53 Periodic boundary conditions (PBCs) were applied in all three dimensions to eliminate any surface effects.
Calculation of the local atomic strain during shear deformation Pure shear loading was imposed on the as-quenched Cu 64 Zr 36 , Cu 50 Zr 50 , Cu 47.5 Zr 47.5 Al 5 , Ni 62 Nb 38 , and Mg 65 Cu 25 Y 10 samples and the annealed Cu 64 Zr 36 samples to different shear strains to induce atomic rearrangements with a rate of 10 7 s −1 . In order to weaken the effect of thermal fluctuation, the deformation temperature was set at 50 K, the sample responses mainly being derived from mechanical activation. To identify the atoms most likely involved in ST zones, the atomic strain was tracked by the local minimum non-affine displacement D 2 min À Á 39 and von Mises strain. 46,54 The D 2 min of central atom i relative to its nearest-neighbor atoms j can be defined as where N i is the nearest neighbors of the ith atoms determined by Voronoi analysis. r i ! ðtÞ is the position of atom ith at time t, and Δt is the time interval for the atomic rearrangement. To achieve the non-affine displacement parameter, we also employed the locally affine transformation matrix, J i , that best maps: d 0 ji n o ! d ji È É ; 8j 2 N i , where d 0 ji and d ji are bond vectors for referenced and current configurations between ith (central) and jth atoms, respectively.

Voronoi tessellation
Structural analysis was implemented using the Voronoi tessellation method to yield information about the coordination environment for atoms in MGs. 5,55 Employing this tessellation, the space is divided into close-packed polyhedra around each atom by constructing bisecting planes along the lines joining the central atom and all its neighbors. The Voronoi index is commonly denoted as <n 3 , n 4 , n 5 , n 6 > with n i as the number of i-edged polygons. For the multicomponent MGs, the bisection should be weighted by the atomic size to make up a Voronoi polyhedron.

Data availability
The data that support the findings of this study are available from the corresponding author, specifically Professor Yanqing Su of the Harbin Institute of Technology (email: suyq@hit.edu.cn) upon reasonable request.