Block scaling in the directed percolation universality class

The universal behaviour of the directed percolation universality class is well understood, both the critical scaling as well as finite size scaling. This article focuses on the block (finite size) scaling of the order parameter and its fluctuations, considering (sub-)blocks of linear size l in systems of linear size L. The scaling depends on the choice of the ensemble, as only the conditional ensemble produces the block-scaling behaviour as established in equilibrium critical phenomena. The dependence on the ensemble can be understood by an additional symmetry present in the unconditional ensemble. The unconventional scaling found in the unconditional ensemble is a reminder of the possibility that scaling functions themselves have a power-law dependence on their arguments.


Introduction
The directed percolation (DP) universality class comprises a huge number of non-equilibrium critical phenomena. Janssen [1] and Grassberger [2] famously conjectured more than 25 years ago that under very general circumstances, all models with a unique absorbing state (AS) belong to the DP universality class. While this conjecture has been confirmed numerically many times, evidence for the presence of the DP universality class in natural systems is still very scarce [3] (but see [4]).
One problem when probing field data for the presence of DP is that field data are more readily obtained in a single measurement rather than as a time series. However, the statistical features to be identified require an entire ensemble of realizations of the process in question. Instead of using a time series, one can resort to sub-sampling, i.e. splitting a large sample of size L d into (L/l) d small ones of size l d . For example, a population pattern obtained by measuring the spatial distribution of species could be split into several distinct blocks and their mutual correlations analysed.
The question 'How does the order parameter of an equilibrium system at the critical point change with the block size it is averaged over?' has been studied in great detail by Binder [5]. In the present work, a corresponding analysis is applied to models belonging to the DP universality class, more specifically to the contact process and to DP itself from the point of view of AS phase transitions. It turned out that the block averaged order parameter needs to be defined very carefully in order to reproduce standard finite size scaling (FSS). In the following, I will present numerical evidence and theoretical arguments that block FSS in DP can be very different from what is expected from equilibrium critical phenomena, depending on the choice of the ensemble. This observation can be readily applied to the analysis of field data, and will be illustrated using surrogate data. 3 cannot escape. In a finite system, the AS is reached with finite probability from anywhere in phase space, so that every (finite) system eventually becomes inactive, lim t→∞ ρ(t) = 0, where t measures the time in the model. However, the order parameter ρ signals a phase transition in a temperature-like tuning parameter p, dividing the parameter space into a region where the decay of the activity ρ is exponentially fast in time, from a region where for sufficiently large systems it is practically impossible to observe ρ = 0.
At least two (seemingly) different methods have been devised to overcome the problem that strictly lim t→∞ ρ(t) = 0 and to obtain the phase transition even in finite systems: one either introduces an external field that induces activity in the system [6] and analyses the model in the limit of arbitrarily small fields or considers the quasi-stationary state [7]. In this latter approach, after initializing and discarding a transient, all averages are taken conditional to activity. The activity entering the observables therefore never vanishes and the order parameter is always nonzero, rendering for example moment ratios well defined. One can show that both methods produce asymptotically equivalent scaling results [8]. In the present study, the quasi-stationary state was used to produce individual samples.
Above the critical point p c , i.e. for p > p c , the ensemble average of the activity ρ in the thermodynamic limit picks up as a power law ρ = A( p − p c ) β , where A is the amplitude, p c is the critical value of the tuning parameter and β is a universal critical index. Similarly, fluctuations of the order parameter scale as where B + and B − are the amplitudes above and below the transition, respectively, L is the linear extent of the system, d is its spatial dimension and γ is an independent critical exponent. Both these power laws are asymptotes and thus acquire corrections [9] away from the critical point.
Another kind of FSS can be explored in addition to the ordinary one just described. Instead of considering the entire system, observables are recorded within small blocks of linear extent l. To this end, I introduce the observables ρ u (l; L) and σ 2 u (ρ)(l; L), which are first and second cumulants of the activity within those blocks of linear size l in a system of size L. The blocks are produced by dividing samples of linear size L generated in the quasi-stationary state into (L/l) d blocks of linear size l each. This length l is conveniently chosen so that it divides L. By construction, each sample of size L d contains at least one active site and, consequently, at least one block of size l d has non-vanishing activity. However, there may be up to (L/l) d − 1 inactive blocks.
In addition to these observables, I introduce ρ c (l; L) and σ 2 c (ρ)(l; L), which are the first and second cumulants of the activity conditional to activity within the respective block. That means that inactive blocks are discarded when averaging. To distinguish the two ensembles, the former (with cumulants ρ u (l; L) and σ 2 u (ρ)(l; L)) will be called 'unconditional' (subscript u) and the latter (with cumulants ρ c (l; L) and σ 2 c (ρ)(l; L)) 'conditional' (subscript c). One can derive the moments of the unconditional ensemble from the corresponding moments in the conditional ensemble and vice versa by re-weighting, because they differ only by a number of samples with vanishing block activity, which are discarded in the conditional ensemble but not in the unconditional one (see [8] for similar considerations for the overall activity). If the fraction of active blocks averaged over the entire unconditional ensemble is a(l; L), then the nth moment of the unconditional ensemble ρ n u is related to the nth moment of the conditional ensemble ρ n c by ρ n u (l; L) = a(l; L) ρ n c (l; L) for n > 0. More than 25 years ago, block scaling was investigated by Binder [5] for the Ising model, where the order parameter is the magnetization density m rather than the activity ρ. In these systems with a symmetric phase space, there is no corresponding distinction between active and inactive blocks. Naïvely transferring these results to DP suggests ρ u,c (l; in leading order of L and in the limit of l being large compared with a lower cutoff l 0 , i.e. l l 0 . Below this constant threshold l 0 , i.e. l < l 0 , the order parameter deviates from the behaviour predicted by equation (2). This phenomenon is known and well understood in classical critical phenomena [12,14]. In the following, the data for l l 0 are not shown; in the two-dimensional (2D) contact process and DP l 0 was estimated to be ≈8 and in the 1D contact process l 0 ≈ 16. The dimensionless scaling functions G u,c (x) in equation (2) are bounded from above and usually also from below (away from 0) for all x ∈ [0, 1]. The scaling function is universal up to a pre-factor, which can always be absorbed into the metric factor C u,c . The latter is required for dimensional consistency.
Scaling behaviour similar to equation (2) is expected to hold for higher order moments and cumulants as well, in particular for the variance of the activity using the scaling relations cited above. Again, the scaling functions F u,c are constrained by being bounded from above and (usually) away from zero from below. So far, I have described what is expected in AS systems as inferred from equilibrium critical phenomena. In the remainder of this paper, I will first introduce the non-equilibrium models and the methods used in the numerical simulations for this study. The results for the different observables (activity, its variance and a moment ratio) in the different ensembles (unconditional and conditional) are then discussed in the light of analytical arguments. The paper finishes with an application to surrogate data, a discussion of the implications and the wider context of the findings and concludes with a brief summary.

Results
Most of the numerics in this work is based on the contact process [13,15] on a 2D square lattice with periodic boundary conditions, but the same results are found for the 1D contact process as well as for 2D site-directed percolation [16]. In fact, it will be argued that they are general features of the DP universality class, if not of all AS phase transitions. In the 2D contact process, occupied sites turn empty with extinction rate 1, whereas empty sites become occupied with rate zλ, with z being the fraction of occupied nearest neighbours. The critical value of λ has been estimated numerically with great accuracy λ c = 1.64877(3) [17] (I used λ = 1.6488 in the present study). The timescale is set by the extinction rate. The 2D contact process belongs to the DP universality class, which is characterized by exponents β = 0.583(3), γ = 0.297 (2) and ν ⊥ = 0.733(4) [18], so that β/ν ⊥ = 0.795 (6).
In site-directed percolation in 2 + 1 dimensions (BCC lattice), the time evolves discretely [16,19]. A site is occupied in the following time step with probability p if at least one of its directed neighbours is occupied, otherwise it is empty. The directed neighbours of a site are four sites in the preceding time step: the site itself, its right and upper nearest neighbours and its upper right next nearest neighbour. The critical value of p in this model has been estimated as p c = 0.34457(1) [18]. This model belongs to the 2D DP universality class as well.
All measurements are taken at the quasi-stationary state: starting from random initial configurations with a small but non-vanishing activity, the systems evolve according to the rules described above, until the observables reach a (quasi-)stationary state. For example, in the 2D contact process (CP) for L = 256, the first 10 5 updates are discarded as transient. Measurements are taken at constant rate after the transient and enter with the same weight until the system reaches the AS (for L = 256 the average lifetime was about 4.8 × 10 4 ) or a maximum time is reached (5 × 10 5 for L = 256). The procedure is repeated until the statistical error is acceptably small; for example 1.42 × 10 8 systems of size L = 256 were started, of which only about 13% survive the transient. Statistical errors were estimated by sub-sampling the ensemble, which copes even with correlated data. In the following, the scaling of the various (block) observables in l and L is analysed.

Order parameter
The first moment of the activity in the unconditional ensemble, ρ u (l; L), does not vary in l at all, because of translational invariance: every (randomly chosen) site is equally likely to be active and therefore ρ u (l; L) is constant in l. In order to establish ordinary FSS complying with equation (2), the scaling function G u (x) must necessarily be a power law itself, otherwise standard FSS, ρ u (l = L; L) ∝ L −β/ν , would not be recovered for l = L, i.e. when a single block covers the entire system. Contrary to what is expected from equilibrium critical phenomena, this scaling function necessarily vanishes at 0, i.e. lim x→0 G u (x) = 0. Figure 1 contains a (trivial) data collapse for ρ u (l; L) according to equation (2), namely ρ u (l; L)l β/ν ⊥ ∝ (l/L) β/ν ⊥ as a function of l/L for various system sizes L, which illustrates the scaling behaviour. The conditional order parameter ρ c (l; L), in contrast, does not suffer from this complication. By construction at least one site per patch is active, ρ c (l; L) l −d , so that lim L→∞ ρ c (l; L) l −d . The latter limit is the thermodynamic limit of a density and its existence is the most basic assumption in statistical mechanics. Provided l is sufficiently large compared with the fixed lower cutoff l 0 , so that equation (2) applies, the limit implies that lim x→0 G c (x) > 0, i.e. G c (x) converges to a nonzero value. This is numerically confirmed by the data collapse of ρ c (l; L) in figure 1.
However, strictly this argument does not apply, because l = 1 cannot be expected to be large compared with the lower cutoff l 0 (and in fact is not in the systems studied numerically in this paper). On the other hand, one might argue that σ 2 u (ρ)(l; L)/σ 2 u (ρ)(l = 1; L) can be expected to remain finite in the thermodynamic limit. While this is not a necessity, the alternative would imply a rather exotic behaviour of the variance, with the dotted line (the 'apparent asymptote') in figure 2 moving further and further away from σ 2 u (ρ)(l = 1; L) (which necessarily scales like L −β/ν ⊥ , the dashed line in figure 2) with increasing L. If σ 2 u (ρ)(l; L)/σ 2 u (ρ)(l = 1; L) converges and does not asymptotically vanish in L, then σ 2 u (ρ)(l; L) inherits the scaling of σ 2 u (ρ)(l = 1; L) in L and F u (x) = F u (1)x β/ν ⊥ for small x, so that σ 2 u (ρ)(l; L) = D u F u (1)l −β/ν ⊥ L −β/ν ⊥ , for small l/L.  (3)) of the conditional ensemble, whereas the sloped dashed line has slope β/ν ⊥ ≈ 0.795, expected to be the asymptote of the unconditional ensemble. The dotted line (apparent asymptote) with slope 0.681, however, fits the numerical data very convincingly, indicating that the asymptotic regime has not yet been reached.
scaling of σ 2 u (ρ)(l; L) in l, at a given, fixed L would then suggest σ 2 u (ρ)(l; L) ∝ l ((−2β/ν ⊥ )+0.681) . This is a reminder that scaling assumptions like equation (3) can numerically be verified only by a data collapse. The effective exponent of 0.681 is of course not a universal quantity and its deviation from 0.795 simply indicates that the asymptote has not been reached.
Not much can be said about the variance in the conditional ensemble. While the second moment has a lower bound (namely ρ 2 c l −d ), indicating that its scaling function does not vanish in the limit of small arguments, no lower bound exists for the variance; in fact, σ 2 c (ρ)(l = 1, L) = 0 by construction. It is, however, reasonable to assume that the variance σ 2 c (ρ)(l, L) is finite in the thermodynamic limit for fixed l, because by construction every patch always retains some activity regardless of the system size. In contrast, in the unconditional ensemble, the moments of the activity within a finite fraction of patches might vanish for a duration which increases with increasing system size, because for samples to continuously contribute to the average only one site needs to be active somewhere in the system. The scaling function F c (x) for the conditional ensemble being asymptotically finite, lim x→0 F c (x) > 0 is in line with the numerical evidence, see figure 2.

Active fraction a(l; L)
The scaling of the various observables is linked by equation (1). The fraction of active blocks, a(l; L), is given by any ratio of moments taken in the unconditional and the conditional ensemble, so that based on the first moments, its scaling is given by see equations (2) and (4). As can be seen from equation (1), if a moment in the conditional ensemble scales like l −nβ/ν , in the unconditional ensemble it will scale like l −(n−1)β/ν L −β/ν . Similarly, if σ 2 u (ρ) = D u F(1)l −β/ν L −β/ν for small l/L, the variance in the conditional ensemble is dominated by a term proportional to (l/L) −2β/ν .
One could consider the fraction of active blocks a(l; L) as the coarse-grained order parameter within a real-space renormalization group scheme [20], so that (L/l) d is the number of coarse-grained sites. The scaling of such an order parameter is proportional to (L/l) −β/ν , consistent with equation (5).

Moment ratio R
All the numerical results shown in figures 1 and 2 are for the 2D contact process. Based on FSS in equilibrium critical phenomena, one would normally expect appropriate moment ratios of the (absolute) order parameter, such as ρ 2 u,c (l; L)/ ρ 2 u,c (l; L), to be universal functions of l/L for l l 0 and to converge to a finite value for l/L → 0. Using equations (2) and (3), the moment ratio is universal assuming universality of the scaling functions and of the amplitude ratios. However, from what has been said earlier, in the unconditional ensemble, R u (l; L) ∝ (l/L) β/ν ⊥ , whereas in the conditional ensemble, R c (l; L) indeed converges to a finite value. Because R c (l; L) is universal, different models belonging to the same universality class, such as the CP and DP, should produce the same values. That is indeed the case, as shown in figure 3. For l = L, the moment ratios based on the unconditional and conditional ensembles coincide by construction and can be compared with the value of 0.7543(4) [21] (table VIII of [21], using R c (L; L) = 1/(1 + (K 2 /m 2 1 )) found by Dickman and Kamphorst Leal da Silva (see the dotted line in figure 3). The deviation of R c (l; L) from R c (L; L) for decreasing l/L does not mean that the former is not universal, just like one would generally expect that the value of the latter depends on various geometrical and topological properties of the lattice, such as its aspect ratio and the type of boundary condition [12].
The numerical situation for R u (l; L) (not shown) remained somewhat unclear. Only very few points for l/L close to 1 seem to overlap within the error. The system sizes simulated did not allow a firm statement as to whether R u (l; L) is actually universal.

Surrogate data
The proposed method can be put to test and compared with others using surrogate data, that is, data generated in a computer implementation of, for example, the contact process, mimicking real-world data as they would be obtained in physical or biological systems. In contrast with the simulation data used above, such data consist of a single realization, as if one were to analyse a satellite image or field data. Figure 4 shows the conditional box scaling of the activity ρ c (l; L) for a single realization of a 2D system of linear size L = 128 as a function of the box size l. The slope of ρ c (l; L) is compared with l −β/ν ⊥ (thick dashed line) and fits very well. No error bar can be given for l = L, as there is only one such sample, the remaining error bars are estimated from the variance of the corresponding sample of size (L/l) 2 , assuming independence 2 .
For comparison, figure 5 is a double logarithmic plot of the two-point correlation function as a function of the absolute distance |r|. Here, I (r) indicates occupation of site r and the sums run over all L 2 sites r . It is taken from the same sample as figure 4, using the translational invariance and the eightfold symmetry of a square. Nevertheless, the data are comparatively noisy. Naïve scaling arguments along the lines of equilibrium phase transitions [22] suggest an asymptote C(r) ∝ |r| −d+2−η = |r| −2β/ν ⊥ . In figure 5, it is shown as a thick dotted line, but this asymptote is compatible with the data only in a narrow, noisy intermediate regime.  The correct scaling behaviour of the correlation function, however, is C(r) ∝ |r| −β/ν ⊥ , which is found by imposing that C(r)/ ρ u does not scale in L [23] (as it would, for example, in equlibrium phase transitions, consistent with [24]). Although still noisy, the data for small r shown in figure 5, indeed, are much better compatible with a scaling exponent of −β/ν ⊥ = −0.795 (6).
Box scaling of the order parameter is due to the correlations captured in the two-point correlation function, as can be seen directly by deriving the variance σ 2 u (ρ)(l; L) from C(r), consistent with σ 2 u (ρ)(l; L) ∝ l −β/ν ⊥ as observed earlier. Box scaling therefore can be regarded as an elegant form of extracting the scaling from the correlation function. From that point of view, the advantage of box scaling over a direct investigation of the two-point correlation function is merely due to its simplicity: box scaling is well understood theoretically and easy to implement in an experiment, in field work or in a computer simulation.
As a final remark, the quality of the surrogate data shows a strong time dependence. Starting from a randomly occupied lattice, choosing too short an equilibration time leads to a lack of correlations, with the system still being dominated by the independent, random initialization. Waiting too long, on the other hand, means that the system is likely to be low in activity and just about to die out completely, producing sparse and biased results.

Discussion
A block scaling analysis could provide a practical method to overcome the problem of limited availability of data in natural systems and allow the analysis of a natural system's scaling behaviour without the need for an entire time series. Block scaling effectively is a form of sub-sampling and the analysis utilizes the universal scaling with changing block size, which is characterized in the present work. As it turns out, in order to obtain block scaling as known from equilibrium critical phenomena, one has to use a conditional ensemble, where moments of the activity in a block enter the average only conditional to non-vanishing activity. The situation in AS phase transitions therefore is very different from what is expected from equilibrium critical phenomena, where no additional condition is needed in order to obtain box scaling corresponding to FSS.
Box scaling effectively measures the correlations between finite boxes throughout the system. In (near-)equilibrium systems, i.e. systems with a Hamiltonian, at the critical point, the probability density function of the order parameter is symmetric around 0 due to the symmetry of the Hamiltonian. As contributions with opposite signs cancel, fluctuations decrease with increasing box size and, consequently, the box averaged (absolute) order parameter and its variance decline.
In AS systems, this mechanism does not exist: the local order parameter is non-negative and therefore cannot cancel out. In an unconditional ensemble the order parameter does not change with block size and fluctuations about the mean are not symmetric. For the order parameter to vanish, its fluctuations must vanish as well. Nevertheless, introducing the conditional ensemble restores the behaviour of (near-)equilibrium systems.
Within the conditional ensemble, the scaling functions of the observables considered converge to a finite value for small arguments, so that the cumulants scale in l with the expected exponents, ρ c ∝ l −β/ν ⊥ and σ 2 c (ρ) ∝ l −2β/ν ⊥ for sufficiently small l/L. Moreover, the moment ratio ρ 2 c / ρ 2 c is universal and converges to a nonzero, universal value for sufficiently small l/L. This is not the case in the unconditional ensemble: the scaling functions are power laws themselves and therefore the moments display an unconventional scaling. As a consequence, the average activity is constant in l, whereas its variance scales like σ 2 u (ρ) ∝ l −β/ν ⊥ L −β/ν ⊥ for small l/L. In addition, the moment ratio ρ 2 u / ρ 2 u vanishes asymptotically in small l/L. The unconventional scaling of the unconditional ensemble is due to the existence of the translational symmetry, which is not present in the conditional ensemble. This symmetry causes the lack of scaling of the first moment, which is connected to that of all other moments through equation (1).
The ratio of any unconditional moment and its conditional counterpart is the fraction of active blocks a(l; L); this quantity itself displays universal behaviour in the ratio l/L. Any unconditional moment can be derived from the corresponding conditional one and vice versa, by multiplying and dividing by a(l; L), respectively.
The key advantage of the conditional ensemble is the scaling of the first moment ρ c , which is usually the easiest to determine, carrying the smallest statistical error. In the unconditional ensemble, the first moment does not scale at all. Moreover, the moment ratio ρ 2 c / ρ 2 c converges to a finite, universal value, which again is a quantity with a comparatively small statistical error. The hope is to use these observables in experimental situations where the system is guaranteed to be at the critical point.