Correlations of non-affine displacements in metallic glasses through the yield transition

We study correlations of non-affine displacement during simple shear deformation of Cu-Zr bulk metallic glasses in molecular dynamics calculations. In the elastic regime, our calculations show exponential correlation with a decay length that we interpret as the size of a shear transformation zone in the elastic regime. This correlation length becomes system-size dependent beyond the yield transition as our calculation develops a shear band, indicative of a diverging length scale. We interpret these observations in the context of a recent proposition of yield as a first-order phase transition.


Introduction
Bulk metallic glasses (BMG), multi-component metals that are kinetically arrested into an amorphous structure, have been suggested for wide range of applications, including as structural materials 1,2 . For practical applications a big problem is their tendency to form shear bands, planar regions in that localize most of the plastic deformation at relatively low strain. These shear bands are the primary mechanism by which metallic glasses fail. Numerous ideas to address this problem have been suggested, such as deliberately introducing pores 3 or creating "nanoglasses" that have an internal microstructure 4 .
The deformation of BMGs is described by the theory of shear transformations or shear transformation zones (STZs) [5][6][7] , localized rearrangements of small regions of atoms. The size of these zones has been estimated to range from a few 8 to many tens of atoms 9,10 . Knowledge of the size of the zones could help to fundamentally understand this class of materials on an atomic level and be used in mesoscale simulations that incorporate STZs [11][12][13] . The size of STZs has been linked to the Poisson ratio 14 as well as the brittle or ductile character of fracture of BMGs [15][16][17] .
Spatial correlation functions of non-affine deformation have recently been employed to quantify the geometry of STZs. Murali et al. 18 looked at the spatial autocorrelation in the non-affine deformation field of deformed bulk metallic glasses in molecular dynamics simulations. They found an exponential decay of the autocorrelation from which they extracted a correlation length ℓ, which they interpreted as the size of an STZ. These findings have been confirmed by similar calculation on Lennard-Glasses 19 . In a similar spirit, Chikkadi et al. 20,21 have discussed the autocorrelation of non-affine deformation in experiments of sheared colloidal glasses. In addition to the global non-affine displacement field, they characterized the local nonaffine deformation through the #$% & measure of Falk & Langer 7 . Their data shows long-range correlations as manifested in a power-law behavior of the autocorrelation function global and local measures for non-affinity. In contrast to Murali et al.'s data 18 , this suggests a scale-free character of the deformation. Calculations of hard-sphere mixtures carried out for the interpretation of these experiments did again yield an exponential decay of the correlation function 21,22 . Varnik et al. 23 argued that this is because of limitation in system size; larger calculation, albeit carried for a 2D soft disk model rather than in 3D, indeed showed power-law correlations. Earlier calculations of the correlation of the vorticity field during deformation of a 2D Lennard-Jones solid showed similar power-law correlations 24,25 .
We here revisit the question of exponential vs. power-law correlations and provide new data on how they evolve through the yield transition. Our molecular dynamics calculations of the deformation of BMGs show the emergence of correlations in the non-affine part of the deformation field of calculation larger than those previously reported. This allows us to extract the correlation function up to distance ~75 times the nearest neighbor distance for the largest systems studied here, similar to previous 2D calculations that showed power-law correlation 23 . While we do find exponential and not power-law correlations, the length-scale ℓ associated with the exponential becomes a function of system size after shear-band nucleation, indicating a divergent length at the nucleation of the band.

Molecular dynamics simulations
We conducted all simulations using Molecular Dynamics (MD) and the second generation of the interatomic Cu-Zr potential by Mendelev et al. 26 . Amorphous sample systems were first obtained by melting and equilibrating systems of Cu50Zr50 (or other stoichiometries, see below) at 2500 K for 100 ps, followed by a linear quench to 750 K at a rate of 6 K ps -1 . This temperature is slightly above the glass transition temperature ( ≈ 600 K, as obtained from the jump in heat capacity when cooling the system through Tg at the same rate. We then aged the system for 1 ns before quenching it to 0 K at a rate of 6 K ps -1 . We used a Berendsen barostat 27 with a relaxation time constant of 10 ps to keep the hydrostatic pressure in the simulation cell at zero and Langevin thermostat 28 with a relaxation time constant of 1 ps to control temperature during quench and equilibration.
To prepare simulations carried out at different temperatures, the amorphous systems were then equilibrated at zero pressure for 200 ps at different temperatures between 0 and 300 K. The cell was subsequently deformed using simple shear deformation at constant volume at an applied shear rate of ̇= 10 2 s 45 up to a maximum of 35% strain. To control temperature, we again used a Langevin thermostat but only thermalize the Cartesian direction normal to the plane of shear to eliminate any drag with respect to some reference velocity field. The bulk of our simulations comprises a cubic cell with an edge length of ≈ 206 Å and 500,000 atoms. The potential influence of finite-size effects was studied using two additional system sizes: A cubic system with twice the edge length and eight times the number of atoms and another cubic system with half the edge length and 1/8 the number of atoms. If not mentioned otherwise, results are reported for the = 206 Å system.

Local strain measure and correlation
To quantify heterogeneous flow of the system, we need measures that identify local deformation events. Falk & Langer 7 introduced a method to determine the local deformation of an atomic system within spheres of radius rcut. The idea is to map for each atom i its atomic neighborhood at time t-∆t to the neighborhood at time t using an affine deformation with deformation gradient : ; , and then find : ; that minimizes the residual error. The final residual error, . Note that for our calculations carried out at constant applied shear rate , we specify the reference frame by its distance in applied strain Δ rather than Δ ≡ Δ /.
To quantify the geometry of these deformation events, we use spatial auto-correlation maps. The auto-correlation map of some field ( ⃗) is defined as Note that Eq. (5) is normalized such that, because of Eq.  perpendicular to the plane where shear is applied (Fig. 2b), shows that the system was supercooled and atoms broke out of their cages at a strain between 0.1% and 1% (see e.g. Refs. 31,32). The shear stress ‚ƒ in the plane of shear (Fig. 3a) initially rose linearly with the strain applied in the xy-plane. At around ≈ 10% the system yielded and the stress ‚ƒ dropped from a peak value to a plateau region where the ‚ƒ remained constant up to an applied strain of = 35%, the maximum strain applied in our calculations. Our five calculations at 0 K, 50 K, 100 K, 200 K and 300 K show that the system softened as temperature increased; from a yield stress of around 1.7 GPa in the athermal limit to 1.2 GPa at 100 K. At small strain where ‚ƒ ( ) is linear, we find localized events (Fig. 3b).

Results
After yield, these localized events coalesce to shear-bands, first vertical (Fig. 3c, see also Ref. 19) but later horizontal (Fig. 3d), developing a clear anisotropic structure. Note that such vertical shearbands occurred only in some of our calculations. At small strain, both shear-band directions are equivalent and the nucleation direction is random. Symmetry breaking at larger strain forces the shear band back into the direction of shear.
After yield (Fig. 4b), Š ( , ) develops a clear anisotropic structure with a band of large correlation parallel to the x-axis. Fig. 4c shows radial averages ( ) of the data of Figs. 4a and b. There are oscillations at small distances that turn into an exponential decay at around 10 Å. Oscillations at small distances are due to the structure of the amorphous solid. We therefore normalize the autocorrelation function and define to remove variations in ( ) due to variations in local atomic density. Figure 4d shows that these oscillations are eliminated in ̅ ( ) for r > 5Å.
In the following, we characterize the exponential decay, by fitting the correlation length ℓ in Eq. (7) over a select section of the correlation function. We distinguish between the behavior at short distances 5 Å < r < 15 Å, which is within the range of  (Fig. 5c). The behavior of ℓ "'%" is different. Its value (Fig. 5d) is independent of QRS used in the computation of #$% & and the data does not collapse when normalized accordingly. The evolution of ℓ ••''S and ℓ "'%" with applied strain clearly show the point where the samples yield (cf. also Fig. 3a). At around 12.5% strain, ℓ "'%" increases dramatically. During further deformation it fluctuates around a value consistently a factor of 3 higher than before yield.
The computation of ̅ ( ) furthermore depends on how the reference frame for the computation of #$% & is chosen. Here, we report results obtained for references frames at constant distance in applied strain, Δ . All results reported above were obtained for Δ = 1%. Fig. 6 demonstrates how the correlation function and ℓ vary as a function of this parameter. Before yield, the correlation function does not depend on Δ and drops exponentially over two decades as a function of distance. Fig. 6a shows this behavior for Δ = 1%, 0.1% and 0.01% which is above, at and below the cagebreaking strain (Fig. 2b). The behavior changes at yield (Fig. 6b), where the initial exponential drop starts to depend on Δ . Fig. 6c shows the influence on the extracted value of ℓ "'%" . ℓ "'%" decreases with decreasing Δ and saturates at ℓ "'%" ≈ 15 Å in the flow region for the lowest Δ = 0.01%.
To clarify the role of Δ on the calculation of the correlation length ℓ, we further test the influence on system size on the correlation functions. Fig. 7a shows that before yield (applied strain = 7%), ̅ ( ) is independent of system size but that a clear size dependence develops when the material flows (Fig. 7c, = 20%). Plotting ℓ "'%" versus applied strain shows that the before yield ℓ "'%" is independent of size but after yield it depends on system size (Fig. 7b). Normalizing distance r or correlation length ℓ "'%" by system size collapses all data in the region where the amorphous solid flows (Fig. 7c and d).
Finally, we test the dependence of the correlation function on temperature and composition. Fig.   8a shows the temperature dependence of ℓ "'%" . Data in the temperature range from 50K to 300K, below the glass transition temperature of ( ≈ 600 of our metallic glass, is superimposed for small strain. It appears that at large strain the higher temperatures lead to a smaller ℓ "'%" but out present data is too noisy to make a firm conclusion. Fig 8b shows ℓ "'%" ( ) for different compositions. Again the data collapses in the elastic regime and there appears to be a slight variation with composition after the sample has yielded.

Discussion
The correlation length ℓ characterizing the exponential decay of the spatial-autocorrelation Δ that determines over how many STZs we average does not affect the results. The situation changes dramatically after the sample has yielded ( > 10%). STZs are now localized within a shear band and it becomes difficult to identify individual STZs ( Fig. 3c and d). The onset of shearbanding is then accompanied by a characteristic length ℓ "'%" proportional to the system size and that depends on strain offset Δ . For small Δ , we find values for ℓ comparable to the ones found in the elastic regime ℓ (Fig. 6c). We hypothesize, that this is because even for the flowing glass we can pick out individual STZs if we look at small enough strain increments, much smaller than the cage-breaking strain (Fig. 2b).
We note that while in the elastic regime our correlation functions look clearly exponential, our system sizes albeit large are yet too small to rule out power-law behavior during flow. Indeed, the fact that our length scale ℓ "'%" depends on system size is indicative of a diverging length or a crossover to a power-law as STZ events become correlated within the shear band. This observation is consistent with a recently proposition that yield in amorphous solids can be interpreted as a first-order phase transitions 34,35 , an interpretation that has a rich history for explaining shear-banding instabilities in non-Newtonian fluids 36 . Jaiswal et al. 34 identify the transition using an order parameter that measures similarity or "overlap" of atomic configuration. The atomic configuration uses overlap with the initial configuration at yield. A central observation is that their "yield" point occurs at larger strains than the overshoot in the stress-strain curve that is typically attributed to yield. This is consistent with our calculations, which show that ℓ rises after the stress has peaked (cf. Fig. 3a and 7b,d).

Summary & Conclusion
We studied the correlation between nonaffine displacements, as characterized by the #$% & measure of Falk & Langer 7 , using molecular dynamics calculations. This multipoint correlation function shows exponential behavior in the elastic regime from which we can extract a length scale ℓ, typically attributed to the size of an STZ. We find that this length scale diverges at yield, as manifested by a size-dependent ℓ in during flow of the material. The diverges of ℓ occurs at strains larger than the peak stress that is typically attributed to the yield point. Our results support a recent proposition that yield in amorphous materials can be interpreted as a first-order phase transition 34,35 .