- Email: [email protected]
“Avalanches” in the ground state of the 3D Gaussian random field Ising model driven by an external field
Computer Physics Communications 147 (2002) 455–458 www.elsevier.com/locate/cpc
“Avalanches” in the ground state of the 3D Gaussian random field Ising model driven by an external field Carlos Frontera a , Eduard Vives b,∗ a Institut de Ciència de Materials de Barcelona, Consell Superior d’Investigacions Científiques Campus UAB, 08193 Bellaterra,
Catalonia, Spain b Departament d’Estructura i Constituents de la Matèria, Facultat de Física, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona,
Catalonia, Spain
Abstract We present a numerical study of the exact ground states of the 3D Gaussian random field Ising model (G-RFIM) with an applied external field B. We combine a max-flow min-cut algorithm with an optimal procedure for determining all the ground states when B is swept from −∞ to ∞. The magnetization of finite lattices (L3 ) is studied as a function of the degree of disorder in the system σ (standard deviation of the Gaussian random fields). The magnetization evolves as a sequence of jumps or “avalanches” with a certain size s. The statistical distribution p(s) becomes a power law p(s) ∼ s −τ for a certain degree of disorder σc (L). The extrapolation of the results to L → ∞ renders σc 2.4 ± 0.1 and τ 1.70 ± 0.07. 2002 Elsevier Science B.V. All rights reserved. PACS: 05.50.+q; 75.10.Nr; 64.60.Fr; 75.60.Ej Keywords: 3D Gaussian random field Ising model; Exact ground states; Avalanches
1. Introduction Avalanche phenomena usually occur in slowly driven systems that exhibit a sudden change of a magnitude with random properties. Their appearance, size and duration are difficult to predict. In many cases power-law distributions of such magnitudes are encountered which indicate a certain degree of criticality, i.e. the absence of characteristic spatial and temporal scales. The most usual physical framework for the description of avalanche phenomena is the so-called selforganized criticality (SOC) theory [1]. The theory ap* Corresponding author.
E-mail address: [email protected] (E. Vives).
plies to non-linear dissipative extended systems (with spatial and temporal degrees of freedom), which are kept in a non-equilibrium steady state by the maintenance of an external input of energy which compensates dissipation. Under such conditions the systems naturally evolve to a state exhibiting avalanches. The theory has been applied not only to the study of many natural phenomena such as piles of granular materials, earthquakes, river networks, biological systems, etc. . . , but also to the study of economics and social behaviour [1]. Nevertheless, SOC is not the only theory that explains avalanches [2]. In many cases avalanches are associated with the (non-equilibrium) metastable evolution through a first-order phase transition. For such
0010-4655/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 2 ) 0 0 3 2 8 - 4
456
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
cases many models have been developed including as main ingredients: (i) a first-order phase transition when driven by an external field, as occurs in the Ising model; (ii) a certain amount of quenched disorder, for instance local random fields; (iii) irrelevance of thermal fluctuations and (iv) a metastable dynamics. The original studies of the T = 0 G-RFIM [3] were followed by many other models [4,5]. When B is swept, the system evolves by avalanches joining metastable states until the whole system has been transformed. When the external field is reversed the system exhibits hysteresis. For a certain degree of disorder the behaviour is critical, i.e. the avalanche sizes and durations become power-law-distributed. This second framework has been applied to the study of avalanches in many first-order phase transitions [6]: Barkhausen noise in magnets, acoustic emission in martensites, superconductivity, etc. The two above scenarios correspond to out-ofequilibrium phenomena. More recently we have suggested that “avalanche”-like phenomena could also occur in equilibrium systems [7]. The models proposed are those presented in the above paragraph but under equilibrium conditions at T = 0, i.e. the system is driven by the external field B through the sequence of exact ground states. Numerical studies for finite lattices (L2 ) of the 2D G-RFIM model have revealed that the magnetization evolves as a sequence of jumps. The amount of disorder in the systems is, in this case, controlled by the standard deviation σ of the Gaussian distribution of random fields. For a certain critical amount of disorder σc (L), the magnetization jumps become power-law-distributed. Nevertheless, this third possibility for the existence of “avalanche”-like phenomena remains unclear since the extrapolation of σc (L) to the thermodynamic limit is difficult. A number of authors [8] have shown that in this 2D case the ferromagnetic phase is unstable for all σ > 0. In order to circumvent this problem, in this paper we have studied the 3D G-RFIM, for which the existence of a ferromagnetic ground state is more generally accepted.
2. Model and numerical details The model is defined on a cubic lattice of size N = L3 . At each site we define a spin variable Si (i =
1, . . . , N ) taking values +1 or −1. The Hamiltonian reads: H =−
n.n. i,j
Si Sj −
N
Si hi − B
i=1
N
Si .
(1)
i=1
The first term, which extends over all nearest-neighbor pairs, stands for a ferromagnetic interaction. The second term stands for the interaction with the local quenched random fields hi , which are independent and Gaussian-distributed with zero mean hi = 0 and variance h2i = σ 2 . The last term accounts for the interaction with the external field B. Efficient algorithms for the study of the exact ground states of (1) have been available for several years. We have used a max-flow min-cut algorithm that for a certain fixed set of the random fields {hi }, and for each value of B in the whole range −∞ < B < ∞ finds the configuration {Si } that minimizes H . The algorithm used [9] determines not only the sequence of ground-state configurations, but also the exact values of B where the magnetization jumps s take place. We have systematically studied systems of sizes L = 8, 12, 16, 20, 24, 28 and 32. In order to compute the histograms, a large number of realizations of the random fields for each value of σ have been analyzed. For instance, for L 28 typical histograms are obtained from several thousands of configurations, whereas for L = 32 only several hundreds are considered. 3. Results Fig. 1 shows the log–log plot of the normalized distribution of “avalanche” sizes p(s) for systems with L = 28 and different values of σ ranging from 2.10 to 2.70. The values of s are discrete ranging from 1 to 283 . Note that data has not been binned so that the discretization of the data in the vertical direction is related to the statistical confidence which, in some cases, is of the order of 10−5 . As can be seen, for the largest σ values the distribution of “avalanche” sizes extends over the small s region whereas for lower values a peak corresponding to very large “avalanches” (with a size of the order of L3 ) appears. In order to model such histograms we have performed maximum-likelihood fits of the two freeparameter expression: p(s) = A(λ, τ, L)e−λs s −τ ,
(2)
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
457
Fig. 3. Numerical estimations of σc (L) versus 1/L. The lines show different fits of Eq. (3) as indicated by the legend. Fig. 1. Distribution of “avalanche” sizes in log–log scale. Data corresponds to a system with L = 28 and different values of σ . The lines correspond to the fits of Eq. (2). The histograms have been consecutively shifted 2 decades upwards in order to clarify the picture.
Fig. 2. Parameters τ (a) and λ (b) obtained from the fit of the expression 2 to the numerical data, for system sizes ranging from L = 8 to L = 32, as indicated by the legend. The dotted lines indicate the values of σ for which λ = 0.
3 −λs −τ where A(λ, τ, L) = 1/ L s is the normals=1 e ization constant. The maximum-likelihood procedure, contrarily to the χ 2 methods, is independent of any binning process and of the type of representation (linear, linear–log or log–log) of the data. The results of the fits are also shown in Fig. 1. In general the fits are good in the small s region but are unable to quantitatively reproduce the peak in the large s region for small σ ’s. The resulting fitted values of the parameters λ and τ are shown in Fig. 2. As can be seen, for
any value of L, there exists a certain σc (L) for which λ = 0 and, therefore, the distribution of “avalanches” becomes power-law. For L = 28, for instance, this occurs for σc (L = 28) 2.6 as can also be seen from Fig. 1 when the peak corresponding to large s values disappears. Estimations of σc (L) can be obtained by interpolation to the λ(σ ) data close to the crossing point, using a linear fit. The obtained values of σc (L) in front of 1/L are shown in Fig. 3. At first sight, it becomes clear that in the thermodynamic limit a critical point with power-law distribution of “avalanches” exists. Before performing a more detailed analysis of the extrapolation to L → ∞ it is convenient to review some of the ground-state properties of the 3D G-RFIM which are not fully understood. In particular, the exact position σc of the T = 0 critical point and the values of the critical exponents characterizing it have been under discussion for a long time. Most of the previous studies of the 3D G-RFIM have been performed at B = 0. After the pioneering work of Ogielski [10] an intense debate focused on whether the transition is discontinuous [11] or continuous [12] with a very small value of the critical exponent β. A very recent and deep analysis [13] has tilted the balance towards the existence of a standard critical point at σc = 2.270 ± 0.004 with β = 0.017 ± 0.005 and ν = 1.37 ± 0.09. To our knowledge only a single study with an applied external field has been performed. MC simulations at finite temperatures [14] with L up to 24 are not totally conclusive, but the authors propose that, instead of a critical point at (Tc (σ ), B = 0), two twin critical points could appear at (Tc (σ ), ±Bc (σ )).
458
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
These are joined by first-order phase transition lines to a triple point at (T3 (σ ), B = 0). The extrapolation to T → 0 gives σ3 = 2.35. According to the standard finite size scaling ansatz σc (L) should behave like: σc (L) = σc + AL−1/ν .
(3)
In Fig. 3 we have compared our data with different fits of Eq. (3): (i) the best linear fit (ν = 1), (ii) the best fit imposing ν = 1.37 [13], (iii) the best three freeparameters fit, and (iv) the best fit imposing ν = 1.37 and σc = 2.27 [13] (restricted to the L 24 data). Although it is difficult to distinguish the quality of the fits (i) and (iii), the fits (ii) and (iv) are qualitatively worse. Moreover, the extrapolated values of σc for the fits (i), (ii) and (iii) are clearly above σc = 2.27 [13]. Our best (and very conservative) estimation is σc = 2.4 ± 0.1. Therefore, we conclude that the “critical” point characterized by a power-law distribution of magnetization jumps when B is swept from −∞ to ∞ occurs at a value above the critical point found from the studies at B = 0. Moreover, from Fig. 2, we can also obtain an estimation of the exponent τ . Although λ exhibits a clear finite size dependence, the value of the exponent τ (at the λ = 0 point) is almost size independent, as indicated by the dashed lines. Our best estimation is τ = 1.70 ± 0.07. Finally, it is also interesting to compare our data with the data from the studies of the same 3D G-RFIM using the metastable dynamics [5] that have found a critical point at σc = 2.16 ± 0.03 with τ = 2.03 ± 0.03. Therefore, our present results do not allow any connection to be established between both critical points.
Acknowledgements The authors acknowledge financial support from CICyT project numbers: MAT1999-0984 (C.F.) and
MAT2001-003251 (E.V.). Part of the data has been obtained by using a free distribution of the maxflow min-cut algorithm from LEDA Research Group (Germany). This research has partially been developed using computational facilities of CESCA (Centre de Supercomputació de Catalunya). References [1] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381; How Nature Works, Oxford. Univ. Press, Oxford, 1997. [2] D. Sornette, J. Phys. I (France) 4 (1993) 209. [3] J.P. Sethna, K. Dahmen, S. Kartha, J.A. Krumhansl, B.W. Roberts, J.D. Shore, Phys. Rev. Lett. 70 (1993) 3347. [4] K. Dahmen, J.P. Sethna, Phys. Rev. Lett. 71 (1993) 3222; E. Vives, A. Planes, Phys. Rev. B 50 (1994) 3839; E. Vives, J. Goicoechea, J. Ortín, A. Planes, Phys. Rev. E 52 (1995) R5; B. Tadi´c, Phys. Rev. Lett. 77 (1996) 3843; J. Kushauer, van Bentum, W. Kleemann, D. Bertand, Phys. Rev. B 53 (1996) 11 647. [5] O. Perkovic, K.A. Dahmen, J.P. Sethna, Phys. Rev. B 59 (1999) 6106. [6] Ll. Carrillo, Ll. Mañosa, J. Ortín, A. Planes, E. Vives, Phys. Rev. Lett. 81 (1998) 1889, and references therein. [7] C. Frontera, E. Vives, Phys. Rev. E 62 (2000) 7470. [8] E.T. Seppälä, M.J. Alava, P.M. Duxbury, Phys. Rev. E. 63 (2001) 066110, and references therein. [9] C. Frontera, J. Goicoechea, J. Ortín, E. Vives, J. Comput. Phys. 160 (2000) 117. [10] A.T. Ogielski, Phys. Rev. Lett. 57 (1986) 1251. [11] J.-C.Anglès d’Auriac, N. Sourlas, Europhys. Lett. 39 (1997) 473; N. Sourlas, Comp. Phys. Comm. 121 (1999) 183. [12] M.E.J. Newman, G.T. Barkema, Phys. Rev. E 53 (1996) 393; U. Nowak, K.D. Usadel, J. Esser, Physica A 250 (1998) 1; A.K. Hartmann, U. Nowak, Eur. Phys. J. B 7 (1999) 105. [13] A.A. Middleton, D.S. Fisher, Phys. Rev. B 65 (2002) 134411. [14] J. Machta, M.E.J. Newman, L.B. Chayes, Phys. Rev. E 62 (2000) 8782.
“Avalanches” in the ground state of the 3D Gaussian random field Ising model driven by an external field Carlos Frontera a , Eduard Vives b,∗ a Institut de Ciència de Materials de Barcelona, Consell Superior d’Investigacions Científiques Campus UAB, 08193 Bellaterra,
Catalonia, Spain b Departament d’Estructura i Constituents de la Matèria, Facultat de Física, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona,
Catalonia, Spain
Abstract We present a numerical study of the exact ground states of the 3D Gaussian random field Ising model (G-RFIM) with an applied external field B. We combine a max-flow min-cut algorithm with an optimal procedure for determining all the ground states when B is swept from −∞ to ∞. The magnetization of finite lattices (L3 ) is studied as a function of the degree of disorder in the system σ (standard deviation of the Gaussian random fields). The magnetization evolves as a sequence of jumps or “avalanches” with a certain size s. The statistical distribution p(s) becomes a power law p(s) ∼ s −τ for a certain degree of disorder σc (L). The extrapolation of the results to L → ∞ renders σc 2.4 ± 0.1 and τ 1.70 ± 0.07. 2002 Elsevier Science B.V. All rights reserved. PACS: 05.50.+q; 75.10.Nr; 64.60.Fr; 75.60.Ej Keywords: 3D Gaussian random field Ising model; Exact ground states; Avalanches
1. Introduction Avalanche phenomena usually occur in slowly driven systems that exhibit a sudden change of a magnitude with random properties. Their appearance, size and duration are difficult to predict. In many cases power-law distributions of such magnitudes are encountered which indicate a certain degree of criticality, i.e. the absence of characteristic spatial and temporal scales. The most usual physical framework for the description of avalanche phenomena is the so-called selforganized criticality (SOC) theory [1]. The theory ap* Corresponding author.
E-mail address: [email protected] (E. Vives).
plies to non-linear dissipative extended systems (with spatial and temporal degrees of freedom), which are kept in a non-equilibrium steady state by the maintenance of an external input of energy which compensates dissipation. Under such conditions the systems naturally evolve to a state exhibiting avalanches. The theory has been applied not only to the study of many natural phenomena such as piles of granular materials, earthquakes, river networks, biological systems, etc. . . , but also to the study of economics and social behaviour [1]. Nevertheless, SOC is not the only theory that explains avalanches [2]. In many cases avalanches are associated with the (non-equilibrium) metastable evolution through a first-order phase transition. For such
0010-4655/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 2 ) 0 0 3 2 8 - 4
456
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
cases many models have been developed including as main ingredients: (i) a first-order phase transition when driven by an external field, as occurs in the Ising model; (ii) a certain amount of quenched disorder, for instance local random fields; (iii) irrelevance of thermal fluctuations and (iv) a metastable dynamics. The original studies of the T = 0 G-RFIM [3] were followed by many other models [4,5]. When B is swept, the system evolves by avalanches joining metastable states until the whole system has been transformed. When the external field is reversed the system exhibits hysteresis. For a certain degree of disorder the behaviour is critical, i.e. the avalanche sizes and durations become power-law-distributed. This second framework has been applied to the study of avalanches in many first-order phase transitions [6]: Barkhausen noise in magnets, acoustic emission in martensites, superconductivity, etc. The two above scenarios correspond to out-ofequilibrium phenomena. More recently we have suggested that “avalanche”-like phenomena could also occur in equilibrium systems [7]. The models proposed are those presented in the above paragraph but under equilibrium conditions at T = 0, i.e. the system is driven by the external field B through the sequence of exact ground states. Numerical studies for finite lattices (L2 ) of the 2D G-RFIM model have revealed that the magnetization evolves as a sequence of jumps. The amount of disorder in the systems is, in this case, controlled by the standard deviation σ of the Gaussian distribution of random fields. For a certain critical amount of disorder σc (L), the magnetization jumps become power-law-distributed. Nevertheless, this third possibility for the existence of “avalanche”-like phenomena remains unclear since the extrapolation of σc (L) to the thermodynamic limit is difficult. A number of authors [8] have shown that in this 2D case the ferromagnetic phase is unstable for all σ > 0. In order to circumvent this problem, in this paper we have studied the 3D G-RFIM, for which the existence of a ferromagnetic ground state is more generally accepted.
2. Model and numerical details The model is defined on a cubic lattice of size N = L3 . At each site we define a spin variable Si (i =
1, . . . , N ) taking values +1 or −1. The Hamiltonian reads: H =−
n.n. i,j
Si Sj −
N
Si hi − B
i=1
N
Si .
(1)
i=1
The first term, which extends over all nearest-neighbor pairs, stands for a ferromagnetic interaction. The second term stands for the interaction with the local quenched random fields hi , which are independent and Gaussian-distributed with zero mean hi = 0 and variance h2i = σ 2 . The last term accounts for the interaction with the external field B. Efficient algorithms for the study of the exact ground states of (1) have been available for several years. We have used a max-flow min-cut algorithm that for a certain fixed set of the random fields {hi }, and for each value of B in the whole range −∞ < B < ∞ finds the configuration {Si } that minimizes H . The algorithm used [9] determines not only the sequence of ground-state configurations, but also the exact values of B where the magnetization jumps s take place. We have systematically studied systems of sizes L = 8, 12, 16, 20, 24, 28 and 32. In order to compute the histograms, a large number of realizations of the random fields for each value of σ have been analyzed. For instance, for L 28 typical histograms are obtained from several thousands of configurations, whereas for L = 32 only several hundreds are considered. 3. Results Fig. 1 shows the log–log plot of the normalized distribution of “avalanche” sizes p(s) for systems with L = 28 and different values of σ ranging from 2.10 to 2.70. The values of s are discrete ranging from 1 to 283 . Note that data has not been binned so that the discretization of the data in the vertical direction is related to the statistical confidence which, in some cases, is of the order of 10−5 . As can be seen, for the largest σ values the distribution of “avalanche” sizes extends over the small s region whereas for lower values a peak corresponding to very large “avalanches” (with a size of the order of L3 ) appears. In order to model such histograms we have performed maximum-likelihood fits of the two freeparameter expression: p(s) = A(λ, τ, L)e−λs s −τ ,
(2)
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
457
Fig. 3. Numerical estimations of σc (L) versus 1/L. The lines show different fits of Eq. (3) as indicated by the legend. Fig. 1. Distribution of “avalanche” sizes in log–log scale. Data corresponds to a system with L = 28 and different values of σ . The lines correspond to the fits of Eq. (2). The histograms have been consecutively shifted 2 decades upwards in order to clarify the picture.
Fig. 2. Parameters τ (a) and λ (b) obtained from the fit of the expression 2 to the numerical data, for system sizes ranging from L = 8 to L = 32, as indicated by the legend. The dotted lines indicate the values of σ for which λ = 0.
3 −λs −τ where A(λ, τ, L) = 1/ L s is the normals=1 e ization constant. The maximum-likelihood procedure, contrarily to the χ 2 methods, is independent of any binning process and of the type of representation (linear, linear–log or log–log) of the data. The results of the fits are also shown in Fig. 1. In general the fits are good in the small s region but are unable to quantitatively reproduce the peak in the large s region for small σ ’s. The resulting fitted values of the parameters λ and τ are shown in Fig. 2. As can be seen, for
any value of L, there exists a certain σc (L) for which λ = 0 and, therefore, the distribution of “avalanches” becomes power-law. For L = 28, for instance, this occurs for σc (L = 28) 2.6 as can also be seen from Fig. 1 when the peak corresponding to large s values disappears. Estimations of σc (L) can be obtained by interpolation to the λ(σ ) data close to the crossing point, using a linear fit. The obtained values of σc (L) in front of 1/L are shown in Fig. 3. At first sight, it becomes clear that in the thermodynamic limit a critical point with power-law distribution of “avalanches” exists. Before performing a more detailed analysis of the extrapolation to L → ∞ it is convenient to review some of the ground-state properties of the 3D G-RFIM which are not fully understood. In particular, the exact position σc of the T = 0 critical point and the values of the critical exponents characterizing it have been under discussion for a long time. Most of the previous studies of the 3D G-RFIM have been performed at B = 0. After the pioneering work of Ogielski [10] an intense debate focused on whether the transition is discontinuous [11] or continuous [12] with a very small value of the critical exponent β. A very recent and deep analysis [13] has tilted the balance towards the existence of a standard critical point at σc = 2.270 ± 0.004 with β = 0.017 ± 0.005 and ν = 1.37 ± 0.09. To our knowledge only a single study with an applied external field has been performed. MC simulations at finite temperatures [14] with L up to 24 are not totally conclusive, but the authors propose that, instead of a critical point at (Tc (σ ), B = 0), two twin critical points could appear at (Tc (σ ), ±Bc (σ )).
458
C. Frontera, E. Vives / Computer Physics Communications 147 (2002) 455–458
These are joined by first-order phase transition lines to a triple point at (T3 (σ ), B = 0). The extrapolation to T → 0 gives σ3 = 2.35. According to the standard finite size scaling ansatz σc (L) should behave like: σc (L) = σc + AL−1/ν .
(3)
In Fig. 3 we have compared our data with different fits of Eq. (3): (i) the best linear fit (ν = 1), (ii) the best fit imposing ν = 1.37 [13], (iii) the best three freeparameters fit, and (iv) the best fit imposing ν = 1.37 and σc = 2.27 [13] (restricted to the L 24 data). Although it is difficult to distinguish the quality of the fits (i) and (iii), the fits (ii) and (iv) are qualitatively worse. Moreover, the extrapolated values of σc for the fits (i), (ii) and (iii) are clearly above σc = 2.27 [13]. Our best (and very conservative) estimation is σc = 2.4 ± 0.1. Therefore, we conclude that the “critical” point characterized by a power-law distribution of magnetization jumps when B is swept from −∞ to ∞ occurs at a value above the critical point found from the studies at B = 0. Moreover, from Fig. 2, we can also obtain an estimation of the exponent τ . Although λ exhibits a clear finite size dependence, the value of the exponent τ (at the λ = 0 point) is almost size independent, as indicated by the dashed lines. Our best estimation is τ = 1.70 ± 0.07. Finally, it is also interesting to compare our data with the data from the studies of the same 3D G-RFIM using the metastable dynamics [5] that have found a critical point at σc = 2.16 ± 0.03 with τ = 2.03 ± 0.03. Therefore, our present results do not allow any connection to be established between both critical points.
Acknowledgements The authors acknowledge financial support from CICyT project numbers: MAT1999-0984 (C.F.) and
MAT2001-003251 (E.V.). Part of the data has been obtained by using a free distribution of the maxflow min-cut algorithm from LEDA Research Group (Germany). This research has partially been developed using computational facilities of CESCA (Centre de Supercomputació de Catalunya). References [1] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381; How Nature Works, Oxford. Univ. Press, Oxford, 1997. [2] D. Sornette, J. Phys. I (France) 4 (1993) 209. [3] J.P. Sethna, K. Dahmen, S. Kartha, J.A. Krumhansl, B.W. Roberts, J.D. Shore, Phys. Rev. Lett. 70 (1993) 3347. [4] K. Dahmen, J.P. Sethna, Phys. Rev. Lett. 71 (1993) 3222; E. Vives, A. Planes, Phys. Rev. B 50 (1994) 3839; E. Vives, J. Goicoechea, J. Ortín, A. Planes, Phys. Rev. E 52 (1995) R5; B. Tadi´c, Phys. Rev. Lett. 77 (1996) 3843; J. Kushauer, van Bentum, W. Kleemann, D. Bertand, Phys. Rev. B 53 (1996) 11 647. [5] O. Perkovic, K.A. Dahmen, J.P. Sethna, Phys. Rev. B 59 (1999) 6106. [6] Ll. Carrillo, Ll. Mañosa, J. Ortín, A. Planes, E. Vives, Phys. Rev. Lett. 81 (1998) 1889, and references therein. [7] C. Frontera, E. Vives, Phys. Rev. E 62 (2000) 7470. [8] E.T. Seppälä, M.J. Alava, P.M. Duxbury, Phys. Rev. E. 63 (2001) 066110, and references therein. [9] C. Frontera, J. Goicoechea, J. Ortín, E. Vives, J. Comput. Phys. 160 (2000) 117. [10] A.T. Ogielski, Phys. Rev. Lett. 57 (1986) 1251. [11] J.-C.Anglès d’Auriac, N. Sourlas, Europhys. Lett. 39 (1997) 473; N. Sourlas, Comp. Phys. Comm. 121 (1999) 183. [12] M.E.J. Newman, G.T. Barkema, Phys. Rev. E 53 (1996) 393; U. Nowak, K.D. Usadel, J. Esser, Physica A 250 (1998) 1; A.K. Hartmann, U. Nowak, Eur. Phys. J. B 7 (1999) 105. [13] A.A. Middleton, D.S. Fisher, Phys. Rev. B 65 (2002) 134411. [14] J. Machta, M.E.J. Newman, L.B. Chayes, Phys. Rev. E 62 (2000) 8782.