Atomic nonaffinity as a predictor of plasticity in amorphous solids
AAtomic nonaffinity as a predictor of plasticity in amorphous solids
Bin Xu,
1, 2
Michael L. Falk, ∗ Sylvain Patinet, and Pengfei Guan † Beijing Computational Science Research Center, Beijing 100193, China Department of Materials Science and Engineering,Johns Hopkins University, Baltimore, Maryland 21218, USA Department of Materials Science and Engineering,Mechanical Engineering, and Physics and Astronomy,Johns Hopkins University, Baltimore, Maryland 21218, USA PMMH, CNRS UMR 7636, ESPCI Paris, PSL University,Sorbonne Universit´e, Universit´e de Paris, F-75005 Paris, France (Dated: January 24, 2020)Structural heterogeneity of amorphous solids present difficult challenges that stymie the predic-tion of plastic events, which are intimately connected to their mechanical behavior. Based on aperturbation analysis of the potential energy landscape, we derive the atomic nonaffinity as an indi-cator with intrinsic orientation, which quantifies the contribution of an individual atom to the totalnonaffine modulus of the system. We find that the atomic nonaffinity can efficiently characterizethe locations of the shear transformation zones, with a predicative capacity comparable to the bestindicators. More importantly, the atomic nonaffinity, combining the sign of third order derivative ofenergy with respect to coordinates, reveals an intrinsic softest shear orientation. By analyzing theangle between this orientation and the shear loading direction, it is possible to predict the protocol-dependent response of plastic events. Employing the new method, the distribution of orientations ofshear transformation zones in a model two-dimensional amorphous solids can be measured. The re-sulting plastic events can be understood from a simple model of independent plastic events occurringat variously oriented shear transformation zones. These results shed light on the characterizationand prediction of the mechanical response of amorphous solids.
Understanding how the heterogeneity of amorphousstructures correlates with mechanical response remains asignificant challenge. Various indicators have been pro-posed to quantitatively predict the locations in the mate-rial susceptible to undergo plastic transformation. Someof these only consider the geometry of structure, such asfree volume [1, 2], five-fold symmetry [3, 4], and local de-viation from sterically favored structures [5], etc. Othersof these take the interaction between particles into con-sideration, like low-frequency normal modes [6–9], po-tential energy [10], local elastic modulus [11], flexibilityvolume [12], mean square vibrational amplitude [13], lo-cal thermal energy [14, 15], local yield stress [16, 17],and saddle points sampling [18], etc. Most of these in-dicators are inherently scalar quantities while deforma-tion mechanism must have a shear-like character. Thisis clearly borne out by the fact that the orientational na-ture of shear transformation zones (STZs), the defectspurported to be associated with plastic rearrangement,can be measured through their high sensitivity to thedeformation protocols. As verified in simulations, underdifferent loading orientations, glasses may exhibit con-trasting mechanical responses, with often entirely differ-ent activated STZs [14, 16–19]. Until now, one wouldneed to test the mechanical response in different orienta-tions to know the protocol-dependent mechanical behav-ior of one STZ [15–17]. While an intrinsic softest orien-tation for a STZ should exist, it is not directly revealed.Information regarding saddle points on the energy land-scape can reveal the intrinsic softest orientation, but this analysis incurs a high computational cost for sampling allof the low-barrier events. Thus, it is imperative to de-velop a low-cost indicator to quantitatively measure theprotocol-dependent behavior of plastic events. In thisletter, based on a perturbation analysis of the energylandscape, we derive a new indicator with intrinsic ori-entation, termed the atomic nonaffinity. We show thatit can efficiently predict the locations of the plastic rear-rangements in a two-dimensional amorphous solid, witha predicative capacity comparable to the best indicators.The relevant orientational information of STZs is natu-rally reflected in this parameter. Moreover, the distribu-tion of orientations of activated STZs is calculated, andunderstood through a simple model of independent STZswith isotropically distributed soft orientations.A two-dimensional glassy system, comprised of 10 particles, was prepared via the same gradual quench andthe same smoothed Lennard-Jones potential describedin Ref. [16, 17]. We start from a special state, in whichthe system is deformed to be very close to the trigger-ing strain of a plastic event via a protocol of athermalquasistatic shear. The spatial distribution of the normalmode with the lowest eigenvalue (refereed as the lowestmode, LM, in the following text) is shown in Fig. 1(a).A plastic event will be triggered in the region (Fig. S1in the supplemental materials (SM)) where the LM is lo-calized if the system was further sheared in this direction(denoted as the reference direction θ L = 0). However, ifthe system was further sheared in other directions, suchas θ L = − π or π , no plastic event would be observed a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n (c) (f)(d)(e) (a) (b) 𝜽 𝐋 = - 𝝅𝟒 𝜽 𝐬 𝝅𝟐𝟑𝝅𝟐 𝝅𝟒 𝝅𝟒𝝅𝟒 𝝅𝟒𝝅 𝑮 𝐋𝐌 ( 𝜽 𝐋 ) 𝜽 𝐬 𝝅𝟐𝟑𝝅𝟐 𝝅𝟒 𝝅𝟒𝝅𝟒 𝝅𝟒𝝅 ∆ 𝜸 𝐜 ( 𝜽 𝐋 ) × − 𝜽 𝐬𝒄 𝝅𝟐𝟑𝝅𝟐 𝝅𝟒 𝝅𝟒𝝅𝟒 𝝅𝟒𝝅 𝑮 𝒄 ( 𝜽 𝐋 ) FIG. 1. (a) The spatial distribution of normal mode withlowest eigenvalue. (b) The D field after shear with strainof ∆ γ = 6 × − in the orientation of θ L = − π/
4. The insetshows the stress-strain curve. (c) The predicted triggeringstrain (line) and triggering strain from simulation (circles) asa function of shear angle θ L . (d) The magnitude of nonaffinemodulus contribution from the lowest mode at different θ L .Blue line represents the range of θ L , where the plastic eventcan be triggered, while the red line represents the range of θ L ,where the plastic event can not be triggered. (e) The spatialdistribution of ˆ G ( θ L ) · sign(∆ γ c ( θ L )) for different orientations.(f) The atomic shear nonaffinity in different orientations forthe atom that has the maximum magnitude of atomic shearnonaffinity in (e). (Fig. 1(b) and Fig. S1 in SM) in the region, even fora large magnitude of strain. Obviously, this protocol-dependent mechanical behavior of amorphous systemscannot be clearly understood solely with scalar indica-tors. Here, we introduce the second and the third deriva-tive of the energy with respect to the vibrational coordi-nate ( q ∗ ) of the LM (denoted as λ ∗ and η ∗ respectively)and the first derivatives of stress of the system with re-spect to q ∗ (denoted as ∂σ xx ∂q ∗ , ∂σ yy ∂q ∗ and ∂τ xy ∂q ∗ respectively).The triggering strain for different shear orientation (de-noted as θ L ) can be derived as (see SM for details ofderivation) ∆ γ c ( θ L ) = λ ∗ η ∗ V ∂τ ( θ L ) ∂q ∗ , (1)where V is the volume of system and ∂τ ( θ L ) ∂q ∗ is the firstderivative of the shear stress with respect to q ∗ at θ L , which is equal to − ( ∂σ xx ∂q ∗ − ∂σ yy ∂q ∗ ) sin 2 θ L + ∂τ xy ∂q ∗ cos 2 θ L .Moreover, a softest shear orientation of the LM, ˜ θ s , as-sociated with the smallest positive triggering strain canbe defined as tan 2˜ θ s = − (cid:18) ∂σxx∂q ∗ − ∂σyy∂q ∗ ∂τxy∂q ∗ (cid:19) , with (2) ∂τ ( θ s ) ∂q ∗ · η ∗ > , and θ s ∈ ( − π , π ] , to take the symmetry of shear into consideration. Thecomputed quantity ˜ θ s corresponds to the orientationalong which the first derivative of shear stress with re-spect to q ∗ is the largest. To verify the validity of thepredictions of Eq. (1) and Eq. (2), further simulationswere performed to directly measure the ∆ γ c ( θ L ) and ˜ θ s .As shown in Fig. 1(c), the predictions agree well withthe simulation results, which suggests that the analysisof LM is successful for fully understanding the protocol-dependent mechanical response of the system close to theinstability. When the system is far from the instability,we suppose that all the modes, especially the ones withsmall eigenvalues, should be taken into consideration.To develop an indicator that takes all modes intoconsideration while maintaining the orientational infor-mation of each mode, we check how different modescontribute to the modulus of system. Following Mal-oney [6, 20], the elastic constants of amorphous solids canbe derived from the second derivative of the total poten-tial energy with respect to strain in athermal quasistaticdeformation. These can be rewritten in the coordinatesof eigenbasis as C ijkl = 1 V (cid:32) ∂ U∂(cid:15) ij ∂(cid:15) kl + (cid:88) m ∂ U∂q m ∂(cid:15) ij · dq m d(cid:15) kl (cid:33) , (3)where U is the potential energy, and q m is the m th coor-dinate of the eigenbasis of the Hessian matrix ( ∂ U∂r i ∂r j ).The first term (Born term) of Eq. 3, accounts for affinedisplacement and is insensitive to the structural stabilityinduced by different cooling rate [21] . The second term,containing the contribution from nonaffine relaxation ineach normal mode, termed the nonaffine modulus ( ˜ C )here, is sensitive to the structural stability. By expressingthe stress as as σ ij = V ∂U∂(cid:15) ij and the nonaffine “velocity”in quasistatic deformation as dq m d(cid:15) kl = − λ m ∂σ kl ∂q m [6], where λ m is the eigenvalue of m th normal mode, the ˜ C ijkl canbe rewritten as˜ C ijkl = (cid:88) m ˜ C ijkl,m = − Vλ m ∂σ ij ∂q m ∂σ kl ∂q m , (4)where ˜ C ijkl,m is the contribution from the m th normalmode, which is always negative. In shear protocols, theshear modulus is the most important elastic constant.Thus, we focus on the nonaffine shear modulus ( ˜ G ) andthe contribution from each mode ( ˜ G m ). The ˜ G m can becalculated by ˜ G m ( θ L ) = − Vλ m (cid:18) ∂τ ( θ L ) ∂q m (cid:19) , (5)which depends on the orientation θ L . The nonaffine shearmodulus contribution of the dominant LM, ˜ G LM , for thestate described in Fig. 1(a)-(c) is shown in Fig. 1(d).The blue line represents the orientational range, where∆ γ c >
0, i.e. ∂τ ( θ L ) ∂q ∗ · η ∗ >
0, and the event can betriggered. The red line represents the orientational range,where the event can not be triggered no matter how largeis the shear strain. Moreover, a softest shear orientationcan also be defined by the largest value of ˜ G LM in theblue range which is and should be consistent with the ˜ θ s derived from Eq. 2.So far our results have been discussed with respect toeigenbasis. We need to transform those results into basisof Cartesian coordinates to develop an indicator in realspace. By expressing the normalized eigenvector in realspace as Ψ m = (cid:80) n,α c mnα e nα , where e nα is a unit vectorcorresponding to the displacement of n th atom in the α (= x or y ) direction, and c mnα is the projection of the m th eigenvector onto the e nα , the ˜ C ijkl can be rewrittenas ˜ C ijkl = (cid:88) n ˆ C ijkl,n = (cid:88) n (cid:88) m,α − Vλ m ∂σ ij ∂q m ∂σ kl ∂q m c mnα . (6)Here ˆ C ijkl,n is the atomic nonaffinity of the n th atom. Asmost local plastic rearrangements are shear-like [22, 23],the atomic shear nonaffinity (ASN) is the most impor-tant component when investigating the STZs and can bewritten asˆ G n ( θ L ) = (cid:88) m,α − Vλ m (cid:18) ∂τ ( θ L ) ∂q m (cid:19) c mnα . (7)Obviously, the value of ˆ G n depends on the orientation θ L . As a result, the spatial distribution of ˆ G n in the pre-vious case shown in Fig. 1(e) exhibits a clear orientation-dependent behavior in the region where the plastic eventis located. More negative values of ˆ G n mean that thecorresponding atom is easier to trigger in the orientation θ L . According to the ˆ G n distribution with different ori-entations, θ L = 0 is the easiest shear direction for theplastic event comparing with θ L = π and π , which isconsistent with the direct loading test in Fig. 1(b) andFig. S1. Moreover, the atom located in the core regionhas the maximum magnitude of atomic shear nonaffinity,denoted as ˆ G c ( θ L ). Fig. 1(f) shows the θ L -dependent ˆ G c ,which has a similar shape as the ˜ G LM . It is attributedto the LM with smallest eigenvalue dominating the vari-ation of ˆ G c , which can be inferred from Eq. 7. Thus, wecan define the softest shear orientation for the n th atom as the softest shear orientation of the mode that domi-nates the variation of ˆ G n . The softest shear orientationof the n th atom is defined as θ n,s = ˜ θ i,s , i = argmax m (cid:88) α | ˜ G m (˜ θ m,s ) | c mnα , (8)and the calculated softest shear orientation of the coreatom ( θ cs ) is shown in Fig. 1(f). The consistency of theproposed softest shear orientations for one STZ from thethree parameters, i.e. the directly calculated triggeringstrain (Fig. 1(c)), nonaffine modulus of the lowest mode(Fig. 1(d)) and the atomic shear nonaffinity (Fig. 1(f)),implies that the θ s defined from atomic shear nonaffinityis effective to for characterizing the orientations of STZs.It is worth noting that the softest shear orientation ob-tained from atomic shear nonaffinity is parameter-freeand can be obtained without triggering the plastic events.We believe it to be consistent to the definition proposedby Nicolas et al. [24], where orientations are obtained byfitting the local strain field induced by the triggering ofplastic events to the Eshelby strain field.Now we have seen the predicative capacity of atomicshear nonaffinity in protocol-dependent mechanical re-sponse of a plastic event close to instability. But whatif the system is not close to instability? And how pre-dictive is this indicator? Predicting plastic events inan amorphous system by analyzing the local indicatorsof initial structure has been extensively studied in theliterature [7–9, 14, 16, 18, 19, 25, 26]. To comparethe reliability of local indicators for predicting plasticevents, one hundred two-dimensional samples preparedwith the same thermal history as the previous samplewere employed for local properties calculations. Theathermal quasistatic shear deformation with a strain stepof ∆ γ xy = 10 − was then applied to each sample, andeach stress drop in the stress-strain curve was associatedwith the resulting atomic rearrangements correspondingto one plastic event. Nonaffine rate [6] was calculatedfor the configurations just before the stress drops, andthe atom with maximum nonaffine rate was identified asthe core atom, whose index is denoted as ID N for N th plastic event. To compare the success of different indica-tors, we transform those indicators to a rank correlation(RC) value following the analysis performed by Patinetet al. [16] as RC Ψ ( n ) = 1 − Ψ ( n ) , (9)where Ψ is one of the indicators, and CDF Ψ ( n ) is the cu-mulative distribution function value of the n th particle forΨ. The spatial distribution field of the calculated RC ˆ G with θ L = 0 is shown in Fig. 2(a). The first ten plasticevents in shear protocols with θ L = 0 are almost all lo-cated at high RC ˆ G regions, which implies the highly pre-dictive power of ˆ G . To quantitatively compare the predic-tive power regarding plastic events for different local in-dicators, the relationship between the locations of plasticevents and the corresponding values of local indicators isdescribed by the average of RC Ψ (ID N ) over 100 samples.The average RC Ψ (ID N ) of investigated local indicators, such as participation fraction (PF) [7–9] in the lowest1% of normal modes, the magnitude of “nonaffine veloc-ity” [6] (nonaffine rate, NR), the mean square vibrationalamplitude (MSVA) [13], the local yield stress (LYS) [16]and the ASN, are shown in Fig. 2(b). The LYS presentsthe highest predictive power in the early stage, since non-linear response to shear is considered. The MSVA, andASN shows comparable predictive power, and the otherindicators have lower predictive power than those three.Worth to note that the predicative power of the indica-tors depends on the stability of configurations and weshow predictive power of those indicators for configura-tions prepared via instantly quenching from high temper-ature liquids in SM, where one can see MSVA and ASNoutperform LYS. We also note that orientational infor-mation in ASN is incomplete and it, as a modulus, hasthe same value for θ L = 0 and π , while local regions gen-erally have different mechanical behavior for those twoprotocols. Thus, taking into account the intrinsic soft-est shear orientation ( θ s ) should effectively predict theprotocol-dependent behavior of plastic events.As discussed in Fig. 1, we expect that the plastic eventsinduced when shearing along direction θ L should be lo-cated at the atoms with | θ s − θ L | < π , and here we testthis expectation in one of the previous samples. We fo-cused on the ”soft” atoms in the sample with RC ˆ G > θ s . The corre-lation between the atom with RC ˆ G > | θ s | < π ( | θ s | > π ) and the first ten plastic events of θ L = 0( θ L = π ) direction is shown in Fig. 2(c) (Fig. 2(d)). Thecorrelation in both Fig. 2(c) and (d) indicate that the pre-dictive power can be increased by screening the high cor-relation regions with intrinsic softest orientation of ˆ G n .(Similar results about protocols of θ L = π and − π areshown in SM.)However, there still exist some number of events thatare not caught by the criterion | θ s − θ L | < π . This canbe attributed to the rotation of θ s during deformation,since the θ s is calculated mainly based on a second-orderperturbation method, and higher order terms and nonlin-ear interactions between different modes can lead to therotation of θ s . To obtain the statistics of the rotation of θ s , the softest shear orientations of core particles of allthe plastic events before shear strain 0 .
12 with θ L = 0in 10 samples are calculated based on the configurationsjust before each event or just after the last event (illus-trated in Fig.3(a)). The distribution of the calculated θ s of those core particles in the configurations just be-fore triggering are shown in Fig. 3(b) and all θ s satisfythe criterion | θ s − θ L | < π , which is what we expectedfor systems close to instability, as discussed previously.Moreover, the peak of probability density is located at θ L implying that the region with the intrinsically soft- | 𝜽 𝒔 | ≤ 𝝅𝟒 𝜽 𝒔 > 𝝅𝟒𝜽 𝐋 = 𝜽 𝐋 = 𝝅𝟐 FIG. 2. (a) The correlation field of atomic shear nonaffinityˆ G ( θ L = 0) and the locations of the first ten plastic events(black circles) triggered in shear protocols with θ L = 0. (b)Correlation between the indicators including local yield stress(LYS), mean square vibrational amplitude (MSVA), atomicshear nonaffinity (ASN), participation fraction (PF), non-affine rate (NR) with the locations of plastic events as a func-tion of the index of the events. Averages are taken over win-dows of five events. Error bar at each window represents thestandard deviation of the mean. (c) Orange circles representthe atoms with RC ˆ G > | θ s | < π . Black circles mark thelocations of the first ten plastic events with θ L = 0. (d) Greencircles represent the atoms with RC ˆ G > | θ s | > π . Tri-angles mark the locations of the first ten plastic events with θ L = π . est orientation closest to the imposed shear orientationis easiest to trigger. However, the distribution is broad-ened as shown in Fig. 3(c) for the calculated θ cs based onthe configurations just after the triggering of the previ-ous event. In this analysis only approximately 75% of theplastic events satisfy the criterion. More statistics abouthow orientations calculated by our perturbation methodchange are presented in SM.The distribution of orientations of the triggered plasticevent shown in Fig. 3(b) is regular. It can be understoodby a simple model of independent plastic events with in-trinsic orientations. In this model, we assume that amor-phous solids are isotropic and that the shear-orientation-dependent triggering strain can be derived from Eq. 1as γ c ( θ L ) = γ c ( θ s )cos[2( θ s − θ L )] , (10)where θ s is the softest shear orientation of a STZ. If weassume that the number density for a particular softestshear direction θ s at different triggering strains γ c ( θ s )(noted as γ c,s ) follows a power law ρ ( γ c,s ) = Aγ αc,s [27,28], the probability density distribution of orientations (cid:1)(cid:2) (cid:1) (cid:3) , (cid:4) (cid:1) (cid:3) , (cid:4) (cid:1) (cid:3) , (cid:4) (cid:5) (cid:1) (cid:3) , (cid:4) (cid:6) (cid:1) (cid:3) , (cid:4) , (cid:5) (cid:1) (cid:3) , (cid:4) , (cid:6) B efore triggering A fter triggering ID (cid:4) ID N : Index of the core atom for N th event (a)(b) (c) FIG. 3. (a) A schematic diagram introducing the notationsused in this figure. (b) The distribution of the softest shearorientations for all the plastic events that are triggered beforeshear strain 0 .
12 in 10 samples. The softest shear orienta-tions are calculated based on the configuration just beforeeach event. Red line follows the function π cos (2∆ θ ). (c)Similar to (b), but the softest shear orientations for each eventare calculated based on the configurations that are just aftertriggering of the previous event. ˆ ρ ( θ s − θ L ) (denoted as ˆ ρ (∆ θ )) will follow (see SM fordetails of derivation)ˆ ρ (∆ θ ) = k cos α +1 (2∆ θ ) (11)The probability density distribution in Fig. 3(b) corre-sponds to α = 1, as it closely fits a distribution function π cos (2∆ θ ) (red line in Fig. 3(b)). These results arealso supported by the probability distribution functionof local yield stress of the samples, in which α ≈ .
1, asshown in Ref.[16].In conclusion, we have derived a general andparameter-free indicator, the atomic nonaffinity. It iswell-defined and is easy to apply in systems beyond the2d Lennard-Jones system discussed here. The atomicnonaffinity has a clear physical meaning in that the sum-mation of atomic nonaffinities corresponds to the totalnonaffine modulus of the system. The softest shear orien-tation of each region is defined based on the atomic shearnonaffinity and stems from anisotropy of the shear stressderivative against the coordinate of the low-frequencymode in different orientations. Combining the sign ofthird order derivative of energy with respect to coordi-nates, it reveals the intrinsic orientation of the plasticrearrangement and directly connects to the anisotropicmechanical response of local regions, which is importantfor understanding aspects of the mechanical behavior ofamorphous solids not directly reflected or defined in other indicators. As atomic nonaffinity is developed based onthe nonaffine response of atoms upon deformation, it nat-urally has a good correlation with the plastic events, out-performing most other soft-mode-based approaches thatneglect the orientational information of local structureswith notable exception of the atomic mean square vi-brational amplitude [13]. Mechanical behavior must becorrelated with structure, and we anticipate that thismethod will be important for elucidating the structuralorigin of the anisotropic mechanical response in specificsystems.B.X. and P.F.G acknowledge financial support bythe National Natural Science Foundation of China(NSFC, Grants No. 51571011/U1930402), the MOST973 program (No. 2015CB856800). M.L.F. ac-knowledges support provided by NSF Grant AwardsNo. 1408685/1409560 and No. 1910066/1909733. We ac-knowledge the computational support from the BeijingComputational Science Research Center. ∗ [email protected] † [email protected][1] F. Spaepen, Acta Metallurgica , 407 (1977).[2] F. Zhu, A. Hirata, P. Liu, S. Song, Y. Tian, J. Han,T. Fujita, and M. Chen, Physical Review Letters ,215501 (2017).[3] M. Wakeda, Y. Shibutani, S. Ogata, and J. Park, Inter-metallics , 139 (2007).[4] Y. Hu, F. Li, M. Li, H. Bai, and W. Wang, Naturecommunications , 8310 (2015).[5] H. Tong and H. Tanaka, Physical Review X , 011041(2018).[6] C. Maloney and A. Lemaˆıtre, Physical Review Letters , 195501 (2004).[7] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R.Reichman, Nature Physics , 711 (2008).[8] M. L. Manning and A. J. Liu, Physical Review Letters , 108302 (2011).[9] J. Ding, S. Patinet, M. L. Falk, Y. Cheng, and E. Ma,Proceedings of the National academy of Sciences of theUnited States of America , 14052 (2014).[10] Y. Shi, M. B. Katz, H. Li, and M. L. Falk, PhysicalReview Letters , 185505 (2007).[11] M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Bar-rat, Physical Review E , 026112 (2009).[12] J. Ding, Y.-Q. Cheng, H. Sheng, M. Asta, R. O. Ritchie,and E. Ma, Nature Communications , 13733 (2016).[13] H. Tong and N. Xu, Physical Review E , 010401(R)(2014).[14] J. Zylberg, E. Lerner, Y. Bar-Sinai, and E. Bouchbinder,Proceedings of the National Academy of Sciences ,7289 (2017).[15] Z. Schwartzman-Nowik, E. Lerner, and E. Bouchbinder,Physical Review E , 060601(R) (2019).[16] S. Patinet, D. Vandembroucq, and M. L. Falk, PhysicalReview Letters , 045501 (2016).[17] A. Barbot, M. Lerbinger, A. Hernandez-Garcia,R. Garc´ıa-Garc´ıa, M. L. Falk, D. Vandembroucq, and S. Patinet, Physical Review E , 033001 (2018).[18] B. Xu, M. L. Falk, J. F. Li, and L. T. Kong, PhysicalReview Letters , 125503 (2018).[19] O. Gendelman, P. K. Jaiswal, I. Procaccia, B. Sen Gupta,and J. Zylberg, Europhys. Lett. , 16002 (2015).[20] T. Barron and M. Klein, Proceedings of the Physical So-ciety , 523 (1965).[21] Y. Q. Cheng and E. Ma, Physical Review B , 064104(2009).[22] A. Argon, Acta Metallurgica , 47 (1979).[23] M. L. Falk and J. S. Langer, Physical Review E , 7192(1998). [24] A. Nicolas and J. Rottler, Physical Review E , 063002(2018).[25] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Mal-one, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu,Physical Review Letters , 108001 (2015).[26] E. Lerner, Physical Review E , 053004 (2016).[27] S. Karmakar, E. Lerner, and I. Procaccia, Physical Re-view E , 055103(R) (2010).[28] H. G. E. Hentschel, P. K. Jaiswal, I. Procaccia, andS. Sastry, Physical Review E92