Quasi-stationary simulation: the subcritical contact process

Abstract

We apply the recently devised quasi-stationary simulation method to study the lifetime and order parameter of the contact process in the subcritical phase. This phase is not accessible to other methods because virtually all realizations of the process fall into the absorbing state before the quasi-stationary regime is attained. With relatively modest simulations, the method yields an estimate of the critical exponent nu|| with a precision of 0.5%.

Contact process; Quasi-stationary simulations; Nonequilibrium phase transitions; Scaling

Quasi-stationary simulation: the subcritical contact process

Marcelo Martins de Oliveira; Ronald Dickman

Departamento de Física, ICEx, Universidade Federal de Minas Gerais, 30123-970, Belo Horizonte, MG, Brazil

ABSTRACT

We apply the recently devised quasi-stationary simulation method to study the lifetime and order parameter of the contact process in the subcritical phase. This phase is not accessible to other methods because virtually all realizations of the process fall into the absorbing state before the quasi-stationary regime is attained. With relatively modest simulations, the method yields an estimate of the critical exponent n_|| with a precision of 0.5%.

Keywords: Contact process; Quasi-stationary simulations; Nonequilibrium phase transitions; Scaling

I. INTRODUCTION

Stochastic processes with an absorbing state arise frequently in statistical physics [1, 2], epidemiology [3], and related fields. Phase transitions to an absorbing state in spatially extended systems, exemplified by the contact process [4, 5], are currently of great interest in connection with self-organized criticality [6], the transition to turbulence [7], and issues of universality in nonequilibrium critical phenomena [8-11].

Systems exhibiting a phase transition to an absorbing state possess (for appropriate values of the control parameter), an active (nonabsorbing) stationary state in the infinite-size limit. But if the system size is finite, the process must eventually end up in the absorbing state. The quasi-stationary (QS) distribution for such a system provides a wealth of information about its behavior. (Since the only true stationary state for a finite system is the absorbing one, simulations of "stationary" properties of models with an absorbing state in fact study the quasi-stationary regime.)

Conventional Monte Carlo methods entail a somewhat complicated procedure for determining QS properties: many independent realizations are performed (using the same initial configuration), and the mean f(t) of some property (for example the order parameter) is evaluated over the surviving realizations at time t. At short times times f(t) exhibits a transient as it relaxes toward the QS regime; at long times it fluctuates wildly as the surviving sample decays. In the supercritical phase (where an infinite system survives forever) one is normally able to identify an intermediate regime free of transients and with limited fluctuations, which can be used to estimate the QS value of f. (Even in this case the method requires careful scrutiny of the data and is not always free of ambiguity [12].) Deep in the subcritical phase, however, conventional simulations are impractical, as nearly all realizations fall into the absorbing state before the quasi-stationary regime is attained. (In other words, one no longer has a separation of the time scales for relaxation to the QS state and for its survival. The effective potential governing the process no longer possesses a local minimum away from the absorbing state.)

We recently devised a simulation method that yields quasi-stationary properties directly [13]. We showed that the method reproduces known scaling properties of the contact process (CP); in [14] it was used to study the static correlation function of the model. In studies of the critical point, the QS simulation method requires an order of magnitude less computer time than conventional simulations, to obtain results of comparable precision. Here we use the method to study the lifetime and the order parameter of the one-dimensional contact process in the subcritical phase.

In the following section we review the method and define the contact process. Then in Sec. III we present our results for the lifetime of the CP on a ring, in the subcritical phase. We summarize our findings in Sec. IV.

II. BACKGROUND

We begin by reviewing the definition of the quasi-stationary distribution. Consider a continuous-time Markov process X_t taking values n = 0, 1, 2,...,S, with the state n = 0 absorbing. (In any realization of the process, if X_t = 0, then X_t' must be zero for all subsequent times, t' > t.) The transition rates w_m,n (from state n to state m) are such that w_m,0 = 0, "m > 0. Let p_n(t) denote the probability of state n at time t, given some initial state X₀. The survival probability P_s(t) = å_{n> 1}p_n(t) is the probability that the process has not become trapped in the absorbing state up to time t. We suppose that as t ® ¥ the p_n, normalized by the survival probability P_s(t), attain a time-independent form. The quasi-stationary distribution _n is then defined via

with

₀º 0. The QS distribution is normalized so:

The basis for our simulation method is the observation [15] that the QS distribution represents a stationary solution to the equation,

where w_n = å_m w_m,n is the total rate of transitions out of state n, and r_n = å_m w_n,m q_m is the flux of probability into this state. To see this, consider the master equation (Eq. (3) without the final term) in the QS regime. Substituting q_n(t) = P_s(t)_n, and noting that in the QS regime dP_s/dt = - ₀ = - P_s å_m w_0,m

_m, we see that the r.h.s. of Eq. (3) is identically zero if q_n = _n for n > 1. The final term in Eq. (3) represents a redistribution of the probability r₀ (transfered to the absorbing state in the original master equation), to the nonabsorbing subspace. Each nonabsorbing state receives a share equal to its probability.

Although Eq. (3) is not a master equation (it is nonlinear in the probabilities q_n), it does suggest a simulation scheme for sampling the QS distribution. In a Monte Carlo simulation one generates a set of realizations of a stochastic process. In what follows we call a simulation of the original process X_t (possessing an absorbing state) a conventional simulation. We define a related process, , whose stationary probability distribution is the quasi-stationary distribution of X_t. (Note that in order to have a nontrivial stationary distribution, cannot possess an absorbing state.) The probability distribution of is governed by Eq. (3), which implies that for n > 0 (i.e., away from the absorbing state), the evolution of is identical to that of X_t. (Since Eq. (3) holds for n > 0, the process must begin in a non-absorbing state.) When X_t enters the absorbing state, however, instead jumps to a nonabsorbing one, and then resumes its "usual" evolution (i.e., with the same transition probabilities as X_t), until such time as another visit to the absorbing state is imminent. The probability that jumps to state j (when a visit to state 0 is imminent), is the QS probability _j.

A subtlety associated with this procedure is that the QS distribution is needed to determine the evolution of when X_t visits the absorbing state. Although one has no prior knowledge of the QS distribution _n, one can, in a simulation, use the history (0 < s < t) up to time t, to estimate the _n. (There is good evidence, after all, that the surviving sample in a conventional simulation converges to the QS state after an initial transient.) In practice this is accomplished by saving (and periodically updating) a sample n₁, n₂, ..., n_M of the states visited. As the evolution progresses, will visit states according to the QS distribution. We therefore update the sample {n₁,n₂,...,n_M} by occasionally replacing one of these configurations with the current one. In this way the distribution for the process (and the sample drawn from it), will converge to the QS distribution (i.e., the stationary solution of Eq. (3)) at long times. Summarizing, the simulation process has the same dynamics as X_t, except that when a transition to the absorbing state is imminent, is placed in a nonabsorbing state, selected at random from a sample over the history of the realization. In effect, the nonlinear term in Eq. (3) is represented as a memory in the simulation.

A. The contact process

To explain how our method works in practice, we detail its application to the contact process (CP) [4, 5, 8]. In the CP, each site i of a lattice is either occupied (s_i (t) = 1), or vacant (s_i (t) = 0). Transitions from s_i = 1 to s_i = 0 occur at a rate of unity, independent of the neighboring sites. The reverse transition is only possible if at least one neighbor is occupied: the transition from s_i = 0 to s_i = 1 occurs at rate lr, where r is the fraction of nearest neighbors of site i that are occupied; thus the state s_i = 0 for all i is absorbing. (l is a control parameter governing the rate of spread of activity.)

Although no exact results are available, the CP has been studied intensively via series expansion and Monte Carlo simulation. The model has attracted much interest as a prototype of a nonequilibrium critical point, a simple representative of the directed percolation (DP) universality class. Since its scaling properties have been discussed extensively [8-10] we review them only briefly. The best estimate for the critical point in one dimension is l_c = 3.297848(20), as determined via series analysis [16]. As the critical point is approached, the correlation length x and correlation time t diverge, following x µ |D|^-n|| and t µ |D|^-n||, where D = (l- l_c)/l_c is the relative distance from the critical point. The order parameter (the fraction of active sites), scales as r µ D^b for D > 0.

Two characteristic times, t_C and t_L, may be identified in the contact process. The first is a relaxation time that governs the decay of temporal correlations in the stationary state: C(t) º ár(t₀) r(t₀+t) ñ- r² ~ e^-t/t_c. The second is a lifetime, determining the asymptotic decay of the survival probability (starting from a spatially homogeneous initial condition) via P(t) ~ e^-t/t_L. The two characteristic times exhibit the same scaling properties in the critical region. In the supercritical or active phase (D > 0), the lifetime grows exponentially with system size and D, while t_C remains finite. Our interest here is in the lifetime (denoted simply as t) of the quasi-stationary state. In the QS probability distribution there is a nonzero flux of probability to the absorbing state,

where is the QS probability of the configuration with exactly one active site and w₀₁ = 1 is the transition rate from this configuration to the absorbing state; the QS lifetime t = 1/r₀. (In QS simulations we take t to be the mean time between successive attempts to visit to the absorbing state.)

An important point in interpreting our simulation results concerns the finite-size scaling (FSS) behavior of t. According to the usual FSS hypothesis [17], finite-size corrections to critical properties are functions of the ratio L/x, or, equivalently, of the quantity DL^{1/n^}. The lifetime is therefore expected to follow (in the subcritical phase, D < 0),

where the scaling function (x) µ xⁿ|| for small x (so that t does not diverge in a finite system), while in the opposite limit ® ₀, a constant. The scaling hypothesis leads to the familiar relation, t(0,L) ~ Lⁿ||^{/n^} at the critical point, and suggests that we attempt to collapse data for diverse system sizes by plotting Dⁿ|| t versus D* º DL^{1/n^}. For the order parameter the expected finite-size scaling form is [8],

In the subcritical regime, the order parameter must decay to zero µ L^-1 as L ® ¥, for any D < 0, so that (x) ~ |x|^{-n^} as x ® -¥. On the other hand, for D = 0 and L finite, r must be nonzero and nonsingular, implying (x) ~ x^-b for x ® 0. Thus r ~ |D|^{b- n^} for D large and negative.

III. SIMULATION RESULTS

In the QS simulations we use a list of size M = 2 ×10³ - 10⁴, depending on the lattice size. The process is simulated in 15 realizations, each of 5 ×10⁸ time steps. As is usual, annihilation events are chosen with probability 1/(1+l) and creation with probability l/(1+l). A site i is chosen from a list of currently occupied sites, and, in the case of annihilation, is vacated. In a creation event, a nearest-neighbor of site i is selected at random and, if it is currently vacant, it becomes occupied. The time increment associated with each event is Dt = 1/N_occ, where N_occ is the number of occupied sites just prior to the attempted transition [8].

In the initial phase of the evolution, the list of saved configurations is augmented whenever the time t increases by one, until a list of M configurations has been accumulated. From then on, we update the list (replacing a randomly selected entry with the current configuration), with a certain probability p_rep, whenever t advances by one unit. A given configuration therefore remains on the list for a mean time of M/p_rep. (Values of p_rep in the range 10^-3 - 10^-4 are used.)

Figure 1 shows the QS lifetime t as a function of D, for lattice sizes L = 20, 40, 80,...,2560. For the larger system sizes, power-law dependence on D is evident, before the lifetime saturates (at very small values of |D|), due, as anticipated, to finite-size effects. In Fig. 2 these data are collapsed using the known values of the DP exponents [18], n_^ = 1.09684(6) and n_|| = 1.73383(3). The data collapse is quite good for the larger system sizes. A least-squares fit to the data in the linear portion of the graph (|D| > 0.02) yields a slope of -1.738(12), in reasonable agreement with the accepted value for n_||. An alternative method for analyzing the data is to estimate, for each value of D, the infinite-size limiting value of the lifetime by plotting t(D,L) versus 1/L (see Fig. 3). Plotting the resulting estimates on log scales, we find a slope of -1.735(9), again in good agreement with the standard value for the exponent n_||.

We verified that the distribution p_T of the lifetime (that is, of the time interval between successive attempts to visit the absorbing state) is exponential:

This is as expected: in the quasi-stationary state the transition rate to the absorbing state is time-independent.

Turning to the order parameter, we see from Eq. (6) that a plot of r* º|D|^-b r versus D* should exhibit a data collapse. The asymptotic behavior of the scaling function implies that r* µ (D*)^{-n^} for large D*. These scaling properties are verified in Fig. 4. (We use the accepted value b = 0.2765. While scaling of the order parameter in the subcritical regime was verified in Ref [8], here we are able to extend the range of D* by an order of magnitude, using the QS simulation technique.)

In Fig. 5 we show the probability distribution p(r) of the order parameter in a large system (L = 2560) for three values of D. As expected, the distribution broadens and shifts toward larger values of r as l approaches the critical value. The probability distribution follows, to a good approximation, the scaling form

where árñ º r(D,L) is the mean value and is a scaling function. This is verified (Fig. 6) by plotting p* = árñp(r) versus r/árñ. We see that the scaling function attains its maximum near r* = 0.6, and that it falls off rapidly as r* ® 0. On the other side of the maximum it exhibits a roughly exponential tail.

IV. SUMMARY

We use the quasi-stationary simulation method to study the lifetime and order parameter of the one-dimensional contact process in the subcritical phase. Our results confirm the expected scaling properties of the lifetime of the quasi-stationary state, and of the order parameter. The QS simulation method is the first to allow such an analysis deep in the subcritical regime. With a rather modest expenditure of computer time, the approach yields an estimate of the critical exponent n_|| that agrees with the accepted value to within uncertainty, with a precision of about 0.5%. Analysis of quasi-stationary properties in the subcritical regime should therefore be a useful tool in the study of absorbing-state phase transitions.

Acknowledgment

This work was supported by CNPq and FAPEMIG, Brazil.

Received on 26 September, 2005

• [1] N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1992).
• [2] C. W. Gardiner, Handbook of Stochastic Methods, (Springer-Verlag, Berlin, 1990).
• [3] M.S. Bartlett, Stochastic Population Models in Ecology and Epidemiology (Methuen, London, 1960).
• [4] T. E. Harris, Ann. Probab. 2, 969 (1974).
• [5] T. Liggett, Interacting Particle Systems (Springer-Verlag, Berlin, 1985).
• [6] R. Dickman, M. A. Muńoz, A. Vespignani, and S. Zapperi, Braz. J. Phys. 30, 27 (2000).
• [7] T. Bohr, M. van Hecke, R. Mikkelsen, and M. Ipsen, Phys. Rev. Lett. 86, 5482 (2001), and references therein.
• [8] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
• [9] H. Hinrichsen, Adv. Phys. 49 815, (2000).
• [10] G. Odor, Rev. Mod. Phys 76, 663 (2004).
• [11] S. Lübeck, Int. J. Mod. Phys. B 18, 3977 (2004).
• [12] S. Lübeck and P. C. Heger, Phys. Rev. E 68, 056102 (2003).
• [13] M. M. de Oliveira and R. Dickman, Phys. Rev. E 71, 016129 (2005).
• [14] R. Dickman and M. M. de Oliveira, Physica A 357, 134 (2005).
• [15] R. Dickman and R. Vidigal, J. Phys. A 35, 1145 (2002).
• [16] I. Jensen and R. Dickman, J. Stat. Phys. 71, 89 (1993).
• [17] M. E. Fisher and M. N. Barber, Phys. Rev. Lett. 28, 1516 (1972).
• [18] I. Jensen, J. Phys. A 32, 5233 (1999).

Publication Dates

• Publication in this collection
  16 Oct 2006
• Date of issue
  Sept 2006

History

• Received
  26 Sept 2005