Interplay of rearrangements, strain, and local structure during avalanche propagation

Jammed soft disks exhibit avalanches of particle rearrangements under quasistatic shear. We follow the avalanches using steepest descent to decompose them into individual localized rearrangements. We characterize the local structural environment of each particle by a machine-learned quantity, softness, designed to be highly correlated with rearrangements, and analyze the interplay between softness, rearrangements and strain. Our findings form the foundation of an augmented elastoplastic model that includes local structure.


I. INTRODUCTION
All disordered solids respond elastically at low strain but flow plastically at sufficiently high strain. As strain increases beyond the elastic regime, disordered solids partially relax via intermittent localized rearrangements until they reach the yield strain, where they begin to flow. Up to the yield strain, disordered solids display surprisingly universal behavior [1]. Beyond the yield strain, however, disordered solids exhibit several different classes of plastic behavior. Foams can flow indefinitely via localized rearrangements without ever fracturing [2]. Many systems exhibit crackling noise or avalanche behavior [3][4][5][6], while still others exhibit shear banding and brittle fracture [7]. Here we focus on avalanche behavior.
An avalanche consists of a series of rearrangements. A class of models known as elastoplastic models describes the courses of avalanches in terms of the interplay of rearrangements and elastic stress [8]: an increase of elastic stress can cause a local region to yield and rearrange, while conversely, a local rearrangement can increase stress elsewhere. It has become increasingly clear, however, that rearrangements and elasticity do not tell the whole story. Systems with identical microscopic interactions can show ductile or brittle behavior depending on preparation history [9,10]. This bolsters approaches that postulate local structures prone to rearrange [11,12], but also points to the need for microscopic understanding of the connection between local structure and the physics included in elasto-plasticity models. While it has been shown that certain local structural environments are much more likely to rearrange than others [13][14][15], the effects of rearrangements on local structure have not been established, even though it is clear that they must exist. It is also clear that elastic stresses can distort the structural environment surrounding a particle. These considerations point to the need for detailed understanding of the interplay of local structure, rearrangements and elasticity.
In this paper, we go back to basics to untangle the interplay of local structure, rearrangements and strain in athermal, quasistatically sheared jammed packings of * ajliu@physics.upenn.edu soft disks. We capture local structure via a machinelearned quantity, softness [1,[15][16][17]. This quantity has been shown to provide useful insight into the dynamics of supercooled liquids and glasses [15,16,18]. Following this earlier approach [15], we describe softness as the weighted sum of a set of structural quantities based on the local pair correlation function, where the weights are chosen to maximize the correlation with rearrangements that occur during avalanches. To tease out the interplay between softness, rearrangements and stress, we look at the effects of each of these factors on the others to develop a "structuro-elasto-plasticity" framework for avalanches in disordered solids.

II. SIMULATION DETAILS AND SOFTNESS
We generate two-dimensional packings of N soft disks in a simulation box with periodic boundary conditions. The disks interact with each other through the pairwise additive Hertzian potential: where σ i is the radius of the ith disk. To avoid crystallization, we use a 1 : 1 mixture of particles with σ = 0.5 and σ = 0.7. Starting from random initial conditions, we minimize the potential energy to find the initial zero-temperature jammed state. We then repeatedly apply a small shearstrain step of δ , minimizing the energy after each step, until the total strain reaches end . The stress-strain relation for a single configuration, shown in Fig. 1, confirms the existence of avalanches. We generated 5 trajectories with N = 10 5 , δ = 10 −5 , and end = 0.1; and 20 trajectories with N = 4000, δ = 10 −4 , and end = 2. This smaller system with N = 4000 is shown in Fig. 1 for visual clarity. It is also used to train the machine-learning algorithm because we need to access larger shear strains, as detailed in the supplementary information (SI) [19]. All of the remaining analysis was carried out on the larger system.
It is well known [20] that during athermal quasistatic shear, energy drops mark rearrangements that can be either localized or extended due to avalanches. In each During an avalanche ...

FIG. 1:
As strain increases, avalanches occur during stress drops. During an avalanche, some constituent particles rearrange, triggering other localized rearrangements far away in the depicted system of N = 4000 particles. Here, the non-affine displacement D 2 min of particles is represented on a black-to-blue-to-red scale with red corresponding to high values of D 2 min . The rightmost plot depicts the cumulative D 2 min measured over the entire stress drop.
step of strain followed by energy minimization, we calculate the final energy to monitor for energy drops. We focus only on energy drops, using steepest descent with adjustable step sizes [19] to capture the over-damped relaxation process from the beginning of the energy drop to the end. We save intermediate configurations that are equidistant in configuration space, so that the sum over particles of particle displacement squared is fixed between successive frames.
To identify rearranging particles, or "rearrangers," we calculate D 2 min [21]: where the sum is over all neighbors of particle k within a distance of R D = 2. Here M k is the number of such neighbors, r ik and r ik are the vector separations between particles i and k at two consecutive frames, respectively, and J k is the "best-fit" local deformation gradient tensor about particle k that minimizes D 2 min . A particle with D 2 min above a certain threshold [19] is a rearranger. The rest of the paper presents results for rearrangers that are small particles in our binary mixture, but we have verified that results for large-particle rearrangers are qualitatively the same.
We draw from the sets of rearrangers and nonrearrangers to train a linear support vector machine to obtain softness, a structural quantity that indicates the propensity of a particle to rearrange [15]. As detailed in the SI [19], the softness of particle i is essentially a weighted integral over the local pair correlation function g i (r). The weight function is inferred by the linear support vector machine to maximize the accuracy of predicting rearrangers. As in Ref. [15], the weighting is highly negative at the first peak of g(r), implying that particles with fewer neighbors have higher softness, consistent with intuition based on the cage picture.

III. THE AVALANCHE PROCESS
In Fig. 1 and the supplemental video [19], we confirm that during avalanches, rearrangements are indeed localized and sequential, as assumed in elastoplastic models [8]. Moreover, consecutive rearrangements can be very far apart.

A. Interplay of rearrangements with strain
We begin our study of the interplay between rearrangements, softness and strain by examining the effect of a rearrangement at the origin on strain at r, averaged over many rearrangements. The local-fit deformation tensor about each particle, J, is already known from calculating D 2 min . From J, which is a second-rank tensor with four degrees of freedom in 2D, we extract three different strain components: the volumetric (isotropic) strain k, total deviatoric strain˜ , and shear strain in the xy direction (the direction of the global shear), xy . We extract these by comparing two consecutive frames when a rearrangement is occurring (see SI).
The near-field behaviors of the local strains depend on microscopic details of how rearrangements locally strain their surroundings, but in the far field we expect the local strains to be well-described by elasticity theory. In the far field, one typically approximates the rearrangement as a point plastic shear strain, equivalent to a pair of point force dipoles. The responses to this source of the three local strains studied at position r and time t following a rearrangement at the origin at t = 0 are given in Eqs. (A5)-(A7). Because the actual shear strain source due to the rearrangement is very long-lived, however, the response to the point plastic shear strain is well-approximated by the infinite-time limit, shown in Eq. (A8).
Analytical derivations [22], numerical measurements [20], and experiments [23] have all found that xy has an r −2 radial dependence and a quadrupolar angular dependence in response to a point force dipole. Indeed our results show that | xy | decays as a power law close to r −2 with a quadrupolar angular dependence, as indicated by the bottom row of Fig. 2. The red solid lines are fits to the dependence expected from continuum elasticity (see Appendix A). Because of the quadrupolar angular dependence, the angular average of xy (r) is zero within statistical noise (almost every error bar crosses the x axis). The deviatoric strain,˜ (middle row in Fig. 2), likewise decays as r −2 (red solid line in left plot) but with an isotropic angular dependence (right plot), as expected from continuum elasticity (see Appendix A).
The volumetric strain k(r) is typically neglected in systems of fixed total volume but as we will show, it plays an important role because softness is strongly dependent on local density. The volumetric strain in response to a shear strain is given in Eq. (A5). The far field response to a rearrangement, however, must also take into account the effect of a point compression source since the rear-rangement can also give rise to local plastic compression. This has a transient effect since the total volume of the system is conserved, but is significant because it gives rise to a contribution to k(r) [Eq. (A9)] that does not angleaverage to zero. As a result, the volumetric strain k(r) is predicted by elasticity theory to be the sum of two terms: a sin(2θ)r −2 term arising from the point shear strain and another term arising from the point compression. The top left plot of Fig. 2 shows that the angular-averaged volumetric strain k(r) is positive at most r and does not exhibit a power-law decay; this corresponds to the second term. In the Appendix we discuss the expectation from elasticity theory that yields the fit (red solid curve) shown. The absolute value of k(r) is dominated by the first term and shows an r −2 decay, consistent with continuum elasticity theory for a point plastic shear strain (red solid line in top middle plot). For r > ∼ 5, k has the expected dipolar angular dependence from the first term (top right plot). In summary, of the two contributions to the local volumetric strain, the term arising from the point shear strain dominates but angle-averages to zero so that the second term is revealed in the angularaveraged k(r).
Although the results shown here are for twodimensional systems, we have confirmed that the expected scalings for volumetric and deviatoric strain are observed in 3 dimensions [19], providing strong evidence in favor of our interpretation of the roles of volumetric, deviatoric and xy-strain.
We now turn to the effect of the induced strain on the next rearrangement. We first compute the frame-dependent pair correlation function of rearrangers g 2 (r, δf ), namely the probability of finding a rearrange- and for rearrangers only (blue dotted). There is a pronounced difference between the two distributions. (b) The probability that a particle is rearranging, P R , as a function of its softness. As the softness increases, P R increases by four orders of magnitude, verifying the high correlation between softness and rearrangements. ment at r δf frames later, given a rearrangement at the origin at frame δf = 0. Results for several values of δf are presented in Fig. 3.
We first focus on the radial dependence. The rearranger pair correlation function g 2 (r, δf ) decays as r −3 , independent of δf . This is consistent with either˜ or xy , which both decay as r −2 , due to the following argument. Two earlier studies of systems with sphericallysymmetric potentials found that the cumulative distribution of the local yield strain has a low-yield-strain tail described by a power law with exponent 1.6 [24,25]. On general grounds this scaling should also apply to our system [24], so the probability that a rearrangement is trig- gered by˜ or xy ∼ r −2 should scale as (r −2 ) 1.6 = r −3.2 , roughly consistent with the scaling we observe in g 2 . Note that the angular dependence of g 2 (r) is nearly isotropic and clearly does not show a quadrupolar dependence. This is consistent with the angular dependence of˜ , not xy (see Fig. 2). We therefore conclude that rearrangement-induced shear strain in any direction can trigger rearrangements equally well. This result contradicts the assumption of many elastoplastic models that xy is solely responsible for triggering rearrangements.

B. Interplay of softness with rearrangements and strain
In training the machine-learning algorithm to obtain softness, we find that 90% of rearrangers have S > 0, while 84% of non-rearrangers have S < 0. Moreover, Fig. 4 shows that the softness distribution for rearrangers is very different from that of the whole population, and that the probability that a particle rearranges increases by four orders of magnitude as softness increases. These results verify that softness affects the propensity to rearrange very strongly.
In turn, rearrangements can affect softness. The average difference in softness of a rearranger immediately before and after the rearrangement is ∆S R = −0.75; the softness of a rearranger drops significantly when it rearranges. Rearrangements can also affect the softness of other particles; we plot the mean softness change ∆S(r) of a particle at r due to a rearrangement at the origin in Fig. 5. Rearrangements make directly contacting neighbors (r < 1) softer and non-contacting nearby particles (1 < r < 5) less soft. Rearrangements also make faraway particles (r > 5) softer or harder depending on the orientation. The distance and angular dependences of the far-field ∆S are consistent with the volumetric strain k (see Fig. 2), suggesting that it is caused by k. This is not surprising since softness is sensitive to density.
To understand the near-field effect of rearrangements on softness, we first note that in a thermal Lennard-Jones system, the mean softness of non-rearranging particles with a given initial softness S 0 evolves toward the global mean value for any S 0 [15] due to rearrangements of other particles. In other words, rearrangements of other particles tend to push softness towards its mean value. Here we study if the same effect exists in our quasistatically sheared system. For particles within a short distance r ≤ 1.6 to a rearranger, we plot the softness change vs. the original softness and perform a linear fit, presented in Fig. 6 (a). We plot the slopes c 1 (r) of such fits at several different r in Fig. 6 (b). For r < 10 and r > 30, c 1 is negative, indicating that softness in our system also has the tendency to approach its mean at these distances. However, c 1 is positive for 10 < r < 30, suggesting the opposite effect. The effect is small and negligible, and is probably because softness tend to increase in this range of r [see Fig. 5 (a)], and the softer a particle is, the floppier its local environment is, and the more tendency it has to deform, even if such deformation generally raises S. More important is the magnitude of c 1 (r): we see that the magnitude of c 1 (r) decays rapidly with r and is well described as a power law: |c 1 (r)| = 0.06r −3.2 . Finally, c 1 (r) appears to be independent of the angle θ.
Overall, our results suggest that the mean softness change of a particle with softness S at r when a particle at the origin rearranges is: where c 1 (r) is given in Fig. 6 (b), and b ≈ 207. To find c 0 , we subtract bk(r) from ∆S(r). Similar to c 1 , we do not find any angular dependence in c 0 . We plot its r-dependence in Fig. 6(b). Clearly, c 0 and c 1 exhibit similar power law decays; we find |c 0 (r)| = 0.3r −3.1 . With the fit, Eq. (3) yields the red curve in Fig. 5(a). Note that the red curve provides an excellent description of the black points (∆S(r)), capturing the sign as well as the magnitude except for two points at small r.

C. How strain and softness induce new rearrangements
We have shown that rearrangements give rise to deviatoric strain that in turn triggers new rearrangements. We have also shown that rearrangers tend to have high softness. Here we study how S and˜ work in tandem to induce rearrangements. When a particle starts rearranging at frame f , we rewind δf frames to calculate the shear strain exerted on this particle between f − δf and f , and the softness S at frame f − δf . As Fig. 7 shows, the amount of shear strain needed to trigger a rearrangement depends strongly on S, but only very weakly on δf . Thus, softer particles require less shear strain to start rearranging. This is consistent with earlier results in thermal systems that found that softer particles have lower activation energies to rearrange [15,17]. Indeed, we have conducted thermal molecular dynamics simulations to find energy barriers comparable to those predicted by Fig. 7 [19].

IV. DISCUSSION
In this paper, we study avalanches that occur during energy drops when a jammed binary Hertzian disk packing is sheared quasistatically, using steepest descent to follow the minimization process. We have confirmed the consistency of our interpretation of the roles of volumetric and deviatoric strain in three-dimensional systems [19]. We expect that the qualitative results of Fig. 8 apply quite generally to both two and threedimensional ductile disordered solids, which generically exhibit avalanche behavior. We find that (1) a rearrangement alters the softness of a nearby particle ac-cording to the difference between its softness and the mean softness. This behavior was first observed for 3D Lennard-Jones systems above the glass transition [15], indicating that it is quite general. (2) A rearrangement alters the softness of distant particles through volumetric strain. The existence of a transient volumetric strain, which has not been considered significant, is a feature of elasticity. The fact that local dilation/compaction increase/decrease the softness is consistent with the previously observed dependence of softness on local density in 3D Lennard-Jones mixtures and with our physical understanding of softness [15], and is therefore also quite general. (3) A rearrangement exerts a deviatoric strain on the rest of the system. This should be generally true for isotropic systems in any dimension. (4) The average yield strain decreases with increasing softness. This is consistent with previous results for 3D Lennard-Jones simulations [15], 2D colloidal glass experiments [26] and 3D aluminum polycrystal simulations [17], showing that the energy barrier for rearrangements decreases with increasing softness.
Due to the generic nature of our findings, the summary of our results in Fig. 8 should hold generally for avalanches in ductile disordered solids. Fig. 8 can be viewed as a structuro-elasto-plastic model that builds upon earlier elasto-plastic models. There are two main differences compared to the earlier models. First, we find that shear strain in any direction due to a rearrangement can trigger the next rearrangement equally well. Elastoplastic models typically focus on the component of the local shear strain with the same orientation as the global shear strain [8,27]. Second, and more significantly, we have elucidated how the local structural environment of a particle affects and is affected by rearrangements and strain.
Our results point to a few factors that may contribute to the ductile behavior observed. Future rearrangements are triggered by the total deviatoric strain rather than the xy-shear strain. As a result, rearrangements trigger successive rearrangements that are isotropically distributed, not concentrated in the direction of maximum xy-strain. In addition, a rearranger lowers the softness of nearby particles, discouraging them from rearranging, while on average raising the softness of distant particles, facilitating their rearrangement. Third, rearrangements tend to push the softness of nearby particles towards the mean, which is quite high for the ductile system. Our approach can be applied directly to systems that exhibit shear-banding and brittle failure to see whether the interplay is different in such systems. Earlier papers have shown that softness is readily identified in experimental systems for which the positions of particles can be tracked with time [1,28,29]. Our analysis for disentangling the interplay of softness, rearrangements and strain can therefore be applied directly to experiments as well as simulations. It is likely that the key to understanding ductile vs. brittle behavior is encapsulated in this interplay.

ACKNOWLEDGMENTS
We thank Hongyi Xiao, Robert J. S. Ivancic, Douglas J. Durian and Robert A. Riggleman for helpful discussions and their careful reading of the manuscript. This work was supported by the UPenn MRSEC NSF-DMR-1720530 (GZ), the Simons Foundation for the collaboration "Cracking the Glass Problem" award #454945 to AJL (GZ,SR,AJL), and Investigator award #327939 (AJL), the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-FG02-05ER46199 (GZ), and NSERC via a PGS-D fellowship (SR).

Appendix A: Continuum-elastic predictions for strain field induced by a rearrangement
The far field of rearrangement events has long been modelled as that of an Eshelby inclusion, which is the elastic response to a point strain source [22].
Elastoplastic models typically only consider σ xy , use an elastic kernel which assumes the medium to be incompressible and take the limit of infinite time (mechanical equilibrium). Since we are interested in understanding the course of avalanches during steepest descent, we need the kernel at finite times with overdamped dynamics. We sketch below the derivation of all components of the continuum strain field.
We begin by considering an infinite elastic medium subject to a point force turning on at t = 0 at the origin.
We wish to find G ik (r, t) such that Taking a Fourier transform in space and a Laplace transform in time gives us with the last equality holding for an isotropic medium in 2d. Heret is the vector normal toq.
We invert the spatial Fourier transform, and then the Laplace transform. The result is where Γ(0, x) ≡ ∞ x ds s −1 e −s is the incomplete gamma function (in this case, also the exponential integral function).
Differentiating this twice and symmetrizing over one of the indices allows us to compute G ijkl , the strain response to a dipole of force.
We obtain [r krl δ ij +r irj δ kl +r irk δ jl +r jrk δ il ] Following previous work, a dipole of xy shear strain at the origin is equivalent to a pair of force dipoles [22]. Assuming this source gives us the elastic strain field (now written in terms of the Poisson ratio ν and the "diffusion constants" D T ≡ µ η and D L ≡ λ+2µ The familiar power law dependences from elastic equilibrium are realized in the large-time limit These results together explain why the volumetric strain is observed to have a sin(2θ) dependence, and why the deviatoric strain magnitude is isotropic.
Notice, however, that dθ k(r, θ, t) = 0 for such a shear strain source. To explain the apparent nonzero value of dθ k(r, θ, t) for short times in our simulations, we must consider the effect of a transient expansion source. The local region surrounding a rearrangement might be expected, on average, to have a different volume than in the initial state.
In an infinite system, the kernel above gives for a point plastic compression at the origin: As long as the Poisson ratio is close to 1, this precisely conserves volume in an infinite system, when added to the point compression at the origin.
We expect that since our system is finite (and the short-time Poisson ratio is far from 1), this kernel would need to be modified near the boundaries of the system to satisfy the periodic boundary conditions and conserve the total volume. We find that it works adequately for the bulk for our data however, and our data at r close to the box size are difficult to resolve -we have chosen the y-range in the top-left box of Fig. 2 to exclude points beyond r = 30 because the error bars are comparable to the absolute value.
The full response to a given event will be a sum of the responses to strain [Eq. (A8)] and compression sources (A9) with appropriate prefactors, although for measurements where its contribution is nonzero we expect the strain source to be dominant.
Appendix B: Compare analytical and numerical k(r), , and xy results Since we have derived a analytical formulae for the strain, Eqs. (A8) and (A9), we can make comparison with our numerical results. We have numerically measured instantaneous elastic constants λ + 2µ = 0.3533 and ν = 0.3408 for our system by applying a small (10 −6 ) strain on the simulation box and measuring the force.
The time interval between frames, t, is not fixed since we record frames that are equidistant in configuration space; see the SI [19]. We plot the distribution of times between frames in supplementary Fig. S1, and find that the most probable time interval is t ≈ 100. The definition of our time implies that η = 1. With these parameters, Eq. (A9) predicts a Gaussian that decays to 0.1% of its peak height at r = 31, roughly consistent with the actual result presented in Fig. 2.
For the total deviatoric strain˜ and xy-strain xy , we have numerically confirmed that they decay as power laws:˜ =c/r 2 and xy = c xy /r 2 (Fig. 2), which matches the prediction in Eq. (A8). The prefactors, i.e., constantsc and c xy , were not predicted in Appendix A since our theory does not take into consideration the average amount of plastic strain caused by a rearranger.
Nevertheless, we can approximately measure this quantity. The strains in equation (A8) are for a plastic strain pl xy = 0 Aδ(r), i.e. the prefactor of the far-field strain is equal to the product of the area A = πr 2 D of the rearrangement and its plastic strain.
Why do our numerical results match the analytical derivations for shear strains produced by a shear source, Eqs. (A6) and (A7), in the infinite-time limit of Eq. (A8), but match that for the volumetric strain produced by a compression source, Eq. (A9), at a finite time? It turns out that at the rearranging site, the plastic shear occurs over a much longer time interval than the plastic compression. We plot these strain components at the rearranging site versus time in supplementary Fig. S5. If we approximate such strain-time curves with Gaussians, then the numerically measured strain at distance r should be the convolution of previously-derived finite-time analytical result and Gaussians, i.e., k(r, numerical) = c k 0 −∞ exp −αt 2 k(r, t − t )dt , xy (r, numerical) = c xy 0 −∞ exp −βt 2 xy (r, t − t )dt , (B4) where k(r, t) and xy (r, t) are given in Eqs. (A9) and (A6), respectively. We numerically compute these integrals for various parameters. For k, the integral fits numerical data well at α = 6.1197 × 10 −5 , as shown in Fig. 2. This indicates that the width of the Gaussian is about α −1/2 = 127.83, roughly consistent with supplementary Fig. S5. For xy , however, it turns out that Eq. (B4) cannot closely fit our numerical result, which decays slightly slower than r −2 (Fig. 2). No matter how small β is, Eq. (B4) gives an xy that decays slightly faster than r −2 . We see two possible reasons for this difference: (1) A finite size effect as r becomes comparable to the box size, or (2) the interference between simultaneous rearrangements in our numerical results. As we discuss in the last paragraph of supplementary Sec. III, we filter out