- Email: [email protected]
Critical exponents for extended dynamical systems with simultaneous updating: the case of the Ising model
Physica D 168–169 (2002) 318–324
Critical exponents for extended dynamical systems with simultaneous updating: the case of the Ising model Gabriel Pérez∗ , Francisco Sastre, Rubén Medina Departamento de F´ısica Aplicada, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Unidad Mérida, Apartado Postal 73 “Cordemex”, 97310 Mérida, Yucatán, Mexico
Abstract The behavior of a simultaneously updated Metropolis simulation of the two-dimensional (2D) Ising model is explored. The configurations generated by such a simulation are not constructed to be equilibrium configurations, and may therefore display a different critical behavior, assuming a phase transition is present. This is relevant to the recently posed question of how the updating scheme affects the critical exponents for chaotic extended systems that show continuous phase transitions. It is found that the simultaneously updated algorithm drives the lattice, for all non-zero temperatures, into blinking maximum energy configurations, making it of little use. A small modification, given by a reduction in the acceptance probabilities for spin flips, restores the ability of the algorithm to generate ordered and disordered states, with a well-defined phase transition between them. The critical temperature becomes a function of the acceptance probability for energy lowering spin flips. The critical exponents still fall into the 2D Ising universality class. This happens even though the simultaneously updated Metropolis algorithm does not really simulate the Ising model in equilibrium. © 2002 Elsevier Science B.V. All rights reserved. PACS: 05.10.Ln; 05.45.Ra; 05.50.+q Keywords: Metropolis algorithm; Ising model; Phase transition
1. Introduction Continuous order–disorder phase transitions have been found recently in extended two-dimensional (2D) systems, with both chaotic [1] and stochastic [2] local dynamics, and diffusive coupling. In all the models looked up until now, the local dynamics is scalar and odd-symmetric, and so they were expected to fall in the 2D Ising universality class [1,3]. However, in most cases it has been found that the critical exponents are close, but not equal, to those of the 2D Ising model [2,4,5]. In practice, one finds that ν ≈ 0.9, ∗
Corresponding author. E-mail address: [email protected] (G. P´erez).
instead of the Ising value ν = 1, while the ratios β/ν and γ /ν remain at their Ising values of 1/8 and 7/4, respectively. It has been conjectured that this may be related to the fact that updating in all these systems is done simultaneously over the whole lattice, since they do go back to the 2D Ising universality class when the updating is changed to asynchronous [4]. This would mean then that the updating scheme is a relevant parameter in the vicinity of the critical point. It is therefore important to ask what happens to standard, well known statistical models if one tries to force them to do simultaneous updating. This proposition in the abstract may seem to make no sense, since equilibrium statistical mechanics is usually studied in terms of ensembles, with no need for a time flow. But for
0167-2789/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 0 2 ) 0 0 5 1 9 - 5
G. P´erez et al. / Physica D 168–169 (2002) 318–324
the purposes of this work, what we are thinking about is the particular simulation algorithms that introduce a “time” in the models, in the sense of an ordered sequence of iteration steps. Of these many algorithms we will concentrate in the Metropolis [6] single-spin-flip Monte Carlo algorithm for the 2D nearest-neighbor Ising model, which is simple and well known [7,8]. In order to give a faithful sampling of the equilibrium distribution, this algorithm requires that the updating be done asynchronously and one can find in the literature many examples where this is accomplished following the order of the lattice (“typewriter updating”), or choosing sites at random. One is also allowed to update any non-interacting part of the lattice at once; for square lattices, half of the lattice can then be updated simultaneously (“checkerboard updating”). By mistake, full simultaneous updating have been tried [9], but it was found soon enough that in this case the Metropolis algorithm drives square lattices into blinking maximum energy patterns, which we will discuss later on. It has been assumed then that simultaneous updating is of no use for the Metropolis algorithm [10]. The rest of this paper is arranged as follows. In Section 2 we will review the Metropolis algorithm, and will discuss why its simultaneous updated version always gets trapped into blinking configurations. We will then introduce a slight modification in the algorithm, so that it will once more show a continuous phase transition as the temperature is changed. In Section 3 we will show the results of the simulations done with this modified algorithm, and in Section 4 we will present our conclusions.
2. The Metropolis and simultaneous Metropolis algorithms for the Ising model The Metropolis algorithm for the simulation of physical systems in the canonical ensemble [6] has been the workhorse of numerical statistical mechanics for almost half a century. For completeness, we will give a short description of the algorithm here, following the very detailed presentation given in [7]. It consists of the construction of a sequence (more exactly, a Markov chain) of configurations, which
319
we will denote with greek indices µ, ν, etc., which appear with frequencies proportional to their Boltzmann weight, Wµ = exp(−βEµ ), where Eµ is the energy of the configuration µ and β ≡ (kB T )−1 . Equilibrium average values of observables are then given by simple averages over this sequence. The Markov chain of configurations one generates has to obey the following two restrictions: (1) it should be ergodic, i.e., for a finite system, one should be able to reach any valid configuration ν starting from any other valid configuration µ and (2) it should exhibit detailed balance, which means that the ratio between the probability of going from one configuration µ to another ν, and the probability for the opposite transition, ν to µ, should be proportional to the ratio of the corresponding Boltzmann weights. Since we will be dealing here only with its application to the Ising model, let us describe the Metropolis algorithm in terms of spin configurations. The steps for one iteration of the algorithm are the following: starting form some arbitrarily generated spin configuration (1), choose one spin. As mentioned in the introduction, this may be done in some order through the lattice, or at random, and one may even choose a group of spins, as long as they do not interact with each other. (2) One then flips the chosen spin and calculates the change in energy E this leads to. (3) If this
E is negative, the spin flip lowers the energy of the lattice, and is always accepted. If not (4), the flip increases the energy of the lattice, and is given a chance of happening proportional to the Boltzmann weight involved in the transition, W = exp(−β E). This is implemented with the generation of a random number η, uniformly distributed between 0 and 1, and accepting the spin flip if η < W . (5) Go back to step (1). It is easy to see the ergodicity of this algorithm for T = 0, since going from any configuration µ to any other configuration ν is always possible, requiring only the generation of a finite string of random numbers above certain thresholds. For the condition of detailed balance, one notices the following: let µ and ν be configurations that differ only in the orientation of one spin. The probability of going from µ to ν, denoted P (µ → ν), is composed of two parts: the probability of generating ν starting from µ, which is just 1/N ,
320
G. P´erez et al. / Physica D 168–169 (2002) 318–324
where N is the number of sites in the lattice, and the probability of accepting the new configuration, A(µ → ν). In terms of these quantities, the condition of detailed balance is P (µ → ν) A(µ → ν)/N = = exp(β(Eµ − Eν )). P (ν → µ) A(ν → µ)/N The Metropolis algorithm is then assigning the maximum possible acceptance probability A(µ → ν) = 1 for energy lowering transitions, and sets A(ν → µ) = exp(−β(Eµ − Eν )) otherwise, giving the correct acceptance ratio for the previous equation. As mentioned, this needs to be done one spin at a time. When the previously described steps are applied on every site of a square lattice simultaneously, one finds that in a few iterations the lattice reaches a state of maximum energy (for ferromagnetic coupling, neighboring spins pointing opposite to each other), and that all spins flip at every iteration. This effect has been called “blinking checkerboard catastrophe” [9], but is not limited to square lattices. In fact, when simultaneous updating is attempted in a triangular lattice, the configuration soon evolves into one of the many local maxima of energy that the frustration of the system offers, and then continues flipping the whole lattice at every iteration. It is not difficult to understand what is actually happening in these cases. Suppose a checkerboard cluster develops by chance. Then, for every site in this cluster, one finds that energy can be lowered flipping the local spin. If one then precedes to do just that, the whole cluster is flipped and therefore goes back in two iterations to the original situation. The same applies to sites in touch with the cluster, which then grows until covering the whole lattice.1 The end reason for this behavior is the fact that spin flips that lower the energy of the system are always accepted. The solution for this anomaly becomes then clear: allow the algorithm to reject a given fraction of the energy lowering spin flips. At the same time, in order to maintain detailed balance in its asynchronous version, a similar extra fraction of the energy-raising 1 The possibility for the formation of large checkerboard-like domains with static walls between them cannot be discarded. A situation like this has not appeared in our simulations.
spin flips should be rejected. The spin-flip acceptance is fixed as A(µ → ν) = A0 , for Eµ > Eν , and A(µ → ν) = A0 exp(β(Eµ − Eν )) otherwise, where A0 is some constant with 0 < A0 < 1. For the asynchronously updated algorithm this change means that a fraction 1 − A0 of all attempted moves are always discarded, and the algorithm is simply a less efficient version of the original Metropolis. Before getting into the results, it is worth to restate that this algorithm is no longer simulating the Ising model in equilibrium. For what we want to study here, it could be considered simply as an stochastic extended dynamical system, which borrows most of its structure from the Metropolis simulation of Ising model, but is not really intended to be a simulation of any physical model in equilibrium. In particular, notice that there is no clear way of establishing detailed balance here. On the other hand, it should be noticed that our evaluation of critical quantities is based in finite size scaling [11], which strictly speaking is justified only for equilibrium models. The motivation for its use here is mostly heuristic.
3. Numerical results for the Ising model simultaneous Metropolis algorithm The simultaneous Metropolis algorithm was implemented in a simulation of the Ising model in a 2D
Fig. 1. Critical temperature as a function of the acceptance factor A0 . The line is provided only as a guide to the eye.
G. P´erez et al. / Physica D 168–169 (2002) 318–324
321
A detailed simulation of the model for A0 = 1/2 was performed, using sides ranging from L = 42 to 128. Periodic boundary conditions were implemented. For each size, the total number of iterations was around 1400 times the correlation time. A dynamical exponent z ≈ 2 was found. First we located the critical point, using again fourth-order cumulants. The results are given in Fig. 2, and out estimation for this parameter is Tc (A0 = 21 ) = 3.7475(10). Fig. 2. Crossing of fourth-order cumulants UL (T ), used to locate the critical point. The actual data contains 20 points for each value of L, covering a larger range. The values of L are, starting from the smallest slope, L = 42, 48, 56, 64, 74, 84, 98, 112 and 128. The lines are second-order polynomial fits to the full data, and give the best polynomial fit in the sense of minimizing χ 2 per degree of freedom.
Here, the digits between parenthesis indicate the uncertainty in the corresponding last digits of the quantity. Using this value, we evaluated the critical
triangular lattice, and order–disorder continuous phase transitions were found in a large range of A0 . All the determinations of critical quantities were done using arguments from finite size scaling. First, we did a preliminary exploration of the phase diagram, and the results are shown in Fig. 1. The critical temperature becomes a function of A0 , and for the values we have checked stays above the correct value in triangular lattices, Tc = 4/ ln (3) ≈ 3.641. This determination was done via crossing of fourth-order cumulants [12] UL (T ) = 1 −
m4L (T )
3(m2L (T ))2
,
where mkL (T ) is the average value found for the kth power of the magnetization density.2 The two end points are not included, since, as already mentioned, A0 = 1 gives the blinking maximum energy patterns, and A0 = 0 would leave any configuration unmodified. The behavior of the model with A0 close to these two extremes is left for future study. 2 As usual when dealing with finite size lattices, it is necessary to consider the absolute value of the instantaneous magnetization, since otherwise the average values of all its odd moments would go to zero for large runs. Since this systems is not in the thermodynamic limit, it cannot display ergodicity breaking.
Fig. 3. (a) Scaling of derivatives of cumulants U and other quantities used to calculate ν. Starting from the lower line, the quantity A corresponds to ∂T UL (T ), ∂T ln ( m ), ∂T ln ( m2 ), ∂T ln ( m3 ) and ∂T ln ( m4 ). The lines have slope 1. (b) Same quantities, after removing the leading term L1/ν = L from their scaling.
322
G. P´erez et al. / Physica D 168–169 (2002) 318–324
exponent ν using the scaling behavior at the critical point of the following quantities: ∂T UL (Tc ) ∝ L1/ν ,
∂T ln mkL (Tc ) ∝ L1/ν .
For the uncertainty levels we have, we may assume that we are already in the scaling regime. Therefore, these equations do not include finite size corrections. The results are shown in Fig. 3(a) and (b), which show the scaling of the derivatives for U and for ln (mk ), with k = 1, 2, 3 and 4. The combined result obtained from these five lines is ν = 1.0070(81) perfectly compatible with the exact Ising model result ν = 1. The magnetization exponent β was calculated
Fig. 4. (a) Scaling of the first four moments of the magnetization density. Starting from the lower line, we have m, m2 , m3 and m4 . The lines have slope −1/8. (b) First four moments of the magnetization density, after removing the leading term L−kβ/ν = L−k/8 from their scaling.
using the scaling relations mkL (Tc ) ∝ L−kβ/ν , where again we run k from 1 to 4, and no finite size corrections are used. The results are given in Fig. 4(a) and (b), and the combined value we find is β = 0.1290(30). ν The Ising value β/ν = 1/8 differs slightly from this result, but is still within its 2σ range. Next, we checked the magnetic susceptibility exponent, using χL (T ) = N (m2L (T ) − (mL (T ))2 ), and its scaling behavior χL (T ) ∝ Lγ /ν .
Fig. 5. (a) Scaling of the specific heat, assuming a power-law dependence. The line corresponds to the best linear fit with slope 0.2232. (b) Scaling of the specific heat, after removing the leading power law term. It is clear that a simple power law does not fit well with the specific heat.
G. P´erez et al. / Physica D 168–169 (2002) 318–324
From the results obtained (not shown) the value found is γ = 1.740(21) ν again, compatible with the Ising value γ /ν = 7/4. Finally, one can explore the behavior of the spe2 (T ) − (e (T ))2 ), where cific heat, given by c = N (eL L k eL (T ) is the average value found for the kth power of the energy density. For the Ising model, this quantity shows a logarithmic divergence as the lattice size grows, which is interpreted as a zero value for the exponent α in a power-law divergence. For the simultaneous Metropolis simulation, and assuming that the scaling relation α + 2β + γ = 2 remains valid, the values already found for ν, β/ν and γ /ν imply that
Fig. 6. (a) Scaling of the specific heat, assuming a logarithmic dependence. The line corresponds to the best linear fit, with slope c0 = 9.88. (b) Scaling of the specific heat after removing the logarithmic term.
323
α ≈ −0.012(26), compatible with the Ising result. The actual scaling found shows that in fact c diverges as ln (L) (see Fig. 5(a) and (b)), and cannot be fitted to a power law (Fig. 6(a) and (b)).
4. Conclusions and final comments As mentioned in Section 1, the phase transitions that appear in diffusively coupled chaotic 2D lattices with scalar order parameter do not fall into the 2D Ising universality class, as was originally expected. Since they do become Ising-like when the updating is changed to asynchronous, the only explanation offered up to now has been to give to the updating scheme the role of a relevant parameter. In order to test this possibility, we have explored a simultaneous version of the well-known Metropolis algorithm for the simulation of the 2D Ising model. It is found that the change from asynchronous to simultaneous updating in this algorithm always drives the lattice into a blinking maximum energy configuration, which are really not relevant to the properties one expects the model to show, even with such a changed simulation. This anomalous result can be avoided via the reduction of the acceptance probability A0 for energy lowering spin flips. This reduction should then be also applied to energy raising flips, in order to preserve detailed balance in the asynchronous version of the algorithm. This modification generated one-parameter family of algorithms, indexed by A0 . For A0 = 1/2, the simultaneous Metropolis simulation of the 2D Ising model gives the same critical behavior of the original model, even if the simulation itself is not designed to explore its equilibrium configurations. This difference with a correct simulation does show in a different critical temperature, which actually becomes a function of A0 . One then has the following main conclusion: for the universality class of the Ising model in the Metropolis simulation, the issue of synchronicity in its updating is irrelevant. It should be mentioned, however, that it could be possible that the reduction in acceptance in energy lowering spin flips may be acting as an effective source of asynchrony in the updating [13]. We believe that
324
G. P´erez et al. / Physica D 168–169 (2002) 318–324
the possible asynchrony thus generated is not enough as to explain the present results, but the possibility does merit some extra checking. The simplest way of doing so would be to raise A0 to some value close to 1. Another option, which is now being studied, is to consider the heat-bath algorithm, which should not need the introduction of an acceptance lowering factor. One is then left with the same puzzle: What makes the Miller–Huse [1] and some other chaotic models [4], and a similar stochastic model [2], show a critical behavior different from that of the Ising model? At the moment, two other options may need to be considered: (1) presence of still not well understood corrections, or (2) existence in parameter space of some close but not yet identified critical point, which would generate cross-over behavior in this dynamical models [13]. These are questions that certainly need to be addressed.
Acknowledgements This work was supported by CONACyT (Mexico) through Grant No. 28383E.
References [1] J. Miller, D. Huse, Phys. Rev. E 48 (1993) 2528. [2] F. Sastre, G. Pérez, Phys. Rev. E. [3] G. Grinstein, C. Jayaprakash, Y. He, Phys. Rev. Lett. 55 (1985) 2527; T. Bohr, G. Grinstein, Y. He, C. Jayaprakash, Phys. Rev. Lett. 58 (1987) 2155. [4] P. Marcq, H. Chaté, P. Manneville, Phys. Rev. Lett. 77 (1996) 4003; P. Marcq, H. Chaté, P. Manneville, Phys. Rev. E 55 (1997) 2606. [5] D. Makowiec, Phys. Rev. E 60 (1999) 3787. [6] N. Metropolis, et al., J. Chem. Phys. 21 (1953) 1087. [7] M.E.J. Newman, G.T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford University Press, Oxford, 1999. [8] D.P. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press, Cambridge, 2000. [9] B. Hayes, Am. Scientist 88 (2000) 384. [10] G. Vichniac, Physica D 10 (1984) 96. [11] V. Privman (Ed.), Finite-size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990. [12] K. Binder, Z. Phys. B 43 (1982) 119; K. Binder, D. Stauffer, in: K. Binder (Ed.), A Simple Introduction to Monte Carlo Simulation and Some Specialized Topics, Applications of the Monte Carlo Method in Statistical Physics, Topics in Current Physics, Vol. 36, Springer, New York, 1984. [13] H. Chaté, Private communication.
Critical exponents for extended dynamical systems with simultaneous updating: the case of the Ising model Gabriel Pérez∗ , Francisco Sastre, Rubén Medina Departamento de F´ısica Aplicada, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Unidad Mérida, Apartado Postal 73 “Cordemex”, 97310 Mérida, Yucatán, Mexico
Abstract The behavior of a simultaneously updated Metropolis simulation of the two-dimensional (2D) Ising model is explored. The configurations generated by such a simulation are not constructed to be equilibrium configurations, and may therefore display a different critical behavior, assuming a phase transition is present. This is relevant to the recently posed question of how the updating scheme affects the critical exponents for chaotic extended systems that show continuous phase transitions. It is found that the simultaneously updated algorithm drives the lattice, for all non-zero temperatures, into blinking maximum energy configurations, making it of little use. A small modification, given by a reduction in the acceptance probabilities for spin flips, restores the ability of the algorithm to generate ordered and disordered states, with a well-defined phase transition between them. The critical temperature becomes a function of the acceptance probability for energy lowering spin flips. The critical exponents still fall into the 2D Ising universality class. This happens even though the simultaneously updated Metropolis algorithm does not really simulate the Ising model in equilibrium. © 2002 Elsevier Science B.V. All rights reserved. PACS: 05.10.Ln; 05.45.Ra; 05.50.+q Keywords: Metropolis algorithm; Ising model; Phase transition
1. Introduction Continuous order–disorder phase transitions have been found recently in extended two-dimensional (2D) systems, with both chaotic [1] and stochastic [2] local dynamics, and diffusive coupling. In all the models looked up until now, the local dynamics is scalar and odd-symmetric, and so they were expected to fall in the 2D Ising universality class [1,3]. However, in most cases it has been found that the critical exponents are close, but not equal, to those of the 2D Ising model [2,4,5]. In practice, one finds that ν ≈ 0.9, ∗
Corresponding author. E-mail address: [email protected] (G. P´erez).
instead of the Ising value ν = 1, while the ratios β/ν and γ /ν remain at their Ising values of 1/8 and 7/4, respectively. It has been conjectured that this may be related to the fact that updating in all these systems is done simultaneously over the whole lattice, since they do go back to the 2D Ising universality class when the updating is changed to asynchronous [4]. This would mean then that the updating scheme is a relevant parameter in the vicinity of the critical point. It is therefore important to ask what happens to standard, well known statistical models if one tries to force them to do simultaneous updating. This proposition in the abstract may seem to make no sense, since equilibrium statistical mechanics is usually studied in terms of ensembles, with no need for a time flow. But for
0167-2789/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 0 2 ) 0 0 5 1 9 - 5
G. P´erez et al. / Physica D 168–169 (2002) 318–324
the purposes of this work, what we are thinking about is the particular simulation algorithms that introduce a “time” in the models, in the sense of an ordered sequence of iteration steps. Of these many algorithms we will concentrate in the Metropolis [6] single-spin-flip Monte Carlo algorithm for the 2D nearest-neighbor Ising model, which is simple and well known [7,8]. In order to give a faithful sampling of the equilibrium distribution, this algorithm requires that the updating be done asynchronously and one can find in the literature many examples where this is accomplished following the order of the lattice (“typewriter updating”), or choosing sites at random. One is also allowed to update any non-interacting part of the lattice at once; for square lattices, half of the lattice can then be updated simultaneously (“checkerboard updating”). By mistake, full simultaneous updating have been tried [9], but it was found soon enough that in this case the Metropolis algorithm drives square lattices into blinking maximum energy patterns, which we will discuss later on. It has been assumed then that simultaneous updating is of no use for the Metropolis algorithm [10]. The rest of this paper is arranged as follows. In Section 2 we will review the Metropolis algorithm, and will discuss why its simultaneous updated version always gets trapped into blinking configurations. We will then introduce a slight modification in the algorithm, so that it will once more show a continuous phase transition as the temperature is changed. In Section 3 we will show the results of the simulations done with this modified algorithm, and in Section 4 we will present our conclusions.
2. The Metropolis and simultaneous Metropolis algorithms for the Ising model The Metropolis algorithm for the simulation of physical systems in the canonical ensemble [6] has been the workhorse of numerical statistical mechanics for almost half a century. For completeness, we will give a short description of the algorithm here, following the very detailed presentation given in [7]. It consists of the construction of a sequence (more exactly, a Markov chain) of configurations, which
319
we will denote with greek indices µ, ν, etc., which appear with frequencies proportional to their Boltzmann weight, Wµ = exp(−βEµ ), where Eµ is the energy of the configuration µ and β ≡ (kB T )−1 . Equilibrium average values of observables are then given by simple averages over this sequence. The Markov chain of configurations one generates has to obey the following two restrictions: (1) it should be ergodic, i.e., for a finite system, one should be able to reach any valid configuration ν starting from any other valid configuration µ and (2) it should exhibit detailed balance, which means that the ratio between the probability of going from one configuration µ to another ν, and the probability for the opposite transition, ν to µ, should be proportional to the ratio of the corresponding Boltzmann weights. Since we will be dealing here only with its application to the Ising model, let us describe the Metropolis algorithm in terms of spin configurations. The steps for one iteration of the algorithm are the following: starting form some arbitrarily generated spin configuration (1), choose one spin. As mentioned in the introduction, this may be done in some order through the lattice, or at random, and one may even choose a group of spins, as long as they do not interact with each other. (2) One then flips the chosen spin and calculates the change in energy E this leads to. (3) If this
E is negative, the spin flip lowers the energy of the lattice, and is always accepted. If not (4), the flip increases the energy of the lattice, and is given a chance of happening proportional to the Boltzmann weight involved in the transition, W = exp(−β E). This is implemented with the generation of a random number η, uniformly distributed between 0 and 1, and accepting the spin flip if η < W . (5) Go back to step (1). It is easy to see the ergodicity of this algorithm for T = 0, since going from any configuration µ to any other configuration ν is always possible, requiring only the generation of a finite string of random numbers above certain thresholds. For the condition of detailed balance, one notices the following: let µ and ν be configurations that differ only in the orientation of one spin. The probability of going from µ to ν, denoted P (µ → ν), is composed of two parts: the probability of generating ν starting from µ, which is just 1/N ,
320
G. P´erez et al. / Physica D 168–169 (2002) 318–324
where N is the number of sites in the lattice, and the probability of accepting the new configuration, A(µ → ν). In terms of these quantities, the condition of detailed balance is P (µ → ν) A(µ → ν)/N = = exp(β(Eµ − Eν )). P (ν → µ) A(ν → µ)/N The Metropolis algorithm is then assigning the maximum possible acceptance probability A(µ → ν) = 1 for energy lowering transitions, and sets A(ν → µ) = exp(−β(Eµ − Eν )) otherwise, giving the correct acceptance ratio for the previous equation. As mentioned, this needs to be done one spin at a time. When the previously described steps are applied on every site of a square lattice simultaneously, one finds that in a few iterations the lattice reaches a state of maximum energy (for ferromagnetic coupling, neighboring spins pointing opposite to each other), and that all spins flip at every iteration. This effect has been called “blinking checkerboard catastrophe” [9], but is not limited to square lattices. In fact, when simultaneous updating is attempted in a triangular lattice, the configuration soon evolves into one of the many local maxima of energy that the frustration of the system offers, and then continues flipping the whole lattice at every iteration. It is not difficult to understand what is actually happening in these cases. Suppose a checkerboard cluster develops by chance. Then, for every site in this cluster, one finds that energy can be lowered flipping the local spin. If one then precedes to do just that, the whole cluster is flipped and therefore goes back in two iterations to the original situation. The same applies to sites in touch with the cluster, which then grows until covering the whole lattice.1 The end reason for this behavior is the fact that spin flips that lower the energy of the system are always accepted. The solution for this anomaly becomes then clear: allow the algorithm to reject a given fraction of the energy lowering spin flips. At the same time, in order to maintain detailed balance in its asynchronous version, a similar extra fraction of the energy-raising 1 The possibility for the formation of large checkerboard-like domains with static walls between them cannot be discarded. A situation like this has not appeared in our simulations.
spin flips should be rejected. The spin-flip acceptance is fixed as A(µ → ν) = A0 , for Eµ > Eν , and A(µ → ν) = A0 exp(β(Eµ − Eν )) otherwise, where A0 is some constant with 0 < A0 < 1. For the asynchronously updated algorithm this change means that a fraction 1 − A0 of all attempted moves are always discarded, and the algorithm is simply a less efficient version of the original Metropolis. Before getting into the results, it is worth to restate that this algorithm is no longer simulating the Ising model in equilibrium. For what we want to study here, it could be considered simply as an stochastic extended dynamical system, which borrows most of its structure from the Metropolis simulation of Ising model, but is not really intended to be a simulation of any physical model in equilibrium. In particular, notice that there is no clear way of establishing detailed balance here. On the other hand, it should be noticed that our evaluation of critical quantities is based in finite size scaling [11], which strictly speaking is justified only for equilibrium models. The motivation for its use here is mostly heuristic.
3. Numerical results for the Ising model simultaneous Metropolis algorithm The simultaneous Metropolis algorithm was implemented in a simulation of the Ising model in a 2D
Fig. 1. Critical temperature as a function of the acceptance factor A0 . The line is provided only as a guide to the eye.
G. P´erez et al. / Physica D 168–169 (2002) 318–324
321
A detailed simulation of the model for A0 = 1/2 was performed, using sides ranging from L = 42 to 128. Periodic boundary conditions were implemented. For each size, the total number of iterations was around 1400 times the correlation time. A dynamical exponent z ≈ 2 was found. First we located the critical point, using again fourth-order cumulants. The results are given in Fig. 2, and out estimation for this parameter is Tc (A0 = 21 ) = 3.7475(10). Fig. 2. Crossing of fourth-order cumulants UL (T ), used to locate the critical point. The actual data contains 20 points for each value of L, covering a larger range. The values of L are, starting from the smallest slope, L = 42, 48, 56, 64, 74, 84, 98, 112 and 128. The lines are second-order polynomial fits to the full data, and give the best polynomial fit in the sense of minimizing χ 2 per degree of freedom.
Here, the digits between parenthesis indicate the uncertainty in the corresponding last digits of the quantity. Using this value, we evaluated the critical
triangular lattice, and order–disorder continuous phase transitions were found in a large range of A0 . All the determinations of critical quantities were done using arguments from finite size scaling. First, we did a preliminary exploration of the phase diagram, and the results are shown in Fig. 1. The critical temperature becomes a function of A0 , and for the values we have checked stays above the correct value in triangular lattices, Tc = 4/ ln (3) ≈ 3.641. This determination was done via crossing of fourth-order cumulants [12] UL (T ) = 1 −
m4L (T )
3(m2L (T ))2
,
where mkL (T ) is the average value found for the kth power of the magnetization density.2 The two end points are not included, since, as already mentioned, A0 = 1 gives the blinking maximum energy patterns, and A0 = 0 would leave any configuration unmodified. The behavior of the model with A0 close to these two extremes is left for future study. 2 As usual when dealing with finite size lattices, it is necessary to consider the absolute value of the instantaneous magnetization, since otherwise the average values of all its odd moments would go to zero for large runs. Since this systems is not in the thermodynamic limit, it cannot display ergodicity breaking.
Fig. 3. (a) Scaling of derivatives of cumulants U and other quantities used to calculate ν. Starting from the lower line, the quantity A corresponds to ∂T UL (T ), ∂T ln ( m ), ∂T ln ( m2 ), ∂T ln ( m3 ) and ∂T ln ( m4 ). The lines have slope 1. (b) Same quantities, after removing the leading term L1/ν = L from their scaling.
322
G. P´erez et al. / Physica D 168–169 (2002) 318–324
exponent ν using the scaling behavior at the critical point of the following quantities: ∂T UL (Tc ) ∝ L1/ν ,
∂T ln mkL (Tc ) ∝ L1/ν .
For the uncertainty levels we have, we may assume that we are already in the scaling regime. Therefore, these equations do not include finite size corrections. The results are shown in Fig. 3(a) and (b), which show the scaling of the derivatives for U and for ln (mk ), with k = 1, 2, 3 and 4. The combined result obtained from these five lines is ν = 1.0070(81) perfectly compatible with the exact Ising model result ν = 1. The magnetization exponent β was calculated
Fig. 4. (a) Scaling of the first four moments of the magnetization density. Starting from the lower line, we have m, m2 , m3 and m4 . The lines have slope −1/8. (b) First four moments of the magnetization density, after removing the leading term L−kβ/ν = L−k/8 from their scaling.
using the scaling relations mkL (Tc ) ∝ L−kβ/ν , where again we run k from 1 to 4, and no finite size corrections are used. The results are given in Fig. 4(a) and (b), and the combined value we find is β = 0.1290(30). ν The Ising value β/ν = 1/8 differs slightly from this result, but is still within its 2σ range. Next, we checked the magnetic susceptibility exponent, using χL (T ) = N (m2L (T ) − (mL (T ))2 ), and its scaling behavior χL (T ) ∝ Lγ /ν .
Fig. 5. (a) Scaling of the specific heat, assuming a power-law dependence. The line corresponds to the best linear fit with slope 0.2232. (b) Scaling of the specific heat, after removing the leading power law term. It is clear that a simple power law does not fit well with the specific heat.
G. P´erez et al. / Physica D 168–169 (2002) 318–324
From the results obtained (not shown) the value found is γ = 1.740(21) ν again, compatible with the Ising value γ /ν = 7/4. Finally, one can explore the behavior of the spe2 (T ) − (e (T ))2 ), where cific heat, given by c = N (eL L k eL (T ) is the average value found for the kth power of the energy density. For the Ising model, this quantity shows a logarithmic divergence as the lattice size grows, which is interpreted as a zero value for the exponent α in a power-law divergence. For the simultaneous Metropolis simulation, and assuming that the scaling relation α + 2β + γ = 2 remains valid, the values already found for ν, β/ν and γ /ν imply that
Fig. 6. (a) Scaling of the specific heat, assuming a logarithmic dependence. The line corresponds to the best linear fit, with slope c0 = 9.88. (b) Scaling of the specific heat after removing the logarithmic term.
323
α ≈ −0.012(26), compatible with the Ising result. The actual scaling found shows that in fact c diverges as ln (L) (see Fig. 5(a) and (b)), and cannot be fitted to a power law (Fig. 6(a) and (b)).
4. Conclusions and final comments As mentioned in Section 1, the phase transitions that appear in diffusively coupled chaotic 2D lattices with scalar order parameter do not fall into the 2D Ising universality class, as was originally expected. Since they do become Ising-like when the updating is changed to asynchronous, the only explanation offered up to now has been to give to the updating scheme the role of a relevant parameter. In order to test this possibility, we have explored a simultaneous version of the well-known Metropolis algorithm for the simulation of the 2D Ising model. It is found that the change from asynchronous to simultaneous updating in this algorithm always drives the lattice into a blinking maximum energy configuration, which are really not relevant to the properties one expects the model to show, even with such a changed simulation. This anomalous result can be avoided via the reduction of the acceptance probability A0 for energy lowering spin flips. This reduction should then be also applied to energy raising flips, in order to preserve detailed balance in the asynchronous version of the algorithm. This modification generated one-parameter family of algorithms, indexed by A0 . For A0 = 1/2, the simultaneous Metropolis simulation of the 2D Ising model gives the same critical behavior of the original model, even if the simulation itself is not designed to explore its equilibrium configurations. This difference with a correct simulation does show in a different critical temperature, which actually becomes a function of A0 . One then has the following main conclusion: for the universality class of the Ising model in the Metropolis simulation, the issue of synchronicity in its updating is irrelevant. It should be mentioned, however, that it could be possible that the reduction in acceptance in energy lowering spin flips may be acting as an effective source of asynchrony in the updating [13]. We believe that
324
G. P´erez et al. / Physica D 168–169 (2002) 318–324
the possible asynchrony thus generated is not enough as to explain the present results, but the possibility does merit some extra checking. The simplest way of doing so would be to raise A0 to some value close to 1. Another option, which is now being studied, is to consider the heat-bath algorithm, which should not need the introduction of an acceptance lowering factor. One is then left with the same puzzle: What makes the Miller–Huse [1] and some other chaotic models [4], and a similar stochastic model [2], show a critical behavior different from that of the Ising model? At the moment, two other options may need to be considered: (1) presence of still not well understood corrections, or (2) existence in parameter space of some close but not yet identified critical point, which would generate cross-over behavior in this dynamical models [13]. These are questions that certainly need to be addressed.
Acknowledgements This work was supported by CONACyT (Mexico) through Grant No. 28383E.
References [1] J. Miller, D. Huse, Phys. Rev. E 48 (1993) 2528. [2] F. Sastre, G. Pérez, Phys. Rev. E. [3] G. Grinstein, C. Jayaprakash, Y. He, Phys. Rev. Lett. 55 (1985) 2527; T. Bohr, G. Grinstein, Y. He, C. Jayaprakash, Phys. Rev. Lett. 58 (1987) 2155. [4] P. Marcq, H. Chaté, P. Manneville, Phys. Rev. Lett. 77 (1996) 4003; P. Marcq, H. Chaté, P. Manneville, Phys. Rev. E 55 (1997) 2606. [5] D. Makowiec, Phys. Rev. E 60 (1999) 3787. [6] N. Metropolis, et al., J. Chem. Phys. 21 (1953) 1087. [7] M.E.J. Newman, G.T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford University Press, Oxford, 1999. [8] D.P. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press, Cambridge, 2000. [9] B. Hayes, Am. Scientist 88 (2000) 384. [10] G. Vichniac, Physica D 10 (1984) 96. [11] V. Privman (Ed.), Finite-size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990. [12] K. Binder, Z. Phys. B 43 (1982) 119; K. Binder, D. Stauffer, in: K. Binder (Ed.), A Simple Introduction to Monte Carlo Simulation and Some Specialized Topics, Applications of the Monte Carlo Method in Statistical Physics, Topics in Current Physics, Vol. 36, Springer, New York, 1984. [13] H. Chaté, Private communication.