Multiperiodic orbits from interacting soft spots in cyclically sheared amorphous solids

Small groups of hysteretic elements can explain how a system returns to a previous state only after multiple driving cycles.


INTRODUCTION
A solid with perfectly elastic behavior deforms reversibly, in the sense that all material points return to their initial positions when a load is removed. Some amorphous solids may be prepared in a reversible plastic state, wherein loading the material in one direction changes its structure through many microscopic events, but loading it in the reverse direction precisely undoes these changes (1)(2)(3)(4)(5). Each microscopic event is localized to a soft spot (6) or shear-transformation zone (Fig. 1A) (7), which resembles a two-level system that switches under forward and reverse shear (4,(7)(8)(9).
Recent simulations using athermal quasi-static shear have revealed an even more remarkable behavior in which the period of particle motions is a multiple of the period of driving (1,10), reminiscent of the familiar action of a retractable pen. Such "multiperiodic" behavior may sound quite tenuous, given the daunting number of mechanically stable configurations and transitions in a packing of even a modest size. Nevertheless, multiperiodicity has been observed in molecular dynamics simulations of amorphous solids in two and three dimensions for several kinds of particle interactions (1,(10)(11)(12)(13)(14)(15). However, the mechanism for this behavior has remained unclear, even as it seems to be associated with an unjamming transition as the confining pressure is decreased (13).
Here, we show how multiperiodicity can arise in a simplified coarse-grained model of interacting soft spots (Fig. 2). We identify how the prevalence of multiperiodicity depends on the spatial arrangement of the soft spots, and we show how to design the behavior on demand. In all cases, the multiperiodic orbits are made possible by frustrated interactions in our model. Our results show that frustrated interactions between soft spots must be considered as an important counterpart to the cooperative interactions that are used to explain avalanches near the yielding transition (5,16,17).

While experiments
have not yet observed multiperiodic behavior, they exhibit the microscopic phenomenology we wish to distill into our model. Figure 1A shows a displacement field from an experiment with two nearby soft spots (see Materials and Methods for details). Each has the characteristics of an Eshelby inclusion-a small region of plastic deformation that is coupled to a quadrupolar elastic deformation of the surrounding material. This extended deformation induces or inhibits the rearrangement of other nearby soft spots, depending on their relative placement (8,16,18). For example, Fig. 1A is suggestive of a frustrated interaction, whereas the arrangement in Fig. 1B suggests cooperative interactions.
Our jumping-off point is to consider the possible behaviors of compact collections of N soft spots by modeling them as interacting hysteretic elements or "hysterons" (8,9). A hysteron has two possible states, s i = ±1; it transitions to the "+" state when the local field-equal to the instantaneous global strain field H plus neighbor interactionsreaches a fixed threshold H i + . Likewise, it transitions to the "−" state at a fixed threshold H i − < H i + . To model the disorder of such packings, these thresholds are set as H i where h i is chosen with uniform probability from the interval [ −1,1], and u i is chosen from [0, 2], for each hysteron independently. Hysteron j imposes a local field on hysteron i equal to J ij s j , where the coupling strength J ij is taken to be symmetric (J ij = J ji ) except where stated otherwise. The magnitude of each J ij (with i ≠ j) is selected with uniform probability so that |J ij | ≤ 1.
To capture the effect of the characteristic quadrupolar elastic deformations of rearranging soft spots, the signs of the J ij are dictated by the spatial configuration of the hysterons. Pairs that are 45° off the shear direction have a frustrated coupling (antiferromagnetic, J ij , J ji < 0), whereas pairs along 0° or 90° have a cooperative coupling (ferromagnetic, J ij , J ji > 0). This rule assumes that all soft spots' displacement fields have approximately the same orientation and polarity relative to the direction of shear, which appears to be true broadly in experiments (2,4,9,18).
Our simulations, available as an open-source Python package (19), probe the system evolution under athermal, quasi-static, oscillatory driving between −H 0 and +H 0 . We initialize the system with H ≪ −1 and all hysterons negative (s i = −1), and we evolve forward using an event-based method. Since flipping one hysteron may prompt a neighbor to flip, we wait for avalanches at fixed field until a stable state is reached; the hysteron farthest past its threshold is flipped first and all the local fields are updated between flips. In extremely rare cases where no stable state can be found or two flips are degenerate, the system is discarded. We continue driving until an absorbing state is reached where the dynamics repeat under further driving.
To search for multiperiodic behavior efficiently given the couplings J ij and thresholds H i ± , we note that increasing the driving amplitude H 0 will not change the dynamics until it is large enough to cause an additional hysteron to flip. Therefore, a finite set of H 0 will exhaust all possible dynamics under symmetric driving. To obtain this set, for each of the 2 N possible states, we compute the two values of H that bound the interval of stability for the state. We then sort the list of absolute values of these H and take the midpoints between successive values as our set of H 0 . We perform a series of simulations starting with the smallest H 0 and continuing until any multiperiodic orbit is found. Such an "amplitude sweep" is likewise an efficient method to search for unfamiliar behavior in experiments. Figure 2 shows an example of a multiperiodic orbit that is achieved for N = 4 hysterons arranged on a square. The system cycles through eight states over two driving periods, repeating this sequence indefinitely thereafter. This is just one possible T = 2 orbit for this spatial arrangement of N = 4 hysterons; it occurs with probability P = 8.37 × 10 −6 (allowing permutation of hysterons and inversion of the H i ± ). Figure 3 shows the prevalence of multiperiodicity in this and other compact arrangements of hysterons. The arrangements labeled a to e show all the unique configurations where N = 4 hysterons are placed within a 2 × 3 lattice that is oriented with the shear direction (up to reflections and rotations by 90°, which do not change the interactions). As before, interactions are between all nearest-neighbor pairs. Arrangement c has the highest probability of T = 2 among this set. These arrangements are some of the simplest ones eliciting multiperiodicity in our model.

Comparing arrangements of hysterons
Arrangements f to h in Fig. 3 show the increasing prevalence of multiperiodicity for N = 6, 8, and 9 hysterons on a square lattice. Arrangement h is multiperiodic with P = 5.3 × 10 −3 , so that if a macroscopic amorphous solid has 20 of these configurations, it will have a ∼10% chance of multiperiodicity. Notably, in contrast to the observed behavior of amorphous systems of many particles (1,3,13,20), small clusters of soft spots reach periodic orbits after very few cycles: For arrangement h, despite the space of 2 9 states, the longest observed transient before a (multiperiodic) limit cycle was just 3 cycles, and it occurred in just 1 of 10 7 systems.
When the lattice is rotated by 45° (exchanging cooperative and frustrated interactions, i.e., J ij → −J ij ), no multiperiodic orbits are observed (arrangements a′ to h′). This curious observation leads us to note another special property of a′ to h′: If we assign a + or − state to any one hysteron, we can then work outward and assign states to all other hysterons, satisfying every interaction. This is because these arrangements are portions of an antiferromagnetic lattice, with ordered ground states. While it is unclear why this property might suppress multiperiodic behavior, it could be a starting point for a deeper understanding of multiperiodicity generally.
The above results demonstrate that multiperiodic orbits can arise in our simple model constructed from coupled hysterons. In the following sections, we identify which attributes of the model are necessary for producing multiperiodicity.

Minimal number of hysterons
Empirically, we find that multiperiodic behavior is impossible for N < 3 hysterons. N = 3 hysterons with symmetric couplings also do not exhibit multiperiodic behavior. However, breaking the symmetry of at least one interaction pair (J ij ≠ J ji ) is enough to allow a T = 3 orbit, if and only if all interactions are frustrated. Under these conditions, we observe T = 3 with P = 4.67 × 10 −3 , with a single unique sequence of states (see the Supplementary Materials). We observe T = 2 with P = 7.80 × 10 −3 , accounting for a variety of different sequences.
Asymmetric couplings in spin systems without external cyclic driving (21-23) have been studied before, but the physical meaning in a driven amorphous solid is unclear. One possible mechanism  might be for soft spots to change states on different time scales, so that when the system is driven at finite frequency, a "slow" hysteron could fail to change in part of the cycle, even when in strict terms it is unstable.

Role of frustration
The observation that all interactions must be frustrated to elicit multiperiodic behavior for N = 3 prompts us to further investigate the role of frustration. In Fig. 4, we vary the fraction of interaction pairs that are randomly chosen to be frustrated (J ij , J ji < 0), and we plot the prevalence of multiperiodicity under these conditions. There is a clear trend across all the data: Multiperiodic behavior becomes exponentially more scarce as the fraction of frustrated pairs is reduced from 2 / 3 down to 0. In all cases, the probability is identically zero in the absence of frustration, a result we have checked up to N = 7. Figure 4 also confirms that the topology of frustrated and cooperative interactions can be just as important as their number: Arrangement a′ in Fig. 3 has N = 4 and 2 / 3 of pairs frustrated, and yet we find no multiperiodic orbits for that specific topology for either symmetric or asymmetric interactions.

Multiperiodicity from nonhysteretic elements
The above results show how coupled hysterons can produce multiperiodic orbits. We now show that hysteresis of the elements is not a necessary ingredient for multiperiodicity. In the absence of hysteresis and when −1 < J ij = J ji < 1, our model of an amorphous solid reduces to a spin glass where each soft spot corresponds to an Ising spin, governed by the Hamiltonian We verified this by writing separate code for such a spin glass and comparing the results with our coupled hysteron code with zero hysteresis. Deutsch and Narayan (24) reported multiperiodic orbits in such spin glasses with as few as five spins, although they focused on larger systems (N ≥ 64). We now elucidate the conditions for multiperiodicity with N = 5, under additional conditions that simplify the interactions even further: All the spin couplings are antiferromagnetic (J ij ≤ 0), and one or more of the couplings are randomly set to zero.
With four couplings set to zero, no multiperiodic orbits were observed in 10 6 systems. With 3 couplings set to 0, of 10 7 systems, we observe multiperiodicity in 1932-all with period T = 3 and a unique topology of interactions. This topology is shown in Fig. 5A and in the inset to Fig. 5B as a portion of a triangular lattice. Without loss of generality, we break the mirror symmetry by requiring |J 34 | < |J 01 | when spins are indexed left to right. This leads to  an additional remarkable uniqueness: At the smallest H 0 for multiperiodicity in each system, there is a unique and highly symmetric steady-state orbit (see the Supplementary Materials).

Regions in parameter space
The evolution of this spin glass model is deterministic given the coupling strengths J ij , an initial condition, and a driving protocol. Working in the reverse direction, a sequence of states may be mapped back to a region of the (high-dimensional) space of J ij that can give this sequence; here, a subset of the unit hypercube [ −1,0] n , where n = 7 is the number of nonzero couplings. Proceeding in this manner, we find a set of 10 inequalities among the J ij that bound the region of parameter space corresponding to this T = 3 orbit, which we list in the Supplementary Materials. The volume of this high-dimensional polygon (i.e., polytope) is found to be 1.86 × 10 −4 . To convert to a probability for multiperiodicity, we multiply this volume by 2 for the indexing degeneracy we lifted, and by 1 / 2 to account for the probability of obtaining the correct network topology. The latter factor may be found by noting that at least two of the three removed edges must share a vertex and enumerating the remaining cases. Thus, the above value is precisely the predicted probability of multiperiodic behavior for some H 0 . It agrees with how often we observe T = 3 in our simulations: P = (1.93 ± 0.04) × 10 −4 .
Such a detailed characterization of the high-dimensional phase space of the J ij is useful for designing systems with robust multiperiodic behavior. For instance, we can compute (25) a Chebyshev center for this polytope-a point that is farthest from its faces, which can thus withstand the largest possible errors in J ij while remaining multiperiodic. We find a Chebyshev center at which is the center of a hypersphere of radius 0.074 that lies entirely within the polytope. We report these coordinates to illustrate that our method can give precise quantitative information about finite regions of phase space that share a common orbit. At these coordinates, T = 3 is attained for any H 0 in the range 1 < H 0 < 1.685 (see the Supplementary Materials). For the general case of normally distributed errors in J ij , Fig. 5B shows that the probability of T = 3 remains high for an SD  up to several hundredths. Thus, the low probability of multiperiodicity in this system stems from the enormity of the parameter space rather than a need for fine-tuning. This same methodology-starting from an orbit and working backward to a region of parameter space-also applies to our model of interacting hysteretic soft spots. For example, setting H 0 = 1, a Chebyshev center for the orbit in which is a distance 0.338 from the nearest face. Note the high degree of symmetry at this Chebyshev center: Each hysteron has identical H + and H − , with identical asymmetric couplings that set up a clear chirality in the system. In the Supplementary Materials, we further characterize all of the above polytopes and list the inequalities that bound them.

DISCUSSION
We have shown how multiperiodicity can arise from the interactions of a small number of localized soft spots with simple, physically motivated interactions. Previous studies of this behavior using molecular dynamics simulations did not consider localization to soft spots (1,(10)(11)(12)(13)(14)(15), while previous attempts to understand it using simplified models (26-28) did not pursue a microscopic picture of the system, e.g., of the sequence or spatial structure of rearrangements. In this work, by focusing on small systems, probing the effect of the spatial structure of the elements, and using an amplitude sweep for the driving field, we have provided a concrete and thorough foundation for addressing the origin of multiperiodicity in amorphous solids, where its robust appearance in simulations has not been well understood. While we have not specialized to particular distributions of parameter values that are a subject of current research (17,29), our general model is nevertheless able to capture a detail of the multiperiodicity found in molecular dynamics simulations: We observe an approximately exponential decay of probability with the period of the limit cycle, T (e.g., in arrangements f to h in Fig. 3), a trend that was reported by Lavrentovich et al. (13) in simulations on jammed solids. The findings of Lavrentovich et al. that multiperiodicity may be associated with an unjamming transition prompt the question of the role of soft spot interactions in this critical transition.
Our results show that frustrated interactions are always necessary for multiperiodic behavior. This comports with existing theory about the random-field Ising model (30), which showed that without frustration, it supports return-point memory-a behavior that precludes a multiperiodic response. These findings suggest that multiperiodicity should be taken as a conspicuous signature of frustration-a counterpart to the yielding and shear-banding behaviors that are often attributed to cooperative interactions (5,16,17).
Our results also offer guidance to experiments searching for multiperiodicity in amorphous solids. Because interactions among soft spots are crucial, strain amplitudes should be large enough to ensure a high density of switching soft spots but small enough to allow a periodic steady state-consistent with results of prior simulations that we can now rationalize with our model. Such experiments also promise to reveal the role of soft spot interactions near yielding (8,17) and to probe the limits of the return-point memory behavior that is incompatible with frustration (8,9,(30)(31)(32). However, experiments must overcome measurement error and a high susceptibility to mechanical noise in this regime (4,9). We have shown that relatively few soft spots are sufficient for multiperiodic behavior, so that localized clusters of soft spots may be the dominant way that multiperiodicity emerges in large systems. Dividing observations of a large experimental system into regions of (10) soft spots could thus enhance sensitivity to multiperiodic orbits while rejecting the effects of mechanical noise or initial conditions playing out elsewhere. Furthermore, it would test the hypothesis that multiperiodic behavior is highly localized rather than being a strictly emergent behavior spread out among many interacting particles. Combinations of small groups with incommensurate periods may be a way for longer-period orbits to arise.
We have also shown that specific multiperiodic behaviors among spins and hysterons correspond to convex regions in high-dimensional parameter space, bounded by systems of inequalities. This both serves as an additional check of our modeling and paves the way for the rational design of systems with these behaviors, for example, as the basis for a digital counter. Most promising are the N = 5 spin configuration (Fig. 5, inset) and the N = 3 and N = 4 hysteron configurations (Figs. 2 to 4), each of which is conducive to a real-space physical implementation, with network topology and bond strengths that might be realized in the lab. Fig. 1 The experimental particle trajectories and micrograph used for Fig. 1 were obtained using methods described in (9), by cyclically shearing a monolayer of bidisperse polystyrene particles adsorbed at an oil-water interface. Because these particles exhibit long-range electrostatic repulsion (33), the material is a disordered, frictionless soft solid. We shear each sample between parallel boundaries that are 1.5 mm apart and 18 mm long; the material extends far beyond the open ends of this working sample. We image an approximately 1.4 × 1.9 mm region within the working sample. In fig. S1, we show micrographs corresponding to Fig. 1(A and B).

Details for
Each material is prepared by combining small and large sulfate latex microspheres (Invitrogen) in suspension, in roughly equal number, and dispersing them at an oil-water interface (9). The small and large particles in Fig. 1A have average diameters of 3.8 and 5.2 m, and the particles in Fig. 1(B and C) have average diameters of 3.5 and 5.4 m, although it is ultimately each particle's electric dipole strength that determines its effective size in the packing (33). Although aggregates of several particles can form during the preparation process, they do not seem to be strongly correlated with the locations of soft spots.
We obtained the plotted displacements (Fig. 1, A and B) by comparing the position of each particle at two different times and subtracting the average motion of the region of surrounding material with radius 16.5a, where a is the mode of the interparticle distance, determined from the pair correlation function g(r) (4,9,34). We chose times when the shear strain  = 0 (the midpoint of shearing), one full cycle apart. Figure 1A shows displacements in a portion of the system upon switching from strain amplitude 0.038 to 0.055. Figure 1B shows displacements in a different experiment upon switching from strain amplitude 0.045 to 0.050.
To obtain the plotted trajectory loops (Fig. 1C), we used positions over a full cycle of shearing at strain amplitude 0.055. Rather than subtracting the average motion within the region shown, we subtract the motion of a set of particles centered ∼35 m below this region, so that the particles in the field of view appear to be displaced horizontally by the global shearing motion.

Inequalities for regions of parameter space
For a system of spins, inequalities that bound regions of parameter space may contain only the parameters H 0 and J ij as variables. To generate such inequalities from a sequence of states, we follow a method that parallels our simulation algorithm. Two examples illustrate our approach. We first consider a spin i that flips to the + state as H is increased. At this instant, the spin has become marginally unstable, so that At this same instant, the other spins are stable, since otherwise they would have flipped before spin i did. For instance, if spin k ≠ i is in the − state, where we use the previous states s j of the spins before spin i flipped. Substituting Eq. 5 into Inequality 6 yields an inequality that contains only the unknowns J ij , as desired. As a second example, we imagine that the flipping of spin i causes another spin l to flip immediately (an avalanche). This tells us not only that spin l is unstable at the same value of H given by Eq. 5 but also that, at that instant, it is farther past its threshold of stability than every other spin. The avalanche ends when all spins are stable; this observation leads to further inequalities by again combining Eq. 5 and Inequality 6 (where Inequality 6 is flipped for spins in the + state). A similar method applies to a system of hysterons, with H i + and H i − as additional unknowns on the righthand side of Eq. 5 and Inequality 6 as needed. In general, additional inequalities are needed to denote that |H | ≤ H 0 at all times.
The resulting inequalities define a high-dimensional polygon (polytope). We use the Python package pycddlib (based on CddLib) to remove any redundant inequalities and the pypoman package (25) to compute Chebyshev centers. The full sets of inequalities for select orbits, further characterizations of the polytopes, and checks of the inequalities against our simulation results are given in the Supplementary Materials.

Period 3 polytope volume in the spin model
To compute the volume of the period 3 polytope for N = 5 spins, we first convert from a set of inequalities to a set of vertices, using the Python package pycddlib. We then compute the volume of the convex hull of these points with the SciPy module spatial.ConvexHull. We find it to be 1.863 × 10 −4 . Measuring this volume using Monte Carlo integration with 10 8 points gives consistent results: (1.857 ± 0.010) × 10 −4 . The full set of inequalities defining the polytope, and the 14 vertices they define, is given in the Supplementary Materials.

Organizing the hysteron simulation orbits
Comparing orbits lets us meaningfully group and count systems with equivalent orbits. We represent each simulation's output as a directed cyclic graph of states and manipulate it with the NetworkX package (35). We obtain the orbit by extracting the longest simple cycle in this graph. This removes trivial excursions: For instance, a system may transition from state + − + to + + + as H is increased and then return directly to + − + as H is decreased; we generally find many other randomly generated systems in which this excursion is missing. To compare these extracted orbits, we then account for all possible permutations of hysterons' identities, reversal of the sequence, and inversion of the system (exchanging all the + and − states).