Contact process with sublattice symmetry breaking
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Contact process with sublattice symmetry breaking
Marcelo Martins de Oliveira ∗ and Ronald Dickman † Departamento de F´ısica e Matem´atica, CAP,Universidade Federal de S˜ao Jo˜ao del Rei,36420-000 Ouro Branco, Minas Gerais - Brazil Departamento de F´ısica and National Instituteof Science and Technology for Complex Systems,ICEx, Universidade Federal de Minas Gerais,C. P. 702, 30123-970 Belo Horizonte, Minas Gerais - Brazil (Dated: November 10, 2018)
Abstract
We study a contact process with creation at first- and second-neighbor sites and inhibition atfirst neighbors, in the form of an annihilation rate that increases with the number of occupied firstneighbors. Mean-field theory predicts three phases: inactive (absorbing), active symmetric, andactive asymmetric, the latter exhibiting distinct sublattice densities on a bipartite lattice. Thesephases are separated by continuous transitions; the phase diagram is reentrant. Monte Carlo simu-lations in two dimensions verify these predictions qualitatively, except for a first-neighbor creationrate of zero. (In the latter case one of the phase transitions is discontinuous .) Our numericalresults confirm that the symmetric-asymmetric transition belongs to the Ising universality class,and that the active-absorbing transition belongs to the directed percolation class, as expected fromsymmetry considerations.
PACS numbers: 05.50.+q,05.70.Ln,05.70.Jk,02.50.EyKeywords: ∗ email: [email protected] † email: dickman@fisica.ufmg.br . INTRODUCTION An absorbing-state phase transition is a far from equilibrium transition between an activeand an inactive phase. When a control parameter such as a creation or annihilation rate isvaried, 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, suchas population dynamics, heterogeneous catalysis, interface growth, and epidemic spreading[1–5]. Interest on this class of phase transitions has been stimulated by recent experimentalrealizations [6, 7].As equilibrium statistical mechanics, absorbing-state phase transitions can be classifiedinto universality classes. However, a complete classification of their critical behavior is stilllacking; much numerical and theoretical work is focused on identifying the factors thatdetermine the universality classes of these transitions.The absorbing-state universality class associated with directed percolation (DP) hasproven to be particularly robust. DP-like behavior appears to be generic for absorbing-state transitions in models with short-range interactions and lacking a conserved density orsymmetry beyond translational invariance [8]. By contrast, models possessing two absorb-ing states linked by particle-hole symmetry belong to the voter model universality class [9].There are also models which, while free of absorbing states, cannot achieve thermal equilib-rium because their transition rates violate the detailed balance principle. An example is themajority vote model. In its ordered phase a Z symmetry is spontaneously broken, leadingto Ising-like behavior in two or more dimensions.In this work we examine whether a spatially structured population model that suffers aphase transition to a single absorbing state can also exhibit a broken-symmetry phase. Ourcandidate for such a model is a modified contact process (CP) with suppression of activity atthe nearest neighbors of active sites. We observe unequal sublattice populations for suitablechoices of the control parameters.The balance of this paper is organized as follows. In the next section we define the modeland analyze its mean-field theory. In Sec. III we present our simulation results; Sec. IV isdevoted to discussion and conclusions. 2 I. MODEL AND MEAN-FIELD THEORY
The CP [10] is a stochastic interacting particle system defined on a lattice, with each site i either occupied by a particle [ σ i ( t ) = 1], or vacant [ σ i ( t ) = 0]. Transitions from σ i = 1 to σ i = 0 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 transitionfrom σ i = 0 to σ i = 1 occurs at rate λr , where r is the fraction of nearest neighbors of site i that are occupied. In the absence of spontaneous creation of particles, the state σ i = 0for all i is absorbing. It can be shown [10] that for a certain critical value λ c the systemundergoes a phase transition between the active and the absorbing state. Although no exactresults are available, the CP has been studied intensively via series expansion and MonteCarlo simulation, and its critical properties are known to high precision [1–5].In order to generate distinct sublattice occupations, we modify the basic contact processin two respects:i) In addition to creation at first neighbors, at rate λ , we allow creation at secondneighbors , at rate λ . For bi-partite lattices such as the square or simple cubic lattice, λ isthe rate of creation in the opposite sublattice, while λ is the rate in the same sublattice asthe replicating particle. Unequal sublattice occupancies are thus favored for λ > λ .ii) We include a nearest-neighbor inhibition effect in the annihilation process. In additionto the intrinsic annihilation rate of unity, there is a contribution of µn , where n is thenumber of occupied first neighbors. Thus if the occupation fraction ρ A of sublattice A ismuch larger than that of sublattice B, any particles created in sublattice B will die outquickly, stabilizing the unequal sublattice occupancies.The one-site mean-field theory (MFT) for this model on a lattice of coordination number q is given by the coupled equations, dρ A dt = − (1 + µq ρ B ) ρ A + ( λ ρ B + λ ρ A )(1 − ρ A ) (1)and dρ B dt = − (1 + µq ρ A ) ρ B + ( λ ρ A + λ ρ B )(1 − ρ B ) (2)which are seen to be symmetric under ρ A ⇌ ρ B . Defining ρ = ρ A + ρ B and φ = ρ A − ρ B ,these equations may be recast as, 3 ρdt = (Λ − ρ − Λ2 ρ − ∆2 φ − µq ( ρ − φ ) ρ (3)and dφdt = (cid:20) ∆ − − λ ρ − µq ( ρ − φ ) (cid:21) φ (4)where Λ ≡ λ + λ and ∆ ≡ λ − λ .From Eq. (3) it is evident that at this level of approximation, the extinction-survivaltransition occurs at Λ = 1, as one would expect. This transition leads to the symmetricactive phase, characterized by ρ > φ = 0. This phase is linearly stable for small ∆. As λ is increased, with λ constant, the symmetric active phase can become unstable, leadingto a phase with | φ | >
0. On further increasing λ , however, ρ can increase to the pointthat the symmetric active phase is again stable. (Perfect sublattice ordering corresponds to ρ = φ , which is only possible for ρ ≤
1; total densities larger than unity tend to suppresssublattice ordering.) Thus the asymmetric active phase should exist for some intermediaterange of λ values.In the symmetric active phase, the stationary activity density is, ρ = 12 κ hp (Λ / + 4 κ (Λ − − (Λ / i (5)with κ ≡ µq /
4. From Eq. (4), this solution is stable with respect to one with φ = 0, when a φ ≡ ∆ − − λ ρ + κρ < . (6)The transition to the asymmetric phase occurs when a φ = 0, which implies(2 λ + Λ) (Λ −
1) = 2( λ −
1) [8 κ ( λ −
1) + (2 λ + Λ)Λ] (7)Eq. (4) may be written as dφ/dt = a φ φ − κφ , showing the expected symmetry under φ →− φ , which in turn suggests that the transition between symmetric- and asymmetric-activephases is Ising-like. The two transitions (active/absorbing and symmetric/asymmetric) co-incide at the point λ = 0, λ = 1. Figure 1 shows the phase diagram for µ = 2, obtainedusing Eq. (7). Of particular note are the observations that (1) the asymmetric-active phaseis observed only for intermediate values of λ , i.e., the phase diagram is reentrant, and (2)4 IG. 1: (Color online) Phase diagram in the λ − λ plane for µ = 2, showing absorbing (ABS), active-symmetric (AS) and active asymmetric (AA) phases. Lines: MFT; symbols: simulation. the asymmetric-active phase does not exist for λ above a certain value, λ ∗ ( µ ). For µ = 2,for example, λ ∗ = 7 . λ ∗ byabout an order magnitude.The line λ = 0 is special, since the subspaces with ρ A > ρ B = 0 (and vice-versa),are also absorbing. MFT yields the following phases: absorbing for λ ≤
1; asymmetric-active-I for 1 < λ < λ I ( µ ); asymmetric-active-II for λ I ( µ ) < λ < λ II ( µ ); and symmetric-active for λ > λ II . The difference between phases asymmetric-active-I and II is that in theformer, ρ = φ (i.e., only one of the sublattices is populated) while in the latter ρ > φ > both sublatticespopulated.) For µ = 2, for example, we find λ I = 30 .
97 and λ II = 41 . linear function of the number of occupied first neighbors, i.e., the mean-field equations dρ A dt = − (1 + µqρ B ) ρ A + ( λ ρ B + λ ρ A )(1 − ρ A ) , (8)and similarly for ρ A ⇌ ρ B . The phase diagram is qualitatively similar to that found above,but the region occupied by the asymmetric-active phase is much reduced. For µ = 2, forexample, sublattice ordering is only found for λ ≤ . III. SIMULATIONS
We performed extensive Monte Carlo simulations of the model on square lattices of linearsize L = 20 , , ...,
320 sites, with periodic boundaries. The simulation algorithm is asfollows. First, a site is selected at random. If the site is occupied, it creates a particle atone of its first-neighbors with a probability p = λ / (1 + λ + λ + µn ), or at one of itssecond-neighbors with a probability p = λ / (1 + λ + λ + µn ). With a complementaryprobability 1 − ( p + p ) the site is vacated. (In order to improve efficiency the sites arechosen from a list which contains the currently N occ occupied sites; we increment the timeby ∆ t = 1 /N occ after each event). For simulations in the subcritical and critical absorbingregime, we employ the quasi-stationary (QS) simulation method detailed in [11], to furtherimprove efficiency.In a series of studies using µ = 2 .
0, we verify that for low values of λ and λ thesystem 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 λ = 0 .
2, the absorbing transition occurs at λ = 1 . λ + λ is somewhat greater than the critical value for the basic process, λ c = 1 . ρ ∼ L − β/ν ⊥ . In Fig. III, we show thatour data follow a power law, with β/ν ⊥ = 0 . β/ν ⊥ = 0 . τ ∼ L ν || /ν ⊥ at criticality, where we find ν || /ν ⊥ = 1 . ν || /ν ⊥ = 1 . m = h ρ i / h ρ i (see Fig. 3) is6nalogous to Binder’s fourth cumulant [14] at an equilibrium critical point, in the sense that m assumes a universal value m c at the critical point. We find = 1 . m = 1 . ln L -4-20246810 l n ρ , l n τ FIG. 2: (Color online) Scaled critical QS density of active sites ln ρ (bottom) and scaled lifetime of the QSstate ln τ (top), versus ln L . Parameters: µ = 2 . λ = 0 . λ = 1 . λ m FIG. 3: (Color online) Quasistationary moment ratio m versus λ , for µ = 2 . λ = 0 .
2. System sizes: L = 20 , , ,
160 (in order of increasing steepness).
In the active phase, we find that as we increase λ at fixed λ , there is a transition toa state with unequal particle densities in the two sublattices ( φ = 0), that is, to the activeasymmetric (AA) phase. A typical configuration in the AA phase is shown in Fig. 4. Ourobservation of a bimodal probability distribution for the order parameter φ (see Fig. III)confirms the existence of the AA phase. Although MFT predicts a symmetry-breaking7
10 15 20 25 30 35 40 45 505101520253035404550
FIG. 4: (Color online) A typical configuration in the active asymmetric phase; particles are represented byblack sites. System size: L = 50. Parameters: µ = 2 . λ = 0 .
2, and λ = 6 . -0.4 -0.2 0 0.2 0.4 φ P ( φ ) FIG. 5:
Probability distribution of the order parameter φ . System size: L = 80. Parameters: µ = 2 . λ = 0 . λ = 2 . . . transition in any number of dimensions, we do not expect such a phase transition in onedimension; simulations on a ring confirm this.In order to classify the AA-AS transition, we perform a finite-size scaling study. First we8nvestigate the behavior of the reduced Binder cumulant [14], given by U = 1 − h φ ih φ i . (9)The intersection points of the cumulants for successive pairs of sizes depend rather weaklyon the sizes, providing a good estimate for the critical point. The value of the cumulant atthe critical point approaches a universal value as L → ∞ . In Fig. 6 we show the results forthe case µ = 2 . λ = 0 .
2. The curves for different sizes intersect at λ = 3 . λ = 31 . L → ∞ .Extrapolating to infinite size, we estimate U ,c = 0 . U ,c = 0 . ... [16] for the two-dimensional Ising model with fully periodic bound-ary conditions. For values of λ between the transitions points, the cumulant approaches2/3, as expected in an ordered phase that breaks a Z (i.e., up-down) symmetry.Further evidence for Ising-like behavior is furnished by the order-parameter variance. Atthe critical point, var( φ ) = h φ i − h φ i , is expected to scale so: L var( φ ) ∝ L γ/ν ; this isconfirmed in Fig. 7. A fit to the data yields the exponent ratio γ/ν = 1 . γ/ν = 7 / λ B i nd e r c u mm u l a n t
30 31 32 33 34 3500.10.20.30.40.50.60.7
FIG. 6: (Color online) Binder cumulant versus λ , for µ = 2 . λ = 0 .
2. System sizes: L = 20 , , , For the λ = 0, simulations show no evidence of the asymmetric-active-II phase predictedby MFT. Simulations reveal only two phase transitions along this line. The first, from9 ln L L v a r φ FIG. 7: (Color online) Scaled order parameter variance ln L var φ versus ln L at the AS-AA transition.System sizes: L = 20 , , ..., µ = 2 . λ = 0 . λ = 3 .
94. The slope of the regression lineis 1.76. the absorbing to the asymmetric active phase, occurs, as expected, at the critical valueof the basic CP on the square lattice. We find λ ,c = 1 . m = 1 . β/ν ⊥ = 0 . ν || /ν ⊥ = 1 . discontinuous . Figure 9 illustrates such transitions for µ = 0 . µ = 2. These resultssignal a qualitative failure of mean-field theory along the line λ = 0. λ m ln L -4-202468 l n ρ , l n τ FIG. 8: (Color online) Quasistationary moment ratio m versus λ , for µ = 2 . λ = 0. System sizes: L = 20 , , ,
160 (in order of increasing steepness). Inset: Scaled critical QS density of active sites ln ρ (bottom) and scaled lifetime of the QS state ln τ (top), versus ln L . Parameters: µ = 2 . λ = 0 and λ = 1 .
10 100 λ ρ , φ FIG. 9: (Color online) Order parameters ρ (solid) and φ (dashed) versus λ , for λ = 0. Parameters: µ = 0 . µ = 2 . IV. CONCLUSIONS
A contact process with creation at both nearest and next-nearest neighbors, and suppres-sion (in the form of an increased annihilation rate) at neighbors of active sites is shown toexhibit a broken-symmetry phase with sublattice ordering. The existence of a continuoussymmetry-breaking transition is predicted by mean-field theory and confirmed in simula-tions on the square lattice. The asymmetric active phase exists for a restricted range of thesecond-neighbor creation rate, λ , when that at first neighbors, λ , is sufficiently small. Fora given value of the suppression parameter µ , there is a certain value, λ ∗ , above which nosublattice order occurs. We note that the possibility of breaking Z symmetry depends onthe fact that creation rates λ and λ apply to different sublattices. We would not expectto observe this symmetry breaking if, for example, there were a creation rate λ for both nearest- and second-neighbors, and a different creation rate at third-neighbor sites.For 0 < λ < λ ∗ , the sequence of phases observed on increasing λ from small values, whileholding µ and λ fixed, is: absorbing, symmetric-active, asymmetric-active, and symmetric-active. The first transition belongs to the universality class of directed percolation, whilethe latter two are Ising-like. For the special case of λ = 0, the first transition is fromabsorbing to asymmetric-active and the second is discontinuous. For both values of µ studiedwe observe continuous transitions for λ >
0. While increasing λ makes the transitionsmoother, we can not exclude the possibility of a discontinuous transition when λ > µ . A complete characterization of the reentrant transition for µ > . λ = 0 marks the meeting of two lines of criticalpoints, one of DP transitions, the other of Ising transitions. The meeting of two such lineshas been studied in the context of a nonequilibrium Potts model by Droz at al. [17], andof the generalized voter model (GVM) by Al Hammal et al. [18]. The former study alsofinds an absorbing state in the DP class, while in the latter, after the Ising and DP linesjoin, the absorbing transition belongs to the GVM class. We note that in these studies themodels have but two (symmetric) absorbing states (i.e., all plus and all minus), whereasfor λ = 0, our model has two (symmetric) absorbing subspaces (activity restricted to onesublattice) and a third, completely inactive absorbing state. In terms of long-time behavior,the absorbing-active phase transition takes the system from the fully inactive state to onewith activity restricted to one sublattice. Staring from all sites active, at the critical point λ = 0, λ = λ c,CP , fluctuations drive local regions into states with activity in only onesublattice, as well as fully inactive regions. We expect that these domains then coarsen, sothat eventually there is activity only in one sublattice; from then on the dynamics is that ofthe basic contact process. Thus it is not surprising that we observe DP-like static behavior.It is nevertheless possible that the short-time critical behavior includes a regime dominatedby domain dynamics as in the GVM. It also seems possible that the discontinuous AA-AStransition observed for λ = 0 belongs to the voter model class. We intend to explore thesequestions in future work.Possible extensions of this work include precise determination of static and dynamic crit-ical properties at the symmetry-breaking transition, studies of the effect of diffusion and/oranisotropy (which may generate modulated phases) and studies of a model in continuousspace. The latter appears to be of particular interest in the context of patterns arising inecosystems or bacterial colonies. We note that the suppression of activity at the nearestneighbors of active sites resembles biological lateral inhibition, known to be important inthe visual system of many animals [19]. Thus extensions of the present study might beuseful in the study of patterned activity on the retina. Finally, if the inhibition rate were adecreasing function (linear or quadratic) of the number of occupied nearest neighbors (i.e., µ < Acknowledgments
This work was supported by CNPq and FAPEMIG, Brazil. We thank Hugues Chat´e forhelpful comments. [1] J. Marro and R. Dickman,
Nonequilibrium Phase Transitions in Lattice Models (CambridgeUniversity Press, Cambridge, 1999).[2] G. ´Odor,
Universality In Nonequilibrium Lattice Systems: Theoretical Foundations (WorldScientific,Singapore, 2007)[3] M. Henkel, H. Hinrichsen and S. Lubeck,
Non-Equilibrium Phase Transitions Volume I: Ab-sorbing Phase Transitions (Springer-Verlag, The Netherlands, 2008).[4] H. Hinrichsen, Adv. Phys. , 815 (2000).[5] G. ´Odor, Rev. Mod. Phys , 663 (2004).[6] K. A. Takeuchi, M. Kuroda, H. Chat´e, and M. Sano, Phys. Rev. Lett. , 234503 (2007).[7] L. Cort´e, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Nature Physics , 420 (2008).[8] H. K. Janssen, Z. Phys. B , 151 (1981);P. Grassberger, Z. Phys. B , 365 (1982).[9] I. Dornic, H. Chat´e, J. Chave and H. Hinrichsen, Phys. Rev. Lett. , 045701 (2001).[10] T. E. Harris, Ann. Probab., , 969 (1974).[11] M. M. de Oliveira and R. Dickman, Phys. Rev. E , 016129 (2005); R. Dickman and M. M.de Oliveira, Physica A , 134 (2005).[12] A. G. Moreira and R. Dickman, Phys. Rev. E , R3090 (1996).[13] R. Dickman, Phys. Rev. E, , R2441 (1999).[14] K. Binder, Phys. Rev. Lett. 47 (1981) 693.[15] R. Dickman and J. Kamphorst Leal da Silva, Phys. Rev. E, , 201 (1993).[17] M. Droz, A. L. Ferreira, and A. Lipowski, Phys. Rev. E , 056108 (2003).[18] O. Al Hammal, H. Chat´e, I. Dornic, and M. A. Mu˜noz, Phys. rev. Lett. , 230601 (2005).
19] W. Lytton,
From Computer to Brain (Springer-Verlag, New York, 2002).[20] D.H. Janzen, Am. Nat. , 501 (1970); J.H. Connell, in
Dynamics of Populations , Ed. P.J.Den Boer and G.R. Gradwell. (Pudoc:Wageningen, 1970)., Ed. P.J.Den Boer and G.R. Gradwell. (Pudoc:Wageningen, 1970).