Stochastic dynamics of coupled systems and damage spreading

Stochastic dynamics of coupled systems and damage spreading

T. Tomé^I; E. Arashiro^II; J. R. Drugowich de Felício^II; M. J. de Oliveira^I

^IInstituto de Física, Universidade de São Paulo, Caixa Postal 66318, 05315-970 São Paulo, São Paulo, Brazil

^IIDepartamento de Física e Matemática, FFCLRP, Universidade de São Paulo, Av. Bandeirantes, 3900, 014040-901 Ribeirão Preto, São Paulo, Brazil

ABSTRACT

We study the damage spreading in the one-dimensional Ising model by means of the stochastic dynamics resulting from coupling the system and its replica by a family of algorithms that interpolate between the heat bath and the Hinrichsen-Domany algorithms. At high temperatures the dynamics is exactly mapped to the Domany-Kinzel probabilistic cellular automaton. Using a mean-field approximation and Monte Carlo simulations we find the critical line that separates the phase where the damage spreads from the one where it does not.

I Introduction

It is well known that most models studied in equilibrium statistical mechanics, such as the Ising model, are defined in a static way through the equilibrium Gibbs probability distribution associated with the Hamiltonian of the model. It is desirable from the theoretical and numerical point of view to assign a dynamics to such models. The stochastic dynamics introduced by Glauber [1] is the prototype example of a dynamics assigned to a static-defined model. The numerous versions of the Monte Carlo method [2], used in statistical mechanics are also examples of dynamics assigned to static-defined models. All of them are markovian processes that have the Gibbs probability distribution as the stationary distribution. In general they are either continuous time process governed by a master equation [3-8] or probabilistic cellular automata [8-13] The latter is defined by a stochastic matrix, whose elements are the transition probabilities, and the former by the evolution matrix, whose nondiagonal elements are the transition rates.

If we wish to simulate, for instance, the Ising model we have to choose one of the possible stochastic dynamics, since there are many. Having decided which dynamics to use, that is, having decided which probabilistic rules to use, we realize that there are several ways of doing the actual simulation corresponding to the chosen probabilistic rules. For instance, for the case of the probabilistic cellular automaton used by Derrida and Weisbuch [10] to simulate the Ising model, and which will concern us here, there are several ways of realizing the dynamics. We may use the so called heat-bath algorithm [14] or the algorithm introduced more recently by Hinrichsen and Domany [15] or any other we may invent. These algorithms govern the movement of the system in phase space and they may be called stochastic equations of motion in phase space. Different algorithms may be the realization of the same probabilistic rule or stochastic dynamics.

The description of a system either by the equation of motion or by the time evolution of the probability are equivalent. An analogy can be made with the Brownian motion which can be described either by the Langevin equation or by its associated Fokker-Planck equation [4,8,16,17]. The first is a stochastic equation of motion of a representative point in phase space whereas the second governs the time evolution of the probability distribution in phase space.

In the study of damage spreading [10,12,15,18-24] it has been realized that algorithms that are realization of the same probalistic rules may yield different results for the damage spreading [21, 15, 22, 24], and they usually do. The damage spreading is a procedure through which we may study the sensibility of the time evolution of systems with respect to the initial conditions. The procedure amounts to coupling the system with a replica of it, each of them following the same equation of motion. The coupling is acomplished by the use of the same sequence of random numbers. The equation of motion for each system together with the use of the same random number defines the equation of motion for the coupled system from which we obtain the joint transition probability [12, 13] for the coupled system.

Suppose one uses an algorithm to couple a system and its replica. This will lead us to a certain joint transition probability. If another algorithm is used, which is also a realization of the same transition probability for a single system, the joint transition probability will be distinct. The correlation between system and replica will also be distinct and, in particular, the Hamming distance which is a measure of the damage spreading will be different. For example, in the one-dimensional Ising model, the heat-bath (HB) algorithm [14] will give no damage spreading whereas the Hinrichsen-Domany (HD) algorithm [15] will exhibit a spreading of damage above a certain temperature [15]. This is an impressive example that the damage spreading is not an intrinsic static property of a given system, but depends on the algorithm, or the stochastic equation of motion, we use to perform the actual simulation [15, 22].

In this paper we introduce a family of algorithms, or equations of motion, spanned by a parameter that interpolates between the HB and HD algorithms. The associated transition probability corresponds, for all values of the parameter, to the Derrida-Weisbush (DW) probabilistic cellular automaton [10]. If we use this family of algorithms to study the spreading of damage, as we will do here, the parameter will have no effect on each system separately since for any possible value of the parameter the algorithm is related to the same transition probability. However, the joint transition probability will depend on the parameter and the properties of the system, including the damage spreading, will also depend on the parameter.

A remarkable property of the dynamics introduced here is that at infinite temperature it is exactly mapped to the Domany-Kinzel (DK) probabilistic cellular automaton [9]. This gives support to a conjecture by Grassberger [22] according to which the generic class of damage spreading transitions is the same as the directed percolation, to which the transition ocurring in the DK probabilistic cellular automaton belongs.

II Single system

Let us consider a one dimensional lattice where at each site one attaches an Ising variable s_i that takes the values +1 or -1 and denote by s = {s_i} the set of all variables of the lattice. The time evolution of the probability P(s) of state s at discrete time is given by

where W(s^¢|s) is the transition probability from state s to state s^¢ which, for a probabilistic cellular automaton is given by [8]

where w_PCA(|s) is the probability that site i will be in state in the next step given that the present state of the system is s. The DW probabilistic cellular automaton [10] for the one dimensional Ising model is defined by

with

and

where

The site i assumes the state -1 with a probability p_i(s) that does not depend on the central site i. If we choose the linear size of the system to be even the dynamics is decomposed into two independent dynamics for each sublattice. It is possible to show [10] that the DW probabilistic cellular automaton has as the stationary probability distribution the Gibbs probability distribuiton associated with the Ising model, namely,

where b = 1/k_BT, so that it defines a stochastic dynamics that can be assigned to the Ising model.

The transition probabilities w_DW(|s_i-1,s_i+1) are shown in Table 1 where we used the parameter p defined by

The actual computer realization of a probabilistic cellular automaton can be made in several ways. Here, we introduce a family of algorithms that are possible realizations of the DW probabilistic cellular automaton. It has a free parameter a that interpolates between the HD and HB algorithms. At each time step all sites of the lattice are updated in a synchronous way by means of the following algorithm, or equation of motion for the spin variables,

if s_i-1 = s_i+1 and

if s_i-1¹ s_i+1 where x_i is a random number uniformly distributed in the interval [0,1].

When a = 0 one recovers the HD algorithm [15]

and when a = 1/2 one recovers the HB algorithm [14,15]

It is straightforward to show that the algorithm defined by Eqs. (9) and (10) yields the one-site transition probability given by Eqs. (4) and (5) for any value of the parameter a.

III Coupled system

Let us denote by s = {s_i} and = {_i} the configurations of the system and its replica, respectively. All sites of the system and its replica are updated in a synchronous way according to the algorithm

if s_i-1 = s_i+1 and

if s_i-1 ¹ s_i+1 and

if

_i-1 = _i+1 and

if

_i-1¹ _i+1. Notice that the random number x_i is the same for both systems.

The coupled system will be described by a four-state probabilistic cellular automaton defined by the time evolution

of the joint probability P(s; ) of state (s; ) at discrete time where W(s^¢; ¢|s; ) is the joint transition probability from state (s; ) to (s^¢; ¢), given by

From the stochastic equation of motion given by Eqs. (14), (15), (16), and (17), we deduce the joint transition probabilities w(;|s_i-i,s_i+1; _i-1, _i+1) that the site i of the system and the replica assume the values and , respectively. The resultant joint transition probabilities are displayed in Table 2 and are valid for 0 < a < p. For p < a < 1, the algorithm yields a joint transition probability which is independent of a and is the one that results by formally replacing, in Table 2, a by p. The joint transition probability fulfill the following properties

which reflects the condition that the system and the replica follow their own dynamics independent of the coupling.

The joint transition probabilities satisfy also the following properties. (a) Reflection symmetry in which the states of sites i-1 and i+1 are interchanged, that is, s_i-1« s_i+1 and _i-1« _i+1. (b) System-replica symmetry in which the states of the system and the replica are interchanged, that is, s_i« _i for all sites. (c) Up-down symmetry defined by the transformation s_i« -s_i and _i« - _i for all sites.

The Hamming distance, which characterizes the damage spreading, is defined by

which is also the order parameter related to the damage spreading phase transition.

IV Relation with the DK automaton

In this section we show an exact relation between the stochastic dynamics defined in Section III and the DK probabilistic cellular automaton [9]. If we let h_i be the occupation variable attached to site i, that is, h_i = 0 or 1 according to whether site i is empty or occupied by one particle, then the transition probabilities w_DK(|h_i-1,h_i+1) of the DK cellular automaton is given by

The DK cellular automaton displays a critical line in the phase diagram p₁ versus p₂ that separates the absorbing state, for which the density of particles is zero, and the active state, for which the density is nonzero.

Now, let us denote by h_i the coupling variable associated to the dynamics of Section III that takes the value 1 or 0 according whether s_i¹ _i or s_i = _i respectively, given by

The relation between the Hamming distance and the coupling variables is just

The joint transition probabilities in the variables h_i and s_i are defined by

where

_i = s_i(1-2h_i)

Summing over the coupling variable we get the following property

which reflects, as in Eq. (20), the condition that the system follows its own dynamics independent of the coupling. The main property we wish to show, however, is that for infinite temperature, that is, for p = 1/2 we have

with the DK transition probabilites defined by p₂ = 0 and p₁ = 1-2a. This means that the subsystem defined by the variables {h_i} follows a dynamics identical to the DK probabilistic cellular automaton. From relation (27) it follows that the Hamming distance coincides with the order parameter of the active state displayed by the DK automaton.

Yet for the case p = 1/2, it is easy to show that the joint transition probability satisfies the property

with the DK transition probabilites defined by p₂ = 0 and p₁ = 1-2a. Therefore, the s-subsystem and the h-subsystem are statistically independent.

V Mean-field solution

The dynamic mean-field approximation has already been used to study systems in nonequilibirum stationary states [6, 7, 12, 27]. Here we set up equations for an approximate solution of the equation that governs the time evolution of the coupled system. We start by writing down the equations that give the time evolution of the one-site and two-site probabilities, namely,

and

From now on we will drop the subscript and use the prime for quantities calculated at time + 1. To obtain a set of closed equations we use the approximation

which defines the dynamic mean-field pair approximation.

The probabilities P(s₁; ₁) and P(s₁, s₃; ₁, ₃) cannot be considered all independent variables. Taking into account that they should have the reflection symmetry and the system-replica symmetry and, in addition, assuming the up-down symmetry the probabilities are related as follows

These seven variables are not yet independent. Only three of them can be considered independent which we choose to be Y, B, and D. The others are related to them by the relations

where P(+) and P(++) are the one-site and two-site probabilities corresponding to a single system. The exact equilibrium solution of the one-dimensional Ising model gives P(+) = 1/2 and P(++) = [1+(tanh bJ)²]/4.

From the time evolution given by Eqs. (32) and (33) and using Eqs. (43), (44), and (45) we obtain the following closed equations for Y, D, and B

where

and

A stationary solution of the evolution equation is such that the stationary probability P(s; ) is zero when s ¹ which corresponds to no damage spreading (Y = 0). From Eqs. (46), (47), and (48) we may obtain solutions with damage spreading (Y ¹ 0). The transition line is obtained by a linear stability analysis of the solution around Y = 0 and by assuming that the variables B and D vanishes linearly with Y. Taking the limit Y ® 0 we find a transition line given by the implicit equation

whose solution is shown in the phase diagram of Fig. 1. In particular, when a = 0 (corresponding to the HD algorithm) we have g = 2(1 – a) which substituted in the equation for the transition line gives

whose solution is a = 0.696173 from which we get p = 0.196173 so that J/k_BT_c = 0.352597 and T_c = 2.83610. When a = 1/2 (correspoding to the HB algorithm) there is no transition.

At infinite temperature, p = 1/2, the mean-field transition line gives a = 1/6. Now, using the relation p₁ = 1 – 2a obtained from the equivalence with the DK automaton, and taking into account the result p₁ = 2/3 obtained in [12] in the pair approximation for the DK automaton we have a = 1/6 in coincidence with our present result.

VI Numerical simulations

Our numerical simulations resulted in the transition line shown in Fig. 1. When a = 0 we have obtained p = 0.285(1) which gives J/k_BT_c = 0.230(1) and T_c = 4.35(2) in agreement with the result by Hinrichsen and Domany [15], namely J/k_BT_c = 0.2305. At infinite temperature, p = 1/2, the numerical results give a transition at a = 0.0955(1). Now, using the relation p₁ = 1 – 2a obtained from the mapping of our model to the DK cellular automaton, we obtain p_1c = 0.809(1) in agreement with previous Monte Carlo numerical results, namely p_1c = 0.8095(5) [28].

The determination of the critical line was obtained by using the time dependent method [26, 7, 22]. We started with two one-dimensional lattices (system and replica) with L = 1000 sites. Both lattices were initialized with completly indepedendent random configurations so that half the spins were damaged at the begining (Y = 1/2). The update was done in a synchronized way by using the algorithm defined by Eqs. (14), (15), (16), and (17), with the same random number for both lattices. The density of damages Y(t), obtained by taking the averages over 2000 samples, were collected from t = 1 to t = 1500 Monte Carlo steps. At the critical point we expect the following asymptotic time behavior [22]

Therefore, a double-log plot of Y versus t will be linear at the critical point. In Fig. 2 we show how the critical value of p was found when a = 0.075. Several values of p, the ones shown in Fig. 2, were checked in order to find a linear behavior in a log-log plot of Y versus t. Our estimate in this case gives p_c = 0.455(1) for a = 0.075. The straight line fitted to the numerical data gives d = 0.16(1) in agreement with a transition belonging to the direct percolation universality class [7]. For other values of a the procedure was the same.

VII Conclusion

We have introduced a family of algorithms to describe the time evolution of the one-dimensional Ising model The family of algorithms interpolates between the HB and the HD algorithms and the resulting stochastic dynamics corresponds to the DW probabilistic cellular automaton. Coupling a system with its replica by using the same sequence of random numbers, we have determined the joint transition probability which defines a four-state probabilistic cellular automaton. By using a dynamic pair mean-field approximation and Monte Carlo simulations we have found that the stochastic dynamics defined by the family of algorithms displays a line of critical points separating a phase where the damage spreads and a phase where it does not. One important feature of the joint stochastic dynamics studied here is that at infinite temperature the joint dynamics is exactly mapped into the DK probabilistic cellular automaton. This result together with the Monte Carlo simulations give support to a conjecture by Grassberger according to which the damage spreading transition is in the universality class of the directed percolation.

Received on 9 April, 2003

• [1] R. J. Glauber, J. Math. Phys. 4, 294 (1963).
• [2] Monte Carlo Methods in Statistical Physics, edited by K. Binder, 2nd. ed. (Springer-Verlag, Berlin, 1986).
• [3] K. Kawasaki in Phase Transitions and Critical Phenomena, edited by C. Domb and M. S. Green (Academic Press, New York, 1972), vol. 2, p. 443.
• [4] N. G. van Kampen, Stochastic Process in Physics and Chemistry (North-Holland, Amsterdam, 1981).
• [5] T. M. Liggett, Interacting Particle Systems (Spinger-Verlag, New York, 1985).
• [6] T. Tomé and M. J. de Oliveira, Phys. Rev. A 40, 6643 (1989).
• [7] J. Marro and R. Dickman, Nonequilibrium Phase Transition in Lattice Models (Cambridge University Press, Cambridge, 1999).
• [8] T. Tomé e M. J. de Oliveira, Dinâmica Estocástica e Irreversibilidade (Editora da Universidade de Săo Paulo, Săo Paulo, 2001).
• [9] E. Domany and W. Kinzel, Phys. Rev. Lett. 53, 447 (1984).
• [10] B. Derrida and G. Weisbuch, Europhys. Lett. 4, 657 (1987).
• [11] J. L. Lebowitz, C. Maes and E. R. Speer, J. Stat. Phys. 59, 117 (1990).
• [12] T. Tomé, Physica A 212, 99 (1994).
• [13] E. P. Gueuvoghlanian and T. Tomé, Int. J. Mod. Phys. B 11, 1245 (1997).
• [14] M. N. Barber and B. Derrida, J. Stat. Phys. 51, 877 (1988).
• [15] H. Hinrichsen and E. Domany, Phys. Rev. E 56, 94 (1997).
• [16] T. Tomé and M. J. de Oliveira, Braz. J. Phys. 27, 525 (1997).
• [17] M. J. de Oliveira, Int. J. Mod. Phys. B 10, 1313 (1996).
• [18] P. Grassberger, Physica A 214, 547 (1995).
• [19] M. Creutz, Ann. Phys. 167, 62 (1986).
• [20] H. Stanley, D. Stauffer, J. Kertesz, and H. Herrmann, Phys. Rev. Lett. 59, 2326 (1987).
• [21] A. M. Mariz, H. J. Hermann, and L. de Arcangelis, J. Stat. Phys. 59, 1043 (1990).
• [22] P. Grassberger, J. Stat. Phys. 79, 13 (1995).
• [23] H. Hinrichsen, J. S. Weitz and E. Domany, J. Stat. Phys. 88, 617 (1997).
• [24] E. Arashiro and J. R. Drugowich de Felício, Braz. J. Phys. 30, 677 (2000).
• [25] M. L. Martins, H. F. Verona de Rezende, C. Tsallis and A. C. N. de Magalhăes, Phys. Rev. Lett. 66, 2045 (1991).
• [26] P. Grassberger and A. de la Torre, Ann. Phys. (NY) 122, 373 (1979).
• [27] R. Dickman, Phys. Rev. A 34, 4246 (1986).
• [28] H. Rieger, A. Schadschneider and M. Schreckenberg, J. Phys. A 27, L423 (1994).

Publication Dates

• Publication in this collection
  11 Nov 2003
• Date of issue
  Sept 2003

History

• Received
  09 Apr 2003
• Accepted
  09 Apr 2003