Evolution of the average avalanche shape with the universality class

A multitude of systems ranging from the Barkhausen effect in ferromagnetic materials to plastic deformation and earthquakes respond to slow external driving by exhibiting intermittent, scale-free avalanche dynamics or crackling noise. The avalanches are power-law distributed in size, and have a typical average shape: these are the two most important signatures of avalanching systems. Here we show how the average avalanche shape evolves with the universality class of the avalanche dynamics by employing a combination of scaling theory, extensive numerical simulations and data from crack propagation experiments. It follows a simple scaling form parameterized by two numbers, the scaling exponent relating the average avalanche size to its duration and a parameter characterizing the temporal asymmetry of the avalanches. The latter reflects a broken time-reversal symmetry in the avalanche dynamics, emerging from the local nature of the interaction kernel mediating the avalanche dynamics.

T he theoretical interpretation of crackling noise 1 , observed in numerous systems, including the Barkhausen effect in ferromagnetic materials 2,3 , plastic deformation 4-6 , structural transitions 7 and fracture 8,9 of solids and earthquakes 10 , has found a formulation in terms of non-equilibrium phase transitions 11 . These transitions separate quiescent and active phases of the system, and naturally give rise to critical scaling 12 . In the vicinity of such a phase transition, the time evolution of the activity signal V(t) or the order parameter of the transition (for example, the interface velocity for a depinning transition) exhibits scale-free bursts or avalanches. Statistical analysis of such fluctuations, together with renormalization group calculations 13 , suggest that in general systems with avalanche dynamics can be classified into universality classes characterized by the values of the critical exponents, depending on, for example, the spatial dimension and the interaction range of the system. The average temporal shape of bursts in a crackling noise signal is a fundamental signature of avalanches, and has been estimated for systems as diverse as plastically deforming crystals 14 , earthquakes 15 and Barkhausen noise [16][17][18] . For the latter, the symmetric average avalanche shape observed in ferromagnetic films of intermediate thickness where the long-range dipolar interactions render the avalanche dynamics mean field-like has been explained within the ABBM model 19 , and shown to be given by an inverted parabola 16,[20][21][22] . In thick enough samples eddy currents induce an effective mass for the propagating domain walls, visible as an asymmetry in the (mean-field-like) average shape of the Barkhausen pulses 17 . In general, clear-cut shape determinations should give strong indications of the underlying physics, such as the kind and range of interactions governing the avalanche dynamics.
Here we present a general scaling form for the average avalanche shapes for non-mean-field systems. It is verified within a large-scale numerical study of avalanches at the depinning transition of driven elastic interfaces in random media. We vary systematically the range of the elastic interaction kernel, and thus the universality class of the avalanche dynamics 23 , showing how the avalanche shape depends on the universality class. The average shape is to a high precision given by a function parameterized by the scaling exponent g characterizing the scaling of the average avalanche size as a function of the avalanche duration, and a parameter a describing the temporal asymmetry of the average avalanche shapes. We find that an inherent asymmetry in the average avalanche shapes is present in systems where the interaction kernel is not fully non-local, reflecting the underlying broken time-reversal symmetry of the avalanche dynamics. Finally, we compare these results with experiments of planar crack front propagation, finding good agreement with the predictions of the scaling theory and the relevant depinning model.

Results
Average avalanche shape scaling function. To obtain a general scaling form for the average shape of the bursts in V(t) corresponding to avalanches of a given duration T, hV(t | T) i, we start from the well-known result that the average avalanche size hs(T) i follows in the scaling regime Thus, the amplitude of the average avalanche shape of avalanches of duration T has with T the relation hV(t | T) ipT g À 1 , and is in general given by the form where f(x) is a scaling function characterizing the average temporal avalanche shape. We make the Ansatz that the average early-time growth of an avalanche is given by a power law of time, that is, hV(t | T) ipt d for t=T b ( 1. Given that equation (2) implies that hV(t | T) ipT g À 1 for any fixed t/T, it follows in particular that hV(bT | T) ipT g À 1 . Therefore, in order for the postulated power law early-time growth to be compatible with this, one needs to have hV(bT | T) ip (bT) d pT g À 1 , that is, d ¼ g À 1. Thus, the average early-time growth scales with t as To arrive at an expression for hV(t | T) i satisfying equations (2) and (3), we multiply equation (3) by (1 À t/T) g À 1 , to describe the (symmetric) deceleration at the end of an avalanche, and obtain In mean-field systems, hV(t | T) i is known to be given in the scaling regime by the inverted parabola, predicted within the ABBM model in the limit of vanishing drive rate and demagnetizing factor, hV(t | T) ipt(1 À t/T)T(t/T)(1 À t/T) (ref. 16). This is in agreement with equation (4), given the mean-field value g ¼ 2. Notice also that random walk bridges obey similar scaling with g ¼ 3/2 (refs 24,25).
Although avalanches in mean-field systems have a symmetric average shape 16 , there is no fundamental reason why this should be true in general. Indeed, it is easy to see from the space-time activity plots of avalanches from, for example, the local quenched Edwards-Wilkinson (qEW) equation 26 used in this study (see Fig. 1a) that the internal dynamics of the avalanches appears to violate time-reversal symmetry: the branching exhibited by the activity pattern at all scales clearly defines a direction of time. This means that the system is not Markovian: a possible signature of that is a temporal asymmetry in the average avalanche shape hV(t | T) i (ref. 24). For an increasing range of the elastic interactions, the time-irreversible nature of the activity pattern becomes less evident (Fig. 1b) and naturally vanishes for the mean-field infinite range model ( Fig. 1c; see below for the definitions of the various interface depinning models): the activity pattern becomes a cloud of points without any internal timeirreversible structure.
Therefore, to account for the time-irreversible avalanche dynamics, we need to allow for a small temporal asymmetry in the average avalanche shape. Thus, we multiply equation (4) by 1 À a(t/T À 1/2), that is, a first order correction term allowing for an asymmetry quantified by a, and obtain If a ¼ 0, equation (5) reduces to equation (4) and corresponds to a symmetrical avalanche shape. With aa0, equation (5) describes a temporally asymmetric avalanche shape, with a positive or negative skewness for a40 and ao0, respectively (see also Supplementary Note 1).
Numerical simulations of interface depinning models. To check equation (5), we perform extensive simulations of a discretized model of a 1d elastic string or interface in a 2d random medium 23 , represented by a set of integer heights h i (t), i ¼ 1 y L, with L the system size. The lateral coordinates x i of the interface are given by x i ¼ i. The total force acting on the interface element i is where the first term on the RHS represents the elastic interactions characterized by the exponent a, Z is uncorrelated quenched disorder and F ext is the external driving force. Notice that in the limit a-N, the elastic interaction term becomes completely local, G 0 (h i þ 1 þ h i À 1 À 2h i )G 0 r 2 h i : thus, equation (6) reduces to the qEW equation. In the opposite limit (a-0) the system loses its spatial structure and we describe it by the mean-field infinite-range model, by replacing the elastic interactions in equation (6) by G 0 (6) reduces to the long-range elastic string expected to describe, for example, planar crack fronts propagating along disordered weak planes between solid blocks 27,28 , contact lines of liquids spreading on solid surfaces 29 and low-angle grain boundaries in plastically deforming crystals 30 . Furthermore, aZ3 belongs to the qEW class, whereas a ¼ 1 is expected to follow mean-field dynamics. The crackling noise signal is given by VðtÞ with y the Heaviside step function. For additional details, see Methods.
As expected, hs(T) i scales with T according to equation (1) (see Supplementary Fig. S1 and Supplementary Note 2). We find three different values of g, that is, g ¼ 2.0 ± 0.01 for ar1, g ¼ 1.79±0.01 for a ¼ 2 and g ¼ 1.56±0.01 for aZ3. These values are in agreement with earlier results, either directly or via scaling relations 12,16,27,31 . Fitting equation (5) to the hV(t | T) i data for various a and different ranges of T in the scaling regime reproduces well these g-exponents (Fig. 2a,b; see also Supplementary Figs S2-S5 and Supplementary Notes 3 and 4). The avalanche shapes exhibit an asymmetry, with the a-parameter evolving continuously from À 0.01 ± 0.02 (for the infinite range model) to 0.081±0.015 (for the qEW equation), as the interaction range is reduced ( Fig. 2c; see also Supplementary  . The corresponding skewness (computed by interpreting hV(t | T) i as a probability density 17 ) of the avalanches exhibits a similar evolution with a ( Supplementary  Fig. S7, Supplementary Note 6). Thus, avalanches whose dynamics is governed by interaction kernels that are not fully non-local are temporally asymmetric, as illustrated by the timeirreversible nature of the corresponding space-time activity patterns (Fig. 1).
Planar crack front propagation experiments. Finally, we consider data from planar crack front propagation experiments 9,32 , as an example of an experimental system with non-mean-field avalanche dynamics, see Methods for details. The scaling of the average size of the avalanches of crack front propagation as a function of their duration is shown in Fig. 3b. In the scaling regime, these are characterized by g ¼ 1.67±0.15, in agreement with the 1d non-local elasticity depinning model with a ¼ 2 (refs 27,28), see also Supplementary Fig. S8. The average avalanche shape is shown in Fig. 3a. Owing to the nonnegligible statistical fluctuations present in the data, it is not possible to detect the small asymmetry predicted by the crack line model: thus, we choose to fit the leading-order behaviour, equation (4), to the data (see also Supplementary Fig. S9). This leads to g ¼ 1.74 ± 0.08, in agreement with the g-value obtained from the fit to the shape obtained from the crack line model 27,28 .
Notice also that the experimental shape clearly differs from both the mean-field inverted parabola and the shape expected for the local qEW equation.

Discussion
We have shown how the average avalanche shape of systems exhibiting crackling noise depends on the universality class of the avalanche dynamics. It is a fundamental fingerprint of an avalanching system and extrapolates when tuning elastic interactions between an inverted parabola for mean-field systems and a shape close to a semicircle for the 1d short-range interface. The broken time-reversal symmetry in the avalanche dynamics emerging from the spatially localized interactions is manifested as a temporal asymmetry in the avalanche shape evolving with the interaction range (see also Supplementary Discussion). Thus, such asymmetries should be looked for in experimental data in systems where the interactions mediating the avalanche dynamics are not fully non-local. These include, for example, domain wall dynamics in magnetic thin films 33 and fluid invasion into disordered media 34,35 .

Methods
Numerical simulations of interface depinning models. We simulate the interface depinning model, equation (6) with periodic boundary conditions. The parallel dynamics of the interface is defined in discrete time t by setting the local velocity where y is the Heaviside step function. The interface is driven with a quasistatic constant velocity drive, where avalanches are triggered by increasing F ext just enough to make exactly one interface element unstable (that is, F i 40 for some i) whenever the previous avalanche has ended. Thus, avalanches can be defined unambiguously without thresholding 28 . During an avalanche, F ext is decreased at a rate proportional to the instantaneous avalanche velocity, k is a parameter analogous to the demagnetizing factor for ferromagnetic domain walls 36 or to the elastic stiffness of the specimen-machine system in a mechanical loading experiment 37 and controls the cut-off of the avalanche distribution, with the cut-off size s 0 obeying s 0 $ k À 1=sk . The crackling noise signal of interest is given by VðtÞ ¼ P i v i ðtÞ. To compute the average avalanche shapes, we collect a large ensemble of avalanches from various duration ranges. For T in the scaling regime, the average shapes corresponding to the various duration ranges fall onto a single curve after normalizing with the maximum amplitude, hV(t | T) i max . The simulations are performed in large system sizes (L ¼ 8,192 ¼ 2 13 for aZ3, L ¼ 32,768 ¼ 2 15 for a ¼ 2, L ¼ 131,078 ¼ 2 17 for a ¼ 1 and L ¼ 8,388,608 ¼ 2 23 for the infinite range model), and by using sufficiently small k-values such that the avalanche cut-off s 0 $ k À 1=sk is large. Examples of the crackling noise signals VðtÞ ¼ P i v i ðtÞ obtained from the model for different interaction ranges are shown in Fig. 4. Owing to the quasistatic driving mechanism, the avalanches can be unambiguously defined without thresholding. Here, for visual clarity a fixed quiet time t q ¼ 100 has been added between each avalanche. All three signals are of the same length.
crack front propagating along the heterogeneous weak plane of a transparent poly (methyl methacrylate) block, made of two sintered rough Plexiglas plates (of dimensions (27, 14 and 1) cm for the top and (30, 12 and 0.4) cm for the bottom plate, respectively) is studied 8,9,32 . We imposed a constant normal displacement d to the bottom plate while the upper one is fixed, resulting in a quasi-mode I creep growth of the crack, see Fig. 5a. The interfacial fracture front was observed in a small central region (to avoid boundary effects) corresponding to around 4.48 mm and 6.72 mm in the direction of propagation and along the front, respectively. The propagation of the crack front is monitored optically by using a camera with a high frame rate compared with the average crack front velocity, recording typically 20,000 images of 2,000 Â 3,000 pixels with a pixel size of r ¼ 2.24 mm at a rate ranging from 1 to 60 fps. From these images, we define the local velocity v along the crack front by measuring the waiting time (wt) the front has spent at a given position (x,y) and setting v(x,y) ¼ r/wt(x,y) (ref. 9), see Fig. 5b. We consider here the 'global' velocity of the crack front V(t), computed (that is, spatially averaged from the local velocities) at a length scale of dx ¼ 200 mm, larger than the correlation length of 100 mm of the local velocities 32 , leading to more than 30 different and independent crackling noise signals V(t) for each creep test (since the lateral size of the imaging area divided by dx is (6720 mm)/(200 mm) ¼ 33.6). Five different creep experiments with low average crack front velocities V ¼ hV(t) i t (V ranging from 0.017 mm s À 1 to 1.13 mm s À 1 ) are considered, resulting in more than 150 crackling noise signals containing in total more than 4,000 bursts or avalanches. An example of a signal V(t) is shown in Fig. 5c.
During a creep test, the average crack front velocity decreases slowly. However, the time intervals studied here are short enough to consider the average crack front velocity constant within such time windows. Avalanches are defined as a continuous occurrence of V(t) above its average value V ¼ hV(t) i t . We have verified that by choosing another threshold values in a reasonable interval around V does not influence the extracted results. The average value is subtracted from V(t) to define the avalanche size (that is, the integral of V(t) À V) and shape. As the different experiments have a different average velocity, we consider the rescaled avalanche durations T 0 ¼ T/T V , where the time scale T V is given by T V ¼ r/V. Also, very short and very long avalanches might affect the scaling behaviour as well as the avalanche shape, owing to a poor resolution and a lack of statistics, respectively. Therefore, we consider only avalanches that last longer than T ¼ 25 dt, (the temporal resolution dt varying between 1/60 s up to 1 s for the various experiments performed). In the inset of Supplementary Fig. S8, we show the probability density functions P(T 0 ) of the normalized durations for all avalanches extracted and for such a subset, that is, P(T 0 |T425). We can observe an exponential cut-off for large duration T 0 42.5, which provides the upper limit for the scaling range we will consider, leading finally to the study of slightly more than 1,500 avalanches. A power-law fit to the avalanche size s as a function of the normalized durations T 0 for those avalanches yields an exponent g ¼ 1.67 ± 0.15 (as shown in Fig. 3b, and in the Supplementary Fig. S8), in good agreement within error bars with the 1d depinning model with a ¼ 2 (ref. 27). We finally normalize the amplitudes of those avalanches by T 0g À 1 and compute their average shape hV(t | T 0 ) i. The resulting average avalanche shape normalized by the maximum amplitude, hV(t | T 0 ) i/hV(t | T 0 ) i max , is shown in Fig. 3a.
Finally, we discuss the effect of the limited statistics available from the experiments. In order to estimate the amount of statistics required to observe the predicted asymmetry, we estimate the statistical error bars d N (V/V max ) of the average avalanche shapes averaged over N avalanches as d N ðV=V max Þ ¼ is the fluctuation magnitude of an individual avalanche. As the typical difference between the normalized avalanche shape of the depinning model with a ¼ 2 and the corresponding symmetrical shape is roughly 0.005 or less (see Fig. 2c), the error bars have to be smaller than that to be able to distinguish the asymmetry from the data. Setting d N (V/V max ) ¼ 0.005 corresponds to N ¼ 40,000, which should be interpreted as a lower limit. Thus, we estimate the order of magnitude of the number of avalanches N required to make conclusions regarding the asymmetry to be of the order of N ¼ 10 5 , that is, two orders of magnitude more than the number of avalanches we have at our disposal from the experiments. Indeed, this can also be seen from the data: in Supplementary Fig. S9 we show again the experimental average avalanche shape, along with a fit of equation (4). The inset displays their difference, showing clearly that the very small correction to the symmetric shape predicted by the crack line model is hidden by statistical noise.