Critical behavior in reaction-diffusion systems exhibiting absorbing phase transitions

Critical behavior in reaction-diffusion systems exhibiting absorbing phase transitions

Géza Ódor

Research Institute for Technical Physics and Materials Science, H-1525 Budapest, P.O.Box 49, Hungary

ABSTRACT

Phase transitions of reaction-diffusion systems with a site occupation restriction, particle creation requiring n > 2 parents, and in which explicit diffusion of single particles (A) is possible, are reviewed. Arguments based on mean-field approximation and simulations are given which support novel kind of nonequilibrium criticality. These are in contradiction with the implications of a suggested phenomenological, multiplicative noise Langevin equation approach and with some recent numerical analyses. Simulation results for one- and two-dimensional binary spreading model, 2A ® 4A, 4A ® 2A, reveal a new type of mean-field criticality characterized by the critical exponents a = 1/3 and b = 1/2, as suggested in a recent preprint [cond-mat/0210615].

I Introduction

The classification of universality classes of second order phase transitions is still one of the most important uncompleted tasks of statistical physics. Recently phase transitions of genuine nonequilibrium systems have been investigated intensively in reaction-diffusion (RD) type models exhibiting absorbing states [1-3]. There has been a hope that in such homogeneous systems symmetries and spatial dimensions are the most significant factors like in equilibrium ones, but gradually it turned out that there may be other relevant factors as well. The best known example is the parity conserving class (PC), which differs from the robust universality class of directed percolation (DP). The DP hypothesis stated by Janssen and Grassberger [4, 5], according to which in one component systems exhibiting continuous phase transitions to a single absorbing state (without extra symmetry and inhomogeneity or disorder) short ranged interactions can generate DP class transitions only. However parity conservation itself proved to be an insufficient condition in many cases [6-9] and rather an underlying A ® 3A, 2A ® Æ (BARW2) process [10] of particles and the Z₂ symmetry of absorbing states is necessary for this class [11]. On the other hand parity conservation in N-component branching and annihilating random walk (N-BARW) systems [10], or by triplet production models [12] was found to be responsible for novel classes again. In one dimensional, multi-component reaction-diffusion systems site restriction turned out to be a relevant, new factor [13, 14]. Global conservation laws by directed percolation and lattice gas models were shown to be irrelevant [15-18], while systems with multiple absorbing states [19] or with multi-components also exhibit DP class scaling behavior [20].

Another important puzzle has been investigated intensively during the past three years that emerges at phase transitions of binary production reaction-diffusion systems [8,9,20-29] (PCPD). In these systems particle production by pairs competes with pair annihilation and single particle diffusion. If production wins steady states with finite particle density appear in (site restricted), while in unrestricted (bosonic) models the density diverges. By lowering the production/annihilation rate a doublet of absorbing states without symmetries emerges. One of such states is completely empty, the other possesses a single wandering particle. In case of site restricted systems the transition to absorbing state is continuous. It is important to note that these models do not break the DP hypothesis, because they exhibit multiple absorbing states which are not frozen, isolated particles may diffuse in them. However no corresponding symmetry or conservation law has been found yet. A non-DP phase transition in a binary production system was already mentioned in the early work of Grassberger [30]. A corresponding bosonic field theoretical model the annihilation fission (AF) process was introduced and studied by Howard and Täuber [21]. These authors claim a non-DP transition in AF, because the action does not contain a linear mass term and the theory is not renormalizible perturbatively, unlike the Reggeon field theory of DP. In field theories of models exhibiting DP class transitions the canonical A ® Æ, A ® 2A reactions are generated by the renormalization transformation unlike here. Further facts opposing DP criticality are the set of different mean-field exponents and the different upper critical dimension of binary production PCPD like models (d_c = 2 vs. d_c = 4)[21, 9].

Subsequent numerical studies reported somewhat different critical exponents, but there has been a consensus for about two years that this model should possess novel, non-DP transition. The first density matrix study by Carlon et al. [22] did not support a DP transition, but reported exponents near to those of the PC class. Since the PCPD does not conserve the particle number modulo 2, neither exhibits Z₂ symmetric absorbing states, PC criticality was unfavored. Simulation studies by Hinrichsen [23] and Ódor [24] and coherent anomaly calculations by Ódor [24, 26] resulted in novel kind of critical behavior, although there was an uncertainty in the precise values of critical exponents. Exponent estimates showed diffusion (D) dependence that was supported by pair mean-field results [24], possessing two distinct classes as a function of D. Recently Park et al. reported a well defined set of critical exponents in different versions of binary production PCPD-like processes [31]. However these simulations were done at a fixed, high diffusion/annihilation rate and agree with Ódor's corresponding results [32]. Kockelkoren and Chaté on the other hand claim another set of critical exponents [33] that agrees with Ódor's low diffusion/annihilation data.

The PCPD model can be mapped onto a two-component model [25] in which pairs are identified as a particle species following DP process and single particles as another, coupled species following annihilating random walk. Simulations of such a two-component system at D = 0.5 showed a continuous phase transition with exponents agreeing with those of the PCPD for high diffusions. This model is similar to another one [34], which exhibits global particle number conservation as well. Field theory [34] and simulations [35, 36] for the latter model reported two different universality classes as the function of D. It would be interesting to see if this conservation law is relevant or not like in case of the DP [18].

Interestingly, higher level cluster mean-field approximations result in a single class behavior by varying D and it turned out that by assuming logarithmic corrections the single class scenario can be supported by simulations too [32]. The origin of such logarithmic corrections may be a marginal perturbation between pairs and single particles in a coupled system description. A field theoretical explanation would be necessary.

Two more recent studies [37, 38] reported non-universality in the dynamical behavior of the PCPD. While in the former one Dickman and Menezes explored different sectors (a reactive and a diffusive one) in the time evolution and gave non-DP exponent estimates, in the latter one Hinrichsen set afloat a speculative conjecture that the ultimate long time behavior might be characterized by DP scaling behavior. In a subsequent preprint [39] Hinrichsen provided a discussion about the possibility of the DP transition based on a series of assumptions. His starting point is a Langevin equation that is mapped onto a wetting process by Cole-Hopf transformation. By analyzing this process within a certain potential he gave arguments for a DP transition. While this Langevin equation with real noise is valid for the bosonic version of PCPD at and above the critical point, its usage in case of site occupancy restricted models is hypothetical, as the noise can be complex at the transition point and may even change sign at the transformation. Furthermore the diffusive field of solitary particles is neglected.

In a very recent preprint [40] Barkema and Carlon continue this line and show that some simulation and density matrix renormalization results may also be interpreted as a signal of a phase transition belonging to the DP class. By assuming correction to scaling exponents that are equal to DP exponents and relevant up to quadratic or 3-rd order in the asymptotic limit they fitted their numerical results in case of two independent exponents. The extrapolated results are close to DP values for D = 0.5. However for smaller D-s and by surface critical exponents this technique gave exponent estimates which are out of the error margin of DP.

Another novel class that may appear in triplet production systems was proposed in [33, 41] (TCPD). This reaction-diffusion model differs from the PCPD in that for new particle generation at least three particles have to meet. For such generalizations Park et al. proposed a phenomenological Langevin equation that exhibits real, multiplicative noise [41]. By simple power counting they found that the triplet model exhibits distinct mean-field exponents and upper critical dimension 4/3 < d_c < 8/3. The simulations in 1d [41] indeed showed non-trivial critical exponents, which do not seem to correspond to any known universality classes. Kockelkoren and Chaté reported similar results in stochastic cellular automata (SCA) versions of general nA ® (n + k)A, mA ® (m – l)A type models [33], where multiple particle creation on a given site is suppressed by an exponentially decreasing creation probability (p^N/2) of the particle number. They claim that their simulation results in 1d are in agreement with the fully occupation number restriction counterparts and set up a general table of universality classes, where as a function of n and m only 4 non-mean-field classes exist, namely the DP class, the PC class, the PCPD and TCPD classes. However more extensive simulations of 1 and 2 dimensional site restricted lattice models [12] do not support some of these results in case of different triplet and quadruplet models. In the 3A ® 4A 3A ® Æ triplet model 1d numerical data can be interpreted as mean-field behavior with logarithmic corrections and in two dimensions clear mean-field exponents appear, hence the upper critical dimension is d_c = 1, which contradicts the Langevin equation prediction. Surprisingly other non-trivial critical behavior were also detected in the 3A ® 5A 3A ® A parity conserving triplet model and in some quadruplet models [12]. The cause of differences between the results of these studies is subject of further investigations. Again proper field theoretical treatment would be important.

The classification of universality classes of nonequilibrium systems by the exponent m of a multiplicative noise in the Langevin equation was suggested some time ago by Grinstein et al. [42]. However it turned out that there may not be particle systems corresponding to real multiplicative noise cases [21] and an imaginary part appears as well if one derives the Langevin equation of a RD system starting from the master equation in a proper way. Furthermore for higher-order processes the emerging nonlinearities in the master equation action do not allow a rewriting in terms of Langevin-type stochastic equations of motion, hence for high-order processes like those of the TCPD a Langevin representation may not exist.

This situation resembles to some extent the decade-long debate over the critical phase transition of driven diffusive systems [43-45]. The latest papers on this topic suggest that the phenomenological Langevin equation originally set up for such systems does not correspond exactly to the lattice models investigated. Simulations for different lattice models show that instead of an external current, anisotropy is the real cause of the critical behavior observed in simulations [46, 47].

II Mean-field classes

In this section I show that mean-field classes of site restricted lattice models with general microscopic processes of the form

with n > 1, m > 1, k > 0, l > 0 and m –l > 0 are different from those of the DP and PC processes backing numerical results which claim novel type of criticality below d_c. The mean-field equation that can be set up for the lattice version of these processes (with creation probability s and annihilation or coagulation probability = 1 – s) is

where r denotes the site occupancy probability and a is a dimension dependent coordination number. Each empty site has a probability (1-r) in mean-field approximation, hence the need for k empty sites at a creation brings in a (1 – r)^k probability factor. By expanding (1 – r)^k and keeping the lowest order contribution one can see that for site restricted lattice systems a rⁿ⁺¹-th order term appears with negative coefficient that regulates eq.(2) in the active phase. In the inactive phase one expects a dynamical behavior dominated by the mA ® Æ process, for which the particle density decay law is known r(t) µ t^1/(m–1) [10]. The steady state solutions were determined in [12] analytically and one can distinguish three different situations at the phase transition: (a) n = m, (b) n > m and (c) n < m .

A. The n = m symmetric case

As discussed in [12] the leading order singularities of steady state solution can be obtained. By approaching the the critical point s_c = l/k+l in the active phase the steady state density vanishes continuously as

with the order parameter exponent exponent b^MF = 1. At the critical point the density decays with a power-law

with a_MF = b^MF/ = 1/n, hence = n, providing a series of mean-field universality classes for n > 1 (besides DP an PC where = 1) and backing the results, which claim novel non-trivial transitions below the critical dimension. Unfortunately determining the the value of d_c is a non-trivial task without a proper Langevin equation. These scaling exponents can be obtained from bosonic, coarse grained formulation too [41], where a rⁿ⁺¹-th order term, with negative coefficient has to be added by hand to suppress multiple site occupancy. It is known however that hard-core particle blocking may be a relevant perturbation in d = 1 dimension [14], so for cases where the upper critical dimension is d_c > 1 the site restricted, N > 1 cluster mean-field approximation that takes into account diffusion would be a more adequate description of the model (see [48]).

B. Then > m case

In this case the mean-field solution provides first order transition (see [12]), hence it does not imply anything with respect to possible classes for models below the critical dimension (d < d_c). Note however, that by higher order cluster mean-field approximations, where diffusion plays a role the transition may turn continuous (see for example [49-51]). The simulations by Kockelkoren and Chaté report DP class transition for such models in one dimension [33].

C. The n < m case

In this case the critical point is at zero branching rate s_c = 0, where the density decays with a^MF = 1/(m – 1) as in case of the n = 1 branching and m = l annihilating models studied by Cardy and Täuber [10] (BkARW classes). However the steady state solution for particle production with n > 1 parents gives different b exponents than those of BkARW classes, namely b^MF = 1/(m – n), defining a new series of mean-field classes [12], for a simplicity I shall call them PkARW classes. It is important to note that one has not found a symmetry or conservation law corresponding to these classes. These mean-field classes imply novel kinds of critical behavior for d < d_c. For n = 2, m = 3 the mean-field exponents are b = 1 and a = 1/2 agreeing with the mean-field exponents of the PCPD class. Indeed for n = 2, m = 3 [33] reports PCPD class dynamical criticality. This supports the expectation that non-mean-field classes follow the distinctions observed in the corresponding mean-field classes. To go further in sections and I investigate the phase transition of the simplest unexplored PkARW classes, in the n = 2, m = 4 model.

D. The role of k and l

In the mean-field approximation k and l do not affect the universal properties, however simulations in one dimension showed [12] that in case of the m = n = k = l = 3 model the critical point was shifted to zero branching rate and a BkARW class transition emerged there contrary to what was expected for n = m. For a stochastic cellular automaton version of these reactions [33] reported a non-trivial critical transition. In general one may expect such effects for large k and l values, for which N-cluster mean-field approximation - that takes into account diffusion - would give a better description.

III Simulations of the 2A ® 4A, 4A ® 2A model in two dimensions

In Sec. II, I introduced PkARW mean-field classes for n < m. Here I explore the phase transition in the simplest model from this class, in the 2A ® 4A, 4A ® 2A model with diffusion rate D = 0.5. Two dimensional simulations were performed on L = 1000 linear sized lattices with periodic boundary conditions. One Monte Carlo step (MCS) - corresponding to dt = 1/P (where P is the number of particles) - is built up from the following processes. A particle and a number x Î (0,1) are selected randomly; if x < D = 0.5 a site exchange is attempted with one of the randomly selected empty nearest neighbors (nn); if x > D = 0.5 two particles are created with probability s at randomly selected empty nn sites provided the number of nn particles is > 2; or if x > 0.5 two particles are removed with probability 1 – s. The simulations were started from fully occupied lattices and the particle density (r(t)) decay was followed up to 4 × 10⁵ MCS.

First the critical point was located by measuring the dynamic behavior of r(t). It turned out that the transition is at zero branching rate (s_c = 0). The density decay was analyzed by the local slopes defined as

where I used m = 4. As Fig.1 shows, the local slope curve for t > 10⁵ MCS extrapolates to the mean-field value a = 0.334(1). This value agrees with the mean-field value a^MF = 1/3.

Density decays for several s-s in the active phase (0.0002 < s < 0.05) were followed on logarithmic time scales and averaging was done over ~ 100 independent runs in a time window which exceeds the level-off time by a decade. The steady state density in the active phase near the critical phase transition point is expected to scale as

Using the local slopes method one can get a precise estimate for b as well as for the corrections to scaling

As one can see in Fig.2 the effective exponent clearly tends to the expected mean-field value b = 0.5 as s ® 0. Assuming a correction to scaling of the form

nonlinear fitting results in d = 0.44(1) correction to scaling exponent.

Besides these corrections to scaling assumptions I also tried to apply different, lowest order logarithmic corrections to the data, but these fits gave exponents slightly away from mean-field values and the corresponding coefficients proved to be very small, therefore I concluded that d_c < 2. In the next section I perform the same analysis in d = 1.

IV Simulations of the 2A ® 4A, 4A ® 2A model in one dimension

The simulations in one dimension were carried out on L = 20000 sized systems with periodic boundary conditions. The initial states were again fully occupied lattices, and the density of particles is followed up to 4×10⁶ MCS. An elementary MCS consists of the following processes:

(a)AÆ « ÆA with probability D,

(b)AAAA ® Æ AAÆ with probability (1 – s)(1 – D),

(c)AAÆÆ ® AAAA or ÆÆAA ® AAAA with probability s(1 – D),

The critical point was again located at s_c = 0. As one can see in Fig. 3 there is a crossover of the local slopes for t > 5×10⁵ MCS and a linear extrapolation for this region results in a = 0.329(5), agreeing with the mean-field value.

The steady state data were analyzed in the active region for 0.0003 < s < 0.5 as in two dimensions, by the local slopes method (eq. 7) and assuming a correction to scaling of the form (8). This resulted in a correction to scaling exponent of d¢ = 0.332. The local slope plotted as a function of s^–d, shown on Fig. 4, extrapolates to b = 0.49(1) in agreement with the mean-field exponent.

V Conclusions

In this paper I reviewed and discussed the state of the art of phase transitions of reaction-diffusion systems exhibiting explicit diffusion and production by n > 1 parents. Arguments are given against the DP criticality recently suggested in some papers. These are supported by a series of mean-field classes that can be classified by the existence of a n = m symmetry in the system and by the n and m values. The need for a proper field theoretical treatment is emphasized. The upper critical dimension in these models is not known. I showed simulation results indicating that for the binary production model, 2A ® 4A, 4A ® 2A, the upper critical dimension is d_c < 1.

Note that while the density decay results in one dimension are in agreement with those of Kockelkoren and Chaté [33], the off-critical order parameter exponent is b^MF = 1/(m – n) which shows that more classes exist at zero branching rate besides the BkARW universality classes. It is still an open question if there is any variant of the PkARW models that exhibits non-mean-field criticality in physical dimensions (d > 1).

Acknowledgements:

Support from Hungarian research funds OTKA (Grant No. T-25286) is acknowledged. The author thanks the trial access to the NIIF Cluster-GRID Hungary.

References

[1] For references see : J. Marro and R. Dickman, Nonequilibrium phase transitions in lattice models, (Cambridge University Press, Cambridge, 1999).

[2] For a review see : H. Hinrichsen, Adv. Phys. 49, 815 (2000).

[3] For a recent review see : G. Ódor, cond-mat/0205644.

[4] H. K. Janssen, Z. Phys. B 42, 151 (1981).

[5] P. Grassberger, Z. Phys. B 47, 365 (1982).

[6] H.Park and H.Park, Physica A 221, 97 (1995).

[7] N. Menyhárd and G. Ódor, J. Phys. A 29, 7739 (1996).

[8] K. Park, H. Hinrichsen, and In-mook Kim, Phys. Rev. E 63, 065103(R) (2001).

[9] G. Ódor, M. A. Santos, M. C. Marques, Phys. Rev. E 65, 056113 (2002).

[10] J. L. Cardy and U. C. Täuber, Phys. Rev. Lett. 77, 4780 (1996); J. L. Cardy and U. C. Täuber, J. Stat. Phys. 90, 1 (1998).

[11] N. Menyhárd, G. Ódor, Braz. J. of Physics, 30, 113 (2000).

[12] G. Ódor, Phys. Rev. E 67, 056114 (2003).

[13] G. Ódor, Phys. Rev. E 63 056108 (2001).

[14] G. Ódor and N. Menyhárd, Physica D 168 305 (2002).

[15] T. Tomé and M. J. de Oliveira, Phys. Rev. Lett. 86, 5643 (2001).

[16] M.M.S. Sabag and M. J. de Oliveira, Phys. Rev. E 66, 036115 (2002).

[17] H. Hilhorst and F. Van Vijland, Phys. Rev.E 65, 035103 (2002).

[18] M. J. de Oliveira, Phys. Rev. E 67, 027104 (2003).

[19] I. Jensen and R. Dickman, Phys. Rev. E 48, 1710 (1993); I. Jensen, Phys. Rev. Lett. 70, 1465 (1993).

[20] G. Ódor, Phys. Rev. E 65, 026121 (2002).

[21] M. J. Howard and U. C. Täuber, J. Phys. A 30, 7721 (1997).

[22] E. Carlon, M. Henkel and U. Schollwöck, Phys. Rev. E 63, 036101-1 (2001).

[23] H. Hinrichsen, Phys. Rev. E 63, 036102-1 (2001).

[24] G. Ódor, Phys. Rev. E 62, R3027 (2000).

[25] H. Hinrichsen, Physica A 291, 275-286 (2001).

[26] G. Ódor, Phys. Rev. E 63, 067104 (2001).

[27] J. D. Noh and H. Park, cond-mat/0109516.

[28] M. Henkel and U. Schollwöck, J. Phys. A 34, 3333 (2001).

[29] M. Henkel and H. Hinrichsen, J. Phys. A 34, 1561 (2001).

[30] P. Grassberger, Z. Phys. B 47, 365 (1982).

[31] K. Park and I. Kim, Phys. Rev E 66, 027106 (2002).

[32] G. Ódor, Phys. Rev E 67, 016111 (2003).

[33] J. Kockelkoren and H. Chaté, Phys. Rev. Lett. 90, 125701 (2003).

[34] F. van Wijland, K. Oerding and H. J. Hilhorst, Physica A 251, 179 (1998).

[35] J. E. de Freitas, L. S. Lucena, L. R. da Silva and H. J. Hilhorst, Phys. Rev. E 61, 6330 (2000).

[36] U. L. Fulco, D. N. Messias and M. L. Lyra, Phys. Rev. E 63, 066118 (2001).

[37] R. Dickman and M. A. F. de Menenzes, Phys. Rev. E 66, 045101 (2002).

[38] H. Hinrichsen, Physica A 320, 249 (2003).

[39] H. Hinrichsen, cond-mat/0302381.

[40] G. T. Barkema and E. Carlon, cond-mat/0302151.

[41] K. Park, H. Hinrichsen and I. Kim, Phys. Rev. E 66, 025101 (2002).

[42] Grinstein G., Muñoz M.A., and Tu Y., Phys. Rev. Lett. 76, 4376 (1996); Tu Y., Grinstein G., and Muñoz M.A. Phys. Rev. Lett. 78, 274 (1997).

[43] J. L. Vallés and J. Marro, J. Stat. Phys. 49, 89 (1989).

[44] K.-t. Leung, Phys. Rev. Lett. 66, 453 (1991).

[45] R.K.P. Zia, L.B.Shaw and B. Schmittman, Physica A 279, 60 (2000).

[46] A. Achahbar, Pedro L. Garrido, J. Marro, Miguel A. Munoz, Phys. Rev. Lett. 87, 195702 (2001).

[47] E. V. Albano and G. Saracco, Phys. Rev. Lett. 88, 145701 (2002).

[48] G. Ódor, Phys. Rev. E 65, 026121 (2002).

[49] G. Ódor, N. Boccara, G. Szabó, Phys. Rev. E 48, 3168 (1993).

[50] G. Ódor and A. Szolnoki, Phys. Rev. E 53, 2231 (1996).

[51] N. Menyhárd and G. Ódor, J. Phys. A 28, 4505 (1995).

Received on 31 March, 2003

Publication Dates

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

History

• Accepted
  31 Mar 2003
• Received
  31 Mar 2003