PHYSICAL REVIEW B 71, 134424 s2005d Hysteresis and avalanches in the T = 0 random-field Ising model with two-spin-flip dynamics Eduard Vives* Departamento d’Estructura i Constituents de la Matèria, Universitat de Barcelona Diagonal 647, Facultat de Física, 08028 Barcelona, Catalonia Martin Luc Rosinberg and Gilles Tarjus Laboratoire de Physique Théorique des Liquides, Université Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris, France sReceived 12 November 2004; published 29 April 2005d We study the nonequilibrium behavior of the three-dimensional Gaussian random-field Ising model at T = 0 in the presence of a uniform external field using a two-spin-flip dynamics. The deterministic, historydependent evolution of the system is compared with the one obtained with the standard one-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 sin number and sized 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. DOI: 10.1103/PhysRevB.71.134424 PACS numberssd: 75.60.Ej, 05.70.Jk, 75.40.Mg, 75.50.Lk I. INTRODUCTION The nonequilibrium random field Ising model sRFIMd 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 sor metastable statesd separated by sizeable barriers, which makes thermally activated processes essentially irrelevant at low enough temperature sthe lifetime of metastable states may then be considered as infinited. As a consequence, these systems remain far from equilibrium on the experimental time scales seven when the driving rate goes to zerod and their response to the external force is made of a series of jumps savalanchesd between neighboring metastable states. This type of behavior is very well modeled by the ferromagnetic RFIM with a zero-temperature single-spin-flip 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 sor “wiping out” effectd 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 weak-disorder one where it has a discontinuous jump. Such a transition has been observed in thin Co/ CoO films2 and Cu-Al-Mn alloys,3 and it has been recently suggested4 that it may also be associated to the change in the adsorption behavior of 4He 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 1098-0121/2005/71s13d/134424s8d/$23.00 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 nonequilibrium 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 and 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 sin the literature on the RFIM, it is implicitly taken for granted that this should not matterd. The single-spin-flip dynamics, however, is not the unique sand may be not the bestd 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,13 that the nonequilibrium 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 nonzero coercive field in the irreversible evolution.14 Since the ground state is stable with respect to the flip of an arbitrary sfinited 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 two-spin-flip dynamics. We compare the results with those obtained with the standard one-spinflip dynamics, in particular those concerning the number and size of the avalanches. We first show in Sec. II that one can indeed define a two-spin-flip algorithm that yields a deterministic evolution of the system with the external field. In particular, the dynamics satisfies the “abelian” property 134424-1 ©2005 The American Physical Society PHYSICAL REVIEW B 71, 134424 s2005d VIVES, ROSINBERG, AND TARJUS which guarantees that the same two-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 return-point memory. In Sec. III we present the results of our numerical simulations on three-dimensional s3Dd lattices and compare them to the behavior of the same system with a one-spin-flip dynamics. The hysteresis loops are significantly reduced when allowing two-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 Sec. 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 one-spin-flip dynamics. The main conclusions of the study are reported in Sec. V. II. MODEL AND DYNAMICS The model is defined on a cubic lattice of linear size L with periodic boundary conditions. On each site si = 1,…,N = L3d there is an Ising spin variable sSi = ± 1d. The Hamiltonian is H=−J o SiS j − oi hiSi − H oi Si , s1d ki,jl where the first sum extends over all distinct pairs of nearestneighbor snnd , H is the external applied field, and hi are quenched random fields drawn independently from a Gaussian distribution with zero mean and standard deviation s. 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 sin the adiabatic limit corresponding to a vanishingly small rate of change of the external fieldd. For this purpose, the Hamiltonian must be supplemented by some dynamical rules. The standard one-spin-flip dynamics used in previous studies consists in minimizing the energy of each spin Hi = − Si f i , s2d where fi = J S j + hi + H o jsid FIG. 1. Stability diagrams showing the state with minimum local energy according to the values of the fields created by the neighborhood sdefined in the textd: sad corresponds to a single spin i and sbd to a pair of neighboring spins i , j. come 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 se.g., by increasing i from 1 to Nd, 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 “nopassing” and abelian properties of this dynamics.1,16 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 two-spin-flip dynamics is rather straightforward. By definition, two-spin-flip stable states are spin configurations whose energy sdefined solely by the Hamiltoniand cannot be lowered by the flip of one or two spins sclearly, new features are only introduced when these two spins are nnd. The slocald energy to be minimized is thus the one associated with a pair ij of nn spins Hij = − Si f 8i − S j f 8j − JSiS j , s3d s4d where is the net field at site i fthe summation in Eq. s3d is over the z nn of ig. From the above expression, it is clear that the minimization of Hi is obtained by aligning each spin with its local field, Si = signsf id, as represented schematically in Fig. 1sad. This provides a stability criterion for any state with respect to this one-spin-flip dynamics. This dynamical rule may be implemented by an algorithm that propagates one avalanche at a time.15 Starting from a stable configuration, one increases sor decreasesd the external field until the local field f i at some site i becomes zero sthis corresponds to the vanishing of the local minimum in which the system was trappedd. The spin Si swhich is uniquely defined because the distribution of the random fields is continuousd is then flipped, which in turn may cause neighboring spins to be- f i8 = J o Sk + hi + H s5d ksidÞj is the field experienced by Si without the influence of the neighbor S j fthe summation in Eq. s5d is over the nn of i excepting jg. One can then think of the dynamics as made of single-spin flips and “irreducibly cooperative” twospin flips. As pictured in Fig. 1sbd, a single-spin flip occurs whenever the net field on Si, f i = f i8 + JS j, or on S j, f j = f 8j + JSi, changes sign. This corresponds to the changes ↓ ↓ ↔ ↑ ↓, ↓ ↓ ↔ ↓ ↑, ↑ ↑ ↔ ↑ ↓, ↑ ↑ ↔ ↓ ↑ in the diagram. An irreducibly cooperative two-spin flip involves a nn pair of spins with the same sign s↑↑ or ↓↓d that cannot flip individually si.e., without the simultaneous flip of the neighbord. This 134424-2 PHYSICAL REVIEW B 71, 134424 s2005d HYSTERESIS AND AVALANCHES IN THE T = 0… occurs whenever the net field on the pair of aligned spins, f ij = f i8 + f 8j , changes sign in the region of the diagram where −J ø f 8i ø J and −J ø f 8j ø J. These two last conditions come from the fact that the spins cannot individually flip shence neither f i nor f j changes sign before the cooperative flip of the paird 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 one-spin-flip dynamics. Starting from a two-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. 1sbd leaves the region where it was originally sassociated with ↑ ↑ , ↑ ↓ , ↓ ↑, or ↓↓d and, depending on the border which is first attained, only a single spin flips or the two spins flip simultaneously.17 It is not hard to show that this dynamics obeys the same properties as the one-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 scorresponding to all the distinct nn pairsd in the diagram of Fig. 1sbd. As the external field is monotonously increased sresp. decreasedd, the local fields can only increase sresp. decreased, the spins can only flip up sresp. downd, and the points can only move up and right sresp. down and leftd 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. sWe have also performed numerical tests in many situations and found no violations of these properties.d The evolution of a system with size L = 30 and s = 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 Si / N sstrictly speaking, e is the enthalpyd. The external field H is varied from a very large initial value where all spins are up to the final value H = −2 shere and after, J is taken as the energy unitd. 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 swith sequential updatingd 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 s = 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. III. NUMERICAL SIMULATIONS We now present the results of numerical simulations performed on 3-dimensional cubic lattices using either the one- FIG. 2. sColor onlined Evolution of the enthalpy and the magnetization sper spind for a system of size L = 30 and s = 2.5 as the external field is changed from H = + ` to H = −2. The trajectories corresponding to sequential and parallel updating algorithms are also compared with the adiabatic one. Points indicate the intermediate unstable states. The inset shows a closeup of the region around the end point. Notice that the end point is the same in all cases. spin-flip or two-spin-flip dynamics with sequential updating. Averaged quantities were obtained with statistics over 103 – 105 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 sespecially the “spanning” avalanchesd, it is more important to perform averages over many disorder realizations than to simulate very large system sizes. A. Hysteresis loops Figure 4 shows the hysteresis loops obtained in a single sample for two different values of the disorder s. For com- FIG. 3. sColor onlined Major hysteresis loop and internal loop s−0.85ø H ø 0.8d obtained with the two-spin-flip dynamics for a system of size L = 30 and s = 2.5. The inset shows the details around the point at H = −0.85 revealing the property of return-point memory. 134424-3 PHYSICAL REVIEW B 71, 134424 s2005d VIVES, ROSINBERG, AND TARJUS FIG. 4. sColor onlined Magnetization curves obtained with the one-spin and two-spin-flip dynamics in a sample of size L = 30 for sad s = 2 and sbd s = 3: dotted-dashed lines correspond to the onespin-flip dynamics and dashed lines to the two-spin-flip dynamics. In addition, the ground-state magnetization is shown as continuous lines. parison, we also display the magnetization curves obtained with the one-spin-flip dynamics and with the algorithm of Ref. 18 which gives the exact ground-state sequilibriumd magnetization. The corresponding behavior of the enthalpy per spin along the ascending branches of the loops is reported in Fig. 5 snote in passing that the ground-state enthalpy does not show any discontinuity as the external field is varied18d. As could be expected, the main effect of the new dynamics is to reduce the size of the hysteresis loops. Specifically, the coercivity si.e., the magnitude of the external FIG. 5. sColor onlined Enthalpy per spin as a function of H corresponding to the loops in Fig. 5 sfor clarity, only the ascending branches are shownd. The continuous lines represent the groundstate behavior. field for which the magnetization is equal to zerod 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-of-equilibrium disorder-induced phase transition under the two-spin-flip dynamics, with a critical value of s at which the discontinuity appears in the thermodynamic limit. It is worth pointing out that the new dynamical rules allow the system to effectively overcome energy barriers of magnitude up to DE = 2J. Indeed, the difference between the two dynamics shows up when a pair of nn 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 one-spin-flip paths, the system has now been able to bypass the higher-energy states, either ↑↓ or ↓↑. By using Eq. s4d, it is easy to show that the relevant barrier height associated with this process is at most 2J. Since cooperative flips of more than two spins do not occur with the chosen dynamics and the system’s trajectory otherwise go through states of decreasing energy, one concludes that DE = 2J is the maximum barrier height that the system may overcome when passing from the one-spin-flip to the two-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 savalanchesd 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 one, two, or three spatial directions sindicated in the following by the index ad. For each individual avalanche, this is a property that can be easily detected during the simulation. Avalanches are thus classified as being nonspanning sa = nsd, one-dimensional s1Dd spanning sa = 1d, two-dimensional s2Dd spanning sa = 2d, or three-dimensional s3Dd spanning sa = 3d. Figure 6 shows the number of 1D, 2D, and 3D spanning avalanches recorded along the descending branch of the hysteresis loops as a function of s. 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 toward larger values of s when the two-spin-flip dynamics is used. The same shift is also found for all studied system s1d sizes. This is a first indication that ss2d c . sc , as will be confirmed by the finite-size scaling analysis presented in the next section. In the case of the one-spin-flip dynamics, a detailed analysis was performed in Refs. 8 and 9, revealing the scenario that occurs in the thermodynamic limit and that is already suggested by the data shown in Fig. 6: when 134424-4 PHYSICAL REVIEW B 71, 134424 s2005d HYSTERESIS AND AVALANCHES IN THE T = 0… FIG. 6. sColor onlined Average numbers of 1D, 2D, and 3Dspanning avalanches as a function of s in a system of size L = 24. Continuous and dashed lines are guides for the eye and correspond to the one-spin-flip and two-spin-flip dynamics, respectively. L → `, N1ssd and N2ssd are expected to display a d singularity at sc and N3ssd a steplike 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 sthere is only one, compact, subcritical avalanche for s , scd whereas the latter only exist at sc shence, the additional d singularity at the edge of the step function whose signature is already visible in Fig. 6d. 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 be- cause of the complexity of the two-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 3D-spanning avalanches. Figure 7 shows the snormalizedd avalanche size distributions Dass ; s , Ld obtained along one branch of the hysteresis loop for three different values of s sfor clarity, the results obtained with the one-spin-flip dynamics are represented by continuous linesd. Surprisingly, one can see that the distribution of non-spanning avalanches in Figs. 7sad–7scd is almost unaltered by the change in the dynamics. The only small difference sbarely visible on the figured induced by the twospin-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 teff < 2.0 in the thermodynamic limit.9 The size distributions D1 and D2 of the 1D and 2Dspanning avalanches shown in Figs. 7sdd–7sgd also appear to be identical with the two dynamics, at least within statistical error bars. sNote that these avalanches do not exist for s = 1.6 because this value is much lower than sc for both dynamics.d The only visible differences between the two dynamics occur in D3, the size distribution of the 3D-spanning avalanches. Specifically, for s = 1.60 and 2.25 si.e., below sc and very close to sc, respectivelyd, the large 3D-spanning avalanches tend to be shifted to even larger sizes. According to Refs. 8 and 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 indis1d cation that ss2d c . sc . In contrast, for s = 2.80 swhich is clearly above scd, the distribution D3 is not affected by the dynamics fFig. 7sjdg: in this case, one expects to detect only critical 3D-spanning avalanches in a finite system. FIG. 7. sColor onlined Avalanche size distributions for s = 1.60 sa,hd, 2.25 sb,d,f,id, and 2.80 sc,e,g,jd. The open symbols and the continuous lines correspond to the two-spin-flip and one-spin-flip dynamics, respectively. The first row shows the distribution Dns of non-spanning avalanches and the other rows show the distributions D1 , D2, and D3 of the 1D, 2D, and 3D-spanning avalanches, respectively. All data have been obtained in a system of size L = 24. In sbd, the dashed line indicates the expected power-law behavior of Dns at criticality with exponent teff = 2.0. 134424-5 PHYSICAL REVIEW B 71, 134424 s2005d VIVES, ROSINBERG, AND TARJUS FIG. 8. sColor onlined Number of avalanches spanning in 1D sad and 2D sbd obtained with the two-spin-flip dynamics for different system sizes. The corresponding scaling plots in scd and sdd have been obtained according to Eqs. s6d–s8d s2d with sc = 2.25, A = −0.2, n = 1.2, and u = 0.1. IV. FINITE-SIZE SCALING ANALYSIS ss2d c and As already mentioned, a precise determination of of some of the critical exponents requires a detailed finitesize scaling analysis which, unfortunately, is not possible with the present algorithm. What can be done, however, is to check whether the set of exponents found with the one-spinflip dynamics can also be used to scale the two-spin-flip data. We first consider N1 and N2, the numbers of 1D and 2Dspanning avalanches, which are quantities that have a simple scaling behavior. Following previous work,8,9 we assume the forms N1ss,Ld = LuÑ1su2L1/nd s6d N2ss,Ld = LuÑ2su2L1/nd, s7d thus indicate that the two-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 kH3ssdl at which the 3D-spanning avalanches occur. This quantity was studied with the one-spinflip dynamics in Ref. 9, which allowed to map out the firstorder line scorresponding to the macroscopic jump in the magnetizationd in the diagram H − s. The dependence of kH3ssdl with the system size is shown in Fig. 9sad. Note that the transition line extends above sc in a finite system. However, the end point, beyond which no 3D-spanning avalanches are found swith the present sampling of disorder realizationsd, becomes closer and closer to the critical point sHc , scd as L is increased. where u and n are critical exponents, Ñ1 and Ñ2 are scaling functions, and u2 is a scaling variable that measures the distance to the critical point. It is defined as u2 = s − ss2d c ss2d c +A F s − ss2d c ss2d c G 2 , s8d 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 n = 1.2, u = 0.1, and A = −0.2, values that we keep here. The only free parameter is thus ss2d c . As shown in Fig. 8, a very good collapse of the new data can be obtained with ss2d 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 ss1d c = 2.21 obtained with the one-spin-flip dynamics8 can be discarded. We have also checked that this conclusion is not modified when the nonuniversal parameter A is allowed to vary. These scaling plots FIG. 9. sColor onlined sad Average field kH3ss , Ldl at which the 3D-spanning avalanches occur as a function of s for different syss2d tem sizes. sbd Scaling plot of the data below sc according to Eq. s2d s9d with Hc = −0.885, B8 = 0.25, and m = 1.5. 134424-6 PHYSICAL REVIEW B 71, 134424 s2005d HYSTERESIS AND AVALANCHES IN THE T = 0… similar analysis for Ñ2ss , Ld shows the same agreement. This is again a strong indication for universality, that goes beyond the equality of the critical exponents.19 V. CONCLUSION FIG. 10. sColor onlined Comparison of the scaling plots of the number of 1D-spanning avalanches obtained with the one-spin-flip ssymbols with error barsd and two-spin-flip ssymbols without error barsd dynamics. The two collapses only differ by the value of the s1d s2d critical disorder fsc = 2.21 and sc = 2.25g. The continuous line is the Gaussian fit proposed in Ref. 8. According to Ref. 9, this set of curves should scale as F 1 − B8 kH3lss,Ld = Hs2d c s − ss2d c ss2d c G − L−1/mĥ3su2L1/nd, s9d where Hs2d c is the critical field, B8 is a nonuniversal tilting constant, m is a critical exponent, and ĥ3 is the corresponding scaling function. Strictly speaking, the scaling should be done separately for the average fields kH3−l and kH3cl 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 ssmalld effect only in the region u2L1/n . 1 where the two kinds of avalanches coexist.9 In Ref. 9, the best collapse of the one-spin-flip data was obtained with m = 1.5, B8 = 0.25, and Hs1d c = −1.425. Here, we keep the same values for m and B8, set ss2d c = 2.25, and consider the critical field as the only free parameter. As shown in Fig. 9sbd, a very good collapse can be obtained with Hs2d c = −0.885 ffor clarity, we only present the collapse for s , ss2d c g. 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 sin absolute valued is in line with the decrease in coercivity illustrated by Fig. 4. Finally, it is interesting to study the influence of the dynamics on the finite-size scaling functions. In Fig. 10, we compare the scaling collapses of the number of 1D-spanning avalanches, Ñ1ss , Ld, obtained with the two dynamics sfor the two-spin-flip dynamics, this is the same curve as in Fig. 8d. As can be seen, the agreement between the two curves is quite remarkable. Even the deviations saround u2L1/n < 0d from the Gaussian fit proposed in Ref. 8 are the same. A We have shown in this paper that the nonequilibrium behavior of the 3D RFIM at zero temperature is not qualitatively altered when going from the standard one-spin-flip metastable dynamics to a two-spin-flip metastable dynamics. The coercivity sthat is essentially the width of the hysteresis loop associated with the evolution of the magnetization with the driving fieldd is significantly reduced by allowing twospin flips, but the main features of the hysteretic behavior, including the presence of a disorder-induced nonequilibrium 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 sexponents and scaling functionsd obtained with the two-spin-flip dynamics is in the same universality class as that obtained with the one-spin-flip dynamics. Changing the dynamical rules used to study the evolution of the model, as done here, helps address several important questions. The first one is the robustness of the hysteretic scenario provided by the T = 0 RFIM with the standard onespin-flip dynamics. As discussed in the 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 two-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 one-spinflip dynamics.1,6–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,13 sAdditional, but less conclusive evidence is provided by exact, but mean-field-like results on the Bethe lattice13 and by perturbation theory near the upper critical dimension d = 6.6d Our present findings suggest that a whole series of T = 0 metastable dynamics involving k-spin flips lead to the same nonequilibrium 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 134424-7 PHYSICAL REVIEW B 71, 134424 s2005d VIVES, ROSINBERG, AND TARJUS of spins. Based on the numerical closeness of the critical exponents and scaling functions, on the fact that the critical s2d s1d eq disorder strength satisfies seq c . sc . sc fwhereas Hc = 0 s2d s1d , Hc , Hc g, 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 In the opposite case, one would expect to see for large enough k a crossover between a behavior dominated by the equilibrium fixed point and the asymptotic behavior controlled by the nonequilibrium fixed point. No such crossover is seen in the present study ssee, e.g., the distribution of avalanches in Fig. *Electronic address: eduard@ecm.ub.es 1 J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. 70, 3347 s1993d. 2 A. Berger, A. Inomata, J. S. Jiang, J. E. Pearson, and S. D. Bader, Phys. Rev. Lett. 85, 4176 s2000d. 3 J. Marcos, E. Vives, Ll. Mañosa, M. Acet, E. Duman, M. Morin, V. Novak, and A. Planes, Phys. Rev. B 67, 224406 s2003d. 4 F. Detcheverry, E. Kierlik, M. L. Rosinberg, and G. Tarjus, Phys. Rev. E 68, 061504 s2003d; Langmuir 20, 8006 s2004d. 5 D. J. Tulimieri, J. Yoon, and M. H. W. Chan, Phys. Rev. Lett. 82, 121 s1999d. 6 K. Dahmen and J. P. Sethna, Phys. Rev. B 53, 14 872 s1996d. 7 O. Perković, K. A. Dahmen, and J. P. Sethna, Phys. Rev. B 59, 6106 s1999d. 8 F. J. Pérez-Reche and E. Vives, Phys. Rev. B 67, 134421 s2003d. 9 F. J. Pérez-Reche and E. Vives, Phys. Rev. B 70, 214422 s2004d. 10 G. Durin and S. Zapperi, in The Science of Hysteresis, edited by G. Bertotti and I. Mayergoyz sElsevier, New York, 2005d. 11 J. P. Sethna, K. A. Dahmen, and O. Perković, in The Science of Hysteresis, edited by G. Bertotti and I. Mayergoyz sElsevier, New York, in pressd. 12 A. Maritan, M. Cieplak, M. R. Swift, and J. R. Banavar, Phys. 7d, but the values of k considered here, k = 1 and 2, are of course too small to detect such effect, if present. ACKNOWLEDGMENTS We acknowledge fruitful discussions with Francisco-José Perez-Reche and Xavier Illa. E.V. also acknowledges the hospitality of the Laboratoire de Physique Théorique des Liquides sUPMC, Parisd during his stay as an invited professor in July 2004. This work has received financial support from CICyT sSpaind, Project No. MAT2004-01291, and CIRIT sCataloniad, Project No. 2001SGR00066. The Laboratoire de Physique Théorique des Liquides is the UMR 7600 of the CNRS. Rev. Lett. 72, 946 s1994d. Colaiori, M. J. Alava, G. Durin, A. Magni, and S. Zapperi, Phys. Rev. Lett. 92, 257203 s2004d; M. J. Alava, V. Basso, F. Colaiori, L. Dante, G. Durin, A. Magni, and S. Zapperi, Phys. Rev. B 71, 064423 s2005d. 14 This objection, however, may be circumvented by considering the remanent magnetization obtained in zero external field after a demagnetization procedure as the order parameter of the transition ssee Ref. 13d. 15 M. C. Kuntz, O. Perković, K. A. Dahmen, B. W. Roberts, and J. P. Sethna, Comput. Sci. Eng. 1, 73 s1999d. 16 D. Dhar, P. Shukla, and J. P. Sethna, J. Phys. A 30, 5259 s1997d. 17 Unfortunately, because of the additional complexity introduced by the cooperative flip of nn pairs, one cannot use anymore the so-called sorted-list algorithm proposed in Ref. 15. This in turn implies that very large system sizes cannot be simulated. 18 C. Frontera, J. Goicoechea, J. Ortin, and E. Vives, J. Comput. Phys. 160, 117 s2000d. 19 Note also that for the scaling collapses it has been possible to use the same values of the a priori nonuniversal constants A and B8 in Eqs. s8d and s9d. 13 F. 134424-8