A general creation-annihilation model with absorbing states

A general creation-annihilation model with absorbing states

Wellington G. Dantas; Armando Ticona^*; Jürgen F. Stilck

Instituto de Física, Universidade Federal Fluminense, Campus da Praia Vermelha, Niterói, RJ, 24.210-340, Brazil

ABSTRACT

A one dimensional non-equilibrium stochastic model is proposed where each site of the lattice is occupied by a particle, which may be of type A or B. The time evolution of the model occurs through three processes: autocatalytic generation of A and B particles and spontaneous conversion A ® B. The two-parameter phase diagram of the model is obtained in one- and two-site mean field approximations, as well as through numerical simulations and exact solution of finite systems extrapolated to the thermodynamic limit. A continuous line of transitions between an active and an absorbing phase is found. This critical line starts at a point where the model is equivalent to the contact process and ends at a point which corresponds to the voter model, where two absorbing states coexist. Thus, the critical line ends at a point where the transition is discontinuous. Estimates of critical exponents are obtained through the simulations and finite-size-scaling extrapolations, and the crossover between universality classes as the voter model transition is approached is studied.

I. INTRODUCTION

The phase transitions exhibited by stochastic models with absorbing states have attracted much attention in recent years, particularly in order to identify and understand the aspects which determine the universality classes in those models. Most of these models have not been solved exactly, but a variety of approximations allow quite conclusive results regarding their critical properties. Stochastic models are, of course, well fitted for simulations, but closed form approximations and other analytical approaches have also been useful in investigating their behavior [1].

One of the simplest and most studied model of this type is the contact process (CP), which was conceived as a simple model for the spreading of an epidemic and proven to display a continuous transition between the absorbing and an active state, even in one dimension [2]. Actually, it was found that the CP is equivalent to other models such as Schlögl's lattice model for autocatalytic chemical reactions [3] and Reggeon Field Theory (RFT) [4]. The CP belongs to the direct percolation (DP) universality class, together with others models such as the Ziff-Gulari-Barshad model of catalysis [5] and branching and annihilating walks with an odd offspring [6]. The DP conjecture states that all phase transitions between an active and an absorbing state in models with a scalar order parameter, short range interactions and no conservation laws belong to this class [7]. This conjecture was verified in all cases studied so far [8].

Here we study a generalization of the CP, with an additional parameter, so that the CP transition point becomes a critical line. This model is similar to the model proposed by Hinrichsen [9] in the particular case where his parameter q is set equal to one. Since the symmetry properties of this generalized model are the same of the CP, it is expected that this critical line should belong to the DP universality class. However, at one point of this line the model is equivalent to the zero temperature Glauber model [10], also called the voter model [11], which displays a spin inversion (or particle-hole) symmetry and therefore belongs to another universality class (the compact directed percolation (CDP) class). Thus the critical line in the phase diagram of the generalized model starts at the CP model and ends at the voter model, a crossover between the two universality classes being observed.

In section II we define the model and show its equivalence to the CP and the voter model in the appropriate limits. The phase diagram of the model is obtained in one- and two-site approximations in section III. In section IV, results of simulations are shown which lead to numerical estimates of the critical line and of dynamic critical exponents. Other estimates of the critical line and static critical exponents are obtained through the exact diagonalization of the time evolution operator of the model for finite lattices, extrapolated to the thermodynamic limit using finite size scaling in section V. Final comments and the conclusion may be found in section VI.

II. DEFINITION OF THE MODEL

Each one of the N sites of a one-dimensional lattice with periodic boundary conditions are occupied by particles A or B. No holes are allowed. Thus, the state of the system at a given time t is described by the occupation variables h = (h₁,h₂,...,h_N), where h_i = 0 or 1 if site i is occupied by particles B or A, respectively.

The time evolution of the system is defined by the following Markovian rules:

1. A site i of the lattice is chosen randomly.

2. If the site is occupied by a particle B, it becomes occupied with a particle A with a transition rate equal to p_an_a/2, where n_a is the number of A particles in first neighbor sites of i.

3. If site i is occupied by a particle A, it may become a particle B through two processes

• Spontaneously, with a transition rate

  p_c

• Through an autocatalytic reaction, with a rate p_bn_b/2, where n_b is the number of B particles in the first neighbors of site i.

We define the time step in such a way that the parameters p_a, p_b, and p_c, which are non-negative, obey the normalization p_a + p_b + p_c = 1, so that only two of them are independent. For convenience, we will discuss the behavior of the model in the (p_a, p_c) plane without loss of generality.

The probability P(h, t) to find the system in state h at time t obeys the master equation

where hⁱ corresponds to the following configuration

and w_i(h) is the transition rate of the model, given by

where l = p_a/(1 – p_a), g = (1 – p_c)/p_a, and the sum is over first neighbors of site i.

The equation for the time evolution of the mean number of A particles at site i, áh_iñ, may be obtained from equations (1) and (3), being given by

and a homogeneous solution, that is, r º áh_iñ and f º áh_ih_i+dñ "i, is given by equation

In principle, we are not able to solve equation (5) because the function f(t) is not known. We may write a differential equation for f(t), but in this equation three variable terms such as áh_ih_jh_kñ will appear, so that a infinite hierarchy of equations will be obtained. A systematic approximate solution is given below. Although it does not provide precise results, it still furnishes a qualitative picture of the behavior of the system.

This model, similarly to what happens in other out of equilibrium systems, displays absorbing states, which are such that once they are reached, they will never be left. The evolution rules of the model define the state where all sites are occupied by B particles (r = 0) as absorbing. Besides this stationary state, others may exist such that r = lim_{t® ¥} r(t) is nonzero. Such states are called active.

It may be useful to remark that this model may be mapped to a spin system if we describe sites occupied by A and B particles by an Ising spin variables s_i = –1 and s_i = –1, respectively. In these variables, the transition rate will be given by

where a = (p_a + p_b + 2p_c)/2, b = (p_a – p_b – 2p_c)/(pa + p_b + 2p_c), e = (p_a + p_b)/(p_a + p_b + 2p_c), and x = (p_a – p_b)/(p_a + p_b+ 2p_c).

Finally, this model corresponds to two known models in particular limits. If we make p_b = 0 or g = 1 the well known contact process is recovered [2],

If now we take p_a = p_b and p_c = 0 in the spin formulation of the model, the zero temperature linear Glauber model is recovered, also known as the voter model,

III. MEAN FIELD SOLUTION

An approximate solution of equation (5) is obtained if we write f as a function of r. This approximation does not account for correlations and is known as the one site mean field approximation, assuming áh_ih_jñ = áh_iñáh_jñ, that is f = r². Within this approximation, equation (5) may be written as

and the stationary solution in terms of the parameters p_a and p_c is

Therefore, in the one-site approximation the absorbing state is separated from the active state by a continuous transition line located at p_a = 1/2 and the behavior of r as a function of p_a for some values of p_c is shown in Fig. (4). As expected, as p_c is increased, the density of A sites decreases, since annihilation of A particles is favored. It may be shown that close to the transition line the order parameter behaves as r ~ D^b, D = (p_a – 1/2) > 0, with b = 1 and, at the critical line, the density decays as r ~ t^–j, with j = 1. These exponents are identical to the mean field results for the contact process.

A somewhat better approximation is obtained if the set of equations for r and f is simultaneously solved. As remarked above, an equation for f = áh_ih_jñ shows mean values of products of three h variables. Therefore, in order to get a closed set of two equations, these mean values should be approximated as functions of f and r. In doing so, we are neglecting higher order correlations and this is known as the two-site mean field approximation.

The following notation is helpful to obtain the second differential equation of the two-site approximation. We call r = P(•), where P(•) is the probability that a site is occupied by an A particle and the probability to find a B particle is P() = 1–P(•). We also define the three probabilities of the configurations of two neighboring sites, and using the relations

together with the evolution rules stated above, the following equation for P(•) is found

As already mentioned the probabilities of three site clusters appear in the equation. These will be written in terms of two-site probabilities through the so called pair approximation

Within this approximation, a closed system of two equations is found

where u º P(•), r º P() = 1 – r - u and s º P(••) = r – u. The first equation above corresponds to equation (5), whereas the second is obtained from equation (13) using the pair approximation (14).

The stationary solutions for the densities are

where . In this approximation the critical value of p_a is a function of p_c. The critical exponents, as expected, have still their mean-field values, only non-universal parameters are different in the two-site approximation as compared to the one-site calculation above. Figure (5) shows r as a function of p_a for some values of p_c in the two-site approximation.

It may be remarked that in the limit p_c ® 0 both approximations show a discontinuous transition between two absorbing states, which happens at p_a = 1/2. In this limit the model corresponds to a spin model where a given spin is reversed with a probability proportional to the number of spins in first neighbor sites which are in the opposite direction. This might be called a biased voter model when p_a is different from p_b.

IV. SIMULATIONS

Although, as seen above, the model exhibits an active state in part of the parameter space in the approximate solutions, numerical simulations are useful providing some confidence that this state is not an artifact of the approximation. To our knowledge, for non-equilibrium models an argument similar to the Peierls construction for equilibrium models is still not available, and therefore no systematic proof of the existence of the active state may be given, although this was accomplished using probabilistic arguments by Harris in his original work on the contact process [2].

The simulation is done for discrete time, and may be described by the following steps:

1. Initially, the system of N sites with periodic boundary conditions is in a state with just one A particle.

2. A list of all sites occupied by A particles is stored, and at each time step one of them is chosen randomly.

3. Once the site is chosen, a random number p uniformly distributed in the interval [0, 1] is generated, if p < p_a a B particle is replaced by an A particle in one of the first neighbors, if possible. Otherwise, the A particle at the chosen site will be turned into a B particle either through the spontaneous or through the autocatalytic process.

4. To define the process switching A to B, another random number q is generated. If q < p_c/(p_b + p_c) the change is spontaneous, otherwise it will happen with a probability proportional to the number of B particles in first neighbors of the chosen site.

5. The time interval associated with the steps above is Dt = 1/N_A, where N_A is the number of the sites occupied by A particles before the step. The process is repeated until either a maximum time t_max is attained or the absorbing state N_A = 0 is reached.

6. Several runs are done and mean values are calculated as a function of time.

The mean number of A particles, áN_A ñ, as a function of time for simulations with N = 10000 sites, t_max = 100000 and 10000 repetitions is displayed in Fig. (6). The number of sites is sufficiently high to ensure that in the simulations the cluster of A particles is much smaller than the system, thus avoiding finite-size effects in the time interval considered. In a region of the (p_a, p_c) plane an active state is found for long times in the simulations. Three curves are shown in Fig. (6). The convex one corresponds to an active stationary state, while the concave curve signals that the system will reach the absorbing state. These two curves are separated by the third, where a power law áN_Añ ~ t^q is seen, corresponding to the critical condition. We adopt, for all critical exponents, the notation proposed in the review article by Hinrichsen [8].

Repeating simulations in the region (p_a + p_c < 1), the critical line may be estimated. The results of these calculations are shown in Fig. (7).

The estimate for the critical line provided by the simulations display a concavity which agrees with the results of the the two-site approximation and, as expected, is situated above the mean field results, which usually overestimate the active region in the parameter space. The value of l = p_a/(1 – p_a) at the critical line attains the limiting value l_c = 3.2945±0.0116 as p_b ® 0, in good agreement with the estimated critical value for the one-dimensional contact process l_c = 3.29785(2) [12], as may be appreciated in Fig. (8).

The simulational result for áN_Añ as a function of time in the limit of the voter model (p_a = 1/2, p_c = 0) is displayed in Fig. (9). One notices that the initial number of particles is almost conserved, showing a narrow dispersion. This aspect may also be appreciated in the histograms shown in Fig. (10), for (p_a = 1/2, p_c = 0) and (p_a = 0.49, p_c = 0.001), with initially one A particle. For the voter model a small dispersion is shown around the maximum value áN_Añ = 1, while in the other case the maximum is shifted to the absorbing state áN_Añ = 0

The simulations allow us to follow the dynamic evolution of mean values that describe the model, such as the mean number of A particles áN_A(t)ñ, the probability of survival P_s(t) that at the time t at least one particle A is present in the system and the mean square radius R²(t) = áå_i i²h_i(t)ñ / áN_A(t)ñ. These variables satisfy the following scaling forms at the critical point [4]:

Through scaling and hyperscaling relations, these critical exponents define all other exponents of the model, and thus their values identify the universality class.

From the simulations, the dynamic critical exponents defined above were estimated as functions of the rate p_c, as is shown in Fig. (11). For most of the values for p_c, particularly those close to the CP (p_c » 0.2327), the estimated exponents are close to the values known for the CP (q = 0.313686, d = 0.159464, and 2/z = 1.265226 [13]). The relative errors are displayed in Fig. (12). Also, the hyperscaling relation 2d/z = 4d + 2q is satisfied numerically with some imprecision, as may be seen in Fig. (13). Some systematic error seems to be present in the estimates for the mean square radius, which propagates to the hyperscaling relation verification. As p_c ® 0, however, when the voter model is approached, a departure of all estimates from the CP values is apparent.

For models with a discontinuous transition with a spin inversion symmetry, such as the voter model, the dynamic critical exponents are given by q_s = 0, d_s = 1/2 e 2/z_s = 1 and the hyperscaling relation is given by d_s + q_s = d/z_s, corresponding to the compact directed percolation (CDP) universality class [8]. An analysis of the exponents in the limit p_c ® 0 may be performed fitting a cubic curve to the estimates, as is shown in Fig. (14) for the exponent d. The values obtained through this procedure are d = 0.44, q = 0.07, and 2/z = 1.01, which are not far from the known values given above. The variation of the estimates for critical exponents with the parameters of the model close to p_c = 0 is probably an apparent effect, due to large fluctuations observed in the simulations in this region, and the whole critical line, with exception of the terminal voter model point, is in the same universality class of the CP, as suggested by the DP conjecture. This situation is quite similar to what is observed in the behavior of the Domany-Kinzel probabilistic cellular automaton (DKPCA) [14], which in a line of its phase diagram corresponds to the CP with discrete time and parallel update, and also exhibits a crossover from the DP to the compact direct percolation (CDP) universality class [15]. It is believed, although to our knowledge not proven, that changing the update process in a stochastic model may affect non-universal quantities only. The point p_a = 1/2,p_c = 0 may be recognized as a multicritical point, and in its neighborhood any stationary density variable should exhibit the scaling form

The critical exponent associated with the density variable g, e_g, should correspond to the CDP universality class, and the scaling function F(z) is singular at a value z₀ of its argument, which corresponds to the critical line. Thus, the critical line is asymptotically given by p_c = z₀(p_a – 1/2)^f and the amplitude z₀ and the crossover exponent f may be estimated through a fit to the simulational results. Due to the large fluctuations mentioned above, these estimates are not very precise. The values we obtained are z₀ = 0.47 ±0.08 and f = 1.80 ±0.03. In the two site approximation, we found f = 2 both for the model described here and for the DKPCA. Although extensive simulational results for the DKPCA are reported in the literature [16], apparently no estimate for the crossover exponent was obtained from them.

V. NUMERICAL DIAGONALIZATION

Another approach to study 1+1 dimensional stochastic systems is the exact solution of models with increasing numbers of sites N followed by extrapolations to the thermodynamic limit [8]. This is accomplished writing the master equation (1) as

where is the time evolution operator of the model. In the representation in which it is diagonal, vanishing eigenvalues µ₀ = 0 correspond to stationary states of the system. For finite systems with absorbing states, only these states are stationary, and no active stationary state is found.

To study the transition between an active and a stationary state, we may consider the behavior of the eigenvalue with the second smallest absolute value G º µ₁ of the operator . This eigenvalue is related to the quasi-stationary [17] state and will eventually become degenerate with µ₀ in the thermodynamic limit, originating the phase transition. Figure (15) shows the behavior of the gap G for some values of the size N. For a fixed value of p_c. the finite size scaling behavior of the gap G is given by [8]

where is the critical value of p_a for a fixed value of p_c. Defining the quantity

we may estimate the critical point (N) finding the intersection of the curves Y_N and Y_N+1 [18]. This procedure resembles the phenomenological renormalization group. The sequence (N) of estimates for a given value of p_c, with N = 4,5,...M, was extrapolated to the thermodynamic limit N ® ¥ using the BST algorithm [19], which even for the limited number of estimates considered lead to rather precise results. Figure (16) shows the extrapolated estimates for the critical line, compared with the results provided by the simulations and the two-site approximation results. The agreement between the simulation and the present results is apparent. The results for the critical line close to the CP point p_b = 0 were extrapolated to the CP limit using a cubic function and the resulting value was l_c = 3.3081±0.0173, which is in reasonable agreement with a more precise estimate in the literature l_c = 3.29785(2) [12]. We also estimate the crossover exponent f, through a fit to the critical line in the neighborhood of the voter model limit and it is results in f = 2.24 ±0.07.

Once the extrapolated estimates for the critical line were obtained, the critical exponents z and x = z – 1/n_^ may be estimated using the asymptotic scaling forms

for a fixed value of p_c and p_a = , on the critical line. An estimate of the exponents is obtained for each size of the system, and the estimates are extrapolated to the thermodynamic limit. Our results for these two static exponents are shown in Fig. (17).

Again a departure of the estimates from the CP values may be observed as p_c becomes smaller. The CDP values for these exponents are z_s = 2 and x_s = 1. Unlike to what was observed in the simulations no systematic trend to these values was observed in the estimates from exact diagonalization as p_c ® 0. We noticed, even considering values of N up to 15, that as p_c becomes smaller the convergence properties are poorer, and larger sizes may be necessary to obtain better estimates in this region. The estimates which were obtained from the simulational data seem to be more reliable that the ones which follow from the numerical diagonalization results.

VI. CONCLUSION

We discussed the transition between an absorbing (r = 0) and an active (r > 0) steady states in a stochastic model of interacting particles. This model displays a line of continuous transition between these states starting at a point which corresponds to the contact process and ends at a point where the model is equivalent to the voter model, showing a discontinuous transition between two absorbing states. The phase diagram thus obtained is similar to what is found in a problem of equilibrium polymerization in a grand-canonical ensemble on an anisotropic square lattice [20], where the polymer is modeled as a self-avoiding walk and the activity of a horizontal link is equal to x and a vertical link is associated with an activity y. A continuous transition is observed in general in the (x, y) plane, between a non-polymerized (r = 0) and a polymerized (r > 0), where r corresponds to the fraction of lattice sites visited by the walk. However, at x = 0 or y = 0 the walk is one-dimensional, and a discontinuous transition is found between the non-polymerized and the fully polymerized (r = 1) phases [21].

The localization of the critical line in the parameter space of the model was found in the one- and two-site approximations, as well as through simulations and numerical diagonalization of the time evolution operator. The evidences obtained from estimates of critical exponents suggest that the whole critical line belongs to the DP universality class, and a crossover to the CDP universality class is observed at the terminal voter model point. From the estimates for the location of the critical line close to the voter model point the crossover exponent f = 1.80 ±0.03 was obtained. Although the estimated value for f is smaller than the mean-field result f = 2, the latter is inside the confidence interval of the estimate, so that it is not possible to decide if the crossover exponent has a non-classical value. Due to the difficulties we had with the simulations and the numerical diagonalization in the multicritical region, we suspect that another approach is necessary to clear this point. We are presently addressing this point through a series expansion approach. Another extension of this work would be the study of the crossover between the universality class in the Domany-Kinzel cellular automaton, which corresponds to the CP with parallel update in a subspace of its parameter space and also has a point of its phase diagram which corresponds to the voter model.

Acknowledgments

We thank Prof. Ronald Dickman for many helpful discussions and a critical reading of the manuscript. This research was partially supported by the Brazilian agencies CAPES, CNPq and FAPERJ, particularly through the project PRONEX-CNPq-FAPERJ/171.168-2003.

REFERENCES

[1] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999)

[2] T.E. Harris, Ann. Probab. 2, 969 (1974).

[3] F. Schlögl, Z. Phys. 253, 147 (1972).

[4] P. Grassberger and A. De La Torre, Ann.Phys. 122, 373 (1979).

[5] R. M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 56, 2553 (1986).

[6] H. Takayasu, A. Tretyakov, and A. Yu, Phys. Rev. Lett. 68, 3060 (1992).

[7] H. K. Janssen, Z. Phys. B, 42, 151 (1981); P. Grassberger, Z. Phys. B, 47, 365 (1982).

[8] H. Hinrichsen, Adv. Phys.49, 815 (2000).

[9] H. Hinrichsen, Phys. Rev. E 55.219 (1997).

[10] R. J. Glauber, J. Math. Phys. 4, 294 (1963).

[11] T. M. Ligget, Interacting particle systems, Springer, Berlin (1985).

[12] R. Dickman and I. Jensen, Phys. Rev. Lett. 67, 2391 (1991); I. Jensen and R. Dickman, J. Stat. Phys. 71, 89 (1993).

[13] I. Jensen, J. Phys A 32, 5233 (1999).

[14] E. Domany and W. Kinzel, Phys. Rev. Lett. 53, 311 (1984); W. Kinzel, Z. Phys. B 58, 229 (1985).

[15] R. Dickman and A. Y. Tretyakov, Phys. Rev. E 52, 3218 (1995).

[16] G. F. Zebende and T. J. P. Penna, J. Stat. Phys. 74, 1273 (1994).

[17] R. Dickman and R. Vidigal, J. Phys. A 35, 1145 (2002).

[18] E. Carlon, M. Henkel and U. Schollwöck, Eur. Phys. J. B 12, 99 (1999).

[19] T.G. Schmaltz et al, J. Phys. A 17, 445 (1984).

[20] J. F. Stilck, Macromol. Symp. 81, 321 (1994).

[21] P. Pfeuty and J. C. Wheeler, Phys. Rev. A. 27, 2178 (1983).

Received on 28 February, 2005

*

Present address: Instituto de Investigaciones Fisicas

Universidad Mayor de San Andres

Casilla 8635, La Paz, Bolivia

Publication Dates

• Publication in this collection
  06 Sept 2005
• Date of issue
  June 2005

History

• Received
  28 Feb 2005