Rigidity-Controlled Crossover: From Spinodal to Critical Failure
RRigidity-controlled crossover: from spinodal to critical failure
Hudson Borja da Rocha
1, 2, ∗ and Lev Truskinovsky † LMS, CNRS-UMR 7649, Ecole Polytechnique, Université Paris-Saclay, 91128 Palaiseau, France PMMH, CNRS-UMR 7636 PSL-ESPCI, 10 Rue Vauquelin, 75005 Paris, France (Dated: January 14, 2020)Failure in disordered solids is accompanied by intermittent fluctuations extending over a broadrange of scales. The implied scaling has been previously associated with either spinodal or criti-cal points. We use an analytically transparent mean-field model to show that both analogies arerelevant near the brittle-to-ductile transition. Our study indicates that in addition to the strengthof quenched disorder, an appropriately chosen global measure of rigidity (connectivity) can be alsoused to tune the system to criticality. By interpreting rigidity as a timelike variable we reveal anintriguing parallel between earthquake-type critical failure and Burgers turbulence.
Failure in disordered solids takes place when elasticity(reversibility) breaks down [1]. The implied abrupt me-chanical degradation can be associated with brittle rup-ture [2], large plastic avalanche [3], or result from othernucleation type event [4]. In strain controlled experi-ments, failure may be accompanied by a dramatic stressdrop, and the challenge is to predict and control suchundesirable events.The mechanism of failure in random elastic systemsis nontrivial because of the intricate interplay betweenthreshold-type nonlinearity, quenched disorder and long-range interactions. While the strength of disorder, thesystem size, and the range of elastic interactions areknown to affect the failure mechanism [5–8], here we fo-cus on the role of system’s rigidity , which has recentlyemerged as another relevant factor in failure-related phe-nomena [9–11].Failure in disordered solids is characterized by scale-free statistics of large events. The associated intermit-tency has been linked to the existence of either spin-odal [3, 12–15] or critical points [7, 16–18]. At largedisorder and infinite system size, failure is known to belinked to percolation [19–23]; however, the physical na-ture of failure at finite disorder remains a subject of de-bate [3, 7, 18, 24].In this Letter, we use an analytically tractable mean-field model to show that both spinodal and critical scalingbehaviors can coexist near the threshold of the brittle-to-ductile transition [25–31]. Ductile response is under-stood here in the sense of stable development of smallavalanches representing micro-failure events [32, 33].Brittle response necessarily involves large events repre-senting system size instabilities [34, 35].Our starting point is the fiber bundle model (FBM)with global stress redistribution [36]. This model wasused to explain a variety of physical phenomena fromfailure of textiles [37], and acoustic emission in loadedcomposites [38] to earthquake dynamics [39]. It is usu-ally studied in the stress control setting, where failure isbrittle and scaling is spinodal [2, 40]. To address failureunder strain control and to be able to tune the systemto criticality, we drive the system differently, using an external harmonic spring [4, 41].In our analysis, brittle failure emerges as a supercrit-ical, while ductile failure as a subcritical phenomenon.The critical behavior can be associated with the brittle-to-ductile transition and we show that due to superuni-versality of the mean-field models [42], the equilibriumand out-of-equilibrium exponents are the same.The main focus of this Letter, however, is the role ofthe system’s rigidity [43] as the regulator of the brittle-to-ductile transition. It is known that rigid, crystal-likesolids subjected to stresses fail catastrophically [44]. In-stead loose, marginally jammed solids fail gradually [9–11]. In view of the minimal nature of our model, wecould construct analytically the rigidity-disorder phasediagram delineating the domain of ductile behavior atlow rigidity and high disorder from the domain of brittlebehavior at low disorder and high rigidity.One of our crucial findings is that in the brittle-to-ductile crossover region, which bridges robust spinodalcriticality with tuned classical criticality, the transitionalexponents are non-universal, depending sensitively onsystem size, disorder, and rigidity. We also show thatwhen rigidity can be conditioned by the system size, fail-ure becomes brittle in the thermodynamic limit, and scal-ing survives only as a finite size effect, cf. [7, 9].Equilibrium (static) avalanches, corresponding tojumps between different globally minimizing configura-tions, have been previously linked to Burgers shocks[45, 46]. Here we extend this analogy showing that ifrigidity is interpreted as "time", and strain as "space",the brittle-to-ductile transition and the associated crit-ical behavior can be viewed as a "finite time" Burgersturbulence [47]. Given that our model is essentially amean-field version of the Burridge-Knopoff model [39],the developed analogy reinforces a conceptual link be-tween earthquakes (fracture) and turbulence [48].Consider a discrete system with dimensionless energy: H = 1 N N (cid:88) i =1 (cid:20) u i ( x i ) + λ X − x i ) (cid:21) + Λ2 ( ε − X ) , (1) a r X i v : . [ c ond - m a t . d i s - nn ] J a n where u i ( x ) is a Lennard-Jones type potential of a break-able element, X is a Weiss-type mean field accountingfor the interaction among breakable elements, and ε isthe controlling parameter representing the harmonic in-teraction of the field X with the loading device, seeFig. 1(a). For determinacy, we assume that the poten-tial u i ( x ) is piece-wise quadratic: u i ( x ) = ( x / l i − x ) + ( l i / x − l i ) , where Θ is the Heaviside function;for x ≤ l i , the element is intact, while for x > l i , itis broken. Here, l i are random numbers drawn fromthe probability distribution f ( l ) . In our numerical il-lustrations, we use Weibull’s distribution with density f ( x ) = ρx ρ − exp ( − x ρ ) ; broad disorder corresponds tosmall ρ [49]. In our generalization of the FBM (1), weintroduced two new parameters: the internal stiffness λ ,and the external stiffness Λ = κ/N , where κ ∼ N is theeffective elasticity of the elastic environment [41]. u ( x ) u N ( x N ) λλN Λ εX (a) (b) (c) . . . εσ Brittle 1 2 3 εσ Ductile
Figure 1. (a) Schematic representation of the system; (b)brittle reponse at
Λ = 0 . ; (c) ductile response at Λ = 5 .Solid black lines: equilibrium path, thin black line: out-of-equilibrium paths; grey lines: metastable states. Parameters: N = 100 , λ = 1 , ρ = 4 . In Fig. 1(b,c), we illustrate the typical behavior of thelocal and global minima of (1) by showing the relationbetween the applied strain ε and the conjugate stress σ = Λ( ε − X ) , see also [50]. Our Fig. 1(b) shows thebrittle behavior, which includes a system size transitionfrom the partially broken to the fully broken state. Incontrast, our Fig. 1(c) illustrates the ductile behavior,characterized by the gradual accumulation of damage.The equilibrium (global minimum) deformation paths areshown in Fig. 1(b,c) by thick black lines. We assumethat failure is reversible, and show by thin black linesthe out-of-equilibrium (marginally stable) paths that aredifferent for loading and unloading.The boundary separating brittle and ductile regimesdepends on the strength of the disorder (our parameter ρ ) and on the dimensionless parameter ν = λ Λ(1 + λ ) , (2)which we interpret as a measure of the structural rigidity of the system [43, 51–53]. When ν is small, meaningthat either Λ is large or λ is small, individual breakableelements interact weakly and the limit λ → can beassociated with the (jamming) threshold beyond which Figure 2. (a) Ensemble averaged brittle-to-ductile transi-tion line; (b) typical averaged stress-strain curves; (c) non-equilibrium avalanche distribution; (d) equilibrium avalanchedistribution. Parameters: N = 10 , λ = 1 . The avalanchedistributions was averaged over realizations. the rigidity is lost [11]. Instead, when ν → ∞ the systemcan be viewed as overconstrained [9, 10].The ensemble averaged brittleness/ductility thresh-old can be found by solving the system of equations − f ( x c ) = f (cid:48) ( x c ) x c and − F ( x c ) = f ( x c ) x c − /ν ,where F ( x ) = (cid:82) x f ( x (cid:48) ) dx (cid:48) is the cumulative distributionof thresholds [50]. In particular, for Weibull-distributedthresholds, the line separating brittle from ductile behav-ior is given by the equation ν ∗ = exp (1 /ρ + 1) /ρ , see Fig.2(a).In the limit N → ∞ the avalanche distribution in themodel (1) can be computed analytically [36, 50] p ( s ) = s s − s ! (cid:90) ∞ (cid:0) − g ( x ) (cid:1) f ( x ) g ( x ) e − h ( x ) s dx, (3)where g ( x ) = f ( x ) x/ (1 − F ( x ) + ν − ) , and h ( x ) = g ( x ) − ln g ( x ) . In the large-event-size asymptotics theuniversal pre-integral multiplier s s − /s ! ∼ s − / rep-resents the classical mean-field contribution, reflectingthe built-in statistics of Brownian return times [39, 54].In the limit s → ∞ the integrated distribution can beobtained by the saddle-point approximation around theglobal minimum, x , of the function h ( x ) [55]. It is a rootof either g ( x ) = 1 or g (cid:48) ( x ) = 0 , and the emergence ofsuch two cases is a general feature of mean-field models[56].Consider first the out-of-equilibrium path (dynamicavalanches). Then, if g (cid:48) ( x ) = 0 while g ( x ) (cid:54) = 1 , weobtain p ( s ) ∼ s − e − s ( h ( x ) − . This is a sub-critical dis-tribution describing the ductile (POP) regime [57], dom-inated by uncorrelated random events. If g ( x ) = 1 but g (cid:48) ( x ) (cid:54) = 0 , the point x is spinodal, and the distribu-tion is super-critical, characterizing the brittle (SNAP)regime [57]. Neglecting the system-size events, we canwrite the corresponding local distribution in the form p ( s, x ) ∼ s − / ( x − x ) e − s g (cid:48) ( x ( x − x ) . The avalanchesize diverges near x , and the integrated distributiontakes the classical form p ( s ) ∼ s − / [36]. Finally, if g ( x ) = 1 and g (cid:48) ( x ) = 0 , the local distribution reads p ( s, x ) ∼ s − / ( x − x ) e − s g (cid:48)(cid:48) ( x ( x − x ) . The character-istic avalanche size again diverges near x and the inte-grated distribution takes the form p ( s ) ∼ s − / . This isthe critical (crackling) regime [57] associated with brittle-to-ductile transition; the exponent / has appeared pre-viously in the context of composite FBM involving break-able and unbreakable springs [58]. Other values of the ex-ponents also appeared in the more complex FBM basedmodels describing richer physics [59].The computed critical exponents coincide with theones known for the mean-field RFIM [14, 60], becausethe energy (1) can be mapped on the soft-spin RFIM. Tothis end we need to minimize out the variable X , whichgives H = − (1 /N ) (cid:80) i,j Jx i x j − (1 /N ) (cid:80) i [ Hx i − v i ( x i )] , where v i ( x ) = u i ( x ) + x + λ Λ ε/ λ + Λ) , see also [50].Note that the Lennard-Jones type potential u i ( x ) wastransformed along the way into the double-well poten-tial v i ( x ) . Other mean-field formulations leading to thesame spinodal and critical exponents that are relevant foramorphous plasticity are discussed in [18, 24]; the sametwo main regimes have been also identified for some sand-pile automata [56]. Interestingly, a numerical analysis ofa non-mean-field model of a structural phase transitionreveals the possibility of a similar coexistence of two scal-ing behaviors [4]. . . ρτ N = 10 N = 10 Figure 3. Finite size crossover associated with brittle to duc-tile transition. Each curve gives the value of the scaling ex-ponent averaged over 250 realizations. Parameters λ = 1 and ν = 0 . (critical value at ρ = 4 ). The exponents andthe uncertainty were computed using the maximum likelihoodmethod [61]. In finite size systems, the crossover from the robust spinodal scaling in the brittle regime (exponent / ) tothe non-robust critical scaling (exponent / ) takes placein an extended transition zone, where the system exhibitsnon-universal exponents, see Fig. 3. The ubiquity of suchtransitional phenomena may explain the large scatter inreported scaling behavior of disordered solids [62–64].The mean-field model (1) can be used to demonstratedirectly the super-universality of the critical regime[42, 65–68]. For instance, one can show that the exponent / is valid for both out-of-equilibrium and equilibriumpaths [50]. Instead, the spinodal criticality, which ex-ists in the out-of-equilibrium model, disappears in theequilibrium model because the SNAP event takes placebefore the spinodal point is reached. Integrating theavalanches should be then performed only up to someMaxwellian x ∗ < x , and since in this case the func-tion h ( x ) attains its minimum at the boundary, we ob-tain p ( s ) ∝ s − / e − s (1 − h ( x ∗ )) . While this distributionhas the same exponent / as in the case of the out-of-equilibrium path, the scaling is now obscured by theexponential cut off. (a) (b) . . . νεσ ε ∗ ν ∗ εν Figure 4. (a) Time (rigidity) evolution of the randomly dis-tributed initial Burgers data σ ( ε ) (red) at ρ = 2 , and λ = 1 ;black line corresponds to ν ∗ = exp (1 /ρ + 1) /ρ . (b) Shockmerging with critical complexity appearing at ν = ν ∗ . Thickblack line shows the shock in the averaged system whichemerges in the limit N → ∞ . We now turn to an intriguing analogy betweenthe equilibrium version of the model (1) and Burg-ers turbulence [47]. If we minimize out the vari-ables x i in (1) and consider the thermodynamic limit N → ∞ [50], the equilibrium problem reducesto finding ˜ H ( ε, ν ) ∼ min X ∈ R (cid:8) ν ( ε − X ) + q ∞ ( X ) (cid:9) , where q ∞ ( z ) = [1 − F ( (cid:112) λ/ ( λ + 1) z )]( z /
2) + (cid:112) λ/ ( λ + 1) (cid:82) z f ( (cid:112) λ/ ( λ + 1) z (cid:48) )( z (cid:48) / dz (cid:48) . We can nowuse the Hopf-Lax formula [69] to turn this variationalproblem into a Cauchy problem for a Hamilton-Jacobiequation ∂ ν ˜ H + ( ∂ ε ˜ H ) = 0 , where the rigidity ν playsthe role of time. This equation must be supplemented bythe initial condition ˜ H ( ε,
0) = q ∞ ( ε ) . Then the tension σ = ∂ ε ˜ H satisfies the inviscid Burgers equation ∂ ν σ + σ∂ ε σ = 0 (4)with initial condition σ = ∂ ε q ∞ ( ε ) . Interestingly, theviscous Burgers equation for σ and the correspondingKPZ equation [70] for H emerge as a finite size effect inthe model (1) with finite temperature.As a result of the reduction of the problem (1) to(4), avalanches become shock waves [45]. In the aver-aged model, the ductile-to-brittle transition can be thenassociated with the shock formation at a finite rigidity ν ∗ = min ε ∈ R (cid:8) − /∂ ε σ ( ε ) (cid:9) , see Fig. 4(b); in the ( ε, ν ) plane this "event" becomes a direct analog of the liquid-vapor critical point.At finite N , the "evolution" equation for the stress re-mains the same as in the case N → ∞ , while the initialcondition changes to σ = ∂ ε q ( ε ) = N − (cid:80) Ni = i ε Θ( l i − (cid:112) [ λ/ ( λ + 1)] ε ) , see [50] for details. In Fig. 4(a) we showhow the increase of rigidity transforms the ductile re-sponse, where avalanches take the form of small Burgersshocks (POP events), into the brittle response with asingle Burgers shock representing a system size SNAPevent. In Fig. 4(b) we track the position of individualshocks and visualize their merging sequence.To highlight the critical nature of the system withrigidity value close to ν ∗ , we studied the ν dependenceof the number of shocks n . In Fig. 5, we show the stan-dard deviation ∆ n = [ K − (cid:80) Ki =1 ( n i − K − (cid:80) Ki =1 n i ) ] / ,where different realizations of disorder are indexed by i = 1 , , ..., K . Note the peak indicating the anomalousbroadening of the distribution around the critical point ν = ν ∗ . The situation is fundamentally different in theconventional decaying Burgers turbulence where the ini-tial data have zero average, which infinitely delays theemergence of scaling. . . . . (a) (b) (c) ν ∆ n √ N N = 10 N = 10 N = 10
890 1090540 740 0 200
Figure 5. Time (rigidity) evolution of the normalized stan-dard deviation for the number of shocks. The statistics wasobtained from K = 1000 realizations of the quenched disorderwith ρ = 3 , and λ = 1 . Inset plots: (a) ν = 0 . (b) ν = 1 ; (c) ν = 6 . So far we were assuming that the rigidity measure ν is finite as N → ∞ . A broader class of elastic envi-ronments can be modeled if we assume that κ ∼ N α ,with ≤ α ≤ . For instance, if the load is transmittedthrough a surface of a 3D body we have α = 2 / and ν ∼ N / . In this setting, small systems would be neces-sarily ductile, while brittle behavior would dominate inthe thermodynamic limit. At a given disorder, the scal-ing will be then seen in a window of system sizes, whilethe (percolation type) critical regime will emerge only atinfinite size and infinite disorder [6, 7, 71].To conclude, we used an analytically transparentmodel to quantify the role of system’s rigidity (globalconnectivity) as a control parameter for the transitionfrom brittle to ductile failure. We showed that this tran-sition can be associated with the crossover from spinodalto classical criticality, generating, in finite size systems, ascaling region with non-universal exponents. Such behav-ior is generic for a broad class of systems, encompassingfracture, plasticity, structural phase transitions, and nowwe established a new link to fluid turbulence.The authors are grateful to R. Garcia-Garcia, K. Dah-men, and M. Mungan for helpful discussions. H.B.R. wassupported by a PhD fellowship from Ecole Polytechnique;L. T. was supported by the grant ANR-10-IDEX-0001-02PSL. ∗ [email protected] † [email protected][1] H. Herrmann and S. Roux, Statistical Models for theFracture of Disordered Media , Random Materials andProcesses (Elsevier Science, 2014).[2] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, Ad-vances in Physics , 349 (2006).[3] I. Procaccia, C. Rainone, and M. Singh, Phys. Rev. E , 032907 (2017).[4] F.-J. Pérez-Reche, L. Truskinovsky, and G. Zanzotto,Phys. Rev. Lett. , 230601 (2008).[5] D. S. Fisher, K. Dahmen, S. Ramanathan, and Y. Ben-Zion, Phys. Rev. Lett. , 4885 (1997).[6] R. Toussaint and A. Hansen, Phys. Rev. E , 046103(2006).[7] A. Shekhawat, S. Zapperi, and J. P. Sethna, Phys. Rev.Lett. , 185505 (2013).[8] S. Roy, S. Biswas, and P. Ray, Phys. Rev. E , 063003(2017).[9] M. M. Driscoll, B. G.-g. Chen, T. H. Beuman, S. Ulrich,S. R. Nagel, and V. Vitelli, Proceedings of the NationalAcademy of Sciences , 10813 (2016).[10] L. Zhang, D. Z. Rocklin, L. M. Sander, and X. Mao,Phys. Rev. Materials , 052602 (2017).[11] C. P. Goodrich, A. J. Liu, and S. R. Nagel, NaturePhysics , 578 EP (2014).[12] R. L. B. Selinger, Z.-G. Wang, W. M. Gelbart, andA. Ben-Shaul, Phys. Rev. A , 4396 (1991).[13] J. B. Rundle and W. Klein, Phys. Rev. Lett. , 171(1989).[14] S. Zapperi, P. Ray, H. E. Stanley, and A. Vespignani,Phys. Rev. Lett. , 1408 (1997).[15] A. Wisitsorasak and P. G. Wolynes, Proceedings of theNational Academy of Sciences , 16068 (2012).[16] Y. Moreno, J. B. Gómez, and A. F. Pacheco, Phys. Rev. Lett. , 2865 (2000).[17] J. V. Andersen, D. Sornette, and K.-t. Leung, Phys. Rev.Lett. , 2140 (1997).[18] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tar-jus, Proceedings of the National Academy of Sciences , 6656 (2018).[19] M. Sahimi and S. Arbabi, Phys. Rev. Lett. , 608(1992).[20] S. Roux, A. Hansen, H. Herrmann, and E. Guyon, Jour-nal of Statistical Physics , 237 (1988).[21] A. Hansen and J. Schmittbuhl, Phys. Rev. Lett. ,045504 (2003).[22] R. Toussaint and S. R. Pride, Phys. Rev. E , 046127(2005).[23] A. A. Moreira, C. L. N. Oliveira, A. Hansen, N. A. M.Araújo, H. J. Herrmann, and J. S. Andrade, Phys. Rev.Lett. , 255701 (2012).[24] M. Popović, T. W. J. de Geus, and M. Wyart, Phys.Rev. E , 040901 (2018).[25] B. Kahng, G. G. Batrouni, S. Redner, L. de Arcangelis,and H. J. Herrmann, Phys. Rev. B , 7625 (1988).[26] F. Yuan and L. Huang, Scientific Reports , 5035 EP(2014).[27] J. Liu, Z. Zhao, W. Wang, J. W. Mays, and S.-Q. Wang,Journal of Polymer Science Part B: Polymer Physics ,758 (2019).[28] D. Şopu, A. Foroughi, M. Stoica, and J. Eckert, NanoLetters , 4467 (2016).[29] G. Subhash, Q. Liu, and X.-L. Gao, International Jour-nal of Impact Engineering , 1113 (2006).[30] J. Zhao, X.-T. Feng, X.-W. Zhang, Y. Zhang, Y.-Y.Zhou, and C.-X. Yang, Engineering Geology , 160(2018).[31] M. Selezneva, Y. Swolfs, A. Katalagarianakis,T. Ichikawa, N. Hirano, I. Taketa, T. Karaki, I. Verpoest,and L. Gorbatikh, Composites Part A: Applied Scienceand Manufacturing , 20 (2018).[32] D. Krajcinovic, S. Mastilovic, and M. Vujosevic, Mecca-nica , 363 (1998).[33] R. Christensen, Z. Li, and H. Gao, Proceedings of theRoyal Society A: Mathematical, Physical and Engineer-ing Sciences , 20180361 (2018).[34] S. Papanikolaou, J. Thibault, C. Woodward, P. Shan-thraj, and F. Roters, arXiv preprint arXiv:1707.04332(2017).[35] E. Berthier, J. E. Kollmer, S. E. Henkes, K. Liu,J. M. Schwarz, and K. E. Daniels, arXiv preprintarXiv:1812.07466 (2018).[36] A. Hansen, P. Hemmer, and S. Pradhan, The FiberBundle Model: Modeling Failure in Materials , Statisti-cal Physics of Fracture and Breakdown (Wiley, 2015).[37] F. T. Peirce, Journal of the Textile Industry , 355(1926).[38] H. Nechad, A. Helmstetter, R. E. Guerjouma, andD. Sornette, Journal of the Mechanics and Physics ofSolids , 1099 (2005).[39] Didier Sornette, J. Phys. I France , 2089 (1992).[40] P. C. Hemmer and A. Hansen, Journal of Applied Me-chanics , 909 (1992).[41] A. Delaplace, S. Roux, and G. P. jaudier Cabot, Interna-tional Journal of Solids and Structures , 1403 (1999).[42] I. Balog, M. Tissier, and G. Tarjus, Phys. Rev. B ,104201 (2014).[43] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Man- ning, Proceedings of the National Academy of Sciences , 6560 (2019).[44] J. R. Rice, Fracture: an advanced treatise , 191 (1968).[45] J.-P. Bouchaud and M. Mézard, Journal of Physics A:Mathematical and General , 7997 (1997).[46] P. Le Doussal and K. J. Wiese, Phys. Rev. E , 051106(2009).[47] J. Bec and K. Khanin, Physics Reports , 1 (2007).[48] A. Basu and B. K. Chakrabarti, Philosophical Transac-tions of the Royal Society A: Mathematical, Physical andEngineering Sciences , 20170387 (2019).[49] If the breakable elements are composed of sub-partslinked in series and if the failure is associated with break-ing of the weakest sub-part, the Weibull distributionemerges rigorously in the thermodynamic limit as thedistribution for the breaking threshold of the whole sys-tem. Here we assume that the distribution of thresholdsfor the sub-parts has a compact support.[50] See Supplemental Material at [URL will be inserted bypublisher] for the study of the structure of the metastablestates, the computation of the statistics of avalanche dis-tribution, mapping on the RFIM and the reduction toBurgers model in the case of finite N .[51] M. F. J. Vermeulen, A. Bose, C. Storm, and W. G.Ellenbroek, Phys. Rev. E , 053003 (2017).[52] H. Crapo, Structural topology, 1979, núm. 1 (1979).[53] J. Z. Kim, Z. Lu, S. H. Strogatz, and D. S. Bassett,Nature Physics , 714 (2019).[54] S. Zapperi, K. B. Lauritsen, and H. E. Stanley, Phys.Rev. Lett. , 4071 (1995).[55] Our asymptotic analysis is valid for arbitrary disorderas long as the function h ( x ) has a minimum. Some longtailed disorders can modify the behavior of the system,for instance, the distribution of thresholds F ( x ) = 0 , for x ≤ , and F ( x ) = 1 − / √ x , for x > , leads to thefunction h ( x ) without a minimum.[56] S. di Santo, R. Burioni, A. Vezzani, and M. A. Muñoz,Phys. Rev. Lett. , 240601 (2016).[57] Here SNAP corresponds to supercritical regimes whichinclude isolated system size avalanches, POP de-scribes subcritical regimes with uncorrelated eventsdistributed around an average size, and the termCRAKLE/crackiling is used for critical regimes withscale-free avalanches [4, 72].[58] R. C. Hidalgo, K. Kovács, I. Pagonabarraga, and F. Kun,EPL (Europhysics Letters) , 54005 (2008).[59] R. C. Hidalgo, F. Kun, K. Kovács, and I. Pagonabarraga,Phys. Rev. E , 051108 (2009).[60] K. Dahmen and J. P. Sethna, Phys. Rev. B , 14872(1996).[61] A. Clauset, C. R. Shalizi, and M. E. J. Newman, SIAMReview , 661 (2009).[62] J. Weiss, W. B. Rhouma, T. Richeton, S. Dechanel,F. Louchet, and L. Truskinovsky, Phys. Rev. Lett. ,105504 (2015).[63] Y. Xu, A. G. Borrego, A. Planes, X. Ding, and E. Vives,Phys. Rev. E , 033001 (2019).[64] G. Sparks and R. Maaß, Phys. Rev. Materials , 120601(2018).[65] F. J. Pérez-Reche and E. Vives, Phys. Rev. B , 214422(2004).[66] A. Maritan, M. Cieplak, M. R. Swift, and J. R. Banavar,Phys. Rev. Lett. , 946 (1994).[67] Y. Liu and K. A. Dahmen, Phys. Rev. E , 061124 (2009).[68] I. Balog, G. Tarjus, and M. Tissier, Phys. Rev. B ,094204 (2018).[69] L. C. Evans, Partial Differential Equations , Graduatestudies in mathematics, Vol. 19 (American Mathemat-ical Society, 2010).[70] M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. , 889 (1986).[71] Z. Olami, H. J. S. Feder, and K. Christensen, Phys. Rev.Lett. , 1244 (1992).[72] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature , 242 (2001). SUPPLEMENTAL MATERIAL
To obtain the avalanche distribution in our generalizedFBM problem with controlled length , we follow the gen-eral methodology largely developed by Hansen and col-laborators in their studies of the classical FBM problemwhich implies control of the force [1–5].
Metastable states.
First, we use the condition ∂ X H =0 to obtain X ( x , ε ) = λ +Λ (cid:16) Λ ε + λ N (cid:80) Ni =1 x i (cid:17) , and thecondition ∂ x i H = 0 to obtain u (cid:48) ( x i ) = λ ( X − x i ) . In viewof permutational invariance, we can characterize the mi-croscopic state by the number of broken bonds, k , whichgives ˆ X ( k, ε ) = (1 + λ )Λ ελ (1 − k/N ) + λ Λ + Λ . (5)For the attached links we have ˆ x ( k, ε ) = λ Λ ελ (1 − k/N ) + λ Λ + Λ , (6)and for the broken links ˆ x ( k, ε ) = (1 + λ )Λ ελ (1 − k/N ) + λ Λ + Λ . (7)The energy of the equilibrium configurations can be writ-ten as H ( k, ε ) = a k ε + S k , (8)where a k = 12 λ Λ( N − k ) λ ( N − k ) + N ( λ Λ + Λ) , and S k is the en-ergy of the broken bonds. If ¯ x i , i = 1 , ..., N is the orderedsequence of failure thresholds, ¯ x ≤ ¯ x ≤ · · · ≤ ¯ x N , wecan write S k = 1 N k (cid:88) i =1 ¯ x i , and S = 0 . We observe that a k is a (strictly) monotonically decreasing sequence while S k is a (strictly) monotonically increasing sequence.The stress-strain relation for a microscopic state char-acterized by the parameter k is σ ( k, ε ) = ∂ H ( k, ε ) ∂ε = λ Λ( N − k ) ελ ( N − k ) + N ( λ Λ + Λ) . (9) Each value of k defines an equilibrium branch extend-ing between the two limits induced by the inequalities ˆ x ( k, ε ) < ¯ x k and ˆ x ( k, ε ) > ¯ x k . For the failure thresh-olds we can then write ε fk = λ + 1 λ (cid:34)(cid:18) − kN (cid:19) ν + 1 (cid:35) ¯ x k , (10)where ≤ k < N . Similar expressions can be obtainedfor the rebuilding thresholds ε rk = (cid:34)(cid:18) − kN (cid:19) ν + 1 (cid:35) ¯ x k , (11)where < k ≤ N . The ensuing equilibrium branches arerepresented by the gray lines in Fig. 1 (b, c) in the maintext.To analyze their (local) stability, we need to study thepositive definiteness of the Hessian matrix for the energy H ( x , X ) M = M . . . − λ . . . . . . ... ...... . . . . . . ... . . . M N − λ − λ . . . . . . − λ N ( λ + Λ) , (12)where M i is either λ + 1 , for ≤ i < N − k , or λ , for N − k ≤ i ≤ N . The sufficient condition for stabilityis that all the principal minors of M are positive. Thefirst N minors are just the product of diagonal thermsare therefore always positive. The last principal minor,the determinant det( M ) = N (cid:89) i =1 M i N (cid:88) i =1 (cid:32) λ + Λ − λ M i (cid:33) . (13)is also positive implying stability of the obtained equi-librium configurations; the unstable configurations mustcontain at least one element in the spinodal state repre-sented in our model by a single point. Equilibrium (global minimum) path.
For large N wecan write S k = 1 N k (cid:88) i =0 ¯ x i ≈ (cid:90) ¯ x k ¯ x x f ( x ) dx, (14)where we used the fact that for ordered distributions wecan use the approximation k/N ∼ F (¯ x k ) [6]. We canthen write the continuous approximation of the discreteenergy in the form H ( x, ε ) = λ Λ(1 − F ( x )) λ (1 − F ( x )) + Λ( λ + 1) ε (cid:90) x f ( x (cid:48) ) x (cid:48) dx (cid:48) . (15)Using the equilibrium condition ∂ H ( x, ε ) /∂x = 0 , andapplying it for the discrete values ¯ x k , we obtain ε gk = 1 √ Λ ν (cid:34)(cid:18) − kN (cid:19) ν + 1 (cid:35) ¯ x k . (16)Note that the three formulas (10), (11) and (16) are dif-ferent only by constant multipliers. Out-of-equilibrium (zero viscosity limit) path.
Eachmicroscopic configuration characterized by parameter k exists in an extended domain of the loading parameter ε between the failure strain ε fk and the rebinding strain ε rk .At large N , we can use the approximation ¯ ε f ( x ) = λ + 1 λ (cid:104)(cid:0) − F ( x ) (cid:1) ν + 1 (cid:105) x. (17)and ¯ σ f ( x ) = [1 − F ( x )] x. Similarly, along the reversepath, ¯ ε r ( x ) = (cid:104)(cid:0) − F ( x ) (cid:1) ν + 1 (cid:105) x. (18)and ¯ σ r ( x ) = λλ +1 [1 − F ( x )] x. Both, equilibrium and outof equilibrium (averaged) stress-strain relations are illus-trated in Fig. 6. . . . εσ Brittle 0 1 2 3 ε Critical 0 1 2 3 ε Ductile0 0 . . x ¯ ε . . x . . x Figure 6. First row: avergaed stress-strain relations , secondrow: averaged strain dependence on the internal variable x ;blue (red) curves correspond to the loading (unloading) out-of-equilibrium paths; black curves correspond to the equilib-rium (global minimum) path. Brittle to ductile transition.
It is easier to see if thesystem is brittle if we consider the out-of-equilibrium(marginal stability) path, even though the actual ductil-ity threshold would be the same if we consider the globalminimum path. All we need to check is the conditionthat the curve ¯ ε ( x ) has a local maximum, which reads [1 − F ( x c )] − f ( x c ) x c + ν − = 0 . To locate the brit-tle to ductile transition we need to find the inflectionpoint on the curve ¯ ε ( x ) characterized by the condition − f ( x c ) − f (cid:48) ( x c ) x c = 0 .In the case of Weibull distribution, we obtain from the first of these two conditions x c = (cid:34) ρ − W (cid:18) − exp (1 /ρ ) ρν (cid:19)(cid:35) /ρ , (19)where W ( x ) is the Lambert function, defined throughthe equation x = W ( x ) e W ( x ) . Then the second condi-tion gives ν = e ρ +1 /ρ, which delineates the boundarybetween brittle and ductile regimes. Statistics of avalanches.
We first compute theavalanche distribution for the case of the out-of-equilibrium loading path; it will be clear that the sameprocedure can be adapted for the out-of-equilibrium un-loading path and for the reversible equilibrium paths.For an avalanche of size s to take place along the out-of-equilibrium loading path and be associated with thefailure of k th (in strength) spring, we must have ε k + j ≤ ε k , for j = 1 , , . . . , s − , and ε k + s > ε k , which, follow-ing [3], we called the forward condition ; to secure that ε k is larger than previous thresholds, we must also requirethat ε j ≤ ε k , for all j < k , which we call the backwardscondition . Given that the rebinding sequence for the un-loading out-of-equilibrium path is ε rk = λλ +1 ε k and for the equilibrium path is ε gk = (cid:113) λλ +1 ε k , the avalanche condi-tion in those two cases and the ensuing avalanche statis-tics will be the same as in the case of out-of-equilibriumloading path, so it is sufficient to deal with this case only.Since we are interested in the asymptotics for theavalanche distribution at large N , we assume that s (cid:28) N . Using the ordered thresholds ¯ x i , and Eq. (10) for thesequence ε k , we can then rewrite ε k + j ≷ ε k in the form ¯ x k + j ≷ ¯ x k (cid:18) jN − k − j + N ν − (cid:19) . (20)Defining δ k = ¯ x k N − k + N ν − , and using the assumptionthat j (cid:28) N − k , we can simplify these relation further ¯ x k + j ≷ ¯ x k + jδ k . (21)Note next that breaking of one spring at the elongation ε k , corresponding to a threshold ¯ x k = x , raises the loadon the remaining fiber by δ k . The average number offibers that breaks as a result of this load increase is equalto the number of thresholds in the interval ( x, x + δ k ) ,which is N f ( x ) δ k . Thus, the average number of fibersbreaking as a result of the failure of the k th fiber is, g ( x ) = f ( x ) x − F ( x ) + ν − , (22)where we again used the approximation k/N ∼ F ( x ) [6].For an avalanche of size s , the increase in load will beapproximately sδ k , which leads to g ( x ) s broken springs.The (forward) probability that the additional s − springsbreak is then given by a Poisson distribution with the rate g ( x ) s , ˜ p f ( s, x ) = ( g ( x ) s ) s − ( s − e − g ( x ) s . (23)To complete this expression, we still need to secure thecondition stating that all the s − inequalities ¯ x k +1
We now focus on the tail of thedistribution p ( s ) assuming that N → ∞ . We use thesaddle-point approximation, which implies that the maincontribution to the integral will come from the vicinity of x = x , where the function h ( x ) = g ( x ) − ln g ( x ) reachesits global minimum. To find x , we need to solve theequation h (cid:48) ( x ) = g (cid:48) ( x ) g ( x ) ( g ( x ) −
1) = 0 . There are threepossibilities,1. g ( x ) (cid:54) = 1 and g (cid:48) ( x ) = 0 (ductile regime),2. g ( x ) = 1 and g (cid:48) ( x ) = 0 (critical regime)3. g ( x ) = 1 and g (cid:48) ( x ) (cid:54) = 0 (brittle regime).If g ( x ) (cid:54) = 1 , and g (cid:48) ( x ) = 0 , we can write, h ( x ) ≈ g ( x ) − ln g ( x ) + g (cid:48)(cid:48) ( x )2 g ( x ) ( g ( x ) − x − x ) . Then usingthe saddle-point approximation in (30), and applying theStirling approximations s ! ≈ s s e − s √ πs , we obtain p ( s ) = s s − s ! e − sh ( x ) φ ( x ) (cid:115) πs | h (cid:48)(cid:48) ( x ) |∼ s − e − s ( h ( x ) − . (31)When simultaneously g ( x ) = 1 and g (cid:48) ( x ) = 0 wehave h (cid:48)(cid:48) ( x ) = 0 , and h (cid:48)(cid:48)(cid:48) ( x ) = 0 ; therefore the Taylorexpansion is h ( x ) ≈ g (cid:48)(cid:48) ( x ) ( x − x ) . We can alsowrite φ ( x ) ≈ − f ( x ) g (cid:48)(cid:48) ( x )2 ( x − x ) , which allows us tore-write the integral (30) in the form, p ( s ) = s s − e − s s ! (cid:90) x − f ( x ) g (cid:48)(cid:48) ( x )( x − x ) × e − s g (cid:48)(cid:48) ( x ( x − x ) dx. (32)Computing the integral explicitly and using Stirling’s ap-proximation we obtain p ( s ) ∼ s − / . In the brittle regime we need to consider separatelyequilibrium and out of equilibrium paths.Consider first the out-of-equilibrium path. We need toexpand the function h ( x ) = g ( x ) − ln g ( x ) up to secondorder to obtain h ( x ) ≈ g (cid:48) ( x )2 ( x − x ) . We canalso expand φ ( x ) to obtain φ ( x ) ≈ − g (cid:48) ( x ) f ( x )( x − x ) .These expansions allow us to approximate the integral(30) by p ( s ) = s s − s ! e − s (cid:90) x g (cid:48) ( x ) f ( x )( x − x ) e − s g (cid:48) ( x )22 ( x − x ) dx. (33)Along the out-of -equilibrium path, the avalanches arecounted up to x = x ; and if we compute the integralexplicitly, and use the Stirling approximations, we obtain p ( s ) ∼ s − / . Consider now the equilibrium path. The actual equi-librium SNAP event takes place at some x ∗ < x , givenby the Maxwell construction. The counting of avalanchesshould be then performed only up to the point x ∗ , andin the integral (30), we must put x c = x ∗ . The function h ( x ) will attain its minimum in the boundary point x ∗ ,which is the upper limit of integration. In such case, thefollowing asymptotic representation holds at N → ∞ [7] (cid:90) x sup x inf e − Nh ( x ) dx → e Nh ( x ∗ ) N h (cid:48) ( x ∗ ) (34)This allows to write, p ( s ) ∼ s − / e − s (1 − h ( x ∗ )) . Mapping on RFIM.
Using the condition ∂ X H = 0 ,we obtain X ( x , ε ) = λ +Λ (cid:16) Λ ε + λ N (cid:80) Ni =1 x i (cid:17) . If we sub-stitute this expression back into H we obtain H = − N (cid:88) i,j Jx i x j − N (cid:88) i [ Hx i − v i ( x i )] , (35)where J = λ λ +Λ) , H = λ Λ ελ +Λ , and v i ( x i ) = u i ( x i ) + x i + λ Λ ε λ + Λ) . Initial condition for the Burgers equation.
In the caseof finite N , the equilibrium condition ∂ x i H = 0 allows us to write H ( X, ε ) = 1 N N (cid:88) i =1 e i ( X ) + Λ2 ( ε − X ) . Here, two metastable branches e i = λλ +1 X Θ( l i − λλ +1 X ) + l i Θ( X − l i ) are defined in each interval X ∈ [ l i , λ +1 λ l i ] . If we choose the branch with the minimal en-ergy, the remaining problem reduces to finding ˜ H ( ε, ν ) = min X ∈ R (cid:26) ν ( ε − X ) + q ( X ) (cid:27) , where q ( X ) = N (cid:80) Ni =1 X Θ( l i − (cid:113) λλ +1 X )+ λ +1 λ l i Θ( X − (cid:113) λλ +1 l i ) . The initial data for the associated Burgersequation are σ = ∂ ε q = 1 N N (cid:88) i =1 ε Θ (cid:32) l i − (cid:114) λλ + 1 ε (cid:33) . In the limit N → ∞ we have N (cid:80) Ni =1 Θ( l i − (cid:113) λλ +1 X ) ∼ (cid:82) ∞ (cid:113) λλ +1 X f ( l ) dl and N (cid:80) Ni =1 l i Θ( X − (cid:113) λλ +1 l i ) ∼ (cid:82) (cid:113) λλ +1 X f ( l ) l dl. Then, in this limit, ˜ H ( ε, ν ) = min X ∈ R (cid:26) ν ( ε − X ) + q ∞ ( X ) (cid:27) , where q ∞ ( X ) = (cid:113) λλ +1 (cid:82) X f ( (cid:113) λλ +1 X (cid:48) )( X (cid:48) / dX (cid:48) +[1 − F ( (cid:113) λλ +1 X )]( X / . The initial condition for the associ-ated Burgers equation is σ ( ε ) = ε [1 − F ( (cid:113) λλ +1 ε )] . ∗ [email protected] † [email protected][1] A. Hansen, P. Hemmer, and S. Pradhan, The FiberBundle Model: Modeling Failure in Materials , Statisti-cal Physics of Fracture and Breakdown (Wiley, 2015).[2] A. Hansen and P. C. Hemmer,
Criticality in fracture: theburst distribution , Tech. Rep. T-TPS-94-8. Trondheim-TPS-8-1994 (Trondheim TU. Inst. Phys., Trondheim,1994).[3] P. C. Hemmer and A. Hansen, Journal of Applied Me-chanics , 909 (1992).[4] A. H. M. Kloster and P. C. Hemmer, Physical Review E , 2615 (1997).[5] S. Pradhan, A. Hansen, and B. K. Chakrabarti, Reviewsof Modern Physics , 499 (2010).[6] B. Arnold, N. Balakrishnan, and H. Nagaraja, A FirstCourse in Order Statistics , Classics in Applied Mathe-matics (Society for Industrial and Applied Mathematics,1992). [7] N. de Bruijn, Asymptotic Methods in Analysis , DoverBooks on Mathematics (Dover Publications, 2014).[8] A. Clauset, C. R. Shalizi, and M. E. J. Newman, SIAMReview , 661 (2009). [9] M. E. Newman, Contemporary physics , 323 (2005).[10] J. Baró and E. Vives, Phys. Rev. E , 066121 (2012).[11] M. L. Goldstein, S. A. Morris, and G. G. Yen, The Euro-pean Physical Journal B - Condensed Matter and Com-plex Systems41