Stretched and non-stretched exponential relaxation in Ising ferromagnets

PHYSICA ELSEVIER

Physica A 232 (19961 171-179

Stretched and non-stretched exponential relaxation in Ising ferromagnets P e t e r G r a s s b e r g e r a,,, D i e t r i c h S t a u f f e r b aphysics Department, University of Wuppertal, 1)-42097 Wuppertal, Germany bInstitute for Theoretical Physics, Cologne University, D-50923 K61n, Germany Received 16 April 1996

Abstract

The decay of the Hamming distance in Monte Carlo simulations of Ising models with model A dynamics is shown to be numerically efficient to investigate the dynamics also below T~.. We observe a simple exponential decay in three and four dimensions below and above the Curie temperature. In two dimensions, a stretched exponential decay is observed only below To. The exponent in the latter is between the two predictions of Takano et al. and Huse and Fisher, but the dynamics of single droplets is clearly found to disagree with Takano et al.

During the last decade the method of damage spreading [ 1-7] has been widely used to investigate time-dependent critical phenomena. Damage spreading is a technique where two identical lattices are simulated with different initial conditions but with the same sequence of random numbers. If the spin difference (the "damage") decreases to zero, then the system is "stable" in the sense that it does not show sensitive dependence with respect to initial conditions for a fixed realization of the random number sequence. We speak in this case of damage healing. In the opposite case where a small initial damage spreads, the system can be considered to be chaotic. The detailed properties with respect to damage depend strongly on the used Monte Carlo scheme. Thus it seems that in the Glauber dynamics of the Ising model there is a damage spreading transition which does not coincide with the critical temperature [8,9], and which is in the directed percolation universality class [10]. The situation is completely different for heat bath dynamics (for the implementation of heat bath dynamics for lattice pairs involved in this statement, see [11]), as explained theoretically * Corresponding author. 0378-4371/96/$15.00 Copyright (~ 1996 Elsevier Science B.V. All rights reserved PH S0378-4371 (96)00206-3

P. Grassberger, D. Stauffer/Physica A 232 (1996) 171-179

172

by Coniglio et al. [12]. Consider two replicas A and B with spins S~A) and S f ), such that initially

S~A)>~Sf )

(t -- O)

(1)

for all sites i. Then it follows that S~A)>~Sf ) for all later times as well. This means that the damage is equal to the magnetization difference for any two replicas which satisfied Eq. (1) initially, and it allows to measure very small magnetization differences [12-14]. This was used in [15-17] to measure the dynamical exponent z, and in [18] to measure the magnetization decay in the Griffiths phase of diluted Ising systems where a slow decay is expected which is between power law and stretched exponential [19, 20]. We use this method here to measure the magnetization decay for pure Ising models at temperatures ~ Tc, with initially fully magnetized configurations, M(0) = 1. In 2 dimensions it is known that the decay of the excess M ( t ) - M ( o c ) is not a simple but a stretched exponential [21-27]

M(t)-M(c~)

~ e -eta,

with 0 < ~ < 1

(T < Tc, d - - 2 ) .

(2)

The same decay is believed to govern also the autocorrelation function C(t) = (Si(t) si(0)) - (s) 2. The qualitative physical reason for Eq. (2) suggested in [21,22,24,25] is the following1: for large t, the only contributions to C(t) are from droplets with reversed magnetization which are large enough to have a lifetime >t t. But there exist conflicting assumptions about the lifetime distributions [21,22, 24, 25] which have lead to different conjectures for the exponent ~. While most authors agree upon ct -- 1 [22, 24, 25], the value ~ = 1 had been proposed in [21]. This difference is even more pronounced in higher dimension. While a stretched exponential with ~ = ( d - l )/(d + 1) was predicted for any d~>2 in [21], a simple exponential decay M ~ e -ct was derived in [22,24,25] for d ~>3 (eventually with correction terms for small t [24], which should however be different for C(t) and M ( t ) - M ( c ~ ) ) . Above Tc, all authors so far have assumed simple exponential decay in all dimensions ~> 1, though we are not aware of any rigorous proof for that. Simulations so far gave simple exponentials above T~ in all investigated dimensions, and stretched exponentials below Tc for d - - 2 [27]. The simulations of Ito [27] also prefered ct = ½ over the alternative prediction ~ = ½, but this was not very significant due to the inherent difficulties in distinguishing different stretched exponentials, in particular if they involve also power law or logarithmic prefactors. For d = 3, the situation was also not very clear. While a simple exponential was prefered in [26], a stretched exponential could not be excluded. 1 As pointed out in [25], it is not really needed that individual droplets live at least t time steps, provided these droplets are correlated. Moreover, the decay of M(t) (as opposed to that of C(t)) follows more easily if the lifetime distribution Pn(t) of individual droplets of size n used in [21,24] is re-interpreted as droplet size distribution at time t. In that case, the Langevin equation used in [22, 24] becomes essentially equivalent to a Langevin formulation of the Fokker-Planck-type equation postulated by Becker and Doering [25].

P. Grassberoer, D. Stauffer l Physica A 232 (1996) 171-179

173

In our simulations, we set up two lattices of identical size with all spins up, and let one of them (called in the following the "reference replica") relax to a state near equilibrium. Then we let both replicas (the reference replica and the one where still all spins are up) relax further using the heat bath method with identical random number sequences, and monitor the Hamming distance ("damage") between them. When this Hamming distance has reached zero, we can stop the iteration since it would remain zero forever. Thus, we reset all spins up in the second replica, keep the reference replica in a state close to equilibrium, and start the iteration again. To "keep the reference replica in a state close to equilibrium", the most straightforward strategy would be just to continue with its last spin configuration. This would however introduce a bias since the damage has died out preferentially when M was large in the reference replica, and the new pair of replicas would start in a biased state. To avoid this bias, we followed two strategies. In the first we used just very large lattices (up to 20000 x 20000 in d = 2 , e.g.). In the second we stored the configuration of the replica to time steps after a new pair was started (with to roughly equal to the average damage lifetime; obviously this meant that at least the iteration of the reference replica has to be followed for ~>t0 time steps, even if the damage had already healed earlier), and used this configuration for the next round. Here time t is measured as usual in Monte Carlo steps per spin. No significant differences were found between these two methods. We used checkerboard updating, but our conclusions do not seem to depend on these programming details. The main advantage of using damage healing is of course that statistical fluctuations are strongly suppressed [12]. In addition, at least two potential sources of systematic errors are eliminated to leading order, as both affect both replicas in the same way: finite size effects and imperfections of the random number generator (RNG). One just should take care to use a RNG with sufficiently long period since the heat bath algorithm is not "chaotic", whence using a RNG with short period would also lead to periodic configurations. We used R250 and a 64-bit linear congruential one, with no discernable difference. For the simulations below To, the lattices must be sufficiently large that the total magnetization n e v e r flips during the entire simulation. We checked that the magnetization of the reference replica always stayed far from zero, even at our highest temperatures. Fig.1 shows data for the damage decay above To, in 2 and 3 dimensions, on a semilogarithmic plot. We see in both cases nice straight lines, indicating simple exponentials with very little corrections. Corrections are of course expected very close to To, but they do not yet show up in Fig. 1. The results shown in Fig. 1 are of course not surprising 2, but they illustrate nicely the accuracy of our method.

2 Stretched eponentials for T > Tc are expected in 1 dimension, where a decay with ~ = ½ holds for finite times [28]. When T ~ 0 (and at~er a suitable renormalization of the time scale), it holds for arbitrarily long times.

P. Grassber#er, D. Stauffer/Physica A 232 (1996) 171-179

174

d = 2, 10000xl0000sites(squares)

at T/'1"c=1.32,

and

d = 3, 80,k3sites (crosses) at T/Tc=l.33

le+10 O

le+09

1e+08

[] + 4-

le+07

+ (3

+ 4-

E "o

4-

1e+06

+ 4+

rrl

÷

100000

+ E]

+ + []

+ +

I0000

+ []

+ + +

1000 0

I 20

I 40

~ 60

f 80 time

I 100

i

I 120

t 140

160

Fig. 1. Simple exponential decays in two and three dimensions above Tc. Data are from 8 runs on lattices of size 100002 (kT/J = 3; squares) and from 27200 runs with 803 (kT/J = 6; crosses).

More interesting is the decay of the damage below To, where it corresponds to the excess magnetization M ( t ) - M ( ~ ) . This is plotted in Fig. 2 again on a semilogarithrnic plot. Data are shown for dimensions 2, 3 and 4, with T ~ 0.9To in all three cases. While the data for d = 3 and d --- 4 are compatible with being simple exponentials (at least on this scale), the data for d = 2 obviously cannot be described in this way. In order to understand better the behavior in d~>3, we formed the ratio [26]: Q = In D(t - 5) - D(t) O(t) - O(t + 5) '

(3)

where D(t) is the damage (Hamming distance) at time t. It should approach a positive constant for a simple exponential decay, and should decay as t "-1 if D(t) ,,~ e -ct~. Results are shown in Fig. 3. They indicate clearly that simple exponentials prevail asymptotically both above and below T~ for d ~>3, though there are large finite-t deviations near and below T~ (notice that Q ~ 0 for T --* To). Similar finite-t corrections are predicted in [24], but for C(t) instead o f M ( t ) - M(cx~). If these corrections were just due to a cross-over to the power decay holding at T = Tc, then they should be described by a power t -#/2v ~ t -°'25. But things cannot be so simple, otherwise this prefactor should appear also in the data for T > Tc where it is not seen [29]. Indeed, a best fit gives for d = 3 below T~ an exponent ~-, 0.6 instead of 0.25. Thus, the effective

P. Grassberger, D. Stauffer l Physica A 232 (1996) 171-179

175

1000"1000 (diamonds), 80*80*80 (+), and 40*40*40*40 (squares) near T/'l'c=0.9

lo+09

i

i

i

i

P le+08

le+07

8 +0 <>

t~

E

[] +

o

lo+06

0

1:7

<> <> <>

100000

o

<>,¢, ,3<>

Itt

,¢, o 0

~><> <>o <> < > o O

10000 + [3

00<>0<} 00 <><>0~<>O0 +

1000

0

I

I

50

100

1

I

150

200

9

time

Fig. 2. Decay below Tc, based on 10000 runs for 10002 (kT/J = 2; diamonds), 20000 rims for 803 (kT/J = 4; crosses) and 5000 runs for 404 (kT/J = 6; squares). (Larger lattices up to 200002 ,4723 and 724 with less statistics confirmed this behavior.) exponent o f the prefactor seems to increase with decreasing T. In spite of this problem we consider our data as very strong evidence for simple exponentials in d >~3, both below and above Tc. Finally, we consider the case of 2 dimensions below Tc. In Fig. 4 we plot log M ( t ) for various values o f T against x/t (panel a), resp. against t 1/3 (panel b). We do not see perfect straight lines in either plot. This might be interpreted as a stretched exponential with ~ between ½ and :. i Alternatively, we could again try a prefactor 1/t a as in d -- 3. If such a prefactor is fitted in combination with ct = 1, this gives rather poor fits with a ~ - 0 . 5 . In contrast, with c¢ --- ½ we find excellent fits for a = 0.25. Neither of these exponents agrees with the a = f l / 2 v ~ 0.06 holding at T = Tc. But this should not be expected after our experience with d -- 3. Comparison with the latter would suggest that also for d = 2 the effective exponent a might increase with decreasing T, which would prefer ct = ½ over ct = 1. This would also be in agreement with [22, 24, 25]. All this prefers the prediction o f these papers over those o f [21], but these arguments are rather weak, and it seems still difficult to decide between these alternatives on the basis o f direct simulations. For indirect evidence in favour o f one or the other prediction, we finally simulated the evolution o f individual droplets (for a similar study see [30]; there, much larger droplets were studied which were too big to play an important role in our relaxation

P. Grassberger, D. Stauffer/Physica A 232 (1996) 171-179

176

Ito ratios for d=3 (kT/J = 4 (+), 4.30 (x), 4.35 (triangles), 6 (diam.)) and d=4 (T
i

I

I

I

I

2.5

..I-

2

o,

x z~

13

0v

""

[]

1.5

~

[]

[]

[]

[]

O

+

o ,.a... to

x ~x x

+

zx x

zxx

=

A× X ~

0.5

X X

A~ XXX X ~A~ A X X X X X X x x x ~ A XXXXXxxxxxxxxxxxxxxxXxxxx

I

I

10

20

I

30 time

X~

I

I

I

40

50

60

Fig. 3. Logarithm of the Ito ratio ln[((M(t-5)-M(t))/((M(t)-M(t+ 5))] for d = 3. This ratio approaches a constant for ct --- 1 and decays as t °-1 for ~t < 1. We use the same data as in Fig. ! above and Fig. 2 below Tc (but omit short times where the deviations from the plateau are stronger). In addition, we also show data slightly below Tc (T = 4.3 and T = 4.35; Tc = 4.511) where Q is much smaller but also seems to converge to a constant > 0. The latter were each from about 106 runs on lattices 23 x 23 × 24. studies, and whose evolution could not be studied for sufficiently long times to see the asymptotic behavior). W e proceeded as before, except that the spins in the second replica were not all set "up" initially. Instead, they were left unchanged except for the spins in a roughly spherical droplet which were all set "down". Thus, we created initial states with droplets o f inverted spins, and let them grow or shrink. According to all theories cited above, the average lifetime o f a droplet with initial radius R should scale as R 2. But while [21] assumed a roughly exponential decay M ( t ) ~ e -at~R2 o f the average magnetization caused by such a droplet, the more detailed analysis o f [24] gives M ( t ) ~ f(R)exp(--ct(d-1)/2), with f ( R ) such that (t) _= f d t t M ( t ) /

f d t M ( t ) ~ R 2. Our results are shown in Fig. 5 for d = 3, and in Fig. 6 for d = 2. Each curve corresponds to a different initial droplet size. We see very clear stretched exponential decays with exponents exactly as predicted in [24]. The decay assumed in [21] is clearly not seen. One might take this as a strong argument in favor o f [22, 24, 25] and against [21], but we should warn the reader that things are not so simple. Our droplets are presumably not equal to the droplets appearing in equilibrium, and thus the connection between the decay o f our artificial droplets and the magnetization decay is not as straightforward as one might hope.

P. Grassberger, D. StaufferlPhysica A 232 (1996) 171-179

177

2-D Ising model, T < T c, heat bath dynamics -1

I

I

I

~

i

i

~,~.~,. iji~-.%.~ ~L. ~ ~l= ~

-2

~

~

2,174 2.128 2.041 2.000 1.961

-e---+-.n.. x ~--

T = 1.887 ~

~,~.~ ~ , ~ ~

-3

I

T = T = T = T = T=

T=175 .....

-4

o

-5

-6

-7

-8

0

(a)

I

I

I

I

5

10

15

20

t1/2

I

I

I

25

30

35

40

2-D Ising model, T < T c, heat bath dynamics -1

i

T

i

i

~-.~'%."~%~

~E*'~'~ , ; ~ . ~

-2

•, ~ , , ~ ii~"~ 'I,~

-3

i

T = 2.174 - e ~ T = 2.128 . . . .

T = 2.041 T = 2.000

~. ~ ~,~ ~

o

T = 1.961 ~ T = 1.887 - ~ T = 1.75 -m--

I

%

-4

"-~

-5

°

\,.

-7

-8

(b)

0

I

I

I

I

I

2

4

6

8

10

t 1/3

12

Fig. 4, Stretched exponential decay in two dimensions below Tc. Data are from very many runs on small lattices (81 × 82), but they are fully compatible with the large lattice simulations shown in Fig. 2, and with 3 runs on a 200002 lattice with T = 2. From bottom to top, the temperatures are 1.75 (8 × 106 runs), 1.887 (7 × 106 runs), 1.961 (7 × 106 runs), 2.0 (6.7 × 106 runs), 2.041 (5 × 106 runs), 2.128 (4.6 × 106 runs), and 2.174 (3 x 106 runs). We plot logarithmically the excess magnetization versus x/t in part a, and versus t 1/3 in part b.

P. Grassberger, D. Stauffer/Physica A 232 (1996,) 171-179

178

single droplets in 2-D Ising model, T = 2.128, heat bath dynamics

-0.5

37 spins J

-1.5 o

-2.5

-3

0

,

,

,

,

2

4

6

81/2

,

,

10

12

14

16

t

Fig. 5. Average magnetization excess due to a single droplet with all spins pointing opposite to the global magnetization, for d = 2. Temperature is T = 2.128; initial droplets consisted of no sites with no = 5, 13, 21, and 37 (from bottom to top). On the horizontal axis is plotted v~. single droplets in 3-D Ising model, T = 4.35, heat bath dynamics 0

i

i

-t .5

i

"',,.

-2.5 7 spins -3

-3.5 20

t 40

= 60

80

100

120

t

Fig. 6. Similar to Fig. 5, but for d : 3. Temperature is T = 4.35, initial droplets consisted of no sites with no = 7, 33, and 81 (from bottom to top). On the horizontal axis is plotted t.

P. Grassberger, D. Stauffer/Physica A 232 (1996) 171 179

179

In s u m m a r y we have studied relaxation in the Ising model up to times where the excess m a g n e t i z a t i o n has decayed to extremely small values. A b o v e Tc (where the situation n e v e r was controversial) we found exponential decay as expected. In the symmetry b r o k e n phase b e l o w Tc our results clearly support the predictions o f [22, 24, 25], at least for d >~3. T h e y speak against the arguments proposed in [21] for an alternative scenario, but they still c a n n o t exclude this scenario directly.

Acknowledgements W e thank H L R Z J/ilich, R R Z K C o l o g n e and R Z der B U G H Wuppertal for c o m p u t e r time, and N.lto for discussions. This work was partially supported b y D F G , SFB 237 and 341.

References [1] M. Creutz, Ann. Phys. (NY) 176 (1986) 62. [2] N.Jan and L. de Arcangelis, in Annual Reviews of Computational Physics, Vol. I, ed., D.Stauffer (World Scientific, Singapore 1994) p.l. [3] 13. Derrida and G. Weisbuch, Europhys. Lett. 4 (1987) 653. [4] H.E. Stanley, D. Stauffer, J. Kert6sz and H. Herrmann, Phys. Rev. Lett.59 (1987) 2326. [5] M.L. Martins, H.F. Verona de Resende, C. Tsallis and A.C.N. de Magalhaes, Phys. Rev. Left. 66 (1991) 2045. [6] I.A. Campbell, Europhys. Lett. 21 (1993) 959. [7] P. Grassberger, J. Statist Phys. 79 (1995) 13. [8] U.M.S. Costa, J. Phys. A 20 (1987) L583. [9] G. le Ca~r, Physica A 159 (1989) 329; J. Phys. A 22 (1989) L647. [10] P. Grassberger, J. Phys. A 28 (1995) L67. [11] A.M. Mariz, H.J. Herrmann and L. de Arcangelis, J. Statist. Phys. 59 (1990) 1043. [12] A. Coniglio, L. de Arcangelis, H.J. Herrmann and N. Jan, Europhys. Lett. 8 (1989) 315. [13] S.C. Glotzer, P.H. Poole and N. Jan, J. Statist. Phys. 68 (1992) 895. [14] R.Matz, D.L.Hunter and N.Jan, J. Statist. Phys. 74 (1994) 903. [15] P. Grassberger, Physica A 214 (1995) 547; Physica A 217 (1995) 227 (erratum). [16] U. Gropengiesser, Physica A 215 (1995) 308. [17] F.Wang and M.Suzuki, Physica A 223 (1996) 34. [18] P. Grassberger, Int. J. Mod. Phys. C (1996), to be published. [19] A.J. Bray, Phys. Rev. Lett. 60 (1988) 720. [20] D. Dhar, M. Randeria and J.P. Sethna, Europhys. Lett. 5 (1988) 485. [21] H.Takano, H.Nakanishi and S.Miyashita, Phys.Rev. B 37, (1988) 3716. [22] D.A.Huse and D.S.Fisher, Phys.Rev. B 45 (1987) 6841. [23] A.T. Ogielski, Phys. Rev. B 36 (1987) 7315. [24] C. Tang, H. Nakanishi and J.S. Langer, Phys. Rev. A 40 (1989) 995. [25] D. Stauffer, Physica A 186 (1992) 197. [26] N. lto, Physica A 192 (1993) 604. [27] N. lto, Physica A 196 (1993) 591. [28] J.J. Brey and A. Prados, Phys. Rev. E 53 (1996) 458. [29] W. Koch, V. Dohm and D. Stauffer, to be published (1996). [30] D. Chowdhury, J. Phys. A 20 (1987) LI171.