Directed percolation conjecture for cellular automata revisited

The directed percolation (DP) hypothesis for stochastic, range-4 cellular automata with acceptance rule $y \le\sum_{j=-4}^4 s_{i-j} \le 6$, in cases of $y<6$ was investigated in one and two dimensions. Simulations, mean-field renormalization group and coherent anomaly calculations show that in one dimension the phase transitions for $y<4$ are continuous and belong to the DP class for $y=4,5$ they are discontinuous. The same rules in two dimensions for $y=1$ show $2+1$ dimensional DP universality; but in cases of $y>1$ the transitions become first order.


I. INTRODUCTION
The phase transition theory of nonequilibrium systems is an area of considerable interest. There is not such a clear view about critical universality classes as there is for equilibrium systems. The cellular automata (CA) approach for modelling nonequilibrium systems can be regarded as a general tool like differential equations. It introduces a space and time discretization -thus making the solution treatable -but preserving the features of interest. As was shown in [1] CA models can describe very rich behavior. Even the simplest one dimensional CA can have various phase transitions, can generate chaotic patterns, or show complex self-organizing or reproductive behavior.
The directed percolation (DP) hypothesis [2] states that all continuous transitions of a scalar order parameter to a unique absorbing state in dimension d belong to the DP universality class, represented by the d + 1 dimensional directed percolation [3] or reggeon field theory. There are some additional conditions -usually not stated explicitly -to be fulfilled as well : short range interaction both in space and time, translational invariance, absence of multicritical points, and the non-vanishing probability for any active state to die locally. Although there is no proof for this conjecture there exists a great number of models, which have been found to belong to this class. Attempts to find counterexamples [4][5][6] by other simulation results [7] have been unsuccessful. The motivation of our investigation is to find a possible exception to the DP conjecture in the cellular automata phase transition of different dimensions. The stochastic cellular automata -investigated by Bidaux, Boccara and Chatè (hereafter referred to as BBC CA) [5] and Jensen [7] using simulations -in one dimension showed controversial results concerning the universality class. Steady state simulation gave non-DP whereas the time-dependent simulation resulted in the DP class. Even though the latter method is considered to be more precise, arguments based on universality suggested that the same rule in two dimensions should produce a transition of 2+1 dimensional DP class. Simulations proved this to be a first order transition [5]. In order to clarify the situation we used various selection of different analytical and simulation methods.
In addition we were interested in the change of universality by varying the lower acceptence value y of the rule : where X ∈ {0, 1} is a two valued random variable such that P rob(X = 1) = p .
By setting up the mean-field steady state equation for the probability of ones (c) : the leading term in the polinomial function f , which determines the critical behavior as p → p c (and c → 0) is c for y = 1, c 2 for y = 2 and c 3 for y = 3. In analogy with the equilibrium Ising model -where the three-spin Ising model has a different universality from the two-spin nearest neighbor Ising model -we expected deviations from the DP behavior.
To check this assumption we extended our investigation to the y = 1, 2, 3 cases in one and two dimensions.

II. STEADY STATE SIMULATIONS
The steady state Monte Carlo simulations (STDMC) were carried out on L = 40000 and L = 256 × 256 sized lattices with periodic boundary conditions. In each case measurements at 30-50 different p values were taken. Following 10 5 -5 · 10 5 relaxation lattice updates at each p the concentration averages and r.m.s. values were determined. The exponent β and p c were fitted together -using the weighted least squares error method -to the critical scaling formula Since the concentration in the immediate neighborhood of the critical point has a relatively large statistical error -owing to the increased relaxation time -data points within the 1% neighborhood of p c were discarded. In order to compensate this truncation of data we performed the fitting in several different neighborhoods of the p c . Since we don't know how the (2) scaling breaks down away from p c , we have checked this tendency of the fitting by varying p cut of the used data points : (p c × 1.01) < {p} < p cut . This kind of analysis enabled us to see how the critical exponent converged as p cut → p c and we were able to extrapolate for the true β exponent in spite of the growing statistical errors in the vicinity of p c .

A. Results in one dimension
The generally accepted 1+1 dimensional DP critical exponets are summarized in Table  I. Figure 1 shows that the estimated critical exponent of the order parameter converges to the value β ≈ 0.28, which is about the DP value determined by Dickman and Jensen by series expansion [8]. The critical point and β estimates ( given in Table II) were determined by linear weighted least-square fitting to the functions of p cut . The results of the BBC [5] (y = 3) model are the most uncertain, but still seem to extrapolate to the DP class value if p cut → p c . Estimates based on fitting with less than 10 data points should not be taken very seriously and are shown for information only.
As a tendency the transitions become more and more steep (and p c → 1) as we increase the lower limit of the CA rule (y). For y = 4 and 5 non-empty steady state exists in the deterministic, p = 1 limit only. One can easily check that only a few oscillating structures can survive. There is a first order transition to the empty state if p < 1, which is quite rare among one dimensional CA. For y > 5 every cell dies out.
To verify our results further simulations and analytical calculations were carried out; these are discussed in the subsequent sections.

B. Results in two dimensions
In two dimensions for y = 1 the transition seems to be continuous and the estimates for the critical exponent converged to the 2 + 1 dimensional DP class value : β ≈ 0.58 [9] (Figure 1/d). For y = 2 and 3 the transitions appear to be discontinuous, see Figures 2. If formula (2) is fitted to the concentration data -as if continuous transitions were assumedresults in p c -s with values lower than p with non-zero densities. Since it is not easy to prove first order transition by simulation, we applied the generalized mean-field method, which has already been found to be very useful for deciding the order of the transition at CA [10]. The results of these calculations are also shown in Figure 2 and they are discussed in section V.
C. Finite-size scaling analysis Finite-size scaling is shown to be applicable to continuous transition to an absorbing state of nonequilibrium systems [11]. At the critical point the steady state density (c) and the fluctuation χ = L d (< c 2 > − < c > 2 ) scale with the system size as : where ν ⊥ is the correlation length exponent in the space direction : Simulations were done in one dimension for lattice sizes : L = 16, 32, 64, ..., 8192 with the corresponding timesteps : t = 200, 400, 800, ..., 102400. Averaging was performed on the 'surviving' samples out of N s = 262144, 131072, 65536, ..., 512 independent runs. The p c values are taken from the time-dependent MC calculations. Figure 3 shows that the c and χ results fall nicely on a straight line on a log-log plot, but for y = 1 the fitting resulted in: for y = 2 : β/ν ⊥ = 0.252(1), γ/ν ⊥ = 0.548 (7), and for y = 3 : which differ slightly from the DP class values [12]: We have shown that the deviations come from systems with small sizes, similarly as for [12]. By omitting points with sizes less than L = 16, 32, 64 the results converge to the DP values (see Tables III, IV, V) and we can see the emergence of a DP like behavior in the large L limit.

III. TIME-DEPENDENT SIMULATIONS
Time-dependent Monte Carlo simulations (TDMC) have been shown to be an efficient method for locating critical points and estimating exponents for models with continuous transition to an absorbing state [12][13][14]. The initial state of this simulation is the smallest a single live 'seed', that can survive at the center of the lattice. The evolution of the population is followed up to several thousands of lattice update steps and the following quantities are measured at each time step : • survival probability s(t) • concentration c(t) • average mean square distance of spreading from the center R 2 (t) The evolution runs are averaged over N s independent runs for each different value of p in the vicinity of p c ( but for R 2 (t) only over the surviving runs ). At the critical point we expect these quantities to behave in accordance with the power law as t → ∞, i.e.
We carried out simulations for the one-dimensional systems up to t = 5000 time steps with N s ≃ 10 6 independent samples. The seed was a single 1 for y = 1, a neighboring pair of 1-s for y = 2, and a triplet for y = 3. We did not follow the evolution of R 2 (t). To estimate the critical exponents we determined the local slopes : using m = 8. Figures 4, 5, 6 show these quantities as functions of t −1 . In the case of powerlaw behavior we should see a straight line as t −1 → 0, when p = p c . The off-critical curves should possess curvature. Curves corresponding to p > p c should veer upward, curves with p < p c should veer down. The critical exponent should be read off as the interception of the critical (straight) curve with the ordinate axis. It is emphasized that it is very difficult to get precise estimates on the critical data by visual inspection of the curves . While for y = 3 the middle curve -which seems to be the nearest to the critical one -one can read off DP exponents, for y = 1, 2 one could conclude that we have non-DP values. The logarithmic c(t) and s(t) data were fitted by parabolas : ln(s(t)) = const. + δ 1 ln(t) + δ 2 ln 2 (t) (12) at different p values. Estimations for p c are based on the conditions η 2 = 0 and δ 2 = 0; the results are shown in Table II. It should be noted, however that even though the difference is small (∆p < 10 −4 ) the conditions η 2 = 0 and δ 2 = 0 do not give exactly the same p c . Also the estimates for the exponents are slightly off the DP class values. The errors must have statistical nature with finite size errors being excluded since the fronts of the populations -over the entire run -cannot exceed the boundary of the lattice.
To get more precise exponents we used the general form -including corrections to scaling -of the local slopes [13] where δ , is a correction to the scaling exponent. As indicated in Table II the exponents agree well with the DP class data.

IV. MEAN-FIELD RENORMALIZATION GROUP ANALYSIS
Renormalization group methods for nonequilibrium systems have been successful when combined with mean-field approximation. At the critical point the scaling behavior of quantities calculated on small clusters of different sizes is exploited. One main branch of such calculations is the Coherent Anomaly Method (CAM) its application to these CA models will be discussed in section VI. Another important approach is the so-called mean-field renormalization group (MFRG) method, originally introduced by Indekeu et al. [15] for equilibrium systems. Later it was shown that the method may also be applied to self-avoiding walk and percolation problems [16]. The application to nonequilibrium systems was demonstrated by [17]. Although this method is not considered to be the most accurate, it is relatively simple and effective. Developments to enhance its performance are under way [18]. Here, we utilize the approach for one-dimensional CA.
First we set up mean-field equations for different sized block probabilities (P i ) in the steady state. Here the index i denotes an n block. However we do not use the Bayesian extension trick as in the case of the generalized mean-field calculations, but take into account the external field with the average concentration (c): Because in the vicinity of the second-order phase transition to the absorbing state c → 0, we expand and keep linear terms only in the equations. We express the probability of a site being occupied (O n (p, c)) by calculating it in terms of different (n) cluster sizes. For the y = 1 model in the n = 1 case the equations are : In view of this, the probability of a site being in state 1 is : Expanding in factors of c we get : For clusters of pairs (n = 2) the equations look like : Now, the probability of 1 after linearization is : At the critical point the scaling requires that the O n (p, c) probability must scale like the external field : This condition enables one to combine the expressions for O n and O n , and to get an estimate for the p c . Comparing the n = 1 and n = 2 results we obtained p c (1, 2) = 0.123106. The exponent ν ⊥ , which measures the critical behavior of the correlation length, can be obtained from the scaling of (p − p c ) resulting from change of the length scale l = n , /n : For n = 1 and n = 2 clusters it is : ν ⊥ (1, 2) = 5.802. This is very far from the DP class ν ⊥ ∼ 1.1 value, but if calculations are done on larger clusters the results converge there. We were able to go up to n = 5 clusters with the help of symbolic MATHEMATICA programming. Table VI shows that there is a slow convergence of the p c and ν ⊥ values. With increasing the cluster sizes we expect more accurate results. Unfortunately we were not able to apply the method for y = 2, 3 models because the mean-field approximations predict first order transitions for small cluster sizes and thereafter no scaling is expected.

V. GENERALIZED MEAN-FIELD RESULTS
The generalized mean-field approximation (GMF) introduced by [19,20] has already been applied to describe phase transitions of different stochastic CA models [21,10] with success. The phase structure and the tricritical point of the one-dimensional BBC model (y = 3) with site-exchange mixing has been explored with the help of GMF [10].
The traditional mean-field approximation has been extended to handle clusters of sizes n > 1 with the help of the Bayesian extension. Correlations of range > n at a given level of approximation (n) are neglected and the convergence of the results is analysed by increasing n; the details of the technique are given in [21]. The extension of the method to the two-dimensional CA model is straightforward. In this case the nearest-and next-nearest neighbor of the lattice point can contribute to the sum defined in (1). Taking into account the probability of a given distribution we used maximal overlap of the clusters. At threepoint level the form of the cluster breaks the x − y symmetry and prevents the configuration from being covered symmetrically. Because of this the result of the equations depends on the way of covering the configuration. To avoid this uncertainty we restricted ourselves to symmetrical clusters (n = 1, 2, 4).

A. Results in one dimension
In the case of the y = 1 model the convergence of the approximations is monotonic. Even the traditional (n = 1) mean-field approximation predicts continuous phase transition. By applying quadratic fitting to the p c (n) results we got an extrapolation value of n → ∞ value, which agrees with the simulation result (see Table VII). For y = 2, 3 the low level approximations gave a discontinuous phase transition, but the gap sizes decrease with increasing n and finally disappear at n = 6 for y = 2, and at n = 3 for y = 3. It is interesting that the convergence is non-monotonic for the y = 3 case in contrast to what we expect for equilibrium systems.

B. Results in two dimensions
As one expects in two dimensions the GMF approximation gives more precise results than in one dimension. For the n = 4 level approximation the p c data differ from the simulations by a few per cent only. For the y = 1 model all approximations predict a second order phase transition, whereas for y = 2, 3 cases the transitions are discontinuous (see Table VIII) and the gap sizes increase with n. We can expect that in the n → ∞ exact limit the transition remains first order in agreement with the STDMC simulations.

VI. COHERENT ANOMALY EXTRAPOLATION
The GMF approximations give slow convergence at the critical point since the correlation length there goes to infinity. Therefore the accuracy of method breaks down in the immediate vicinity of p c . To obtain estimates for the critical exponents one needs to apply an extrapolation technique. This technique usually is based on finite size scaling theory. The mean-field renormalization technique is one such approach.
The coherent anomaly method (CAM) introduced by Suzuki [22] is based on the scaling relation that the solution for singular quantities at a given (n) level of approximation (Q n (p)) in the vicinity of the critical point is the product of the classical singular behavior multiplied by an anomaly factor a(n), which becomes extremely large as n → ∞ (and p n c → p c ): where p is the control parameter and ω cl is the classical critical index. The divergence of this anomaly factor scales as thereby permitting the estimation of the true critical exponent ω, given a set of GMF approximation solutions. The method has been tested on a great number of solved and unsolved systems [23]. The application to cellular automata was shown by [24]. We followed the technique described there in the case of the present CA models. A new parametrization is used suggested by [25] that has an invariance property: p ↔ p −1 . We took into account a correction term to the anomaly scaling: The GMF calculation gave enough data to do the CAM extrapolations for the one dimensional y = 1 case (Table IX). In a similar way to the MFRG case for y = 2, 3 models the first order predictions of the low level GMF results excluded CAM calculations. Table X shows that β results are stable in contrast to neglecting points from the nonlinear fitting procedure by Eq.(26) and indicating DP class behavior.

VII. CONCLUSIONS
A deatiled analysis for a family of totalistic, critical CA in one and two dimensions has been carried out in order to check the DP conjecture and to clarify the former contradictory simulation results of the BBC (y = 3) model. Different critical exponents were determined by simple steady state, finite size scaling and time dependent simulations. The simulation results are more extensive than the previous ones and none of them seems to break the DP rule. Rather they all give estimates for different critical exponents that are close to the DP universality values. We wish to emphasize that careful analysis is necessary for each kind of simulation to get unambiguous results. We have shown that the systematic calculation of GMF provides a reliable way to distinguish the order of phase transitions. Mean-field renormalization group and coherent anomaly method extrapolations to GMF provided estimates for critical exponents. These results are again in accordance with the DP conjecture.