Mapping out the glassy landscape of a mesoscopic elastoplastic model

Understanding the response of a disordered solid to an externally imposed forcing, such as stress or strain, is important in order to characterize the transitions between rigid and flowing states in a wide variety of soft matter systems. Examples for such behavior include the jamming transition in granular materials,¹ the yielding transition in amorphous solids,^2,3 and the depinning transition of a pinned elastic interface, such as flux-lines in type II superconductors.⁴ 

The interplay between the deformation energy cost and gain, as the disordered solid adapts to the imposed forcing by deforming, gives rise to rich dynamics on a complex energy landscape. For small loading, the response of the solid is largely elastic, characterized by few plastic deformation events. However, as the loading is increased, plastic deformations start to proliferate and eventually lead to yielding and flow. The manner in which the transition to yielding occurs has been found to strongly depend on the degree of initial annealing, e.g., aging, of the sample.^5–10 

Over the last several years, a large body of experimental^11–15 and numerical^{6,9,10,16–30} work has been carried out to understand the nature of the yielding transition in amorphous solids. These results reveal an intriguingly complex and dynamical spectrum of response that, besides its dependence on the degree of annealing and the amount of loading, also shows dependence on history, system size, and dimensionality.

Of special recent interest has been the response of amorphous solids to cyclic shear, in particular, under athermal quasistatic (AQS) conditions.^{9,10,18,19,23,26,27,31} Experiments and simulations show^{18,19,31–36} that for small oscillatory strain amplitudes, the solid settles into a cyclic response after just a few driving cycles. As the strain amplitude is increased, the transients to cyclic response become increasingly longer, and multi-periodic response, i.e., cycles that repeat every T > 1 driving periods, starts to emerge. This behavior continues until a critical strain amplitude is reached, beyond which cyclic response is no longer attainable and particles start to diffuse across the sample. The transition from cyclic to diffusive behavior has been found to be rather sharp and is called the irreversibility transition.

At the same time, a cyclic response to periodic loading can also be regarded as a form of memory that encodes information about the forcing that produced the response.³⁷ Such memory effects have been observed experimentally as well as numerically in periodically sheared amorphous solids, colloidal suspensions, as well as other soft-condensed matter systems.^{15,32–34,38–41}

Along with atomistic (AS) models of amorphous solids, spatially coarse-grained mesoscopic elastoplastic models³ have been introduced. Due to their conceptual simplicity, mesoscopic models are appealing both from a numerical as well as a theoretical perspective. Initially, the main goal of these models was to capture the response under monotonous loading by shear strain. However, more recently, mesoscopic models have been constructed that study the response under oscillatory shear.^42,43 In order to be able to realistically capture cyclic response, particularly key features of the irreversibility transition, a prescription for replacing mesoscopic elements once they have yielded has to be provided. We will refer to the available choices generated by a given replacement prescription as the “landscape” of the mesoscopic model.

Of special interest are the two recently introduced elastoplastic models that study the response under cyclic shear: The model by Liu et al.⁴² assumes that mesoscopic elements that yielded are replaced at random and hence irreversible, while the model of Khirallah et al.⁴³ is complementary in that it is fully deterministic—elements that yielded are replaced by ones with identical, i.e., non-random, yield stresses. The only source of randomness is the initial internal stress configuration. Thus, in terms of the landscape terminology, the model of Liu et al. has a totally random disorder-landscape, while the one of Khirallah et al. is totally ordered. Despite these differences, both models nevertheless recover key features of the response of amorphous solids to cyclic shear, such as the irreversibility transition and divergence of lengths of transients as the transition is approached. Let us finally note that these types of mesoscopic models have been used as a starting point for developing even further coarse-grained models, such as the recently introduced stochastic mesostate models.^44–46 These models also capture key features of the irreversibility transition of amorphous solids under oscillatory shear.

Here, we present a depinning-like mesoscopic elastoplastic model with a quenched disorder landscape. Our model, therefore, interpolates between the two types of landscapes considered before. Specifically, the model we consider has two features: (i) a local yielding protocol that allows us to mimic thermalization and aging, and thereby to tune the history of our samples; (ii) the quenched disordered landscape that allows us to capture, in rather great detail, the transients and the evolution to cyclic response in terms of the localized plastic events.

As in previous works,^42,43 we first focus on the stress response under monotonous loading by an externally applied shear strain. Our model recovers the brittle-to-ductile transition: as our initial glass is increasingly aged better, the stress response exhibits a stress peak that gets more pronounced with the duration of aging.

We next focus on the irreversibility transition under oscillatory shear and its dependence on both the degree of annealing and system size. We find that, for poorly and moderately aged (PA and MA) samples, the transient times to cyclic response diverge as the irreversibility transition is approached. In the case of poorly aged samples, this divergence follows a power-law with an exponent that is comparable with estimates obtained in recent studies.^19,22,23,43

We finally turn to a more detailed comparison between the disorder landscape of mesoscopic and atomistic models. To this end, we make use of the fact that the AQS dynamics of driven disordered systems has a natural representation in terms of a transition graph, the t-graph.^39,41,47,48 The AQS dynamics is thereby encoded into the topology of the t-graph and provides a unified setting within which we can compare in great detail the properties of the disorder landscapes underlying our mesoscopic and atomistic models.

We perform such comparisons by focusing on a particular topological feature of the t-graph, its strongly connected components (SCCs). An SCC is a collection of mechanically stable configurations, actually elastic branches, which are connected in a bi-directional manner by plastic deformation pathways: a pair of configurations belongs to the same SCC if there is a deformation pathway that leads from one to the other and back. Hence, the plastic events triggered by transitions between states belonging to the same SCC are mechanically reversible, while transitions connecting different SCCs are irreversible.⁴¹ Any periodic response must necessarily be confined to a single SCC, and therefore the size of the SCCs and their dynamic accessibility is a limiting factor for the length of transients and the cyclic strain amplitudes at which cyclic behavior can be attained.⁴¹ 

We have organized the manuscript in two main sections and sub-divided these into subsections. In Sec. II, we present the mesoscopic elastoplastic model, its behavior under monotonous and cyclic loading, and the dependence upon glass preparation. In Sec. III, we present a characterization of the underlying disorder landscape of the model based on the use of transition graphs. Throughout the two parts, when applicable, we also compare our results qualitatively with atomistic simulations we have carried out. We conclude with a discussion of our results. We conclude this introduction by providing the plan of the paper along with a summary of our main results.

We present in Sec. II A the mesoscopic model developed for the present study. In Sec. II B, we detail the protocols of preparation that allow us to mimic annealing at high temperature and aging at vanishing temperature, respectively. In Sec. II C, we show that varying the level of aging allows us to recover upon monotonous loading either a ductile response or a brittle one, where a stress peak is followed by a softening branch. In Sec. II D, we focus on the analysis of the irreversibility transition upon cyclic driving. In particular, we discuss the dependence of the transition on sample size and preparation by aging. The poorly aged glasses show a power-law divergence of the transient time to cyclic response as the strain amplitude approaches the irreversibility transition. We find that this behavior changes qualitatively as the samples are aged better and their sizes are sufficiently large. In this case, the irreversibility transition appears to be discontinuous. Cyclic response is attained rather quickly for amplitudes below a critical strain, or not at all.

In Sec. III, a characterization of disorder landscape is given with an emphasis on the prevalence of limit cycles. Our analysis is based on the transition graphs, as recently proposed in the context of atomistic simulations.³⁹ The representation of Athermal Quasi-Static (AQS) dynamics via t-graphs extracted from simulations is introduced in Sec. III A. Their extraction from simulations of sheared amorphous solids is presented in Sec. III B. The graph topology and the crucial role of the strongly connected components (SCCs) in the context of cyclic loading are discussed in Sec. III C. In Sec. III D, we discuss the effect of glass preparation on the topological properties of the t-graphs. In particular, we show the scale free character of the size distribution of SCCs. In Sec. III E, we show that in order to understand better the disorder landscape underlying the differently aged glasses, one has to combine the topological properties of the t-graph with structural information, such as the stability range of the SCCs. This combined information, as well as the study of the evolution of the plastic strain, is discussed in Sec. III F, and we report on an interesting age-induced symmetry breaking transition associated with a phase-separation-like process in the disorder landscape.

Note that throughout the manuscript, we present qualitative comparisons between results obtained from our mesoscopic and atomistic simulations. Let us stress that our main aim in showing such comparisons is not to reach a quantitative agreement. Rather, our intention has been to use the molecular dynamic results as a guide and a qualitative reference against which the different results obtained in the present work are compared.

The paper concludes with a discussion followed by a series of appendices. Methodological details of atomistic and mesoscopic simulations are summarized in  Appendixes A and  B, respectively. An estimate of the effective sizes of the atomistic and mesoscopic simulations, in order to allow us to compare these, is given in  Appendix C. Finite-size effects on the stress response upon monotonous loading are discussed in  Appendix D. Details on the properties of the catalogs used for assembling the transition graphs, which were extracted from simulations of the atomistic model and the mesoscopic models, are presented in  Appendix E. A discussion of the strip-like arrangement of SCCs on the plane of exit strains is given in  Appendix F, while in  Appendix G, we provide a transition graph perspective for the dependence of the irreversibility transition on finite-size and aging of the glasses.

We consider a scalar 2D lattice-based mesoscale elastoplastic model. The physics of this class of models relies on the coupling between a threshold dynamics and an elastic interaction induced by the incremental local plastic slip that arises as a result of a mechanical instability.³ We use here a variant of the model introduced in Refs. 49–51. A detailed presentation is given in  Appendix B. Here, we give a brief introduction to the model and emphasize its novel features, focusing on the properties of the stress landscape.

We consider a square grid of N × N cells of size a × a. The model is scalar so that we account for one and only one shear direction, along which we can shear the system forward (FW) and backward (BW). We assume a uniform shear modulus μ. Each individual cell (i, j) is characterized by a stack of local elastic branches indexed by a variable ℓ, each of which relates the local stress σ_ij to the local strain ɛ_ij, as shown in Fig. 1. The stability of each such local elastic branch ℓ is limited by two bounds: a maximum stress threshold $σij,ℓ+$ and a minimum stress threshold $−σij,ℓ−$⁠. Note that in order to ease notation, whenever no explicit reference to a particular branch number ℓ is made, we will omit it in the following. The two thresholds $σij+$ and $σij−$ are drawn from a random distribution with support in R⁺ so as to ensure $−σij−<σij+$⁠, i.e., the existence of a stability domain for the cell (i, j). Whenever the local stress σ_ij overcomes one of the two bounds, the cell experiences a plastic event and its stability domain is shifted to a neighboring elastic branch. Since the cell is surrounded by other cells and can be seen as an Eshelby inclusion within an elastic matrix,⁵² this plastic event induces a stress redistribution in the system so that other cells can get destabilized.

Details about the implementation and the driving of the model are given in  Appendix B, but we summarize below the main novel features of the present model with respect to the previous variants presented in Refs. 49 and 51:

Bidirectionality—Since cyclic loading (in addition to a simple monotonous loading) is considered here, two local thresholds are defined instead of only one for each cell: one threshold $σij+$ in the forward direction and the other one $σij−$ in the backward direction.

Annealed vs quenched disorder—In the case of a monotonous loading, every time a cell experiences a plastic deformation, its threshold is renewed. Independently of the particular method chosen to draw random thresholds, the quenched character of the disorder is automatically obtained since a unique value of the threshold $σij,ℓ±$ is assigned to each triplet $(i,j,εijpl)$⁠. The possibility of back and forth motions requires more care in the definition of the threshold disorder. Here, when performing cyclic loading tests, we use a quenched disorder. In practice, we resort to a counter-based random number generator (see Ref. 54 for a pedagogical introduction) to assign efficiently a unique pair of random thresholds to each triplet without having to store the full sequence in memory.

Preparation of the system before cyclic loading—The structure and the mechanical behavior of glasses do depend on their thermo-mechanical history. In order to account for this preparation dependence in mesoscopic elastoplastic models, one usually specifies a particular distribution of local thresholds and/or internal stress in the initial configuration.^5,53 In contrast to atomistic simulations, these distributions do not derive from a well-defined quench protocol but have to be introduced by hand. Here, we propose two simplistic protocols of preparation allowing us to mimic (i) an instant quench from a high temperature liquid and (ii) aging at vanishing temperature. Although they are un-realistic caricatures of actual glass preparations, the combination of these two protocols allow us to tune continuously the state of the system from a very disordered fresh soft glass to a very aged hard glass in a (statistically) reproducible way.

As explained earlier, the present model is stress based and relies on threshold dynamics: plasticity sets in at cell (i, j), if and only if the local stress overcomes one of the two thresholds in the positive or negative shear directions: $σij>σij+$ or $σij<−σij−$⁠. Despite the absence of an explicit energy landscape, which would allow us to equilibrate the system at finite temperature and to subsequently perform a quench to zero temperature,⁵⁴ it is possible to implement two limit-cases of glass preparation: instant quench from a high temperature liquid and aging at vanishing temperature, respectively.

At a high temperature, the local energy barriers associated with the stress thresholds are very low with respect to the available thermal energy so that, in the T → ∞ limit, all plastic rearrangements are equally probable. We then define a thermal step by selecting a site uniformly in space at random and choosing one of the two directions with probability 1/2. The chosen site, thus, experiences a plastic slip and jumps onto a new elastic branch, which is characterized by two new plastic thresholds. Next, the stress field is updated to account for the stress redistribution. The stress redistribution can make some other sites mechanically unstable and thereby induce an avalanche. Updates are then performed until the avalanche stops and the system is stable again. The system is subjected in this manner to a sequence of thermal steps until it reaches a stationary state. In Fig. 2(a), we show, for different system sizes N, how the mean stress-threshold $σ+̄$ of our samples evolves with the number of thermal steps. We see that when plotted against the average number of thermal steps per site, the curves for the different sizes collapse and $σ+̄$ reaches a stationary value rather quickly, after about 4–5 thermal events per site. The inset of the figure shows the corresponding evolution of the standard deviations δσ⁺, δσ, and δx⁺ of the stress-threshold, the internal stress, and the plastic strength, respectively. When plotted against the average number of thermal steps per site, we find again little size dependence. In Fig. 2(b), we show the stationary distributions of the stress-thresholds, internal stress, and local plastic strengths for our N = 64 sample.

We now turn to the other limit, namely aging at very low temperature, T → 0. In the framework of activated behavior, the activity at low temperature is restricted to overcoming the lowest barriers. Moreover, in the limit of vanishing temperature, the lowest barrier becomes dominant. We define an extremal aging step as follows: recall that for each site (i, j) its plastic strength in the positive and negative directions are given as $xij+=σij+−σij$ and $xij−=σij+σij−$⁠, respectively. We identify the site and direction with lowest plastic strength and let it experience a local slip so that stresses are redistributed, and new stress thresholds are assigned to the yielded site. As in the case of the “thermal” procedure with randomly selected sites, a stability check is performed after each slip. If one or more sites get unstable, they are updated in turn and with the most unstable sites updated first, as explained before. The procedure is iterated until the avalanche triggered by the initial extremal step terminates. Then, the next site and direction of lowest plastic strength is identified and allowed to slip.

The present “aging” procedure is thus similar to the “thermal” procedure, differing only in the choice of the initial site to be slipped: in the case of “aging” an extremal site is selected for slip, i.e., the cell and direction with least plastic strength, while in the thermal case the selection of site and direction is random. This difference drastically alters the dynamics, since it induces a systematic statistical bias. When a site yields, it acquires a new pair of thresholds. The latter are drawn from a prescribed distribution. However, in the framework of the aging procedure, this takes place at an extremal site, which is characterized by a very low plastic strength (either in the positive or in the negative direction). We thus get a typical exhaustion phenomenon: low thresholds get replaced by “normal” ones. This systematic bias induces a drift in the threshold distributions and thus a systematic plastic hardening.^49,55

Starting from an initial state corresponding to the inherent state obtained from a “high temperature liquid,” as described in Sec. II B 1, we thus “age” the system by slipping a number of least stable sites. As shown in Fig. 3(a), we observe a logarithmic growth of the mean thresholds $σ+̄$ with the number of aging steps. Again, the dependence of this evolution on system size becomes negligible when we consider the average number of aging steps per site, instead of the total number of steps. We find that after about 10³ aging steps per site, the mean threshold doubles in value.

The inset of the figure shows the evolution of the standard deviation of the stress-threshold, internal stress, and plastic strength. The standard deviation of thresholds shows a slow decrease (about 20% over 10³ aging steps per site). Together with the doubling of the mean thresholds over the same range of 10³ steps, this corresponds to a significant narrowing of the threshold distributions upon aging.

Interestingly, after a fast decrease in the early stage of the aging protocol (less than one aging step per site), the standard deviation of internal stress remains almost constant upon aging. In recent studies on the dependence of plastic behavior of amorphous solids on glass preparation,^6,53 the width of the stress fluctuation distribution has been used as a proxy for the level of stability of the amorphous solids while keeping the value of the plastic threshold constant (actually uniform). We get here a different situation: an increase in the mean threshold and stability of the stress fluctuations upon aging. A way to reconcile these contrasting observations is to consider the fluctuations of the local plastic strength x^± = σ^± ∓ σ and to note that in the case of uniform thresholds the standard deviation of plastic strength equals that of internal stress δx^± = δσ. Upon aging, we indeed observe a continuous decrease of δx^± that gets halved after about 10³ aging steps per site.

In Fig. 3(b), we display distributions of the stress thresholds σ⁺ for our N = 64 samples, which were either not aged at all (thermal), or aged at 0.8, 15, and 150 aging steps per site, for N = 64. These aging levels have been indicated by the appropriately colored circles on the graph showing the evolution of mean stress-thresholds with aging in panel (a). Henceforth, we will refer to these levels of aging as poorly aged (PA), moderately aged (MA), and well-aged (WA).

The effect of our aging procedure is dramatic: it opens a growing gap in the distribution of stress-thresholds σ⁺. In spirit, we recover here a phenomenology that is close to that of ultrastable glasses obtained via swap Monte Carlo methods.⁵⁶ The opening of a gap will induce a perfect elastic behavior over a finite range of strains that contrasts with the quasi-elastic behavior (short elastic branches punctuated by plastic events) typically observed in less equilibrated glasses.

Depending on glass preparation, the stress–strain curves show either a monotonous behavior up to a plateau or exhibit a stress peak followed by a softening branch that slowly approaches the stress plateau at a steady-state stress Σ_ss. The existence of a stress peak is usually associated with shear-banding behavior.

In Sec. II B, we proposed a glass preparation protocol for our mesoscopic model that mimics aging at vanishing temperature. While tuning an aging duration is very different from tuning a quench rate from the liquid state, both methods allow us to transit continuously from a soft/poorly equilibrated glass to a hard/well equilibrated glass. Our protocol actually allows us to obtain in this way very different glassy states. In Fig. 4, we show stress–strain curves corresponding to a poorly aged, a medium-aged, and well-aged glass, aged at an average of 0.8, 15, 150 number of steps per site. The system size is N = 32, and the curves were obtained by averaging over 500 realizations. While the poorly aged glass does not exhibit a stress peak, such a peak emerges and becomes more pronounced as the samples are aged more. Thus, by tuning the duration of aging, we are able to transit from a poorly aged to a well-aged glass. We checked that all curves do converge to the same stress plateau for large enough shear strains. More details on the size dependence of these stress–strain curves are shown in Fig. 14 of  Appendix D.

Comparison with atomistic simulations—In the inset of Fig. 4, we show for reference two stress–strain curves obtained by atomistic simulations under athermal quasi-static shear for a slow and fast quench, respectively. The details of the simulations are provided in  Appendix A. The slow quench curve shows a distinct stress peak, while apart from fluctuations, the fast quench curve is almost monotonous. Due to computational time limitations, it is difficult to obtain strongly contrasting quenches and consequently stress–strain curves when using molecular dynamics for the glass preparation. The recently developed swap Monte Carlo methods give access to a wider range of glass preparation although they are more restrictive with respect to the nature of the model glasses.⁵⁶ 

Let us emphasize that it has not been attempted here to adjust the parameters of the elastoplastic model to quantitatively reproduce the stress–strain curve obtained by atomistic simulations. Rather, our goal is to compare generic features, such as the brittle to ductile transition under monotonous loading, and how the behavior upon cyclic loading depends on the soft/hard nature of a glass. Recent analyses of coarse-graining atomistic simulation to be used to feed mesoscopic elastoplastic models with realistic parameters can be found in Ref. 57.

In this section, we consider the irreversibility transition and, in particular, the response to cyclic shear of our poorly aged (PA) and moderately aged (MA) mesoscopic glasses whose preparation was described in Sec. II B. The well-aged (WA) mesoscopic glasses yield a response to cyclic shear that is qualitatively similar to that of the (MA) glasses and will therefore not be considered in this section.

When subjected to cyclic shear loading, amorphous solids tend to either evolve into periodic response or reach a diffusive regime depending on the value of the amplitude ɛ_amp of the loading cycles. This transition presents typical features of a critical transition. In particular, power-law divergence of the number τ of loading cycles to reach the periodic response below the transition, as well as the power law dependence of the diffusivity above the transition have been observed both for atomistic and mesoscopic models.^{9,10,18,19,23,26,27,31} The features of the irreversibility transition depend on glass preparation.^9,10,29,30 Figure 5 shows the response of an N = 32 sample of a poorly aged mesoscopic glass to cyclic shear at strain amplitude below, panel (a), and above, panel (b), of the irreversibility transition. In the former case, a cyclic response was obtained after τ = 158 cycles, while for the latter, cyclic response was still absent after 10⁴ driving cycles. The main plots show the evolution of stress and strain over the last 30 cycles, each of which have been color-coded in increasing shades of red, as indicated by the legend in (a). The lower insets show a detail from the main plot, while the upper insets show the evolution of the difference of stresses at the beginning of two consecutive driving cycles. Below the irreversibility transition, panel (a), this stress difference eventually vanishes (after τ = 158) while above the transition in panel (b) it keeps showing finite fluctuations at least until τ = 10⁴. In the following, we show results for the size dependence of the irreversibility transition in our PA and MA mesoscopic glasses. Specifically, we consider systems of size N = 16(7500), 32(3000), 64(400), where the numbers in parenthesis indicate the number of realizations used to obtain our results.

We first focus on the poorly aged (PA) systems. Figure 6 shows the mean success-rate p_succ, i.e., the fraction of PA systems (circles) within our ensemble of realizations that reach a limit cycle when subject to a given number τ_max of symmetric loading cycles at amplitude ɛ_amp: 0 → ɛ_amp → 0 → −ɛ_amp → 0. The different colors correspond to the system sizes, as indicated in the legend of the figure. For system sizes N = 16, 32, and 64, we used a cut-off of τ_max = 10⁴ driving cycles so that if cyclic response had not been established at that point we considered the run to be unsuccessful. A clear transition can be observed between a low amplitude regime with convergence to a limit cycle and a high amplitude regime with no limit cycle. The transition between these two regimes gets increasingly sharper with system size. A clear size dependence is also observed in the location of the transition that tends to occur at lower strain amplitudes for larger systems. The size effect exhibited by our poorly aged glasses is all the more striking as it turns out to be completely absent in the response to monotonous loading, and only weakly present in the case of our moderately and well-aged glasses (Fig. 14 in  Appendix D).

For each size N, we estimate the strain ɛ_irr(N) at which the irreversibility transition occurs, as the loading amplitude such that 50% of the realizations reach a limit cycle, i.e., p_succ = 1/2, as indicated by the pink horizontal line in Fig. 6. The inset of Fig. 6 shows the size and aging dependence of ɛ_irr(N) for N = 16, 32, and 64 for the PA, MA, and WA glasses. We see that for a given degree of aging, ɛ_irr(N) decreases with increasing system size. Moreover, a dependence of ɛ_irr on aging at fixed system size is clearly visible, in particular, for the larger sizes N = 32 and 64. At these sizes, the MA glasses have slightly larger ɛ_irr than the PA ones, while the WA glasses have overall larger values of ɛ_irr for all system sizes considered. The behavior of ɛ_irr with aging is consistent with atomistic simulations of cyclically sheared amorphous solids that show that the strain marking the onset of the irreversibility transition is largely independent of aging for sufficiently poorly aged samples, but that it starts to increase as the samples are better aged.^9,10

We turn next to the response of our moderately aged (MA) glasses to cyclic shear. The diamond symbols in Fig. 6 show the fraction p_succ of MA glasses in our ensembles of realizations that reach a limit cycle when subject to cyclic loading of amplitude ɛ_amp. Similarly to the poorly aged samples, as the system size is increased, the irreversibility transition exhibits an increasingly sharper decline of the success-rate from one to zero. However, for a given system size, the rapid fall-off of the success rate in the MA glasses occurs at consistently larger strain values than for the PA glasses, which is in agreement with the behavior of ɛ_irr discussed earlier.

Another feature of the irreversibility transition is the divergence of the duration of the transient regime: atomistic simulations show that the number of loading cycles needed to reach the limit cycles diverges as a power-law according to τ(ɛ_amp) ∝ |ɛ_irr − ɛ_amp|^−α, as shown in Refs. 19, 22, 23, and 43.

In Fig. 7(a), we plot τ(ɛ_amp) against |ɛ_irr(N) − ɛ_amp| for our poorly aged glasses and different system sizes N. Here, ɛ_irr(N) is the loading amplitude at which half of the realizations reach limit cycle, as defined previously. Once again, a significant size effect is observed: for a given |ɛ_irr(N) − ɛ_amp|, the larger the system size, the shorter the transient regime. An indicative power-law behavior of exponent α = 2.8 is plotted as a dashed line. We see that the results obtained for N = 16, 32, and 64 are reasonably consistent with this trend over roughly one decade for the larger samples. Note that the value α = 2.8 is close to the estimate of α ≈ 2.7, recently reported in Ref. 43, as well as α ≈ 2.6, which was obtained using atomistic simulations by Regev et al.^19,23 It is also close to the value α ≈ 2.66 obtained by Corté et al.³³ for a simplified model of interacting particles under flow. The saturation observed for large values of τ(ɛ_amp) naturally stems from the hard limit associated with the finite number of loading cycles τ_max = 10⁴ that we used in our numerical simulations for N = 16, 32, and N = 64.

In Fig. 7(b), we also plot the period of the limit cycle T(ɛ_amp) against |ɛ_irr(N) − ɛ_amp| for our poorly aged glasses and different system sizes N. As already observed in Ref. 43, we see that the limit cycles get more and more complex, with an increasing period when the amplitude ɛ_amp of the cyclic loading approaches the irreversibility transition ɛ_irr. For illustrative purposes, we show that the fast increase of the period is consistent with a power law behavior T(ɛ_amp) ∝ |ɛ_irr(N) − ɛ_amp|^−β with β = 1.5 plotted as a dashed line in Fig. 7(b).

In Fig. 8, we show the same observables τ(ɛ_amp) and T(ɛ_amp) close to the irreversibility transition, now for the moderately aged glasses. For small systems sizes (N = 16, 32), we again observe a diverging trend in the transient duration and the limit cycles period. It appears actually that the larger the system size, the narrower the range of amplitudes over which this diverging behavior holds. Another behavior gradually becomes dominant: for large system sizes, a limit cycle is reached after just a few loading cycles, and the response is mainly elastic. Moreover, as can be seen for the N = 64 glass in Figs. 6 and 8, the transition to irreversibility is rather abrupt and discontinuous. The system either reaches a T = 1 cyclic response rather quickly or not. These findings are consistent with results reported in the literature, e.g., the work by Bhaumik et al.¹⁰ where the authors consider simulations of a 3D amorphous solid subject to cyclic shear. Figure S5(a) of in their paper’s supplement shows the evolution of the transient for a well-aged sample by monitoring the average potential energy per particle. Depending on the shear amplitude, the transient is either very short or a cyclic response is not attained at all.  Appendix G contains additional details on the dependence of the irreversibility transition on system size and aging of our mesoscopic glasses, which is complemented by the results of the transition graph description of the disorder landscape and which we turn next to.

In order to characterize better the disorder landscape underlying our mesoscopic model, we turn next to the transition graph (t-graph) representation of the dynamics under AQS shear.⁵⁸ As was shown recently,^39,41 such t-graphs can be extracted from atomistic simulations of sheared amorphous solids. Features of the AQS dynamics, such as yielding and return point memory, are thereby encoded in the topology of the t-graph.^39,41,47 Thus, t-graphs provide useful information about the underlying disorder landscape. At the same time, the representation of AQS dynamics via t-graphs extracted from simulations provides a unified framework within which we can compare the dynamics of atomistic as well as mesoscopic models in a rather direct and comprehensive manner.

To fix ideas, we consider first the sheared amorphous solid in an atomistic setting. Under AQS conditions, a given mechanically stable particle configuration can be sheared in the positive and negative direction until a mechanical instability occurs. Denoting by ɛ^± the critical values of the external shear strain at which the instability sets in, for shear strains between ɛ⁻ and ɛ⁺, the configuration of particles deforms smoothly and reversibly in response to the applied shear strain. These sets of mechanically stable particle configurations constitute an elastic branch of the system, which we simply refer to as a mesostate.³⁹ We will use capital letters to label mesostates and denote the critical strain values of a mesostate A by ɛ^±[A]. When ɛ = ɛ⁺[A] (or ɛ = ɛ⁻[A]), a fast relaxation to a new mechanically stable particle configuration occurs. This particle configuration must necessarily be part of another mesostate, i.e., belong to a different elastic branch, say B. Thus, the instability at ɛ = ɛ⁺[A] triggers a transition from mesostate A to B. A similar transition occurs when ɛ = ɛ⁻[A]. The transition between mesostates can therefore be represented in terms of a directed graph, the AQS transition graph or simply t-graph. The vertices of the t-graph are the mesostates, while from each mesostate, we have two outgoing transitions that constitute the directed edges of the graph. We shall denote the transitions when ɛ = ɛ⁺[A] or ɛ = ɛ⁻[A] as the U-, respectively, D-transition out of A. We will refer to the states that these transitions lead to as UA and DA.

The t-graph, along with the critical strains ɛ^±[A] associated with each mesostate, forms a complete representation of the AQS dynamics under arbitrary shearing protocols.⁴⁷ Given an initial mesostate A and a shear protocol, the sequence of mesostates visited can be read off by following the corresponding U- and D-edges, while checking each time whether the critical strains needed to trigger the transition have been exceeded or not.

Note that since UA and DA are mesostates reached from A, their stability ranges must contain the strains ɛ^±[A] at which these transitions were triggered, i.e., we have the AQS conditions⁴⁷ 

Consequently, ɛ⁺[A] < ɛ⁺[UA] < ɛ⁺[U²A] < ⋯, and thus the upper critical strains are monotonously increasing with repeated U-transitions. An analogous result holds for the lower strain threshold under D-transitions. An immediate consequence of this observation is that the sub t-graphs, which are obtained by considering only transitions under U (or D), are necessarily acyclic, i.e., they cannot contain any cycles. Thus, any cyclic behavior must arise from an interplay of the U- and D-transitions.

The numerical algorithm of extracting t-graphs from simulations of sheared amorphous solids has been described in detail in the supplementary material of Ref. 39. Here, we will sketch out the main idea. We start with an initial configuration that is part of a mesostate O, which we call the reference state. We assign to O the generation number g = 0. Next, we execute the U- and D-mesostate transitions out of O, leading to the mesostates UO and DO, and we assign these to generation g = 1. Every time we reach a new mesostate, we compare it to the catalog of mesostates we have obtained so far to see whether it has been encountered before. If not, we add it to our catalog. By proceeding generation by generation, we acquire in this manner a catalog of mesostates: each mesostate A is assigned an ID, its critical strains ɛ^±[A] and the IDs of the mesostates it transits into under a U- or D-transitions are determined. The t-graph is then assembled from such catalogs. In our mesoscopic models, each mesostate corresponds to a configuration of the local elastic branches associated with each of the cells. The event based nature of their simulations facilitates the identification of mesostates and their transitions.

We obtain catalogs from ten realizations each of the N = 32 poorly, moderately, and well-aged glasses, as described in Sec. III A. Tables I, II, and V (in  Appendix E) detail various properties, such as the number of generations and mesostates contained in them. The ranges of strains that these catalogs sample are indicated in Fig. 4, showing how far these catalogs reach out in strain relative to the yield strain under monotonous loading. In addition, we produced ten catalogs from samples of an ultra-stable glass aged by an average of 4 × 10³ steps per-site. For comparison purposes, we also extracted catalogs from our atomistic simulations, using a set of 8 soft and 30 moderately hard reference configurations, that were obtained via fast and slow quenches from a high-temperature liquid. The description of these atomistic catalogs is given in  Appendix E.

Figure 9 shows sample t-graphs from each of the four sets of samples: fast quenched atomistic glass (AS Fast No. 4), slow quenched atomistic glass (AS Slow No. 2), poorly aged mesoscopic glass (Meso PA No. 1), and the moderately aged mesoscopic glass (Meso MA No. 4). The numbers specify the particular realization of the glass, as listed in Tables I, II, VII, and VIII. The placement of the vertices of the graph is arbitrary. The mesostate corresponding to the initially prepared glass, i.e., the reference state, is indicated by the large vertex in red and labeled O. Note the general tree-like structures in all four t-graphs that appear to be qualitatively similar, despite the different underlying model (atomistic vs mesoscopic) and also the different degrees of glass preparation. The color of each vertex indicates the SCC that it belongs to, as we discuss in Sec. III C.

A note of caution when comparing simulations of atomistic and mesoscopic models is in order. As we argue in  Appendix C, our atomistic simulations correspond to an elastoplastic model with size approximately between N = 5 and 10. Thus, our atomistic simulations involve systems of smaller size and possibly suffer more from finite-size effects. Moreover, the way the atomistic systems have been aged is different from the aging protocol used for our mesoscopic systems. All of these features make detailed comparisons difficult and we want to stress again that our main aim in presenting our mesoscopic model is not to quantitatively reproduce the results of atomistic simulations. We will return to this point when comparing the SCC size distributions in Sec. III D.

We will probe the topology of the t-graphs more deeply by focusing on their SCCs to which any cyclic response must be confined,⁴¹ as we explain now. Two mesostates A and B are connected, if on the t-graph there is a directed path of U- and D-transitions that leads from A to B. Physically, this implies that there is some shearing protocol that, when applied to A, gives rise to a deformation pathway terminating in B. We say that two mesostates A and B are mutually reachable, if there is a deformation pathway from A to B as well as one from B to A. Mutual reachability is an equivalence relation (in particular, if the pairs A, B and B, C are mutually reachable, so must be the pair A, C). Therefore, the vertices of the t-graph can be partitioned into equivalence classes under mutual reachability and these classes form its SCCs.⁵⁹ Numerical details on how to extract SCCs from t-graphs have been provided in Ref. 41.

By construction, transitions between any two mesostates belonging to different SCCs are irreversible: there may be a deformation pathway from one to the other, but not back, since otherwise the pair of states would have been mutually reachable. Thus, mutual reachability also partitions the set of transitions between mesostates into reversible ones, i.e., those connecting a pair of mesostates within the same SCC, and irreversible ones, where the two mesostates must belong to different SCCs. Any periodic and hence reversible response to some shear protocol must therefore be confined to a single SCC. The SCCs are thus the “containers” of reversible behavior.⁴¹ 

Tables I and II show the properties of the ten catalogs with N = 32, which were obtained by taking the poorly and moderately aged mesoscopic glasses as reference states. The second column lists the number of generations g_comp up to which all outgoing mesostate transitions were identified. Thus, g_comp = 39 means that we have identified every mesostate that can be reached from the reference configuration by a sequence of 39 U- and D-transitions. Next, N₀ and N_SCC list the number of mesostates and SCCs contained in the catalog. The last row of each table provides the cumulative totals. We will discuss the results shown in the last four columns later in this section.

SCC size distributions—In Fig. 10, we compare the size distribution of the SCCs found in these catalogs. The blue boxes and black circles show the size distribution of SCCs extracted from all ten catalogs of the N = 32 mesoscopic glasses. All curves are normalized but have been vertically offset for clarity. Observe that the size distributions are broad and that the moderately aged catalogs contain larger SCCs. Nevertheless, power-law fits using the method of Clauset et al.⁶⁰ yield a comparable power-law exponent of about 2.3 ± 0.3 for both distributions.⁶¹ For comparison purposes, we also show the SCC size distributions obtained from our atomistic simulations under slow and fast quench, labeled as AS slow (triangles) and AS fast (diamonds), corresponding to moderately and poorly aged glasses. These catalogs reveal similarly broad distributions, with the moderately aged catalogs containing again larger SCCs, while the fitted power-law exponents 2.7 ± 0.3 are comparable.

Note the presence of a finite-size cut-off around SCC sizes of about 30 and 100 for the mesoscopic PA and MA catalogs, respectively. The SCC size distributions obtained from the atomistic simulations do not feature such a cut-off. In  Appendix C, we argue that the mesoscopic equivalent size N corresponding to our atomistic simulations is approximately between N = 4 and N = 10. Thus, the atomistic samples are in effect smaller. Figure 13 of  Appendix C shows the corresponding SCC size distributions when we compare the atomistic simulations with the size distributions obtained from N = 8 mesoscopic catalogs. While still there, the finite-size cut-off appears to be less prominent in the distributions of the mesoscopic models, particularly for the MA samples. We think that the suppression of the cut-off is a finite-size effect. In fact, among the eight AS fast catalogs there are considerable fluctuations in the size s_max of the largest SCC found in each of them. Ordered from smallest to largest, we find s_max = 106, 243, 244, 259, 379, 413, 458, and 929. Among these, the smallest value s_max = 106 is realized in catalog No. 4, whose transition graph is shown in Fig. 9, while the largest value s_max = 929 is observed in catalog No. 2, whose transition graph is given in Fig. 1(b) of Ref. 41. In fact, we checked for AS fast that the data points for the largest SCC sizes in Fig. 10 are singletons corresponding to the largest SCCs found in the catalogs.

Note that while the t-graphs and SCC size distributions obtained from our mesoscopic and atomistic model are qualitatively similar, the dependence of these on the degree of aging is rather weak. In other words, the topology of the t-graphs alone does not appear to contain features that are directly linked to the different amount of aging these samples have been subjected to. As we will show next, the effect of aging on the samples reveals itself when we combine the topological features of the t-graphs with additional physical properties, such as the prevalence of cycles, the strain stability ranges and plastic strains associated with the mesostates and their SCCs.

Prevalence of cycles—We next turn to the population of cycles in our catalogs. We are again interested in cycles that can be traversed under a symmetric cyclic shear protocol: 0 → ɛ_amp → −ɛ_amp → 0 with some shear amplitude ɛ_amp. We consider every mesostate in our catalog that is stable at zero strain and apply this cyclic shear protocol, checking whether a cyclic response has set in or not. The column labeled n_cycles of Tables I and II lists the total number of distinct cycles found in our catalogs obtained from our moderately and poorly aged mesoscopic glasses. We find that the poorly aged catalogs contain a significantly larger number of cycles although the total number of mesostates in these catalogs is comparable (836 082 and 676 760 mesostates, respectively).

As we have noted before, the mesostates forming a cyclic response must all be confined to a single SCC, i.e., a cycle cannot traverse multiple SCCs. We, therefore, ask next how the cycles found in the catalogs are distributed across the available SCCs. In particular, we ask for the number of SCCs that support at least one symmetric cycle, which we define as $NSCCsupp$ and list in Tables I and II. For ease of comparison, we have put together in Table III the cumulative totals listed in the last lines of these tables along with the corresponding data from our atomistic simulations.

Starting with the mesoscopic glasses, there is again a stark contrast between catalogs obtained from poorly aged (PA) and moderately aged (MA) samples (first two rows of Table III). In the MA glasses, the symmetric cycles are contained in a relatively small fraction of SCCs (10 120 out of a total of 218 942 available ones), while for the poorly aged catalogs, a large fraction of SCCs supports at least one such cycle (134 563 SCCs that support symmetric cycles out of a total of 163 587). From Tables I and II, we see that this is true also for the individual catalogs. It is thus apparent that in the moderately aged catalogs a relatively small fraction of SCCs support most of the cycles found, while in the poorly aged catalogs, the opposite is the case and almost every SCC supports at least one cycle. A similar, albeit less pronounced behavior is seen also in our atomistic simulations, cf. the last two rows of Table III. Note that the cumulative data for poorly aged (moderately aged) initial states have been sampled from 8 (30) catalogs (Tables VII and VIII in  Appendix E) so that it is hard to compare the overall number of cycles. Nevertheless, we observe also in our atomistic simulations that the number of cycle supporting SCCs in the poorly annealed catalogs appears to be disproportionally larger.

We finally consider the largest SCCs that support symmetric cycles, comparing their sizes $ssuppSCCmax$ and the number of cycles they contain $ncyclesmaxSCC$⁠. These numbers are shown in the last column of Tables I and II. Again, we find contrasting behavior. The largest cycle supporting SCCs found in the moderately aged catalogs are generally larger than those in the poorly aged ones, but despite this, they contain fewer cycles.

Our results for the prevalence of symmetric cycles can be summarized as follows: while the poorly aged catalogs contain a large number of such cycles that are distributed across a large number of SCCs of various sizes, we find that the opposite is true for the catalogs obtained from the moderately aged samples. For the latter, the number of symmetric cycles contained is far less and these cycles are confined to a small subset of available SCCs.

In order to understand better the difference of the disorder landscape arising from moderately aged and poorly aged samples, we coarse-grain the t-graph to the level of SCCs, since—as we have shown—any cyclic response must be confined to a single SCC. Every SCC has at least one outgoing U- and one outgoing D-transition. Let us denote the states from which these outgoing transitions originate as the U- and D-exits of the SCC. Suppose now that the SCC has only one U- and one D-exit and denote the threshold strains triggering these exiting transitions as $ESCC±$⁠. Consequently, given any mesostate A belonging to that SCC and applying strains confined to the interval $ESCC−<ε<ESCC+$⁠, the resulting sequence of mesostates must remain confined to the SCC. This follows from the observation made before, namely that for any mesostate A, ɛ⁺[A] < ɛ⁺[UA] and ɛ⁻[DA] < ɛ⁻[A].

In the case of multiple U- or D-exits from an SCC, we define $ESCC+$ and $ESCC−$ as the largest, respectively, lowest, strain triggering the outgoing transitions. It actually turns out that for the SCCs considered in our catalogs only a very small fraction of SCCs have multiple U- or D-exits.⁶² Therefore, assuming that each SCC has exactly one outgoing U- and D-transition, consequently in order for the SCC to support cyclic response under the strain protocol 0 → ɛ_amp → −ɛ_amp → 0⋯, we must require that $ESCC+>εamp$ and $ESCC−<−εamp$⁠. In particular, this implies that

Distinguishing the SCCs by (i) their size and (ii) whether they support a symmetric cycle or not, we now ask how these SCCs are scattered in the plane spanned by $ESCC−$ and $ESCC+$⁠. Panels (a) and (b) of Fig. 11 show the SCC scatter plots obtained from single catalogs of our atomistic poorly aged and moderately aged samples, while panels (c) and (d) show the same for catalogs obtained from our mesoscopic poorly aged and moderately aged N = 32 samples. Panels (e) and (f) show SCC scatter plots obtained from even further aged mesoscopic samples, with an average of 150 and 4000 aging steps per site, respectively (details of these catalogs are provided in  Appendix E). In each panel of the figure, the number indicates the particular sample from which the data shown came from. The size of the symbols represent the size of the SCCs, as indicated in the legend, while the boxed symbol shape indicates that the SCC actually supports a limit-cycle, as determined by inspecting our catalogs. The highlighted upper left quadrant of each plot corresponds to the region where the inequality (2) holds. Since this is the region where any SCC that supports cyclic response under symmetric oscillatory shear must be located, we will refer to it as the cycle-quadrant.

We start with a comparison of the poorly aged (PA) and moderately aged (MA) SCC scatter plots obtained from our atomistic and mesoscopic glasses [panels (a)–(d)]. Comparing the catalogs obtained from the PA samples [panels (a) and (c)] with those of the MA samples [panel (b) and (d)], we see that in all cases the cycle supporting SCCs (boxes) are indeed confined to the cycle-quadrant, i.e., the highlighted region in the top left part of the figure, as they should. Moreover, note the relative sparsity of cycle-supporting SCCs in the atomistic (b) and mesoscopic (d) MA samples, when compared with their poorly aged counterparts, panels (a) and (c). This is consistent with our earlier observation, namely that relative to the poorly aged samples, in the MA catalogs only a small fraction of SCCs actually support symmetric cycles.

Plotting the SCCs against their exit strains $ESCC−,ESCC+$ also visualizes possible correlations in the locations of cycle supporting SCCs. For the poorly aged samples, panels (a) and (c), these SCCs fill out the cycle-quadrant rather uniformly and the extent to which this region is filled seems to be limited mainly by the size of the catalog we have sampled, i.e., the number of generations we tracked.⁶³ This is in contrast to the case of the moderately aged samples, panels (b) and (d): not only are there fewer SCCs in the cycle-quadrant, but these SCCs tend to cluster around its boundaries, $ESCC+=0$ and $ESCC−=0$⁠, thereby implying that these SCCs can only support cycles with low amplitudes of a symmetrical shear protocol. In fact, for the mesoscopic samples, we find that the scarcity of SCCs within the cycle-quadrant and their clustering near its boundary becomes even more pronounced when the samples are aged more, as shown in the SCC scatter plots of panels (e) and (f) that were generated from samples that underwent 150 and 4000 aging steps per site, respectively.

Thus, panels (a)–(d) reveal that the SCC scatter plots obtained from our mesoscopic model are qualitatively very similar to their atomistic counterparts: our mesoscopic model captures rather well the difference of the samples due to their aging as well the spatial distribution of the SCCs in the plane of exit strains $ESCC−,ESCC+$⁠.

Before proceeding, we should note that there are sample-to-sample fluctuations in the scatter plots obtained from the individual glasses. This is also apparent in the variation of catalog properties listed in Tables I and II, as well as in the tables for the other catalogs given in  Appendix E. In particular, the spatial population of SCCs in the cycle-quadrant varies from sample to sample. Moreover, within a given sample, the populations of SCCs in the $ESCC−,ESCC+$-plane does not perfectly display the statistical $ESCC±→−ESCC∓$ symmetry that arises under interchange of the forward and reverse shearing directions, even though the number of SCCs shown in these plots are rather large. Nevertheless, the features we have been discussing so far and in the following are typical and appear to be robust from sample to sample.  Appendix F contains a discussion of the strip-like arrangement of the SCCs in the plain of exit strains, which is highly pronounced in the case of the atomistic samples, as well as the better-aged mesoscopic ones.

Having established that the SCC scatter plots are a good proxy to probe topological features of the disorder landscape, we next look more closely at the effect of aging on our mesoscopic glasses. Panels (c)–(f) of Fig. 11 show SCC scatter plots obtained from increasingly better aged samples of our mesoscopic glass which—apart from the PA, MA, and WA samples we considered so far—now includes also an ultra-aged (UA) glass, obtained from a treatment with 4000 aging steps per site. The properties of the 10 catalogs extracted from these glasses are listed in Table VI of  Appendix E.

Note that the moderately aged (MA), well-aged (WA), and ultra-aged (UA) samples each display distinct outlier SCCs in the cyclic quadrant. For the MA sample, these SCCs are located around $(ESCC−,ESCC+)=(−0.05,0.05)$⁠, while for the WA samples, these are found at larger strains. These SCCs turn out to be formed by mesostates that can be reached from the initially prepared glass by strain deformation protocols that do not go beyond the stress-peak and hence do not suffer the subsequent large stress-drop.

To understand why with increased aging the cyclic quadrant becomes less densely populated by SCCs and why these tend to cluster near its boundaries, we consider next the plastic strains. Recall that with each mesostate A we associate an elastic branch in the stress–strain plane. In the case of our mesoscopic model, this branch is by construction linear and the plastic strain ɛ^pl[A] associated with the branch is the (extrapolated) value of the strain where the stress vanishes. By averaging over the plastic strains of the mesostates that belong to an SCC, we obtain a coarse-grained plastic strain for each SCC. The colors of the plot symbols shown in Fig. 11 represent the plastic strains of the SCCs, as indicated by the color table legends. Note that for the mesoscopic samples, panels (c)–(f), we have color-coded the same range of plastic strains. Thus, the shift of colors toward blue and red as the samples get better aged indicates that the magnitudes of typical plastic strains increase with aging.

Moreover, we see that the distribution of plastic strains across SCCs is strikingly different for the differently aged samples. The well- and ultra-aged samples reveal a clear bi-modal distribution of plastic strains, characterized by very few SCCs that have vanishing plastic strains.⁶⁴ For the poorly aged sample [panel (c)], the distribution of plastic strains appears to be unimodal, with a large number of SCCs, particularly those in the cycle quadrant, having plastic strains of very small magnitude.

To understand better the segregation of SCCs by plastic strain, we turn to the mesostates and their deformation histories. Given a mesostate A, we consider the deformation path that leads to it from the initially prepared glass state O. In particular, we are interested in mesostates whose deformation path experiences the stress-peaks and subsequent stress-drops that are encountered during monotonous shearing in the forward or backward direction. We distinguish such states by whether the stress-drop in the forward or backward direction was experienced first, and call these FW, respectively, BW mesostates. In Fig. 12(a), we plot for each level of aging (MA, WA, or UA) the plastic strain distribution extracted from the elastic branches associated with the FW mesostates found in all ten catalogs. We do the same for BW mesostates. For each level of aging, we have thus two distributions of plastic strains: one associated with FW and the other associated with BW mesostates. We see that these distributions are peaked and the location of the peaks move away from each other as the samples are better aged. This can also be seen in the inset of Fig. 12(a) where we plot the averages of these distributions against the aging level. The bars accompanying each symbol indicate the standard deviation of these distributions.

Panels (b) and (c) of Fig. 12 show the deformation path of one particular FW mesostate A from the WA catalog. In both panels, the horizontal axis labels the mesostate transitions starting from the initial glass state O and leading to A. Panel (b) shows the evolution of strain, while (c) depicts the evolution of stress. Black and red circles mark the values at which each transition occurs, respectively, indicating whether the transition happened as a result of a strain increase (black) or decrease (red). The blue vertical lines indicated the extent in strain (b) and stress (c) for each of the elastic branches associated with the mesostates encountered along the deformation path.

The protocol of applied strain that leads from O to A has an initial segment where the strain is monotonously increased to about ɛ = 0.08 (black dashed horizontal line), giving rise to 15 U-transitions. Subsequently, the strain is monotonously decreased to about zero (red dashed horizontal line) over 25 D-transitions. The large FW stress drop of ΔΣ ≈ 0.4 is seen to occur at step g = 7, while the strain is still increasing. The green curve superimposed in panel (b) shows the plastic strain associated with each of the mesostate elastic branches along the deformation path. We see that the stress-drop at step g = 7 is accompanied by a large increase in the plastic strain. As we keep on increasing the driving strain, the plastic strain continues to increase with g, though much more slowly. Remarkably, once we start decreasing the strain again, the plastic strain does not change appreciably. This shows that the changes in plastic strain accrued as a result of experiencing the stress drop are subsequently very difficult to undo, since even a monotonous and prolonged decrease of strain does not seem to change the plastic strain value very much.

We verified that for the better aged glasses, MA, WA, and UA, that all mesostates with an appreciable plastic strain have a deformation history that experiences the stress drop. We thus are able to link the bi-modal nature of the SCC plastic strain distribution to the passage through the corresponding stress peak that is then accompanied by a stress-drop, as demonstrated in Fig. 12(b) (see also the t-graph excerpts shown in Fig. 15 of  Appendix E, where transitions accompanied by larger stress drops have been marked). The better the aging, the larger the stress drops, and hence the larger the jumps in plastic strain, and the more separated are the peaks of the FW and BW distributions.

More importantly, as the evolution of plastic strains in Fig. 12(b) clearly shows, these gains in plastic strain due to the experienced stress drop are apparently very hard to undo by subsequently shearing in the reverse direction. We find that under shearing in the reverse direction the sample has now been significantly softened, i.e., it has become more plastic (the last five transitions from g = 35 to 40 happen at nearly constant strain), indicating a rejuvenation of the sample.²⁸ Thus for the well-aged and ultra-well-aged samples the diagonal $ESCC−+ESCC+=0$ divides the plane of exit strains into an upper and lower half. SCCs located in the upper (lower) half of the plot are SCCs whose mesostates were reached by passing through the forward (reverse) stress peak. This is also consistent with the excerpts from the corresponding transition graphs shown in Fig. 15 of  Appendix E.

We have introduced a depinning-like mesoscopic model of amorphous plasticity characterized by a tunable aging and embedded in a quenched disorder landscape. When driven by an externally applied shear, the model recovers many phenomena exhibited by sheared amorphous solids: a brittle-to-ductile transition under monotonous strain loading, as well as an irreversibility transition under symmetric oscillatory shear, i.e., of the form 0 → ɛ_amp → 0 → −ɛ_amp → 0, and its dependence on the extent the sample has been aged.

We find that the irreversibility transition exhibits a strong dependence on system size as well as on the extent of prior aging of our mesoscopic glasses. Close to the transition, poorly aged systems show a power-law behavior for both the duration of the transient (number of loading cycles needed to reach the limit cycle) and the mean period of the cyclic response (measured in units of the number of driving cycles). Moreover, with increasing system size, the strain value at which the irreversibility transition occurs seems to converge to a well-defined value in the infinite system size limit.

In the case of the better aged samples, we find that cyclic response under oscillatory shear emerges after only a few loading cycles. The dependence on system size is more pronounced in this case. Samples of small size exhibit a cyclic response containing many plastic events and this response continues up to strain amplitudes at which the system would have yielded under monotonous shear loading. However, as the sample size increases, the cyclic response becomes more elastic and the range of strain amplitudes at which it is exhibited shrinks. Changing the system size in our moderately and well-aged samples allows us to gradually transition from a cyclic response whose phenomenology is characteristic of poorly aged glasses to one where this cyclic response is dominantly elastic at larger sizes.

In order to better understand the dependence of the dynamics of our mesoscopic model on the prior aging, we turned next to the study of the transition graphs (t-graphs) that capture the transitions between accessible elastic branches via plastic events. The topology of the t-graphs encodes the dynamics under arbitrary shear loading protocols and thus provides a complementary tool to characterize the disorder landscapes underlying our differently aged systems. We considered a particular topological quantity characterizing the t-graphs, its strongly connected components (SCCs), since any cyclic response has to be confined to a single SCCs. The size distribution of SCCs sampled from both atomistic and mesoscopic simulations of differently aged samples all follow a power-law with an exponent that varies little with the extent of aging but is slightly smaller for the mesoscopic systems than for the atomistic ones.

A closer inspection that also takes into account physical properties associated with the SCCs, in particular their range of stability and typical plastic strains, turns out to be extremely informative. We find that the sample age induces a gradual phase separation between domains of stability centered either on the initial state or at a finite positive or negative plastic strain. The complex age-dependence of the interplay between the amplitude of the center of the domain and the width of the stability ranges has important consequences on the accessibility of limit cycles depending on the particular parameters of the cycling protocols.

Our findings have also implications for memory formation in amorphous solids. Cyclic response under oscillatory shear can encode information and thus form a “memory” about the forcing that caused the response.³⁷ Viewed within the framework of the t-graphs any periodic response must be confined to one of its SCCs. Thus, the evolution under oscillatory shear is primarily a search for a confining SCC. In fact, such SCCs not only contain the cycle forming the cyclic response, but a hierarchy of nested cycles, one of which forms the cyclic response. A hierarchical organization of cycles is typically associated with return point memory.^39,41,47 In particular, the size of an SCC, i.e., the number of configurations they contain, can be regarded as a proxy for memory capacity.⁴¹ 

Since we find that the distribution of SCC sizes is broad, irrespective of the thermal histories of the glasses from which these distributions were sampled, this suggest a high memory capacity even for well-aged glasses. However, a closer look at the stability ranges of the SCCs found in these glasses reveals that only the poorly aged samples have a large abundance of SCCs that can support symmetric cyclic shearing protocols. Contrastingly, in the case of the well-aged glasses, very few SCCs support cyclic response to such oscillatory shear protocols. We find that those that do are characterized by loading/driving histories that did not experience the stress-peak and subsequent stress drop. Consequently, their cyclic response is largely elastic and confined to few and relatively small SCCs.

On the other hand, loading histories in which a stress peak and subsequent stress drop are encountered invariably give rise to rejuvenation of the sample, which is also accompanied by a jump of the plastic strain to non-zero values. As a result, a large number of SCCs become dynamically accessible. However, due to the jumps in plastic strain, these SCCs will only support cyclic response to oscillatory shear if the shear strain is centered on the value of their plastic strain.

Having demonstrated that our mesoscopic model reproduces key features of amorphous solids under variable athermal quasistatic loading, we conclude with a discussion of possible directions for future research. Compared to atomistic models, the computational cost of simulation of mesoscale models is rather low, allowing us to perform extensive numerical computations as well as probing system sizes not accessible to atomistic simulations. In this context, it would be nice to understand better the complex interplay between finite-size effects and the degree of aging that we have observed under oscillatory shear. In the same vein, a detailed statistical analysis of the spatial structure and correlations of sites that undergo plastic activity will be of interest both near the yielding transition and also in the evolution of the transients toward cyclic response under oscillatory shear. In this context, it would be relevant to understand how the spatial structure of sites of plastic activity associated with transitions within an SCC correlates with the size of the SCC and its stability range. In fact, one can regard the set of such active sites as a fingerprint of its SCC and ask how this set changes under transitions to neighboring SCCs, thereby defining an overlap function between SCCs. Since SCCs are containers of periodic response, the strength of such overlaps will have implications for memory formation. Strong overlaps would imply that similar cyclic responses can be realized in neighboring SCCs. At the same time, such overlaps can also be used to characterize in greater detail the topology of the disorder landscape and its possible hierarchical organization.

The authors would like to thank two anonymous referees for their insightful comments and suggestions. M.M. would like to thank Monoj Adhikari and Srikanth Sastry for many useful discussions and for sharing their insights on cyclic shearing and the evolution of plastic strain. I.R. and M.M. acknowledge the kind hospitality of ESPCI, which was made possible in part by Chaire Joliot awards. The authors would like to thank Lila Sarfati and Gaël Tejedor for a careful reading of the manuscript and useful comments. C.E.M., M.M., and D.V. acknowledge the KITP Program Memory Formation in Matter organized in winter 2018 supported, in part, by Grant No. NSFPHY-1748958. This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Grant Agreement No. 754387. M.M. was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Project No. 398962893, the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)— Project No. 211504053—SFB 1060, and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—GZ 2047/1, Project No. 390685813. M.M. also acknowledges partial support by the TRA Modeling (University of Bonn) as part of the Excellence Strategy of the federal and state governments. I.R. was funded by the Israel Science Foundation Grant No. (ISF) through Grant No. 1301/17.

The authors have no conflicts to disclose.

D. Kumar: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). S. Patinet: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). C. E. Maloney: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). I. Regev: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). D. Vandembroucq: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). M. Mungan: Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Atomistic simulations were performed on a two-dimensional binary system with N = 1024 particles of two sizes, where half the particles are 1.4 times larger than the other half. We used a two-body radially symmetric interaction introduced in Ref. 65 and used in Ref. 19, employing the same units of temperature and time discussed there. The initial sample is prepared by first simulating the system at a high temperature in a liquid state and then quenching the liquid to zero temperature. We used two different preparation protocols to obtain soft and hard glasses. To obtain a soft glass, starting from T = 1 we equilibrated the system for t = 20 simulation time units and then reduced the temperature to T = 0.1 and equilibrated for another t = 50. To obtain a hard glass, starting from T = 1 we cooled the system to T = 0.1 in steps of ΔT = 0.025, where at each step the system was equilibrated for t = 10. Once an initial solid sample was prepared, it was sheared quasistatically using a standard AQS protocol: at each strain step, the system is sheared using the Lees–Edwards boundary conditions⁶⁶ such that the total strain increases by 10⁻⁴. Immediately after strain is applied, the energy is minimized using the FIRE minimization algorithm.⁶⁷ 

We consider a scalar 2D lattice-based mesoscale elastoplastic model. The physics of this class of models relies on the coupling between a threshold dynamics and an elastic interaction induced by the incremental local plastic slip that arises as a result of a mechanical instability.³ 

More specifically, we consider a square grid of N × N cells of size a × a. The model is scalar, so that we account for one and only one shear direction, along which we can shear the system forward and backward. We assume a uniform shear modulus μ. Each individual cell (i, j) is characterized by a stack of local elastic branches indexed by a variable ℓ, each of which relates the local stress σ_ij to the local strain ɛ_ij, as shown in Fig. 1. The stability of each such local elastic branch ℓ is limited by two bounds: a maximum stress threshold $σij,ℓ+$ and a minimum stress threshold $−σij,ℓ−$⁠. Note that in order to ease notation, whenever no explicit reference to a particular branch number ℓ is made, we will omit it in the following. The two thresholds $σij+$ and $σij−$ are drawn from a random distribution with support in R⁺ so as to ensure $−σij−<σij+$⁠, i.e., the existence of a stability domain for the cell (i, j).

In the present model, the local stress σ_ij experienced by the cell (i, j) originates from two distinct contributions: a global stress Σ due to the external loading, and an internal stress associated with the interactions with other cells, so that $σij=Σ+σijint$⁠. The latter contribution fluctuates spatially and is by definition of zero average so that we have $σijint̄=0$⁠, and therefore $σij̄=Σ$⁠. Here, $Ā$ denotes the spatial average of the observable A.

Due to the external loading and the stress interactions, the local stress σ_ij is in general non-zero so that the amount of (external) stress that needs to be applied in order to reach one of the boundaries of the elastic branch is not a priori equal to the stress thresholds $σij+,σij−$⁠. Instead, it is given by the local plastic strengths in the positive and negative directions, which we define as $xij+=σij+−σij$ and $xij−=σij+σij−$⁠, respectively. Note that for a mechanically stable configuration we require that $−σij−<σij<σij+$ so that the quantities $xij±$ must be positive in that case.

The separation Δɛ between two neighboring local elastic branches that belong to a given cell (i, j) defines the local plastic strain $εij,ℓpl$ experienced by the cell after the local stress has reached threshold in one or the other direction.

Stress interaction—Local plastic strains are generated within an elastic matrix (the other cells of the lattice). This incompatibility induces an internal Eshelby stress field of quadrupolar symmetry.⁵² Since we assume homogeneous elasticity, the elastic response to a unit plastic slip can be computed once and for all. The internal stress thus directly arises from the convolution of the field of plastic strain with the Green function of Eshelby stresses. The latter is computed from the discrete Fourier transform of the analytical solution in the reciprocal space. Details on the implementation and a discussion can be found in Refs. 49 and 68.

The typical stress drop associated with a rearrangement of plastic strain Δɛ is of order μΔɛ. For the sake of comparisons with atomistic simulations, we consider here μ = 10, a typical value observed in Lennard-Jones binary model glasses.^19,68

Random landscape—The stress thresholds are drawn from a random distribution P(σ^±). Here, we consider a Weibull distribution of parameters λ = 1.0, k = 2.0, where λ and k are constants in the cumulative density function given by $1−e−(σ±/λ)k$⁠. The plastic strain increment $Δε=εij,ℓ+1pl−εij,ℓpl$ between two neighbor elastic branches ℓ and ℓ + 1 is also a random variable, cf. Fig. 1. We choose it to be correlated with the two plastic thresholds associated with the transition ℓ → ℓ + 1, i.e., $σij,ℓ+$ in the forward direction and $σij,ℓ+1−$ in the backward direction. More specifically, we choose Δɛ from a uniform distribution in [0, Δɛ_max] with $Δεmax=η(σij,ℓ++σij,ℓ+1−)/(2μ)$⁠, where η is a tunable parameter. Note that the parameter η thus controls the strength of the elastic interaction:^49,69 the larger η, the larger the short range stress kicks that trigger the avalanche, but also the larger the amplitude of the mechanical noise arising from the small positive and negative contributions of the long range stress interaction. We have set η = 1 in our simulations.

Nature of disorder—In the following, we will consider two different cases: (i) an annealed disorder where after a plastic slip new values of the thresholds $σij+$⁠, $σij−$ are computed in the absence of any memory; (ii) a quenched disorder, as a result of which the stress landscape of any given cell remains fixed so that the very same elastic branches are revisited in the course of a back and forth motion.

The landscape with quenched disorder is implemented through the use of a counter-based random number generator (CRNG)⁷⁰ so that the value of a threshold at the local elastic branch ℓ only depends on the index ℓ of that branch and on a previously defined key κ. In this way, the access to, say, $σℓ+=fκ(ℓ)$ requires just a simple call to the generator without the need of storing a full sequence of random numbers.

In the following, we will use an annealed disorder throughout the glass preparation step and a quenched disorder throughout the quasi-static shear driving steps. More specifically, we first “fabricate” our glasses using a two-step process, which mimics a thermalization step at high T and a subsequent aging step at vanishing temperature. We control the degree of aging of our glasses in this manner. Further details are given in Sec. II B. At the end of this preparation protocol, the different fields (thresholds in the forward/backward directions $σij±$ and internal stress σ_ij) are stored; the plastic strain field is reinitialized at zero and this initial configuration is inserted as the slice of index ℓ = 0 of a stack of quenched disorder thresholds at each cell (i, j). This quenched configuration is then used to perform mechanical loading.

Driving—Two kinds of mechanical loading are considered in this study: monotonous shear loading and cyclic loading. In both cases, the driving is strain controlled and changed quasi-statically. The elementary steps consist in (i) identifying the first site⁷¹ (i*, j*), which becomes unstable in the shear loading direction, i.e., the extremal site; (ii) incrementing the external strain ɛ up to the point where the extremal site (i*, j*) becomes unstable; (iii) incrementing the plastic strain of (i*, j*) by Δɛ to trigger the transition ℓ → ℓ ± 1 to the next elastic branch by the instability (plastification); (iv) updating the internal stresses of all sites; (v) identifying any site that has in turn become unstable due to the internal stress update, plastifying these sites as well, updating the internal stress, etc., until the end of the avalanche, i.e., until all sites have become stable again; (vi) repeat steps (i)–(v) as needed.

Avalanches—The precise treatment of step (v), i.e., the avalanche, deserves more detail. Once a list of unstable sites has been identified, the question remains about the order in which these sites will be updated. Indeed, since the elastic interaction can induce both positive and negative stress kicks, an unstable site can be healed and get stable again after another one has been plastified and the resulting internal stresses at the other sites have been updated, steps (iii) and (iv). Hence, the order of the updates matters. The effect of the ordering of updates on the dynamical properties has been recently discussed by Ferrero and Jagla.⁷² Some of us opted for a synchronous update:⁵¹ all unstable sites are plastified simultaneously in parallel; the internal stress is updated afterward; after this first sweep, a new configuration is reached, a stability test is performed, if all sites are stable, the avalanche is over, otherwise a new list of unstable sites is identified and the process is iterated until a stable configuration is reached. Here we make a different choice and perform a sequential update: the most unstable site, i.e., the extremal site, is updated first (plastic slip followed by an update of the associated elastic stress field) and we repeat this procedure until all sites become stable again. This choice of updating protocol happens to be very close to the extremal driving proposed in Ref. 55.

Although our goal is not to quantitatively map the elastoplastic model onto atomistic simulations, we must ensure that the disordered landscape statistics are comparable between the two types of models. From this point of view, being able to estimate the number of simulated elements of the mesoscopic model, i.e., the system’s number of degrees of freedom, is essential for a reasonable comparison that takes also into account possible finite-size effects.

To estimate the equivalent number of simulated elements, one must first determine the element size of the elastoplastic model below which the mechanical description is unresolved. This size corresponds to an upper limit of the characteristic plastic rearrangement size. Several experimental approaches have been performed to estimate the size of rearrangements ranging from direct observations in colloidal systems⁷³ to indirect estimations from strain rate sensitivity analysis in metallic glasses.⁷⁴ In all of these cases, the results show that plastic rearrangement cores contain a few dozen particles, so that the overall sizes of these cores range from about two to three particle diameters.

The determination of this length scale in atomistic simulations poses several difficulties. First, the presence of avalanches makes it challenging to identify the individual rearrangements. Second, there is no method yet to spatially distinguish between the non-linear and non-affine elastically strained zones from the non-reversible plastic responses. Finally, another complication arises from the fact that the same zone can contain several slip directions in a realistic particle system,^75,76 resulting in an effective higher density of potential rearrangements than that of a scalar description. Several approaches have been implemented to deal with these difficulties. They rely on the analysis of long-range elastic fields,^75,77 the quantitative calibration of elastoplastic models,⁵⁷ the calculation of the spatial extension of rearrangements,⁷⁸ the strain’s spatial correlations,^79,80 and the reproduction of the mechanical response from the spatial density of barriers.⁶⁸ These approaches, particularly those using a two-dimensional system under AQS loading like ours, lead to a consistent estimate of the linear size of plastic rearrangements lying between three and seven particle diameters. For our atomistic system containing 1024 atoms, these bounds lead to an equivalent mesoscopic system size between N = 5 and 10.

We conclude this section with a comparison of the SCC distributions obtained from our atomistic simulations with those obtained from catalogs of its mesoscopic size-equivalent with N = 8. The properties of our N = 8 mesoscopic catalogs are given in Table IV. The distributions are broad, and a prominent size cut-off at largest SCC sizes is now less prominent for the mesoscopic distributions when compared with the SCC size distributions obtained from the N = 32 mesoscopic catalogs, as shown in Fig. 10.

Figure 14 shows the dependence of stress under monotonous strain loading on system size and aging. Different colors correspond to different system sizes, as indicated in the legend, while the line shapes correspond to the different degrees of aging. The curves have been obtained at various extent of aging and for systems of size N = 16(4000), 32(1600), 64(750), where the numbers in parenthesis indicate the number of realizations used to obtain our results. While the poorly aged samples (PA with 0.8 aging steps per site) show no discernible size-dependence, with increasing amount of aging a rather weak system size dependence emerges, particularly near the stress peak, as shown in the inset.

In addition to Tables I and II in the main text, which describe the properties of catalogs extracted from moderately and poorly aged glasses of our mesoscopic model with N = 32, we also list here the properties of (i) two increasingly better-aged mesoscopic catalogs, prepared from glasses subjected to 150 and 4000 aging steps per site, which we will refer to as the well-aged (WA) and ultra-aged (UA) glasses, respectively, and (ii) atomistic catalogs obtained from 8 poorly aged and 30 well-aged reference configurations. The aging of the atomistic glasses is controlled by the rate of quenching to zero temperature from a high temperature liquid, as described in  Appendix A. We refer to these as fast (AS Fast) and slowly (AS Slow) quenched atomistic glasses, respectively.

Figure 15 shows an excerpt of the transition graph extracted from samples of sample No. 8 of our WA glass [panel (a)] and sample No. 3 of the ultra-stable UA glass [panel (b)]. The number of mesostates displayed in the t-graph excerpts shown are 1665 and 4610, respectively. We have obtained the graphs shown in Fig. 9 of the main text as well as in Fig. 15 by starting out in the reference configuration and following SCCs and the transitions between them until at least 1500 mesostates have been collected. For every SCC reached in this way, we also added the remaining mesostates belonging to that SCC so that the total number of vertices constituting the graph excerpt is typically larger than 1500. The number of SCCs shown in the excerpts of the two graphs in Fig. 15 are 216 (WA) and 19 (UA). The sizes of SCCs seen in the WA excerpt are small (s_SCC ≤ 56), while the UA excerpt has three very large SCC with sizes s_SCC = 3173, 1271, and 82, shown in pink, yellow, and green, respectively. These findings are consistent with the SCC scatter plots shown in panels (e) and (f) of Fig. 11. We believe that the emergence of the giant SCCs in the ultra-stable sample is a finite-size effect.

The stress–strain curves of the well-aged samples under uniform shear exhibit large stress changes across the yielding transition. For the WA and UA samples shown in Fig. 15, the magnitude of these stress-jumps under shear in the forward and reverse directions are ΔΣ = 0.46, 0.45 for the WA glass and ΔΣ = 0.62, 0.70 for the ultra-stable UA glass. In the graphs shown in Fig. 15, we have highlighted transitions that involve stress-jumps with a magnitude of at least 0.1 by fat black (U-transition) and red arrows (D-transition). Despite the relatively low threshold value chosen for these jumps, only very few transitions in the two graphs shown experience large stress changes. Note that for both the WA and UA samples, the transitions involving the large stress-jumps under forward and reverse shear tend to partition the graph into two halves (at least to the resolution of the number of vertices shown). This effect is even more dramatic for the ultra-stable glass sample where the transitions with large stress jumps immediately leads to giant SCCs.

We discuss the spatial arrangement of SCCs along a strip-like region in the $ESCC−,ESCC+$⁠, which is clearly evident for the atomistic systems as well as the WA and UA mesoscopic samples in the SCC scatter plots of Fig. 11. The diagonal dashed line in each of the plots corresponds to

As one would expect, the larger the strain range $ΔESCC$ over which mesostates are trapped within an SCC, the larger the size of the SCC itself. This trend is clearly visible in all six panels of the plots. The smallest (and most numerous) SCCs are clustered around small values $ΔESCC$⁠, while particularly for the WA and UA samples there appears to be a value $ΔESCC=ΔEmax$ beyond which it is unlikely to find SCCs, except for the outlier SCCs that we have associated with mesostates not having experienced the stress-peak. A naive estimate for $ΔEmax$ can be made as follows. Denote by Σ_ss the steady-state yield stress reached under monotonous loading (cf. Fig. 4). Assuming that between −Σ_ss and Σ_ss the system responds purely elastically, we obtain the estimate $ΔEmax=Σss/μ$⁠. From Fig. 4, we find that for the mesoscopic samples Σ_ss ≈ 0.85, while for the atomistic samples Σ_ss ≈ 2.4. The dashed lines shown in Fig. 11 correspond to these choices, i.e., $ESCC+−ESCC−=Σss/μ$⁠.

As our cyclic shear simulations show, at strain amplitudes close to but below the irreversibility transition, cyclic response may eventually be attained, but after a long transient. In particular, for the better-aged MA and WA samples we find that with increasing system size the transition to irreversibility becomes abrupt, meaning that we either reach cyclic response after a few driving cycles (typically one or two cycles) or not at all, implying a rather sharp and possibly discontinuous transition from reversibility to irreversibility. Our simulations indicate that this transition becomes smoother when the system size is fixed and the samples are less aged or when we reduce the system size at fixed aging steps per site. Thus, for example, for system sizes N = 16 and N = 32, the moderately aged samples are able to attain limit-cycles [cf. Fig. 8(a)], and even at strain amplitudes that are well beyond the location of the stress peak under uniform shear, which is thought of as marking the onset of yielding. To demonstrate this, we have used red triangles to mark on the monotonous loading curves of Fig. 4, the strain amplitudes beyond which the probability of finding a cycle is less than 2%. For the MA and WA samples, there are located beyond the stress peak.

These observations are consistent with findings in recent work by one of us on periodically sheared 3d atomistic glass formers.⁸¹ There it was found that small samples that were moderately or well-aged exhibit cyclic response at amplitudes well beyond the value of the strain at the stress peak. As the size of the samples increases, a sharp irreversibility transition at the stress peak is recovered. We should note, however, that in Ref. 81 such behavior was found to be the case only for totally asymmetric shear protocols of the form 0 → ɛ_amp → 0 → ɛ_amp→ ⋯.

Summarizing all of these findings: (i) we think that for the better-aged samples, the cycles reached after relatively long transients and at amplitudes beyond the stress-peak, meaning that the cyclic driving must have passed through it at least once, are an artifact of the system’s finite size. (ii) Related to this, we also find that as the samples get increasingly better aged, finite-size effects are not necessarily characterized by long transients. A case in point is the response to cyclic shear for the N = 32 UA sample. Here, it turns out that a cyclic response after only a few cycles is reached even at strain amplitudes as large as 0.4, which is well beyond the location of the stress peak. Consequently, such cycles are composed of a large number of plastic events, and hence mesostates, meaning that the SCCs confining them must be quite large as well. This finding is consistent with the “giant” SCCs already visible in the transition graph excerpt of the N = 32 UA sample, as shown in Fig. 15(b). A better understanding of such finite size effects, particularly in the ultra-stable glasses is clearly desirable and the subject of future work.