Effect of surface and internal defects on the mechanical properties of metallic glasses

Despite the significance of surface effects on the deformation behaviours of small-scale metallic glasses, systematic investigations on surface states are lacking. In this work, by employing atomistic simulations, we characterise the distributions of local inhomogeneity near surfaces created by casting and cutting, along with internal distributions in pristine and irradiated bulk specimens, and investigate the effects of inhomogeneity on the mechanical properties. The cast surface shows enhanced yield strength and degrees of shear localisation, while the cut surface shows the opposite effects, although the fraction of vibrational soft spots, known to indicate low-energy barriers for local rearrangement, is high near both surfaces. Correspondingly, plastic deformation is initiated near the cut surface, but far from the cast surface. We reveal that improved local orientational symmetry promotes strengthening in cast surfaces and originates from the effectively lower quenching rate due to faster diffusion near the surface. However, a significant correlation among vibrational soft spots, local symmetries, and the degree of shear localisation is found for the pristine and irradiated bulk materials. Our findings reveal the sensitivity of the surface state to the surface preparation methods, and indicate that particular care must be taken when studying metallic glasses containing free surfaces.

vibrational modes is used to capture vibrational soft spots, which are fragile local clusters indicating potential STZ nucleation sites [16][17][18] . Irreversible plastic jumps are quantified by defining the atomic deformation gradient tensor 11,19 . Such plastic rearrangements accumulate until the local stress concentration field is strong enough for a cascade of STZs to evolve throughout the specimen. Recently, chemical heterogeneity has been used to explain the reaction barriers to the coalescence of STZs 20 . However, although surface effects have been shown to intensify as the specimen size is decreased in both computer simulations 21,22 and experiments 23,24 , the spatial distribution of structural disorder near the specimen surface and its effect on the mechanical properties of metallic glass have not been investigated in detail.
Here, by employing atomistic simulations, we systematically study the distribution of various structural defects on the surface and within the bulk, and the effects of defect distribution on the mechanical properties, by examining metallic glasses prepared by four different methods. We consider the pristine bulk, the specimen with cut surfaces, the specimen with cast surfaces, and irradiated bulk, referred to as 'pristine' , 'cut' , 'cast' , and 'irradiated' henceforth. For the clarification, we note that the 'bulk' in the present study represents nanoscale specimen with full PBCs. We carry out athermal mechanical simulations to characterise the inelastic transformation driven by external shear strain 11,25 , although realistic deformations in MGs are determined by the intricate coupling between thermal activation and external loading 8,12,26 . We find that cast surfaces show enhanced yield strength and degrees of shear localisation, while cut surfaces show the opposite effects, although the fractions of vibrational soft spots are increased near both surfaces. In addition, by analysing the spatial distribution of local stress fields, we find that the cascade of STZs, or shear localisation, is initiated near the surface for the cut at a smaller strain than in the pristine, while plastic deformation is initiated far from the surface in the cast. We reveal that the improved local orientation symmetry promotes the strengthening effect of the cast, and originates from the effectively lower quenching rate by faster diffusion near the surface. This implies that vibrational soft spots do not necessarily indicate potential STZ nucleation sites for specimens with free surfaces. On the contrary, a strong correlation among vibrational soft spots, local symmetries, and degrees of shear localisation is found for the pristine and irradiated. Our findings reveal the sensitivity of the surface state to the surface preparation method, and indicate that special care must be used when studying metallic glasses containing free surfaces. Figure 1a shows the stress-strain curves under simple shear tests of the four samples. The strength of the cut is smaller than that of the pristine, although the initial elastic responses are almost identical. The cut has a central microstructure identical to that of the pristine, but the surface region is expected to be softer because of the abrupt surface formation by the removal of the PBCs. The strength of the cast, however, is slightly larger than that of the pristine. The atoms near the surface have higher mobility, allowing more time to rearrange into lower-energy configurations, leading to a stiffening effect, as shown later. The irradiated shows significantly reduced stiffness and strength compared to the other samples, and exhibits more frequent activation steps (i.e. sudden discrete drops in stress) early in the shear test than the others. However, the flow stress at large strain is identical for all four specimen after the chains of STZs percolate the specimen, which is consistent with the literature 27 .

Results and Discussion
Previous studies have shown that the macroscopic strength is closely related to the order of the microscopic shear localisation 2 . Here, we represent the atomic strain tensor η i as J J I ( ) henceforth referred to as the atomic shear strain. The localised atomic shear strain of the samples at 20% applied shear strain (inset of Fig. 1a) shows that the cast exhibits a higher degree of shear localisation than the pristine does. However, the cut shows a much lower degree of shear localisation than the pristine, both centrally and near the free surface. This result implies the importance of an accurate consideration of the free-surface effects in predicting the deformation modes of metallic glasses. The surface states of the cut and cast significantly affect the degrees of shear localisation in the central planes, located ~27 Å away from the surfaces (Fig. 1c,d). In comparison, the irradiated has the least shear localisation. Our results show that the deformation mode of nanoscale metallic glass samples can be significantly affected by defects at both the exterior (surface treatment) and interior (irradiation processes) despite having identical chemical compositions and quenching rates. For a more detailed analysis of the plastic deformation, we investigate the spatial correlation function of atomic strain for all samples, as depicted in Fig. 2. The shear stress around local isolated inclusions with different moduli has four-fold symmetry. For example, the stress distribution of an elastic medium with a circular hole (i.e. void inclusion) with the radius R under the pure shear stress S is given as  29 . Hence, in the presence of multiple isolated inclusions considered as STZs, the spatial correlation for the atomic strain has four-fold symmetry. The spatial correlation function for the atomic strain tensor component η xy is defined as and represent the averages over all atoms in the simulation cell 30 . The formation of STZs can also be visualised using the correlation function for the non-affine displacement D min 2 (see Supplementary  Information). Figure 2a-d represent the spatial correlation of xy η for the pristine, cut, cast, and irradiated obtained for → = d r dx dy ( , , 0). From top to bottom, the correlation fields are shown at 5%, 10%, and 20% engineering shear strain. The correlation functions for the pristine and cast maintain clear four-fold symmetry up to 10% strain, as depicted in Fig. 2a and c, while the correlation fields of the non-affine displacement D min 2 are isotropic, as shown in Supplementary Figs 1-8. The four-fold symmetry is abruptly broken when the STZ cascade is initiated, and the correlation functions form a narrow band at 20% strain. On the contrary, as shown in Fig. 2b and d, the cut and irradiated have distorted four-fold symmetries along the dx direction at 5% strain, indicating the occurrence of some STZ cascades at this low strain. Interestingly, the shear localisation proceeds along the same axis for the irradiated, while the orientation of localisation takes the direction of the dy-axis in the cut. Such behaviour can be directly visualised from the apparently shear-transformed atoms (atoms with large plastic jumps, as defined formally later) in the Supplementary Figs 3-4. The STZ cascade typically proceeds in a consistent direction once initiated, because the stress concentration around an ellipsoidal inclusion (such as an STZ chain) is highest near the most prolonged tip of the inclusion. However, the cut sample has two distinct regions of the stiff centre and soft surface. While the initial plastic deformation and shear localisation occur near the surface at low applied stress, the STZ cascade cannot be maintained because the required stress for STZ formation is significantly higher in the central region. A new STZ cascade is nucleated throughout the central region at a higher stress; this may proceed along a different direction. The detailed evolution of the correlation function and the apparently shear-transformed atoms can be found in Supplementary Figs 1-8.
In order to investigate the structural origin of the observed plastic deformation, we examine the local stiffness values of the samples based on the vibrational density of states (V-DOS), as represented in Fig. 3. Sharper first peaks in the V-DOS generally correlate to the involvement of fewer low-frequency modes in the sample 17 . As expected from the strengths of the four samples, the first peak of the cast is the sharpest, as shown in Fig. 3a. However, the internal structural defects in the irradiated blunt the first peak, corresponding to a significant increase in the fraction of low-frequency modes. The increasing order of defects seems to have effects similar to those of the increasing quenching rate reported in the literature 17 the total number of atoms), e ( ) i λ → represents the eigenvectors from the Hessian matrix constructed from atomic , Ω is the set of selected elements whose vibration energy w  is smaller than the lowest 1% value of 3.71 meV in the pristine, and n is the cardinality of the set Ω 18 . We use the same cutoff value of w  for the four samples to accommodate the relative differences in vibrational soft spots among the samples. In order to consider the spatial distribution of vibrational soft spots along the aperiodic axes, we subdivide the samples into layers orthogonal to the surface, as depicted in Fig. 1b. The central plane of the layer direction (z-axis) is set to 0 and the distances to each edge are normalised to ±1. Figure 3b presents the variation of the participation fraction for the 1% lowest vibrational frequency modes through the layers. We find that the averaged participation ratio rises near the surface in for both the cast and cut. In addition, the cut shows a higher participation ratio near the surface than the cast does, originating from the abrupt surface creation. The pristine and irradiated show homogeneous distributions of soft spots, as expected. In the absence of a free surface, we can state that, because of the larger fraction of vibrational soft spots, STZ nucleation occurs more uniformly in the irradiated than in the pristine, eventually reducing the degree of shear localisation. However, although the vibrational soft spot fraction is higher at the casted surface, plastic deformation is not initiated near the surface and the degree of shear localisation is higher than that in the pristine. This implies that the vibrational soft spots do not necessarily indicate potential STZ nucleation sites for specimens with free surfaces. 3D views of the participation fraction distributions of the samples are illustrated in Fig. 3c-f.
Next, we examined the spatial distribution of local orientation symmetry to investigate the origins of different mechanical properties in the four specimens. The Voronoi tessellation, a typical method for determining the Kasper polyhedra, cannot be clearly defined for aperiodic boundaries, because of ambiguities in defining the surface encompassing the total volume. Instead, we adopted the BOO parameter Q 6 , which can be defined only with neighbouring atomic arrangements 14 . We first examined the correlation between the Cu-centred Q 6 and various Cu-centred locally favoured motifs as well as geometrically unfavoured motifs (GUMs) known for the Cu 64 Zr 36 system for the pristine (Fig. 4a) 6 . Because the Cu-centred icosahedron is known as the most stable local cluster in the selected chemical composition, we attempted to capture the local structures from the BOO parameter Q 6 = 0.663, corresponding to perfect icosahedral symmetry 15 . We define the range of BOO parameters to include slightly distorted Cu-centred icosahedra and test the correlation between the icosahedral order found from both Voronoi tessellation and the Q 6 parameter, as shown in Fig. 4a-c and Supplementary Fig. 9a,b. The range of Cu-centred Q 6 is determined to include the same amount of atoms as the total number of Cu-centred full icosahedra from the Voronoi tessellation in the pristine. In the range considered (Q 057174 6 > . ), the atoms selected by Cu-centred icosahedra from Voronoi tessellation and the atoms in the selected Q 6 range match at 85% of sites, as shown in Fig. 4c. We also find that Cu atoms with higher Q 6 parameters, illustrated by white circles, are located at the low-participation-fraction sites, i.e. sites with higher local structural rigidity, in the pristine, as shown in Fig. 4b and in all samples in Supplementary Fig. 10. The correlation between Zr-centred Q 6 and the participation fraction is depicted in Supplementary Fig. 9a, where no strong correlation is seen. Having established the validity of using the Q 6 parameter to identify local icosahedral order, we calculate the fraction of selected Q 6 parameters found along the layer of the specimen containing free surfaces (Fig. 4d). We find that more stable microstructures, i.e. higher local symmetries, are developed near the cast surface, except in the outermost layers. In the cast, the favourable SROs developed near the surface reduce the surface nucleation of STZs, as shown in Supplementary Figs 5 and 6.
In order to understand the origin of enhanced symmetry near the cast surface, we calculated the self-diffusivity constant D for the cast at 700 K, below the glass transition temperature (~745 K) of the selected chemical composition 2 . The self-diffusivity D is defined as where the numerator is the definition of the mean square displacement (MSD), 〈 〉 symbolises the average over all atoms in the cast, and the initial time is t 0 31 . Figure 5a presents the MSDs during the 200 ns of NVT simulations at 700 K obtained at the outermost and central layers of the cast. The atoms near the surface travel significantly higher diffusion distances than those near the central region. D along the aperiodic layer is then calculated, using the last 20 ns of the simulation time, as shown in Fig. 5b. Interestingly, D in the outermost layer is ~5 times larger than that in the central layer. This implies that the atoms near the surface remain mobile even below the glass transition temperature. Hence, the atoms near the surface experience much longer effective structural relaxation times, and form higher-symmetry structures after the quenching process. We note that, when the cut sample is annealed at 700 K for 200 ns, structural properties and mechanical responses become almost identical to the cast sample, which also proves that the facile atomic diffusion at surfaces is responsible for the higher structural symmetry and stiffer mechanical response of the cast surface as shown in Supplementary Fig. 11. Because such relaxation process occurs much slower at 300 K, we observe significantly less changes on structural and mechanical properties after 300 K annealing. Figure 6 presents the correlation between the event of irreversible atomic rearrangement and the soft spots, based on the participation fraction analysis. Supplementary Fig. 12a represents the frequency of the non-affine displacement accumulated for up to 10% engineering shear strain. The non-affine displacement is defined from J i as ( ) where N i is the number of nearest neighbours 5 and the reference configuration is taken from the previous strain step. We define the threshold of non-affine displacement at the value of 0.1 Å 2 , shown as a dotted purple line, as indicated in Supplementary Fig. 12. Figure 6a represents the distribution of the number of plastic events along the aperiodic axis. In the cut, most highly vibrational soft spots, coloured red, undergo the plastic event, while a  smaller fraction of highly vibrational soft spots does so in the cast. While the plastic events are first localised near the surfaces before proceeding to the central region for the cut, the cast exhibits dominant plastic events far from the surface because of the lower structural symmetry (lower BOO) in the central region (see Supplementary Figs 3-6 for details). As depicted in Fig. 6b, the cut exhibits frequent surface-driven plastic rearrangement early in strain application. The irradiated, on the other hand, exhibits uniform plastic jump distribution along all layers. The time evolutions of the plastic jump distributions for the four samples are illustrated in Supplementary Fig. 13. We note that the width of the apparently shear-transformed atoms near the cut surface remains almost the same at ~15 Å, until the STZ cascade completely percolates the central region. In other words, the softening from the abrupt surface formation extends ~15 Å from the surface. Hence, for metallic glasses in atomistic simulation with sizes comparable to the softening depth, the degree of shear localisation would be affected significantly.
We expect a similar effect by surface softening for the metallic glass samples tested in experiments. For nanopillars experimentally prepared by the focused ion beam (FIB) method, the depth of softening would be significantly higher, which explains the enhanced plasticity of FIB sample at 100 nm scale 32 . In contrast, 300 nm or larger diameter structures do not show significant changes in deformation behaviour under compression or tension regardless of FIB damage, suggesting that the FIB damage affects the local symmetry in the region much shallower than 300 nm 23,32,33 . In contrast, significant plasticity at 400 nm scale was achieved when the interior of the specimen is also damaged via a series of proton irradiation processes 10 . Our annealing simulations on the cut sample ( Supplementary Fig. 11) imply that the surface damage by FIB may be recovered after long time relaxation process at high temperature, and the critical size for the intermittent-to-homogeneous deformation 23,33 can be increased.

Summary and Conclusion
We have investigated defect distributions and the effects thereof on the mechanical properties for pristine, cut, cast, and irradiated Cu 64 Zr 36 metallic glasses. We find that vibrational soft spots are concentrated near the surfaces, regardless of the preparation methods, because of broken bonds at the surfaces. Meanwhile, local orientational order is enhanced only near the cast surface, because of the longer effective relaxation time originating from the higher atomic mobility near the surface. Hence, the rigidity and the degree of shear localisation are enhanced for the cast compared to the pristine, while they are diminished for the cut. This implies that vibrational soft spots do not necessarily indicate low orientational order or preferential STZ nucleation sites in the presence of free surfaces. The irradiated experiences homogeneous damage throughout the volume; the reduced rigidity and lower degree of shear localisation can be understood from the correlation among vibrational soft spots, low orientational order, and preferential STZ nucleation sites. Our findings imply that metallic glass specimens with sizes comparable to the surface damage depth (15 Å for the cut sample in the present study) experience lower degrees of shear localisation and higher ductility. In real experiments where surface damage from FIB or other cutting methods is deeper than that in the present study, abrupt changes in the deformation modes, such as brittle-to-ductile transitions, may occur at larger scales. In the future work, by employing the advanced sampling techniques such as autonomous basin climbing (ABC) methods for accessing experimental strain rate, we plan to investigate the rate-dependent intermittent-to-homogeneous deformation transition observed in experiments 8,34 .

Methods
We studied the Cu 64 Zr 36 metallic glass because it has a relatively high density of five-fold locally favoured motifs called icosahedra 2 . The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) was employed to perform molecular statics/dynamics simulations 35 . The embedded-atom-method (EAM) interatomic potential, developed by Sheng et al., was selected to describe interatomic interactions in the CuZr binary system 2 . In order to successfully examine the surface effects, we used high surface-area-to-volume-ratio cubic specimens with dimensions of ~54 Å × 54 Å × 54 Å containing 10,000 atoms. Each sample in the study was equilibrated for 10 ns at 2000 K and quenched to 300 K at the rate of 10 9 K/s. We chose the quenching rate of 10 9 K/s as a computationally limiting case to have well-relaxed glass structure 36 . All stress components along the periodic boundaries were relaxed during the quenching processes.
The four types of metallic glass samples of the pristine, cut, cast, and irradiated are depicted in Fig. 1. Four samples represent the pristine, with PBCs along all three axes; two samples (cut and cast) have free boundaries along the z-axis and PBCs along the x-and y-axes; another pristine bulk sample incorporates irradiation damage. The cut is prepared by abruptly turning off the PBC along the z-axis from the pristine after quenching. The cast is produced by cutting the sample immediately after the 2000 K equilibration process and quenching it in the presence of the external harmonic spring potential wall represented by k(z−r c ) 2 along the aperiodic direction where z is the distance from the central x-y plane of the specimen. The potential wall is activated when z is larger than the cutoff distance r c . The effective spring constant k is set to 0 5eV/Å 2 .
and r c is defined as one-half the vertical length of the sample at ~28 Å to confine the atoms during the quenching period. The irradiated is constructed by applying a series of atomic bombardments to the pristine bulk. For a bombardment event, a velocity corresponding to 100 eV in energy is assigned to a randomly chosen atom in a random direction. Then, the sample is equilibrated for 1.5 ps with the microcanonical ensemble (NVE) and 3.0 ps with the isothermal-isobaric ensemble (NPT) sequentially. This bombardment process is repeated 200 times to reproduce a certain level of structural transitions 10 . The short-range pairwise Lennard-Jones potential (ε = . 2 0 LJ σ = . Å, 1871 LJ eV) is augmented during the irradiation process to resolve the soft core problem of the EAM potential. The quenched samples at 300 K are then relaxed using a conjugate gradient. Athermal quasi-static simple shear tests (AQS) are performed in the x-y plane by gradually changing the shape of the simulation cell with incremental steps of 0.001 of the engineering strain ( ENG xy  ).
Scientific RepoRts | 7: 13472 | DOI:10.1038/s41598-017-13410-3 The structural order was characterised by various measures. We obtained the entire normal mode spectrum from the full Hessian matrix and identified soft vibrational modes from the participation ratio of each atom at the vibration energy corresponding to the lowest 1% normal-mode frequency of the pristine 17 . With Voronoi tessellation, we obtained the Kasper polyhedron for each atom, from which we could identify local icosahedral structures in the pristine. Because of the ambiguity of defining boundaries by Voronoi tessellation, the local icosahedral order for the cut and cast were analysed by the Q 6 parameter from the bond orientational order (BOO) analysis 14,15 , after establishing its correlation with the Kasper polyhedron in the pristine. The local deformation gradient tensor J i was used to define the atomic strain tensor η i and the non-affine displacement D min 2 5 . J i mapped the relative displacement against the reference, i.e. { } { } J d d : where N i is the number of nearest neighbours of the i th atom. The cutoff distance for assigning N i was chosen as the first minimum position of the radial distribution function of the pristine bulk sample of 3.5 Å. In this study, the reference state for the atomic strain tensor was fixed at the initial zero-strain state, while the reference for non-affine displacement was the configuration at the previous strain step. We defined the apparently shear-transformed atoms as those with accumulated D min 2 values over the threshold of 0.1 Å 2 until reaching 10% of the applied (engineering) shear strain. These atoms were likely within STZs. The total number of irreversible jumps was captured by counting the number of apparently shear-transformed atoms 37 .