Contact process with sublattice symmetry breaking
Abstract
We study a contact process with creation at first and secondneighbor sites and inhibition at first neighbors, in the form of an annihilation rate that increases with the number of occupied first neighbors. Meanfield theory predicts three phases: inactive (absorbing), active symmetric, and active asymmetric, the latter exhibiting distinct sublattice densities on a bipartite lattice. These phases are separated by continuous transitions; the phase diagram is reentrant. Monte Carlo simulations in two dimensions verify these predictions qualitatively, except for a firstneighbor creation rate of zero. (In the latter case one of the phase transitions is discontinuous.) Our numerical results confirm that the symmetricasymmetric transition belongs to the Ising universality class, and that the activeabsorbing transition belongs to the directed percolation class, as expected from symmetry considerations.
pacs:
05.50.+q,05.70.Ln,05.70.Jk,02.50.EyI Introduction
An absorbingstate phase transition is a far from equilibrium transition between an active and an inactive phase. When a control parameter such as a creation or annihilation rate is varied, the system undergoes a phase transition from a fluctuating state to an absorbing one. Such transitions are currently of great interest, arising in a wide variety of problems, such as population dynamics, heterogeneous catalysis, interface growth, and epidemic spreading marro ; henkel ; odor07 ; hinrichsen ; odor04 . Interest on this class of phase transitions has been stimulated by recent experimental realizations take07 ; pine .
As equilibrium statistical mechanics, absorbingstate phase transitions can be classified into universality classes. However, a complete classification of their critical behavior is still lacking; much numerical and theoretical work is focused on identifying the factors that determine the universality classes of these transitions.
The absorbingstate universality class associated with directed percolation (DP) has proven to be particularly robust. DPlike behavior appears to be generic for absorbingstate transitions in models with shortrange interactions and lacking a conserved density or symmetry beyond translational invariance grasjans . By contrast, models possessing two absorbing states linked by particlehole symmetry belong to the voter model universality class voter . There are also models which, while free of absorbing states, cannot achieve thermal equilibrium because their transition rates violate the detailed balance principle. An example is the majority vote model. In its ordered phase a symmetry is spontaneously broken, leading to Isinglike behavior in two or more dimensions.
In this work we examine whether a spatially structured population model that suffers a phase transition to a single absorbing state can also exhibit a brokensymmetry phase. Our candidate for such a model is a modified contact process (CP) with suppression of activity at the nearest neighbors of active sites. We observe unequal sublattice populations for suitable choices of the control parameters.
The balance of this paper is organized as follows. In the next section we define the model and analyze its meanfield theory. In Sec. III we present our simulation results; Sec. IV is devoted to discussion and conclusions.
Ii Model and MeanField Theory
The CP harrisCP is a stochastic interacting particle system defined on a lattice, with each site either occupied by a particle [], or vacant []. Transitions from to occur at a rate of unity, independent of the neighboring sites. The reverse transition, at a vacant site, is only possible if at least one of its neighbors is occupied: the transition from to occurs at rate , where is the fraction of nearest neighbors of site that are occupied. In the absence of spontaneous creation of particles, the state for all is absorbing. It can be shown harrisCP that for a certain critical value the system undergoes a phase transition between the active and the absorbing state. Although no exact results are available, the CP has been studied intensively via series expansion and Monte Carlo simulation, and its critical properties are known to high precision marro ; henkel ; odor07 ; hinrichsen ; odor04 .
In order to generate distinct sublattice occupations, we modify the basic contact process in two respects:
i) In addition to creation at first neighbors, at rate , we allow creation at second neighbors, at rate . For bipartite lattices such as the square or simple cubic lattice, is the rate of creation in the opposite sublattice, while is the rate in the same sublattice as the replicating particle. Unequal sublattice occupancies are thus favored for .
ii) We include a nearestneighbor inhibition effect in the annihilation process. In addition to the intrinsic annihilation rate of unity, there is a contribution of , where is the number of occupied first neighbors. Thus if the occupation fraction of sublattice A is much larger than that of sublattice B, any particles created in sublattice B will die out quickly, stabilizing the unequal sublattice occupancies.
The onesite meanfield theory (MFT) for this model on a lattice of coordination number is given by the coupled equations,
(1) 
and
(2) 
which are seen to be symmetric under . Defining and , these equations may be recast as,
(3) 
and
(4) 
where and .
From Eq. (3) it is evident that at this level of approximation, the extinctionsurvival transition occurs at , as one would expect. This transition leads to the symmetric active phase, characterized by and . This phase is linearly stable for small . As is increased, with constant, the symmetric active phase can become unstable, leading to a phase with . On further increasing , however, can increase to the point that the symmetric active phase is again stable. (Perfect sublattice ordering corresponds to , which is only possible for ; total densities larger than unity tend to suppress sublattice ordering.) Thus the asymmetric active phase should exist for some intermediate range of values.
In the symmetric active phase, the stationary activity density is,
(5) 
with . From Eq. (4), this solution is stable with respect to one with , when
(6) 
The transition to the asymmetric phase occurs when , which implies
(7) 
Eq. (4) may be written as , showing the expected symmetry under , which in turn suggests that the transition between symmetric and asymmetricactive phases is Isinglike. The two transitions (active/absorbing and symmetric/asymmetric) coincide at the point , . Figure 1 shows the phase diagram for , obtained using Eq. (7). Of particular note are the observations that (1) the asymmetricactive phase is observed only for intermediate values of , i.e., the phase diagram is reentrant, and (2) the asymmetricactive phase does not exist for above a certain value, . For , for example, . Also shown in Fig. 1 are simulation results (details of our simulation method are discussed in the following section). MFT and simulation are in qualitative agreement, but as is usually the case, MFT overestimates the regions in parameter space corresponding to the active and ordered phases. In particular, MFT overestimates by about an order magnitude.
The line is special, since the subspaces with and (and viceversa), are also absorbing. MFT yields the following phases: absorbing for ; asymmetricactiveI for ; asymmetricactiveII for ; and symmetricactive for . The difference between phases asymmetricactiveI and II is that in the former, (i.e., only one of the sublattices is populated) while in the latter . The transitions between these phases are continuous. (Naturally, the asymmetricactiveII and symmetricactive phases can only be reached from initial conditions with both sublattices populated.) For , for example, we find and .
We also studied a version in which the suppression (enhanced annihilation rate) is a linear function of the number of occupied first neighbors, i.e., the meanfield equations
(8) 
and similarly for . The phase diagram is qualitatively similar to that found above, but the region occupied by the asymmetricactive phase is much reduced. For , for example, sublattice ordering is only found for . Since meanfield theory offers at best a qualitative description of the model, we turn, in the following section, to simulation.
Iii Simulations
We performed extensive Monte Carlo simulations of the model on square lattices of linear size sites, with periodic boundaries. The simulation algorithm is as follows. First, a site is selected at random. If the site is occupied, it creates a particle at one of its firstneighbors with a probability , or at one of its secondneighbors with a probability . With a complementary probability the site is vacated. (In order to improve efficiency the sites are chosen from a list which contains the currently occupied sites; we increment the time by after each event). For simulations in the subcritical and critical absorbing regime, we employ the quasistationary (QS) simulation method detailed in qssim , to further improve efficiency.
In a series of studies using , we verify that for low values of and the system becomes trapped in the absorbing (inactive) state. When we increase and/or the system undergoes a phase transition from the inactive to the active symmetric (AS) phase. For example, for , the absorbing transition occurs at . (Note that at this point is somewhat greater than the critical value for the basic process, , on the square lattice agmrd .) At the critical point, the quasistaionary order parameter is expect to decay as a power law, . In Fig. 2, we show that our data follow a power law, with , in good agreement with the DP value dickman99 . The lifetime of the QS state at criticality, where we find , also shown in Fig. 2. This value again is in very good agreement with the best known DP value . The moment ratio (see Fig. 3) is analogous to Binder’s fourth cumulant bind81 at an equilibrium critical point, in the sense that assumes a universal value at the critical point. We find , in comparison to the known DP value dicjaf . These results permit us to state that, as expected, the absorbing transition belongs to the universality class of directed percolation.
In the active phase, we find that as we increase at fixed , there is a transition to a state with unequal particle densities in the two sublattices (), that is, to the active asymmetric (AA) phase. A typical configuration in the AA phase is shown in Fig. 4. Our observation of a bimodal probability distribution for the order parameter (see Fig. 5) confirms the existence of the AA phase. Although MFT predicts a symmetrybreaking transition in any number of dimensions, we do not expect such a phase transition in one dimension; simulations on a ring confirm this.
In order to classify the AAAS transition, we perform a finitesize scaling study. First we investigate the behavior of the reduced Binder cumulant bind81 , given by
(9) 
The intersection points of the cumulants for successive pairs of sizes depend rather weakly on the sizes, providing a good estimate for the critical point. The value of the cumulant at the critical point approaches a universal value as . In Fig. 6 we show the results for the case and . The curves for different sizes intersect at , and again at when .
Extrapolating to infinite size, we estimate at the transitions, consistent with the kami93 for the twodimensional Ising model with fully periodic boundary conditions. For values of between the transitions points, the cumulant approaches 2/3, as expected in an ordered phase that breaks a (i.e., updown) symmetry.
Further evidence for Isinglike behavior is furnished by the orderparameter variance. At the critical point, var, is expected to scale so: var; this is confirmed in Fig. 7. A fit to the data yields the exponent ratio , again in good agreement with the known value of for the Ising model in two dimensions.
For the , simulations show no evidence of the asymmetricactiveII phase predicted by MFT. Simulations reveal only two phase transitions along this line. The first, from the absorbing to the asymmetric active phase, occurs, as expected, at the critical value of the basic CP on the square lattice. We find , and confirm that this transition belongs to the DP class (see Fig. 8, where we find , and , consistent with the values for DP in 2+1 dimensions, cited above). The second transition, between the activeasymmetric and the active symmetric phases, is discontinuous. Figure 9 illustrates such transitions for and . These results signal a qualitative failure of meanfield theory along the line .
Iv Conclusions
A contact process with creation at both nearest and nextnearest neighbors, and suppression (in the form of an increased annihilation rate) at neighbors of active sites is shown to exhibit a brokensymmetry phase with sublattice ordering. The existence of a continuous symmetrybreaking transition is predicted by meanfield theory and confirmed in simulations on the square lattice. The asymmetric active phase exists for a restricted range of the secondneighbor creation rate, , when that at first neighbors, , is sufficiently small. For a given value of the suppression parameter , there is a certain value, , above which no sublattice order occurs. We note that the possibility of breaking symmetry depends on the fact that creation rates and apply to different sublattices. We would not expect to observe this symmetry breaking if, for example, there were a creation rate for both nearest and secondneighbors, and a different creation rate at thirdneighbor sites.
For , the sequence of phases observed on increasing from small values, while holding and fixed, is: absorbing, symmetricactive, asymmetricactive, and symmetricactive. The first transition belongs to the universality class of directed percolation, while the latter two are Isinglike. For the special case of , the first transition is from absorbing to asymmetricactive and the second is discontinuous. For both values of studied we observe continuous transitions for . While increasing makes the transition smoother, we can not exclude the possibility of a discontinuous transition when , for larger values of . A complete characterization of the reentrant transition for , including the location of the possible tricritical points separating the line of second order transitions from that of first order transitions is left for future work.
The absorbingactive transition at marks the meeting of two lines of critical points, one of DP transitions, the other of Ising transitions. The meeting of two such lines has been studied in the context of a nonequilibrium Potts model by Droz at al. droz , and of the generalized voter model (GVM) by Al Hammal et al. alhammal . The former study also finds an absorbing state in the DP class, while in the latter, after the Ising and DP lines join, the absorbing transition belongs to the GVM class. We note that in these studies the models have but two (symmetric) absorbing states (i.e., all plus and all minus), whereas for , our model has two (symmetric) absorbing subspaces (activity restricted to one sublattice) and a third, completely inactive absorbing state. In terms of longtime behavior, the absorbingactive phase transition takes the system from the fully inactive state to one with activity restricted to one sublattice. Staring from all sites active, at the critical point , , fluctuations drive local regions into states with activity in only one sublattice, as well as fully inactive regions. We expect that these domains then coarsen, so that eventually there is activity only in one sublattice; from then on the dynamics is that of the basic contact process. Thus it is not surprising that we observe DPlike static behavior. It is nevertheless possible that the shorttime critical behavior includes a regime dominated by domain dynamics as in the GVM. It also seems possible that the discontinuous AAAS transition observed for belongs to the voter model class. We intend to explore these questions in future work.
Possible extensions of this work include precise determination of static and dynamic critical properties at the symmetrybreaking transition, studies of the effect of diffusion and/or anisotropy (which may generate modulated phases) and studies of a model in continuous space. The latter appears to be of particular interest in the context of patterns arising in ecosystems or bacterial colonies. We note that the suppression of activity at the nearest neighbors of active sites resembles biological lateral inhibition, known to be important in the visual system of many animals lytton . Thus extensions of the present study might be useful in the study of patterned activity on the retina. Finally, if the inhibition rate were a decreasing function (linear or quadratic) of the number of occupied nearest neighbors (i.e., ), one should observe symbiotic coexistence of activity in the two sublattices. The latter is reminiscent of the JanzenConnel hypothesis for the maintenance of tree species biodiversity in tropical rainforests janzenconnell .
Acknowledgments
This work was supported by CNPq and FAPEMIG, Brazil. We thank Hugues Chaté for helpful comments.
References
 (1) J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
 (2) G. Ódor, Universality In Nonequilibrium Lattice Systems: Theoretical Foundations (World Scientific,Singapore, 2007)
 (3) M. Henkel, H. Hinrichsen and S. Lubeck, NonEquilibrium Phase Transitions Volume I: Absorbing Phase Transitions (SpringerVerlag, The Netherlands, 2008).
 (4) H. Hinrichsen, Adv. Phys. 49, 815 (2000).
 (5) G. Ódor, Rev. Mod. Phys 76, 663 (2004).
 (6) K. A. Takeuchi, M. Kuroda, H. Chaté, and M. Sano, Phys. Rev. Lett. 99, 234503 (2007).
 (7) L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Nature Physics 4, 420 (2008).

(8)
H. K. Janssen, Z. Phys. B 42, 151 (1981);
P. Grassberger, Z. Phys. B 47, 365 (1982).  (9) I. Dornic, H. Chaté, J. Chave and H. Hinrichsen, Phys. Rev. Lett. 87, 045701 (2001).
 (10) T. E. Harris, Ann. Probab., 2, 969 (1974).
 (11) M. M. de Oliveira and R. Dickman, Phys. Rev. E 71, 016129 (2005); R. Dickman and M. M. de Oliveira, Physica A 357, 134 (2005).
 (12) A. G. Moreira and R. Dickman, Phys. Rev. E 54, R3090 (1996).
 (13) R. Dickman, Phys. Rev. E, 60, R2441 (1999).
 (14) K. Binder, Phys. Rev. Lett. 47 (1981) 693.
 (15) R. Dickman and J. Kamphorst Leal da Silva, Phys. Rev. E, 58 4266 (1998).
 (16) G. Kamieniarz and H.W.J. Blöte, J. Phys. A: Math. Gen. 26, 201 (1993).
 (17) M. Droz, A. L. Ferreira, and A. Lipowski, Phys. Rev. E 67, 056108 (2003).
 (18) O. Al Hammal, H. Chaté, I. Dornic, and M. A. Muñoz, Phys. rev. Lett. 94, 230601 (2005).
 (19) W. Lytton, From Computer to Brain (SpringerVerlag, New York, 2002).
 (20) D.H. Janzen, Am. Nat. 104 , 501 (1970); J.H. Connell, in Dynamics of Populations, Ed. P.J. Den Boer and G.R. Gradwell. (Pudoc:Wageningen, 1970).