text:
# Quasi-stationary distributions for stochastic processes with an absorbing state

## I Introduction

Stochastic processes with an absorbing state arise frequently in statistical physics , epidemiology , and related fields. On one hand we have autocatalytic processes (in lasers, chemical reactions, or surface catalysis, for example), and population models (in genetic or epidemiological modeling), in which a population of $`n0`$ individuals goes permanently extinct when the absorbing state, $`n=0`$, is reached. On the other hand, phase transitions to an absorbing state in spatially extended systems , exemplified by the contact process , are currently of great interest in connection with self-organized criticality , the transition to turbulence , and, more generally, issues of universality in nonequilibrium critical phenomena . Such models are fequently studied using deterministic mean-field equations, Monte Carlo simulation, and renormalization group analyses.

While fluctuations play an essential role in the vicinity of the absorbing state, exact solution of the master equation is not, in general, feasible. A useful tool in this situation is van Kampen’s $`\mathrm{\Omega }`$-expansion in inverse powers of the system size, about the deterministic, macroscopic solution. It is clearly desirable to develop additional methods for analyzing processes with absorbing states. In this work we study models whose macroscopic equation admits a nonabsorbing (or active) stationary state, but which, for a finite system size, must always end up in the absorbing state. We show that the quasi-stationary distribution for such a system is readily constructed, and that it provides a wealth of information about its behavior. In fact, simulations of “stationary” properties of lattice models with an absorbing state actually study the quasi-stationary regime, since the only true stationary state (for a finite system) is the absorbing one.

For models without spatial structure, such as uniformly distributed populations, well stirred chemical reactors, or networks in which each node communicates equally with all others, our results provide a complete description of long-time properties, conditioned on survival. In models with spatial structure, typified by nearest-neighbor interactions on a lattice, as in the contact process, the description in terms of one or a few variables corresponds to a mean-field theory that cannot adequately capture critical fluctuations. Our approach does, nonetheless, allow one to put some of the fluctations “back into” mean-field theory, and so we are able to study finite size effects and moment ratios that are beyond the grasp of simpler approaches. For one-step models involving a single variable, the computational demands of our approach are trivial; more complicated processes can be analyzed using straightforward numerical procedures.

The idea of the quasi-stationary state or distribution is quite simple. Consider a Markov chain $`n_t0`$ in continuous time, with $`n=0`$ absorbing, and such that the macroscopic or rate-equation description (i.e., neglecting fluctuations) admits an active stationary state ($`n_{st}>0`$). Then in many cases one expects the probability distribution $`p_n(t)`$ to bifurcate, after some initial transient, into two components, one with all weight on the absorbing state, the other, $`q_n(t)`$, concentrated near the macroscopic value $`n_{st}`$, such that $`p_n0`$ for $`n0`$. The survival probability (for the process not to have fallen into the absorbing state), is $`P(t)=_{n1}q_n(t)`$. We study the case in which $`q_n`$ has attained a time-independent form, in which its only time dependence is through the overall factor $`P_s(t)`$. Since the probability distribution is supposed to bifurcate, we may write

$$p_n(t)=Q(t)\delta _{n,0}+q_n(t)$$
(1)
where $`Q(t)=1P(t)`$, and by definition, $`q_0(t)=0`$. The quasi-stationary distribution (QS) is then defined via the condition:

$$p_n(t)=P(t)\overline{p}_n,(n1),$$
(2)
where the $`\overline{p}_n`$ are time independent. (Again, $`\overline{p}_00`$.) Since $`P(t)`$ is the survival probability we have the normalization:

$$\underset{n1}{}\overline{p}_n=1.$$
(3)
Clearly, not every process with an absorbing state admits a QS distribution. The exponential decay process (with transition rates $`W_{n,m}=m\delta _{n,m1}`$) and the birth-and-death process ($`W_{n,m}=m\delta _{n,m1}+m\lambda \delta _{n,m+1}`$), for example, do not. While we will not try to determine rigorously the conditions under which a QS distribution exists, it seems reasonable to suppose that the macroscopic equation for the process in question admits an active stationary solution. All of the examples considered satisfy this condition. The voter model represents a special case, in which the macroscopic equation provides no information on the fate of the process, and the QS distribution bears no resemblance to the stationary (absorbing) distribution.

The balance of this paper is organized as follows. In Sec. 2 we illustrate the evaluation of the QS distribution, and the analysis of associated statistical properties, for the simple case of the contact process on a complete graph. Sec. 3 is concerned with the closely related Malthus-Verhulst process (MVP). Both the contract process and the MVP exhibit a continuous phase transition in the infinite-size limit. In Secs. 4 and 5 we study a models with a discontinuous phase transition: the second Schlögl process, and the voter model. In Sec. 6 we return to the contact process, this time in a two-site approximation, and in Sec. 7 study a pair of competing MVPs. Sec. 8 contains a summary and discussion of our results.

## II Contact process on a complete graph

In the contact process (CP) , each site of a lattice is either occupied ($`\sigma _i(t)=1`$), or vacant ($`\sigma _i(t)=0`$). Transitions from $`\sigma _i=1`$ to $`\sigma _i=0`$ occur at a rate of unity, independent of the neighboring sites. The reverse can only occur if at least one neighbor is occupied: the transition from $`\sigma _i=0`$ to $`\sigma _i=1`$ occurs at a rate of $`\lambda m`$, where $`m`$ is the fraction of nearest neighbors of site $`i`$ that are occupied; thus the state $`\sigma _i=0`$ for all $`i`$ is absorbing. The stationary order parameter $`\rho `$ (the fraction of occupied sites), is zero for $`\lambda <\lambda _c`$. It is easy to show that in mean-field approximation (i.e., treating occupancy of sites as statistically independent events, and supposing that the density $`\rho (t)`$ is spatially uniform), the density follows the rate equation:

$$\frac{d\rho }{dt}=(\lambda 1)\rho \lambda \rho ^2.$$
(4)
This equation predicts a continuous phase transition (from $`\rho 0`$ to $`\rho =1\lambda ^1`$ in the stationary state) at $`\lambda _c=1`$. While this is qualitatively correct, $`\lambda _c`$ is in general larger than 1; $`\lambda _c=3.29785(2)`$ in one dimension, for example.

A mean-field description such as Eq. (4) does not, of course, yield correct critical exponents for dimension $`d<d_c=4`$, but that is not of concern here. Instead, we would like to study the simplest stochastic process for which Eq. (4) represents the macroscopic limit, (i.e., the limit of an infinite system, in which $`\rho (t)`$ is a deterministic variable). By studying such a process we can recover the fluctuations ignored in the macroscopic equation. Since the process has an absorbing state, its limiting ($`t\mathrm{}`$) distribution is trivial ($`p_n(t)\delta _{n,0}`$), but we can study the long-time properties, conditioned on survival, using the QS distribution.

The stochastic birth-and-death process whose macroscopic limit is given by Eq. (4) is the contact process on a complete graph, in which the rate for a vacant site to turn occupied is $`\lambda `$ times the fraction of all sites that are occupied, rather than the fraction of nearest neighbors. Since each site interacts equally with all others, all pairs of sites are neighbors, defining a complete graph. Given its connections with epidemic modeling and percolation, this model has also been called “percolitis” .

The state of the process is specified by a single variable $`n`$: the number of occupied or infected sites. This is a one-step process with nonzero transition rates

$$W_{n1,n}=n$$
(5)
$$W_{n+1,n}=\lambda \frac{n}{\mathrm{\Omega }}(\mathrm{\Omega }n)$$
(6)
on a graph of $`\mathrm{\Omega }`$ sites. These expressions yield the master equation

$$\frac{dp_n}{dt}=(n+1)p_{n+1}+(n1)\overline{\lambda }(\mathrm{\Omega }n+1)p_{n1}\left[1+\overline{\lambda }(\mathrm{\Omega }n)\right]np_n,$$
(7)
where $`\overline{\lambda }\lambda /\mathrm{\Omega }`$. Since $`P(t)=_{n1}p_n(t)`$, we have $`dP/dt=p_1(t)`$ in the contact process, and under the QS hypothesis,

$$\frac{1}{P}\frac{dP}{dt}=\overline{p}_1.$$
(8)
If we insert Eq. (2) in the master equation, then, dividing through by $`P(t)`$ we obtain

$$\overline{\lambda }(n1)[\mathrm{\Omega }n+1]\overline{p}_{n1}+(n+1)\overline{p}_{n+1}[\overline{\lambda }(\mathrm{\Omega }n)+1]n\overline{p}_n+\overline{p}_1\overline{p}_n=0,$$
(9)
for $`n1`$. For a given $`n`$, Eq. (9) furnishes a relation for $`\overline{p}_{n+1}`$ in terms of $`\overline{p}_n`$ and $`\overline{p}_{n1}`$. Letting $`u_n=n[\overline{\lambda }(\mathrm{\Omega }n)+1]`$, we have

$$\overline{p}_2=\frac{1}{2}\overline{p}_1(u_1\overline{p}_1)$$
(10)
and

$$\overline{p}_n=\frac{1}{n}\left[(u_{n1}\overline{p}_1)\overline{p}_{n1}+(n2)(\mathrm{\Omega }2)\overline{\lambda }\overline{p}_{n2}\right]$$
(11)
for $`n=2,\mathrm{},\mathrm{\Omega }`$. Thus the QS distribution is completely determined once $`\overline{p}_1`$ is known; the latter is determined via the normalization condition, Eq. (3).

In practice the following iterative scheme converges very quickly. Starting with a guess for $`\overline{p}_1`$ (unity, for instance), one uses the recurrence relations to find the corresponding $`\overline{p}_n`$, and then evaluates $`S=_{n=1}^\mathrm{\Omega }p_n`$. The procedure is then repeated using $`\overline{p}_1^{}=\overline{p}_1/S`$; after a modest number of steps, $`S`$ converges to unity, at which point the $`\overline{p}_n`$ represent the QS distribution.

We have confirmed that the master equation for the contact process on a complete graph does in fact attain the QS distribution at long times. We integrate the master equation using a fourth-order Runge-Kutta scheme , and stop the integration when the mean population conditioned on survival, $`n_s=_nnp_n(t)/P(t)`$, changes by less than $`10^{10}`$ per time step. The resulting distribution, $`p_n^{}=p_n(t)/P(t)`$ is compared with the QS distribution furnished via the recurrence relations in figure 1; they are identical.

Having constructed the QS distribution and verified that it indeed represents the long-time distribution (conditioned on survival), as given by the master equation, we now examine some properties of the QS state. Figure 2 shows the quasi-stationary population density $`\rho =n/\mathrm{\Omega }`$ versus $`\lambda `$ for various system sizes. This plot looks strikingly similar to finite-size plots of the order parameter at a continuous phase transition. In particular, it is evident that the $`\lambda `$-dependence is smooth for any finite system, but that on increasing $`\mathrm{\Omega }`$, the function $`\rho (\lambda ,\mathrm{\Omega })`$ approaches a singular limit.

On finite-dimensional lattices the value of the order parameter at the critical point is expected to scale as $`\rho (\lambda _c,\mathrm{\Omega })\mathrm{\Omega }^{\beta /\nu _{}}`$, where $`\beta `$ and $`\nu _{}`$ are the critical exponents governing the order parameter and the divergence of the correlation length . Here we find $`\rho (\lambda _c,\mathrm{\Omega })0.685\mathrm{\Omega }^{1/2}`$ (from results for $`\mathrm{\Omega }=50010^5`$). We note that the mean-field exponent values $`\beta =1`$ and $`\nu _{}=1/2`$ for the contact process would lead one to expect $`\rho (\lambda _c,\mathrm{\Omega })\mathrm{\Omega }^2`$. The reason for the difference would seem to be that on a complete graph (where each site is, in effect, a unit distance from all others), the correlation length is undefined, and so finite-size scaling ideas do not apply. Instead, we can understand the exponent -1/2 by noting that the microscopic equation (4), gives $`\rho =0`$ at the critical point, so that $`n`$ is purely a fluctuation. The central limit theorem then requires $`n\mathrm{\Omega }^{1/2}`$. Away from the critical point, $`\lambda =1`$, the finite-size correction to the mean-field solution, $`\rho =1\lambda ^1`$, is $`𝒪(\mathrm{\Omega }^1)`$.

Another property of interest at a continuous phase transition is the moment ratio $`m=n^2/n^2`$. This quantity is analogous to Binder’s reduced fourth cumulant , at an equilibrium critical point: the curves $`m(\lambda ,\mathrm{\Omega })`$ for various $`\mathrm{\Omega }`$ cross near $`\lambda _c`$ (the crossings approach $`\lambda _c`$), so that $`m`$ assumes a universal value $`m_c`$ at the critical point. For the basic contact process (and other models in the universality class of directed percolation), $`m_c=1.174`$ and 1.326 in one and two dimensions, respectively . In figure 3 we plot $`m`$ for the quasi-stationary contact process on a complete graph. Qualitatively, the behavior is similar to that observed on finite-dimensional lattices, with $`m(\lambda )`$ becoming ever steeper with increasing system size; $`|dm/d\lambda |_{\lambda _c}\mathrm{\Omega }^{1/2}`$ for large $`\mathrm{\Omega }`$. The crossings approach $`\lambda _c`$, as expected. We find that $`m(1,\mathrm{\Omega })1.660+𝒪(\mathrm{\Omega }^{1/2})`$ for large $`\mathrm{\Omega }`$.

As noted above, the survival probability decays as $`dP(t)/dt=p_1(t)`$. We may therefore interpret $`\overline{p}_1`$ as the quasi-stationary decay rate. The QS lifetime, $`\tau =1/\overline{p}_1`$, grows as $`\mathrm{exp}[\text{const.}\mathrm{\Omega }(\lambda \lambda _c)]`$ for $`\lambda >\lambda _c`$, as can be seen from figure 4, where we plot $`\mathrm{\Omega }^1\mathrm{ln}\tau `$ versus $`\lambda `$ for several system sizes. At the critical point, $`\tau \mathrm{\Omega }^{1/2}`$.

It may appear surprising that we are able to discuss quasi-stationary properties for $`\lambda <\lambda _c`$, where an active stationary state does not exist, even in the thermodynamic limit. The existence of nontrivial QS properties for $`\lambda <\lambda _c`$ is in fact a finite-size effect, analogous to a nonzero magnetization in the Ising model above the critical temperature, on a finite lattice. For $`\lambda >\lambda _c`$ the QS state converges, as $`\mathrm{\Omega }\mathrm{}`$, to the true stationary state, while for $`\lambda <\lambda _c`$ it converges to the absorbing state. But since the properties of any finite system are nonsingular, we must expect a smooth decay of, for example, the density, as $`\lambda 0`$.

### A Pseudo-stationary distribution

The contact process on a finite graph (be it a regular lattice, a complete graph, or small-world network) has only the absorbing state as its true stationary state. It is nonetheless instructive to examine the result of forcing a time-independent solution on the master equation. This is easily done by introducing the generating function $`F(z,t)_np_n(t)z^n`$. The master equation, Eq.(7), is equivalent to the partial differential equation,

$$\frac{F}{t}=(1z)\left[(1\overline{\lambda }z)\frac{F}{z}+\overline{\lambda }z\frac{}{z}\left(z\frac{F}{z}\right)\right].$$
(12)
Setting $`F/t=0`$ and letting $`G(z)=F/z`$, a first integration yields

$$G(z)=Cz^{\mathrm{\Omega }1}e^{1/\overline{\lambda }z},$$
(13)
where $`C`$ is a constant. Integrating once more gives

$$F_{ps}(z)=C\underset{n=0}{\overset{\mathrm{}}{}}\frac{z^{\mathrm{\Omega }n}}{\overline{\lambda }^nn!(\mathrm{\Omega }n)}.$$
(14)
While the corresponding distribution,

$$p_{ps,m}=\frac{C\overline{\lambda }^m}{(\mathrm{\Omega }m)!m}$$
(15)
yields $`dp_m/dt=0`$ when inserted in the master equation, it is not a valid solution since it cannot be normalized, and does not satisfy the correct boundary condition at $`m=0`$. For this reason we call $`p_{ps,m}`$ a pseudo-stationary distribution. Similar pseudo-stationary distributions have been derived by Muñoz from the Fokker-Planck equation for continuous versions of the contact process and of a model with multiplicative noise .

If we treat $`p_{ps,m}`$ as a proper probability distribution, restricting it to $`m=1,\mathrm{},\mathrm{\Omega }`$, and normalizing it on this set, then the resulting distribution turns out to be very close to the quasi-stationary one when $`p_m`$ is negligible in the vicinity of $`m=0`$. Figure 5 (for $`\mathrm{\Omega }=100`$) shows that the two distributions are essentially identical for $`\lambda =2`$ (as they are for larger $`\lambda `$ values), and for $`\lambda =1.5`$, except for small $`n`$. For $`\lambda =1.2`$, however, the distributions are radically different; in particular, we see that the pseudo-stationary distribution does not respect the absorbing boundary at $`n=0`$.

In summary, the pseudo-stationary distribution provides a simple and useful approximation to the true QS distribution for $`\lambda `$ (and $`\mathrm{\Omega }`$) sufficiently large that $`p_m0`$ for $`m1`$. (In this case the survival probability decays extremely slowly, rendering the hypothesis of a stationary distribution - while strictly incorrect - at least reasonable.) As we approach the critical value, however, the region near $`n=0`$ assumes dominance, and $`p_{ps}`$ bears little relation to the quasi-stationary distribution.

## III Malthus-Verhulst process

The Malthus-Verhulst process (MVP) is a birth-and-death type process, $`n(t)`$, very similar to the contact process on a complete graph, but with the difference that $`n`$ may assume any nonnegative integer value. The nonzero transition rates for this one-step process are:

$$W_{n1,n}=\mu n+\frac{\nu }{\mathrm{\Omega }}n(n1),$$
(16)
$$W_{n+1,n}=\lambda n.$$
(17)
In the first expression, $`\mathrm{\Omega }`$ represents the system size, as discussed by van Kampen . In the MVP there is no lattice structure, and no fixed limit to the population size, as there is in the contact process; rather, growth is limited by the term $`\nu `$, which represents competition between individuals for access to finite resources. (Thus $`\nu `$ represents an intrinsic competition parameter, and $`(\nu /\mathrm{\Omega })n(n1)`$ is the rate of interactions between pairs of individuals in a system of size $`\mathrm{\Omega }`$.) The saturation effect that permits a nontrivial quasi-stationary state is imposed in the death term, whereas it appears in the birth term, in the contact process.

Letting $`n=\mathrm{\Omega }\rho `$, the macroscopic equation for the MVP is

$$\frac{d\rho }{dt}=(\lambda \mu )\rho \nu \rho ^2$$
(18)
which has the stationary solution $`\overline{\rho }=(\lambda \mu )/\nu `$. The master equation is

$$\frac{dp_n}{dt}=[\mu +\overline{\nu }n](n+1)p_{n+1}+(n1)\lambda p_{n1}[\lambda +\mu +\overline{\nu }(n1)]np_n,$$
(19)
with $`\overline{\nu }=\nu /\mathrm{\Omega }`$. Inserting the quasistationary form of the probability distribution, one readily obtains the recurrence relation

$$\overline{p}_n=\frac{(q_{n1}\mu \overline{p}_1)\overline{p}_{n1}\lambda (n2)\overline{p}_{n2}}{n[\mu +\overline{\nu }(n1)]},$$
(20)
where

$$q_n=[\lambda +\mu +\overline{\nu }(n1)]n$$
(21)
The QS distribution can be found using the same iterative procedure as for the contact process. (We cut off the distribution at a certain $`M`$ such that $`\overline{p}_{}10^{20}`$ for $`n>M`$.)

The quasi-stationary properties of the MVP are quite similar to those of the contact process on a complete graph. Figure 6 shows the dependence of $`\rho `$ on $`\lambda `$ for system sizes $`\mathrm{\Omega }=100`$, 10<sup>3</sup> and 10<sup>4</sup>. (We fix $`\mu =\nu =1`$ in all numerical studies of the MVP.) The density grows linearly with $`\lambda `$ for $`\lambda >\lambda _c=\mu `$, unlike in the contact process, where the macroscopic saturation term is also proportional to the birth rate, $`\lambda `$. As in the contact process on a complete graph, the lifetime $`\tau \mathrm{\Omega }^{1/2}`$ at $`\lambda _c`$, and the QS density decays $`\mathrm{\Omega }^{1/2}`$.

The moment ratios $`m`$ (figure 6, inset) show the same qualitative trend as in the CP. The asymptotic value of $`m`$ at the critical point, moreover, is the same as in the CP: we find $`m1.6600.814\mathrm{\Omega }^{1/2}`$ for large $`\mathrm{\Omega }`$. This identity of $`m`$ values suggests that the limiting QS probability distributions at $`\lambda _c`$ are the same for the two processes. This is verified numerically; for $`\mathrm{\Omega }=10^3`$, for example, the distributions differ by less than about $`2\times 10^4`$.

In fact, the critical QS probability distribution for the MVP (and by extension for the contact process) enjoys the scaling property

$$\overline{p}_n=\frac{1}{\sqrt{\mathrm{\Omega }}}f(n/\sqrt{\mathrm{\Omega }}).$$
(22)
The scaling function $`f`$ may be found using a method that parallels van Kampen’s analysis of the master equation . First, we treat $`n`$ in Eq. (20) as a continuous variable, expanding $`\overline{p}_{n\pm 1}`$ to second order. Setting $`\lambda =\mu `$ at the critical point, we obtain

$$\frac{1}{2}[(2\mu +\overline{\nu })n+\overline{\nu }n^2]\frac{d^2\overline{p}}{dn^2}+[2\mu +\overline{\nu }n(1+n)]\frac{d\overline{p}}{dn}+2\overline{\nu }n\overline{p}=\mu \overline{p}_1\overline{p}.$$
(23)
Next we make the change of variables $`y=n/\mathrm{\Omega }^{1/2}`$, and $`f(y)=\mathrm{\Omega }^{1/2}\overline{p}_n`$. This yields an expansion in inverse powers of $`\mathrm{\Omega }^{1/2}`$; at lowest order ($`\mathrm{\Omega }^1`$) we find

$$\mu y\frac{d^2f}{dy^2}+[2\mu +\nu y^2]\frac{df}{dy}+2\nu yf=\mu f(0)f(y).$$
(24)
The initial value $`f(0)`$ must be chosen such that $`_0^{\mathrm{}}f𝑑y=1`$.

To integrate Eq. (24) numerically we let $`g(y)=df/dy`$, and note that $`g(0)=[f(0)]^2/2`$, and $`dg/dy|_{y=0}=(2\nu /\mu )f(0)`$. Setting $`\mu =\nu =1`$, numerical integration yields the scaling function shown in figure 7, which agrees well with the distributions for the MVP and the contact process on a complete graph, in the limit of large system size. \[Note that the change of variables used to derive Eq. (24) yields the very same equation when applied to the recursion relation for the contact process, Eq. (9), at its critical point, $`\lambda =1`$.\] The asymptotic scaling function $`f(y)`$ yields the moment ratio $`m=1.660063`$, again in agreement with results from the recurrence relations. (Note that $`m=f(0)/y`$, as may be verified directly from Eq. (24).)

Despite some superficial differences, the MVP and the contact process on a complete graph show the same type of phase transition (continuous, with mean-field-like critical exponents), in the limit $`\mathrm{\Omega }\mathrm{}`$, and the same scaling properties at their respective critical points.

As in the case of the contact process, one may define a pseudo-stationary distribution by demanding that the r.h.s. of the master equation be zero. The resulting expression,

$$p_{ps,n}=C\frac{(\lambda /\overline{\nu })^n}{n(n+\mu /\overline{\nu }1)!}$$
(25)
is very similar to that of the contact process on a complete graph, and again represents a good approximation to the quasi-stationary distribution if $`\overline{p}_n0`$ for $`n1`$.

## IV Schlögl’s Second Model

The processes studied in the preceding sections exhibit a continuous phase transition; next we consider a process with a discontinuous transition. Schlögl’s second model  describes the set of autocatalytic reactions $`A0`$, $`2A3A`$, and $`3A2A`$ in a well stirred system. Since the growth process is quadratic in the density (rather than linear, as in the contact and MV processes), a low-density active state is unstable, and we expect a discontinuous transition to a state with a nonzero density as the growth rate is increased. The nonvanishing transition rates are:

$$W_{n1,n}=\mu n+\frac{\nu }{\mathrm{\Omega }^2}n(n1)(n2),$$
(26)
$$W_{n+1,n}=\frac{\lambda }{\mathrm{\Omega }}n(n1).$$
(27)
These yield the macroscopic equation:

$$\frac{d\rho }{dt}=\mu \rho +\lambda \rho ^2\nu \rho ^3,$$
(28)
which admits an active stationary solution,

$$\rho _s=\frac{1}{2}\left[\lambda +\sqrt{\lambda ^24\mu \nu }\right],$$
(29)
for $`\lambda 2\sqrt{\mu \nu }`$ (The stationary density jumps from zero to $`\sqrt{\mu \nu }`$ at the transition. It is known that in a spatially extended system the transition is actually continuous for $`d<4`$ , but here we treat the well stirred system, which exhibits a discontinuous phase transition.)

Since this is a one-step process, it is a simple matter to derive the recurrence relations for the quasi-stationary probabilities:

$$\overline{p}_n=\frac{(q_{n1}\mu \overline{p}_1)\overline{p}_{n1}\overline{\lambda }(n2)(n3)\overline{p}_{n2}}{n[\mu +\overline{\nu }(n1)(n2)]}$$
(30)
where

$$q_n=[\mu +\overline{\lambda }(n1)+\overline{\nu }(n1)(n2)]n,$$
(31)
with $`\overline{\lambda }=\lambda /\mathrm{\Omega }`$ and $`\overline{\nu }=\nu /\mathrm{\Omega }^2`$. These may be solved using the same method as for the contact process or the MVP.

In the following numerical example, we fix $`\mu =\nu =1`$, and study the neighborhood of the transition at $`\lambda =2`$. Figure 8 shows the quasi-stationary density near the transition; with increasing $`\mathrm{\Omega }`$ the density approaches the discontinuous stationary solution, Eq. (29). The density goes from $`𝒪(1/\mathrm{\Omega })`$ to $`𝒪(1)`$, in a transition region of width $`\mathrm{\Omega }^{1/2}`$. The moment ratio $`m`$, at the transition, appears to approach unity (from above) slowly, as $`\mathrm{\Omega }\mathrm{}`$; data for $`\mathrm{\Omega }5\times 10^4`$ yield $`m1+2.64\mathrm{\Omega }^{0.265}`$. Similarly, the density at the transition approaches its limiting value from below: $`\rho 11.5\mathrm{\Omega }^{0.211}`$. The lifetime at the transition appears to grows as $`\tau \mathrm{\Omega }^{0.35}`$. Essentially the same power laws (but with different prefactors) are found at the transition for $`\mu =2`$ and $`\nu =1/2`$.

The probability distribution is bimodal in the neighborhood of the transition, a hallmark of a discontinuous transition. This is evident in figure 9, where the relative amplitudes of the two peaks, one at $`n=1`$, the other near $`n1000`$, shift rapidly near the transition. The bimodal distribution of course presages a hysteresis loop when $`\mathrm{\Omega }\mathrm{}`$.

## V Voter model on a complete graph

A somewhat simpler model exhibiting a discontinuous phase transition is the voter model . It differs from the systems considered previously in that it possesses two absorbing states. In the voter model each site on a lattice is in one of two states, $`A`$ or $`B`$; each $`AB`$ nearest-neighbor pair has a unit rate to change to the state $`AA`$ (with probability $`\omega `$) or to $`BB`$ (with probability $`1\omega `$). Let $`n(t)`$ denote the number of sites in state $`A`$; both $`n=0`$ and $`n=\mathrm{\Omega }`$ are absorbing, on a lattice of $`\mathrm{\Omega }`$ sites. Starting from a random initial condition with equal densities of $`A`$ and $`B`$, the stationary probability distribution switches discontinuously from $`\delta _{n,0}`$ to $`\delta _{n,\mathrm{\Omega }}`$ as $`\omega `$ is increased through 1/2. Of particular interest is the coarsening dynamics for $`\omega =1/2`$ . In this case $`n`$ is constant, and the probabilities of the two absorbing states are determined by $`n(0)`$.

In this brief section we study the model on a complete graph, for the case $`\omega =1/2`$. The master equation is

$$\frac{dp_n}{dt}=\frac{n1}{2\mathrm{\Omega }}(\mathrm{\Omega }n+1)p_{n1}+\frac{n+1}{2\mathrm{\Omega }}(\mathrm{\Omega }n1)p_{n+1}\frac{n}{2\mathrm{\Omega }}(\mathrm{\Omega }n)p_n.$$
(32)
This yields the macroscopic equation $`d\rho /dt=0`$, reflecting the dominance of fluctuations in this process. From the terms for $`n=0`$ and $`n=\mathrm{\Omega }`$ we find that in the QS state, the survival probability obeys

$$\frac{1}{P}\frac{dP}{dt}=\frac{1}{2}(1\mathrm{\Omega }^1)(\overline{p}_1+\overline{p}_\mathrm{\Omega }),$$
(33)
allowing us to write the following relation for the QS distribution:

$$\overline{p}_n(1\mathrm{\Omega }^1)(\overline{p}_1+\overline{p}_{\mathrm{\Omega }1})=\frac{n1}{\mathrm{\Omega }}(\mathrm{\Omega }n+1)\overline{p}_{n1}+\frac{n+1}{\mathrm{\Omega }}(\mathrm{\Omega }n1)\overline{p}_{n+1}2\frac{n}{\mathrm{\Omega }}(\mathrm{\Omega }n)\overline{p}_n.$$
(34)
The solution is the uniform distribution, $`\overline{p}_n=1/(\mathrm{\Omega }1)`$ for $`n=1,\mathrm{},\mathrm{\Omega }1`$ (here, $`\overline{p}_0=\overline{p}_\mathrm{\Omega }=0`$).

In the case of the voter model, the QS distribution looks very different from the stationary distribution, even in the infinite-$`\mathrm{\Omega }`$ limit. This is because the process has no active stationary state: on a complete graph it must always fluctuate into one of the absorbing states. If we exclude all trials that hit one of the boundaries, we are in effect looking at a random walker confined to an interval, whose distribution is asymptotically uniform.

## VI Two-site approximation for the contact process

In this section we return to the contact process, constructing a stochastic description in terms of nearest-neighbor pairs of sites, on a ring of $`\mathrm{\Omega }`$ sites. The idea of the approximation is similar to that of two-site or pair mean-field dynamic mean-field theory , except that in the present case we treat the number of occupied sites, $`n`$, and the number of doubly-occupied nearest-neighbor pairs, $`z`$, as stochastic variables. When necessary, we invoke the usual factorization of three-site probabilities in terms of those for two sites, in order to construct the transition rates.

We begin by establishing the range of allowed values for $`z`$. Using ‘1’ and ‘0’ to represent occupied and vacant sites, respectively, $`z`$ is the number of (11) nearest-neighbor (NN) pairs. Let $`w`$ represent the number of (10) NN pairs. (For obvious reasons, the number of (01) pairs is again $`w`$.) Note that $`w`$ is not an independent variable, since each 1 is followed wither by a 0 or another 1, yielding $`n=z+w`$. Similarly, the number of (00) pairs $`v`$ is given by $`v=\mathrm{\Omega }2n+z`$. (Here we used $`n+2w+z=\mathrm{\Omega }`$.) Evidently $`v`$ is restricted to the set $`\{0,\mathrm{},\mathrm{\Omega }\}`$ while $`w`$ must be in $`\{0,\mathrm{},\mathrm{\Omega }/2\}`$ (We assume $`\mathrm{\Omega }`$ even.) These considerations imply certain limits for $`z`$ on a ring of $`\mathrm{\Omega }`$ sites, as listed in table 1.

Table 1. Allowed values for $`z`$ in the CP on a ring.    $`n`$ $`z`$    0, 1 0  $`2,\mathrm{},\mathrm{\Omega }/2`$ $`0,\mathrm{},n1`$  $`\mathrm{\Omega }/2+1,\mathrm{},\mathrm{\Omega }1`$ $`2n\mathrm{\Omega },\mathrm{},n1`$  $`\mathrm{\Omega }`$ $`\mathrm{\Omega }`$

Next we construct the transition rates $`W_{n^{},z^{};n,z}`$. From a given state $`(n,z)`$ there are at most five possible transitions, viz., to the states $`(n^{},z^{})=(n1,z)`$, $`(n1,z1)`$, $`(n1,z2)`$, $`(n+1,z+1)`$, and $`(n+1,z+2)`$. The first three represent annihilation of a particle with, resp. zero, one or two occupied nearest neighbors, while the last two represent creation of a particle at a site with one, or two occupied neighbors, resp. The rate of the transition $`(n,z)(n1,z)`$ is the number of isolated particles, that is, triplets of the form (010). This number is not determined completely by $`n`$ and $`z`$; we estimate it as the number of (01) pairs times the conditional probability of a vacant site given an occupied neighbor: $`w\times (w/n)`$. This estimate, however, is obviously wrong when $`z=n1`$ and $`n2`$, in which case there are no (010) triplets. Thus we have

$$W_{n1,z;n,z}=\{\begin{array}{cc}0,\hfill & z=n1,n>1\hfill \\ \frac{(nz)^2}{n},\hfill & \text{ otherwise }\hfill \end{array}$$
(35)
By similar arguments, we find:

$`W_{n1,z1;n,z}`$ $`=`$ $`\{\begin{array}{cc}2,\hfill & z=n1,n>1\hfill \\ (nz)\frac{2z}{n},\hfill & \text{ otherwise }\hfill \end{array}`$ (38)
$`W_{n1,z2;n,z}`$ $`=`$ $`\{\begin{array}{cc}0,\hfill & z=1\hfill \\ n2,\hfill & z=n1,n=3,\mathrm{},\mathrm{\Omega }1\hfill \\ \frac{z^2}{n},\hfill & \text{ otherwise }\hfill \end{array}`$ (42)
$`W_{n+1,z+1;n,z}`$ $`=`$ $`\{\begin{array}{cc}\lambda ,\hfill & z=n1\hfill \\ \lambda (nz)\frac{\mathrm{\Omega }2n+z}{\mathrm{\Omega }n},\hfill & \text{ otherwise }\hfill \end{array}`$ (45)
$`W_{n+1,z+2;n,z}`$ $`=`$ $`\{\begin{array}{cc}0,\hfill & z=n1\hfill \\ \lambda \frac{(nz)^2}{\mathrm{\Omega }n},\hfill & \text{ otherwise }\hfill \end{array}`$ (48)
In the macroscopic limit, the master equation defined by these rates yields a pair of coupled equations for the particle density $`\rho `$ and the pair density $`\zeta =z/\mathrm{\Omega }`$ . The critical value $`\lambda _c=2`$ in this approximation.

Since this is not a one-step process, we cannot obtain the QS distribution via recurrence relations, as was done in the previous examples. The most obvious approach - direct integration of the master equation - works reasonably well for modest system sizes ($`\mathrm{\Omega }<200`$ or so), but becomes very slow for large $`\mathrm{\Omega }`$. This is due in part to the large number of equations ($`\mathrm{\Omega }^2`$), and to the fact that the transition rates are $`𝒪(\mathrm{\Omega })`$, forcing one to use a time step $`\mathrm{\Omega }^1`$, to avoid numerical instability. A much more efficient iterative scheme that converges directly to the QS distribution is described in Ref. . The basic idea of this approach is that in the QS regime,

$$\overline{p}_n=\frac{r_n}{w_nr_0},$$
(49)
for any nonabsorbing state $`n`$, where $`r_n=_n^{}W_{n,n^{}}p_n^{}(t)`$ is the flux of probability into state $`n`$ ($`r_0`$ is the flux into the absorbing state), and $`w_n=_n^{}W_{n^{},n}`$ is the total transition rate exiting $`n`$. The iterative scheme  starts with an arbitrary distribution $`p_n`$ normalized on $`n1`$ (or, in the multivariate case, on the set of nonabsorbing states). Next one evaluates

$$p_n^{}=ap_n+(1a)\frac{r_n}{w_nr_0},$$
(50)
where $`a`$ is a parameter between zero and one. The new distribution $`p_n^{}`$ is then normalized, and inserted in the r.h.s. of Eq. (50), and the process repeated until it converges. For large $`\mathrm{\Omega }`$ the time to reach the QS distribution is three or more orders of magnitude smaller than that required for numerical integration. (In the present case convergence is enhanced if one uses a small value for $`a`$, e.g., 0.01. A further economy is realized by noting that, near $`\lambda _c`$, $`\overline{p}_{nz}`$ is extremely small for $`n\mathrm{\Omega }/2`$, allowing one to truncate the distribution. For $`\mathrm{\Omega }=1280`$ and $`\lambda =2`$, for example, $`\overline{p}_{nz}10^{17}`$ for $`n=350`$.)

We studied the QS properties of the contact process on a ring of $`\mathrm{\Omega }=10`$, 20, 40,…, 2560 sites, using the pair approximation. Figure 10 shows that the dependence of the particle density and the moment ratio on $`\lambda `$ and on system size is qualitatively similar to that found in the complete graph case. The scaling laws for the density ($`\mathrm{\Omega }^{1/2}`$) and the lifetime ($`\mathrm{\Omega }^{1/2}`$) at the critical point are the same as for the CP on a complete graph. The critical values of the moment ratio $`m`$ at $`\lambda _c`$, extrapolated to infinite system size, yield $`m_c=1.653`$, nearly the same as for the CP on a complete graph. The critical probability distribution $`\overline{p_n}`$ again shows a scaling collapse (figure 11), but the scaling function differs somewhat from that for the complete graph. Figure 12 shows the joint QS distribution $`\overline{p}_{n,z}`$ at the critical point, for $`\mathrm{\Omega }=40`$. For each $`n`$ there is a most probable number, $`z^{}(n)`$, of pairs; $`z^{}`$ grows linearly with $`n`$.

## VII Competing Malthus-Verhulst processes

As a final example we consider a population consisting of two types, $`A`$ and $`B`$, that evolve according to the MVP analyzed in Sec. 3. (The two subpopulations may, for example, possess different alleles of a certain gene, in a species with asexual reproduction.) $`A`$ and $`B`$ interact via two mechanisms: (1) each individual competes with all others (regardless of type) for access to resources; (2) on reproducing, $`A`$ mutates to $`B`$ with probability $`q`$, and vice-versa. Thus we have the following set of transition rates:

$$W_{n_A1,n_B;n_A,n_B}=\mu _An_A+\frac{\nu }{\mathrm{\Omega }}n_A(n_A+n_B1),$$
(51)
$$W_{n_A,n_B1;n_A,n_B}=\mu _Bn_B+\frac{\nu }{\mathrm{\Omega }}n_B(n_A+n_B1),$$
(52)
$$W_{n_A+1,n_B;n_A,n_B}=\lambda _Apn_A+\lambda _Bqn_B,$$
(53)
$$W_{n_A,n_B+1;n_A,n_B}=\lambda _Bpn_B+\lambda _Aqn_A,$$
(54)
where $`p=1q`$.

Evidently the state $`n_A=n_B=0`$ is absorbing. For mutation probability $`q=0`$, $`n_A=0`$ is also absorbing (similarly $`n_B=0`$), and in the QS state only one type is present, i.e., $`\overline{p}_{n_A,n_B}=0`$ if both $`n_A`$ and $`n_B`$ are nonzero. This prompts us to investigate the effect of a nonzero mutation rate: Will the predominance of one species persist when mutation is possible? To address this issue we study the order parameter

$$\mathrm{\Delta }\frac{|n_An_B|}{n_A+n_B}$$
(55)
in the QS state. Figure 13 shows how $`\mathrm{\Delta }`$ varies with $`q`$, for various system sizes. The competing MVPs are both at the critical point, $`\lambda _i=\mu _i`$. (For simplicity we take all rates, $`\lambda _i`$, $`\mu _i`$, and $`\nu `$, equal to unity). The order parameter decreases with increasing mutation probability, as expected, and with increasing system size. From analysis of $`\mathrm{\Delta }`$ at various $`q`$-values, as a function of $`\mathrm{\Omega }`$, it appears that $`lim_\mathrm{\Omega }\mathrm{}\mathrm{\Delta }(q,\mathrm{\Omega };\lambda _c)=0`$ for any nonzero $`q`$. In other words, an arbitrarily small mutation rate restores the symmetry of the process. We observe qualitatively similar behavior in $`\mathrm{\Delta }(q,\mathrm{\Omega })`$ above the critical point ($`\lambda >1`$). At the critical point, the decay of $`\mathrm{\Delta }(q,\mathrm{\Omega })`$ with $`\mathrm{\Omega }`$ is slower than exponential, but faster than a power law; it can be fit approximately by a stretched exponential $`\mathrm{exp}[c\mathrm{\Omega }^{0.2}]`$, where the factor $`c`$ depends on $`q`$ (see figure 13 inset).

## VIII Discussion

We study quasi-stationary probability distributions for Markov processes exhibing a phase transition between an active and an absorbing state. The QS distribution, which represents the long-time limit of the probability distribution, conditioned on survival, yields a complete statistical description of the process in this limit. While the QS distribution can always, in principle, be found by integrating the master equation, we present much more efficient methods for its generation: recurrence relations for one-step processes, and an iterative scheme for multi-step or multivariate processes.

Some insight into the nature of the quasi-stationary distribution is afforded by the considering the equation (in the notation of Sec. VI), for $`n1`$,

$$\frac{dq_n}{dt}=w_nq_n+r_n+r_0q_n,$$
(56)
where $`r_n=_n^{}W_{n,n^{}}q_n^{}`$ is the flux of probability into state $`n0`$. (As usual we take state 0 as absorbing.) Without the final term on the right-hand side, this is simply the master equation. Including this term, however, the equation has as its stationary solution the QS distribution, $`\overline{p}_n`$. In other words, the QS distribution corresponds to the stationary state of a “process” in which all flux to the absorbing state is immediately reinserted into the non-absorbing subspace. The portion alloted to state $`n`$ is equal to its probability, $`q_n`$.

Thus the QS distribution does not result from imposing a simple boundary condition (such as periodic, or reflection at the origin), on the original master equation. Since Eq. (56) is nonlinear, it is not a master equation. It may happen, nevertheless, that the QS distribution for a certain process (possessing an absorbing state), is also the stationary distribution for some other process, as in fact was found for the voter model.

In this work we study QS distributions for processes with an absorbing state, but without spatial structure. We find that the Malthus-Verhulst process and the contact process on a complete graph have the same scaling properties near their respective critical points. A study of Schlögl’s second model illustruates how a discontinuous phase transition emerges in the infinite-size limit, associated with a rapid change in the QS distribution, which is bimodal in the vicinity of the transition. The voter model on a complete graph serves as a negative example: it has no active stationary state, and the QS distribution, which is uniform on the set of allowed density values, is very different from the true long-time distribution, even in the infinite-size limit. Our two final examples are bivariate processes. The pair approximation to the contact process shows scaling properties that are very similar to those of the complete-graph version. A study of a pair of competing Malthus-Verhulst processes reveals that mutation causes a radical change in the QS distribution.

Several directions for future work on quasi-stationary distributions can be mentioned. First, it should be possible to derive continuum equations, analogous to the Fokker-Planck equation, to describe the QS distribution in the limit of large system size. Second, analysis of QS properties for a series of system sizes promises to be a useful method for studying lattice models possessing an absorbing state. Through finite-size scaling analyses and extrapolation procedures, it should be possible to extract useful estimates of critical properties for such models. Finally, numerical study of more elaborate models used in epidemiology and population dynamics may yield a better understanding of such processes.

Acknowledgment

We are grateful to Miguel Angel Muñoz for helpful discussions. This work is supported by CNPq, Brazil.

FIGURE CAPTIONS

Figure 1. Probability distributions for the contact process on a complete graph, $`L=100`$ sites. The curves (connecting values at integer points) represent the long-time limit of the master equation, normalized to the survival probability $`P(t)`$. Points represent the quasi-stationary distribution determined by the recurrence relations. Left curve: $`\lambda =1.2`$; right: $`\lambda =2`$.

Figure 2. Quasi-stationary population density versus $`\lambda `$ for the contact process on a complete graph. System sizes 100, 200, 500, 1000, and 2000, top to bottom on the left-hand side. The inset is a detail of the critical region, including the curve for the infinite-size limit, $`\rho 1\lambda ^1`$.

Figure 3. Quasi-stationary moment ratio $`m`$ versus $`\lambda `$ for the contact process on a complete graph. System sizes 20, 50, 100, 200, 500, and 1000, in order of increasing steepness.

Figure 4. Quasi-stationary lifetime $`\tau `$ versus $`\lambda `$ for the contact process on a complete graph. System sizes 100, 500, 1000, and 5000, upper to lower on left-hand side.

Figure 5. Quasi-stationary distribution (solid line) and pseudo-stationary distribution (points) for the contact process on a complete graph, $`L=100`$ sites for $`\lambda =2`$ (upper), 1.5 (middle), and 1.2 (lower). The inset in the middle panel is a detail of the region near $`n=0`$.

Figure 6. Quasi-stationary density versus $`\lambda `$ for the Malthus-Verhulst process. System sizes 100, 10<sup>3</sup>, and 10<sup>4</sup>, upper to lower on the left-hand side. The inset shows moment ratio $`m`$ versus $`\lambda `$ for the same system sizes.

Figure 7. Scaling plot of the quasi-stationary distribution, $`f=\mathrm{\Omega }^{1/2}\overline{p}_n`$ versus $`y=n/\mathrm{\Omega }^{1/2}`$, at the critical point: $`+`$: MVP, $`\mathrm{\Omega }=10^4`$; $``$: contact process on complete graph, $`\mathrm{\Omega }=10^3`$; solid curve: asymptotic scaling function $`f(y)`$ obtained via numerical integration.

Figure 8. Quasi-stationary density versus $`\lambda `$ for the second Schlögl process with $`\mu =\nu =1`$. System sizes 100, 500, 1000, and 5000, (in order of increasing steepness). The dashed curve is the mean-field solution, which is discontinuous at $`\lambda =2`$.

Figure 9. Quasi-stationary distribution for the second Schlögl process with $`\mu =\nu =1`$, system size $`\mathrm{\Omega }=1000`$. $`\lambda =1.98`$, 2, and 2.02 (upper to lower on left-hand side).

Figure 10. Quasi-stationary density versus $`\lambda `$ for the contact process in the pair approximation. System sizes 20, 40, 80, 160, 320, and 640. The inset shows moment ratio $`m`$ versus $`\lambda `$ for the same system sizes.

Figure 11. Scaling plot of the quasi-stationary distribution: $`P^{}=\mathrm{\Omega }^{1/2}\overline{p}_n`$ versus $`n^{}=n/\mathrm{\Omega }^{1/2}`$, at the critical point of the contact process. Open circles: pair approximation, $`\mathrm{\Omega }=640`$; solid line: pair approximation, $`\mathrm{\Omega }=2560`$; filled circles: contact process on complete graph, $`\mathrm{\Omega }=1000`$.

Figure 12. Quasi-stationary probability distribution for the critical contact process in the pair approximation, system size $`\mathrm{\Omega }=40`$.

Figure 13. Quasi-stationary order parameter $`\mathrm{\Delta }`$ versus mutation rate $`q`$ for a pair of competing MVPs at the critical point. System sizes $`\mathrm{\Omega }`$ = 20, 50, 100, 200, 500, 1000 (upper to lower). The inset shows $`\mathrm{ln}\mathrm{\Delta }`$ versus $`\mathrm{\Omega }^{0.2}`$ for (upper to lower) $`q=0.02`$, 0.05, and 0.1.


format: text

text:
# Photometric validation of a model independent procedure to extract galaxy clustersBased on observation collected at the European Southern Observatory, Chile, ESO N∘ 62.O-0230 and 64.O-0317(A).

## 1 INTRODUCTION

Clusters of galaxies are the largest virialized structures in the Universe, and accurate knowledge of their global properties is needed to constrain models of galaxy formation and evolution.
The first step in this direction requires the construction of a statistically well-defined sample of clusters in the nearby universe to be used as a ”local template”. Unfortunately, due to historical and observational reasons, and in spite of much effort, existing samples cannot be considered ideal. Existing cluster catalogs, in fact, fall into three main categories: i) large catalogs derived from photographic surveys (POSS-I, UKST) by visual inspection and covering wide portions of the sky (Abell (1958); Abell et al. (1989); Zwicky et al. 1961- (68)) but missing the needed depth, homogeneity and completeness (Postman et al. (1986); Sutherland (1988)); ii) catalogs machine extracted with objective criteria from photographic plates (cf. Dodd & MacGillivray (1986); Dalton et al. (1992); Lumsen et al. (1992)), reaching, in some cases, limiting magnitudes fainter than (i), but not covering equally wide areas of the sky and so far available only for the Southern hemisphere (only UKST plates); iii) accurate and deeper catalogs usually derived from CCD data and selected on the basis of objective criteria but covering much smaller regions of the sky and containing only a small number of objects (cf. Postman et al. (1996); Olsen et al. (1999)). Excluding the ones, which will be derived from the Sloan Digital Sky Survey (SDSS) (Kim et al. (2000); Kepner et al. (1999)), for the Northern sky no automatically extracted catalogs of galaxies (nor of clusters) could be produced until the recent completion of the Digital Palomar Sky Survey (DPOSS).
DPOSS, which covers three bands (J,F,N), is characterized by deeper limiting magnitudes than other catalogs extracted from photographic surveys and therefore offers an important opportunity to investigate galaxy clusters in the nearby (i.e. $`z<0.4`$) Universe. Recently, DPOSS photometric calibration (Weir, Djorgovski & Fayyad (1995)) and extraction of the Palomar Norris Sky Catalog (PNSC) were completed at Caltech in collaboration with the observatories of Rome, Naples and Rio de Janeiro, partners in the CRoNaRio project (Djorgovski et al. 1998a ; Andreon et al. (1997)). This catalog contains astrometric, photometric and morphological information for all objects detected down to limiting magnitudes of $`g_J21.5`$, $`r_F20.5`$ and $`i_N20.0`$ in the Gunn & Thuan photometric system.
With the availability of these new data, various methods to build galaxy cluster catalogs based on color information have been proposed (Gal et al. (1999); Gal et al. 2000a ).

In order to exploit the scientific potential of the DPOSS, we developed a model independent procedure to search for cluster candidates. The procedure is detailed elsewhere and therefore we briefly summarize its main characteristics here (Puddu et al. (2000); Puddu et al., in preparation). First, for a given region of the sky, we extract the individual catalogs obtained from the corresponding J, F and N DPOSS plates and calibrate them to the g, r and i Gunn-Thuan System using the procedure described in Weir et al. (1995). Then, after correcting for misclassified objects in at least one of the three filters (see Puddu et al. (2000)), we produce a matched (in the three bands) catalog complete to r $`19.75`$ and g $`20.2`$ and, in order to ensure completeness, we disregard all objects fainter then r $`=19.5`$. The spatial galaxy distribution in the matched catalog is then binned into equal area square bins of $`1.2`$ $`arcmin`$ to produce a density map. S-Extractor (Bertin & Arnouts (1996)) is run on the resulting image in order to identify and extract all the overdensities having number density $`2`$ $`\sigma `$ above the mean background and covering a minimum detection area of $`4`$ pixels (equivalent to $`4.8`$ $`arcmin^2`$) on the map convolved with a 3 x 3 Gaussian filter. We wish to stress that, as in the Schectman (1985) approach, we are not assuming any a priori cluster model.
The procedure is run on a given plate catalog to produce a preliminary candidate cluster catalog. All previously known Abell and Zwicky clusters are recovered together with many new cluster candidates (Fig. 1 shows, as an example, DPOSS field 610) which need to be confirmed.
The best way to validate cluster finding algorithms would be to apply them to a field for which redshifts are available for a fairly deep magnitude limited sample of galaxies. Since such a sample does not exist, we were obliged to follow an alternative route. It is well known that early-type galaxies are preferentially located in the cores of rich clusters and that they present a small scatter in the color-magnitude diagram (Dressler (1980); Dressler & Gunn (1992); Stanford et al. (1998)). These properties turn the color-magnitude diagram into a powerful tool to disentangle true clusters from overdensities of galaxies caused by random object alignment along the line of sight.
In this paper we use the early type sequences detected in the color-magnitude diagrams obtained from multiband optical photometry to confirm a sample of 12 candidates. These candidate clusters were selected from some fully reduced DPOSS plates available at the time the test was performed. The selection was performed randomly in order to include overdensities covering a wide range of S/N ratios in the detection maps. We also use a sample of 8 X-ray selected clusters at nearby and intermediate redshift as templates to calibrate the photometric redshift estimate.

The paper is structured as follows: in section 2 we describe the data and the data reduction strategy, in section 3 we show the color-magnitude diagrams for the calibration sample and illustrate the procedure implemented to derive the photometric redshift estimate, while in section 4 we apply this method to the candidates. In section 5 conclusions and future developments are discussed.

## 2 Observations and Data Reduction

### 2.1 The observations and the sample

All data used in this paper were obtained in imaging mode with DFOSC at the ESO 1.54m Danish telescope (La Silla - Chile) during two observing runs (March 1999 and March 2000) blessed by dark time and photometric conditions. The CCD (LORAL/LESSER C1 W7) has 2052 x 2052 pixels, each pixel covering $`0.39^{\prime \prime }`$, corresponding to a field of $``$ $`13.3^{}`$ x $`13.3^{}`$. Data were taken in the g, r and i filters of the Thuan & Gunn system (Thuan & Gunn (1976); Wade, Hoessel & Elias (1979)). Seeing averaged $`1.2^{\prime \prime }`$ and was always better than $`1.5^{\prime \prime }`$. The slight difference between our setup and the original Thuan & Gunn filters resulted in a significant color correction which had to be taken into account in the calibration procedure.
Exposure times ranged from 40 to 50 minutes in the g band and 20 to 30 minutes in the r and i bands, depending on the target. In order to obtain higher S/N and more accurate photometric measurements, exposures for those clusters at intermediate redshift were usually repeated in two or three slightly offset frames.
The observed fields (all located in a region with $`0^{}<\delta <+12^{}`$ and therefore observable from both hemispheres) included 12 candidate clusters plus 8 clusters with known redshifts to be used as comparison sample. In order to use the same material to both calibrate the corresponding DPOSS fields and to validate our algorithm, we selected galaxy overdensities in such a way that we had at least two (up to four) candidates and/or clusters in each DPOSS field. One candidate, observed on two different nights, also provided an independent check of the photometric accuracy of the second run. For the comparison sample we selected 8 clusters from the X-ray selected sample of Gioia & Luppino (1994) and Ebeling et al. (1996).
The equatorial coordinates of the observed fields are given in Tables 1 and Tables 2.
We wish to stress that one of the main problems encountered in our work was the well known lack of a suitable set of photometric standards for the Gunn-Thuan system which, along with the lack of faint stars suitable for CCD observations, very often prevents good coverage of the airmass-color plane. The problem is even worse for observers in the Southern hemisphere where the number of available standards is uncomfortably small. We succeeded, however, in observing an average of 4-5 standard stars per observing night.

### 2.2 Data Reduction and Photometric Calibration

The raw images (both scientific and calibration) were prereduced using the standard procedures available in the IRAF package. First, the frames were corrected for instrumental effects (overscan and bias) and flat fielded. Individual dome and sky flats in each filter were median stacked to increase the S/N ratio. For the first run, flatfielding was performed using sky flats only, but the experience gained in this run suggested a slightly different procedure for the second run, using dome flats to achieve better correction of the small scale pixel-to-pixel variations. Dome flats were first used to correct the sky flat frames for the higher frequency fluctuations, and the resulting frames were then smoothed and stacked to map the lower frequency fluctuations and combined with the average dome flat to produce the final master flats.
In the first run, we divided the exposures into two or three frames for the same field; in these cases, the images were combined in each filter by medianing (three exposures) or averaging (two exposures) the aligned frames. Standard star photometry was performed using the apphot package in IRAF. Due to the need to defocus most of the stars to avoid saturation, stars were measured through 10 apertures with diameters up to $`90`$ pixels ($`35.1^{\prime \prime }`$), and the local sky was determined using a $`10`$ pixel wide annulus outside of the largest aperture. In order to determine the zero-point offset and the airmass and color terms we used the $`40`$ pixels ($`15^{\prime \prime }`$) aperture, for the focused and unsaturated stars, and the asymptotic magnitudes for the defocused ones.

The IRAF fitparams task was used to fit the data with the relation

$$m_{true}=m_{inst}+Z_p+K_eX+CIcolor_{inst}$$
(1)
where $`Z_p`$ is the zero point of the magnitude scale; $`K_e`$ is the extinction coefficient and $`X`$ the airmass; $`CI`$ is the instrumental color term coefficient. For the first run (March 1999), due to the paucity of standard stars, we could determine single night coefficients only for the r and i filters. The coefficients were consistent from one night to the other and we used a mean fit for the g-r color.
For the second run (March 2000) we instead derived the coefficients for each night and in each band (using the g-r and r-i colors). The resulting calibration coefficients of the various nights are consistent within the errors and, therefore, in order to improve the quality of the fit, we adopted a unique pair of extinction coefficient and color term for the whole run. These constants were then used to derive the zero points for each night. In Fig. 2 we show for the g and r filters, the fit residuals as a function of the estimated magnitude, using different symbols for different nights.

### 2.3 Object Detection and Photometry

The object catalogs were produced individually for each band using S-Extractor: all objects larger than $`4`$ pixels and $`2\sigma `$ above the background counts were included and their photometric and morphological features measured. We used a photometric reference aperture with a diameter $`3`$ times larger than the average seeing.
For each CCD field, the three single band catalogs where matched taking into account the shifts between pointings (measured using the geomap and geotrans IRAF tasks).

To obtain an estimate of the external photometric errors, the candidate cluster 64\_781 was observed on two different nights. In this way we could evaluate possible night-to-night magnitude offsets in both the g and r filters. The typical weighted mean values for these offsets are $`0.05`$ for the g filter and $`0.007`$ for the r filter, i.e. they are of the same order as the rms errors from the three parameter calibration fit (Fig. 3).

Since our goals require high accuracy for the color determination, we further checked the photometric calibration, using the following test: in the color-color diagram (Fig. 4) we plotted the linear sequence of the Gunn-Thuan standards (open circles) together with all the unsaturated stars (S-Extractor stellarity index $`<0.8`$) within the limiting magnitude, selected from some CCD cluster candidate fields (in Fig. 4 we show the 5\_778 field). For all of these fields, the sequence of the selected stars is linear (excluding the very red sources, dominated by stars of spectral type M, which have a constant g-r color while r-i depends on the spectral subtype; see Fukugita et al. (1996) and Finlator et al. (2000)) and overlap quite well with the standard sequence. This means that the colors of the main sequence stars are well determined, since the relation between g-r and r-i is the same for the cluster field stars and for the standards.

## 3 Validation method via color-magnitude diagrams

From the g and r matched catalogs we excluded obvious stars (stellarity index $`>0.95`$) and then selected a box of $`300`$ pixels in size (corresponding to a typical core cluster diameter of $`750`$ kpc at $`z0.2`$), centered on the approximate cluster center and a second box of equal size located as far as possible from the cluster, to be used for the evaluation of the background contribution.

Our procedure is summarized in Fig. 5, which shows the results for four sample candidates which are representative of the various morphologies encountered.
For all candidates in our sample we first obtained the color-magnitude diagrams for both the cluster and background objects. Then, in order to enhance the early-type sequence, we performed the statistical subtraction of the background contribution by eliminating for each object in the background diagram the corresponding nearest galaxy in the cluster+background diagram. If we then isolate the objects contained within a narrow strip of the color-magnitude diagram centered around the early-type sequence, the galaxy overdensities become more evident in both the spatial distribution and in the number counts radial profile (Fig. 5; the radial profile was calculated by choosing as cluster center the barycenter of the density distribution).
The plots in Fig. 5 can be used as a criterion to distinguish true clusters (candidates (a) and (b)), even if they are difficult to detect. In some cases (usually candidates which are either too distant or too poor), despite the presence of an apparent sequence in the color-magnitude diagram, the objects do not form a physical overdensity, but turn out to be uniformly distributed on the sky. In these cases it is more difficult to reach any definite conclusion about the physical nature of the candidate.

### 3.1 Color-magnitude diagrams for the calibration cluster sample

We chose as templates a sample of X-ray clusters for which, at least in principle, the early-type sequence in the color-magnitude diagrams should be easily detectable. This sample was also used to investigate whether or not it was possible to derive an acceptable estimate of the redshift from the g-r color of the early-type sequences.

In Fig. 6 we plot the color-magnitude diagrams for the 8 clusters in the X-ray sample including only the cluster contribution (i.e., after the statistical subtraction of the background). It is quite evident that some of the early-type sequences are only broadly outlined (MS1401, MS1426, MS1532), which may be caused either by the intrinsic faintness of the cluster members or by cluster structural features (poorness, looseness and presence of interacting systems; see the comments in Gioia & Luppino (1994) about these three clusters). For each cluster we derived a median g-r color using only the 5 brightest galaxies after the background subtraction (continuous line), from which we also estimated the redshift. The crosses represent the g-r colors corresponding to the literature redshifts.

### 3.2 Photometric Redshift Estimate

Some techniques for deriving redshifts from broadband photometry consist of matching observed elliptical galaxy colors with those predicted from the Spectral Energy Distributions (SEDs) (Visvanathan & Sandage (1977)) of a template elliptical galaxy at zero-redshift and corrected according to the redshift: since ellipticals become redder as their redshift increases and since the redshift dependent correction (k-correction), is monotonically increasing in the near and intermediate redshift Universe, colors can be used to infer the cluster redshift.
There is no agreement in the literature for the g-r color of ellipticals at zero-redshift: $`0.47^m`$, according to Schneider, Gunn & Hoessel (1983); $`0.38`$ mag for Frei & Gunn (1994); $`0.40^m`$ to $`0.49^m`$ according to Fukugita et al. (1995). Differences in these values likely depend on the galaxy spectrum template adopted for the ellipticals and on the use of a synthetic or an observed spectrum for the standard stars defining the photometric system. To a lesser extent, differences are due to the variations in the actual shape of the Gunn g and r filters (convolved with the atmosphere, mirror and glass transmissions, CCD quantum efficiency, etc.) and possibly also to the way in which the colors are computed.

In the absence of a definite value, we left the zero-redshift color of ellipticals as the unique free parameter and constrained it with our own observations by (robustly) fitting the relation between color and redshift. Fig. 7 shows (filled dots) the observed colors from the color-magnitude relation vs. known spectroscopic redshift for our X-ray cluster sample. The expected color of ellipticals (continous line) were computed using the Schneider, Gunn & Hoessel (1983) k-correction curve and our own determination of the elliptical colors. The average g-r color of ellipticals of zero redshift turns out to be $`0.44^m`$, i.e. the average of the four previously quoted literature values. Errors on the colors are given as one third of the interquartile range, which roughly corresponds, for a Gaussian distribution of five points to the error on the mean. We prefer these to the standard error since they are more robustly determined. Figure 7 shows that all points are compatible with the curve within $`1\sigma `$, excluding two points , which are within $`2\sigma `$. The agreement is good, provided that there is only one free parameter (the rest-frame elliptical color). Figure 8 compares the photometric redshift, estimated from the color-magnitude diagrams and the spectroscopic redshift. The agreement is good, and the error (interquartile range) on the redshift is, on average, $`\delta z=0.01`$, i.e. $`3000`$ km/s. Tab. 3 lists the estimated photometric redshifts, with the errors computed as previously defined, for the putative clusters; since these clusters are fairly rich systems, this error is likely to be a lower limit.

## 4 The validation of the candidate cluster sample

In Fig. 9 we show the early-type sequences obtained (after subtracting the background) for the 12 candidate clusters in our sample. We confirm 9 of the 12 candidates as true clusters, one is a fortuitous detection of a cluster at $`z0.4`$, one is a false detection and the last is undecidable on the basis of the available data.

The better definition of the early-type sequences observed in the DPOSS confirmed clusters sample with respect to the X-ray sample is likely due to the different specific properties of the two samples: at a given z, our optically selected clusters are on the average richer and more centrally concentrated then the X-ray selected ones.

For the cluster OV27\_694, the early-type sequence is less defined due to the overlap of two independent clusters/groups along the line of sight.
The case of OV21\_694 merits special attention. As a visual inspection of the corresponding POSS-II F plate shows, the marginal detection of an early type sequence (Fig. 5 and 9) seems to be due to the chance alignment of a distant cluster (at redshift $`z0.4`$) with a rich galaxy field. The existence of such a foreground rich field has therefore triggered the search algorithm. As far as OV24\_694 and OV26\_727 are concerned, the marginal evidence for an early-type sequence does not correspond to a defined overdensity in the number count radial profiles (using galaxies in the strip centered on the mean g-r color). The case of OV26\_727 is a false cluster detection, since it has a low S/N ratio and low isophotal richness (see Tab. 1; it may be a group). Visual inspection of the POSS-II plate shows that OV24\_694 lies in a crowded field rich with galaxies; from the sky diagram (Fig. 5, case (d)) it is also evident that a large fraction of these foreground galaxies have the same color. Thus, this field could be part of a larger loose cluster or a cluster in a region with variable background.

## 5 Summary and conclusions

The aim of our work was to test the validity of a model independent cluster finding algorithm, implemented to extract a statistically well defined sample of cluster candidates from photometrically calibrated DPOSS data (see Paper I for details). The advantages of a model independent approach are that i) the program does not assume any a priori knowledge about the clusters, and ii) it objectively looks for statistically meaningful overdensities in the galaxy density field. The main problem in validating any cluster finding algorithm is the lack of a suitable data set to use as a template, i.e., the lack of a region of the sky containing a large sample of clusters with well defined redshifts and properties. In the absence of such a data set, we adopted a photometric approach based on the use of the sequence defined in the color-magnitude diagrams of clusters by bright early-type galaxies as a diagnostic tool. We obtained deep multiband CCD photometry for a sample of 12 candidate clusters extracted from the DPOSS data, plus an additional sample of 8 X-ray clusters with known redshifts to be used as a template to calibrate the photometric redshift procedure. Results may be summarizied as follows: among the 12 clusters candidates, 10 are confirmed clusters, 1 is false and 1 is uncertain. The X-ray selected cluster sample was then used both to check the accuracy ($`\mathrm{\Delta }z0.01`$) and to find the zero point (i.e., the average zero-redshift $`gr`$ color for elliptical galaxies in our system:$`(gr)_0=0.44`$) for the photometric redshift procedure. This procedure is being applied to a larger sample of clusters derived from both DPOSS calibration data and from other archive datasets. Future papers will deal with the analysis of a larger sample of clusters ($`200`$) identified on both DPOSS and archive data and will focus on the derivation and analysis of luminosity functions (both individual and cumulative) and of radial number count profiles (Strazzullo (2001)).