Micro-plasticity in a fragile model binary glass

Atomistic deformation simulations in the nominally elastic regime are performed for a model binary glass with strain rates as low as $10^{4}$/sec (corresponding to 0.01 shear strain per 1$\mu$sec). A robust elasticity is revealed that exhibits only minor elastic softening, despite quite different degrees of structural relaxation occurring over the four orders of magnitude strain rates considered. A closer inspection of the atomic-scale structure indicates the material response is distinctly different for two types of local atomic environments. A system spanning iscosahedrally coordinated substructure responds purely elastically, whereas the remaining substructure also admits microplastic evolution. This leads to a heterogeneous internal stress distribution which, upon unloading, results in negative creep and complete residual-strain recovery. A detailed structural analysis in terms of local stress, atomic displacement, and SU(2) local bonding topology shows such microscopic processes can result in large changes in local stress and are more likely to occur in geometrically frustrated regions characterized by higher free volume and softer elastic stiffness. These insights shed atomistic light onto the structural origins that may govern recent experimental observations of significant structural evolution in response to elastic loading protocols.


I. INTRODUCTION
The amorphous solid is a material strongly out of equilibrium and very much dependent on its thermal and loading history [1,2]. Despite this, structural glasses have a robust and extended elastic regime [3], and a well defined and high yield stress [4,5], giving them realworld industrial and engineering applications. This robustness and reproducibility implies a strongly constrained microstructure and therefore a spatially-correlated structural disorder from which structural length-scales can emerge [6,7]. Such emergent structural heterogeneity, which is difficult to characterize experimentally, strongly influences the macroscopic plasticity of the amorphous solid [8][9][10].
Classical thermally-activated micro-plasticity during the nominally elastic region, offers direct insight into a class of fundamental excitations of the solid that exist also under zeroload conditions [11]. Signatures of such microplastic activity can, for example, be indirectly seen in micro-tensile experiments [12], creep and damping [13,14], isostatic macroscopic loading [15][16][17], fatigue [18], or stress-relaxation experiments [18]. The attribute "indirect" is used here since deformation in the absence of resolvable shear-band formation, is resolved as either non-linear pre-yield behavior, a permanent residual strain upon unloading, or an increase in the stored excess enthalpy. More specific examples are the rejuvenation and a local increase in nano-scale elastic heterogeneity after cyclic loading [18], or the residual creep strain after applying a constant load at low homologous temperature [15,16].
Such studies reveal changes in measurable quantities that imply structural changes which, within the precision of experiments, can hardly be resolved. This inability to track detailed structural changes can be somewhat overcome with novel scattering methods using both x-rays and/or electrons, where for example fluctuation microscopy inside a transition electron microscope [19] or speckle patterns from coherent x-ray scattering [20] can shed light on the atomic scale structural change or dynamics. The present authors pursued such an approach in very recent work by tracking the structural dynamics during nominally elastic loading [21] and during cryogenic protocols [22] that activate structural excitations underlying the Gamma-mode [23]. Both these efforts reveal a rich and non-trival signature of structural dynamics quantified with time-dependent relaxation times that hint towards thermally-activated microplastic processes heterogeneous in both space and time. However, the quantified signals remain an average of a large number of atoms, and a detailed picture of the spatio-temporal characteristics of individual structural excitations is lacking. This is where atomistic simulations have for a long time played a pivotal role by providing model based understanding of the microscopic structure and processes of metallic glasses, with a primary focus on binary and ternary model glass systems. Such efforts strongly support the early ideas of Frank [24], Chaudhari and Turnbull [25] and Nelson [26,27], which asserted a short-range order dominated by 5-fold coordinated bonds characterized by localized structural motives with (or close to) an icosahedral symmetry. Indeed, the fraction of local icosahedral environments is a strong measure of the degree to which a glass is relaxed [28][29][30][31][32][33][34], and the constraints associated with their connectivity is one avenue to understand medium range order and the associated emergence of structural length-scales [34,35]. Simulation has also given fundamental insight into the underlying microscopic processes during stress-driven athermal plasticity, demonstrating that mechanisms which mediate the collective athermal macroscopic response are fundamentally local [36][37][38]. Such localized activity is compatible with early thermal plasticity theories [39,40] and their coarse-grained simulation implementations [41][42][43][44][45][46]. More recently, potential energy landscape exploration methods also revealed a local picture for thermally-activated structural excitations under both zero [47][48][49][50][51] and applied load [52] -a result confirmed by finite temperature molecular dynamics simulations [31,32,53] suggesting the possibility of glass plasticity driven by thermal fluctuations at sufficiently long deformation time scales.
Here we employ finite temperature molecular dynamics at the micro-second timescale, to demonstrate that a well relaxed binary model glass system exhibits a robust elastic regime largely insensitive to any structural relaxation occurring during the elastic loading. We do this for a well known Lennard-Jones system [54], whose resulting glassy structure forms a thermally-driven negative creep that leads to full recovery of the residual strain at a timescale comparable to the initial deformation.

II. SIMULATION METHODOLOGY
A. Choice of Lennard-Jones potential, 50:50 chemical composition and physical units The binary Lennard-Jones potential of Wahnström [54] has been used extensively in atomistic simulations of model binary glass systems [31-33, 49, 50, 52-54, 57, 58] and is widely known as a model fragile glass former. Since it is a pair potential, it is unable to describe the unsaturated nature of the metallic bond, and therefore the quantitative aspects of material specific binary model alloys. For this, embedded atom [59] or secondmoment [60] empirical many-body potentials are needed. Despite this aspect, the potential captures the essential structural physics of bulk binary metallic glasses as for example CuZr, ZrNi, NiNb, HfCu, CaAl etc., and the competition between liquid-like and crystal-like bonds which is believed to underlie such metallic glasses and its medium range order [24,25,28,35].
Fundamentally this is due to the scale of local density variations, which whilst important for characterizing the structural state of the amorphous solid, are too small to distinguish the explicit many-body aspect of the material specific potential. For systems in which there are large changes in coordination, such as that at a surface or an internal cavitation, such Lennard Jones potentials would not be able to capture the bond stiffening expected for a metal. These aspects have been discussed in more detail in Refs. [31][32][33]53].
The present work employs a 50:50 chemical composition of small to large atoms, which is the potential's eutectic composition, therefore producing a homogeneous undercooled liquid and a glass with little chemical segregation.
The Lennard-Jones potential has the following form: where σ 22 = 5/6σ 11 and σ 12 = σ 21 = 11/12σ 11 for the Wahnström parametization. The atoms of type 1 may be considered as the larger atom type. The atomic masses of the two atom types are arbitrarily chosen such that m 1 /m 2 = 2. For a molecular dynamics iteration, a time step of 0.002778σ 11 m 1 /ε is used. The distance unit is taken as σ = σ 11 and the energy unit as ε, with stresses in the units of ε/σ 3 . For this work, the potential is truncated to a distance 2.5σ.
When metallic units (representative of say, CuZr) are taken for ε and σ, an MD timestep is of the order of a femto-second. Thus one billion MD iterations corresponds to approximately one micro-second. Throughout this paper, simulation times will be measured with respect to one billion MD iterations, which in turn is approximated as one micro-second.
Absolute temperatures are expressed as an energy, k B T , using a value of k B = 8.617 × 10 −5 ε where ε is assumed to have units of electron-Volt.

B. Sample preparation and deformation simulations
The sample preparation protocol exploits an initial NVT ensemble to produce a glassy structure at a local minimum in the PEL [31-33, 50, 52]. This sample has undergone 2  [32]. sponse for all strain rates, and at shear strains above 0.02-0.03, an increasing strain-rate effect is seen, which leads to a strain-rate dependent peak-stress response and a transition to plastic flow. This plastic regime will be considered in a forthcoming publication. A closer inspection of the elastic region reveals a weak shear softening (reduction of effective elastic shear modulus) with decreasing strain rate. Fig. 1b plots the internal energy per atom as a function of shear-strain. Unlike Fig. 1a, a strong strain-rate dependence is now seen in the nominally elastic regime, where for the lowest strain-rate, the internal energy significantly decreases. For increasing strain-rates, the internal energy curves tend to approach the expected elastic response, 1/2Gγ 2 where G is the shear modulus derived from the 10 7 strain data (derived from strain values below 0.01).
Atoms are now classified into two groups as a function of the shear strain: those icosahedrally and non-icosahedrally coordinated using Voro++ [61] within the OVITO software package [62]. For our well-relaxed sample, approximately 25% of the smaller atoms are icosahedrally coordinated. A connectivity analysis of the icosahedrally coordinated atoms reveals a dominant system spanning cluster consisting of connected fragments of C15 Laves backbone polyhedra. This is a known structural feature of the model LJ glass [33,57] and also of small clusters in a model CuZr glass [55,56,63]. For each group of atoms, an average potential energy per atom is calculated as a function of shear strain. These potential energies are obtained by performing a congjugate gradient relaxation to a local inherent state con- Also shown in these figures is the expected elastic response, 1/2Gγ 2 . Fig. 1c shows that the average energy of the icosahedral atoms follow more closely the elastic response at all shear-rates with only a weak dependence on strain rate, whereas Fig. 1d confirms that the non-icosahedrally coordinated atoms underlie the strain-rate dependent energy relaxation seen in Fig. 1b. The red data set in panels b, c, and d describe the evolution of a later discussed zero-load configuration (see Sec. III C).
To investigate the degree of micro-plasticity occurring during the deformation, samples are unloaded to zero shear stress at various strains. The unloadings were performed at the upper strain rate of 10 7 /sec, to ensure minimal structural relaxation and therefore maximal probing of the internal stress state at a particular strain. Fig. 2a displays the resulting residual strain for all the strain rates. The figure reveals that micro-plasticity occurs within the nominally elastic regime for all strain rates, and increases (in terms of its residual strain signature) as the strain rate reduces. Within this elastic regime, the residual strain increases with strain and for the lowest strain rate the residual strain begins to increase more rapidly beyond a strain of 0.02 suggesting the onset of the flow regime at these micro-second timescales. The observation of micro-plasticity is compatible with the findings entailed in Fig. 1a, suggesting that the icosahedral regions respond purely elastically, whereas the non-icosahedral regions undergo some form of plasticity that surprisingly results in structural relaxation. To gain further insight into the origin of this micro-plasticity, the spatial loadbearing properties of the glassy structure are now investigated.
To determine the spatial origin of the evolving global shear stress, the virial stress σ µν = 1/(2V ) ij F µ ij r ν ij (where F ij and r ij are the force and distance vectors between atoms i and j, V is the volume of the simulation cell, and the factor of 1/2 is used to correct for double counting) is written as two terms, one involving a summation of i over icosahedrally coordinated atoms and the other over non-icosahedrally coordinated atoms. These two terms are used individually and divided respectively by the fraction of icosahedrally (Ω I ) and nonicosahedrally (Ω NI ) coordinated atoms. This allows determining the average local stress for each of the two atomic environments individually. Doing so, modifies the 1/V pre-factor to give a representative volume of the icosahedral and non-icosahedral regions, and therefore an estimate of the corresponding stresses. An average of the two quantities, weighted according to Ω I and Ω NI , gives the exact global stress of the system. Fig. 2b shows their evolution as a function of strain for the 10 4 /sec strain-rate sample derived from the corresponding cg-quenched inherent-state configurations. The data reveals that different shear stresses are carried by both regions of the material, with the non-icosahedral regions carrying a lower stress. The observed difference between the icosahedral and non-icosahedral stress responses decreases with increasing strain rate (data for other strain rates not shown), as suggested in Figs. 1c and d which show similar elastic responses at the higher strain rates. Thus the difference in shear stress between the two regions increases with decreasing strain rate, and is therefore the origin of the elastic softening seen in Fig. 1a.

B. Negative creep and strain recovery
The picture which therefore emerges is that, upon loading, the micro-plasticity occurring within the non-icosahedrally coordinated regions facilitates a reduction in stress (and energy) within these regions. Thus, upon unloading, a back-stress develops and produces the residual strains seen in with decreasing strain rate is strong evidence for thermally-activated micro-plasticity. If this is indeed the case, the unloaded heterogeneous micro-structure should experience negative creep, with the non-icosahedral regions facilitating micro-plasticity to reduce the negative stress which in turn will reduce the positive stress of the non-icosahedral regions, homogenizing the internal stress to eventually exhaust the negative creep. This has the experimental implication that time-dependent strain recovery may be present after a constant force experiment within the nominally elastic regime, evidence of which can be found in Refs. [13,17]. C. Comparison to zero-load evolution The observations that 1) significant structural relaxation occurs during deformation at low enough strain rates, 2) this occurs in non-icosahedral regions which also undergo a reduction in shear, 3) residual strain exists upon unloading indicating micro-plasticity and 4) negative creep is seen resulting in recovery towards the pre-loading geometry, all point towards a strong thermally-activated micro-plasticity component.
Because of the inherently non-equilibrium nature of the glass structure, structural relaxation under zero-load also occurs via the thermal activation of microscopic scale processes.
These localized structural changes are mediated by localized structural excitations (LSE).
LSEs, can be observed in atomistic simulations at finite temperature within the dynamical heterogeneities of the under-cooled liquid [64][65][66][67][68][69], and also well below the glass transition temperature of the amorphous solid both [31,33,50,53]. They involve a core region in which atoms rearrange, often in a string-like fashion in which one atom replaces the position of a neighbouring atom. This localized activity is accommodated by a far-field Eschelby strain field which has a similar energy-scale as the LSE core [52]. Thermally-activated LSEs should be distinguished from athermally stress-driven shear transformation zones [36][37][38]70] arising from inflections in the PEL giving core structures significantly more localized than that of LSEs.
These considerations motivate the questions also eluded to in the introduction: i) are there differences between the structural relaxation occurring under load and the relaxation processes that would normally occur under zero load? and ii) would such zero-load structural relaxation be related to the thermally-activated micro-plasticity observed in Fig. 1?
To investigate the nature of structural relaxation without load, isothermal simulations are performed for 5 micro-seconds at zero fixed pressure and at 0. This structural relaxation is quantified in Fig. 5, which plots the evolution of the icosahedral content during loading as a function of shear strain for all strain rates. The data has been corrected for a general decrease in icosahedral content as a function of shear strain.
This smooth reduction in icosahedral content, which is reversible upon unloading, also exists in purely athermal deformation and originates from the Voronoi based detection algorithm not being able to identify some strongly strained local icoasahedral environments. For the highest 10 7 /sec strain rate, Fig. 5 shows that the corrected icosahedral content slightly decreases, indicating weak rejuvenation and a likely reason for why the energy per atom increases more rapidly than that expected from elasticity (see Figs. 1b and d). For the two intermediate strain rates, the icosahedral content remains approximately constant. However for the 10 4 /sec strain rate, the icosahedral content increases until a strain of 0.01-0.02, and then plateaus at higher strains. Also shown is the icosahedral evolution for the zero-load sample plotted as a function ofε × t for the 10 4 /sec strain rate, displaying a similar evolution in icosahedral content. Thus the initial rise in icosahedral content is due to the initial transient relaxation regime seen in both zero-load and loading simulations (Fig. 4).
Returning to questions i) and ii), structural evolution in terms of energy relaxation is similar for a zero of finite applied shear stress, with both being mediated by thermally activated LSE activity. For the case of an applied load, the LSE activity clearly results in a net shear stress reduction with the non-icosahedral regions. How this similarity and difference might arise will be discussed in Sec. V.

A. Local stress measure distribution
The results of the loading and unloading simulations reveal a developing internal stress heterogeneity that correlates with the icosahedral and non-icosahedral atomic environments.
How is this stress heterogeneity manifested at the local microscopic level? Insight into this question may be seen by rewriting the expression for the global stress as where N is the number of atoms and V atom is the average volume per atom. σ µν i should be viewed as a measure of atom i's contribution to the global stress, where the average value gives the global stress. It may also bee seen as an estimate of the actual local atomic stress.   Fig. 6a. This is done for the zero-load and unloaded sample by performing the summation in Eqn. 2 but only including terms within the range ±γ, and plotting the summation as a function of γ. Fig. 6b demonstrates that for small values of γ, which includes only the most probable part of the distribution, the global stress remains quite small. Deviations away from this begin in the tails of the The results of this section indicate that the stress heterogeneity associated with the icosahedral and non-icosahedral regions manifests itself in differences in the local atomic stress distributions, where the distribution is consistently narrower for the icosahedral regions.
On the other hand, the average stress difference of these regions originates from large, but rare, local stress contributions at the tails of these distributions. This result suggests that much of the internal stress state of the material remains inactive, with only certain regions undergoing structural change associated with relaxation and/or micro-plasticity.

B. Localized structural excitations
The present data, indicates that a similar thermally-activated relaxation scenario occurs for both zero-and finite-load conditions. Fig. 7a  To quantitatively confirm this similarity in LSE structure, atomic displacements occurring between a time interval, ∆t, are binned to obtain the van Hove self-correlation function, G(r, ∆t) = δ(r−|r i (∆t)−r i (0)|) , a quantity that is experimentally accessible using coherent photon correlation spectroscopy methods [21]. Fig. 7c plots the resulting correlation function obtained from the configurations prior to load and at a shear strain of 0.02 for all strain rates corresponding to a respective ∆t of 2.0 nsec, 20.0 nsec, 200.0 nsec and 2 µsec. All correlation functions demonstrate two structural features: 1) a narrow peak whose maximum is close to zero with a width that spans displacements below the characteristic bond-length (here ∼ 1σ), and 2) a peak at this characteristic bond-length which exhibits an exponential tail at larger displacements. The intensity of this latter peak increases with decreasing strain rate (increasing ∆t). Fig. 7c also includes the van Hove self-correlation function at zero-load for ∆t = 2µsec which due to a similar time interval may be directly compared to the 10 4 strain rate data, and which demonstrates a close overlap with the lowest strain rate. A similar overlap can be seen for the higher strain rates when using the appropriate ∆t (not shown).
In addition, comparison of the second peak associated with bond-length displacements, for the different strain rates show a similar overall structure -the main effect being that the intensity of the second peak increases with decreasing strain rate. More generally, the character of thermally-activated LSE activity is the same, whether under load or not, where the second peak is associated with LSE core atom bond length (or greater) displacements. The observed exponential tail at larger distances is a characteristic of collective behaviour and conjectured to be a universal feature of the glass transition regime [71]. The present work extends this property to the thermally-activated LSE processes below the glass transition regime. The sub-bond-length displacements associated with the first peak correspond to the accommodating LSE Eschelby strain field indicating an intrinsic landscape of small barrier energies which has been viewed as a collective generalization of rattling-in-the-cage excitations [72].

C. Spatial visualization
To visualize the spatial region defined by the icosahedrally coordinated atoms, the initial configuration is shown in Fig. 8a, which displays only the icosahedral atoms. Here, blue/red atoms correspond to the small/large atoms. A surface mesh encloses these atoms. Inspection of the figure demonstrates the dominance of the smaller atom and that the connected cluster of atoms is system spanning. To gain further insight into the structural nature of the non-icosahedral regions, a local SU(2) topological analysis is performed, following a similar approach as that done in Ref. [34].
This procedure, uses a modified radical Voronoi tessellation procedure to determine the bond-order of each nearest neighbour pair, defined as the number of common neighbours of the two atoms. The relevance of this description of glassy structure has its origins in Refs [24][25][26][27], where the liquid like 5-fold bond order represents the minimally frustrated packing of atoms around the bond. As a result, this bond is referred to as a defect free bond. On the other hand, 4-fold and 6-fold bonds are not minimally frustrated atomic packings, and are referred to as defected bond orders. From this perspective, the icosahedral environment is an atom with 12 defect-free bonds and thus represents a minimally frustrated short-range-order structure. The geometrical relationship between such defected bonds, how they may occur and be arranged in terms of bond orientation, has been predicted in the work of Nelson [26,27] which demonstrated via an SU(2) algebraic structure, that only certain combinations of such defect bonds could exist at a particular atomic site. The recent work of Ref. [34] shows that a large proportion of atomic environments of the model Lennard-Jones binary glass presently considered follow the rules dictated by the SU(2) algebraic structure.  (2) topologies. [34]. These later topologies entail disclination nodes in which both green and red disclinations intersect according to the rules of the underlying SU(2) algebra (see Refs. [26,27]). Their presence indicates highly bond-frustrated environments. Indeed the work of Ref. [34] found that minimization of 4-fold disclination content indicates a reduction in the bond frustration of the glassy structure (See Fig. 5 of Ref. [34] and the related discussion) and that such environments are characterized by higher free volume and locally softy (possibly negative) elastic stiffness shear moduli, both of which facilitate LSE activity [34].

V. DISCUSSION
The present atomistic simulation work, which employs samples relaxed and deformed at the micro-second timescale, has revealed the emergence of thermally-activated microplasticity within the elastic regime as the strain rate reduces. This micro-plasticity is heterogeneously distributed in regions characterized by enhanced geometric frustration, whereas regions characterized by minimal geometrical frustration represented by a high local icosahedral content deform almost entirely elastically. This heterogeneous micro-plasticity induces a corresponding heterogeneity in the internal stress, which upon unloading, drives negative creep to homogenize the internal stress and entirely recover the residual strain induced by the micro-plasticity. The thermally-activated local structural excitations (LSEs) underlying these deformation mechanisms also results in the relaxation of the glassy structure. Such relaxation also occurs without load, but without the developing internal stress heterogeneity and may be understood from a classical thermal-activation picture of plasticity, as was applied by Spaepen [39], Argon [40] and Bulatov and Argon [41] which all rely on the stochastic occurrence of thermally-activated localised plastic events, biased by the applied load. These works, in turn, exploit the seminal work of Eshelby [74].
Most generally, a change in energy due to thermal activation of an LSE may be written as a sum of a core contribution, E core , and an Eshelby elastic energy contribution [74]: Here V 0 is the volume of the Eshelby inclusion appropriate for the LSE, σ µν T and ε µν T are the rigid far-field elastic stress and strain due to the Eshelby construction, σ µν C is the correction due to internal relaxation around the Eshelby inclusion, and σ µν L is the characteristic external stress at the LSE. Due to its strongly local character, E core is expected to have a much weaker dependence on σ µν L than the Eshelby contribution. Taking σ µν L as a simple shear in the xy plane of magnitude σ, it's contraction with ε µν T gives ε µν T σ µν L = 2σγ T cos 2θ T and an energy which depends on the LSE's Eshelby shear-strain magnitude γ T and orientation θ T within the xy plane. A similar expression for the corresponding saddle-point energy may also be obtained [41] resulting in a barrier energy whose leading order stress dependence will be ∼ σγ T cos 2θ T , giving the general result that thermally-activated LSEs with a far-field strain signature that reduces the applied load (i.e. cos 2θ T < 0) are more likely to occur, reducing the global elastic energy through the global plastic strain increment ε µν T (here γ T ). This picture of a barrier energy dependent on the applied stress has been confirmed through ARTn simulations which investigated how particular barrier energies varied with load, finding that in terms of a shear load, LSE barrier energies would either increase or decrease depending on their far field shear-stress signature [52].
The above heuristic picture of thermally-activated plasticity gives a simple reason as to why similar degrees of energy relaxation and displacement activity are seen when comparing the loaded sample to the zero-load sample -a certain degree of thermal activation will occur irrespective of the presence of an applied load. Under zero load, over time, this activity will result in structural relaxation with no net deformation of the sample. However with an applied load, thermal activity will be biased towards those LSEs which reduce the applied elastic load. The observation that this bias does not affect the degree of thermallyactivated structural relaxation (over a given time interval), infers that the glass structure is sufficiently complex to support a very large number of LSEs with a broad range of far-field strain signatures.
The observation of a developing heterogeneity and how this can affect plasticity has been investigated in past simulation work [75]. Indeed, in the work of Baumer and Demcowitz [75], which considered a model CuNb binary glass at a concentration away from its eutectic, and therefore with strong chemical segregation, the authors found a system-spanning icosahedral network at the interface between the regions rich in Cu and rich in Nb. This occurred at the the glass transition temperature regime, giving a strong analogy to the general phenomenon of gelation. Detailed analysis of their shear loading simulations performed at the glass transition temperature and at a strain rate of 10 9 /sec revealed the non-icosahdral regions (referred to as liquid-like) plasticised. This resulted in a characteristic cohesive energy that did not depend strongly on the applied strain, whereas the icosahedral network responded approximately elastically. Moreover, upon unloading, no residual strain was observed indicating the plastic processes within the liquid-like non-icosahedral regions were entirely reversible. The present work establishes a related phenomenon but at lower temperatures and slower rates.
In particular, the significantly enhanced time scales enables the observation of recovery via classical thermally-activated micro-plasticity and creep. This is nicely demonstrated in Fig. 3, where time-dependent strain recovery (negative creep) proceeds after unloading. Due to the slow loading rate and the ability of the nonicosahedral network to admit thermally-activated structural micro-plasticity, an incompatibility stress (and therefore strain) between both structural components emerges upon unloading. The magnitude of the residual strain must, if true thermally-activated micro-plasticity occurs during loading, depend on temperature, strain rate, and the initial structural state of the glass. These observations provide an atomistic understanding of various experimental observations of inelasticity and strain recovery after both static and cyclic loading conditions at stresses well below yielding [12][13][14][15]17], which typically are described with viscoelastic Kelvin-Voigt or Maxwell models. The applicability of such phenomenological laws to the here observed atomic-scale strain recovery at the micro-second time scale will be reported elsewhere.
After the elastic energy component has been subtracted from the cohesive energy, the degree of energy relaxation of the loading and zero-load curves is quite similar for the lower strain rates (Fig. 4b). From this, and the observation of strain recovery, one might conclude that little relative rejuventation occurs due to the micro-plasticity. Experimental changes in enthalpy due to rejuventation occurring within the nominally elastic regime are typically a few 100 J/mol (see Refs. [15,16,76]) and can be as high as ∼ 1 kJ/mol [23]. 100J/mol is approximately 0.001 eV per atom and assuming the Lennard Jones parameter ε is at the scale of one electron volt, detecting such relative energy changes in the simulation data becomes difficult because of the fluctuations associated with the finite size of the system. It is however encouraging that the present simulations could entail relative rejuvenation energies not much larger than this typical energy scale of 100 J/mol. In terms of absolute change (not relative to the zero-load sample) significant relaxation has of course occurred over the timescale of the deformation, which is certainly not the case during a typical experimental deformation. Fig.4b reveals a structural relaxation involving a reduction in energy of about 0.02ε per atom. This corresponds to ∼ 2 kJ/mol when ε is taken as an electron volt, which is comparable to what is seen in the fast quenching of a ribbon [77]. Together, these aspects motivate future deformation simulations on more relaxed configurations. For example, the configuration at 5 µsec in Fig. 4a for which all µsec transients have disappeared appears as an ideal candidate for such an endeavor.

Inspection of
Beyond a shear strain of 0.02, Fig. 1a indicates a greater strain rate dependence with an increased shear softening with decreasing strain rate. In addition Fig. 2a shows that with respect to shear strain, this regime entails a more rapid increase of the residual strain upon unloading. These observations, at the micro-second timescale, could be a very early signature of the onset of plastic flow at the shear strain levels where bulk metallic glasses characteristically deform in experiments.

VI. CONCLUDING SUMMARY
The nominal elastic regime has been studied for a well known fragile model binary glass using strain rates spanning four orders of magnitude. The work finds: 1. a robust elastic response with minimal softening for all strain rates considered, despite large differences in the degree of structural relaxation occurring during the loading.
2. the underlying reason for this is a system spanning spatial region which responds elastically, with the remaining parts of the material experiencing structural relaxation and micro-plasticity. These latter processes are responsible for the observed minimal elastic softening.
3. upon unloading, a residual strain remains that fully recovers through negative creep. 4. this strain recovery is driven by a back-stress within the elastically responding spatial regions, eventually leading to the homogenization of the internal stress.
5. zero-load simulations reveal similar structural relaxation, indicating very similar thermally activated processes underlie relaxation and micro-plasticity 6. the microscopic structural features distinguishing the elastic region are a high density of icosahedral coordinated environments and a negligible 4-fold disclination content, both of which indicate a minimally bond frustrated system. 7. thermally activated localized structural excitations (LSEs), mediate both the relaxation and micro-plasticity, and occur mainly in the non-icosaheral regions as evidenced by enhanced atomic displacement during loading and increased disclination content (strong bond frustration).
8. these LSEs result in large but very local changes in internal stress, that are in fact statistically quite rare, reflecting an underlying spatio-temporal heterogeneity that emerges at the micro-second time scale.
The above should motivate new experiments on bulk metallic glasses that specifically probe the temperature, time and strain rate dependence of micro-plasticity and its relation to rejuvenation, as well as further atomistic simulations done at micro-second (and indeed the sub-millisecond) timescales. As pointed out in a recent overview article by the present authors [11], micro-plasticity probes the fundamental thermally activated unit-scale plastic processes that are characteristic of the zero-load structural state. A proper understanding of these local processes is essential when entering the high-stress regime of collective phenomena that ultimately lead to rejuvenation, yield and plastic flow.