When subject to cyclic forcing, amorphous solids can reach periodic, repetitive states, where the system behaves plastically, but the particles return to their initial positions after one or more forcing cycles, where the latter response is called multi-periodic. It is known that plasticity in amorphous materials is mediated by local rearrangements called “soft spots” or “shear transformation zones.” Experiments and simulations indicate that soft spots can be modeled as hysteretic two-state entities interacting via quadrupolar displacement fields generated when they switch states and that these interactions can give rise to multi-periodic behavior. However, how interactions facilitate multi-periodicity is unknown. Here, we show, using a model of random interacting two-state systems and molecular dynamics simulations, that multi-periodicity arises from oscillations in the magnitudes of the switching field of soft spots, which cause soft spots to be active during some forcing cycles and idle during others. We demonstrate that these oscillations result from cooperative effects facilitated by the frustrated interactions between the soft spots. The presence of such mechanisms has implications for manipulating memory in frustrated hysteretic systems.

## I. INTRODUCTION

Understanding the response of amorphous materials to time-varying mechanical loads is one of the main research topics in different fields, such as mechanical engineering, chemical engineering, materials science, geophysics, and condensed matter physics. Although the traditional notion of plasticity is that of irreversible deformation^{1–4} in recent years, it was discovered that under athermal quasi-static (AQS) conditions, plastic deformation can also be reversed.^{5–10} This form of deformation thus presents an intermediate type of deformation—*reversible plasticity*. One way in which reversible plasticity arises is as a result of oscillatory shear. When oscillatory shear of moderate amplitude is applied to an amorphous material, it experiences transient and irreversible dynamics until it responds periodically to the cyclic forcing so that the same microscopic configurations reappear periodically. Measured in units of the number of driving cycles, one denotes the length of the transient leading to cyclic behavior as *τ*, while the period of the response is *T*. In the case of *T* > 1, this type of periodic response is called *multi-periodic* (in the case of driven spin-glasses, this type of response has also been called *subharmonic*^{11}). It is known that as the amplitude of an applied oscillatory shear approaches a critical value, the transients get increasingly longer, while multi-periodic cycles start to appear. Once the critical amplitude is reached, a transition to diffusive behavior occurs, and a periodic response is no longer attainable. This transition has been called the *irreversibility transition*.^{5–7,9,12–18} Understanding the mechanism leading to multi-periodic response is of interest, as it emerges at strain amplitudes close to the yielding transition.^{5,8,10,19–21} At the same time, a *T* > 1 response is also of interest due to the possibility of encoding mechanical memory.^{22–26} In this work, we focus on the emergence of *T* = 2 cycles that turned out to be easier to analyze than cycles of larger periods.

To uncover the origin of multi-periodicity, we first use a variant of a recently introduced model of interacting hysteretic elements (hysterons), where the interactions are frustrated, which was shown to exhibit multi-periodicity.^{22–24} Contrary to previous studies, we focus on large systems and show that multi-periodicity results from oscillations of the switching strains of the hysterons. For a system with a large number of hysterons, these oscillations are a result of complex cooperative effects that involve several hysterons.

Next, by identifying the switching strain threshold of soft spots in molecular dynamics simulations of amorphous solids, we show that the switching strain oscillations are also responsible for multi-periodicity in amorphous solids. We further provide evidence for the existence of cooperative effects similar to the ones that we observed in the hysteron model. The relation between cooperative effects and multi-periodicity may explain why multi-periodic cycles emerge in amorphous solids at strain amplitudes close to the irreversibility transition, where the number of active soft spots becomes large.

## II. *T* = 2 IN A SYSTEM OF HYSTERONS

Experiments and simulations have shown that plasticity in amorphous materials is mediated by local rearrangement events. The regions where these rearrangements occur are called “soft spots” or “shear transformation zones.”^{27–31} An isolated soft spot has hysteresis properties: when sheared in the positive direction, its constituent particles undergo a localized rearrangement at a strain *γ*^{+}, and when subsequently strained in the opposite direction, this rearrangement is reverted at some strain *γ*^{−} < *γ*^{+}.^{32–36} Due to the long-range elastic deformations accompanying these state changes, soft spots interact with each other.

One can, therefore, model an amorphous solid as a set of interacting *hysterons*, where a hysteron corresponds to a soft spot with maximal and minimal strain thresholds. The *interacting Preisach model*, which was originally introduced by Hovorka and Friedman,^{37} and recently revisited within the context of memory formation in driven mechanical systems,^{22–25} demonstrated some important consequences of the interactions. In the following, we will refer to models of interacting hysterons as *iPreisach* models. It was shown that interactions in which the coupling can have both positive (“ferromagnetic”) and negative values (“anti-ferromagnetic”), and are thus frustrated, can induce switching order changes—scrambling,^{23} changes in the stability of hysterons,^{24} and enable the emergence of multi-periodic cycles.^{22}

To understand the subtle differences between a *T* = 1 cycle and a *T* > 1 cycle, we use a version of the iPreisach model with random interactions, similar to the model presented in Refs. 22 and 24 (see Appendix A for details). We consider a system of *N* = 100 hysterons that is driven by cyclic shear at some strain amplitude Γ_{max}. The strain thresholds $\Gamma i\xb1$ of the hysterons are given as

where $\gamma i\xb1$ are the strain thresholds without interactions, *ɛ*_{j} is a random strength factor, *s*_{j} = 0, 1 is the state of the *j*th hysteron, and *A*^{+} and *A*^{−} are the random, symmetric interaction matrices that determine how the upper and lower switching fields, respectively, of one hysteron change as a result of the state change of other hysterons. For our iPreisach model, a phenomenon that was not reported before, as much as we are aware, is that in some cases, after a hysteron changes its state, it can happen that for one of the hysterons Γ^{+} becomes smaller than Γ^{−}. If we use the usual notion of a soft spot as a two-level system with two minima,^{28} we can interpret this as the merging of the two stable minima of the soft spot (see Appendix A for further details). In our algorithm, we assume that this means that the state of this hysteron is *frozen* and cannot switch until Γ^{+} becomes larger than Γ^{−} again; only then the hysteron is unfrozen and regains its previous state.

To study the mechanisms leading to multi-periodicity, we start in a stable configuration with random *s*_{i} values (0 or 1) and a realization of the thresholds and interaction matrices. We focus on realizations that exhibit multi-periodic response *T* > 1 at some shear amplitude Γ_{max} and applied strain protocol −Γ_{max} ↗Γ_{max} ↘ − Γ_{max} ↗…, and let the system evolve until it reaches the multi-periodic cycle. Next, with the system in its zero-strain configuration **S** at the beginning of the cycle, we manually deactivate all the hysterons active in this cycle by fixing their states *s*_{i}. We then activate the hysterons one by one and in a random order. After each activation, we let the system evolve to attain a cyclic response. The cycles found after the activation of the first two hysterons must necessarily have *T* = 1,^{22–24} but at some point, the activation of another hysteron must lead to a response with *T* > 1. Each activation order yields a multi-periodic cycle with a different set of hysterons. Among those, we choose the cycle with the lowest number of hysterons for simplicity reasons. In the particular realization of the iPreisach model that we will consider here, multi-periodic response sets in for the first time after 29 hysterons have been activated. With the activation of hysteron 23, a multi-periodic cycle with *T* = 2 is attained. Adding 23 is, therefore, the *tipping point*, allowing the system of hysterons to switch from *T* = 1 to *T* = 2. Table I compares the switching patterns of the hysterons in the *T* = 1 and *T* = 2 cycles. We see that although hysteron 23 triggered the transition to the *T* = 2 cycle, once it has been activated it switches from 1 to 0, and then remains in state 0. Therefore, we can think of the triggering effect of hysteron 23 as a perturbation to the strain thresholds $\Gamma i\xb1,1T$ of the hysterons active in the *T* = 1 cycle. From Eq. (1), we find that the strain thresholds $\Gamma i\xb1,2T$ of the *i*th hysteron in the *T* = 2 cycle are given in terms of $\Gamma i\xb1,1T$ as

This shift in the strain thresholds causes four hysterons, 16, 93, 9, and 1, to have a switching pattern with period *T* = 2. To understand why, we perform next a careful analysis of their dynamics.

The *switching order* of the hysterons that are active in the *T* = 1 and *T* = 2 cycles, i.e., the order in which the individual hysterons change state during the different segments of these cycles, is shown in Table II. Multi-periodicity is manifested by the fact that some of the hysterons, in this case, 16, 93, 9, and 1, repeat their switching pattern after more than one forcing cycle. Note that hysteron 1 is the last hysteron to switch its state from 0 to 1 during the ↗ segment of the first driving cycle but also the first one to switch it back to 0 when the strain starts to decrease in the second driving cycle. Thus, its switching pattern does not affect the dynamics of the other hysterons. For that reason, and since hysteron 23, which caused the transition from *T* = 1 to *T* = 2 is not active during the latter, the set of hysterons active during the *T* = 1 and *T* = 2 cycles is effectively the same.

In Fig. 1(a), we schematically represent the dynamics of three hysterons; 93, 16, and 9, and show how they induce oscillations upon each other to form a multi-periodic cycle. We can see that 16 causes the activity of 93 to oscillate, 93 causes the activity of 9 to oscillate, and 9 causes the activity of 16 to oscillate, closing a loop. In Fig. 1(b), we see the switching strain thresholds of 16 and 9 during the first and second forcing cycles (we assert that the threshold is above Γ_{max} by continuing to increase the strain until the hysteron changes its state), and in Figs. 2(a) and 2(b), we can observe the value of their strain thresholds for every simulation step. From these figures, it is clear that the multi-periodicity of these hysterons arises due to oscillations of the threshold around Γ_{max}. When $\Gamma i+>\Gamma max$, the hysteron cannot switch, and when $\Gamma i+<\Gamma max$, it is able to switch (here, *i* = 16, 9).

In principle, we can think of a configuration where the threshold of all three soft spots oscillate around Γ_{max} and the switching or unswitching of each one causes Γ_{±} of one of the others to shift above or below Γ_{max}. In this scenario, all of the period-1 hysterons serve as a kind of “background” that does not affect the switching order of the period-2 hysterons. In Refs. 22 and 24, it was pointed out that a three-hysteron system is the minimal size that allows multi-periodicity, and a realization of a specially designed three-hysteron system that exhibits multi-periodicity due to strain threshold oscillations is demonstrated in Appendix B. However, as we discuss in the following, in a setting with many hysterons, the background hysterons play an important role, as they can cooperatively alter the switching behavior of other hysterons from one cycle to the next.

In the case of 16, it is clear why threshold oscillations are possible—the value of $\Gamma 16+$ is always close to Γ_{max} and small perturbations due to interactions can easily drive it above and below Γ_{max}. By inspecting Table II, we can also see that the shift of $\Gamma 16+$ above Γ_{max} is caused directly by the switching of 9. However, when we look carefully at the dynamics of 9, we can see that it does not follow this line. The value of $\Gamma 9+$ during activation is quite far from Γ_{max}, and the switching of 16 and 93 alone cannot cause it to oscillate around Γ_{max}.

The question is then how does $\Gamma 9+$ perform oscillations around Γ_{max}? By studying Fig. 2(b) and Table II, we can see that these oscillations are not caused directly by the switching of a single hysteron but by a mechanism that involves a change in the switching order of several hysterons, including background period-1 hysterons. In the first forcing cycle, hysteron 9 switches directly after 87 and before 93. However, because 93 did not switch in the ↘ segment, in the second forcing cycle, 69 and 21 switch before 41, which then causes $\Gamma 9+$ to become larger than Γ_{max}, and for that reason, hysteron 9 does not switch at all during the second forcing cycle. Therefore, the fact that 93 stayed frozen caused a cascade of events that indirectly causes $\Gamma 9+$ to become temporarily larger than Γ_{max}. We observe that the switching order of 69, 41, and 21 changes periodically between the first and second ↗ segments. In other words, the switching order of these three hysterons has a period of *T* = 2. The different state of hysteron 93, therefore, scrambles the switching order of these three hysterons, where in the first ↗ segment, these hysterons switch in the order of 41, 69, and 21, while in the second, the switching order is 69, 21, and 41. The combined effect that we call *cascade-scrambling*, thus, causes the threshold of hysteron 9 to oscillate around Γ_{max} with period *T* = 2. The multi-periodicity of hysteron 93 comes from an altogether different mechanism. In Figs. 1(c) and 2(c), we can see that it is a result of $\Gamma 93+$ becoming smaller than $\Gamma 93\u2212$ during the first forcing cycle, where Fig. 1(c) shows the range of values that they get during the first forcing cycle after 87 switches. As we discussed above, when $\Gamma 93+<\Gamma 93\u2212$, we define the soft spot to become *frozen*, i.e., unable to switch, as the hysteretic property vanishes. We assume that this corresponds to an elastic state of the soft spot, which is attained when its microscopic structure is distorted (see Appendix A for further detail).

Next, we studied the dynamics of the switching thresholds of soft spots in sheared amorphous solids and observed behavior resembling the hysteron model. This will be discussed in Sec. III.

## III. *T* = 2 IN MOLECULAR DYNAMICS SIMULATIONS

To study the response of an amorphous solid to cyclic shear, we use a two-dimensional molecular dynamics simulation. We simulated 1024 particles in two dimensions, interacting with a radially symmetric potential and integrated with a leap-frog solver.^{38} The potential that we use has a repulsive part identical to the Lennard-Jones potential and a minimum followed by a small hump.^{5,39} The main difference from the Lennard-Jones potential is the hump, which represents a potential barrier for breaking bonds, but it seems that this difference does not affect the existence of multi-periodic cycles observed in Lennard-Jones systems as well as in other types of potentials.^{10} In our simulations, 50% of the particles were 1.4 times larger than the other 50% to prevent crystallization. We prepared initial configurations by performing molecular dynamics simulations with an Andersen thermostat in which the system was kept at a high temperature *T* = 1 (in Lennard-Jones energy units) in a liquid state. The simulation ran for 20 simulation time units, after which the temperature was reduced to *T* = 0.1 and was run for another 50 simulation time units. The final configuration was then quenched to zero temperature using the fast inertial relaxation engine (FIRE) minimization algorithm that uses the molecular dynamics integrator to find a nearby local minimum.^{40} The result was a relatively soft initial configuration.

This initial configuration was then used to generate a state-transition graph (*t*-graph) that encodes the possible plastic deformations of the system, in a manner previously introduced by two of us^{32} (see Appendix C for details). In this graph, each node represents a stable configuration of the system, and black or gray (red or orange) arrows represent plastic events due to an increase (decrease) in the strain. The *t*-graph shown in Fig. 3 was extracted from a large graph obtained from molecular dynamics simulations. It corresponds to a *T* = 2 cycle with the lowest number of mesostates. The label “Start” marks the beginning of the driving cycle. To study the microscopic origin of the dynamics, we have created an algorithm that identifies all the soft spots in a cycle which has perfect or near-perfect loop return-point memory (*ℓ*RPM) property^{32,33,35,41–43} (see Appendix D for details). Multi-periodic cycles are more complex since they do not have the *ℓ*RPM property and thus using the topology is less straightforward. However, the mesostates belonging to such a cycle can be separated into clusters of *ℓ*RPM-like hierarchically organized cycles and sub-cycles. These clusters are sets of mesostates that are all part of one periodic cycle that obeys a near-perfect RPM, meaning that transitions out of the cycle or its sub-cycles are exclusively from the endpoints of the cycle only. Such a decomposition is shown in Fig. 3(b) where mesostates belonging to the same clusters are drawn as vertices of the same color. Table III depicts the switching patterns of the active soft spots. We can see that there are three period-1 soft spots—3, 4, and 6—that switch in all four segments, and seven period-2 soft spots—1, 2, 5, 7, 8, 9, and 10. Figure 3(c) shows the difference in Γ^{+} (black) and Γ^{−} (red) values of selected soft spots in the first and second driving cycles. We can see that for most of the soft spots that have period-2, e.g., soft spots 1, 5, 7, 8, and 9, the strain thresholds $\Gamma i\xb1$ oscillate around ±Γ_{max} in a manner similar to the threshold oscillations of hysterons 9 and 16 in the interacting hysteron example given in Sec. II. More specifically, for soft spot 1, $\Gamma 1+>\Gamma max$ in the first driving cycle but $\Gamma 1+<\Gamma max$ in the second driving cycle. Consequently, soft spot 1 switches only in the second driving cycle. Similarly, in the case of soft spots 5, 7, 8, and 9, their lower strain threshold $\Gamma i\u2212$ oscillates around −Γ_{max} between the first and the second driving cycles.

The dynamics of soft spots 2 and 10 are somewhat different. The measured values of $\Gamma 10+=\u22120.014$ and $\Gamma 2+=0.003$ are much smaller than Γ_{max} = 0.032, which indicates that it is not likely that they perform oscillations around Γ_{max}. Furthermore, when we increased the strain to values much larger than Γ_{max}, we could not find an activation threshold for these soft spots in the second forcing cycle. This indicates that both 2 and 10 are effectively inactive in one of the two forcing cycles. On the other hand, we observe that the measured value of $\Gamma 10+=\u22120.014$ is significantly lower than the measured value of $\Gamma 10\u2212=\u22120.0076$. This suggests that these soft spots may be frozen due to a *supercritical pitchfork bifurcation*, in which the two minima of an underlying effective potential merge into one, as explained in Appendix A. The effect of such a bifurcation would be to freeze hysterons 2 and 10 temporarily, in a manner similar to the freezing mechanism of hysteron 93 in the interacting hysteron example given in Sec. II.

In Fig. 3(a), we see that most of the period-2 soft spots are clustered together, which is consistent with the hypothesis by Keim and Paulsen^{22} that multi-periodic dynamics arises from interactions between localized clusters of strongly coupled soft spots. However, soft spot 5, which also has period-2, resides outside of this cluster and is thus less likely to be strongly coupled to other soft spots. On the other hand, when following the activation pattern of soft spot 5, and soft spots 3 and 4 that are activated before it, we can see that when 5 is idle, the activation pattern of 3 and 4 is scrambled: while 3 switches before 4 in the ↘ segment of the first cycle, 4 and 5 switch before 3 in the ↘ segment of the second cycle. This suggests that the period-2 of soft spot 5 may be a result of an indirect coupling with another period-2 soft spot via a cascade-scrambling mechanism.

Another mechanism observed in molecular dynamics simulations is the overlapping of soft spots. This is demonstrated in Fig. 3(a) as overlapping of the ellipses that represent soft spots. An overlap, in this case, means that two (or more) soft spots share particles. One implication is the existence of conditional switching patterns, meaning that one soft spot cannot switch before an overlapping one does. For example, we see that soft spots 7 and 9 overlap and, indeed, 7 switches before 9 in the ↘ segment of the first cycle and then switches after 9 in the ↗ segment of the second cycle. In the current simulation results, we did not find any indication that overlaps affect multi-periodicity. However, we expect that in a system with many soft spots, overlap between soft spots can cause some soft spots to become disabled, which may require modification of the iPreisach model.

## IV. DISCUSSION

In this work, we studied the mechanisms that lead to a multi-periodic response in sheared amorphous solids using a model of interacting hysteretic elements (hysterons) and molecular dynamics simulations of amorphous solids. In both models, we have shown that some hysterons (soft spots in the actual amorphous solid) exhibit multi-periodic responses during a multi-periodic cycle, while others remain uni-periodic and act as a “background.” We have shown that multi-periodicity of an individual hysteron results from strain threshold oscillations—a periodic change in the strains at which the hysteron switches, which causes the hysterons not to be able to switch during part of the cycle. The analysis of the interacting hysteron model reveals that multi-periodicity emerges when three or more hysterons interact in such a way that during alternating cycles, one or more of them are inactive but become activated in the next cycle. Careful analysis revealed three mechanisms that the alternating activity of a specific hysteron can arise from. In the first mechanism, the switching threshold of the hysteron oscillates around the forcing magnitude Γ_{max} due to direct interactions with other multi-periodic hysterons. These oscillations cause the hysteron to become active and inactive periodically, giving thereby rise to a *T* > 1 switching pattern. In this case, the background period-1 hysterons do not play an active role in the multi-periodicity of the hysteron. In the second mechanism, multi-periodicity is also caused by oscillations of the switching threshold around Γ_{max}. However, the change in the thresholds is not caused by direct interaction with a multi-periodic hysteron but is mediated by the period-1 hysterons that now play an active role in the multi-periodic dynamics. Oscillations are then formed when the switching order is scrambled in a way that prevents the cooperative effect from deactivating the hysteron. In the third mechanism, the hysteron becomes active and inactive or “frozen” periodically when the minima of the double-well potential merge into a single one, and its stability region thus becomes null. An analysis of the switching patterns of soft spots in a molecular dynamics simulation of a periodically sheared amorphous solid, using a specially developed algorithm, finds signatures of all three mechanisms, which suggests that they, indeed, play a role in the multi-periodic dynamics of actual amorphous solids.

The nature of the mechanisms that we uncovered has significant implications for the understanding, modeling, and uses of plasticity in amorphous solids. It has long been recognized that interactions between soft spots play an essential role in the plasticity of amorphous solids, but up until now, the focus has been on the role that interactions play in the size and duration of avalanches. The order in which these soft spots switch was not thought to play an important role in the dynamics. Here, we have shown that changes in the switching order of soft spots are significant for the dynamics under periodic forcing. Therefore, a theoretical or numerical model of the dynamics of amorphous solids under oscillatory shear should include switching order effects. In addition to the effects of soft-spot interactions on the order of switching, we have also found evidence suggesting that interactions may cause a soft spot to be temporarily inactive, which may also be important for understanding plasticity in amorphous solids. A more detailed model and the study of higher-period cycles might unveil additional mechanisms of multi-periodicity.

The mechanisms that we unveiled may be useful for manipulating memory encoding and readout in amorphous solids. Furthermore, since the interacting hysteron model is rather generic, our results may also be relevant to other mechanical systems exhibiting memory effects such as crumpled and corrugated sheets^{25,44} where in the latter scrambling effects have been observed.^{25} Therefore, it is possible that manipulating such effects can allow the creation and manipulation of cycles with periods larger than one also in these or related systems.

## ACKNOWLEDGMENTS

I.R. and A.S. were supported by the Israel Science Foundation through Grant No. 1301/17. M.M. was supported 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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX A: THE iPREISACH MODEL—SIMULATION DETAIL

The iPreisach model can be thought of as an abstract representation of a set of interacting soft spots in an amorphous solid. As we pointed out in the main text, each soft spot can be regarded as a two-level system, and its potential energy, thus, can be approximated by a double-well potential, i.e., a soft spin, as considered by Puglisi and Truskinovsky.^{45} Here, the two minima of the potential correspond to the states *s*_{i} = 0, 1 of the *i*th soft spot. In Fig. 4, we see an example of a double-well potential of the form

as well as the effect of changing the external driving Γ, where *x* is a coordinate representing the state of the soft spot, Λ and Δ are constants. Starting from a situation where *x* is at the left minimum and increasing Γ, the potential is distorted. When a critical value Γ^{+} is reached, the minimum becomes unstable (a saddle-node bifurcation), and the system transitions to the other minimum. Similarly, starting from the other minimum and decreasing Γ, at a critical value Γ^{−}, the minimum is switched again [Fig. 4(a)]. Λ and Δ can also change due to interactions between hysterons and they can thereby alter the limits of stability Γ^{±}, which then can also cause a bifurcation into a state in which the two minima merge into a single minimum (a supercritical pitchfork bifurcation^{46}). This new minimum thus represents an intermediate soft spot configuration that cannot switch [Fig. 4(b)]. We will refer to such hysterons as temporarily frozen since their response is purely elastic.

To model the conditions in an amorphous solid, we simulate a system of *N* interacting two-state hysteretic elements, hysterons, similar to the models presented in Refs. 22–24. As in the amorphous solid, the system is subject to a cyclically varying external field Γ with a given maximal strain amplitude Γ_{max}. Each hysteron is assigned a strength coefficient *ɛ*_{i} drawn independently from a uniform distribution $U[0.1,1.0]$, representing a variation in its properties. In the absence of interactions, the switching fields of each hysteron are given by $\gamma i+$ and $\gamma i\u2212$. We call $\gamma i\xb1$ the *bare* switching strains or strain thresholds. Once interactions are introduced, the thresholds become the *dressed* thresholds $\Gamma i\xb1$ and are given by Eq. (1). The matrices *A*^{±} are symmetric and their elements are drawn independently and randomly from a normal distribution with a zero mean $N(0,0.0064)$, which favors the weaker interaction values. In this way, we represent the stronger interaction between the fewer proximate soft spots and weaker interactions with the distant ones.

The initial state *s*_{i} = 0, 1 of each hysteron is assigned independently and randomly with equal probability. Choosing $\gamma i+>0$ and $\gamma i\u2212<0$, which are drawn independently from the uniform distributions $U[0,1]$ and $U[\u22121,0]$, respectively, guarantees that the initial configuration is mechanically stable at the start of the simulation run when Γ = 0. However, there are no limitations on the values of $\Gamma i\xb1$ once the system evolves. The dynamics is event-driven: in the increasing part of the drive cycle, at each step, the strain Γ is incremented to the value of the smallest Γ^{+} from the set of hysterons that are in the state *s*_{i} = 0. When Γ^{+} > Γ_{max}, we set Γ = Γ_{max}, switch the direction of the strain, and repeat the process where now we switch the strain to the largest Γ^{−} from the set of hysterons with *s*_{i} = 1. We repeat this until Γ^{−} < −Γ_{max} at which point we reverse the drive direction again and so on.

To take into account avalanches, after each update of a hysteron *i* and before proceeding to the next strain threshold, we check whether other hysterons *j* ≠ *i* become unstable, i.e., $\Gamma j+<\Gamma $ when *s*_{j} = 0 or $\Gamma j\u2212>\Gamma $ when *s*_{j} = 1, in which case we flip these hysterons as well in the same simulation step. In the molecular dynamics simulations, the dynamics during an avalanche is purely relaxational. To be consistent with such dynamics, we chose the switching order during an avalanche such that the hysteron that “overshot” Γ by the largest extent will switch first. The switching order is, thus, determined according to the *instability* $\Delta Ej$ of each hysteron, defined as

We thus switch the hysteron with the largest $\Delta Ej$. After every switch, we update the values of all $\Gamma i\xb1$ and repeat the process until all hysterons are stable and then we continue increasing (decreasing) Γ.

### APPENDIX B: STRAIN THRESHOLD IN A SIMPLE CONFIGURATION WITH THREE HYSTERONS

As was pointed out in Refs. 22 and 24, the smallest number of hysterons required to produce a multi-periodic cycle is three. Therefore, we used our iPreisach model to study a simple *T* = 2 cycle, generated using three interacting hysterons with strain thresholds set to be *γ*^{+} ≈ Γ_{max} and *γ*^{−} ≈ −Γ_{max} similar to the case of hysteron 16 in the main text. The results are summarized in Fig. 5, where the first driving cycle starts from a state **S** = (0, 0, 0). During this cycle, hysteron 3 switches two times (once in each segment), hysteron 1 switches once, and hysteron 2 does not switch at all. Therefore, at the end of the first driving cycle, the system is in the state **S** = (1, 0, 0). In the second driving cycle, hysteron 2 switches twice, hysteron 1 switches once, and hysteron 3 does not switch at all, which brings the system back to the **S** = (0, 0, 0) state. When looking at the dressed strain thresholds $\Gamma i\xb1$ of the three hysterons, we can see that all three are period-2 hysterons that exhibit the kind of threshold oscillations that were observed in hysterons 16 and 9 (Sec. II) and in soft spots 5, 7, 8, and 9 in the molecular dynamics simulations (Sec. III). We also see that each hysteron switches in different segments of the cycle and, thus, causes the threshold oscillations of the other hysterons. This means that the interactions cause a collective periodic behavior of out-of-phase threshold oscillations. More specifically, during the first driving cycle, $\Gamma 1+$ and $\Gamma 3+$ are smaller than the maximal applied strain Γ_{max}, which causes them to switch, while $\Gamma 2+$ is larger than Γ_{max} and for that reason, it does not switch during the first cycle. In this example, it is clear what is causing the oscillations: once *s*_{3} = 1, $\Gamma 1\u2212$ changes from its bare value $\gamma 1\u2212$ to its dressed value,

and for that reason, 1 remains in *s*_{1} = 1 until the second driving cycle. When 3 switches back to *s*_{3} = 0, the 1 negative threshold goes back to its bare value $\gamma 1\u2212$, which is larger than −Γ_{max} and, thus, can switch back to *s*_{1} = 0 during the second driving cycle. Similarly, when *s*_{3} = 1, the strain threshold of 2 is pushed to a value that is larger than Γ_{max} and, thus, *s*_{2} = 0 during the first driving cycle. But when *s*_{1} = 1 and *s*_{3} = 0, $\Gamma 2+$ is pushed to a value that is smaller than Γ_{max} and can therefore switch in the second driving cycle.

### APPENDIX C: CONSTRUCTING *t*-GRAPHS FROM MOLECULAR DYNAMICS SIMULATIONS

To construct the graph, we start with a stable configuration at some strain Γ (in this case, a state that is part of the periodic cycle of interest). We identify its strain thresholds, Γ^{+} and Γ^{−}, where a plastic event occurs, including one or more soft spots. All the configurations in the range Γ^{−} < Γ < Γ^{+} are considered part of a single mesostate *O*, represented as a vertex on the graph. The plastic events at Γ^{+}[*O*] and Γ^{−}[*O*] lead to new configurations that are part of mesostates *A* (that is reached after Γ^{+}[*O*]) and *B* (that is reached after Γ^{−}[*O*]) that have different strain thresholds Γ^{+}[*A*], Γ^{−}[*A*] and Γ^{+}[*B*], Γ^{−}[*B*], respectively. The mesostates *A* and *B* are represented as two separate vertices on the *t*-graph, and the plastic transition leading from *O* to *A* is represented as a black (or gray) arrow and the plastic transition from *O* to *B* is shown as a red (or orange) arrow. This process is repeated iteratively for each new mesostate *M*, where plastic transitions at Γ^{+}[*M*] and Γ^{−}[*M*] can lead to new mesostates or back to mesostates that were encountered before. When this happens, a graph loop is closed. Detailed explanations and examples are given in Ref. 32 and its supplementary material; here, we have followed the notation used in related studies.^{22–25,32,41,47,48}

### APPENDIX D: SOFT SPOT IDENTIFICATION METHODOLOGY

For periodic cycles that have the *ℓ*RPM property, there is a well-defined hierarchy of sub-cycles. The main cycle is composed out of a lower and an upper mesostate, the endpoints of the cycle, which we will denote by **S**_{zero} and **S**_{one}, respectively. Starting in the lower endpoint **S**_{zero}, the upper endpoint is reached by increasing the strain to some value Γ_{max}, while a subsequent decrease of strain to some value Γ_{min} brings the system back to **S**_{zero}. Thus, under repeated applications of Γ_{min} → Γ_{max} → Γ_{min}, the graph cycle with endpoints **S**_{zero} and **S**_{one} is traversed. The *ℓ*RPM property assures that if starting again from **S**_{zero}, but increasing the strain this time up to some value Γ < Γ_{max}, so that some state **S** is reached, then a sufficiently large decrease of the strain will always bring the system back to **S**_{zero}. Thus, **S**_{zero} and **S** are the endpoints of a sub-cycle. Similarly, starting from the upper endpoint **S**_{one} and decreasing the strain to some value Γ > Γ_{min}, a subsequent increase of strain will eventually return the system to **S**_{one}.

In terms of the soft spots, one can associate with a cycle a set of active soft spots that change their states as the cycle is traversed, so that at the lower endpoint **S**_{zero} of the cycle, all active soft spots are in their 0 state, **S**_{zero} = (0, 0, …, 0), while at the upper endpoint, they are in their 1 state, **S**_{one} = (1, 1, …, 1). Thus, when *ℓ*RPM is present in the corresponding sub-cycles, only a subset of these soft spots is active. Consequently, a given soft spot will be active in multiple sub-cycles and its switching behavior will depend on the states of the other active soft spots. By using the *t*-graph to decompose a cycle into its sub-cycle and combining this information with the displacement fields of the plastic events, we can identify the soft spots that are active in each of the transitions. Moreover, by appropriate book-keeping combined with a method to compare the displacement fields **v** of different plastic events, we can identify the soft spots participating in the transitions of interest. To identify whether two different displacement fields represent the same plastic event, we compare **v**, the 2*N*-vector of displacements during one plastic event (*N* is the number of particles),

to the vector **u** of the other plastic event. Here, $(vxi,vyi)$ is the displacement of the *i*th particle. To compare the fields **v** and **u**, we calculate the scalar product of the normalized displacement vectors $v\u0302$ and $u\u0302$:

If the plastic events involve the same set of soft spots, the result will be $\u223c1$ or ∼−1 (we used a threshold of 0.98), where the latter corresponds to cases when the compared fields are of the same soft spots in transitions of increasing and decreasing strain (in the case of a single soft spot, this is a case where we are comparing the displacement fields generated by the same soft spot *i* in a transition *s*_{i} = 0 → 1 to the displacement field generated by the same soft spot in the transition *s*_{i} = 1 → 0). In plastic events where more than one soft spot are switched—an avalanche—we assume that the fields can be linearly superimposed. This assumption can break down when the soft spots are too close to each other but still worked in most of the cases that we have studied. When this assumption is not satisfied, we have to intervene manually.

The identification of soft spots involved in transitions between clusters has to be handled separately. In this case, special care has to be given to transitions involving avalanches, since identifying all the individual soft spots of the avalanche may require searching for “nearby” loops, i.e., graph loops that can be reached by a few transitions from the cycle of interest, and in which they are active individually, that is, they are not part of an avalanche. An example where all soft spots involved in the relevant transitions of the *t*-graph have been identified is given in Fig. 3 of the main text.