Inertia and universality of avalanche statistics: The case of slowly deformed amorphous solids

By means of a finite elements technique we solve numerically the dynamics of an amorphous solid under deformation in the quasistatic driving limit. We study the noise statistics of the stress-strain signal in the steady state plastic flow, focusing on systems with low internal dissipation. We analyze the distributions of avalanche sizes and durations and the density of shear transformations when varying the damping strength. In contrast to avalanches in the overdamped case, dominated by the yielding point universal exponents, inertial avalanches are controlled by a non-universal damping dependent feedback mechanism; eventually turning negligible the role of correlations. Still, some general properties of avalanches persist and new scaling relations can be proposed.


I. INTRODUCTION
Avalanche behavior is ubiquitous in nature. Several systems respond to a slow constant driving with intermittent dynamics. Examples are: Barkhausen noise in ferromagnets [1], stick-slip dynamics in the plasticity of solids [2], earthquakes [3], creep of magnetic domain walls [4], serrated stress of driven foams [5] and crack propagation [6]. The common phenomenon is that different regions of a heterogeneous system are loaded towards instability thresholds; a region that first gets unstable yields, destabilizes others, and this goes on producing an "avalanche". Recently, such avalanche dynamics has been evidenced in the time series of the stress tensor in deformation experiments of amorphous systems, such as grains, foams or metallic glasses [7][8][9], and has attracted considerable theoretical interest [10][11][12].
Global quantities linked to such collective behavior are usually power law distributed and allow for the understanding of the system parameter dependencies in terms of scaling functions. This is a gift from the self-organized criticality (SOC) paradigm, in which stationary states with "critical" behavior are attractors of the dynamics [13][14][15] and dominate the scenario as soon as a balance between branching and killing probabilities is fulfilled. Among the factors that may break-down this balance in SOC we count inertia [16][17][18][19].
Various works have addressed inertial-like effects in the context of SOC, like sand-pile models with threshold weakening [17,19] or depinning-like models that include stress overshoots [18,20] or softening [21,22]. Moreover, inertia has been explicitly considered on the Burridge-Knopoff model [23][24][25] in the context of seismic faults. In all cases, avalanche size distributions were found to deviate from the critical power-law scaling. It should be said, nevertheless, that none of this models was intended to illustrate the effect of inertia in the deformation of amorphous solids.
In this respect, recent studies on atomistic simulations of glasses under deformation [26,27] argue, on the contrary, that inertial effects drive the system to a "new underdamped universality class", rather than taking it away from criticality. On the other hand, the same kind of atomistic approach [28], as well as more coarse grained method [29], have signaled a strong contrast of the finiteshear-rate rheology between the overdamped and the underdamped cases; with no signs of universal behavior in the latter. Strongly inertial underdamped systems tend to produce, in particular, a non-monotonic flow curve and the associated localization of the deformation [28,29].
In this work, we study the influence of inertia in the statistics of avalanches at the yielding point. By means of a finite-elements based elasto-plastic model we analyze the stress time series of a sheared system and relate the stress-drops there observed with avalanches of several sites yielding in a collective process. At high damping we recover the critical avalanche statistics found in models of amorphous solids with overdamped dynamics [10][11][12]30]. When damping is decreased, we observe that the avalanche statistics smoothly deviates from the critical distribution. Even when a power-law shape is observed for a range of avalanche sizes, the distribution ceases to be scale-free, developing a bump at large values set by the system size. Also the exponent of the power-law regime systematically changes as damping is decreased, indicating a departure from universal behavior.
We understand the avalanche statistics of systems with low dissipation as a combination of two distinct kinds of events. On one hand, avalanches dominated by inertial effects, typically large in size, frequently system-spanning and with a rather well defined size, that populates the bump or peak of the distributions. On the other hand, smaller and more localized avalanches, which mostly contribute to the power-law regimes of the distributions, remain controlled by the critical point attractor of the overdamped dynamics and show little influence of inertia.
We analyze the emergence of this dichotomy by a careful study of the finite system-size scaling of the avalanche distributions, their dependence on damping, and the changes on avalanches geometry. Within this new scenario, we justify the observed alteration by inertia of the probability gap for the density of shear transformations and propose a scaling relation that links it to the avalanches geometry.

II. MODEL AND PROTOCOL
We perform a finite-element-based simulation of amorphous systems in d = 2 as in [29]. Consider a continuous medium with a displacement field u(r) from an equilibrium configuration with position r. The equation of motion for the displacement field u(r, t) in such a continuum reads where ρ is the uniform mass density and σ is the stress. Conventionally, the state of a stress tensor is expressed as σ = −pI + σ dev with p = −tr(σ)/d the hydrostatic pressure, σ dev the deviatoric stress, and I the identity tensor. We further define the effective shear stress as σ = ( 1 2 σ dev : σ dev ) 1 2 . The contribution to the stress at r is assumed to depend only on the gradients of u(r) and/oru(r), more accurately on their symmetric parts and/or˙ . Similarly, the deformation rate may be also decomposed as˙ =˙ v I +˙ dev with˙ v = tr(˙ ) the volumetric strain rate and˙ dev the deviatioric strain rate. In this case, the effective shear rate can be expressed aṡ = ( 1 2˙ dev :˙ dev ) 1 2 . Each finite-size element represents microscopically rearranging zones (as in real particulate packings) at a coarse-grained level. In their rigid phase, these zones are modeled by a Kelvin-type solid, while past a yielding threshold σ y (r, t), a simple Newtonian fluid is employed. Given this qualitative picture, the constitutive equations are determined by with K and µ the bulk and shear moduli, respectively. The first terms on the right-hand side of the equations mimic the elastic contribution, while the second term, only present for the deviatoric piece, represents viscous dissipation with the viscosity coefficient η. Here n(r, t) is equal to one during the fluid phase -whose life time is limited by a threshold τ p -and zero otherwise. We set τ p = τ min v with τ min v defined here later on. Fixing this time gives a distribution of (plastic) strains in the range %0.5 − 1. Realistically, Eshelby zones undergo strong shear deformations, rather than dilation, to relax the stress, and are, therefore, nearly incompressible. Incompressibility is enforced in the elastic regime by setting K/µ ≈ 10. As in [29], yield stresses are drawn randomly from an exponential distribution with a mean valueσ y and a lower cut-off σ min y , reminiscent of structural disorder in glassy dynamics. Also, we dynamically assign a new random threshold to each element after yielding.
Having defined a proper set of constitutive equations, an irregular set of grid-points and linear plane-strain elasto-plastic triangular elements is then employed in a L × L periodic cell with typical grid-size of h to discretize Eq.(1) in space (see Fig. 1). Periodic boundaries are implemented by carefully assigning (based on images at the "borders" of the simulation cell) a lists of neighboring nodes for each node of the irregular grid, list that remains fix during the simulation.
An area-preserving (simple) shear is implemented to drive the system in a quasi-static manner; that is, an infinitesimal strain-step [31] ∆γ ≈ 10 −5 is initially applied to the simulation box each time, followed immediately by a stress quench thanks to the dissipative terms present in both Kelvin and Newtonian descriptions. The quench runs until the maximum net force on any grid-point is less than 10 −6 times the average force. The damping rate for the viscous term is τ −1 d = ηq 2 /ρ, where the vibrational frequency is τ −1 v = c s q with c s = µ/ρ the transverse wave speed. Here q is the wave number (spatial frequency). The largest q corresponds to the inverse element size h −1 , from which the the shortest vibrational time-scale τ min v = h/c s . The numerical time integration during quench is performed by means of the velocity Verlet algorithm with ∆t min(τ d , τ v ). The dimensionless quantity Γ = τ −1 d /τ −1 v , called damping ratio or damping coefficient hereafter, quantifies the relative impact of dissipation. We expect to capture an overdamped dynamics with Γ 1, while Γ ≤ 1 should in principle lead to an underdamped regime.

III. RESULTS
As for MD simulations and scalar elasto-plastic models, imposing a large enough shear deformation results, after a transient regime, into a steady-state plastic flow. The standard observable is the mean stress σ xy (γ) averaged across the sample. The stress signal is characterized by a series of elastic loading periods, where the global stress grows linearly with the applied deformation, interrupted by stress drops ∆σ, which are due to plastic activity and mark an abrupt release of stored elastic energy. Fig. 2(a) illustrates this common scenario in both high and low damping simulations. A clear qualitative difference between the high and low damping cases is evident from the data: while overdamped fluctuations tend to span a broad range of scales (absence of a characteristic size), inertial signals seem to be better described by a typical, large, characteristic stress drop, appearing in a quasi-periodic fashion. The averaged stress drop value ∆σ as a function of the inverse damping coefficient Γ −1 (bigger the more inertial is the system) is shown in Fig. 2(b). At low values of Γ −1 the mean stress drop saturates, which can be used to define the overdamped limit.

A. Avalanche size distributions
The magnitude of the stress drop reflects the number of sites involved in a correlated sequence of yielding events or avalanche. Therefore, in agreement with earlier studies [26,27], we simply define an extensive avalanche size as S = σ xy ∆σL d . By using the mean stress value σ xy as a multiplicative factor that depends on Γ, S acquires (within a constant factor) the dimensions of an energy drop, that better quantifies the collective behavior and allows for a comparison among systems with different damping. Figure 3 shows distributions of avalanche sizes S for different inverse damping coefficients at a fixed system size L = 40 built from a statistics of 10000 events. We first notice that an overdamped system (Γ −1 = 0.2) displays, after a characteristic lower cutoff, a power-law decay in S that is later on suppressed by an exponential decay. In fact, based on previous results, we expect to observe in this limit a behavior P S ∼ S −τ G(S/S c ) with G(·) a rapidly decaying function and the cutoff set by S c . The inset of Fig. 3 shows curves for such a high damping at different system sizes rescaled as P S S τ c vs S/S c , where we introduce the size dependent cutoff S c (L) ∝ L d f with d f the fractal dimension of the avalanches. A collapse onto a unique master curve is obtained for large S. The power-law fitting and collapse found are characterized by the values τ = 1.27 ± 0.05 and d f = 0.95 ± 0.05, in agreement with previous results of various elasto-plastic models [10,12,30] and molecular dynamics simulations in the overdamped limit [12,26,27].
When the damping is decreased and inertia starts to be relevant, two main new features become apparent: On the one hand the P S distributions deviate from a pure power-law shape; a shoulder develops at large values of S and evolves into a clear local maximum at very low damping. On the other hand, an apparently robust powerlaw regime survives at smaller values of S; nevertheless the exponent τ characterizing the decay of P S systematically deviates from its overdamped value as the damping is decreased, reaching values of τ 1.5 for the lower damping displayed, with no signs of saturation. Furthermore, a closer inspection of the distributions obtained at low damping evidences the existence of an intermediate regime in which a probability depletion is produced; simultaneously showing how the inertial peak at large S degrades the free-scale regime and suggesting a separation between two different kind of events. Let us recall that, indeed, a differentiation of two groups of events has been proposed for the Burridge-Knopoff model [23][24][25].
In our case, this distinctive feature of the avalanche size distribution can be better understood by analyzing the distribution of avalanche durations, P T . Figure 4 shows this distribution in units of τ v , which is the relevant microscopic time for inertial systems [28]; similar to P S , the distribution shows at large damping a power-law regime, exponentially suppressed for long durations. As damping is decreased, a peak at large values of T appears. A vertical dashed line marks the size dependent value T L = L/c s . Figure 4 suggests a simple interpretation of the appearance of the inertial peak in the distributions. At the local level, the immediate consequence of the lack of dissipation is a stronger effect of a yield event in its surrounding, that could be seen as an effective reduction in the local thresholds. Nevertheless, the most notorious inertial feature, as we understand now, comes from a long-range action. When energy is not rapidly dissipated it keeps traveling around the lattice in the form of shear elastic waves. Inertial effects become dominant when an avalanche duration is long enough such that the avalanche is able to reinforce itself through its own elastic waves (traveling around the system or possibly being reflected at the boundaries of a non-periodic system). We expect this time threshold to be proportional to T L , the time needed by the elastic waves to propagate across one of the system's main axis. This is confirmed by the behavior of P T at different damping coefficients. In the still overdamped case (Γ −1 = 0.5, putting the limit at Γ = 1)), P T is already exponentially decaying by T T L meaning that a big majority of avalanches are not affected by this finite-size dependent inertial effect. On the other hand, under-damped systems allow for larger durations T > T L , meaning that a growing number of avalanches are long-lived enough to sustain their activ- Increasing the system size does not diminish the inertia fingerprint in the avalanche size distribution. On the contrary, the characteristic "bump" becomes further evident as more massive avalanches are allowed by a bigger simulation box.
ity by using the energy stored in the elastic waves they emitted.
Although the wave speed does not depend on Γ, another time scale relevant for the efficiency of this feedback mechanism is the one governing the damping of the elastic waves emitted by one event, T r . This time scale will increase as damping decreases, and also be affected by the density of active sites that may scatter the wave. The efficiency of the feedback effect will become stronger and create longer lasting avalanches as T r increases compared to T L and the elastic waves are able to cross the system repeatedly. This is illustrated by the drift of the peak in P T towards larger values of T as Γ decreases.
In such a context, it is logic to suspect that the inertial effect observed may be a simple consequence of a finite system size. It is appealing to imagine that the characteristic bimodal distribution of P S may vanish when increasing the system size and that the scale-free scaling would persist longer and longer. Interestingly, this does not happen. On the contrary, at a fixed low damping and increasing the system size we observe that the effect of inertia gets more and more marked.
In Fig. 5 we show the avalanche size distributions at Γ −1 = 1000 and different system sizes. As we increase L an occasional persistence of the power-law behavior is not evidenced. Instead, a depletion is caused in between the scale-free regime and the inertial peak and the peak becomes sharper. Bigger systems, allow for bigger avalanches to develop. The larger the avalanche, the bigger its mass and increased its chances to be pulled out from a scale-free size statistics by inertia.
Finally, let us emphasize that the particular role played by the system size in this discussion is intrinsically related to the use of a quasistatic protocol. At most one avalanche is taking place in the system at any given in- stant. On the other hand, for a system driven at a finite but small strain rate, the size of the system should be replaced by the strain rate dependent distance between avalanches.

B. Local distances to threshold
A key quantity that characterizes the specific avalanche behavior in plastic solids is the distribution P x of local distances from instability x i = σ yi − σ i , with i an index running over the N blocks that compose the system. At each loading stage in the stress time series, the minimum x value in the lattice denoted by x min determines the global instability threshold for the next event. The extent of an avalanche will be controlled by the full distribution P x and its evolution during the stress drop.
It has been shown, following stability arguments [32], that the functional form of P x near the yielding transition should be a power-law P x ∼ x θ with θ > 0 for small x, a pseudo-gap. Moreover, assuming the independence of the x i , a simple scaling argument links the exponent θ with the exponents τ and d f controlling the distribution of avalanche sizes, This prediction was found to hold in numerical simulations of overdamped systems, both in quasi-static protocols [10,32] and approaching the limit of vanishing strain-rate [12] from finite values, giving consistent sets of exponents {θ, τ, d f } in two and three dimensions. We now report the influence of the damping on P x by analyzing configurations right after a stress drop has occurred. Figure 6 shows P x for systems with different damping Γ. Our estimation of a power law P x ∼ x θ in the over-damped case (Γ −1 = 0.2) agrees within error bars with previous estimates of the exponent θ in 2D elasto-plastic models [10,12]. The estimate θ 2D ≈ 0.52 [12] is displayed as a dash line in the inset of Fig. 6, which also illustrates the size dependence of P x affecting the cutoff at small values of x.
When lowering the dissipation rate Γ a systematic change in the behavior of P x is observed, with the development of a steeper gap in P x at small values of x. Basically, as inertial avalanches start to dominate, the probability of surviving (not yielding) with a small value of x (i.e., very close to the local threshold) decreases. Small barriers will be overcome during a massive avalanche. An apparent tendency to preserve a behavior of the form P x ∼ x θ at low x is observed, with an exponent θ that would increase as damping lowers. This seems to contrast with the behavior of P S and P T that show a clear characteristic peak appearing as inertia becomes relevant (making possible a distinction between two kinds of avalanches). However, the situation is a bit more complex. In fact, the steeper growth at the smallest values of x in each curve is indeed a trace of the presence of two kind of events. The fact that P x considers the full set of local values of a configuration in contrast to S or T that only provide a global avalanche property, just renders more difficult the possibility of visualizing the avalanche heterogeneity at this point.
If one insists on thinking the largest power-law growth for each curve Fig.6, as a damping-modified marginal stability scenario where only the exponent θ changes, depending on Γ, we can find no simple explanation for that exponent θ(Γ). In fact, due to the emergence of the upper peak in P S , the scaling relation linking τ , d f and θ derived in [10] for the overdamped case, i.e. τ = 2 − θ θ+1 d d f , is no longer expected to be valid in the inertial case. Instead, we show below that the behavior displayed by P x is the combined result of an alternation between events roughly classified in two different kinds, those similar to the well-known overdamped avalanches and those dominated by inertial effects.

C. Identifying and splitting two types of avalanches
In order to better understand the response of P S , P T and P x to changes in the damping coefficient, we analyze x min ≡ min{x i }, the quantity that determines the loading needed to trigger the next stress drop after an avalanche ends. Figure 7 shows the distributions of such x min values during the stationary plastic flow for systems with different damping Γ. In the overdamped limit we find that P x (x min ) seems to initially grow as x θ min with θ ≈ 0.52, the same exponent as for the initial grow of the full distribution P x (x) [33]. This being anticipated by Karmakar et al. [34], was not seen instead in other works [32], were P x (x min ) seems flat at small x min . As inertia becomes important P x (x min ) acquires a bimodal shape; this is, two characteristic local maxima separated by an intermediate local minimum. Denoting by x cross min the position of the minimum in P x (x min ) (a saddle point for Γ = 1), we find a logarithmic dependence of x cross min with Γ −1 . x cross min roughly separates two kinds of events: On the one hand, those that have been presumably massive, affecting many sites, and have left the system far away from instability, accumulating on a peak of "large" x min values. On the other hand, those events that we can consider to be similar to the ones in the overdamped case, more localized avalanches, involving a small number of sites compared to the full system, and leaving back some relatively weak spots without yielding. A priori, we could link the characteristic "large" x min values, with the typically big stress drops observed in Fig. 2-left for the underdamped systems.
For the time being, let us take the proposed separation as an ansatz, and test its validity with a self-consistency criterion. We then use x cross min to discriminate avalanches according to the x min value that they yield. Figure 8 shows in the main plot the P S distribution for L = 40 and the lowest damping simulated Γ −1 = 1000. The same curve as in Fig.3 is shown by open triangles, and has been built with the contribution of all avalanches in the set. Now we discriminate avalanches yielding small values of x min and avalanches yielding large values of x min and plot their distributions with filled light color circles and filled dark color squares respectively. After rescaling these contributions according to their weight in the full PS versus S shown to be composed by two contributions, one coming from the overdamped-like events and other one coming from the "inertial" events. Data correspond to a linear system size L = 40 (2046 blocks) and Γ −1 = 1000. Upper Insets: Activity maps of selected avalanches. Red dots represent points on the grid that were activated at least once during the avalanche. On the left, a typical avalanche corresponding to the overdamped sub-set. On the right, an inertial avalanche. Lower Inset: PS shown for Γ −1 = 1000 and different system sizes L = 20, 40, 80. An alignment for the position of the rightmost peak in PS is obtained when rescaling S by L d f , with a damping dependent fractal dimension d f (Γ −1 = 1000) ≈ 1.75. P S , we see that the discrimination between two kinds of avalanches, provided by x cross min is quite accurate. Other values of Γ (not shown) display the same beautiful splitting of P S contributions, while, of course, the low-x min contribution earns more and more weight as the overdamped limit is approached.
As it happens in studies of avalanches in spin systems [35], by accepting the idea of a separation between different kinds of avalanches, one opens the door to go further and analyze scaling properties of their respective distributions. In the following we focus on some scaling features of inertial avalanches that contrast with the overdamped case.

D. Growing fractal dimension and incipient shear localization
Let us first briefly recall the finite-size scaling discussion of Sec.III A. The lower inset of Fig.8 shows the same data as in Fig. 5, corresponding to Γ −1 = 1000 and L = 20, 40, 80, now rescaled as P S S d f τ vs S/L d f , where we have used d f = 1.75 and τ = 1.5. Notice that the inertial effect seems, in fact, stronger for bigger systems; the relative peak height increases with L. Of course, this does not happen for the equivalent scaling in the overdamped limit (Fig. 3-inset), where the distribution cutoffs identically collapse. This tells us that the basic ingredient for the emergence of two different kinds of collective events is, indeed, the lack of dissipation, and not the finiteness of the box. Now we turn our attention to the scaling chosen in the abscissas axis of Fig.8-inset. Here we have attempted to "align" the positions of the inertial peaks. When doing so, we find a dependence S peak ∼ L d f , where d f is a damping dependent fractal dimension for the avalanches and governs the "new" cutoff of P S . In all cases d f is bigger than the fractal dimension of the overdamped case. An effective fractal dimension bigger than d f suggests that the geometry should be very different between overdamped and inertial avalanches. Indeed, this is clearly manifested in the spatial coverage of a typical avalanche. As insets of Fig. 8 we show two images representing a system sample, where we have depicted activity maps of selected avalanches. Red dots represent points on the grid that were activated at least once during the avalanche. On the left, a typical avalanche corresponding to the overdamped sub-set shows a sparse quasi-1d arrangement of active sites, consistent with a fractal dimension d f ≈ 0.95 as the one obtained from the finite size scaling of overdamped systems [12,27]. On the right, an inertial avalanche belonging to the rightmost peak of P S , shows a much more dense and broad structure, compatible with a fractal dimension closer to the dimension of the system (d = 2). When we inspect the largest avalanches for different damping coefficients, we find that the generic form their pattern in space changes systematically, becoming denser as d f evolves from 0.95 to 1.75 in the studied regime.
Inertia, therefore, is found to modify the geometry and fractal dimension of the avalanches at the yielding point, creating much denser events. The geometry of these system spanning and broad avalanches are suggestive of incipient shear bands. In fact, very recent works both using finite elements methods [29] and molecular dynamics [28], have associated the absence of dissipation with a non-monotonic flow-curve in the rheology of the system, and therefore mechanical instability and the emergence of strain localization (see also [22]). Although we are working with a quasi-static protocol, the self-sustained effect of inertial waves generating the inertial avalanches may be the same effect that generates the shear bands in finite strain rate driving protocols.

E. Px contributions and stability-geometry scaling relation
We now turn to the analysis of the different contributions to P x . Again, using the splitting criterion between overdamped-like avalanches, that leave a x min < x cross min , and inertial avalanches, that yield a x min > x cross min , we plot the full-set distributions of distances to yielding P x together with two sub-set distributions. Figure 9 show such a superposition of curves for our as a function of Γ −1 .
lowest damping system (Γ −1 = 1000). We first notice that the distribution P x for the sub-set of configurations after an inertial avalanche has a much sharper gap than the global distribution. This is consistent, of course, with the fact that these events typically have a larger x min value than the average x min . In addition, a remarkable finding is that the P x of inertial events display also a power-law growth at small x. In other words, despite the massive events that push all local sites to have a relatively large value of x after the avalanche, these values are still arranged in a marginally stable fashion, as P x ∼ x θ , but now with a steeper exponent θ > θ overdamped . In fact, assuming that the very general finite-size property x min ∼ L − d θ +1 [32,36] also holds for this sub-set, and asking x min ∆σ during stationary plastic flow, we expect where d f is the exponent ruling the finite-size dependency of the peak position in P S (see Fig.8 inset), or equivalently, the scaling of the sub-set of inertial avalanches. The above relation between θ and d f holds well for all low-damping systems studied in this work, as long as ∆σ is controlled by the inertial peak. This can be seen by fitting θ in the inertial sub-set of the P x distribution. For example, for Γ −1 = 1000 we have θ 6.91 and d f 1.75, for Γ −1 = 10, θ 3.12 and d f 1.5. Another interesting observation is that, even the subset of configurations that are left behind by an overdamped-like avalanche show at the smaller values of x a growth that tends to be compatible with the same exponent θ (see dotted line, parallel to the dashed fit in Fig.8). We interpret this feature as a fingerprint of inertial avalanches, with a depletion (θ > θ) in the amount of sites close to yielding even after one or several smaller avalanches has taken place. This fingerprint of inertial avalanches could also explain the dependence on damping of the exponent τ in the "overdamped-like" subsets of avalanche size distributions. Even when short-duration avalanches do not directly feel the effect of inertia, they have to deploy correlations on a particularly heterogeneous landscape left behind by inertial avalanches. As a result, the separation in two classes appears as a convenient classification, but does not fully account for the complexity introduced by inertial effects.

IV. SUMMARY
We have analyzed the noise statistics of stress signals produced by an amorphous solid under quasi-static deformation, through numerical simulations of a realistic continuum model treated with classical finite-element techniques. In particular we have focused our analysis on the dependence of such statistics with the ability of the system to dissipate energy, spanning a wide range of damping values. We validate our model by comparing our results for the distributions of different observables in the overdamped limit with previously reported numerical results in a variety of different models and techniques.
Both the distributions of avalanche sizes P S and durations P T display a growing peak at high values of S and T , respectively, when we lower the damping coefficient. We have associated these characteristic peaks with a "new" kind of avalanches, peculiar to inertial systems, which are triggered and amplified by elastic waves generated by other avalanches. Such events are, in nature, more related to effective thermal heating than to the deployment of large spatial correlations, in line with [28]. These system spanning inertial avalanches are also characterized by a geometry that is reminiscent of shear bands of strain localization, although our protocol is quasi-static.
The distribution of minimal distances to yielding P x (x min ) allowed us to formulate an ad hoc criterion to discriminate between these inertial avalanches and overdamped like events that are much more localized and not affected by the elastic wave propagation. Using this criterion we have found an explanation for the behavior at different damping of the full-set P x distribution, that is considered as the steady state property that controls the stability of the system and indirectly all the properties of its noise signal. In particular, following very general arguments, we propose a scaling relation between the distribution P x of inertial avalanches and their fractal dimension (Eq. 3).

About the upper size cutoff
As mentioned in the introduction, the irruption of inertial effects in the -otherwise overdamped-driven dynamics has been found to break down the universal avalanche statistics of several systems. It is worth stressing, though, an important difference between the sand-pile problem with local yield thresholds and invariably positive load redistributions (depinning models included), and the elasto-plastic models that describe plastic flow in solids: Systems with a positive load redistribution produce avalanches that are compact objects in space (d f ≥ d). The system size L controls the cutoff for the avalanche sizes. An overdamped system will explore this upper boundary displaying some system-size spanning avalanches, but there is no possibility to observe bigger avalanches than S c ∼ L d f . When the overdamped condition is released in such systems, the avalanche size distribution is modified [17]. The original power-law of the overdamped case deforms into a shorter scale-free regime followed by a kink or peak that depends on the value of the damping. However, for all damping choices, the largest avalanche accessible for a given system size remains the same.
In the deformation of amorphous solids, instead, we have avalanches with d f < d due to the heterogeneous (Eshelby) redistribution of stresses. The cutoff of the avalanche distribution in overdamped systems is far from being a massive avalanche involving all sites of the system. In fact, both for 2D and 3D systems, d f is found to be close to one [12,27]; meaning that even a systemspanning avalanche leaves most of the sites untouched. Therefore, when inertia comes into play, the system still has room to make avalanches grow further. This can be seen in Fig. 3. As we lower the damping starting in the overdamped limit, first a plateau and then a peak develops in P S ; all this happening to the right of the overdamped system-size cutoff. May this explain why, for moderate damping, inertial effects are not strongly evident in yield stress systems and avalanche distributions remain quite similar to the one of the overdamped case, with occasionally the added features of a small plateau at large sizes and a weak change in the τ exponent. A picture that makes it appealing to conclude about a still universal scenario, slightly modified by inertia [26,27].
However, going deeper in the inertial regime while understanding the system size role, we realize that indeed critical behavior breaks down while a damping-dependent typical event size emerges in the form of a very clear peak in P S . Even though a scale-free power-law regime remains observable and is only weakly affected by inertia, the inertial peak increasingly dominates the statistics of events. Furthermore, since systems that dissipate energy at different rates show dissimilar avalanche distributions, even when characterized by some scaling exponents, we find it inaccurate to talk about universality.

Connections to other non universal statistics
We have seen that the characteristic avalanche distribution of underdamped systems in our prescription is a superposition of two populations. Inertial and overdamped-like events interleave in time producing a unique statistics that we cannot discriminate beforehand. Both kinds of events are results of the same dynamical rules and boundary conditions. It simply happens that from the competition of different time scales present in the dynamics, both kinds emerge. Even more, sometimes we cannot tell if one event is of one kind or the other.
Such a situation is not unique to the introduction of inertial effects. It can also be observed in cases where a second time scale is introduced by viscoelasticity [37] or retardation [22]. Indeed, a quasi-periodic oscillation of the stress field was argued in [37] to be a possible explanation for deviations from a pure power law (Gutenberg-Richter (GR) like) in the distribution of earthquake magnitudes. In fact, some years ago, a discussion arose in the seismology community contrasting opposite models of earthquake statistics: On one hand the famous GR power-law decay kind of distribution; on the other hand, the "characteristic earthquake" hypothesis predicting a time recurrence of typically big earthquakes of a characteristic magnitude on each individual fault (see [38,39] and references therein). In spite of a marked predilection for a pure GR picture, the discussion remains somewhat open [39][40][41], and specific models with such characteristics are developed for single faults [42]. We do not pretend to accredit here our simple model with relevance to the phenomenology of earthquakes, but to highlight the ubiquity of the "characteristic event" feature. In line with previous works [23][24][25], our model shows how a single set of dynamical rules can lead to the emergence of distinct kind of events, ones with no characteristic size and a Gutenberg-Richter distribution, and others with a magnitude that fluctuates around a typical value in a peaked fashion. For example, the role played by inertia in our model can be compared to the idea of "dynamic triggering" of earthquakes [43], which has attracted considerable attention in the last ten years. In the present case we see that dynamical amplification through sound waves, rather than triggering, appears as a dominant mechanism. Our large avalanches could however be understood as consisting of an initial event dynamically triggering a series of aftershocks at remote distances through the wave propagation, the total stress drop magnitude being the outcome of the whole series.
Finally, we note that laboratory systems such as metallic glasses also display statistics that deviate from a pure power law with exponential cutoff, as can be seen from the inspection of cumulative distributions shown for example in reference [7]. Making similar studies in granular suspensions, in which inertial effects can be controlled by using solvents of various viscosities, would be of great interest. Interestingly, recent experiments on granular layers sheared between elastic plates [44], intended as a model for an isolated strike-slip earthquake fault, also indicate a separation between two classes of events, with large, system spanning events emerging in the tail of a broad continuum power-law spectrum.

Concluding remarks
Overall, inertia and amplification through sound waves appears as a possible mechanism for enriching the statistics of intermittent phenomena in deformed solids, with deviations from universal power law statistics and large system spanning events that should be connected with the possibility of shear band formation in these systems. We hope that this work may stimulate experimental studies of intermittent behavior in systems that may display strong inertial effects and/or strain localization, and we also consider generalizing this study to finite strain rates in the future.