text:
# Looking for Dust and Molecules in Nova V4743 Sagittarii Based on observations collected with the Swedish ESO Submm Telescope at the European Southern Observatory, La Silla, Chile

## 1 Introduction

V4743~Sgr was discovered as a possible nova by Katsumi Haseda (hase02 (2002)) at about 5 mag on September 20th, 2002. West (west02 (2002)) measured its position as $`\mathrm{R}.\mathrm{A}.(\mathrm{J2000})=19^\mathrm{h}01^\mathrm{m}09\stackrel{s}{.}38`$, $`\mathrm{Dec}.(\mathrm{J2000})=22^{}00^{}05\stackrel{}{.}9`$ with an uncertainty of $`0\stackrel{}{.}75`$. Kato et al. (kato+02 (2002)) confirmed that the object is an Fe ii-class nova, and measured the FWHM of the H$`\alpha `$ emission line of 2400 km s-1. The light curves of AAVSO 2002 reveal it to be a very fast nova undergoing a steep decline from its maximum and reaching three magnitudes below maximum within 15 days. This is usually understood as the transition phase where dust formation might begin.

The formation of dust has been suggested in many novae because of a) the behaviour of their visual light curves which indicate the increase of interstellar extinction at the beginning of the transition phase, and b) the development of the infrared excess at the same time (see Bode & Evans bode+89 (1989) for an overview). The dust that forms during nova eruptions consists of a variety of mineralogies, and in some novae, simultaneous evidence is found for probably co-existing carbon–rich and oxygen–rich condensates (Evans et al. evan+97 (1997)).

While most recent novae have been well observed in the optical and even infrared, hardly any observation exists at sub-millimetre or millimetre wavelengths. Still, these measurements are valuable, because they allow to investigate i) the expansion of the ionised shell from emission of free-free radiation as a progenitor of the colder neutral matter, ii) the formation of molecules like CO and other abundant species of the interstellar medium (ISM), iii) the subsequent formation of dust, which has a deep impact on our understanding of the enrichment of the ISM and hence the formation of stars.

The James Clerk Maxwell Telescope nova–monitoring group detected V4743 Sgr with SCUBA (Ivison ivis02 (2002)). On October 3rd, about 12 days after visual maximum, they measured a flux density of $`39\pm 4`$ mJy at 850 $`\mu `$m and determined an upper limit of 225 mJy at 450 $`\mu `$m.

In this paper, we present observations of the 1.2 mm continuum emission and the spectral line emission of CO(1–0) and SiO(3–2) of V4743 Sgr taken in a range between about 20 and 60 days after the visual maximum. We compare the 1.2 mm emission with Ivison’s 450 and 850 $`\mu `$m data and discuss the implications for the origin of the millimetre emission.

## 2 Observations and data reduction

The observations were carried out at the SEST on La Silla, Chile. Millimetre continuum observations at 1.2 mm were taken using the SEST Imaging Bolometer Array (SIMBA) operating at 250 GHz with a FWHM bandwidth of 90 GHz. The maps were produced with the fast-scanning mode (Reichertz et al. reic+01 (2001)) without a wobbling secondary mirror. The correlated sky noise was eliminated during the data reduction process with MOPSI111MOPSI has been developed by Dr. R. Zylka, IRAM, Grenoble, France. The zenith opacity at 1.2 mm was determined with skydips, and ranged between $`\tau _0=0.18`$ and $`0.44`$. Flux calibration was achieved with a series of Uranus maps. We estimate it to be accurate to 10%.

Spectra of the rotational transitions of CO(1–0) at 115.3 GHz and SiO(3–2) at 127.1 GHz were acquired using the SESIS 100/150 heterodyne receiver providing a FWHM beam size of about 60″. We used the acousto-optical spectrometer (AOS) at low resolution (LRS) in order to achieve the largest possible velocity coverage. Two spectra for each molecule were produced, one at the centre velocity of 0 km s-1, and another with an offset of +800 km s-1 (LSR). The separation of each two of the 1440 spectrometer channels is 0.7 MHz being equivalent to 1.8 km s-1 and 1.6 km s-1, respectively. Dual beam switching was chosen for the rejection of the sky emission and system imbalances. A summary of the observational parameters is given in Tables 1 and 2.

## 3 Results

V4743 Sgr was detected in all 1.2 mm continuum measurements with a S/N of at least 4. A confusion with another source is unlikely. Apart from rare phenomena like novae and supernovae, only quasars and similar background sources show a strong variation as witnessed in our study. 30 per $`\mathit{}^{}`$ are expected for the detection limit of our observations (Blain et al. blain98 (1998)). This gives about 0.02 sources for our map coverage.

The individual 1.2 mm maps are presented in Fig. 1. Note that the structures inside the source are small in comparison to the beam size (lower left corner) and are therefore probably not real. The flux densities have been determined in two ways: a) by integrating over a 1′ aperture while subtracting the background ($`S_{\mathrm{ap}}`$), and b) by fitting a two–dimensional Gaussian to the source and integrating the flux within ($`S_\mathrm{G}`$). Table 3 lists both values for all epochs as well as the 1$`\sigma `$ rms noise as the photometric uncertainty. A millimetre light curve is plotted in Fig. 2. We used $`S_\mathrm{G}`$ for a further analysis. The point-like source should have a Gaussian shape, but atmospheric fluctuations distort the images of weak sources. Therefore, a simple aperture photometry might be misleading. In addition, the integration over a Gaussian is more comparable to on-off measurements.

The spectral line observations of the CO(1–0) and the SiO(3–2) rotational transitions only provided upper limits of $`T_\mathrm{A}^{}=0.06`$ K for the CO and $`T_\mathrm{A}^{}=0.03`$ K for the SiO measurements within the covered velocity range.

## 4 Discussion

Millimetre continuum radiation from novae can arise from thermal dust emission (e.g. Evans et al. evan+97 (1997)) or as free–free emission from hot ionised gas (e.g. Ivison et al. ivis+93 (1993)). The decision on which is the dominant emission source is made on the shape of the spectral energy distribution (SED), i.e. the slope $`\alpha `$ of its Rayleigh–Jeans tail. The SED of dust emission is usually described by a modified blackbody law. The flux density $`S_\nu `$ at a frequency $`\nu `$ is given by

$$S_\nu =\mathrm{\Omega }B_\nu (T)\left(1\mathrm{exp}\left(\tau _\nu \right)\right)\tau _\nu =\left(\nu /\nu _\mathrm{c}\right)^\beta $$
(1)
where $`\mathrm{\Omega }`$ is the solid angle of the emitting source, $`T`$ the dust temperature, $`\tau _\nu `$ the optical depth, and $`B_\nu (T)`$ the Planck function. The optical depth is a function of the frequency with an exponent $`\beta =\alpha 2`$, which often is referred to as the submm emissivity index, and $`\nu _\mathrm{c}`$ is the critical frequency, i. e. the frequency where $`\tau _\nu =1`$. Values between $`\beta =1`$ and $`\beta =2.5`$ are found for different dust grains. Free–free emission instead yields a flat SED ($`\alpha =0.1`$) for optically thin and $`\alpha =2`$ for optically thick media (e.g. Chen et al. chen95 (1995)).

Unfortunately, no simultaneous measurements at mm and submm exist for V4743 Sgr. Therefore, we calculated 1.2 mm values for the time of the SCUBA measurements for the cases of free–free and dust emission, and included them in Fig. 2. The steepest submm slope that is consistent with the SCUBA data is $`\alpha =2.6`$. If the mm emission were due to free–free radiation from an expanding ionised shell, it would be difficult to explain the decreasing flux density in the beginning. Hence, free–free emission seems to be an unlikely source, and we favour the dust interpretation. Additional evidence comes from the optical light-curves provided by AAVSO. Coincidently with the increase in the submm, they show a rise of the $`R`$-band intensity, while values at shorter wavelengths remain constant.

The formation of dust in novae is usually observed by an increasing IR emission accompanied by a prominent dip in the visual light curve of the nova which can be explained by extinction introduced by newly formed dust. Such an extinction dip has not been observed for V4743 Sgr. Together with the lack of any emission from CO and SiO molecules, we conclude that the millimetre emission does not originate from newly formed but old dust. It might be either a remnant from the pre-CV phase or a former nova outburst. In that case, V4743 Sgr would actually be a recurrent nova. A similar interpretation was made for V4444 Sgr which showed a rise in the IR without any dip in the optical light curve (Venturini et al. vent+02 (2002)). If we follow the argument of Kawabata et al. (kawa00 (2000)), we can use Eq. (6) of Clayton & Wickramasinghe (clay76 (1976)) to verify that dust cannot be present inside a radius of $``$50 AU from the nova. The largest possible radius of the expanding nova shell of V4743 Sgr can be derived by extrapolating the expansion speed measured from the H$`\alpha `$ line width of 2400 km s-1 for 23 days to the first millimetre maximum. This yields 4.8 109 km or 32 AU. Hence, the dust cannot originate from the recent outburst. The presence of pre-existing dust around novae has been shown for several more cases (e.g. Malakpur mala77 (1977)).

$$M_{\mathrm{dust}}=\frac{S_\nu D^2}{\kappa _\nu B_\nu (T)}$$
(2)
Eq. (2) allows to estimate the mass of the surrounding dust, if $`D`$ the distance to the source, $`\kappa _\nu `$ the mass absorption coefficient at the measured frequency $`\nu `$, and $`T`$ the dust temperature are known. $`\kappa _\nu `$ can differ at least by a factor of 10, depending on the local environment. In our case, we adopt a value of $`\kappa _{250}=2`$ cm2 g-1, which is typical for circumstellar dust (Krügel & Siebenmorgen kruegel94 (1994)) consisting of relatively small grains (Kawabata et al. kawa00 (2000)).

The determination of the distance $`D`$ of CVs is generally quite difficult, although some methods exist to estimate the distance modulus of novae. Firstly, the absolute visual magnitude of a nova 15 days after the maximum is supposed to be similar for all novae, i.e. $`M_V(15)=5.69\pm 0.14`$ mag (Capaccioli et al. capa+90 (1990), and others herein). Secondly, the absolute visual magnitude of a nova at its maximum is related to the time $`t_2`$, in which the nova descends by 2 mag (Cohen cohe88 (1988)). We estimated the distance modulus of V4743 Sgr with both methods and get $`(mM)=13.9\pm 0.4`$ which yields a maximum possible distance of 6 kpc. However, the crucial point is to measure the interstellar extinction towards the nova. The region looks very patchy and HD~176537, the closest star with a measured extinction, still lies 40′ away from the nova itself and hence its value $`A_V=1.42`$ mag (Guarinos guar92 (1992)) might differ quite a lot from that of V4743 Sgr. Additionally, its heliocentric distance is only about 100 pc. Assuming that HD 176537 and the nova are sharing the same interstellar environment, and scaling the extinction value accordingly, we estimate the distance of V4743 Sgr to be around $`D=400`$ pc with $`A_V=5.7`$ mag. Using the $`JH`$ and $`HK`$ colours of the 2Mass point sources in a vicinity of 1′ in comparison with intrinsic colours from Bessel & Brett (bess+88 (1988)), we estimate the average extinction as $`A_V=3.4\pm 0.5`$ mag which yields $`D=1200`$ pc for the nova. We want to point out that these values are very speculative and have to be regarded with caution. For this reason, we relate the mass estimates to a standard distance of 100 pc.

A well established estimate of the dust temperature is crucial for deriving a proper dust mass from millimetre observations. The present data are insufficient to achieve this, and temperature determinations are sparse throughout the literature. Evans et al. (evan+97 (1997)) quote a dust temperature of some 750 K for V705~Cas after 250–300 days after the nova eruption, and Bode & Evans (bode+89 (1989)) give values between 700 and 800 K for NQ~Vul between 80 and 240 days after the outburst. However, these values have been derived from infrared measurements during the ionisation phase and therefore do not trace the coolest dust involved. Consequently, they may not be appropriate for our purposes.

On the other hand, Lynch et al. (lynch+01 (2001)) derive a temperature of 280 K for the cool dust component in V445~Pup. Beck et al. (beck+90 (1990)) present models of nova environments, in which one of the most important ingredients for dust formation, C i, attains a temperature of about 100 K already soon after the nova eruption. We therefore use typical dust temperatures of 100 and 300 K for the mass determination. The so derived dust masses are listed in Table 4 for 100 and 300 K, respectively. As discussed above, they are scaled to a distance of 100 pc.

We determined column densities of the surrounding dust from the derived masses. The projected area in the line of sight was calculated from the lowest possible dust radius, i.e. 50 AU (see above). Since the area can be larger than this, these column densities are only an upper limit and are also listed in Table 4.

If we assume the distance estimates from above, we get average dust masses of about $`\mathrm{5\hspace{0.17em}10}^5M_{}`$ for $`D=400`$ pc and $`\mathrm{5\hspace{0.17em}10}^4M_{}`$ for $`D=1200`$ pc. The referring column densities are $`\mathrm{1\hspace{0.17em}10}^{22}`$ cm-2 and $`\mathrm{1\hspace{0.17em}10}^{23}`$ cm-2, respectively. We recall that a distance of 400 pc is probably a gross underestimate. We therefore think that the actual values are more likely to be found at the higher end.

If the flux variation were due to forming dust, the dust formation rate $`\dot{M}_{\mathrm{dust}}`$ could be estimated from the dust masses between the first two measurements (2.9 days). For the distance estimates from above we find $`\mathrm{3.5\hspace{0.17em}10}^{23}`$ and $`\mathrm{3.1\hspace{0.17em}10}^{24}`$ g s-1, respectively. We compare these dust formation rates with the results from two other novae, Nova Serpentis 1978 (LW~Ser) and Nova Cassiopeiae 1993 (V705~Cas). Gehrz et al. (gehrz+80 (1980)) measured the mass of the dust produced during the outburst of LW Ser from infrared observations. The infrared data show a steep rise to a maximum within 41 days, which we take as the time of the dust formation. The average dust formation rate is therefore $`\mathrm{1.8\hspace{0.17em}10}^{20}`$ g s-1. Evans et al. (evan+96 (1996), evan+97 (1997)) show a visual light curve of V705 Cas with a prominent extinction dip of 29 days, which we define as the dust formation phase. The average dust formation rate in this case yields $`\mathrm{7.9\hspace{0.17em}10}^{18}`$ g s-1. Our rates, if real, would be at least $`2000`$ times higher. A factor of about $`17000`$, however, seems more likely. We interpret this discrepancy again as an indication that the radiating dust was most likely not produced after the nova eruption.

$$\tau _\nu =Nm_\mathrm{d}\kappa _\nu $$
(3)
From Eq. (3) we determine an average optical depth of the medium at 1.2 mm of $`<\tau _{250}>=\mathrm{9\hspace{0.17em}10}^3`$, where $`N`$ is the column density and $`m_\mathrm{d}=\mathrm{4.5\hspace{0.17em}10}^{26}`$ g is the mass of a dust particle. This is taken as 1% of the particle mass of the ISM, where we correct for the typical mass abundance of helium. With Eq. (1) we can calculate the frequency for which the medium becomes optically thick. Transforming this frequency into a wavelength we get $``$175 $`\mu `$m for $`\beta =0.6`$. The size of the bulk of the involved scattering dust particles should be of the order of $`a=\lambda /2\pi `$, i.e. $`38\mu `$m.

## 5 Summary

V4743 Sgr was detected at 1.2 mm, the flux density varying strongly in time. A local intensity maximum has been found around day 23 after the visual maximum. It continued to rise steeply afterwards. The spectral line observations yielded only upper limits for CO(1–0) and SiO(3–2). The dominant emission source is found to be heated dust rather than free–free emission. We have roughly estimated the distance of the nova as $`1200\pm 300`$ pc. Using this value, the dust mass at the first maximum ranges between $`\mathrm{3\hspace{0.17em}10}^4M_{}`$ and $`\mathrm{1\hspace{0.17em}10}^3M_{}`$ depending on its temperature. The high dust mass, the accordingly high formation rate, the absence of molecules, the lack of any dip in the visual light curve while the $`R`$-band magnitude rises simultaneously with the mm flux, and the small shell radius favour the idea that the observed dust has already been present in the system as a remnant of either the pre-CV common envelope phase or a former nova outburst. Dust grains were estimated to have sizes around 38 $`\mu `$m.

###### Acknowledgements.

In this research, we have used, and acknowledge with thanks, data from the AAVSO International Database, based on observations submitted to the AAVSO by variable star observers worldwide. This publication makes furthermore use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We also thank Andreas Lundgren and Dr. Lars-Åke Nyman for carrying out the November observation.


format: text

text:
# Non-deterministic density classification with diffusive probabilistic cellular automata

## I Introduction

Cellular automata (CA) and other spatially-extended discrete dynamical systems are often used as models of complex systems with large number of locally interacting components. One of the primary problems encountered in constructing such models is the inverse problem: the question how to find a local CA rule which would exhibit the desired global behavior.

As a typical representative of the inverse problem, the so-called density classification task Gacs et al. (1987) has been extensively studied in recent years. The CA performing this task should converge to a fixed point of all 1’s if the initial configuration contains more 1’s than 0’s, and to a fixed point of all 0’s if the converse is true. While it has been proved Land and Belew (1995) that the two-state rule performing this task does not exist, solutions of modified tasks are possible if one allows more than one CA rule Fukś (1997), modifies specifications for the final configuration Capcarrère et al. (1996), or assumes different boundary condition Sipper et al. (1998). Approximate solutions have been studied in the context of genetic algorithms in one Mitchell et al. (1994) and two dimensions Morales et al. (2001).

In what follows, we will define a probabilistic CA which solves the density classification problem in the stochastic sense, meaning that the probability that all sites become eventually occupied is equal to the density of occupied sites in the initial string.

We will assume that the dynamics takes place on a one-dimensional lattice with periodic boundary conditions. Let $`s_i(k)`$ denotes the state of the lattice site $`i`$ at time $`k`$, where $`i𝐙`$, $`k𝐍`$. All operations on spatial indices $`i`$ are assumed to be modulo $`L`$, where $`L`$ is the length of the lattice. We will further assume that $`s_i(k)\{0,1\}`$, and we will say that the site $`i`$ is occupied (empty) at time $`k`$ if $`s_i(k)=1`$ ($`s_i(k)=0`$).

The dynamics of the system can be described as follows: empty sites become occupied with a probability proportional to the number of occupied sites in the neighborhood, while occupied sites become empty with a probability proportional to the number of empty sites in the neighborhood, with all lattice sites updated simultaneously at each time step. To be more precise, let us denote by $`P(s_i(k+1))|s_{i1}(k),s_i(k),s_{i+1}(k))`$ the probability that the site $`s_i(k)`$ with nearest neighbors $`s_{i1}(k),s_{i+1}(k)`$ changes its state to $`s_i(k+1)`$ in a single time step. The following set of transition probabilities defines the aforementioned CA rule:

$`P(1|0,0,0)`$ $`=0`$ $`P(1|0,0,1)`$ $`=p`$
$`P(1|0,1,0)`$ $`=12p`$ $`P(1|0,1,1)`$ $`=1p`$
$`P(1|1,0,0)`$ $`=p`$ $`P(1|1,0,1)`$ $`=2p`$
$`P(1|1,1,0)`$ $`=1p`$ $`P(1|1,1,1)`$ $`=1,`$ (1)
where $`p(0,1/2]`$ (the remaining eight transition probabilities can be obtained using $`P(0|a,b,c)=1P(1|a,b,c)`$ for $`a,b,c\{0,1\}`$). The probabilistic CA defined by (I) can be defined explicitly if we introduce a set of iid random variables $`\{X_i\}_{i=0}^L`$ with probability distribution $`P(X_i=1)=p`$, $`P(X_i=0)=1p`$, and another set $`\{Y_i\}_{i=0}^L`$ with probability distribution $`P(Y_i=1)=2p`$, $`P(Y_i=0)=12p`$. Dynamics of the rule (I) can then be described as

$`s_i(k+1)=X_i(1s_{i1})(1s_i)s_{i+1}`$
$`+(1Y_i)(1s_{i1})s_i(1s_{i+1})`$
$`+(1X_i)(1s_{i1})s_is_{i+1}`$
$`+X_is_{i1}(1s_i)(1s_{i+1})`$
$`+Y_is_{i1}(1s_i)s_{i+1}`$
$`+(1X_i)s_{i1}s_i(1s_{i+1})`$
$`+s_{i1}s_is_{i+1}.`$ (2)
To make the above formula easier to read, we omitted the time argument, denoting $`s_i(k)`$ by $`s_i`$. After simplification and reordering of terms, we obtain

$`s_i(k+1)=s_is_iY_i+X_is_{i1}+X_is_{i+1}`$
$`+`$ $`(s_{i1}s_i+s_is_{i+1}2s_{i1}s_is_{i+1}+s_{i1}s_{i+1})(Y_i2X_i).`$
## II Difference and differential equations

The state of the system at the time $`k`$ is determined by the states of all lattice sites and is described by the Boolean random field $`𝐬(k)=\{s_i(k):i=0\mathrm{}L\}`$. The Boolean field $`\{𝐬(k):k=0,1,2,\mathrm{}\}`$ is then a Markov stochastic process. Denoting by $`E_{𝐬(0)}`$ the expectation of this Markov process when the initial configuration is $`𝐬(0)`$ we will now define the expected local density of occupied sites by $`\rho _i(k)=E_{𝐬(0)}\left[s_i(k)\right]`$. The expected global density will be defined as

$$\rho (k)=L^1\underset{i=0}{\overset{L}{}}\rho _i(k).$$
(4)
While both $`\rho _i(k)`$ and $`\rho (k)`$ depend on the initial configuration $`𝐬(0)`$, we will drop this dependence to simplify notation. We will assume that the initial configuration is exactly known (non-stochastic), hence $`\rho (0)=_{i=0}^Ls_i(0)`$ is the fraction of initially occupied sites.

Taking expectation value of both sides of (I), and taking into account that $`E_{𝐬(0)}\left[Y_i2X_i\right]=0`$, we obtain the following difference equation

$$\rho _i(k+1)=\rho _i(k)+p\left(\rho _{i+1}(k)+\rho _{i1}(k)2\rho _i(k)\right).$$
(5)
After summing over all lattice sites this yields

$$\rho (k+1)=\rho (k),$$
(6)
which means that the expected global density should be constant, independently of the value of parameter $`p`$ and independently of the initial configuration $`𝐬(0)`$. We can therefore say that the probabilistic CA defined in (I) is analogous to *conservative CA*, i.e., deterministic CA which preserve the number of occupied sites Hattori and Takesue (1991); Boccara and Fukś (1998); Fukś (2000); Pivato (2001).

Note that up to now we have not made any approximations, i.e., both (5) and (6) are exact. We can, however, consider limiting behaviour of (5) when the physical distance between lattice sites and the size of the time step simultaneously go to zero, using a similar procedure as described in Lawniczak (2000). Let $`x=ϵi`$ and $`t=\tau k`$. Now in (5) we can replace $`\rho (i,k)`$ by $`\rho (x,t)`$, $`\rho (i\pm 1,k)`$ by $`\rho (x\pm ϵ,t)`$ and $`\rho (i,k+1)`$ by $`\rho (x,t+\tau )`$, which results in the following equation:

$$\rho (r,t+\tau )=\rho (x,t)+p\left(\rho (x+ϵ,t)+\rho (xϵ,t)2\rho (x,t)\right).$$

We will consider diffusive scaling in which time scales as a square of the spatial length, meaning that $`\tau =ϵ^2`$. Taking Taylor expansion of the above equation in powers of $`ϵ`$ up to the second order we obtain

$$_t\rho =p_x^2\rho ,$$
(7)
i.e., the standard diffusion equation. Due to the form of (7), in what follows we will refer to the process defined in (I) as *diffusive probabilistic cellular automaton (DPCA).*

## III Absorption probability

We will now present some simulation results illustrating dynamics of DPCA. Since main features of DPCA remain qualitatively the same for all values of the parameter $`p`$ in the interval $`(0,1/2]`$, we have chosen $`p=0.25`$ as a representative value to perform all subsequent simulations.

Let $`N(k)=_{i=1}^Ls_i(k)`$ be the number of occupied sites at time $`k`$. If we start with $`N(0)=0`$, then $`N(k)=0`$ for all $`k>0`$. Similarly, if $`N(0)=L`$, then $`N(k)=L`$ for all $`k>0`$. The DPCA has thus two absorbing states, corresponding to all empty sites (to be referred to as $`\mathrm{𝟎}`$) and to all occupied sites (to be referred to as $`\mathrm{𝟏}`$). If we start with $`0<N(0)<L`$, then the graph of $`N(k)`$ resembles a random walk, as shown in Figure 1.

Both sample trajectories shown there eventually end in the absorbing state, one of them in $`\mathrm{𝟎}`$, another one in $`\mathrm{𝟏}`$. This is a general property of the DPCA: regardless of the initial configuration, the system sooner or later ends up in one of the two absorbing states. Although the time required to reach the absorbing state can be large for a given realization of the process, the expected value of the number of time steps required to reach the absorbing state is finite, as it is the case for all finite absorbing Markov chains Kemeny and Snell (1960). Figure 2 illustrates this property for $`L=100`$ and the initial configuration with $`30`$ occupied sites clustered around the center, ie., located at $`i=35,36,\mathrm{},64`$. All other sites are empty. We start with an assembly of $`200`$ of such initial configurations, all plotted as vertical lines in which black pixels represent occupied sites, while white pixels represent empty sites, as in Figure 2a. Each of these initial configurations evolves according to the DPCA rule, and after $`k=100`$ ($`k=1000`$) iterations they are again plotted as $`200`$ vertical lines in Figure 2b (2c). After $`12000`$ iterations all $`200`$ configurations reach absorbing states, as illustrated in Figure 2d. Obviously, some reach the state $`\mathrm{𝟏}`$, while others $`\mathrm{𝟎}`$, yet it turns out that the fraction of configurations which ended up in the state $`\mathrm{𝟏}`$ is very close to $`30\%`$, the same as the fraction of occupied sites at $`k=0`$.

To explain this phenomenon, let us define $`u_N(k)`$ to be the probability that the number of occupied sites at time $`k`$ is $`N`$. Since the Markov process $`\{𝐬(k):k=0,1,2,\mathrm{}\}`$ is finite and absorbing, no matter where the process starts, the probability that after $`k`$ steps it is in an absorbing state tends to $`1`$ as $`k`$ tends to infinity Kemeny and Snell (1960). This implies that

$`\underset{k\mathrm{}}{lim}u_N(k)=0\text{ if }N0\text{ and }NL,`$ (8)
$`\underset{k\mathrm{}}{lim}\left(u_L(k)+u_0(k)\right)=1.`$ (9)
The expected global density, as defined in (4), is independent of $`k`$, hence

$$\rho (0)=L^1E_{𝐬(0)}\left[N(k)\right]=L^1\underset{N=1}{\overset{L}{}}Nu_N(k).$$
(10)
Taking the limit $`k\mathrm{}`$ of both sides of the above equation, and using (8) and (9), we obtain

$`\underset{k\mathrm{}}{lim}u_L(k)=\rho (0),`$ (11)
$`\underset{k\mathrm{}}{lim}u_0(k)=1\rho (0).`$ (12)
We have shown that *the probability that the DPCA reaches the absorbing state $`\mathrm{𝟏}`$ is equal to the initial fraction of occupied sites $`\rho (0)`$. The probability that it reaches $`\mathrm{𝟎}`$ is $`1\rho (0)`$.* This is in agreement with the behavior observed in Figure 2.

The above can be viewed as a probabilistic generalization of the density classification process. In the standard (deterministic) version of the density classification problem we seek a CA rule which would converge to $`\mathrm{𝟏}`$ ($`\mathrm{𝟎}`$) if the fraction of occupied sites in the initial string is greater (less) than $`1/2`$, i.e.,

$`\underset{k\mathrm{}}{lim}u_L(k)=\mathrm{\Theta }(\rho (0)),`$ (13)
$`\underset{k\mathrm{}}{lim}u_0(k)=\mathrm{\Theta }(1\rho (0)),`$ (14)
where $`\mathrm{\Theta }()`$ is the step function defined as $`\mathrm{\Theta }(x)=0`$ if $`x1/2`$, and $`\mathrm{\Theta }(x)=1`$ if $`x>1/2`$. Thus the difference between the deterministic and the probabilistic density classification introduced here is the replacement of the step function in (13-14) by the identity function in (11-12).

As opposed to deterministically determined outcome in the standard density classification process, in DPCA it is just *more probable* that the system reaches $`\mathrm{𝟏}`$ then $`\mathrm{𝟎}`$ if the fraction of occupied sites in the initial string is greater than $`1/2`$, and it is *more probable* that it reaches $`\mathrm{𝟎}`$ then $`\mathrm{𝟏}`$ if the converse is true. Additionally, DPCA can in some sense *measure* concentration of occupied sites in the initial string. If we want to know what is the initial density of occupied sites, we need to run DPCA many times with the same initial condition until it reaches the absorbing state, and observe how frequently it reaches $`\mathrm{𝟏}`$. This frequency will approximate $`N(0)/L`$, with accuracy increasing with the number of experiments.

## IV Time to absorption

Simulations shown in Figure 1 indicate that for large $`k`$ the system is typically in a state in which blocks of both empty and occupied sites are relatively long. We can use this observation to obtain the approximate dependence of the time required to reach the absorbing state on the density of initial configuration.

If we assume that in a given configuration all occupied sites are grouped in a few long continuous blocks, then the value of $`N(k)`$ cannot change too much in a single time step. For simplicity, let us assume that the only allowed values of $`\mathrm{\Delta }N(k)=N(k+1)N(k)`$ are $`\{2,1,0,1,2\}`$.

Since $`\rho (k)`$ is time-independent, the probability that $`N(k)`$ increases by a given amount should be equal to the probability that it decreases by the same amount in a single time step. In the agreement with the above, let us define $`p_j`$ as the probability that $`\mathrm{\Delta }N(k)`$ takes the value $`j`$, so that $`p_j`$ is non-zero for $`j\{2,1,0,1,2\}`$, where $`p_2=p_2`$, $`p_1=p_1`$, and $`2p_2+2p_1+p_0=1`$. If $`T_z`$ denotes the expected time to reach the state $`N(k)=0`$ or $`N(k)=L`$ starting from $`N(0)=z`$, a simple argument Feller (1968) yields the difference equation which $`T_z`$ must satisfy

$$T_z=p_2T_{z+2}+p_1T_{z+1}+p_0T_z+p_1T_{z1}+p_2T_{z2}+1.$$
(15)
The solution of this equation satisfying boundary conditions $`T_0=0`$ and $`T_L=0`$ is given by

$$T_z=\alpha z(Lz)=\alpha L^2\rho (0)(1\rho (0)),$$
(16)
where $`\alpha =1/(8p_2+4p_1)`$, meaning that the mean time to absorption scales with lattice length as $`O(L^2)`$. The above result would remain valid even if we allowed further jumps than $`\pm 2`$ (although the form of the coefficient $`\alpha `$ would be different).

In order to verify if this result holds for the DPCA, we performed a series of numerical experiments, computing the average time required to reach the absorbing state for 1000 realizations of the DPCA process, for a range of initial densities. Results are shown if Figure 3. One can clearly see that data points are aligned along a curve of parabolic shape, as expected from (16).

## V Conclusion

The probabilistic CA introduced in this article solves the density classification problem in a non-deterministic sense. It is interesting to note that the DPCA conserves the average number of occupied sites, similarly as deterministic rules employed in solutions of related problems mentioned in the introduction. Indeed, conservation of the number of occupied sites is a necessary condition for density classification by CA if one allows modified output configuration, as recently shown in Capcarrère and Sipper (2001). This suggests that a wider class of probabilistic CA conserving $`\rho (k)`$ might be an useful paradigm in studying how locally interacting systems compute global properties, and certainly deserves further attention.

Acknowledgements

The author acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada.