Hysteresis and avalanches in the T=0 random-field Ising model with 2-spin-flip dynamics

We study the non-equilibrium behavior of the three-dimensional Gaussian random-field Ising model at T=0 in the presence of a uniform external field using a 2-spin-flip dynamics. The deterministic, history-dependent evolution of the system is compared with the one obtained with the standard 1-spin-flip dynamics used in previous studies of the model. The change in the dynamics yields a significant suppression of coercivity, but the distribution of avalanches (in number and size) stays remarkably similar, except for the largest ones that are responsible for the jump in the saturation magnetization curve at low disorder in the thermodynamic limit. By performing a finite-size scaling study, we find strong evidence that the change in the dynamics does not modify the universality class of the disorder-induced phase transition.


I. INTRODUCTION
The non-equilibrium random field Ising model (RFIM) was introduced by Sethna et al. [1] as a model for the Barkhausen effect in ferromagnets and more generally as a prototype for many experimental systems that show hysteretic and jerky behavior when driven by an external force. Because of the presence of disorder, these systems have a "complex free energy landscape" with a multitude of local minima (or metastable states) separated by sizeable barriers, which makes thermally activated processes essentially irrelevant at low enough temperature (the lifetime of metastable states may then be considered as infinite). As a consequence, these systems remain far from equilibrium on the experimental time scales (even when the driving rate goes to zero) and their response to the external force is made of a series of jumps (avalanches) between neighboring metastable states. This type of behavior is very well modeled by the ferromagnetic RFIM with a zero-temperature single-spinflip dynamics in which a spin flips only if this lowers its energy. The local character of the energy minimization is then at the origin of irreversibility. With this dynamics, the RFIM satisfies the property of return-point memory (or "wiping out" effect) which is a feature observed in several experimental systems with good approximation. Moreover, in dimension d ≥ 3, the model is known to exhibit an out-of-equilibrium phase transition between a strong-disorder regime where the magnetization hysteresis loop is smooth on the macroscopic scale and a weakdisorder one where it has a discontinuous jump. Such a transition has been observed in thin Co/CoO films [2] and * Electronic address: eduard@ecm.ub.es Cu-Al-Mn alloys [3], and it has been recently suggested [4] that it may also be associated to the change in the adsorption behavior of 4 He in dilute silica aerogels [5]. The two regimes, strong and weak disorder, are separated by a critical point characterized by universal exponents and scaling laws which have been extensively studied by analytical and numerical methods [6,7]. In particular, much effort has been recently devoted to analyze the number of avalanches, their size, and their geometrical properties above, below, and at criticality [8,9]. These results provide a comprehensive, though rather complex scenario for the phase diagram of the non-equilibrium RFIM with a metastable dynamics in the thermodynamic limit. A recent discussion of the relevance of the model to the description of the Barkhausen effect in real magnets can be found in Refs. [10,11].
An issue that so far has not been studied is the robustness of this theoretical description with respect to a change in the dynamics (in the literature on the RFIM, it is implicitely taken for granted that this should not matter). The single-spin-flip dynamics, however, is not the unique (and may be not the best) way of simulating hysteretic dynamical processes in actual systems. It is clear for instance that the hysteresis loop will shrink if the dynamics allows for a better equilibration of the system by employing multiple-spin flips. Then, what will be the avalanche properties ? Will there still be a phase transition ? If so, will the critical behavior be the same as with the single-spin-flip dynamics ? There is in fact the intriguing possibility, supported by numerical simulations and analytical arguments [6,9,12], that the non-equilibrium and equilibrium transitions of the T = 0 RFIM belong to the same universality class, even if criticality occurs in zero external field at equilibrium and at a non-zero coercive field in the irreversible evolution [13]. Since the ground state is stable with respect to the flip of an arbitrary (finite) number of spins, this may indicate that the disorder-induced transition has a universal character at criticality which does not depend on the specific choice of the dynamics [9].
In order to shed some light on this issue and check the robustness of the transition, we study here the nonequilibrium T = 0 RFIM with a 2-spin-flip dynamics. We compare the results with those obtained with the standard 1-spin-flip dynamics, in particular those concerning the number and size of the avalanches. We first show in section II that one can indeed define a 2-spin-flip algorithm that yields a deterministic evolution of the system with the external field. In particular, the dynamics satisfies the "abelian" property which guarantees that the same 2-spin-flip stable configuration is attained, whatever the order in which the unstable spins are relaxed during an avalanche. It also has the property of returnpoint memory. In section III we present the results of our numerical simulations on 3-d lattices and compare them to the behavior of the same system with a 1-spinflip dynamics. The hysteresis loops are significantly reduced when allowing 2-spin flips, but the main features, including the presence of a disorder-induced transition with an associated critical point, are not significantly altered. We then perform in section IV a finite-size scaling analysis, which allows for a determination of the critical properties. We find that, within statistical uncertainty, the exponents and scaling functions are identical to those obtained with the standard 1-spin-flip dynamics. The main conclusions of the study are reported in section V.

II. MODEL AND DYNAMICS
The model is defined on a cubic lattice of linear size L with periodic boundary conditions. On each site (i = 1, . . . , N = L 3 ) there is an Ising spin variable (S i = ±1).

The Hamiltonian is
where the first sum extends over all distinct pairs of nearest-neighbor (n.n), H is the external applied field, and h i are quenched random fields drawn independently from a Gaussian distribution with zero mean and standard deviation σ. We are interested in studying the sequence of states along irreversible paths at T = 0 when the system is driven by the external field H (in the adiabatic limit corresponding to a vanishingly small rate of change of the external field). For this purpose, the Hamiltonian must be supplemented by some dynamical rules.
The standard 1-spin-flip dynamics used in previous studies consists in minimizing the energy of each spin, where is the net field at site i (the summation in Eq. (3) is over the z n.n. of i). From the above expression, it is clear that the minimization of H i is obtained by aligning each spin with its local field, S i = sign(f i ), as represented schematically in Fig. 1(a). This provides a stability criterium for any state with respect ot this 1-spin-flip dynamics. This dynamical rule may be implemented by an algorithm that propagates one avalanche at a time [14]. Starting from a stable configuration, one increases (or decreases) the external field until the local field f i at some site i becomes zero (this corresponds to the vanishing of the local minimum in which the system was trapped). The spin S i (which is uniquely defined because the distribution of the random fields is continuous) is then flipped, which in turn may cause neighboring spins to become unstable and thus initiate an avalanche. The avalanche stops when a new metastable state is reached. The external field is then changed again, and so on. When the spins that become unstable during the avalanche are sequentially reversed (e.g., by increasing i from 1 to N ), it is of course crucial that the final state does not depend on the sequential order. Thanks to the ferromagnetic nature of the couplings, this is indeed the case as a result of the so-called "no-passing" and "abelian" properties of this dynamics [1,15]. Moreover, the same state is also reached when all unstable spins are flipped in parallel, which allows to measure the "time" it takes an avalanche to occur [7]. Setting the rules for a 2-spin-flip dynamics is rather straightforward. By definition, 2-spin-flip stable states are spin configurations whose energy (defined solely by the Hamiltonian) cannot be lowered by the flip of one or two spins (clearly, new features are only introduced when these two spins are n.n.). The (local) energy to be minimized is thus the one associated with a pair ij of n.n. spins, where is the field experienced by S i without the influence of the neighbor S j (the summation in Eq. (5) is over the n.n. of i excepting j). One can then think of the dynamics as made of single-spin flips and "irreducibly cooperative" 2-spin flips. As pictured in Fig. 1(b), a single-spin flip occurs whenever the net field on This corresponds to the changes ↓↓↔↑↓, ↓↓↔↓↑, ↑↑↔↑↓, ↑↑↔↓↑ in the diagram. An irreducibly cooperative 2-spin flip involves a n.n. pair of spins with the same sign (↑↑ or ↓↓) that cannot flip individually (i.e. without the simultaneous flip of the neighbor). This occurs whenever the net field on the pair of aligned spins, These two last conditions come from the fact that the spins cannot individually flip (hence neither f i nor f j changes sign before the cooperative flip of the pair) and the condition that once a pair has flipped, none of the spins can flip back individually.
The corresponding algorithm is a simple extension of the one described above for the 1-spin-flip dynamics. Starting from a 2-spin-flip stable configuration, the external field is varied until one finds a pair of spins that becomes marginally stable: the representative point of this pair in the diagram of Fig. 1(b) leaves the region where it was originally (associated with ↑↑, ↑↓, ↓↑, or ↓↓) and, depending on the border which is first attained, only one spin flips or the two spins flip simultaneously [16].
It is not hard to show that this new dynamics obeys the same properties as the 1-spin-flip dynamics, in particular the crucial abelian property. This is again a consequence of the ferromagnetic nature of the interactions. One only needs to note that the state of the system can be represented by a set of zN/2 points (corresponding to all the distinct n.n. pairs) in the diagram of Fig. 1(b). As the external field is monotonously increased (resp. decreased), the local fields can only increase (resp. decrease), the spins can only flip up (resp. down), and the points can only move up and right (resp. down and left) in the diagram. By slightly modifying the arguments of Ref. [1], one then can prove the no-passing rule, the abelian property and the existence of return-point memory. Instead of paraphrasing the demonstrations given in Ref. [1], we choose here to illustrate these properties by a numerical example. (We have also performed numerical tests in many situations and found no violations of these properties.) The evolution of a system with size L = 30 and σ = 2.5 is shown in Fig. 2 where the energy per spin, e = H/N , is plotted as a function of the magnetization m = S i /N (strictly speaking, e is the enthalpy). The external field H is varied from a very large initial value where all spins are up to the final value H = −2 (here and after, J is taken as the energy unit). Two of the curves display the sequence of unstable states that are obtained after a sudden change of the external field using either a sequential or a parallel updating algorithm. We also show the metastable evolution corresponding to the adiabatic driving (with sequential updating) along the hysteresis loop. In all cases the final state is the same. This is true even when the intermediate states are distinct, for instance when the spins are chosen sequentially in a different order.
The property of return-point memory property is illustrated in Fig. 3, again for a system with size L = 30 and σ = 2.5. The minor loop is obtained by reversing the evolution of the external field first in the decreasing branch at H = −0.85 and then at H = +0.80. As can be seen from the inset, the internal loop closes before the return point, so that the evolution follows that of the major loop for a small region of H −0.85. We now present the results of numerical simulations performed on 3-dimensional cubic lattices using either the 1-spin-flip or 2-spin-flip dynamics with sequential updating. Averaged quantities were obtained with statistics over 10 3 to 10 5 different realizations of the random field distribution and system sizes ranging from L = 8 up to L = 48. As emphasized in Ref. [8], in order to describe properly avalanche properties (especially the "spanning" avalanches), it is more important to perform averages over many disorder realizations than to simulate very large system sizes.  Fig. 4 shows the hysteresis loops obtained in a single sample for two different values of the disorder σ. For comparison, we also display the magnetization curves obtained with the 1-spin-flip dynamics and with the algorithm of Ref. [17] which gives the exact ground-state (equilibrium) magnetization. The corresponding behavior of the enthalpy per spin along the ascending branches of the loops is reported in Fig. 5 (note in passing that the ground-state enthalpy does not show any discontinuity as the external field is varied [17]). As could be expected, the main effect of the new dynamics is to reduce the size of the hysteresis loops. Specifically, the coercivity (i.e., the magnitude of the external field for which the magnetization is equal to zero) is decreased by more than 30% when allowing pairs of spins to flip together. Accordingly, the enthalpy difference between the ground state and the metastable states which are visited along the loops is also reduced. Nevertheless, the loops display the same key feature, that is a change from a discontinuous to a continuous behavior as the disorder is increased. This suggests that there is also an out-ofequilibrium disorder-induced phase transition under the 2-spin-flip dynamics, with a critical value of σ at which the discontinuity appears for the first time in the thermodynamic limit.

III. NUMERICAL SIMULATIONS
It is worth pointing out that the new dynamical rules allows the system to effectively overcome energy barriers of magnitude up to ∆E = 2J. Indeed, the difference between the two dynamics shows up when a pair of n.n. spins with same sign can cooperatively flip, say from ↓↓ to ↑↑ when the external field is increased, whereas each of its spin cannot individually flip. This means that along 1-spin-flip paths, the system has now been able to bypass the higher-energy states, either ↑↓ or ↓↑. By using Eq. (4), it is easy to show that the relevant barrier height associated with this process is at most 2J. Since cooperative flips of more than 2 spins do not occur with the chosen dynamics and the system's trajectory otherwise go through states of decreasing energy, one concludes that ∆E = 2J is the maximum barrier height that the system may overcome when passing from the 1-spin-flip to the 2-spin-flip dynamics.

B. Avalanches
As shown in recent studies [8,9], a good characterization of the disorder-induced critical point can be reached by analyzing the number and size distribution of the magnetization jumps (avalanches) that compose the hysteresis loops in finite systems. For that purpose, it is necessary to classify the avalanches in several categories, according to their behavior as the system size L is increased. One first has to distinguish whether or not an avalanche spans the system from one side to the other, in 1, 2 or 3 spatial directions (indicated in the following by the index α). For each individual avalanche, this is a property that can be easily detected during the simulation. Avalanches are thus classified as being nonspanning (α =ns), 1d-spanning (α = 1), 2d-spanning (α = 2), or 3d-spanning (α = 3). Fig. 6 shows the number of 1d, 2d and 3d-spanning avalanches recorded along the descending branch of the hysteresis loops as a function of σ. The data, averaged over disorder, correspond to a system of size L = 24. It can be seen that the behavior of the three quantities is completely equivalent under the two dynamics. The only difference is a shift towards larger values of σ when the 2-spin-flip dynamics is used. The same shift is also found for all studied system sizes. This is a first indication that σ (2) c > σ (1) c , as will be confirmed by the finite-size scaling analysis presented in the next section. In the case of the 1-spin-flip dynamics, a detailed analysis was performed in Refs. [8,9], revealing the scenario that occurs in the thermodynamic limit and that is already suggested by the data shown in Fig. 6: when L → ∞, N 1 (σ) and N 2 (σ) are expected to display a δ-singularity at σ c and N 3 (σ) a step-like behavior. It was shown, moreover, that there are two types of 3d-spanning avalanches, subcritical and critical, which scale with different exponents. The former are responsible for the discontinuity in the magnetization curve in the thermodynamic limit (there is only one, compact, subcritical avalanche for σ < σ c ) whereas the latter only exist at σ c (hence the additional δ-singularity at the edge of the step function whose signature is already visible in Fig. 6). In a finite system, however, all kinds of avalanches may exist close enough to the critical point, and it is quite difficult to discriminate subcritical from critical avalanches. In Ref. [9], an elaborate analysis was needed to show that these avalanches have different fractal dimensions at criticality. This study is impossible here because of the complexity of the 2-spin-flip algorithm that forbids the use of large systems with good enough statistics. Therefore, in the following, we shall not distinguish between these subcritical and critical 3dspanning avalanches.   figure) induced by the 2-spin-flip dynamics is that there are a little less avalanches of size s = 1 and a little more avalanches of size s = 2, but the rest of the distribution is almost the same. In particular, with both dynamics, the expected power-law behavior of the distribution will be characterized by the same exponent τ ef f ≈ 2.0 in the thermodynamic limit [9].
The size distributions D 1 and D 2 of the 1d and 2dspanning avalanches shown in Figs. 7(d) to 7(g) also appear to be identical with the two dynamics, at least within statistical error bars. (Note that these avalanches do not exist for σ = 1.6 because this value is much lower than σ c for both dynamics.) The only visible differences between the two dynamics occur in D 3 , the size distribution of the 3d-spanning avalanches. Specifically, for σ = 1.60 and 2.25 (i.e. below σ c and very close to σ c , respectively), the large 3d-spanning avalanches tend to be shifted to even larger sizes. According to Refs. [8,9], these avalanches are probably subcritical spanning avalanches and their average size is thus a measure of the order parameter. Therefore, this result is another indication that σ c . In contrast, for σ = 2.80 (which is clearly above σ c ), the distribution D 3 is not affected by the dynamics (Fig.7(j)): in this case, one expects to detect only critical 3d-spanning avalanches in a finite system.

IV. FINITE-SIZE SCALING ANALYSIS
As already mentioned, a precise determination of σ (2) c and of some of the critical exponents requires a detailed finite-size scaling analysis which, unfortunately, is not posssible with the present algorithm. What can be done, however, is to check whether the set of exponents found with the 1-spin-flip dynamics can also be used to scale the new 2-spin-flip data. We first consider N 1 and N 2 , the numbers of 1d and 2d-spanning avalanches, which are quantities that have a simple scaling behavior. Following previous work [8,9], we assume the forms where θ and ν are critical exponents,Ñ 1 andÑ 2 are scaling functions, and u 2 is a scaling variable that measures the distance to the critical point. It is defined as where the parameter A accounts for a second order correction that plays a role when the studied systems are not very large, as is the case here. In Ref. [8], the best choice for the collapse of the scaling plots was obtained with ν = 1.2, θ = 0.1, and A = −0.2, values that we keep here. The only free parameter is thus σ (2) c . As shown in Fig. 8, a very good collapse of the new data can be obtained with σ (2) c = 2.25. Taking into account the quality of the plots, there is an uncertainty of ±0.01 on this value, but, clearly, the value σ (1) c = 2.21 obtained with the 1-spin-flip dynamics [8] can be discarded. We have also checked that this conclusion is not modified when the non-universal parameter A is allowed to vary. These scaling plots thus indicate that the 2-spin-flip dynamics only induces a shift of the critical value of the disorder but does not change the universality class of the transition.
A second check of this result may be obtained by measuring the average field H 3 (σ) at which the 3d-spanning avalanches occur. This quantity was studied with the 1-spin-flip dynamics in Ref. [9], which allowed to map out the first-order line (corresponding to the macroscopic jump in the magnetization) in the diagram H − σ. The dependence of H 3 (σ) with the system size is shown in Fig. 9(a). Note that the transition line extends above σ c in a finite system. However, the end-point, beyond which no 3d-spanning avalanches are found (with the present sampling of disorder realizations), becomes closer and closer to the critical point (H c , σ c ) as L is increased.
According to Ref. [9], this set of curves should scale as is the critical field, B ′ a non-universal tilting constant, µ a critical exponent, andĥ 3 the corresponding scaling function. Strictly speaking, the scaling should be done separately for the average fields H 3− and H 3c at which the subcritical and critical 3d-spanning avalanches occur [9]. Indeed, the number of these avalanches scales differently with L. However we expect the lack of scaling to have a (small) effect only in the region u 2 L 1/ν ≃ 1 where the two kinds of avalanches coexist [9].
In Ref. [9], the best collapse of the 1-spin-flip data was obtained with µ = 1.5, B ′ = 0.25, and H (2) c = 2.25, and consider the critical field as the only free parameter. As shown in Fig. 9(b), a very good collapse can be obtained with H (2) c = −0.885 (for clarity, we only present the collapse for σ < σ (2) c ). Again, this result is consistent with the assumption that the new dynamics does not change the universality class of the transition. On the other hand, the significant decrease in the critical field (in absolute value) is in line with the decrease in coercivity illustrated by Fig. 4.
Finally, it is interesting to study the influence of the new dynamics on the finite-size scaling functions. In Fig. 10, we compare the scaling collapses of the number of 1d-spanning avalanches,Ñ 1 (σ, L), obtained with the two dynamics (for the 2-spin-flip dynamics, this is the same curve as in Fig. 8). As can be seen, the agreement between the two curves is quite remarkable. Even the deviations (around u 2 L 1/ν ≈ 0) from the Gaussian fit proposed in Ref. [8] are the same. A similar analysis forÑ 2 (σ, L) shows the same agreement. This is again a strong indication for universality, that goes beyond the equality of the critical exponents [18].

V. CONCLUSION
We have shown in this paper that the non-equilibrium behavior of the 3-d RFIM at zero temperature is not qualitatively altered when going from the standard 1spin-flip metastable dynamics to a 2-spin-flip metastable dynamics. The coercivity (that is essentially the width of the hysteresis loop associated with the evolution of the magnetization with the driving field) is significantly reduced by allowing 2-spin flips, but the main features of the hysteretic behavior, including the presence of a disorder-induced non-equilibrium transition and the distribution of avalanches, remain similar. By using a finite- size scaling analysis, focused on the number and size of the avalanches, we have furthermore provided strong evidence that the critical behavior (exponents and scaling functions) obtained with the 2-spin-flip dynamics is in the same universality class as that obtained with the 1spin-flip dynamics.
Changing the dynamical rules used to study the evolution of the model, as done here, helps adressing several important questions. The first one is the robustness of the hysteretic scenario provided by the T = 0 RFIM with the standard 1-spin-flip dynamics. As discussed in introduction, a basic assumption underlying the theoretical description is that the system gets trapped in metastable states on the experimental time scale and can only escape when a change in the external field makes the relevant state looses its stability. However, real systems are not at T = 0 and some partial, local equilibration, due for instance to thermally activated processes, may take place even though the system remains far from equilibrium on the experimental time scale. Introducing cooperative 2-spin flips is a simple way to check the effect of partial equilibration processes. Our results clearly point towards the robustness of the whole theoretical picture drawn from previous studies of the model using the T = 0 1-spin-flip dynamics [1,6,7,8,9].
A second question concerns the relation between the disorder-induced critical properties observed in the nonequilibrium behavior of the RFIM at T = 0 and the equilibrium critical behavior associated with the paramagnetic to ferromagnetic transition. There is numerical evidence that, at least to a good approximation, critical exponents and scaling functions associated with the two kinds of criticality are the same [6,9,12]. (Addi-tional, but less conclusive evidence is provided by exact, but mean-field-like results on the Bethe lattice [12] and by perturbation theory near the upper critical dimension d = 6 [6].) Our present findings suggest that a whole series of T = 0 metastable dynamics involving k-spin flips lead to the same non-equilibrium criticality, with the critical disorder strength increasing with k and the critical coercive field decreasing with k. Equilibrium behavior at T = 0 as a function of the external field involves the system's ground state, i.e. a state stable to flips of any arbitrary finite number of spins. Based on the numerical closeness of the critical exponents and scaling functions, on the fact that the critical disorder strength satisfies σ eq c > σ c ), and on the similarity of the underlying physics at T = 0, it is tempting to speculate that the equilibrium behavior can be obtained as the limit of a series of k-spin-flip metastable dynamics with increasing k, with the critical properties of the whole series belonging to the same universality class and governed by the same fixed point [9].