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

Bulk metallic glasses (BMGs), 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 [3], planar regions in that localize most of the plastic deformation at relatively low strain [4, 5]. These shear bands are the primary mechanism by which metallic glasses fail in tensile [6] or cyclic (fatigue) loading [7, 8] or during indentation [9]. Numerous ideas to address this problem have been suggested, such as deliberately introducing heterogeneities, such as pores [10], nanocrystals [11] or internal interfaces [12].

The deformation of BMGs is described by the theory of shear transformations or shear transformation zones (STZs) [13–15], localized rearrangements of small regions of atoms. The size of these zones has been estimated to range from a few [16] to many tens of atoms [17, 18]. 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 [19–21]. The size of STZs has been linked to the Poisson ratio [22] as well as the brittle or ductile character of fracture of BMGs [23–25].

Spatial correlation functions of non-affine deformation have recently been employed to quantify the geometry of STZs. Murali et al [26] looked at the spatial autocorrelation in the non-affine deformation field of deformed BMGs in molecular dynamics (MD) simulations. They found an exponential decay of the autocorrelation from which they extracted a correlation length ${\ell },$ which they interpreted as the size of an STZ. These findings have been confirmed by similar calculations on Lennard-Jones-Glasses [27]. In a similar spirit, Chikkadi et al [28, 29] 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 non-affine deformation through the ${D}_{{\rm{\min }}}^{2}$ measure of Falk and Langer [15]. Their data shows long-range correlations as manifested in a power-law behavior of the autocorrelation function in both global and local measures for non-affinity. In contrast to Murali et al's data [26], 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 [29, 30]. Varnik et al [31] 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. In a large study comparing different simulation models and experiments, Cubuk et al [32] appear to find again only exponential correlation in ${D}_{{\rm{\min }}}^{2}.$ Note that correlation of the vorticity (strain) field during deformation of a 2D Lennard-Jones solid showed similar power-law correlations [33, 34], and we have recently found power-law correlation in the displacement-difference field [35]. The question of power-law versus exponential correlation clearly hinges on the details of the quantity whose correlations are studied.

We here revisit the question of exponential versus power-law correlations in ${D}_{{\rm{\min }}}^{2},$ provide new data on how they evolve through the yield transition and show that the question of exponential versus power-law correlation depends on how ${D}_{{\rm{\min }}}^{2}$ is calculated. Our MD calculations of the deformation of BMGs show the emergence of correlations in the non-affine part of the deformation field of calculations 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 [31]. While we do find exponential and not power-law correlations, the length-scale ${\ell }$ associated with the exponential becomes a function of system size after shear-band nucleation, indicating a divergent length at the nucleation of the band, but only if the 'strain window' used to calculate ${D}_{{\rm{\min }}}^{2}$ is chosen large enough.

2.1. MD simulations

2.2. Local strain measure and correlation

To quantify heterogeneous flow of the system, we need measures that identify local deformation events. Falk and Langer [15] introduced a method to determine the local deformation of an atomic system within spheres of radius r_cut. 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 ${\bar{F}}_{i},$ and then find ${\bar{F}}_{i}$ that minimizes the residual error. The final residual error

$$\begin{eqnarray}&&{D}_{{\rm{\min }},i}^{2}={\min }_{{\bar{F}}_{i}}\left\{\displaystyle \frac{1}{N}\displaystyle \sum _{k}^{N}{\left[{{\boldsymbol{r}}}_{ik}\left(t\right)-{\bar{F}}_{i}\cdot {{\boldsymbol{r}}}_{ik}\left(t-{\rm{\Delta }}t\right)\right]}^{2}\theta \left({r}_{{\rm{cut}}}-{r}_{ik}\right)\right\},\end{eqnarray} \tag{ 1 }$$

is a measure for the non-affine component of local deformation. Here ${\rm{\Theta }}(r)$ is the Heaviside step function. STZs and shear bands can be identified by looking for regions with high values of ${D}_{\min ,i}^{2}.$ Note that for our calculations carried out at constant applied shear rate $\dot{\varepsilon },$ we specify the reference frame by its distance in applied strain ${\rm{\Delta }}\varepsilon ,$ the strain window, rather than ${\rm{\Delta }}t\equiv {\rm{\Delta }}\varepsilon /\dot{\varepsilon }.$

To quantify the geometry of these deformation events, we use spatial auto-correlation maps. The auto-correlation map of some field $Q(\vec{r})$ is defined as

$$\begin{eqnarray}&&{\mathscr{A}}\left[Q\right](\vec{r})=V\displaystyle \int {d}^{3}r^{\prime} Q\left(\vec{r}^{\prime} \right){Q}^{* }\left(\vec{r}-{\vec{r}}^{\mbox{'}}\right)=\displaystyle \frac{V}{{N}_{p}^{2}}\displaystyle \sum _{i}\displaystyle \sum _{j}{Q}_{i}{Q}_{j}^{* }\,\delta \left(\vec{r}-{\vec{r}}_{ij}\right),\end{eqnarray} \tag{ 2 }$$

where the last identity is the expression obtained for ${N}_{p}$ point particles for which $Q\left(\vec{r}\right)=\displaystyle \sum _{i}{Q}_{i}\delta \left(\vec{r}-{\vec{r}}_{i}\right)$ where ${Q}_{i}$ is the value of quantity $Q$ on atom $i.$ Note that for any quantity $Q$ this autocorrelation map obeys the following sum rules

$$\begin{eqnarray}&&{\mathscr{A}}\left[Q\right]\left(\vec{r}\to 0\right)=\left\langle {Q}_{i}^{2}\right\rangle \,{\rm{and}}\,{\mathscr{A}}\left[Q\right]\left(\vec{r}\to \infty \right)={{Q}_{i}}^{2},\end{eqnarray} \tag{ 3 }$$

where the average $\left\langle {Q}_{i}\right\rangle =\displaystyle \sum _{i}{Q}_{i}/{N}_{p}.$

The radial average of this auto-correlation map gives a function ${\mathscr{A}}[Q](r),$ depending only on the distance and not the direction between the two atoms. The auto-correlation function of unity is the radial distribution (or pair correlation) function

$$\begin{eqnarray}&&{g}_{2}\left(r\right)={\mathscr{A}}\left[1\right]\left(r\right).\end{eqnarray} \tag{ 4 }$$

By virtue of equation (3) it is normalized such that ${g}_{2}\left(r\right)\to 1$ as $r\to \infty .$ We are specifically interested in the correlations of ${D}_{{\rm{\min }}}^{2}$

$$\begin{eqnarray}&&C\left(\vec{r}\right)=\displaystyle \frac{{\mathscr{A}}\left[{D}_{\min }^{2}\right]\left(\vec{r}\right)-{\left\langle {D}_{\min ,i}^{2}\right\rangle }^{2}}{\left\langle {\left({D}_{\min ,i}^{2}\right)}^{2}\right\rangle -{\left\langle {D}_{\min ,i}^{2}\right\rangle }^{2}}.\end{eqnarray} \tag{ 5 }$$

Note that equation (5) is normalized such that, because of equation (3), $C\left(\vec{r}\to 0\right)=1$ and $C\left(\vec{r}\to \infty \right)=0.$ We compute ${\mathscr{A}}\left[Q\right](\vec{r})$ at short distances by direct evaluation of equation (2) and at long distances using a fast Fourier transform to speed up the convolution in equation (2), allowing us to efficiently compute the correlation function up to half the size of our systems. Details of this algorithm are given in the appendix.

Figure 1(a) shows a snapshot of the quenched system before shear. The radial distribution function ${g}_{2}\left(r\right)$ (figure 2(a)) is indicative of a disordered structure with broad nearest and second-nearest neighbor peaks and non-zero probability for finding a neighbor between them. The first neighbor peak is located at ${r}_{{\rm{NN}}}=2.8\,\mathring{\rm{A}}$ and indicated by a vertical dashed line and we use r_NN to normalize all distances reported below. The value of the non-affine displacement ${D}_{\min ,i}^{2}$ depends on the cutoff distance ${r}_{{\rm{cut}}}$ for identifying neighbors of an atom and on the distance ${\rm{\Delta }}\varepsilon$ between current and reference configuration in the time domain. In the following, we will show results obtained for ${r}_{{\rm{cut}}}$ being an integer multiple of the nearest-neighbor distance ${r}_{{\rm{NN}}}$ as given by the position of the first peak in ${g}_{2}\left(r\right).$ These distances are indicated by the vertical dashed lines in figure 2(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After equilibrating, we deformed our metallic glass under simple shear. Figures 1(b) and (c) show exemplary snapshots of these calculations. As an important dynamical quantity, we analyzed the mean square displacement (MSD) as a function of applied strain $\varepsilon ,$ where we only considered the component of the displacement in the z-direction, perpendicular to the plane where shear is applied (figure 2(b)). This is necessary in order not to contaminate the displacement by the applied shear; we could have subtracted the streaming velocity alternatively. Note that at constant strain rate $\dot{\varepsilon },$ the applied strain $\varepsilon =\dot{\varepsilon }t$ is the same as the usual parameter time t. The MSD shows the typical behavior for a glass [39]: initially, ${\rm{MSD}}\propto {\varepsilon }^{2},$ showing the ballistic motion of each atom within its local environment. The MSD then saturates at a value ≈0.01 r_NN, corresponding to a distance of around 0.1 r_NN² much smaller than the nearest neighbor distance (see figure 2(a)) but around the Lindemann criterion for melting [40]. Atoms do not diffusive in this regime but are trapped within their local cages. Finally, at large $\varepsilon ,$ the system entered a diffusive regime where ${\rm{MSD}}\propto \varepsilon .$ This diffusive regime is accompanied by a breaking out of the individual cage, because the mean distance traveled by the atoms now exceeds the nearest neighbor distance. Without mechanical agitation, this breaking out of the cage happens after what is known as the $\beta$-relaxation time. Under mechanical agitation, it happens at the cage-breaking strain ${\varepsilon }_{c},$ here between ${\varepsilon }_{c}\approx 0.1 \%$ and 1% (see also [41, 42]).

The shear stress ${\sigma }_{xy}$ in the plane of shear (figure 3(a)) initially rose linearly with the strain $\varepsilon$ applied in the xy-plane. At around $\varepsilon \approx 10 \%$ the system yielded and the stress ${\sigma }_{xy}$ dropped from a peak value to a plateau region where ${\sigma }_{xy}$ remained constant up to an applied strain of $\varepsilon =35 \% ,$ the maximum strain applied in our calculations. Our five calculations at 0, 50, 100, 200 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 0.8 GPa at 300 K.

Zoom In Zoom Out Reset image size

Figures 3(b)–(d) show a map of ${D}_{{\rm{\min }}}^{2}$ during deformation, here computed for a cutoff ${r}_{{\rm{cut}}}=3\,{r}_{{\rm{NN}}}$ and a strain window of ${\rm{\Delta }}\varepsilon =1 \% ,$ about the cage-breaking strain ${\varepsilon }_{c},$ before the frame shown in the figure. At small strain $\varepsilon$ where ${\sigma }_{xy}\left(\varepsilon \right)$ is linear, we find localized events (figure 3(b)). After yield, these localized events coalesce to shear-bands, first vertical (figure 3(c), see also [27, 43–45]) but later horizontal (figure 3(d)), developing a clear anisotropic structure. Note that such vertical shear-bands 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. We note that the shear-bands are only visible for strain windows ${\rm{\Delta }}\varepsilon \geqslant {\varepsilon }_{c}.$ Figures 3(e) and (f) show maps of ${D}_{{\rm{\min }}}^{2}$ for ${\rm{\Delta }}\varepsilon =0.1 \%$ and ${\rm{\Delta }}\varepsilon =0.01 \% ,$ respectively. Observed at these small strain windows, the shear band dissolves into individual disconnected localized events.

To statistically quantify this (random) structure we computed the ${D}_{{\rm{\min }}}^{2}$ auto-correlation maps, $C\left(\vec{r}\right).$ Figures 4(a) and (b) show a slice $C\left(x,y,z=0\right)\equiv {C}_{0}(x,y).$ Before yield (figure 4(a)), ${C}_{0}\left(x,y\right)$ shows a rotationally symmetric structure with a visible ring at the nearest-neighbor distance ${r}_{{\rm{NN}}}.$ After yield (figure 4(b)), ${C}_{0}\left(x,y\right)$ develops an anisotropic structure with a band of correlation parallel to the x-axis. Figure 4(c) shows radial averages $C(r)$ of the data of figures 4(a) 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

$$\begin{eqnarray}&&\bar{C}\left(r\right)=C\left(r\right)/{g}_{2}\left(r\right)\end{eqnarray} \tag{ 6 }$$

to remove variations in $C\left(r\right)$ due to variations in local atomic density. The normalized correlations are also shown in figure 4(c). The oscillations are eliminated in $\bar{C}\left(r\right)$ for r > 5 Å. Note that beyond yield, the radial symmetry of $C\left(r\right)$ is lost (see figure 4(b)). Figure 4(d) shows $C\left(r\right)$ parallel and perpendicular to the shear band. The positions where the data is taken from is marked with the black and red line in figure 4(b). This shows that the correlation is practically constant in the shear band direction, while it decays perpendicular to it.

Zoom In Zoom Out Reset image size

The radially averaged function $\bar{C}\left(r\right)$ decays exponentially, as visible by a constant slope in the log-linear plots of figures 4(c), 5(a) and later. We find that there are two regions of exponential decay with different correlation lengths, clearly visible in figure 5(b). We characterize the exponential decay

$$\begin{eqnarray}&&\bar{C}\left(r\right)=A\exp \left(-r/{\ell }\right),\end{eqnarray} \tag{ 7 }$$

by fitting the correlation length ${\ell }$ in equation (7) over a select section of the correlation function. At short distance $r\,\lt \,{r}_{{\rm{cut}}},$ the characteristic length ${{\ell }}_{{\rm{short}}}$ appears to be affected by the choice of ${r}_{{\rm{cut}}}$ within which the non-affine part of the local deformation field is computed, as presented in more detail below. The initial decay crosses over to a second exponential at distances $r\gt {r}_{{\rm{cut}}}$ with a characteristic length ${{\ell }}_{{\rm{long}}}$ that does not depend on the specific choice of ${r}_{{\rm{cut}}}$ and reference frame and is a characteristic of the material under investigation.

Zoom In Zoom Out Reset image size

We first focus on the behavior at short scales. The computation of $\bar{C}\left(r\right)$ involves the cutoff radius r_cut as a parameter. r_cut determines the local atomic neighborhood within which ${D}_{{\rm{\min }}}^{2}$ is calculated. To test whether the length scales ${{\ell }}_{{\rm{short}}}$ and ${{\ell }}_{{\rm{long}}}$ depend on this length, we parametrically vary r_cut between r_cut = 2 r_NN = 5.6 Å and r_cut = 5 r_NN = 14 Å. The resulting correlation functions at 7% and 12% applied strain are shown in figures 5(a) and (b), respectively. The radius r_cut varies by a factor of 2.5 while the individual correlation functions move systematically upwards. As a consequence, the extracted value ${{\ell }}_{{\rm{short}}}$ depends systematically on r_cut. Indeed, we can collapse all ${{\ell }}_{{\rm{short}}}$ values on a single curve when normalizing by ${r}_{{\rm{cut}}},$ ${{\ell }}_{{\rm{short}}}/{r}_{{\rm{cut}}}$ (figure 5(c)). The behavior of ${{\ell }}_{{\rm{long}}}$ is different. Its value (figure 5(d)) is independent of ${r}_{{\rm{cut}}}$ used in the computation of ${D}_{{\rm{\min }}}^{2}$ and the data does not collapse when normalized accordingly. The evolution of ${{\ell }}_{{\rm{short}}}$ and ${{\ell }}_{{\rm{long}}}$ with applied strain $\varepsilon$ clearly show the point where the samples yield (see also figure 3(a)). At around 12.5% strain, ${{\ell }}_{{\rm{long}}}$ increases dramatically. During further deformation it fluctuates around a value consistently higher by a factor of 4 than before yield.

The computation of $\bar{C}\left(r\right)$ furthermore depends on how the reference frame for the computation of ${D}_{{\rm{\min }}}^{2}$ is chosen. Here, we report results obtained for reference frames at constant distance in applied strain, ${\rm{\Delta }}\varepsilon .$ All results reported above were obtained for ${\rm{\Delta }}\varepsilon =1 \% .$ Figure 6 demonstrates how the correlation function and ${\ell }$ vary as a function of this parameter. Before yield, the correlation function does not depend on ${\rm{\Delta }}\varepsilon$ and drops exponentially over two decades as a function of distance. Figure 6(a) shows this behavior for ${\rm{\Delta }}\varepsilon =1 \% ,$ $0.1 \%$ and $0.01 \%$ which is above, at and below the cage-breaking strain ${\varepsilon }_{c}$ (figure 2(b)). The behavior changes at yield (figure 6(b)), where the initial exponential drop starts to depend on ${\rm{\Delta }}\varepsilon .$ Figure 6(c) shows the influence on the extracted value of ${{\ell }}_{{\rm{long}}}.$ ${{\ell }}_{{\rm{long}}}$ decreases with decreasing ${\rm{\Delta }}\varepsilon$ and saturates at ${{\ell }}_{{\rm{long}}}\approx 15\,\AA$ in the flow region for the lowest ${\rm{\Delta }}\varepsilon =0.01 \% .$

Zoom In Zoom Out Reset image size

To clarify the role of ${\rm{\Delta }}\varepsilon$ on the calculation of the correlation length ${\ell },$ we further tested the influence on system size on the correlation functions. Figure 7(a) shows that before yield (applied strain $\varepsilon =7 \%$), $\bar{C}\left(r\right)$ is independent of system size but that a clear size dependence develops when the material flows (figure 7(c), $\varepsilon =20 \%$). Plotting ${{\ell }}_{{\rm{long}}}$ versus applied strain $\varepsilon$ shows that the before yield ${{\ell }}_{{\rm{long}}}$ is independent of size but after yield it depends on system size (figure 7(b)). Normalizing distance r or correlation length ${{\ell }}_{{\rm{long}}}$ by system size collapses all data in the region where the amorphous solid flows (figures 7(c) and (d)).

Zoom In Zoom Out Reset image size

To gain further insights into the system size dependence of ${{\ell }}_{{\rm{long}}},$ we took one of our medium-sized simulations at 20% applied strain and replicated it parallel (A) and perpendicular to the shear band (B, see insets in figure 8), creating simulation boxes with aspect ratios ≈2. We continued straining these 'supercell' systems up to 50% applied strain. System A continued to shear along the replicated shear band. System B had initially two shear bands, but one of the shear bands disappeared during shear in favor of the other, such that the final system had only a single shear band (inset in figure 8). While ${{\ell }}_{{\rm{long}}}$ of system A remains unaffected, ${{\ell }}_{{\rm{long}}}$ of system B rises to twice its initial value, like the cubic systems with twice the edge length.

Zoom In Zoom Out Reset image size

Finally, we test the dependence of the correlation function on temperature and composition. Figure 9(a) shows the temperature dependence of ${{\ell }}_{{\rm{long}}}.$ Data in the temperature range from 50 to 300 K, below the glass transition temperature of ${T}_{g}\approx 600\,{\rm{K}}$ of our metallic glass, is superimposed for small strain. It appears that at large strain the highest temperature leads to a smaller ${{\ell }}_{{\rm{long}}},$ while the values for the other temperatures are approximately the same. Figure 9(b) shows ${{\ell }}_{{\rm{long}}}\left(\varepsilon \right)$ 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.

Zoom In Zoom Out Reset image size

The correlation length ${\ell }$ characterizing the exponential decay of the spatial-autocorrelation functions $\bar{C}\left(r\right)$ of ${D}_{{\rm{\min }}}^{2}$ have in the past been interpreted as giving the size of the STZ [26]. Our results clearly show that the decay of $\bar{C}\left(r\right)$ with distance $r$ is exponential in MD calculations of BMGs, confirming other results obtained for EAM [26], Lennard-Jones [27] and hard-sphere glasses [29, 30]. However, there are two regions of exponential decay with different correlation lengths. At short distance $r\,\lt \,{r}_{{\rm{cut}}},$ the characteristic length ${{\ell }}_{{\rm{short}}}$ is strongly affected by the choice of ${r}_{{\rm{cut}}}$ within which the non-affine part of the local deformation field is computed. Our results indicate ${{\ell }}_{{\rm{short}}}\propto {r}_{{\rm{cut}}}$ such that ${{\ell }}_{{\rm{short}}}$ does not characterize any intrinsic material scale. The initial decay crosses over to a second exponential at distances $r\gt {r}_{{\rm{cut}}}$ with a characteristic length ${{\ell }}_{{\rm{long}}}$ that does not depend on the specific choice of ${r}_{{\rm{cut}}}$ and reference frame and is a characteristic of the material under investigation. For the CuZr glasses investigated here we find ${{\ell }}_{{\rm{long}}}\sim 5\mbox{--}10\,\AA .$ This is on the order of the values reported for FeP in [26] $\left({\ell }=8.5\,\mathring{\rm{A}}\right)$ but smaller than the values for MgAl $\left({\ell }=11.1\,\mathring{\rm{A}}\right)$ and CuZr $\left({\ell }=15.0\,\mathring{\rm{A}}\right)$ reported there at an applied strain of $\varepsilon =4 \%$ for simulations carried out with an earlier version of the EAM potential used here [46]. Additionally, the work described in [26] used the initial configuration at $\varepsilon =0$ as reference (and hence ${\rm{\Delta }}\varepsilon =\varepsilon =4 \%$) and looked at correlations of global non-affine displacements rather than ${D}_{{\rm{\min }}}^{2}.$ Recent work using a Lennard-Jones model for CuZr reports ${\ell }=5\,\AA$ [27]. This appears to indicate that the actual value of the correlation length is highly model-dependent and may also depend on the preparation of the glass. As an extreme example, for a poorly tempered system that does not show shear bands, we would not expect the correlation length to increase suddenly and become system size dependent at the onset of yield. For our systems, we find that the values extracted from our calculations are robust to variations of temperature and stoichiometry.

The situation in the pseudoelastic regime before yield is characterized by individual regions of large ${D}_{{\rm{\min }}}^{2}$ (figure 1(b)) that are typically attributed to individual STZs. Therefore, $\bar{C}\left(r\right)$ measures the autocorrelation of the deformation field of an individual STZ. Since the overall density of STZs is low, the strain offset ${\rm{\Delta }}\varepsilon$ that determines over how many STZs we average does not affect the results. The situation changes dramatically after the sample has yielded ($\varepsilon \gt 10 \%$). STZs are now localized within a shear band and it becomes difficult to identify individual STZs (figures 3(c) and (d)). The onset of shear-banding is then accompanied by a characteristic length ${{\ell }}_{{\rm{long}}}$ proportional to the system size $L$ that depends on strain window ${\rm{\Delta }}\varepsilon .$ For strain windows smaller than the cage-breaking strain, ${\rm{\Delta }}\varepsilon \lt {\varepsilon }_{c},$ we find values for ${\ell }$ comparable to the ones found in the elastic regime (figure 6(c)). This is because even for the flowing glass with a shear band we can identify individual STZs if we look at small enough strain windows. Figures 3(e) and (f) show examples of the distribution of regions of large ${D}_{{\rm{\min }}}^{2}$ for ${\rm{\Delta }}\varepsilon \lt {\varepsilon }_{c}$ that clearly show individual STZs. The correlation length ${{\ell }}_{{\rm{long}}}$ computed for ${\rm{\Delta }}\varepsilon \lt {\varepsilon }_{c}$ therefore, like in the pseudoelastic regime, characterizes the size of an individual STZ rather than the correlation along the shear band. Note that Cubuk et al [32] have chosen to evaluate ${{\ell }}_{{\rm{long}}}$ at a value of ${\rm{\Delta }}\varepsilon$ that corresponds to the minimum in the ${{\ell }}_{{\rm{long}}}({\rm{\Delta }}\varepsilon )\,$behavior obtained during flow of the material. Figure 6(c) indicates, that ${{\ell }}_{{\rm{long}}}\left({\rm{\Delta }}\varepsilon \right)={\rm{const}}$ for ${\rm{\Delta }}\varepsilon \lt {\varepsilon }_{c}$ and a system size dependence for ${\rm{\Delta }}\varepsilon \gt {\varepsilon }_{c}$ in our calculations. Cubuk et al are therefore likely in the region where ${\rm{\Delta }}\varepsilon \lt {\varepsilon }_{c}.$ They find exponential correlation in all cases and their extracted scale is independent of analysis protocol and system size.

The actual value of ${{\ell }}_{{\rm{long}}}$ depends on the size perpendicular to the shear band. Our 'supercell' calculations show that duplicating a shearing simulation (with an existing shear band) parallel to the shear band (system A in figure 8) does not affect ${{\ell }}_{{\rm{long}}},$ while duplicating it perpendicular to the shear band (system B in figure 8) also doubles ${{\ell }}_{{\rm{long}}}.$ Note that in the latter case only one of the two shear bands survived during consecutive strain. The scaling ${{\ell }}_{{\rm{long}}}\propto L$ is therefore simply a consequence of the ratio between the volume occupied by the band (it is of a thickness on the order of the size of the STZ that does not depend on system size) and the cell volume, or in other words, ${{\ell }}_{{\rm{long}}}$ measures the distance between individual shear bands in the superlattice of shear bands created by the periodic boundary conditions. For an experimental system of a size much larger than the STZ size, this length scale would look divergent.

We would like to note that our system sizes, although large, are yet too small to rule out power-law behavior during flow. Indeed, the fact that for ${\rm{\Delta }}\varepsilon \gt {\varepsilon }_{c}$ our length scale ${{\ell }}_{{\rm{long}}}$ depends on system size for is indicative of a diverging length scale. This could be a signature of a cross-over to a power-law as STZ events become correlated within the shear band. This observation is consistent with a recent proposition that yield in amorphous solids can be interpreted as a first-order phase transitions [47, 48], an interpretation that has a rich history for explaining shear-banding instabilities in non-Newtonian fluids [49]. Jaiswal et al [47] identified the transition using an order parameter that measures similarity or 'overlap' of atomic configuration. The atomic configuration loses 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 ${\ell }$ rises after the stress has peaked. As a guide to the eye, the gray vertical band in the stress–strain curve (figure 3(a)) indicates the applied strain where ${{\ell }}_{{\rm{long}}}$ rises rapidly (see also figures 7(b), (d)).

We note that power-law correlations can be found in other measures for the non-affine displacement field rather than ${D}_{{\rm{\min }}}^{2},$ for example by looking at the global non-affine displacement field rather than at the non-affine displacements in augmentation spheres of radius r_cut as employed for the computation of ${D}_{{\rm{\min }}}^{2}.$ In the elastic regime, the displacement-difference correlation function of a disordered body shows power-law behavior that describes the elastic Green's function [50, 51]. In the plastic regime, power-law scaling persists albeit with a different exponent that is compatible with a self-affine geometry for the deformation field [35]. Similar scaling has been observed for correlation of the strain field [34]. However, this type of scaling cannot be detected in the bare correlation function of the non-affine displacement (only in the difference correlation); and it appears to be not manifested in the local deviation from affinity as described by ${D}_{{\rm{\min }}}^{2}.$

We studied the correlation between non-affine displacements, as characterized by the ${D}_{{\rm{\min }}}^{2}$ measure of Falk and Langer [15], using MD calculations. This multipoint correlation function shows exponential behavior in the elastic regime from which we can extract a length scale ${\ell },$ typically attributed to the size of an STZ. We find that this length scale diverges at yield, as manifested by a size-dependent ${\ell }$ in during flow of the material. The divergence of ${\ell }$ 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 [47, 48].

We thank Suzhi Li and Jan Mees for useful discussion and the Deutsche Forschungsgemeinschaft (DFG) for funding this work through grant PA 2023/2. All simulations were conducted with LAMMPS [52] on JURECA at the Jülich Supercomputing Center (grant 'hfr13') and on NEMO at the University of Freiburg (DFG grant INST 39/963-1 FUGG). Post-processing and visualization was carried out with ASE [53], matscipy [54] and Ovito [55].

Appendix. 

Spatial correlations can be computed straightforwardly by implementing equation (2) directly, with the $\delta$-function broadened into discrete bins. This calculation, however, becomes prohibitive for large distances, since the number of pairs of atoms scales with the distance squared. In order to speed-up the calculation (for large distances), we map our field Q onto a regular grid

$$\begin{eqnarray}&&{Q}_{x,y,z}=\displaystyle \sum _{i}{Q}_{x,y,z;i},\end{eqnarray} \tag{ 8 }$$

where x, y and z are integers denoting points on a grid of N_DFT = N_xN_yN_z grid points. We chose Q_x,y,z;i = Q_i for $\vec{R}=x{\vec{a}}_{1}+y{\vec{a}}_{2}+z{\vec{a}}_{3}$ closest to the position of atom i, ${\vec{r}}_{i},$ and zero otherwise. Other choices like a Gaussian interpolation or a distribution over all neighboring grid points are possible, but not used here. They only affect the results at small distances. Note that ${\vec{a}}_{i}$ are the cell vectors that describe the tilted cell.

The discrete Fourier transform, and the inverse transform, of Q_x,y,z are given by

$$\begin{eqnarray}&&{\rm{D}}{\rm{F}}{\rm{T}}\{{Q}_{x,y,z}\}=\displaystyle \sum _{x,y,z}{Q}_{x,y,z}\exp \left(-i{k}_{x}x-i{k}_{y}y-i{k}_{z}z\right),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rm{I}}{\rm{D}}{\rm{F}}{\rm{T}}\{{\tilde{Q}}_{x,y,z}\}=\displaystyle \frac{1}{{N}_{{\rm{D}}{\rm{F}}{\rm{T}}}}\displaystyle \sum _{x,y,z}{\tilde{Q}}_{x,y,z}\exp \left(i{k}_{x}x+i{k}_{y}y+i{k}_{z}z\right),\end{eqnarray} \tag{ 10 }$$

with ${k}_{i}=2\pi /{N}_{i}.$ The autocorrelation map in reciprocal space is obtained from the transform of equation (2)

$$\begin{eqnarray}&&\tilde{{\mathscr{A}}}\left[Q\right]\left(\vec{k}\right)=\displaystyle \frac{V}{{N}_{p}^{2}}{\rm{D}}{\rm{F}}{\rm{T}}\left\{{Q}_{x,y,z}\right\}{\rm{D}}{\rm{F}}{{\rm{T}}}^{* }\left\{{Q}_{x,y,z}\right\},\end{eqnarray} \tag{ 11 }$$

which in real space yields

$$\begin{eqnarray}&&{\mathscr{A}}\left[Q\right](\vec{R})=\displaystyle \frac{{N}_{{\rm{D}}{\rm{F}}{\rm{T}}}}{{N}_{p}^{2}}{\rm{I}}{\rm{D}}{\rm{F}}{\rm{T}}\left\{{\rm{D}}{\rm{F}}{\rm{T}}\left\{{Q}_{x,y,z}\right\}{\rm{D}}{\rm{F}}{{\rm{T}}}^{* }\left\{{Q}_{x,y,z}\right\}\right\}.\end{eqnarray} \tag{ 12 }$$

The tilt of simulation box is accounted for automatically, by using the appropriate cell vectors in the Fourier transformation. Specifically, the phase of equations (9) and (10), ${k}_{x}x+{k}_{y}y+{k}_{z}z=\vec{R}\cdot \left({k\vec{{\rm{g}}}}_{1}+{l\vec{{\rm{g}}}}_{2}+{m\vec{{\rm{g}}}}_{3}\right)$ where ${\vec{g}}_{i}$ are the reciprocal lattice vectors. Our reciprocal lattice vectors are normalized such that ${\vec{a}}_{k}\cdot {\vec{g}}_{l}=2\pi {\delta }_{kl}.$ Care needs to be taken that when computing radial averages of ${\mathscr{A}}\left[Q\right](\vec{{R}}),$ the tilt of the box vectors is appropriately accounted for.

The periodic boundaries of the system, in reciprocal space, also pose no problem as all pairs of grid points at a certain (real space) distance one from another collapse onto the same reciprocal grid point. Note that the fidelity of the autocorrelation at short distances depends on the choice of the size of the grid and the interpolation prescription. We here use ${\rm{\Delta }}x=1\,\AA$ for the radially averaged correlation functions and use ${\rm{\Delta }}x=0.5\,\AA$ for the correlation maps shown in figure 4. We systematically checked implementation of equation (12) against a brute-force estimate of equation (2) to check that as ${\rm{\Delta }}x\to 0,$ we obtain the brute force estimate at short distances. We have implemented this algorithm in matscipy [54] and Ovito [55] and it is available in the public repositories of these codes.