text:
# Broken scaling in the Forest Fire Model

## I Introduction

The Drossel-Schwabl Forest Fire Model Drossel and Schwabl (1992) (DS-FFM) is one of the paradigms of non-conservative SOC Jensen (1998). Its importance comes primarily from the fact that the model has non-conservative microdynamics. It therefore answered the question whether conservation is necessary for criticality in driven systems.

The claim that the DS-FFM is critical comes from the fact that it shows powerlaw-like behavior for several geometrical properties of the dissipation events. The average size of these is divergent in the so-called SOC limit, where all timescales get separated so that the rate of the external drive becomes infinitely slow and the total inflow diverges. If one assumes stationarity this is trivial to prove Drossel and Schwabl (1992); Clar et al. (1994, 1996). However, as usual in numerical simulations, it is not possible to investigate the model in the limit of divergent drive ($`\theta 0`$ in the notation below), as finite size limits the correlation length and therefore destroys any possible criticality Clar et al. (1994). It is remarkable that most of the literature available for this model is mainly concerned with finding critical exponents and identifying supposedly critical quantities. It seems that no authors question whether the model is critical at all and if so in which sense. In this paper we carefully investigate the “scaling function” of the cluster size distribution and show that it is indeed an open question whether the model is truly critical: Not only is there no way to prove its criticality, there is also numerical evidence that the model does not become scale free.

## II Definition of the model and methods

The model has been described several times and in great detail elsewhere Drossel and Schwabl (1992); Clar et al. (1996, 1994). Therefore, the description presented here is rather succinct. The model is defined on a d-dimensional lattice of linear size $`L`$, where each lattice site can be in one of two states, “occupied” or “empty”. The lattice is then updated according to the following rules:

* Driving: Choose randomly $`(1/\theta )`$ sites, one after the other. If its state is “empty” turn it into “occupied”.
* Relaxation: Choose one site at random. If it is empty continue with the first step. Otherwise make “empty” all the sites in the cluster the site chosen belongs to. In this case the update is considered to be successful. Continue with the first step.

Here a cluster is defined in the usual fashion as the set of occupied sites which are connected via nearest neighbor interactions, i.e. two sites belong to the same cluster if they are nearest neighbors or there is a path between them along sites, which belong to the same cluster. We have applied periodic boundaries in all our simulations and restrict ourselves to the two dimensional square lattice.

The cluster removed in the second step is called the “burnt cluster”. To measure the overall distribution of clusters within the system, one usually measures the size of the burnt cluster Clar et al. (1994), the distribution of which is biased by a factor $`s`$. To see that, we define $`n(s)`$ to be the ensemble-averaged, site-normalized density of clusters of size $`s`$ in the system. Then, the probability that a randomly chosen site is connected to a cluster of size $`s`$ is $`sn(s)`$, as in standard percolation Stauffer and Aharony (1994). This distribution is probably the most important in the model. Other quantities are the distribution of the burning time, which is defined as the maximum Manhattan distance (shortest path on the square lattice) from the initially chosen site of a burnt cluster to all other sites in the same cluster, and the correlation functions as defined and discussed in Honecker and Peschel (1997). In this paper we solely concentrate on the distribution $`n(s)`$.

Using a new implementation of the model Pruessner and Jensen  we are able to simulate the system on very large scales and at the same time keep track of the entire distribution $`n(s)`$, instead of measuring the biased distribution $`sn(s)`$, as done usually Clar et al. (1994). Between two updates the changes in $`n(s)`$ are only of the order $`(1/\theta )`$, so it is a highly correlated quantity. However, by using standard cluster labeling techniques Hoshen and Kopelman (1976), it is possible to calculate the full histogram $`n(s)`$ essentially without increasing the computing time, which depends almost exclusively on $`\theta `$ and is essentially independent of the system size. Compared to the standard simulation, we gain up to two orders of magnitude in performance 111 As the standard deviation (actually the estimator of the standard deviation of the estimated mean; this quantity includes the correlation time) decreases with the square root of the computing time, one has to compare the products of the square root of the computing time and the standard deviation. In the algorithm presented here, the ratio of the computing time varies between $`1.4`$ and $`2.1`$ . A similar method was recently introduced for standard percolation Newman and Ziff (2000). Using the same amount of computing time the results are significantly less noisy than those of the standard implementation (for example Honecker (1997), which we have used as reference to check the validity of our results). Large system sizes enable us to rule out any finite size effects as described below. The results have been partly cross-checked using a different random number generator (all results presented here make use of ran2 from Press et al. (1992), and for checking ran1 from Press et al. (1992) has been used).

Finite size effects have been ruled out by the following direct method: For each value of $`(1/\theta )`$ a significantly larger system was simulated with exactly the same value for $`(1/\theta )`$. The linear size $`L`$ was typically increased by a factor $`2`$. The smallest systems we used were $`L=1000`$, the largest $`L=32000`$. By comparing the histograms of different system sizes in conjunction with the standard deviation calculated for them, it is possible to decide whether a system size is affected by finite size effects or not. Compared to other simulations published Schenk et al. (2001, 2000), which also claim not to suffer from finite size effects, our system sizes are huge. This simple method of “redoing” all simulations and using lattices which are much larger than actually needed has the obvious disadvantage of being inefficient, but there is probably no more a direct way of identifying finite size effects Pruessner et al. (2001). This waste of computing power is overcompensated by the efficiency of the algorithm and self-averaging Ferrenberg et al. (1991).

## III Results

The focus of this paper is the presumably universal scaling function of the distribution $`n(s)`$. Similar to the correlation function one expects

$$n(s;\theta )=s^\tau 𝒢(s/s_0(\theta ))$$
(1)
if “simple scaling” applies, which is already known not to be the case in the presence of finite size effects Schenk et al. (2000). The $`L`$-dependence of this quantity is omitted in the following wherever the context allows it. It is worthwhile to note that this is usually the definition of the exponent $`\tau `$. The function $`𝒢`$ is the (presumably) universal scaling function, which depends only on the ratio $`s/s_0(\theta )`$, where $`s_0(\theta )`$ is the only scale of the distribution. This scale depends only on the external parameters, in our case $`\theta `$. Simple scaling allows another scale, namely the lower cutoff, but this is fixed or at least bounded.

The function $`𝒢(x)`$ is usually smooth for small values of $`x`$, therefore it does not make a big difference whether we investigate $`n(s)`$ or

$`\stackrel{~}{n}(s;\theta )`$ $`=`$ $`{\displaystyle \frac{n(s;\theta )}{n(1;\theta )}}`$ (2)
$`=`$ $`{\displaystyle \frac{𝒢(s/s_0(\theta ))}{𝒢(1/s_0(\theta ))}}s^\tau .`$
That this particular choice of the normalization does not affect the overall results can be seen in Tab. 1, where the absolute value of $`n(1,\theta )`$ is listed for different values of $`\theta `$. Also shown in this table is the first moment of the distribution or the average density of trees, which is defined as

$$\rho =\underset{s}{}sn(s;\theta )$$
(3)
and the second moment of the distribution,

$$s=\frac{_ss^2n(s;\theta )}{_ssn(s;\theta )},$$
(4)
which is the average size of the cluster connected to a randomly chosen occupied site.

Before presenting the actual results, we first discuss the numerical quality of the results.

### III.1 Avoiding finite size effects

Throughout this paper we initially performed $`510^6`$ successful updates (as defined in Sec. II) as transient (and therefore rejected them) and the same number for producing statistics, apart from runs for calculating error-bars, where only $`10^6`$ updates has been used for statistics, see below. It is known that the transient can be very long Honecker and Peschel (1997) (note that the time unit in Honecker and Peschel (1997) is expressed in our units by multiplying it with $`(1/\theta )/(\rho L^2)`$), but in all cases presented the number of initial steps seemed to be more than sufficient. Numerical checks indicate that the cluster size distribution is very stable against the size of the transient, i.e. even a transient, which is presumably too short, still produces reasonable results for $`n(s)`$.

All systems have been initialized by a random independent distribution of trees with density $`0.41`$.

The standard deviation of the binned histogram is not completely trivial to calculate. In particular, its computation requires a significant amount of CPU time, and was therefore only calculated for the smaller system sizes (up to $`L=8000`$) and in shorter runs (only $`10^6`$ updates for statistics, but $`510^6`$ for transient). We resorted to visual examination for the larger systems when comparing $`\stackrel{~}{n}(s;\theta ,L)`$ for different system sizes. Fig. 1a and b show the ratio of $`\stackrel{~}{n}(s;\theta ,L)`$ for two different system sizes. A deviation of this ratio from $`1`$ indicates a difference in the statistics and therefore the presence of finite size effects. Fig. 1a shows a typical case we accepted as reasonable agreement. Here $`L_1=4000`$ and $`L_2=8000`$ do not seem to differ for $`(1/\theta )=2000`$. Fig. 1b shows a case of finite size corrections we have dismissed (note the different scales in the two graphs). It differs from Fig. 1a only by $`L_1=1000`$.

Fig. 2 illustrates the strong agreement of $`\stackrel{~}{n}(s;\theta )`$ at the same value of $`\theta `$ for the same two different sizes $`L`$ as in Fig. 1a. The two sets of data are virtually indistinguishable, but in this kind of plot it is also almost impossible to see a difference between the data of $`L_1=1000`$ and $`L_2=8000`$, as shown in the inset of Fig. 2. This is also the case with the rescaled data below, and the use of very large systems throughout this paper might therefore be “overcautious” in avoiding finite size effects, although such large sizes are obviously required for an accurate *quantitative* analysis of this model. However, when it comes only to qualitative analysis, such a judgment seems to be justified. On the other hand, an increase in system size hardly increases the computing time and affects “only” the memory requirements, which forced us to implement the algorithm for parallel machines. The side effect of using multiple CPUs at the same time is a significant reduction of the simulation time especially for large values of $`(1/\theta )`$, a fact which compensates the complications of parallel coding.

Another indicator for the absence of finite size effects is the scaling of the standard deviation of $`\rho `$: If the lattice can be split into independent parts, i.e. if subsets of the lattice can be considered as independent, the standard deviation of $`\rho `$ should scale like $`1/L`$ for different values of $`L`$ at given $`(1/\theta )`$. Such a behavior can be seen in Tab. 1, although the standard deviation of $`\rho `$ could be calculated only roughly. This might explain the slight mismatch for $`(1/\theta )=32000`$, $`L=16000,32000`$.

For the highest values of $`(1/\theta )`$ we could not yet do the comparison to another system, so the curve for the largest value of $`(1/\theta )`$ in Fig. 3 is dotted, as their quality is not known. However, it is reasonable to assume that it is not affected by finite size scaling.

### III.2 The scaling function

Comparing the different histograms $`\stackrel{~}{n}(s;\theta )`$ for different values of $`(1/\theta )`$ in a plot enables us not only to find the exponent $`\tau `$, but also to find the universal function $`𝒢`$ as defined in Eqn. 1. A rough, naive estimate of $`\tau `$ is given by $`\stackrel{~}{n}(s;\theta )`$ fitted against $`s^\tau `$, which gives a value of $`\tau ^{}2.1`$ in our case. Plotting now $`\stackrel{~}{n}(s;\theta )s^\tau ^{}`$ double logarithmically should allow us to find the “true” value of $`\tau `$ by performing a data collapse, i.e. choosing $`\tau ^{}`$ in such a way that horizontal shifts (corresponding to the choice of the scale $`s_0(\theta )`$ in the scaling function) make all curves collapse. This is shown in Fig. 3, where $`\tau ^{}=2.1`$ was chosen so that the maxima for the second bumps are almost equally high: denoting their position on the abscissa for each value of $`\theta `$ by $`s_{\text{max}}(\theta )`$, we have chosen $`\tau ^{}`$ such that

$$\stackrel{~}{n}(s_{\text{max}}(\theta );\theta )s_{\text{max}}^\tau ^{}(\theta )\text{const.}.$$
(5)
According to (1) the constant is simply the maximum value of $`𝒢`$, namely $`𝒢(s_{\text{max}}(\theta )/s_0(\theta ))`$, where the value of the argument is therefore the same for all $`\theta `$.

The value of $`\tau ^{}`$ is close to (but not within the error of) the exponent found in the literature, $`\tau =2.14(3)`$ Clar et al. (1994, 1996) ($`\tau =2.15(2)`$ in Grassberger (1993), $`\tau =2.159(6)`$ in Honecker and Peschel (1997)), which is shown in the same figure for comparison. However, it is impossible to force the minima (see the down pointing marks in Fig. 3) to the same height while maintaining the constraint that the maxima remain aligned, i.e. these minima cannot be a feature of the same universal scaling function. Otherwise (1) would hold and the quantity

$$\stackrel{~}{n}(s_{\text{min}}(\theta );\theta )s_{\text{min}}^\tau ^{}(\theta ),$$
(6)
where $`s_{\text{min}}(\theta )`$ denote the position of the minima, would assume the same value for all $`\theta `$, because they are local minima of $`𝒢`$, which is supposed to be the same for all $`\theta `$.

Since these minima cannot be included in the simple scaling defined in (1), they must be explicitly excluded by introducing a lower cutoff, so that simple scaling supposedly sets in only above these cutoffs, excluding especially the minima. However, such a lower cutoff would apparently have to diverge for $`(1/\theta )\mathrm{}`$ – something that is certainly beyond any established concept of scaling. Even when accepting this peculiar scaling behavior, a data collapse for the second bump still seems to be unsatisfactory, as shown in Fig. 4.

If one accepts a divergent lower cutoff of the scaling function, one has to face the fact that this would describe the behavior of $`\stackrel{~}{n}`$ in a region, which becomes physically less and less interesting in the limit $`(1/\theta )\mathrm{}`$, because the vast majority of events are situated at small $`s`$ and as the second bump moves out to infinity, the scaling function hence covers a smaller and smaller part of $`\stackrel{~}{n}`$. However, only a region of $`\stackrel{~}{n}`$ which covers a non-vanishing fraction of events can be physically relevant.

Concentrating now on the behavior of $`\stackrel{~}{n}`$ up to the minimum (see arrows pointing downwards in Fig. 3), one finds that this region is also badly described by a function like (1). First of all, the question of which region is supposedly described by the function needs to be answered. A unique lower cutoff and a $`\theta `$ dependent upper cutoff needs to be found. At first view it looks appealing to choose these two marks such that they cover the set of data, where the curves fall on top of each other. In this case the lower cutoff would be $`1`$ and the upper cutoff, $`s_{\text{naive}}`$, would have a value smaller than the minima marked by downwards pointing arrows in Fig. 3. However, this would be described by a function like

$$\stackrel{~}{n}(s;\theta )=f(s)𝒢(s/s_{\text{naive}})$$
(7)
rather than (1). Note the parameter independent function $`f(s)`$ describing the shape of the curve, while $`𝒢(s/s_{\text{naive}})`$ is a sharp cutoff function. Eqn. 7 does not allow for an exponent, $`f(s)`$ is an arbitrary function. Writing it as

$$f(s)=s^\tau (a_0+\text{higher order corrections})$$
(8)
defines $`\tau `$ to be the steepest descent of this part of the curve and gives a value between $`\tau _{\text{stp.}}=2.22`$ and $`\tau _{\text{stp.}}=2.19`$ (see Fig. 3).

This concept appears to be rather naive – on the other hand, it is hard to assume that (1) can still hold: it would correspond to (7) with $`f(s)`$ replaced by $`s^\tau `$, which is a straight line in a double logarithmic plot. Therefore (1) can apply only to a region in Fig. 3 where the data that fall on top of each other form a straight line. Those features not already collapsing would then collapse when properly tilted (choosing the right $`\tau `$) and shifted (choosing the right $`s_0`$). Introducing a lower cutoff at $`s=10`$ and discarding the data for $`(1/\theta )2000`$ then leads to a data collapse in a narrow range as shown in Fig. 5. It is worthwhile mentioning that even for some $`10<s<200`$, namely for values of $`s`$ between the squares and the filled circles, none of the data collapse. The exponent used in this “collapse” is $`\tau _{\text{stp.}}=2.19`$, as mentioned above.

By considering the function $`f(s)`$ it becomes apparent that $`\stackrel{~}{n}`$, and therefore the model, cannot be scale free: it depends on the fixed, microscopic scale $`s=1`$. This entails that it is always possible to tell $`(1/\theta )`$ by looking only at the shape of $`\stackrel{~}{n}`$; a diagram showing only this shape, without any scales on the axes, reveals $`(1/\theta )`$, since a scale is intrinsically given by the features of $`f(s)`$. One would only need to rescale and tilt it until it fits the plot Fig. 3 and one could identify $`(1/\theta )`$. Only if $`f(s)`$ were scale free, i.e. a straight line in a double logarithmic plot, would this not be possible.

### III.3 Two length scales

That $`\stackrel{~}{n}`$ contains features to define at least two scales, which apparently diverge in $`(1/\theta )`$ with different exponents, becomes clear when analyzing the scaling of the minima and maxima as marked in Fig. 3, using the definitions

$`s_{\text{min}}(\theta )`$ $``$ $`(1/\theta )^{x_{\text{min}}}`$ (9)
$`s_{\text{max}}(\theta )`$ $``$ $`(1/\theta )^{x_{\text{max}}}.`$ (10)
Of course, the exact position of the extrema of $`\stackrel{~}{n}(s;\theta )s^\tau ^{}`$ depends on its tilt, i.e. on the choice of $`\tau ^{}`$. However, their *scaling* in $`(1/\theta )`$ does not depend strongly on this choice. In particular $`x_{\text{min}}`$ and $`x_{\text{max}}`$ are different for all choices of $`\tau ^{}`$. A plot of $`s_{\text{min}}(\theta )`$ versus $`(1/\theta )`$ for different values of $`\tau ^{}`$ is shown in Fig. 6. For small values of $`(1/\theta )`$ the minimum is not pronounced enough to survive for large values of $`\tau ^{}`$, so these curves do not give a data point. Using a linear fit of $`\mathrm{log}s_{\text{min}}(\theta )`$ versus $`\mathrm{log}(1/\theta )`$ of the minimum as found in the rescaled ($`\tau ^{}`$) and binned histogram, gives an “exponent” between $`x_{\text{min}}=0.93`$ and $`x_{\text{min}}=0.98`$. The same procedure applied to the maxima gives an “exponent” in the range $`x_{\text{max}}=1.18`$ and $`x_{\text{max}}=1.22`$, shown in Fig. 7. One may expect that $`x_{\text{min}}`$ tends towards $`x_{\text{max}}`$ for decreasing $`\tau ^{}`$, as $`s_{\text{min}}`$ increases and might enter the scaling region of $`s_{\text{max}}`$, but neither “exponent” exhibits a systematic variation, and the quality of the fit certainly suffers from the rough procedure that searches for the extrema in the *binned* histogram. This is unfortunately necessary because of statistical fluctuations, in conjunction with the absence of error-bars for all data points.

The scale of the clusters, $`s_{\text{min/max}}`$ is related to the correlation length $`\xi `$ by the fractal dimension $`\mu `$, i.e. (see Clar et al. (1994))

$$s_{\text{min/max}}\xi ^{\mu _{\text{min/max}}}.$$
(11)
Since $`\xi (1/\theta )^\nu `$, one should expect $`\nu =x_{\text{min/max}}/\mu _{\text{min/max}}`$. The minima are supposed to be dominated by smaller, fractal events (see Schenk et al. (2001)), so $`\mu _{\text{min}}=1.96(1)`$ Clar et al. (1994) and therefore $`\nu _{\text{min}}[0.47,0.50]`$. The maxima are more likely to be dominated by compact fires, so $`\nu _{\text{max}}[0.59,0.61]`$. It is unclear how the two exponents $`\nu _{\text{min/max}}`$ are related exactly to the exponents of the two correlation lengths found by Honecker and Peschel Honecker and Peschel (1997) for the connected correlation function $`\nu =0.576(3)`$ and for the tree-tree correlation function $`\nu =0.541(4)`$.

## IV Discussion

Prima facie the DS-FFM looks like a percolation process, and one might naively think that it is indeed a percolation process which organizes itself to the critical density: sites are occupied randomly and independently and (at least in the thermodynamic limit) there is only one cluster which is removed with non-vanishing probability, namely the largest. In this way the density of occupied sites is automatically reduced below the percolation threshold whenever the threshold is reached. It is puzzling how remarkably close the tree density in the DS-FFM is to the density of empty sites in critical site percolation on a square lattice ($`\rho _{\text{FFM}}0.4078`$ and $`1\rho _{\text{perc}}=0.40725379(13)`$ Newman and Ziff (2000) respectively). However, the removal process involved in the DS-FFM introduces spatial correlations which are not present in standard percolation. These correlations are expressed, for example, in the form of a patchy tree density distribution Schenk et al. (2001).

The purpose of this paper is *not* to add yet another model to the enormous zoo of SOC models. However, in order to investigate certain features of the given model and identify underlying mechanisms, it makes sense to modify it slightly. The outcome for the histogram of the DS-FFM modified such that the *largest* cluster is removed after each driving step is shown for a few values of $`(1/\theta )`$ in Fig. 8. The distinctive feature of a minimum which scales differently from the maximum is again present, as the peaks of the maxima have approximately the same height, while the height of the local minima varies among different values of $`\theta `$. The inset of this figure shows the histogram on the same scale as Fig. 3 together with the data of the original model (dotted) with the corresponding values of $`\theta `$. One can understand that they do not fall on top of each other, because the relaxation rule in the modified model erases much larger clusters than in the original model.

Fig. 9 shows a second modification of the model, where again the largest cluster is removed during relaxation and in addition the driving is changed such that the density of trees, $`\rho `$, is the same before each relaxation; the trees removed during the relaxation are just filled in randomly afterwards. This model differs from standard percolation only by its updating scheme 222Actually, it also differs from standard percolation because it fixes the number of occupied sites rather than simply the probability of being occupied. However, this difference becomes irrelevant for sufficiently large systems.. In order to compare the outcome with the original model, the values of $`\rho `$ have been chosen close to the values given in Tab. 1. Indeed, the feature of different scaling of the extrema is still present, but it disappears completely if the density is increased to $`\rho _{\text{perc}}=0.592746`$ Newman and Ziff (2000), which is shown in the same figure as the large bump. This curve does not vary much if a much smaller system size is simulated at this density, so we expect it essentially to be free of finite size corrections. Since it represents a *correlated* percolation process, it is just consistent that this bump does not cover the exact results for the lattice animals of standard percolation Stauffer and Aharony (1994); Sykes and Glen (1976); Mertens (1990) at $`\rho =\rho _{\text{perc}}`$ shown as filled circles in Fig. 9. The dotted graphs in the figure show the corresponding data of the original model. Again they do not match apart from the region of very small $`s`$. Unfortunately the simulations of the so-modified model are very expensive in CPU time, because the mass of the largest cluster needs to be refilled after each relaxation, so that only $`50.000`$ updates for transient and statistics could be done.

Since the feature of different scaling survives the modifications described above, it seems reasonable to assume that any relaxation rule that favors the largest cluster leads to the peculiar behavior. Its disappearance at high densities can be explained by the extremely small cutoff in the distribution, which leads to a domination of the statistics by very small clusters, while a single, enormously large one dominates the burning (the average size of the burnt cluster for $`\rho =\rho _{\text{perc}}`$ was $`355811`$). However, much more careful and detailed investigations of models like the modification described above are required to gain a full understanding of the underlying mechanisms. In particular, this should include a modification of the rules such that the feature disappears.

Honecker and Peschel Honecker and Peschel (1997) have calculated the correlation length not only for the probability that two sites belong to the same cluster, but also for the probability that two sites are occupied at all. The correlation function for the latter is of course a $`\delta `$ peak in ordinary percolation, as there are no spatial correlations for the distribution of occupied sites by construction. However, in the DS-FFM the correlation length for this quantity, $`\xi `$, is finite and seems to diverges when approaching the critical point. It is highly remarkable that Honecker and Peschel conclude from their simulations that this correlation length diverges slightly slower  than the correlation length of the probability for two sites to belong to the same cluster, $`\xi _s`$. This seems to indicate that for sufficiently large scales the spatial correlation of the occupation probability becomes arbitrarily small, so that on sufficiently large scales the DS-FFM occupation is uncorrelated and therefore tends to standard percolation. In other words, it seems to be possible to rescale or “renormalize” the DS-FFM to make the occupation correlation arbitrarily small. This would introduce higher order interactions, as known from standard real space renormalization group and would explain the difference in critical density between the rescaled DS-FFM and standard percolation. However, if this “mapping” is valid, one should find the exponent for the divergence of $`\xi _s/\xi `$ to be the same as in standard percolation, but this is precluded by numerics.

It has been suggested at least twice Honecker and Peschel (1997); Schenk et al. (2001), that the DS-FFM is a superposition of cluster distributions $`n_{\text{perc}(s,p)}`$ of standard percolation for a whole range of concentrations $`p`$, weighted by a certain distribution function $`w(p)`$, i.e. $`_0^1𝑑pw(p)n(s,p)`$. Obviously such an assumption neglects spatial correlations. We recall the following result from standard percolation theory Stauffer and Aharony (1994),

$$n(s,p)s^\tau 𝒞(s/(pp_c)^{1/\sigma }),$$
(12)
where $`𝒞`$ denotes the cutoff function and the exponents $`\sigma `$ and $`\tau _{\text{perc}}`$ have their standard definitions. Assuming that the weighting function $`w(p)`$ is analytic around the critical concentration in standard percolation, $`p_c`$, (12) leads to

$$_0^1𝑑pw(p)n(s,p)s^{(\tau _{\text{perc}}+\sigma )}.$$
(13)
This gives rise to an exponent $`\tau =223/912.45`$, however, this is definitely not supported by numerics (see Fig. 3).

It remains completely unclear how to characterize the scaling of the DS-FFM in two dimensions. Apparently it is not a mere superposition of two simple scalings, as recently speculated Schenk et al. (2001). Moreover the model does not seem to be scale free as described above and it does not seem to be possible to identify a unique power law behavior of the cluster size distribution. Nevertheless *effective* power law behavior over restricted regions has clearly been produced by the model, making it potentially relevant to observation.

All we can conclude is that *the DS-FFM is not critical* in the sense of simple scaling. It reminds us that a divergent moment (here $`s`$, the second moment) can be regarded as a unique sign of emergent scale invariance only if we are certain that one single scale is sufficient to characterize the system. If there is more than one relevant scale, different properties of the system might depend on different scales which may or may not diverge.

###### Acknowledgements.

The large scale simulations represent the backbone of this paper. Andy Thomas at the Department of Mathematics at Imperial College was eager to support us technically, making machines available and upgrading the systems in our department such that parallel computing became a pleasure. A great deal of the results in this paper was possible only because of his work. We thank him very much. This work partly relies on resources provided by the Imperial College Parallel Computing Centre. We want to especially thank K. M. Sephton for his support. Another part of this work was possible only because of the generous donation made by “I-D Media AG, Application Servers & Distributed Applications Architectures, Berlin”. We especially thank M. Kaulke and O. Kilian for their support. G.P. wishes to thank A. Honecker, I. Peschel and K. Schenk for very helpful communication, as well as N. R. Moloney for providing the lattice animal data and for proofreading. The authors gratefully acknowledge the support of EPSRC.


format: text

text:
# Coherent switching of semiconductor resonator solitons

\[

## Abstract

We demonstrate switching on and off of spatial solitons in a semiconductor microresonator by injection of light coherent with the background illumination. Evidence results that the formation of the solitons and their switching does not involve thermal processes.

\] Spatial solitons in semiconductor resonators have recently been considered for applications as mobile binary information carriers for certain optical parallel information processing tasks . We have recently found that such dark and bright spatial solitons can exist in nonlinear passive semiconductor resonators . We demonstrated in  that bright solitons can be switched on (and also temporarily be switched off) by addressing with a light beam incoherent with the background light field, which sustains the solitons. In this case a light-induced increase of temperature  is involved in both switching on and off which limits the switching speed to typically 1 $`\mu `$s.

An alternative mechanism of switching such solitons on or off is by changing the carrier density in the resonator at the location of the soliton. This can be effected by addressing the position of the soliton (to be created or to be switched off) with a sharply focused light beam coherent with the background light. The latter’s intensity is to be within the coexistence range of soliton and low transmission resonator state. If the address beam is in phase with the background light the intensity at the soliton location will be increased when admitting the address beam. This increases the carrier density locally and therefore can switch a soliton on. Vice versa, an existing soliton can be switched off when the address beam is in counterphase to the background light. By destructive interference the light intensity at the location of the soliton can decrease below the existence range of solitons so that a soliton may be switched off.

In this coherent scheme primarily the changing of the local carrier density causes the switching, whereas in the incoherent scheme reported in  the increasing of the temperature causes the switching. Nonthermal incoherent switching can also be done - e.g. by application of local optical pumping  pulses which increase the carrier density. Nonetheless the coherent scheme is the only one permitting on- and off-switching in a straightforward manner . We conducted therefore an experiment to demonstrate this mechanism.

Fig. 1 shows the optical arrangement used for the observations. Light of single frequency Ti:Al203-laser at 850 nm wavelength illuminates an area on the resonator sample of $``$50 $`\mu `$m diameter. Part of the light is split away from the main beam and then superimposed with the main beam by a Mach-Zehnder interferometer arrangement, to serve as the address beam. In the address path of the interferometer the light is expanded so that the focus of the address beam is narrow ($``$ 8 $`\mu `$m). An electro-optic modulator switches the address beam on for a short time ($``$ 40 ns). One of the interferometer mirrors can be moved by a piezo-electric element to control the phase difference between the background light and the address light.

In order to avoid long-term thermal drifts the measurements are carried out within a few $`\mu `$s. The laser light is admitted through an acousto-optic modulator for durations of a few $`\mu `$s and repeated every ms. For observation a CCD-camera combined with an electro-optic modulator serves to take snapshots on a nanosecond timescale. The modulator opens with a variable delay after the start of the illumination so that snapshots can be taken at different times during the illumination. The dynamics is followed more directly locally by imaging a fast photodiode onto the center of one soliton on the resonator sample. The light intensity at the soliton location is then recorded as a function of time.

The semiconductor sample was a passive resonator containing 12 GaAs quantum wells cladded by two Bragg reflector stacks consisting of 23 and 17 pairs of AlGaAs/AlAs quarter-wave layers, respectively. The Bragg mirrors were designed to be non-absorptive above 810 nm. The resonator itself was a 3/2 $`\lambda `$ resonator of the ’cos’ type and the width of the GaAs quantum wells was 15 nm. Suitable spacers between the QW material and the Bragg stacks place the resonance wavelength of this resonator in the vicinity of the band gap at 850 nm. The radial thickness variation of the layers from the center to the edge of the circular 2 inch wafer causes a radial variation of the resonance wavelength of only 7 nm. This very good large-scale homogeneity allows to maintain the relative detuning of laser light to resonance center over large distances on the wafer. This kind of sample should therefore be well suited for experiments on motion of spatial solitons. The resonance wavelength everywhere on the wafer is close to (below) the band edge wavelength. Thus the nonlinearity is strongly absorptive.

Working at $`\lambda `$ = 850 nm, bright spatial solitons are readily observable (Fig. 2); they form spontaneously or with addressing pulses. Figs. 3 and 4 show the results of the switching experiments. The dotted lines show the incident intensity and the solid line the intensity reflected from the sample. Note that the observation is in reflection such that a bright soliton (high transmission) shows up as a reduction in reflected intensity.

Fig. 3a shows the switch-on. The background intensity is chosen slightly below the spontaneous switching threshold. With the application of the in-phase address pulse the reflected light intensity increases (constructive interference) and the soliton is switched on (evidenced by the reduced reflectivity). Evidently the soliton remains switched on until the background light is reduced to the lower limit of the existence range of the soliton.

Fig. 3b shows a ”failed switch” in which the intensity of the address beam is too low to switch. Equally to Fig. 3a the in-phase address pulse shows up as increased reflected light intensity, however, the total intensity at the address position does not reach the upper limit of the soliton existence range so that no soliton forms.

Fig. 3c shows the switching off of a soliton. The background light intensity is chosen above the coexistence range of soliton and low transmission resonator state, so that a soliton forms spontaneously. The address beam is then applied in counter-phase to the background field (evident by a reduction of reflected light due to the destructive interference between background- and address light), the soliton disappears, and the resonator returns to the low transmission state. Fig. 3d shows a ”failed switch-off” analogously to Fig. 3b.

In Fig. 4 we demonstrate that the switching here is really primarily controlled by the light field in the resonator varying the carrier density rather than by temperature increase or generating of phonons as in the incoherent case .

In  a soliton could be switched off by increasing locally the light intensity. The mechanism being that the additional carrier generation increased the lattice temperature, which shifts the band edge. With this the existence range of solitons shifted. The soliton therefore switched off and remained off until the lattice had cooled again. The address pulse was perpendicularly polarized to the background field, thus incoherent with it.

In Fig. 4 the intensity is also increased, but in a coherent way, by constructive interference with the (in-phase) address pulse. As opposed to  the intensity increase does not switch the soliton off, showing that a thermal shift of the existence range does not noticeably occur. Rather here the switching is controlled through the (fast) electronic nonlinearity.

This absence or negligibility of thermal effects in the formation of solitons is equally evidenced by the fast spontaneous switchings in Fig. 3. This nonthermal switching is in agreement with the experiments  where the working wavelength of the structure used was (equal to the structure grown for there present observations) slightly below the band edge.

The absence of the thermal effects in the soliton formation was related in  to the predominantly absorptive character of the nonlinearity. Here essentially the same conditions prevail. From the working wavelength of 850 nm we would conclude that the nonlinearity here is even more purely absorptive than in  as the dispersive nonlinearity should vanish at the precise band edge wavelength. Then the solitons observed here should rely solely on the absorptive nonlinearity and the otherwise important mechanism of nonlinear resonance  would be absent.

Fig. 5 shows fast repeated switching on and off of a soliton. The phase of the address beam was modulated rapidly here by an additional electro-optic phase modulator in the address path at a frequency of 700 kHz and with a phase excursion of somewhat larger than $`\pi `$. Fig. 5 proves again that no slow thermal process is involved in soliton formation or its switch off.

We conclude by noting that these experiments demonstrating the switching of solitons constitute another step towards application of semiconductor solitons to optical information processing e.g. in certain telecom tasks.

Acknowledgments
This work was supported by the European ESPRIT Project ”PIANOS” and by Deutsche Forschungsgemeinschaft under grant We743/12-1. Help from C.O.Weiss and helpful discussions with R.Kuszelewicz are acknowledged.


text:
# Update of Discovery Limits for Extra Neutral Gauge Bosons at Hadron Colliders

## Abstract

We study and compare the discovery potential for heavy neutral gauge bosons ($`Z^{}`$) at the various hadron colliders under discussion at Snowmass 2001 which range in $`\sqrt{s}`$ from 14 TeV for the LHC to 200 TeV for a variant of the VLHC. Typical search limits for $`pp`$ colliders are $`0.250.30\times \sqrt{s}`$ assuming $`100`$ fb-1 to $`1`$ ab-1 of integrated luminosity with some variation due to differences of fermion couplings in the different models. Discovery limits at the Tevatron are $`1`$ TeV for $`15`$ fb-1, approximately 30–50% higher than this rough guideline, due to the higher $`q\overline{q}`$ luminosities in the $`p\overline{p}`$ beams.

Extended gauge symmetries and the associated heavy neutral gauge bosons, $`Z^{}`$, are a feature of many extensions of the standard model such as grand unified theories, Left-Right symmetric models, and superstring theories. If a $`Z^{}`$ were discovered it would have important implications for what lies beyond the standard model. It is therefore important to study and compare the discovery reach for extra gauge bosons at the various facilities that are under consideration for the future cvetic ; godfrey ; capstick ; leike ; riemann ; rizzo . Included in the list of proposed facilities considered at the Snowmass’01 workshop are high energy hadron colliders. In this report we update previous studies cvetic ; godfrey ; capstick ; rizzo  to include the high energy hadron colliders discussed at this meeting which range in $`\sqrt{s}`$ from 14 TeV to 200 TeV.

Many models that predict extra gauge bosons exist in the literature. We present search limits for several of these models that have received recent attention. Although far from exhaustive, the list forms a representative set for the purposes of comparison. The Effective Rank 5 E6 Model starts with the GUT group $`E_6`$ and breaks via the chain $`E_6SO(10)\times U(1)_\psi SU(5)\times U(1)_\chi \times U(1)_\psi `$ with the $`Z^{}`$ charges given by linear combinations of the $`U(1)_\chi `$ and $`U(1)_\psi `$ charges. Specific models are $`Z_\chi `$ corresponding to the extra $`Z^{}`$ of $`SO(10)`$, $`Z_\psi `$ corresponding to the extra $`Z^{}`$ of $`E_6`$, and $`Z_\eta `$ corresponding to the extra $`Z^{}`$ arising in some superstring theories. The Left-Right symmetric model (LRM) is based on the gauge group $`SU(2)_L\times SU(2)_R\times U(1)_{BL}`$, which has right-handed charged currents and restores parity at high energy. The Alternative Left-Right Symmetric model (ALRM) is based on the same gauge group as the LRM but now arising from $`E_6`$ where the fermion assignments are different from those of the LRM due to an ambiguity in how they are embedded in the 27 representation.. The Un-Unified model (UUM) is based on the gauge group $`SU(2)_q\times SU(2)_l\times U(1)_Y`$ where the quarks and leptons each transform under their own $`SU(2)`$ and the KK model (KK) contains the Kaluza-Klein excitations of the SM gauge bosons that are a possible consequence of theories with large extra dimensions. We also consider a $`Z^{}`$ with SM couplings (SSM). Details of these models and references are given in Ref. cvetic .

The signal for a $`Z^{}`$ at a hadron collider consists of Drell-Yan production of lepton pairs capstick ; godfrey ; guts ; pp ; ehlq  with high invariant mass via $`ppZ^{}l^+l^{}`$. The cross section for the production of on-shell $`Z^{}`$s is given by capstick :

$$\frac{\text{d}\sigma (ppf\overline{f})}{\text{d}y}=\frac{x_Ax_B\pi ^2\alpha _{em}^2(g_Z^{}/g_{Z^0})^4}{9M_Z^{}\mathrm{\Gamma }_Z^{}}\left(C_{L}^{f}{}_{}{}^{2}+C_{R}^{f}{}_{}{}^{2}\right)\underset{q}{}\left(C_{L}^{q}{}_{}{}^{2}+C_{R}^{q}{}_{}{}^{2}\right)G_q^+(x_A,x_B,Q^2)$$
(1)
where

$$G_q^+(x_A,x_B,Q^2)=\underset{q}{}\left[f_{q/A}(x_A)f_{\overline{q}/B}(x_B)+f_{\overline{q}/A}(x_A)f_{q/B}(x_B)\right]$$
(2)
The cross section for $`Z^{}`$ production at hadron colliders is inversely proportional to the $`Z^{}`$ width. If exotic decay modes are kinematically allowed, the $`Z^{}`$ width will become larger and more significantly, the branching ratios to conventional fermions smaller. We will only consider the case that no new decay modes are allowed. The partial widths are given (at tree level) by

$$\mathrm{\Gamma }_{Z^{}f\overline{f}}=M_Z^{}g_Z^{}^2(C_{f_L}^2+C_{f_R}^2)/24\pi $$
(3)
We obtain the discovery limits for this process based on 10 events in the $`e^+e^{}+\mu ^+\mu ^{}`$ channels using the EHLQ quark distribution functions ehlq  set 1, taking $`\alpha =1/128.5`$, $`\mathrm{sin}^2\theta _w=0.23`$, and including a 1-loop $`K`$-factor in the $`Z^{}`$ production k-factor . We include 2-loop QCD radiative corrections and 1-loop QED radiative corrections in calculating the $`Z^{}`$ width. Using different quark distribution functions results in a roughly 10% variation in the $`Z^{}`$ cross sections rizzo  with the subsequent change in discovery limits. Detailed detector simulations for the Tevatron and LHC validated our approximations as a good estimator of the true search reach. Furthermore, the results of our previous studies following this approach are totally consistent with subsequent experimental limits obtained at the Tevatron.

Lowering the number of events in the $`e^+e^{}+\mu ^+\mu ^{}`$ channels to 6 raises the discovery reach about $`10\%`$ while lowering the luminosity by a factor of ten reduces the reach by about a factor of 3 capstick .

In our calculations we assumed that the $`Z^{}`$ only decays into the three conventional fermion families. If other decay channels were possible, such as to exotic fermions filling out larger fermion representations or supersymmetric partners, the $`Z^{}`$ width would be larger, lowering the discovery limits. On the other hand, if decays to exotic fermions were kinematically allowed, the discovery of exotic fermions would be an important discovery in itself; the study of the corresponding decay modes would provide additional information on the nature of the extended gauge structure.

The discovery limits for various models at hadron colliders are shown in Fig. 1. These bounds are relatively insensitive to specific models. In addition, since they are based on a distinct signal with little background they are relatively robust limits. Typical search limits for $`pp`$ colliders are $`0.250.30\times \sqrt{s}`$ assuming $`100`$ fb-1 to $`1`$ ab-1 of integrated luminosity with some variation due to differences of fermion couplings in the different models. The Tevatron, a $`p\overline{p}`$ collider, has a 50% higher discovery reach than this rough guideline, indicating that valence quark contributions to the Drell-Yan production process are still important at these energies.

###### Acknowledgements.

This work was supported by the Natural Sciences and Engineering Research Council of Canada. The author thanks JoAnne Hewett for her encouragement.


text:
# Calibration of the distance scale from galactic Cepheids:I

## 1 Introduction: Discussion of the problems related to Cepheids

As an extension of our study of the kinematics of the local universe (KLUN+) we need an accurate value for the global Hubble constant and accurate distances of individual galaxies. The calibration of the distance scale is thus a fundamental step in this process. The aim of this work was to calibrate the distance scale from nearby galactic Cepheids for which the HIPPARCOS satellite measured geometrical parallaxes. This should avoid the step of calibrating the distance scale by assuming a given distance to the Large Magellanic Cloud (LMC). Unfortunatelly, it turns out that these measurements are very difficult to use due to a statistical bias (Lutz and Kelker, 1973). The difficulties can be solved by proper treatment, like the one proposed by Feast and Catchpole (1997). It has been shown that this leads to unbiased results (Pont et al., 1997; Lanoix et al. 1999),

On the other hand, individual measurements of Cepheids from HIPPARCOS are relatively inaccurate because of the distance of galactic Cepheids. Excluding $`\alpha `$ UMi which does not pulsate in the fundamental mode, the best geometrical parallax of an individual Cepheid obtained from HIPPARCOS is 3.32$`\pm `$ 0.58 marcsec for $`\delta `$ Cephee. This leads to an uncertainty in the distance modulus of 0.38 magnitude. In comparison, the quasi-geometrical method of Barnes-Evans applied to Cepheids (Gieren, Fouqué, Gomez, 1998; hereafter GFG), gives distance moduli with a typical uncertainty less than 0.1 magnitude (the external error can be estimated to about 0.2 magnitude according to Table 7 in GFG). We call this method quasi-geometrical because it requires only a few assumptions. The method is independent of any determination of the LMC distance and has a relatively small systematic error (about 0.2 magnitude). Thus, we decided to calibrate the distance scale using the work done by Gieren, Fouqué and Gomez (1998).

Nevertheless, other difficulties appear. The slope of the Period-Luminosity relation (hereafter, PL relation) determined from the adopted calibrating galactic Cepheids differs from the slope obtained for the LMC by the same authors (GFG) (Table 1).. For the LMC, the slopes in V and I bands are now confirmed by the OGLE survey (Udalski. et al., 1999). What slope should we adopt?

The true physical relation is actually a Period-Luminosity-Color (hereafter, PLC) relation written as $`M=\alpha \mathrm{log}P+\beta C_o+\gamma `$, where $`M`$ is the absolute magnitude and $`C_o`$ the intrinsic color. The PL relation is simply the projection of the PLC onto the P-L plane. In the PLC relation the slope $`M/logP`$ is constant. However, the observed slope of the PL relation depends on the distribution of observed Cepheids in the PLC plane (i.e., on the color distribution of the sample). Hence, the slope in a given photometric band may partially depend on the metallicity, because it affects the intrinsic color. Linear non-adiabatic models do predict that the slope is constant when one uses bolometric magnitudes (Baraffe et al., private communication), whereas non-linear models predict that the slope depends on the metallicity also for the bolometric magnitudes (Bono et al., 2000 and references therein) and predict that the slope in a given band depends on the metallicity. Because the metallicity of the LMC differs from the metallicity in the Solar neighbourhood, the choice of slopes in different bands is difficult. In order to avoid this dilemma we decided to apply the method of ’sosie’ (Paturel, 1984) because it does not require knowledge of the slope and zero point of the PL relation 111this method was first introduced to solve the same kind of problems for the Tully-Fisher relation (1977).

The correction for extinction produced by interstellar matter is another difficulty. It can be solved by assuming that the extinction law is universal. We will thus assume that the extinction on an apparent magnitude is proportional to the color excess ($`A_\lambda =R_\lambda (CC_o)`$, where $`C`$ is the reddened color). The factor of proportionality $`R_\lambda `$ is taken from tabulations (e.g., Cardelli, Clayton & Mathis, 1989 ; Caldwell & Coulson, 1987 ; Laney & Stobie, 1993). It depends on both the considered band and color. With such an assumption it is possible to use the Freedman and Madore (1990) precepts of de-reddening. Two bands are needed in order to calculate a color. Because most extragalactic Cepheids are measured in V- and I-band from The Hubble Space Telescope (hereafter, HST), we will use these two bands. Thus, the Freedman and Madore (1990) de-reddening method will be adapted to the sosie method, used in V and I photometric bands.

Finally, an ultimate difficulty comes from the incompleteness bias. This bias was first studied by Teerikorpi (1987) for application to galaxy clusters (Bottinelli et al., 1987). It was first denounced by Sandage (1988) in application to the PL relation and re-discussed later by Lanoix, Paturel and Garnier (1999a). The sample to which we are applying the PL relation must be statistically representative of the calibrators themselves. Indeed, due to the intrinsic scatter of the PL relation, there is a given distribution of absolute magnitudes at a given period. At increasing distances the fainter end of this distribution is progressively missed and the distribution of the actual sample changes. Restricting the sample to Cepheids with a period larger than a given limiting period reduces this bias. The limiting period depends on a first estimate of the distance, on the apparent limiting magnitude and on the characteristics of the PL relation (dispersion, slope and zero-point). In fact, the full theory of Teerikorpi is applicable. The method is much more complete than the rough rule of thumb used as a quick approach in an application in which a detailed treatment was not needed. However, we want to derive final distance moduli and the precise bias correction must be used. Note that the slope and zero point of the PL relation are needed but only as second order terms and thus, the uncertainties mentioned about their choice do not present any significant difficulty (this will be confirmed in section 4.3). The incompleteness bias will be corrected using the precepts given by Teerikorpi (1987).

In section 2 we will describe the material used for this study: the calibrating sample by GFG and our extragalactic Cepheid database (Lanoix et al., 1999b).

In section 3 we describe the ’sosie’ method and give the basic equation for the calculation of the distance modulus of an extragalactic Cepheid.

In section 4 we give the results obtained for 1840 Cepheids belonging to 36 nearby galaxies described in the previous section. We also discuss these results and compare them with those recently published by Freedman et al. (2001).

## 2 Observational material

The guideline in the constitution of the observational material is the selection of the most secure observations. This leads us to reject some data, as explained below, both galactic and extragalactic.

### 2.1 The list of galactic Cepheids

The starting point of our study is the choice of the galactic Cepheids used for the calibration. We adopt the list given in Gieren, Fouqué and Gomez (Table 3 in GFG) but we rejected three Cepheids (EV Sct, SZ Tau and QZ Nor) because they do not pulsate in the fundamental mode (they are overtone Cepheids). They correspond to the three lowest periods of the list. Because we use only the V and I photometric bands, three Cepheids are also rejected (CS Vel, GY Sge and S Vul) because they do not have I-band magnitude. Thus, 28 Cepheids remain. Their distance moduli are adopted directly from Table 5 given by GFG. Only three Cepheids have a mean error in their distance modulus larger than $`0.1`$ magnitude. We give in Table 2 the adopted calibrating sample of galactic Cepheids.

### 2.2 The list of extragalactic Cepheids

In 1999 we have constructed an Extragalactic Cepheid database (Lanoix et al., 1999b) by collecting 3031 photometric measurements of 1061 Cepheids located in 33 galaxies. This list has been updated. Especially, the V and I band measurements by Udalski et al. (OGLE survey, 1999) were added for the LMC from the data available through ’astro/ph9908317’. The new database contains 6685 measurements for 2449 Cepheids in 46 galaxies. In order to make this compilation available, the full contents of the extragalactic part will be published in electronic form for the A&A archives at CDS. A description is given in the Annex.

In this database, each light curve has been inspected in order to describe the main features. In the present study only light curves considered as ’Normal’ are used 222Lanoix et al. give eight classes of light curves: ’Normal’, ’Symetrical’, ’Bumpy’, ’Scattered’, ’Overtone’, ’Low amplitude’, ’Peculiar’, ’No curve’. A ’Normal’ light curve is characterized by a non-symetrical variation: a fast increase and a slower decrease.. We reject all peculiar light curves including light curves classified as ’low amplitude’ because they are often associated with overtone Cepheids.

Only the mean V and I band magnitudes are kept. When several magnitudes are averaged from different sources we keep the mean only if the mean error is less than 0.05 magnitude. It is to be noted that HST measurements of seven galaxies 333IC4182, NGC3368, NGC3627, NGC4496A, NGC4536, NGC4639 and NGC5253 have been analyzed by two independent groups. This leads to two different sets of magnitudes. Independent treatment of both sets shows that the distance modulus differs by less than 0.1 magnitude, except for IC4182 for which the difference is 0.28 magnitude (Lanoix, private communication). Because we have no means to decide which set is the best we decided to keep them both.

The final catalogue (Table 3) results in 1840 extragalactic Cepheids. They belong to 36 galaxies, 27 of which come from HST observations and 9 from ground-based observations. The full Table is available in electronic form in the A&A archives at CDS.

## 3 Method of sosie

The method of ’sosie’ was introduced (Paturel, 1984) to avoid the problem encountered in the practical use of the Tully-Fisher relation (Tully and Fisher, 1977), a linear relation between the absolute magnitude of a galaxy and its 21-cm line width. Here we are in similar conditions with a linear relationship between the absolute magnitude and an observable parameter, the logarithm of the period. In French, the word ’sosie’ refers to someone who looks very similar to someone else without being necessarily genetically related. Here two Cepheids will be considered as ’sosie’ if their light curves have the same shape and if they have the same period (within a given error). Because of the selection based on the shape of the light curve we will consider that all Cepheids of our sample pulsate in the fundamental mode. They all obey the same P-L relation.

We write the distance modulus of a calibrating Cepheid and of an extragalactic Cepheid through a universal PL relation. The calibrating Cepheid is identified with subscripts ’c’ and no subscript for the extragalactic one. Presently, we assume that both stars have the same metallicity and the same intrinsic color. We will see how to bypass this problem, later.

$$\mu _c=m_c^oa\mathrm{log}P_cb,$$
(1)
$$\mu =m^oa\mathrm{log}Pb$$
(2)
$`m^o`$ is the apparent mean magnitude in a given band. The superscript ’$`o`$’ means ’corrected for extinction’. If one selects an extragalactic Cepheid having the same period as the calibrating one, i.e., $`\mathrm{log}P=\mathrm{log}P_c`$, the distance modulus of the extragalactic Cepheid is then

$$\mu =\mu _c+m^om_c^o$$
(3)
The distance modulus of the extragalactic Cepheid is deduced without having to know the slope and zero-point of the PL relation.

In order to correct for extinction we apply the previous equation to two different bands and express the extinction term as a function of the color excess $`E=E(BV)`$. In order to make the notations clearer we note the apparent magnitudes $`V`$ and $`I`$ for the two considered bands. From Eq. 3 one has:

$$\mu =\mu _c+VR_VEV_c+R_VE_c$$
(4)
$$\mu =\mu _c+IR_IEI_c+R_IE_c$$
(5)
which can be written as

$$\mu =\mu _c+VV_cR_V(EE_c)$$
(6)
$$\mu =\mu _c+II_cR_I(EE_c)$$
(7)
Then, eliminating $`EE_c`$ between the two previous equations we obtain:

$$\mu =\mu _c+\frac{(VV_c)\frac{R_V}{R_I}(II_c)}{1\frac{R_V}{R_I}}$$
(8)
This is the desired equation. It can be written in a more elegant manner by using the reddening-free Wesenheit function (Van den Bergh, 1968):

$$W=\frac{V(R_V/R_I).I}{1(R_V/R_I)}$$
(9)
$$\mu =\mu _c+WW_c$$
(10)
In practice, the intrinsic color is not known and this equation is valid only for Cepheids of the same intrinsic color and metallicity. Thus, for a true sample, we will write it as (see the discussion below):

$$\mu =\mu _c+WW_c$$
(11)
$`W`$ is an observable quantity. Then, the mean distance modulus of a sample of Cepheids which have the same period of pulsation as a calibrating Cepheid can be obtained directly from Eq. 11.

The physical relationship in this result is a Period-Luminosity-Color relation. This means that we should search for sosie of calibrators by considering both their similarity in $`\mathrm{log}P`$ and intrinsic color $`C_o`$. But the intrinsic color is not observable. Thus, equation 11 must be considered as a statistical relation exactly as the PL relation. Because of the statistical relation between $`C_o`$ and $`\mathrm{log}P`$, the selection in $`\mathrm{log}P`$ will guarantee that a calibrator Cepheid and its sosies have, on average, the same intrinsic color. So, the problem of the intrinsic color is partially bypassed. For the metallicity problem, the solution is to consider that the method is valid only for galaxies having nearly the same metallicity as the calibrating Cepheids. In the present paper this means that, stricto sensu, only galaxies with a nearly Solar metallicity can be considered as valid. In practice, we applied the method to different kinds of galaxies without noting strong metallicity dependence.

As a test, we apply the method to the calibrating sample itself. Indeed, some galaxies of the sample can be considered as sosie of another. Note that each calibrating Cepheid has at least itself as a sosie. Obviously, we will not consider this special case. We will accept two Cepheids as sosie when the difference of their $`\mathrm{log}P`$ is smaller than $`0.07`$. With a PL slope of $`3`$, this will give an uncertainty $`0.2`$ mag. in the distance modulus. We adopt the ratio $`R_V/R_I=1.69`$ because it corresponds to the most widely accepted one (it corresponds to a ratio of total-to-selective absorption $`A_V/(A_VA_I)=R_V/(R_VR_I)=2.45`$).

In table 4 we give the distance moduli obtained with equation 11 for 23 Cepheids which are sosie of another calibrator. In Figure 1 the comparison of the calculated distance moduli with the calibrating ones is given.

From a direct regression we find that the slope is not different from one ($`1.00\pm 0.03`$). The observed mean difference between the calculated distance modulus and its standard value is obtained together with its standard deviation:

$$\mu _{sosie}\mu _{standard}=0.00\pm 0.24.$$
(12)
The method does not introduce any systematic shift in the zero point. This means that the calibrating Cepheids constitute a coherent system (at least for the 28 Cepheids used in the test). The observed standard deviation (0.24) is in agrement with the expected standard deviation 0.2.

## 4 Application to extragalactic Cepheids

### 4.1 Preliminary determination of extragalactic distance moduli

The method is applied to the 1840 Cepheids of Table 3. To accept two Cepheids as sosie, we still adopt the criterion $`|\mathrm{log}P\mathrm{log}P_c|<0.07`$ which guarantees that the standard deviation is about $`0.2`$ mag., assuming a PL slope of $`3`$. We adopt the ratio $`R_V/R_I=1.69`$ which corresponds to the first order terms proposed by Caldwell & Coulson (1987) and Laney & Stobie (1993). This is also the value adopted by Freedman et al. (2001), following Cardelli et al. (1987), for their HST key project about Cepheids 444This value corresponds to a ratio of total-to-selective absorption $`A_V/(A_VA_I)=R_V/(R_VR_I)=2.45`$.. For each of the 36 host galaxies we plot the different distance moduli given by Eq. 11 as a function of $`\mathrm{log}P`$. This result appears in Figure 3.

The most important feature to point out is a significant trend leading to higher distance moduli for long period Cepheids. This trend is visible for almost all the host galaxies. This is visible even for nearby galaxies if short periods are observed. For distant galaxies the trend is visible also at long periods. This was expected from the incompleteness bias we discussed elsewhere (e.g., Lanoix et al., 1999a). Another signature of the bias comes from the fact that only nearby galaxies (IC1613, IC4182, LMC, NGC224, NGC3109; NGC5253) have Cepheids with short periods. This clearly depends on the limiting magnitude of the considered host galaxy. This important question is discussed in the following subsection.

### 4.2 Correction for the incompleteness bias

In order to get the proper distance moduli we have to correct for the incompleteness bias. In a previous paper (Lanoix et al. 1999a) we suggested using a rule of thumb to avoid this bias. The rule consists of using only $`logP`$ values larger than a given limit $`logP_l`$. This limit is expressed as:

$$logP_l=\frac{V_{lim}\mu b2\sigma }{a}$$
(13)
Unfortunatelly, this method does not take into account the pieces of information contained in smaller periods. The detailed theory of this incompleteness bias was given by Terrikorpi (1987) in the study of galaxy clusters. The bias for extragalactic Cepheids is of the same nature because the Cepheids of a given galaxy are all at the same distance from us, like the galaxies of a cluster. Assuming that the dispersion $`\sigma `$ at a given $`\mathrm{log}P`$ is constant, the basic equations adapted to the problem of extragalactic Cepheids are the following (for the sake of simplicity we will consider only the V band):

The observed distance modulus $`\mu `$ will appear smaller than the true one. The bias $`\mathrm{\Delta }\mu `$ at a given $`\mathrm{log}P`$ is:

$$\mathrm{\Delta }\mu =\sigma \sqrt{\frac{2}{\pi }}\frac{e^{A^2}}{1+erf(A)}$$
(14)
where

$$A=\frac{V_{lim}\mu a_v\mathrm{log}Pb_v}{\sigma \sqrt{2}}$$
(15)
and

$$erf(x)=\frac{2}{\sqrt{\pi }}_0^xe^{t^2}𝑑t.$$
(16)
In these equations $`a_v`$ and $`b_v`$ are the slope and zero point of the PL relation. We adopt the values found by GFG from their galactic sample, i.e., $`a_v=3.037`$ and $`b_v=1.021`$ 555GFG give $`b_v=4.058`$ because they consider the zero-point at $`\mathrm{log}P=1`$. Note that this requirement seems to reduce the interest of the sosie method because the slope and zero-point are needed anyway. In fact, the slope and zero-point appear only as parameters in a second order correction.

Two additional quantities are required to apply these equations:

* The limiting magnitude $`V_{lim}`$.
* The standard deviation $`\sigma `$ of the PL relation at a constant $`logP`$.

The first quantity is derived from the histograms of $`V`$ presented in Figure 2 for each galaxy. We adopt for $`V_{lim}`$ the fainter edge of the most populated class. In a few cases where the histogram has no dominant class, we move the value by $`\pm 0.5`$ magnitude. $`V_{lim}`$ values do not change significantly when one changes the binning size. Only one galaxy (NGC5457) changed by more than the binning size, but its histogram shows two classes with almost the same population. Nevertheless, the global influence of a change in $`V_{lim}`$ is discussed in section 4.3 (Table 5) and we show its influence on each individual galaxy in Table 6. The second quantity ($`\sigma `$) is derived by a direct linear regression on each plot of Figure 3. The adopted quantities $`V_{lim}`$ and $`\sigma `$ are listed in columns 2 and 3 of Table 6.

These parameters being fixed, there is no free parameter to adjust the bias curve to the plot of Figure 3 except the distance modulus $`\mu `$ itself which is then determined through an iterative process. The final bias curves are plotted in Fig.3 for each host galaxy. In column 9 of Table 6 we give the number of remaining sosies after the cut-off at $`V_{lim}`$. In Fig.3 the points which are rejected by the cut-off are represented by crosses.

### 4.3 Analysis of the results

Freedman et al. (2001) recently published their final study of their HST keyproject (HSTKP). They publish distance moduli calculated differently to those used here. They calibrate the PL relation with the LMC distance modulus, assumed to be $`\mu (LMC)=18.5`$. They, adopt the V- and I-band PL relations and an extinction law giving $`R_V/R_I=1.69`$. In order to avoid bias, they cut their sample at a given limiting period $`\mathrm{log}P_l`$ as explained above and they apply a small (but still uncertain) correction for metallicity effect.

The comparisons between the HSTKP results and our solution is shown in Fig.4 for 31 galaxies in common. There is a fair agreement. A direct regression between HSTKP distance moduli and ours leads to a slope which is not significantly different from one ($`1.017\pm 0.010`$) and a zero point difference which is not significantly different from zero ($`0.11\pm 0.16`$). Assuming both determinations carry the same uncertainty, this means that our distances are good within $`0.16/\sqrt{2}=0.1`$ magnitude. This is the internal uncertainty.

From a detailed check of Fig.4 one sees a slight departure from a slope of one at large distances. The effect is, on average, 0.17 mag. for $`\mu `$ larger than 30 mag. Two possibilities can explain this discrepancy:

* The PL relation of the GFG sample shows a departure from linearity for large $`\mathrm{log}P`$. This effect is visible (see for instance Figure 4 in GFG) even when one excludes the three overtone Cepheids ($`\mathrm{log}P<0.6`$). Judging by the error bars of individual points, this non-linearity seems real.
* The distance moduli of Freedman et al. may suffer from a small residual incompleteness bias. Using a simulation we have shown that it is difficult to remove the bias just by cutting the sample at a given $`\mathrm{log}P_l`$. If we refer to our Fig. 7 in Lanoix et al. 1999a, one can see that at large distances ($`\mu >32`$) the bias may reach $`0.17`$ magnitude after the $`\mathrm{log}P`$ cutoff. At intermediate distances ($`29<\mu <32`$) the bias may still reach 0.08 mag.

Three external sources of uncertainty come from: (i) the adopted ratio $`R_V/R_I`$, (ii) the adopted limiting magnitude $`V_{lim}`$ and (iii) from the adopted PL relation used for second order bias correction. In order to check the stability of the solution, we repeated the previous calculations with another PL relation (the one found by GFG for LMC), with a variation of $`R_V/R_I`$ by $`\pm `$ 0.2 and a variation of $`V_{lim}`$ by $`\pm 0.5mag`$. The results are summarized in Table 5, where we give the mean shift between distance moduli from different solutions and the adopted mean distance moduli (reference solution). One can see that the choice of the PL relation has no actual influence on the result. However, a change of $`R_V/R_I`$ by $`\pm `$ 0.1 may change the mean distance modulus by nearly 0.2 magnitude and a change of $`V_{lim}`$ by $`0.5mag`$ may produce similar change. The influence of $`V_{lim}`$ depends clearly on the actual distribution of magnitudes. For some galaxies the effect is negligible while it is large for some others. In order to give a better judgement of the stability of the distance modulus with respect to the adopted $`V_{lim}`$, we give the changes $`\mathrm{\Delta }\mu ^{}`$ when $`V_{lim}`$ is reduced by 0.5 mag (respectively, $`\mathrm{\Delta }\mu ^+`$ when $`V_{lim}`$ is augmented by 0.5 mag).

The actual uncertainty (internal plus external) can thus reach 0.3 magnitude and may be more if our actual sources of uncertainty act in the same sense.

## 5 Conclusion

The distance scale can be calibrated using galactic Cepheids. LMC provides us with numerous Cepheids located at the same distance. This gives a way to derive an accurate slope for the Cepheid PL relation. But its low metallicity (with respect to most of the galaxies of the sample) is a cause of suspicion; we are not sure that this slope can be applied to all kinds of metallicity.

So, we preferred, in a first step, to calibrate the distance scale by using accurate distances of galactic Cepheids published by Gieren, Fouqué and Gomez (1998). These distances are based on the geometrical Barnes-Evans method.

Further, we applied the concept of ’sosie’ (Paturel, 1984) to extend distance determinations to extragalactic Cepheids without having to know either the slope or the zero-point of the PL relation. The distance moduli are obtained in a straightforward way. For the calibrating galactic Cepheids we checked the internal coherence from the same method.

The correction for the extinction is made by using two bands ($`V`$ and $`I`$) according to the principles introduced by Freedman and Madore (1990). There is no need for color excess estimation.

Finally, the incompleteness bias is corrected according to the precepts introduced by Teerikorpi (1987). Without any free parameters (except the distance modulus itself), the bias curve calculated for each individual host galaxy fits very well the observed distance moduli. This gives us confidence in our final distance moduli. Nevertheless, the small departure from the measurements published recently by Freedman et al. (2001) at distances larger than 10Mpc ($`\mu =30`$) must be clarified.

In order to bypass the uncertainty due to metallicity effects it is suggested to consider only galaxies having nearly the same metallicity as the calibrating Cepheids (i.e. Solar metallicity). In Table 6 the distance moduli that can be considered as more secure are noted with an asterisk ($``$). Galaxies with $`\mathrm{\Delta }\mu `$ larger than $`0.3`$ mag. or with small $`n`$ do not receive this flag. For a given ratio $`R_V/R_I`$, the uncertainty of the distances is about 0.1 magnitude but the total uncertainty may be about 0.3 magnitude. The choice of a given $`R_V/R_I`$ ratio is a first source of uncertainty. The actual ratio depends on the extinction law in our Galaxy, on the extinction law in the host galaxy and on the color of the considered Cepheid. For the future it would be interesting to search for a clue allowing us to decide which value is the best in a given direction for a Cepheid in a given host galaxy. The proper determination of the limiting magnitude of the sample is a second source of uncertainty. It can be accurately determined only when a large number of Cepheids is available to provide us with good statistics.

Presently, the calibration of the distance scale can barely be better than $`\sigma _\mu =0.3`$ magnitude. Thus, the uncertainty on the Hubble constant, $`\sigma (H)\sigma _\mu H/5\mathrm{log}e`$, cannot be better than about $`10km.s^1.Mpc^1`$.

## Appendix A The extragalactic Cepheid database

The description of this database was given by Lanoix et al. (1999b). Because the database is no longer available on the world-wide-web the present data are published in electronic form in the A&A archives at CDS. All the data are made available, even when they are not used in the present paper, where only Normal Cepheids in V and I-bands are considered. Additional measurements were collected including the LMC ones by Udalski et al. (1999)666Note that we kept 720 normal Cepheids among the 1182 available with $`\mathrm{log}P>0.5`$. and those by Gibson et al. (1998, 1999). Data are now available for 2449 Cepheids of 46 galaxies (instead of 1061 Cepheids of 33 galaxies).

The identification of a Cepheid is given on a first line as follows:

* the name of the host galaxy,
* the name of the Cepheid and the bibliographic code from where this name is taken,
* the adopted period (in log)
* The classification of the shape of the light curve, following Lanoix et al. (1999b)

On this first line we also give the number of measurements attached to this Cepheid. Note that the Cepheid name for LMC is simply the Cepheid number from Udalski et al., without the field number (SC), that was not needed here (only three Cepheids appear with the same number in different fields: 1, 16 and 19, but they are not in our list). We tried to keep the Cepheid name of the first discovery. This was not always done, e.g., the names given by Graham (1984) are referenced as Mad87 because of the renumbering adopted by Madore (1987).

On the following lines, individual measurements are given:

* the magnitude,
* the type of magnitude (mean, maximum, minimum, average) coded according to Lanoix et al.,
* the photometric bands (B,V,R,I …) coded according to Lanoix et al.(1999b),
* the reference code. The full reference and the associated code appears in the references.

A sample is given below to show how the data are organized.

```

IC1613  V1         Fr88a  0.7480     N    8
 21.36     mea       B   Fr88a
 20.79     mea       V   Fr88a
 20.36     mea       R   Fr88a
 20.14     mea       I   Fr88a
 20.50     max       B   Sa88a
 22.03     min       B   Sa88a
 21.27     ave       B   Sa88a
 21.39     mea       B   Sa88a
IC1613  V20        Fr88a  1.6220     B    5
 16.66               H   Ala84
 18.98     max       B   Sa88a
 20.71     min       B   Sa88a
 19.85     ave       B   Sa88a
 19.90     mea       B   Sa88a
IC1613  V22        Fr88a  2.1650     S    9
 15.47               H   Ala84
 19.10     mea       B   Fr88a
 17.75     mea       V   Fr88a
 17.14     mea       R   Fr88a
 16.62     mea       I   Fr88a
 17.74     max       B   Sa88a
 20.44     min       B   Sa88a
 19.09     ave       B   Sa88a
 19.09     mea       B   Sa88a
IC1613  V25        Fr88a  0.9600     B+   5
 18.62               H   Ala84
 20.10     max       B   Sa88a
 21.84     min       B   Sa88a
 20.97     ave       B   Sa88a
 20.87     mea       B   Sa88a
IC1613  V53        Fr88a  0.5900     O    3
 21.13     max       B   Car90
 21.70     min       B   Car90
 21.46     mea       B   Car90
...
...
...

```

###### Acknowledgements.

We thank the HST teams for making their data available in the literature prior to the end of the project. We thank R. Garnier, J. Rousseau and P. Lanoix for having participated to the maintenance of our Cepheid database. We thank P. Teerikorpi for his comments and the anonymous referee for very constructive remarks.