Damage spreading and critical exponents for “model A” Ising dynamics

PHY$1CA ELSEVIER

Physica A 214 (1995) 547-559

Damage spreading and critical exponents for "model A" Ising dynamics Peter Grassberger Physics Department, University of Wuppertal, D-42097 Wuppertal, Germany Received 24 October 1994; revised 15 November 1994

Abstract Using damage spreading and heat bath dynamics, we study the Ising model in 2 and 3 dimensions with non-conservative dynamics. Our algorithm differs in some important points from previous ones, which makes it rather efficient. We give estimates for the exponent z which seem to be the most precise published so far (2.172-t- 0.006 for d = 2, 2.032-t- 0.004 for d = 3). We also give precise estimates of the exponent 0' introduced by Janssen et al. (Z. Phys. B 73 (1989) 539) and of analogous but in principle independent exponents. We find surprisingly that some of the latter agree with 0', and give an explanation for this.

1. Introduction In this paper we shall be concerned with the critical dynamics of the Ising model when neither energy nor the order parameter is conserved ("model A" according to [ 1 ] ). After some conflicting numerical results (see [2] for a review), this was considered as essentially understood when high precision Monte Carlo measurements in 3 dimensions [ 2 ] agreed very well with field theoretic predictions [ 3 ]. But more recently the situation has become somewhat confused again. In 2 dimensions the exponent z (which is defined as r ~ ~:z; T and sc are temporal and spatial correlation lengths) was found in several simulations [ 4 - 8 ] to be considerably higher (2.165-2.5) than the value 2.126 predicted in [3]. In 3 dimensions the same is true. There values 2.05 - 2.35 were measured in [9,5,10,11,7,8], while z = 2.013 had been predicted in [3]. Another interesting aspect is that Janssen et al. [12] showed that there are new universal critical exponents if one starts with non-equilibrium initial states which are suddenly quenched to T = Tc. Finally, using the concept o f "damage spreading" [ 140378-4371/95/$09.50 (~) 1995 Elsevier Science B.V. All rights reserved SSDI 0378-437 1 (94)00285-1

P. Grassberger/Physica A 214 (1995) 547-559

548

16] (which will be explained below), it was shown [ 17] that there are new critical phenomena in Ising models simulated by means of the Glauber algorithm which do not seem to be related to the standard Ising transition. No such new phenomena are seen with heat bath dynamics. In the present paper we show that the situation becomes simpler if the above problems are tackled simultaneously. In particular, we show that very precise estimates of z and of the new exponents of [ 12] can be obtained by using heat bath damage spreading with non-equilibrium initial conditions. We should point out that damage spreading has been used already in [4,13,5,11,18] for measuring z, but with considerably different details. These authors started with different initial states, they fixed some of the spins during the simulations, they measured different observables, and they used different multispin coding or no multispin coding at all. All this together leads to our results being more precise in spite of modest CPU time. In particular, they are not affected by any logarithmic corrections like those which had plagued some of the earlier methods. But our simulations will also bring forward a new surprising result. Let us first recall the arguments of [ 12]. Consider a system which is equilibrized at T = cx~ in a strong magnetic field with B / T small but finite. Thus all spins are uncorrelated, but the net magnetization M0 is small and non-zero. At time t = 0 the magnetic field is switched off and the system is quenched suddenly to T = Tc where Tc is the Curie temperature. As shown in [12], this will lead to an increase of the magnetization which follows a power law,

M ( t ) = Mot ° ' ,

0'>0.

(1)

This persists up to some time t*, after which M ( t ) decays with the equilibrium exponent

M ( t ) oct -#/~z.

(2)

Requiring both expressions to be satisfied (with a proportionality constant O(1) in Eq. ( 2 ) ) gives the cross-over time t* which thus scales with M0 as

t* ~ M o 1/[#/~z+°'].

(3)

While 0' was calculated in [ 12] up to O ( e 2) for d = 4 - e, it was also estimated for d = 3 by Monte Carlo Simulations in [ 19]. We finally mention that another exponent (called 0 in [ 12] ) describes the initial growth of M if B = 0, i.e. if M0 = 0 except for statistical fluctuations. It was shown in [ 12] that 0 is not independent of 8'. Let us now see how these results are related with damage spreading. In the following we shall exclusively deal with the heat bath algorithm (damage spreading with Glauber dynamics will be treated in [20] ).

P. Grassberger/Physica A 214 (1995) 547-559

549

As pointed out by Coniglio et al. [21 ], heat bath dynamics has the following important property called "monotonicity" in [22] 1. Assume we have two replicas of the same lattice, let us call them A and B, and assume that we have at time t = 0 Sx(i) (0) >_ S(a) (0)

(4)

for each site x. Thus each "up" spin in B is also "up" in A, while the opposite is not necessarily true. It is straightforward to check that Eq. (4) remains true for all subsequent times, if we use the heat bath algorithm with identical sequences of random numbers for updating A and B. For a precise description of the heat bath algorithm, see [23]. In the following we implement it such that we alternatively update simultaneously all spins on even and odd checkerboard sublattices. Let us now assume that the initial configurations of A and B are identical except for a single spin at x = 0 for which we assume s(A) ( 0 ) = + l / 2 , S0~a)(0) = - 1 / 2 (only So is "damaged"). Then the response function is defined as C ( x , t) = (Sx~A)(t) - S~B) (t)). Due to monotonicity, we have however Sx~A)(t) > S~xB)(t) for all x and t, whence we have also C ( x , t) = (IS~(A)(t) - S~B) (t)[).

(5)

Thus the response function is obtained by observing the "damage" at all positive times. The connection with the critical exponents is established by observing that Eq. ( 1 ) must also hold for the magnetization difference (M(g)(t)

-

-

M(B)(t)) = Z

C(x,t),

(6)

X

and making a standard scaling ansatz (d = spatial dimension) C ( x , t ) -,~ ta'-d/z~b(rZ/t)

,

r = Ixl.

(7)

The exponent of the prefactor follows directly from Eqs. (1) and (6). Less trivial and crucial for the following - is that the exponent z appearing in the argument of the scaling function is the same as for equilibrium dynamics. This is an important result of the field theoretic treatment of [ 12]. It is very similar to the well known result that the correlation length in critical systems with "ordinary" fiat boundaries is the same parallel and perpendicular to the boundary [24]. Our first strategy ("strategy I") will thus consist in creating a single damage and observing the growth of the total amount of damage and of its average distance from the origin, R~ = ((s~(A)(t) -- S(~a)(t))r2).

(8)

The latter should scale as R2t ~ t 2/z. l In contrast to what was suggested by an anonymousreferee to Coniglio et al., it has nothing to do with the FKG property.

550

P. Grassberger/Physica A 214 (1995) 547-559

Unfortunately, this strategy is not very efficient. In order to see the scaling behavior without being molested by finite size corrections we need rather large lattices. But in each lattice only the region near the origin is really used, while the updating of the overwhelming rest is just a nuisance. There are several ways how to overcome this problem. The first ("strategy II") consists in damaging not a single site at t = 0 but a randomly chosen finite fraction of sites. In this way it becomes difficult to measure R 2 (we would have to remember which was the "ancestor" of every damaged site), but it is a convenient way to measure 0'. Our last strategy (III) consists in damaging entire linear subspaces of the lattice. In 2D this means that we take S~A) (0) = 1 = -S~ B) (0) on an entire line, while in 3D we can either damage a line or a plane. In each case we should expect R 2 ~ t 2/z for the spreading away from the subspace, but we cannot any longer expect that the total damage increases with the exponent 0'. This is again completely analogous to surface exponents in static critical systems. Different exponents were e.g. measured in [25] for percolation with different boundary conditions.

2. Implementation We studied the square and the simple cubic lattices, both with helical boundary conditions. The latter were used because they allow easily to label each site by a single integer (instead of d integers). Since we divided the lattices into two checkerboard sublattices for updating, we used lattices of size L x (L + 1 ) resp. L 2 x (L + 1 ), with odd values of L up to 229 in d = 3 and L < 1771 in d = 2. For d = 2 we used the analytically known value of Tc, while we used J/kBTc = 0.221655 [2,26,27] for d = 3. All simulations were done on 64 bit machines (DEC ALPHA workstations) where a convenient random number generator is obtained by ni+l = 1313 ni mod 264 [28]. We used this because of its speed and long periodicity, and because of recent problems with shift register algorithms [ 29,30]. We used multispin coding where each site corresponds to a different word in a 1D array of size L d - I x ( L + 1 ), and each bit of such a word corresponds to a different lattice (for a similar code and for references to earlier papers, see [6] ). This is particularly convenient for damage spreading since we can then use the same sequence of random numbers for each successive pair of bits in each word. Due to storage limitations most simulations were not done with 64 simultaneous lattices but with 32, i.e. with 16 pairs for which the mutual damage was measured by counting the number of discordant bits. In order to speed up the simulations we indeed used the same random numbers for all pairs at those sites which were not damaged in any pair, while we used different random numbers at those sites where at least one pair was damaged. Since we used independent random initial conditions for each pair, the evolution of different pairs was not identical even where there was no damage. The set of damaged sites was monitored by means of lists of damage sites similar to the lists of growth sites used for efficient simulations of percolation growth [25]. With this algorithm we could update 107 - 4 x 107 sites per

P. Grassberger/Physica A 214 (1995) 547-559

551

second on a DEC A L P H A 300AXP work station, depending on the word length and on the fraction of damaged sites. We count each updating of one sublattice as one time step (and also count the damage only on this sublattice). The number of so defined time steps was 103 in d = 3 and 2 x 103 in d = 2 for most simulations. We made no corrections for finite size effects. The latter should become important only for t ~ L z which was far from the values used in our simulations. But we also checked for the absence of finite size corrections by making runs on smaller lattices, and by monitoring the m a x i m a l distance of damaged sites. The latter was always less than L / 2 , indicating that the finiteness of the lattice could not have made any effect. This showed clearly that any finite size correction, if present at all, should be negligible compared to our statistical errors 2. In this respect our simulations differ most from those of [ 11,19] where the finite size behavior was an essential ingredient in the analysis. The total amount of CPU time for all simulations described below was ~ 1000 h.

3.

Results

In Fig. 1 we show a log-log plot of R 2 / t against t for an initial line damage on a 2D lattice. The data are averaged over > 3000 runs on square lattices with 263 < L < 723. For z = 2 we would have a horizontal line. Otherwise, its slope should be - ( z - 2 ) / z . We see that definitely z > 2, but there are considerable corrections to scaling. In order to take these into account, we define local effective exponents 2 - Zeff -

-

Zeff

ln[R2(2t)/4R2(t/2)r] -

In 4

(9)

and plot them against a suitable power of 1/t. This power should be chosen such that we obtain a straight line, as this allows the easiest extrapolation to t ~ oo and also gives an estimate o f the leading confluent singularity. Such a plot (with abscissa 1/t °'52) is shown in Fig. 2 and gives z = 2.172 + 0.006

( d = 2).

(10)

We verified that this was compatible with spreading from a single damaged site, though the error there was larger by about one order of magnitude. Similar results for spreading from a plane in d = 3 are shown in Fig. 3. This time the lattice sizes were 179 < L < 229, and the data are averaged over 70 runs. We now see virtually no corrections to scaling. Correspondingly, our error bars are now somewhat smaller in spite of somewhat worse statistics: z =2.032-t-0.004

(d=3).

(11)

2 In [6] finite size effects were supposedly seen for 500 time steps on lattices of size > 1000 x 1000 with helical boundary conditions. It is easily seen that such finite size effects are strictly impossible (they are identically zero for L _> 2t), whence these were most likely statistical fluctuations.

P. Grassberger/Physica A 214 (1995) 547-559

552

2-D 0.55

.

.

.

.

.

.

.

.

Ising m o d e l , i

.

.

.

.

h e a t bath,

line d a m a g e

.

I

.

.

.

.

.

.

.

.

.

.

.

I

0.5

0.45

A

0.4 V

0.35

0.3

. . . . . . . .

' 10

. . . . . . . .

' 100

. . . . . . . .

' 1000

t

Fig. 1. Log-log plot of R2/t against t for the spreading of damage away from a line embedded in a 2D lattice.

2-D -0.045

Ising m o d e l , ,

line d a m a g e ,

,

,

i 0.2

i 0.25

-0.05

-0.055

rr

-0.06

E

o

-o.065

>=

g

-0.07

-0.075

-0.08

-0.085 0

i 0.05

i 0.1

i 0.15 1 / t °'52

Fig. 2. Effective exponents for the ratio R2(t)/t defined in Eq. (9) plotted against 1/t °'52. The straightness of the data suggests that the leading correction to scaling has an exponent ~ 0.52.

P. Grassberger/Physica A 214 (1995) 547-559

553

3-D IsJngmodel, heat bath, plane damage

0.36

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

i

0.34

A rrV

0.33

0.32

0.31

0.3

.

.

.

.

.

.

.

.

'

.

~

,

10

,

,

, , , I

100

~

,

.

.

.

.

.

.

i

IO00

t

Fig. 3. Same as Fig. 1, but for damage spreading from a plane embedded in 3D.

Here we also made runs with single sites and with lines initially damaged. Both are compatible with these results but much less significant statistically. While it is not surprising that the single site runs gave very big errors, this was much less obvious for the runs with line damage. Part of the reason for the latter became clear when looking at the amount of damage, i.e. at the total difference in magnetization. If we study the spreading of damage from a d'-dimensional subspace of a ddimensional lattice, we could have guessed that the total damage contains a factor Re -d' ~ t ( d - d ' ) / z for purely geometric reasons. If we start from a high dimensional subspace, we start already from a large damage which can grow only into fewer directions. We should thus expect it to grow slower than if we had started, say, from a point. Certainly this argument is not rigorous. It is not quantitatively exact, e.g., for spreading of percolation away from a fiat boundary [25], but there it is at least qualitatively correct. Surprisingly, this seems not to be true for damage spreading. Both in 2 and in 3 dimensions we found that the total damage increases with the same exponent O', independent of the dimensionality of the subspace on which the initial damage was localized! For d = 2 this is shown in Fig. 4, while it is shown in Fig. 5 for d = 4. While even the corrections to scaling seem independent of the type of initial damage in d = 2, at least they differ in d = 3 and lead to a slower growth away from a damaged plane. Indeed, even our best estimates for O' are slightly different for point and plane

554

P. Grassberger/Physica A 214 (1995) 547-559

2-D Ising model, heat bath, point/line damage 5

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

I

4.5 4 3.5

A

2.5

E v

2

1.5

1

I

.

.

.

.

1

.

.

.

.

i

i

,

,

,

10

,

, , , I

100

.

.

.

.

.

.

.

.

I

1000

t

Fig. 4. Total amount of damage (defined by Eq. (6)) divided by the initial amount of damage, for d = 2. Full line: initially damaged region was a line; long dashed line: randomly placed multiple point damages; short dashed line: single point damage.

damage (0.105 vs. 0.102), but we believe that this is due to uncertainties in taking into account these corrections. An explanation of this universality of 0 ~ will be given in the next section. Indeed, in Figs. 4 and 5 the "point damage" data were not obtained from lattices with single damaged points, but from lattices where a small but finite fraction of randomly chosen sites was damaged (strategy II). Notice that this is indeed much closer to the original scenario of Janssen et al. [ 12]. We found no difference beyond the statistical errors when varying this initial density between 2 -13 and 2 -6 (notice that the latter is still smaller than the initial densities used in [ 19] ). The bulk of the "point damage" data in Figs. 4 and 5 is obtained with initial density 2 -8. The exponents obtained from Figs. 4 and 5 are O' = [ 0.191 5:0.003 L 0.104 ± 0.003

(d = 2) (d = 3).

(12)

These values agree perfectly with the less precise estimates of [31,32] for d = 3, and of [19] for d = 2 . To convince ourselves that the independence of O~ of the initial damage is not trivial we performed finally a last set of runs with completely different damage. This time both lattices were initialized completely independently and randomly. Thus in average half of the spins are damaged, but no attempt is made to control the sign of the damage. Thus

555

P. Grassberger/Physica A 214 (1995) 547-559

3-D Ising model, heat bath, point/line/plane damage

2.4

.

.

.

.

.

.

.

.

i

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

i

2.2 2 1.8

1.6

A

8~ 0~ "o

1.4

v

1.2

,

i

,

,

l l l , l

h

,

. . . . . .

10

I

100

.

.

.

.

.

.

.

.

i

1000

t

Fig. 5. Similar to Fig. 4, bat for d = 3 and for plane damage (full line), randomly placed multiple point damages (long dashed line), and single point damage (short dashed line).

monotonicity is of no immediate help. We expect that the damage should decrease in this case with some power,

Is~A)(t) -- S~B)(t)I ~

t -°''

.

(13)

x

The data shown in Fig. 6 indicate that this is indeed correct for d = 3. They moreover demonstrate that scaling is only observed exactly at the Curie temperature. There we find an exponent

O" = 0.43 + 0.02.

(14)

We have not found any way how to correlate this with other critical exponents. In particular, it differs from the exponent fl/z~' = 0.25 which governs the decay of the magnetization form a completely ordered state. Fig. 5 verifies also that the threshold for damage spreading agrees with the Curie point (which it does not for Glauber updating [ 17] ), and that the value of ~'t - defined as 7- ~ IT - T c I . . . . agrees with the equilibrium value.

556

P. Grassberger/Physica A 214 (1995) 547-559

2-D Ising model, Glauber dynamics, partially ordered start with M = 3/4

01 .

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

.

I

.

.

.

.

.

o,o1 0.001

0.9882

t

,

i

. . . . .

i

10

.

.

.

.

.

.

,,t

100 t

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

1000

Fig. 6. Decay of damage with time for completely random and independent initial configurations, at 7 different temperatures close to Tc.

4. Discussion One purpose of our work was to provide accurate estimates for z. In spite of a relatively modest amount of CPU time we believe that we have achieved this goal. While our value for d = 2 agrees with the most precise previous Monte Carlo estimate [6] (for earlier estimates, see the references in that paper), our value for d = 3 is somewhat lower than most previous MC estimates. In spite of this, our values are still slightly larger than those predicted by field theory [3], though not to any disturbing extent. Technically, this was made possible by several ingredients: ( 1) Using damage spreading with the heat bath algorithm to reduce statistical fluctuations. That this can reduce errors considerably was pointed out first in [21]. (2) Using multi-spin coding which is particularly useful for damage spreading where we want to use the same random numbers for more than one lattice. Efficiently providing independent random numbers for different lattices is just the main problem in multi-spin coding when N lattices (N= number of bits per word) are simulated simultaneously. (3) Starting with configurations far from equilibrium. In this way we not only reduce CPU time, but we also eliminate a possible source of systematic error. The latter would arise if we wanted to start with equilibrium configurations, but would not wait long enough. Technically, using far from equilibrium initial configurations is made possible by the observation of [ 12] that this does not affect z.

P. Grassberger/Physica A 214 (1995) 547-559

557

(4) Observing directly the average distance by which the damage spreads. In most previous simulations z was not measured directly, but only through the combination f l / v z which governs the relaxation of magnetization. When distances were measured directly [13,5,4,11,18], it was usually through more complicated variables which in some cases involved also logarithmic corrections [ 18]. It it hard to say which of the above points was the most important. All of them were used in one way or another in previous analyses, but it seems that they were never combined together. Indeed, our estimate of z is so precise that measurements of f l / v z through the decay of magnetization from an initially ordered state could give, together with Eqs. (9,10), the most reliable estimates of the static exponent f l / v and of Tc. Such measurements should also be done with the damaging technique, i.e. by measuring the magnetization difference between two lattices running with the same random number sequence. A somewhat different method for obtaining static exponents entirely from simulations of systems which are far from equilibrium was recently proposed in [33]. We were helped by the lucky accident that damage spreading is very fast when the original damage is on an entire hyperplane, with random configuration off this hyperplane. This is demonstrated by the fact that the increase of damage for this initial condition is as fast as for damages concentrated on lower dimensional subspaces - in spite of the fact that it starts off already with a much higher level. It is also demonstrated by simulations (not mentioned above) where the undamaged parts of the lattice were ordered. In that case, we observed much slower initial spreading and larger fluctuations. This is in agreement with analogous observations by Stauffer [ 18]. One of our most surprising results is that the increase of the amount of damage is governed by the same exponents O', independent of the dimensionality of the initially damaged region - provided just that it was not the entire lattice. For an explanation, let us consider the probability P ( t ) that a point damage is not yet healed after t time steps. We assume that it scales with some power of t, i.e. P ( t ) ~ t -p. T h e cross over time t* (see Eq. ( 3 ) ) can also be estimated as that time where the distance between the not yet healed damage clusters is of the same order as the correlation length sc N tl/z. With the above ansatz for P ( t ) this gives

p = id-

3/v)/z

- 0'~

0.9 1.1

(d=2) (d=3).

(15)

The density of damaged sites within a not yet healed cluster decays then as ~ t -/z/"z. Thus the damage clusters originating from single sites have rather small chances of survival, but those which do survive have the same magnetization density as would have a lattice which was completely ordered at t = 0. Now let us look at damage from hyperplanes. If M ~ t o' and R t ~ t 1/z , the average density of damaged sites in a layer of width 2Rt around the hyperplane is m ~ t ° ' - I / z . Both in 2 and 3 dimensions this exponent is negative and smaller than f l / v z . Thus the damage from a hyperplane is not distributed uniformly, but breaks up into isolated patches. These can grow into all

P. Grassberger/Physica A 214 (1995) 547-559

558

possible directions, invalidating our argument that the damage growth should be hindered when the original damage is too large. Notice that this argument can also be applied to other growth phenomena. For percolation, in particular, it had predicted correctly that the number of growth sites from an entire wetted edge of a cube increases as fast as if the growth had started from a single point on the edge [25]. After completing this work I learned that a similar argument has been used before in [34].

Acknowledgements I am indebted to N. Jan, H.K. Janssen, U. Ritschel and D. Stauffer for correspondence and discussions which I found extremely helpful. This work was supported by Deutsche Forschunggemeinschaft, SFB 237.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [ 16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

B.I. Halperin, EC. Hohenbetg and S.-K. Ma, Phys. Rcv. B 10 (1974) 139. S. Wansleben and D.E Landau, Phys. Rev. B 43 (1991) 6006. R. Bausch, V. Dohm, H.K. Janssen and R.K.E Zia, Phys. Rcv. Lett. 47 (1981) 1837. K. Maclsaac and N. Jan, J. Phys. A 25 (1992) 2139. D.L. Hunter, L. de Arcangelis, R. Matz, P.H. Poole and N. Jan, Physica A 196 (1993) 188. N. Ito, Physica A 196 (1993) 591. C. Miinkel, D.W. Heermann, J. Adler, M. Gofman and D. Stauffer, Physica A 193 (1993) 540. M. Siegert and D. Stauffer, Physica A 208 (1994) 1. H.-O. Heuer, J. Phys. A 25 (1992) L567. N. Ito, Physica A 192 (1993) 604. R. Matz, D.L. Hunter and N. Jan, J. Stat. Phys. 74 (1994) 903. H.K. Janssen, B. Schaub and B. Schmittmann, Z. Phys. B 73 (1989) 539. S.C. Glotzer, P.H. Poole and N. Jan, J. Stat. Phys. 68 (1992) 895. M. Creutz, Ann. Phys. 167 (1986) 62. H.E. Stanley, D. Stauffer, J. KertEsz and H. Herrmann, Phys. Rev. Lett.59 (1987) 2326. N. Jan and L. de Arcangelis, in: Annual Reviews of Computational Physics, vol. 1, ed. D. Stauffer (World Scientific, Singapore 1994) p. 1. G. le CaEr, Physica A 159 (1989) 329. D. Stauffer, J. Phys. A 26 (1993) L599. Z.-B. Li, U. Ritschel and B. Zheng, preprint SI-94-11 (Univ. Siegen, 1994). P. Grassberger, to be published. A. Coniglio, L. de Areangelis, H.J. Herrmann and N. Jan, Europhys. Lett. 8 (1989) 315. T.M. Liggett, Interacting Particle Systems (Springer, New York, 1985). A.M. Mariz, H.J. Herrmann and L. de Areangelis, J. Stat. Phys. 59 (1990) 1043. H.W. Diehl, in: Phase Transitions and Critical Phenomena, voi. 10, C. Domb and J.L. Lebowitz, eds. (Academic Press, NY, 1986) p. 75. E Grassberger, J. Phys. A 26 (1993) 1023. C.E Balllie, R. Gupta, K.A. Hawick and C.S. Pawley, Phys. Rev. B 45 (1992) 10438. A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44 (1991) 5081. NAG Fortran manual mark 15, vol. 9 (1991).

P Grassberger/Physica A 214 (1995) 547-559

[29] [30] [31] [32] [33] [34]

A.M. Ferrenberg, D.P. Landau and Y.J. Wong, Phys. Rev. Lett. 69 (1992) 3382. P. Grassberger, Phys. Lett. A 181 (1993) 43. D.A. Huse, Phys. Rev. B 40 (1989) 304. K. Humayun and A.J. Bray, J. Phys. A 24 (1991) 1915. Z.B. Li, L. Schiilke and B. Zheng, Siegen preprint Si-94-13 (1994). A.J. Bray, K. Humayun and T.J. Newman, Phys. Rev. B 43 (1991) 3699.

559