Finite Size Scaling analysis of the avalanches in the 3d Gaussian Random Field Ising Model with metastable dynamics

A numerical study is presented of the 3d Gaussian Random Field Ising Model at T=0 driven by an external field. Standard synchronous relaxation dynamics is employed to obtain the magnetization versus field hysteresis loops. The focus is on the analysis of the number and size distribution of the magnetization avalanches. They are classified as being non-spanning, 1d-spanning, 2d-spanning or 3d-spanning depending on whether or not they span the whole lattice in the different space directions. Moreover, finite-size scaling analysis enables identification of two different types of non-spanning avalanches (critical and supercritical) and two different types of 3d-spanning avalanches (critical and subcritical), whose numbers increase with L as a power-law with different exponents. We conclude by giving a scenario for the avalanches behaviour in the thermodynamic limit.


I. INTRODUCTION
Systems with first-order phase transitions exhibit a discontinuous change of their properties when driven through the transition point. Sometimes, due to the existence of energy barriers larger than thermal fluctuations, such systems evolve following a path of metastable states and exhibit hysteresis. Metastable phenomena develop more often in the case of systems at low temperature and with quenched disorder. In many cases the first-order phase transition occurs, instead of at a certain transition point, in a broad range of the driving parameter and the discontinuity is split into a sequence of jumps or avalanches between metastable states. Moreover, under certain conditions such avalanches do not show any characteristic spatial or time scale: the distribution of their size and duration becomes a power-law. This framework, which has sometimes been called fluctuationless first-order phase transitions [1,2], is one of the basic mechanisms responsible for power-laws in nature [3]. Experimental examples have been found in a broad set of physical systems: magnetic transitions [4], adsorption [5], superconductivity [6], martensitic transformations [7], etc.
A paradigmatic model for such fluctuationless firstorder phase transitions in disordered system is the Gaussian Random Field Ising Model (GRFIM) at T = 0 driven by an external field H. The amount of quenched disorder is controlled by the standard deviation σ of the Gaussian distribution of independent random fields acting on each spin. Metastable evolution is obtained with appropriate local relaxation dynamics which assumes a separation of time scales between the driving field rate dH/dt and the avalanche duration. The response of the * Electronic address: jperez@ecm.ub.es † Electronic address: eduard@ecm.ub.es system to the driving field can be followed by measuring the total magnetization m(H). The response exhibits the above-mentioned metastable phenomena: hysteresis and avalanches.
Since the model was introduced some years ago [8,9], different studies (numerical and analytical) have been carried out in order to characterize the hysteresis loops m(H) and the magnetization avalanches [10,11,12,13,14,15]. Two of the most well-studied properties are the number of avalanches N (σ) and the distribution D(s; σ) of avalanche sizes s along half a hysteresis loop. For large amounts of disorder (σ > σ c ) the loops look smooth and continuous. They consist of a sequence of a large number of tiny avalanches whose size distribution D(s; σ > σ c ) decays exponentially with s. On the other hand, for small amounts of disorder (σ < σ c ), besides a certain number of small avalanches, one or several large avalanches produce a discontinuity ∆m in the hysteresis loop. For an intermediate critical value σ c the distribution D(s, σ c ) of avalanche sizes s can be approximated by a power-law: D(s; σ c ) ∼ s −τ .
Many of the properties of the GRFIM have been understood by assuming the existence of a T = 0 critical point (σ c , H c ) on the metastable phase diagram. The more recent estimation [13] renders: σ c = 2.16 ± 0.03 and H c = 1.435 ± 0.004. Although partial agreement on the values of the critical exponents has been reached, other features are still controversial.
One of the fundamental problems is the definition of the order parameter. From a thermodynamic point of view the discontinuity of the hysteresis loop ∆m seems to be an appropriate order parameter if ∆m > 0 for σ < σ c and ∆m = 0 for σ > σ c . Nevertheless, in the T = 0 numerical simulations, due to the finite size of the system and for a given realization of disorder, all the magnetization changes are discontinuous. Note that this does not occur for standard thermal numerical simulations in which, due to thermal averaging, magnetization is continuous for finite systems. Only finite-size scaling analysis will reveal which are the "large avalanches" and whether or not avalanches become vanishingly small in the thermodynamic limit. It is thus very important to study the properties of the "spanning" avalanches. These are avalanches that, for a finite system with periodic boundary conditions, cross the system from one side to another. In particular it would be interesting to measure the number N s (σ) of spanning avalanches and their size distribution D s (s; σ).
A second unsolved question, related to the previous one, is the spatial structure of the avalanches. It has been suggested that they are not compact [10,16]. A fractal dimension (d f = 1/0.34 < 3) has been estimated from the avalanche size distribution [11]. It would be interesting to understand how such a fractal behaviour may, in the thermodynamic limit, represent a magnetization discontinuity.
A third problem is the definition of the scaling variables in order to characterize the critical properties close to the critical point (σ c , H c ). When focusing on the study of avalanche properties, it should be pointed out that the scaling analysis is performed by using quantities (N (σ) and D(s, σ)) measured recording all the avalanches along half a hysteresis loop. The measurement of non-integrated distributions, i.e. around a certain value of H, will require large amounts computing effort in order to reach good statistics for large enough systems. Therefore, the dependence on the field H is integrated out and the distance to the critical point σ c is measured by a single scaling variable u(σ). Although in pioneering papers [8,9] the most usual scaling variable u 1 = (σ − σ c ) /σ c was used in order to scale the avalanche size distribution, forthcoming studies [10,11,13] changed the definition to u 3 = (σ − σ c ) /σ. Apparently both definitions are equivalently close to the critical point, but it can be checked that the "phenomenological" scaling of the distributions D(s; σ) using u 3 (with u 3 > 0.04) as suggested in the inset of Fig. 1 in Ref. 10, is not possible when using u 1 .
Finite-size scaling analysis has been carried out [13,16] for the number of spanning avalanches N s (σ; L). Nevertheless, such finite-size scaling has not been presented either for the avalanche size distributions D(s; σ, L) or for the number of non-spanning avalanches N ns (σ; L). Most of the studies [10,13] have proposed collapses by neglecting the fact that simulated systems are finite. There is an exception [12,17] for which the scaling of the avalanche distributions with L has been studied. In this case, nevertheless, the dependence on the distance to the critical point has been neglected and, consequently, parameterdependent exponents have been obtained. In our opinion, scaling of the avalanche distribution must be studied on a two dimensional plane, including a scaling variable which accounts for the finite-size L and another which accounts for the distance to the critical point.
Previous studies have provided simulations of very large system sizes (up to L = 1000) [14]. This has been advantageous for the study of self-averaging quantities.
Nevertheless, the properties of the spanning avalanches are non-self averaging. This is because, as will be shown, the number of spanning avalanches per loop does not grow as L 3 . This means that, in order to obtain better accuracy, it is more important to perform averages over different disorder configurations (which will be indicated by · ) than to simulate very large system sizes.
In this paper we present intensive numerical studies of the metastable 3d-GRFIM and focus on analysis of the spanning avalanches. In section II the model, the definition of a spanning avalanche and the details of the numerical simulations are presented. In section III raw numerical results are given. In section IV some of the Renormalization Group (RG) ideas will be reviewed, which will be taken into account for the analysis of the critical point. A finite-size scaling analysis of the avalanche numbers is presented in section V. The same analysis for size distributions and their k-moments are presented in sections VI and VII respectively. Section VIII presents a discussion on the behaviour of magnetization. The discussion of the results in relation with previous works is presented in section IX. Finally in section X a full summary and conclusions is given.

II. MODEL
The 3d-GRFIM is defined on a cubic lattice of size L × L × L. On each lattice site (i = 1, . . . , L 3 ) there is a spin variable S i taking values ±1. The Hamiltonian is: where the first sum extends over all nearest-neighbour (n.n) pairs, H is the external applied field and h i are quenched random fields, which are independent and are distributed according to a Gaussian probability density where the standard deviation σ, is the parameter that controls the amount of disorder in the system. Note that h i = 0 and h 2 i = σ 2 . The system is driven at T = 0 by the external field H. For H = +∞ the state of the system which minimizes H is the state with maximum magnetization m = When the external field H is decreased, the system evolves following local relaxation dynamics. The spins flip according to the sign of the local field: where the sum extends over the 6 nearest-neighbouring spins of s i . Avalanches occur when a spin flip changes the sign of the local field of some of the neighbours. This averaged number notation avalanches N (σ, L) spanning avalanches Ns(σ, L) non-spanning avalanches Nns(σ, L) critical non-spanning avalanches Nnsc(σ, L) non-critical non-spanning avalanches Nns0(σ, L) 1d-spanning avalanches N1(σ, L) 2d-spanning avalanches N2(σ, L) 3d-spanning avalanches N3(σ, L) critical 3d-spanning avalanches N3c(σ, L) subcritical 3d-spanning avalanches N3−(σ, L) normalized size distribution notation avalanches D(s; σ, L) spanning avalanches Ds(s; σ, L) non-spanning avalanches Dns(s; σ, L) critical non-spanning avalanches Dnsc(s; σ, L) non-critical non-spanning avalanches Dns0(s; σ, L) 1d-spanning avalanches D1(s; σ, L) 2d-spanning avalanches D2(s; σ, L) 3d-spanning avalanches D3(s; σ, L) critical 3d-spanning avalanches D3c(s; σ, L) subcritical 3d-spanning avalanches D3−(s; σ, L) may start a sequence of spin flips which occur at a fixed value of the external field H, until a new stable situation is reached. H is then decreased again. This "adiabatic" evolution corresponds to the limit for which avalanches are much faster than the decreasing field rate. Note that, once the local random fields are fixed, the metastable evolution is completely deterministic, no inverse avalanches may occur and the hysteresis loops exhibit the return point memory property [8].
The size of the avalanche s corresponds to the number of spins flipped until a new stable situation is reached. Note that the corresponding magnetization change is ∆m = 2s/L 3 .
For a certain realization of the random fields, corresponding to a given value of σ, we have recorded the sequence of avalanche sizes during half a hysteresis loop, i.e. decreasing H from +∞ to −∞. The two main quantities (see Table I) that are measured after averaging over different realizations of disorder, are the total number of avalanches per loop N (σ, L) and the distribution of avalanche sizes D(s; σ, L), normalized so that: Note that given this normalization condition and the fact that s is a natural number, then D(s; σ, L) ≤ 1.
The numerical algorithm we have used is the so-called brute force algorithm propagating one avalanche at a time [14]. We have studied system sizes ranging from L = 5(L 3 = 125) to L = 48(L 3 = 110592). The mea-Closure relations N = Ns + Nns Nns = Nnsc + Nns0 Ns = N1 + N2 + N3 N3 = N3c + N3− Normalization condition   Table I. The dependence on σ, L and s has been suppressed in order to clarify the Table. The subscript α stands for all the possible sub-indices in Table I. sured properties are always averaged over a large number of realizations of the random field configuration for each value of σ. Typical averages are performed over a number of configurations that ranges between 10 5 for L ≤ 16 to 2000 for L = 48. We have used periodic boundary conditions: the numerical simulations correspond, in fact, to a periodic infinite system. Therefore, strictly speaking, all avalanches are infinite. Nevertheless, we need to identify which avalanches will become important in the thermodynamic limit. The definition that best matches this idea is the concept of spanning avalanches: those avalanches that, at least in one of the x, y or z directions, extend over the length L. This definition is very easy to implement numerically in the brute force algorithm. Spanning avalanches are detected by using three (x, y, z) mask vectors of size L whose elements are set to 0 at the beginning of each avalanche. During the evolution of the avalanche the mask vectors record the shade of the flipping spins along the three perpendicular directions (by changing the 0's to 1's). When the avalanche finishes, it can be classified as being non-spanning, 1d-spanning , 2d-spanning or 3d-spanning depending on the number of such mask vectors that have been totally converted to 1. The number and size distribution of 1d, 2d and 3d-spanning avalanches is also studied and averaged over different realizations of disorder. Table I shows the definitions of avalanche numbers and distributions that will be used throughout the paper. In Table II the list of mathematical relations between the avalanche numbers and distributions is given. We will use the subscript α to indicate any of the avalanche numbers or distributions in Table I.
It should be mentioned that, although the definition of spanning avalanches used in this paper is equivalent to the definition in previous works [13,14,16], the average number of spanning avalanches N s , in some cases, does not coincide with the previous estimations. We guess that the reason is because, in previous works, the method used to count spanning avalanches was averaging twice the 2d-spanning avalanches and was averaging three times the 3d-spanning avalanches. There-fore, in order to compare, for instance with Ref. 13, one should take into account that their number of spanning avalanches N is not equal to the present N s but satisfies: N = (N 1 + 2N 2 + 3N 3 )/3. Moreover, we should point out the following remark before presenting the data. As a consequence of the numerical analysis, several "kinds" of avalanches will be identified (see Table I). Such a separation in different kinds will, in some cases, be justified by the measurement of different physical properties (such as whether the avalanche spans the lattice or not) but, in other cases, will be an "a priori" phenomenological hypothesis to reach a good description of the data. Although some authors will prefer to identify such new "kinds" of avalanches as "corrections to scaling", it will turn out that after the finite size scaling analysis we will be able to identify which different physical properties characterize each "kind" of avalanche. Fig. 1 shows an example of the distribution of avalanches D(s; σ, L) on a log-log scale for three values of σ corresponding to a system with size L = 24. The qualitative behaviour of D(s; σ, L) is that already described in the introduction: when σ is decreased the distribution changes from being approximately exponentially damped (σ > σ c ) to a distribution exhibiting a peak for large values of s (σ < σ c ). Therefore, one can suggest that at the critical value σ c the distribution exhibits power-law behaviour. Nevertheless, it is also evident from Fig. 1 that the finite size of the system masks this excessively simplistic description. Only after convenient finite-size scaling analysis shall we discover which features remain in the thermodynamic limit.

III. NUMERICAL RESULTS
The peak occurring for σ < σ c is basically caused by the existence of spanning avalanches. This is shown in As can be seen the distribution of spanning avalanches D s (s; σ, L) is far from simple. It exhibits a multipeak structure caused by the contributions from D 1 (s; σ, L) , D 2 (s; σ, L) and D 3 (s; σ, L) shown in Fig. 2c. Moreover, D 3 (s; σ, L) itself also exhibits two peaks suggesting that the 3d-spanning avalanches may be of two different kinds. We shall denote critical 3d-spanning avalanches (indicated by the subscript 3c) as those corresponding to the peak on the left and subcritical 3d-spanning avalanches (indicated by the subscript 3−) as those corresponding to the peak on the right. As will be explained below, the 1dspanning avalanches, the 2d-spanning avalanches and the critical 3d-spanning avalanches do not exist in the thermodynamic limit except when σ = σ c . This is the reason for having chosen the word "critical" for this kind of 3d spanning avalanche. It will also be shown that, in the thermodynamic limit, subcritical 3d-spanning avalanches only exist for σ ≤ σ c . As regards the non-spanning avalanches, they will also be classified into two types at the end of this section, although this separation cannot be deduced from the behaviour in Fig. 2 Fig. 3 shows the evolution of D 1 (s, σ, L), D 2 (s, σ, L) and D 3 (s, σ, L) when σ is increased. Note that the righthand peak of D 3 (s; σ, L) shifts to smaller values of s and becomes flat, indicating that the mean size of these subcritical 3d-spanning avalanches decreases. Moreover, above σ c [ Fig. 3(d)] this right-hand peak disappears and a peak on the left emerges.
Besides the normalized distributions, it is also interesting to analyze the actual average numbers of spanning avalanches N 1 (σ, L), N 2 (σ, L) and N 3 (σ, L), which also exhibit singular behaviour at σ c as shown in Fig. 4.
From the direct extrapolation of the data corresponding to different system sizes to L → ∞, we can make the following assumptions: in the thermodynamic limit N 1 (σ) and N 2 (σ) will display a δ-function discontinuity at σ c . N 3 (σ) will display step-like behaviour: for σ < σ c there is only one 3d-spanning avalanche, for σ > σ c there are no 3d-spanning avalanches and at σ = σ c the data supports the assumption that N 3 will also display a δ-function singularity at the edge of the step function. This reinforces the suggestion that there are two different types of 3d-spanning avalanches: as will be shown, in the thermodynamic limit, the number of subcritical 3d spanning avalanches N 3− behaves as a step function, whereas the number of critical avalanches N 3c exhibits divergence at σ c .
The total number of spanning avalanches N s (σ, L) and non-spanning avalanches N ns (σ, L), are displayed in Figs. 5(a) and 5(b) respectively. N s (σ, L) shows, as a result of the divergence of N 3c , N 1 and N 2 , a δ-function singularity at σ c when L → ∞ suggesting that the critical point is characterized by the existence of ∞ spanning avalanches. We would like to point out that previous studies have not clarified this result for the 3d-GRFIM [13] The analysis of N ns is more intricate. Fig. 5(b) shows that N ns (σ, L) grows with σ and L. For large amounts of disorder (σ → ∞) one expects that the hysteresis loop consists of a sequence of non-spanning avalanches of size 1. Therefore, their number will equal L 3 . To reveal this behaviour Fig. 6 shows the dependence of N ns (σ, L)/L 3 as a function of σ. One expects that these lines tend to 1 when σ → ∞. Moreover, a closer look reveals that at σ c ≃ 2.21, there is a contribution to N ns (σ, L)/L 3 which decreases with system size. For low values of σ one expects that non-spanning avalanches always exist, except at σ = 0. This last statement can easily be understood by noticing that an approximate lower bound to the number of non-spanning avalanches can be computed by analyzing how many of the spins S i will flip by themselves, independently of their neighbours, due to the fact that the local field h i is either larger than 6 or smaller than −6. This analysis renders N ns /L 3 > Φ err (6/σ) where Φ err is the error function. From these considerations, we expect that for L → ∞ the curves in Fig. 6 tend to a certain limiting behaviour which increases smoothly from 0 to 1. This can also be appreciated in the inset in Fig. 6, which shows the behaviour of N ns (σ, L)/L 3 as a function of L for four different values of the amount of disorder: σ = 1.7, σ = 2.21 ≃ σ c , σ = 2.5 and σ = 2.7. The four curves exhibit a tendency to extrapolate to a plateau when L → ∞. For the case of σ ≃ σ c an estimation of the extrapolated value is N ns (σ c , L)/L 3 → 0.028.
Consequently, it is necessary to consider the existence of, at least, two kinds of non-spanning avalanches. Those whose number N ns0 increases as L 3 will be denoted as non-critical non-spanning avalanches (with the subscript ns0), and those whose number N nsc increases with L with a smaller exponent will be called critical non-spanning avalanches (with the subscript nsc). In fact, a log-log plot  of N ns (σ c , L)/L 3 − 0.028 versus L provides an estimation for this exponent N nsc (σ c , L) ∼ 0.085L 2.02 . All the assumptions that have been presented, corresponding to behaviour in the thermodynamic limit, will be confirmed by the finite-size scaling analysis presented in the following sections.

IV. RENORMALIZATION GROUP AND SCALING VARIABLES
The basic hypothesis for the analysis of the above results using RG techniques is the existence of a fixed point in the multidimensional space of Hamiltonian parameters. This fixed point sits on a critical surface which extends along all the irrelevant directions. By changing the two tuneable parameters σ and H, the critical surface can be crossed at the critical point (σ c , H c ). As has been explained in the Introduction concerning the analysis of the avalanche number and size distributions, the dependence along the external field direction H has been integrated out. One expects that such integration may distort some of the exponents and the shape of scaling functions, but not the possibility of an RG analysis. This is because the integration range crosses the critical surface where the divergences occur. For a L → ∞ system we assume a unique scaling variable u(σ) which measures the distance to σ c . The dependence of u on σ should be smooth, but its proper form is unknown [18]. We will discuss three different possibilities: 1. The standard choice is to use a dimensionless first approximation by expanding u(σ) as: Nevertheless, in general, the correct scaling variables may have a different dependence on σ. For instance, this may be due to the existence of other relevant parameters, such as the external field, which has been integrated out.

2.
A second choice is to extend the expansion of u(σ) to second order by including a fitting amplitude A: 3. A third choice, which has been used in previous analyses and may be "phenomenologically" justified is: Note that the Taylor expansion of this function is: Fig. 7, shows the behaviour of the three scaling variables u 1 (σ), u 2 (σ) and u 3 (σ). For the representation of u 2 we have chosen A = −0.2, which is the result that we will fit in the following sections. The three choices are equivalently close enough to the critical point. Nevertheless, the amplitude of the critical zone, where the scaling relations are valid, may be quite different. Since A < 0, the variable u 2 cannot be used for σ ≫ σ c since u 2 (σ) shows a maximum at σ = 7.735 = 3.5σ c . A similar problem occurs with u 3 since, due to its asymptotic behaviour (u 3 → 1 for σ → ∞), systems with a large value of σ cannot be distinguished one from another. For the finite system, the magnitudes presented in Table I, depend on σ, L and, in the case of the size distributions, on s. In order to identify the scaling variables, let us consider a renormalization step of a factor b close to the critical point [19,20], such that lengths behave as: (The variables with the b subscript correspond to the renormalized system). We expect that after re-scaling the variable u, measuring the distance between σ and σ c changes as: which is the standard definition of the exponent ν which characterizes the divergence of the correlation length when σ → σ c . Under the same renormalization step we assume that: This latter equation introduces an exponent d α (which has been called 1/νσ by other authors [8]) and can be interpreted as the fractal dimension of the avalanches. As mentioned in the previous section, we expect to find different types of avalanches. As will be shown numerically from the scaling plots in the following sections, it is possible to assume that the different types of avalanches behave with the same fractal dimension d α = d f , except for subcritical 3d-spanning avalanches (for which d 3− = d f ) and non-critical non-spanning avalanches.
Close to the critical point the system exhibits invariance under re-scaling. Therefore, in order to propose a scaling hypothesis of the numbers of avalanches N α and the avalanche size distributions D α , it is important to construct combinations of the variables u, L and s, which remain invariant after renormalization. We find: Note that these three invariant quantities are not independent since equation (12) corresponds to equation (14) multiplied by equation (13) to the power of −1/νd α .

V. SCALING OF THE NUMBERS OF AVALANCHES Nα(σ, L)
The discussion in the previous section, enables us to propose the following scaling hypothesis: The exponent θ α characterizes the divergence of the avalanche numbers at the critical point when L → ∞.
Note that this definition of θ α (which is the same used in previous works [13]) is not consistent with the standard finite-size scaling criterion for which the magnitudes grow with exponents divided by ν. [19,20,21]. As will be shown, the behaviour of the number of 1dspanning avalanches, 2d-spanning avalanches and critical 3d-spanning avalanches can be described with the same value of θ 1 = θ 2 = θ 3c = θ, so that: N 3c (σ, L) = L θÑ 3c uL 1/ν .
We have tried, without success, to scale the number of critical non-spanning avalanches with the same exponent θ. We therefore need to define a different exponent θ nsc , so that: As regards the number of N 3− avalanches, which is different from zero away from the critical point in the thermodynamic limit, we propose a scaling hypothesis that is compatible with the limiting behaviour at σ = 0 and σ = ∞. This leads us to the following assumption: since in the absence of disorder we expect that the hysteresis loop displays a single avalanche of size L 3 , and, consequently the number of avalanches must be N 3− = 1 independent of the value of L. As regards N ns0 it has already been discussed that such avalanches will exist in the thermodynamic limit for all values of σ. Moreover, they are probably not related to critical phenomena at σ c . For this reason we propose the following non-critical dependence: In particular, as already mentioned, for large values of disorder (σ → +∞) these avalanches will be of size s = 1, and their number will be N ns0 (∞) = L 3 . It should also be mentioned that the scaling equations (15) admit alternate expressions by extracting the variable uL −1/ν with the appropriate power so that it cancels out the dependence on L: Nevertheless, such expressions are not very useful for the scaling analysis close to σ c since they will display a large statistical error due to the fact that u → 0 when σ → σ c . Figs. 8 and 9 show the best collapses corresponding to equations (16) and (17) with the three different choices for the variable u, explained in section IV. Data corresponding to L = 5, 8, 10, 12, 16, 24, 32 and 48 have been used. The quality of the collapses close to σ c is quite good in the three cases. The values of the free parameters that optimize each collapse are indicated on the plots. By visual comparison one can see that u 2 is the best choice since it allows the smaller sizes to collapse too. Of course, this is because the collapses in this case have an extra free-parameter A. As regards the quality of the overlaps, no remarkable differences are observed between the choices u 1 and u 3 . In the following collapses we will use u 2 with A = −0.2. Thus, the best estimations of the free parameters are: σ c = 2.21 ± 0.02, ν = 1.2 ± 0.1 and θ = 0.10 ± 0.02.
The procedure for improving the collapse of the data corresponding to different system sizes, which will be used many times throughout this paper, renders what we will call "the best values" of the free parameters. Error bars represent the estimated range of values for which the collapses are satisfactory. We would like to note that the obtained value of σ c (for the three choices of the variable u) is slightly higher than the value σ c = 2.16 ± 0.03 proposed in Ref. 13.
From the fact that the scaling functions in Figs. 8(b) and 9(b) are bounded and go exponentially to zero for u 2 L −1/ν → ±∞ (as can also be checked from a log-linear  plot) one can deduce that, in the thermodynamic limit, 1d-spanning avalanches and 2d-spanning avalanches only exist at σ = σ c . Their numbers increase as L 0.10 with amplitudesÑ 1 (0) = 0.12 ± 0.01 andÑ 2 (0) = 0.07 ± 0.01. Moreover, the peaks of the scaling functionsÑ 1 andÑ 2 which are displaced from u 2 = 0, account for the fact that for a finite system the maximum number of 1d and 2d spanning avalanches occurs for a certain σ c (L) which shifts towards σ c from above.
As regards the 3d spanning avalanches, according to the previous discussions one must consider the contributions from N 3c and N 3− . From the scaling assumptions (18) and (20) and the last closure relation in Table II   can write: This equation indicates that N 3 (σ, L) cannot be collapsed in a straightforward way. We propose here a method to separate the two contributions in equation (23). This method, which we will call double finite-size scaling (DFSS), will be used several times throughout the paper for the analysis of similar equations. By choosing two systems with sizes L 1 and L 2 and amounts of disorders σ 1 and σ 2 so that u(σ 1 )L Thus, we can check for the collapse of data corresponding to different pairs of (L 1 , L 2 ). From the numerical point of view, the DFSS method works quite well. An analysis of error propagation reveals that the scaling function corresponding to the contribution with a smaller exponent will display more statistical errors. Fig. 10 shows the results of the DFSS analysis of N 3 according to equation (23). The different symbols, in this case, indicate the values of L 1 and L 2 used for each set of data. Fig. 10(a) corresponds toÑ 3− (u 2 L 1/ν ) and Fig. 10(b) corresponds toÑ 3c (u 2 L 1/ν ). It should be emphasised that such collapses are obtained without any free parameter. The values of θ, σ c , ν and A are taken from the previous collapses of N 1 and N 2 .
Again, from the shape of the scaling functions we can deduce the behaviour in the thermodynamic limit: from the crossing points of the scaling functions with the u 2 = 0 axis, we find that N 3c (σ c , L) = (0.16 ± 0.02)L 0.10 and N 3− (σ c , L) = 0.79 ± 0.02. As occurred previously with the number of 2d and 1d spanning avalanches,Ñ 3c can also be very well approximated with a Gaussian function with amplitude a 3c = 0.706 ± 0.005, peak position x 3c = 1.244 ± 0.007 and width w 3c = 0.802 ± 0.009. The fact thatÑ 3c (u 2 L 1/ν ) vanishes exponentially for u 2 L 1/ν → ±∞ confirms that, in the thermodynamic limit, such avalanches only exist at the critical point. Furthermore, from the fact thatÑ 3− tends to 1 and to 0 exponentially fast when u 2 L 1/ν → ±∞ we deduce that one subcritical 3d-spanning avalanche will exist for σ < σ c and there will be none above this value.
To end with the analysis of the number of avalanches, we will separate the two contributions to N ns : In this case the DFSS method cannot be applied sincẽ N nsc andÑ ns0 depend on different variables. A first check of the validity of this hypothesis has already been presented in section III. The fit of equation (26) to the data corresponding to σ = σ c (u = 0) , shown in the inset of Fig. 6, gives estimations of θ nsc ≃ 2.02 ,Ñ ns0 (σ c ) = 0.028, andÑ nsc (0) = 0.085. Furthermore, we can also check that the derivative with respect to σ behaves as:  Fig. 11(a) demonstrates that the data (estimated using a two-point derivative formula) is compatible with this behaviour. The line shows the best fit (with two free parameters:Ñ ′ nsc (0) andÑ ′ ns0 (σ c )) of the function (27) with θ nsc + 1/ν − 3 = −0.15. One obtains N ′ nsc (0) = −0.136 ± 0.011 andÑ ′ ns0 (σ c ) = 0.102 ± 0.003. The good agreement is a test of the dependence with the variables uL 1/ν and σ of the functionsÑ nsc andÑ ns0 respectively. To go further into the analysis of N ns , one must provide some extra hypothesis on the shape of the scaling functions. Given the fact that we have found almost a perfect Gaussian dependence of the scaling func-tionsÑ 1 ,Ñ 2 andÑ 3c one can guess thatÑ nsc will also have a Gaussian dependence. By forcing the Gaussian function to satisfyÑ nsc (0) = 0.085 and the fact that thatÑ ′ nsc (0) = −0.136 (from previous estimations) we end up with a trial function with a single free parame-ter that should be enough to satisfactorily scale the data from Fig. 6. The best collapse is shown in Fig. 11(b) which corresponds toÑ ns0 (σ). The functionÑ nsc used for the collapse is shown in the inset and corresponds to a Gaussian function with amplitude a nsc = 0.085, peak position x nsc = −0.6 and width w nsc = 1.485. It is interesting to note that the peak position of this scaling function occurs at a value u 2 L 1/ν = x nsc < 0 as opposed to the case of the previous scaling functionsÑ 1 ,Ñ 2 andÑ 3c for which the peak position was at u 2 L 1/ν > 0. This indicates that the properties of the 1d, 2d and 3c critical avalanches have opposite shifts with finite size L compared to the nsc critical avalanches.
To end with the analysis of the number of nonspanning avalanches it is interesting to compare the functionÑ ns0 (σ) with the approximate lower bound (Φ err (6/σ)) discussed in section III, which is represented by a continuous line in Fig. 11(b). The difference between the two curves, which becomes bigger when σ increases, is due to the existence of clusters of several spins (not considered in the extremely facile analysis presented here) that flip independently of their neighbours contributing to the number of non-critical non-spanning avalanches.

VI. SCALING OF THE DISTRIBUTIONS OF SIZES Dα(s; σ, L)
Close to the critical point there are different ways to express the invariance of the size distributions corresponding to different choices of a pair of invariants among the three exponents proposed in equations (12), (13) and (14). For any generic distribution D α (s; σ, L) one can write the following nine generic expressions: Although we have used the generic index α, it is evident that such assumptions can only be proposed for the distributions of avalanches of a single kind, i.e. D 1 , D 2 , D 3c , D 3− and D nsc . For the composite distributions D 3 , D s , D ns and D, one expects mixed behaviour, and concerning D ns0 we cannot expect a dependence on uL 1/ν . The exponents τ α could also be different for the different kinds of avalanches, but as will be discussed in the following paragraphs, in all cases τ α = 1 except for τ nsc , which will take a larger value. As argued before, when scaling the numbers of avalanches, the last three expressions (34), (35) and (36) are not very useful for the numerical collapses because they introduce large statistical errors. Moreover, when trying to check the collapses expressed by equations (29) and (32), the two independent variables of the scaling function converge to zero when the critical point is approached. Thus, such a collapse cannot be checked for u = 0. Therefore, the interesting scaling equations are (28), (30), (31) and (33).
The behaviour of the scaling functions is, in some cases, restricted by the normalization conditions. If scaling holds for the whole range of s = 1, · · · , L 3 , from equation (28), one can write: If 0 < d α < 3, by defining a new variable x = sL −dα , for large L, the above expression is transformed into the following integral: For those distributions for which the integral converges, it is necessary that τ α = 1. We expect that this condition can be applied to the cases of D 1 , D 2 , D 3c and D 3− . In these four cases, as can be seen in Fig. 2(c), the distributions show a marked decay in the two limits of s → 0 and s → L 3 . (Note that the plots have logarithmic scales and that D 3c and D 3− correspond to the left-hand and right-hand peaks in D 3 respectively). For the distribution D nsc the exponent τ nsc can be larger than 1 since this distribution may extend into the small s region and convergence of the integral in (38) cannot be ensured. Fig. 12 shows a 3d view of the collapses corresponding toD 1 sL −d f , u 2 L 1/ν . The lines show three cuts of the scaling surface corresponding to u 2 L 1/ν = 1.21, u 2 L 1/ν = 0 and u 2 L 1/ν = −0.56. The collapses of the curves corresponding to the different sizes are satisfactory within statistical error. The only free parameter in this case is d f . The best estimation renders a fractal dimension d f = 2.78 ± 0.05 for such 1d-spanning avalanches. Similar behaviour is obtained forD 2 sL −d f , u 2 L 1/ν . Although, in principle, we have considered d f as a free parameter, the best collapses are obtained with the same value d f = 2.78 as that obtained for the 1d-spanning avalanches.
The analysis of D 3c and D 3− is more difficult. According to the corresponding distribution relation (see Table II), and assuming the scaling hypothesis (18), (20) and (28), one can write: where we have taken into account the fact that for the subcritical 3d-spanning avalanches τ 3d− = 1 and they have a fractal dimension d 3− . Although it is possible to conceive a DFSS treatment to separate the two contributions in (39), the hard numerical effort needed as well as the associated statistical uncertainties make it very difficult. In the next section we will show that it is enough to analyze the behaviour of the k-moments of the distributions to obtain the critical exponents.

VII. SCALING OF THE k-MOMENTS OF THE DISTRIBUTIONS
Besides the scaling of the entire distributions D α (s; σ, L) that exhibit large statistical errors, it is also useful to analyze the behaviour of their k-moments. For those distributions for which the integral in Eq.(38) converges (and, therefore, τ α = 1), we can check the corresponding scaling functions. By using a similar argument as that used for deriving equation (38), we get: (40) As an example of such collapses, we have indicated the behaviour of the scaled first moment of the distribution D 1 (s; σ, L) on the horizontal plane of Fig. 12. In this case the collapses are obtained without any free parameter.
As will be seen later, it is more convenient to analyze the scaling behaviour of the products N α s k α . By using equations (15) and (40), one gets: Fig. 13 shows the collapses corresponding to the first moment (average size) of D 1 and D 2 . No free parameters are used in this case. Similar scaling plots can be obtained from the analysis of the second moments with the same set of scaling exponents.  As regards the scaling of N 3 s 3 , multiplying expression (39) by s, summing over the whole s-range, and imposing condition (37), one obtains: This equation can be separated by a DFSS analysis. Figs. 14(a) and 14(b) show the collapses corresponding toÑ 3c Ψ 1 3c andÑ 3− Ψ 1 3− respectively. The only free parameter in this scaling plot is the fractal dimension of the subcritical 3d-spanning avalanches. The best value is d 3− = 2.98 ± 0.02. Note that the shape of the scaling function in Fig. 14(b) indicates that, in the thermodynamic limit, the critical 3d-spanning avalanches only contribute to the first-moment for σ = σ c On the other hand, the shape of the scaling function in Fig. 14(a) indicates that, in the thermodynamic limit, the subcritical 3d-spanning avalanches may contribute to the first moment in the whole u 2 < 0 range. Note that, as revealed by the inset in Fig. 14(a), the behaviour in the region of negative values of u 2 L 1/ν isÑ 3− Ψ 1 3− ∼ (|u 2 |L 1/ν ) β3− with β 3− = 0.024 ± 0.012. This numerical value is compatible with the equation: This hyperscaling relation, when introduced into equation (42), results in a second term that grows with L 3 . As will be analyzed in the next section, such a term will be responsible for the order parameter behaviour in the thermodynamic limit.
The analysis of the moments of the non-spanning avalanches presents extra difficulties, as occurred in the analysis of their number. The expected behaviour is: As explained previously, the DFSS cannot be applied, given the different dependence on uL 1/ν and σ of the two terms in (44). The possibility of using a trial function is, now, more difficult since we cannot make a straighforward hypothesis on the shape of Ψ k nsc . In order to fit the value of τ nsc and d nsc we can analyze the dependence of the k-moment (for k = 2 and k = 3) and its derivatives with respect to σ at σ = σ c (u = 0). Data is shown in Fig. 15(a) and Fig. 15(b) with log-log scales. The almost perfect power-law behaviour for different values of k and for the derivatives, indicates that the second term in (44) plays no role in σ c . This is because the exponent of the first term is much larger than 3. Indeed, the best fits are obtained with d nsc = d f = 2.78 ± 0.05 and τ nsc = 1.65 ± 0.02 which render large values of the exponent of the first term (> 5.8). Similar fits can be obtained from higher moments with the same values of the exponents d nsc and τ nsc .

VIII. MAGNETIZATION DISCONTINUITY
In this section we discuss the behaviour of the discontinuity ∆m in the magnetization of the hysteresis loop. We would like it to behave as an order parameter. For large systems, it is clear that only spanning avalanches may produce a discontinuity in the magnetization. We can evaluate the total average magnetization jump ∆m s due to the contribution of all the spanning avalanches (1d, 2d, 3c and 3−). Fig. 16(a) shows the behaviour of ∆m s versus σ for different system sizes. According to the scaling analysis in the previous section, ∆m s will behave as: This equation tells us that ∆m s will display a mixed scaling behaviour. The first term in (46) accounts for the contributions of the 1d-spanning, 2d-spanning and critical 3d-spanning avalanches. We can define an exponent β c so that: This relation is the same relation that other authors have called "violation of hyperscaling" [11,13,22]. From our best estimations of θ, ν and d f , we obtain β c = 0.15± 0.08.
At this point, it is interesting to compare equations (43) and (47.) We would like to note that we could also have introduced an exponent θ ′ that would transform Eq. (43) into an equation similar to (47). Nevertheless, the quality of the scalings of the numbers of 3d-spanning avalanches in Fig. 10 shows that such an exponent θ ′ is either zero or very small. Moreover, an analysis of the behaviour ofÑ 3− for uL 1/ν → −∞ reveals an exponential drift versusÑ 3− = 1 which reinforces the idea that there is no need for an hyperscaling exponent θ ′ . Note that a value θ ′ > 0 implies that the number of subcritical 3d-spanning avalanches (3−) will be infinite, in the thermodynamic limit. On the other hand, our assumption that θ ′ = 0 indicates that N 3− behaves as a step function in the thermodynamic limit.
By inserting equations (43) and (47) into (46) one can easily read the mixed scaling behaviour of ∆m s : where, β c /ν = 0.12 and β 3− /ν = 0.02. Φ is a scaling function with a peaked shape (it corresponds to the sum of the scaling functions in Figs. 13(a), 13(b) and 14(b)) and Φ ′ is twice the scaling function in Fig. 14(b). Consequently, in the thermodynamic limit, only the second term associated to the subcritical 3d-spanning avalanches will contribute to the magnetization jump (order parameter). For finite systems, the first term may affect the scaling of the data close to σ c given the peaked shape of Φ.
This behaviour can be observed in Figs. 16(b) and 16(c), where the two possible scalings show the breakdown of the collapse for u 2 L 1/ν < 0 when using the exponent β c /ν and the breakdown of the collapse for u 2 L 1/ν ≃ 0 when using the exponent β 3− /ν. The larger the system, the better will be the data collapse in Fig.  16(c) and the worse will be the collapse in Fig. 16  In both cases, the lines show the best fits of Eq.44 and its derivative at σ = σc. Table III shows a summary of the exponents that characterize the avalanche numbers and distributions obtained from our numerical simulations. We would like to point out that such exponents are independent of σ and L in a very large region around the critical point both for σ > σ c and σ < σ c simultaneously. Such an achievement has not been possible in previous analyses, even with larger system sizes. The reason is that some of the contributions we have identified (namely 3− and ns0), which reduce finite size effects, were previously neglected.

IX. DISCUSSION
In Table III we also indicate previous estimations of the exponents found in the literature [13]. The comparison is quite satisfactory. Let us analyze the eight exponents: • Although the value of ν does not fall within the error bars in Ref. 13  for the collapses may introduce some deviations in this value. By using u 1 we obtain ν = 1.14 and using u 3 we obtain ν = 1.4.
• As regards θ our value is in agreement with the value previously reported [13] (We would like to note that in Ref. 13, the authors also report a value of 0.015 ± 0.015 probably due to a misprint). The fact that this exponent is non-zero implies that there are infinite spanning avalanches at the critical point in the thermodynamic limit.
• As regards θ nsc , to our knowledge there are no previous finite-size scaling analyses of the number of non-spanning avalanches.
• Concerning d f and d 3− , the numerical values are consistent with the value d f = 2.98 ± 0.43 estimated previously [10]. We shall note that this previous estimation was obtained from the analysis of the distributions of non-spanning avalanches. It should therefore correspond to our d f and not to our d 3− (which corresponds to the subcritical 3dspanning avalanches). Note also that the difference between d f and d 3− suggests that there might be real physical differences between such two kinds of avalanches. The possibility of distinguishing them in the numerical simulation will be studied in a future work.
• The exponent τ nsc , according to our definitions, describes the scaling behaviour of the distribution of critical non-spanning avalanches. Previous measurements of a similar exponent, have analyzed N ns without distinguishing between critical (nsc) and non-critical (ns0) non-spanning avalanches and have not considered the fact that the system is finite. We can estimate what the value of an effective exponent τ ef f will be for the distribution of nonspanning avalanches for very large systems. From equation (44), taking k > 1, and large values of L, only the first term in the sum survives, so that: N ns (σ) s k ns (σ) = = L θnsc+(1+k−τnsc)d fÑ nsc uL 1/ν Ψ k nsc uL 1/ν .
On the other hand, in the same limit, the analysis of equation (26) renders: Combining the last two equations, we get an estimation for the pseudo-scaling behaviour of the k-moment of the non-spanning avalanches: If one approximatesÑ ns0 (σ) byÑ ns0 (σ c ) and imposes s k ns ∼ L −(τ ef f −k−1)d f S k (uL 1/ν ) it is possible to deduce that the effective exponent is τ ef f = τ nsc + (3 − θ nsc ) /d f . From our numerical estimations of the different exponents in Table III, one obtains τ ef f = 2.00 ± 0.06. This value is in very good agreement with the value τ ef f = 2.03 ± 0.03 found previously [13]. Nevertheless, we would like to point out that according to our analysis, such an exponent is not a real critical exponent and, therefore, will depend on σ for finite systems as has been found previously [16].
• As regards the values of β c and β 3− we would like to note that previous analyses have not identified the two contributions to ∆m s . It is therefore not strange that different values have been obtained previously: 0.17 ± 0.07 [8], 0.0 ± 0.43 [11], 0.035±0.028 [13]. The larger the system, the closer the effective exponent becomes to β 3− .
Finally, it is interesting to compare the behaviour of spanning avalanches, with the problem of percolation [23]. In percolation, the number of percolating clusters behaves as a step function, in the thermodynamic limit for d < 6, exactly as N 3− . The order parameter is, in this case, the probability for a site to belong to the percolating cluster. However, this is precisely what we are evaluating by the function N 3− s 3− /L 3 which is the second term in (46) and is the only relevant term in the thermodynamic limit. As occurs in percolation, the hyperscaling relation (43) between β 3− , ν and d 3− is fulfilled since only one infinite avalanche contributes to the order parameter for σ → σ c from below. Moreover, in the percolation problem for d > 6 [24], the number of percolating clusters exhibits, besides the step function, an extra δfunction singularity at the percolation threshold. In our case (3d-GRFIM) we also have such a contribution at σ = σ c , which we have identified as 1d, 2d and critical 3d-spanning avalanches [the first term in (46)]. The existence of such an infinite number of avalanches exactly at σ c (the number of which grows as L θ ) implies the breakdown of the hyperscaling relation β c = ν[3 − (θ + d f )].

X. SUMMARY AND CONCLUSIONS
In this paper we have presented finite-size scaling analysis of the avalanche numbers and avalanche distributions in the 3d-GRFIM with metastable dynamics. After proposing a number of plausible scaling hypotheses, we have confirmed them by obtaining very good collapses of the numerical data corresponding to systems with sizes up to L = 48.
The first result is that, in order to obtain a good description of the numerical data, one needs to distinguish between different kinds of avalanches which behave differently when the system size is increased. Avalanches are classified as being: non-spanning, 1d-spanning, 2dspanning or 3d-spanning. Furthermore, we have shown that the 3d-spanning avalanches must be separated into two classes: subcritical 3d-spanning avalanches with fractal dimension d 3− = 2.98 and critical 3d-spanning avalanches with fractal dimension d f = 2.78, as the 1d− and 2d−spanning avalanches. Non-spanning avalanches occur for the whole range of σ. We have also proposed a separation between critical non-spanning avalanches and non-critical non-spanning avalanches in order to obtain good finite-size scaling collapses. The non-critical nonspanning avalanches are those whose size is independent of the system size and whose number scales as L 3 . The critical non-spanning avalanches also have a fractal dimension d f = 2.78.
The second important result, is the scenario for the behaviour in the thermodynamic limit: below the critical point, there is only one subcritical 3d-spanning avalanche, which is responsible for the discontinuity of the hysteresis loop. Furthermore, at the critical point there are an infinite number of 1d-, 2d-, and 3d-critical spanning avalanches.
For finite systems, the six different kinds of avalanches can exist above, exactly at and below σ c . The finitesize scaling analysis we have performed has also enabled us to compare different scaling variables u, which measure the distance between the amount of disorder in the system σ and the critical amount of disorder σ c . The best collapses are obtained using the variable u 2 = (σ − σ c )/σ c + A [(σ − σ c )/σ c ] 2 with A = −0.2. So far the analysis presented in this paper, is restricted to the analysis of the numbers and distributions of avalanches integrated along half a hysteresis loop. Our analysis of the average magnetization discontinuity ∆m starts from the hypothesis that only spanning avalanches may contribute to such a discontinuity. However, as a future study, we suggest that the measurement of correlations in the sequence of avalanches and the analysis of non-integrated distributions, may reveal details of the singular behaviour at the critical field H c . For instance, non-spanning avalanches could show a tendency to accumulate in H c , in the thermodynamic limit. This could change some of the conclusions reached in this work.
As a final general conclusion we have shown that it is not necessary to simulate very large system sizes to estimate the critical exponents for this model. In order to identify the different kinds of avalanches it may even be better to analyze small systems with larger statistics.

XI. ACKNOWLEDGEMENTS
We acknowledge fruitful discussions with Ll.Mañosa and A.Planes. We acknowledge an anonymous referee for helping us to understand the difference between the present method of counting avalanches and the method used in previous works. This work has received financial support from CICyT (Spain), project MAT2001-3251 and CIRIT (Catalonia) , project 2000SGR00025. F.J. P. also acknowledges financial support from DGICyT. Scaling of ∆ms by considering the 1d, 2d, and the critical 3dspanning avalanches. Note the lack of collapse for the region u2L 1/ν 0. (c) Scaling of ∆ms by considering the subcritical 3d-spanning avalanches. Note the lack of collapse for the region u2L 1/ν ∼ 0. Symbols indicate the system sizes according to the legend.