text:
# (Strange) Meson Interferometry at RHIC

## 1 Introduction

We discuss predictions for the kaon interferometry measurements in Au+Au collisions at the Relativistic Heavy Ion Collider (RHIC) that accelerates the nuclei up to nucleon-nucleon center-of-mass energies of $`\sqrt{s}_{NN}=200`$GeV. Correlations of identical particle pairs, sometimes also called HBT interferometry, provide important information on the space-time extension of the particle emitting source as produced for example in ultrarelativistic heavy ion collisions . In this case, QCD lattice calculations have predicted a transition from quark-gluon matter to hadronic matter at high temperatures. For a first-order phase transition, large hadronization times have been expected due to the associated large latent heat as compared to a purely hadronic scenario . Entropy has to be conserved while the number of degrees of freedom is reduced throughout the phase transition. Thus, one has expected a considerable jump in the magnitude of the interferometry radius parameters and the emission duration once the energy density is large enough to produce quark-gluon matter. Two alternative scenarios of the space-time evolution, with and without a phase transition, are illustrated in Fig. 1 in the $`z`$-$`t`$-diagram. After the collision of the two nuclei, each with nucleon number $`A`$, the system is formed at some eigen-time $`\tau `$ (indicated by the hyperbola) and the initial expansion proceeds either in a hadronic state (left-hand side) or in a state dominated by partonic degrees of freedom, for example a quark-gluon plasma (QGP) (right-hand side). In the latter case, the formation of a mixed phase, leads to large hadronization times and thus to rather long emission durations. The freeze-out is defined as the decoupling of the particles, i.e., the space-time coordinates of their last (strong) interactions. Summarizing this idea illustrated in Fig. 1, the interferometry of identical particle pairs and in particular the excitation function of the interferomtry parameters have been considered as an ideal tool to detect the existence and the properties of a transition from a thermalized quark-gluon plasma to hadrons. Quantities of great interest are for example the critical temperature $`T_c`$ or the latent heat of the phase transition. The prolonged emission duration should lead in particular to an increase of the effective source size in the outward direction, i.e., parallel to the transverse pair velocity. One also expects correlation lengths depending on the specific entropy of the collision. For recent reviews see for example .

Figure 1: Illustration of the space-time evolution in the $`z`$-$`t`$-diagram with (right) and without (left) a phase transition. Proceeding through a first-order phase transition with a large latent heat should lead to large hadronization times, thus yielding eventually large interferometry radius parameters and emission duration.

## 2 Space-time evolution and kaon interferometry

Here, we discuss calculations based on a relativistic two-phase dynamical transport model that describes the early quark-gluon plasma phase by hydrodynamics and the later stages, after hadronization from the phase boundary of the mixed phase, by microscopic transport of the hadrons. For the initial dense (hydrodynamical) phase of a QGP a bag model equation of state exhibiting a first order phase transition is employed . Hence, a phase transition in local equilibrium that proceeds through the formation of a mixed phase, is considered. A crossover  or a rapid out-of-equilibrium phase transition similar to spinodal decomposition  may yield smaller radii and emission times. In the hadronic phase, resonance (de)excitations and binary collisions are modeled based on cross sections and resonance properties as measured in vacuum . A detailed description of this relativistic hybrid transport model and its predictions can be found elsewhere . The particular differences and advantages compared to pure hydrodynamical calculations are for example due to the explicit calculation of the freeze-out phase-space distributions. Moreover, the system evolution in the later dilute stages is based on cross sections and resonance (de)excitations are possible. The assumption of ideal hydrodynamics might be not valid anymore in these later stages close to a freeze-out criterion like a fixed temperature $`T_{\mathrm{fr}}`$. Another important example will be shown below (Fig. 5), demonstrating the possibility of direct emission from the early phase in this transport model used here. The mean freeze-out times are strongly depending on the transverse momentum of the particles as illustrated in Fig. 2 for kaons and anti-kaons. Obviously, large transverse momentum kaons are correlated to early mean freeze-out times. The difference between kaons and anti-kaons is almost negligible, only at very small $`p_t`$ anti-kaons may be emitted a tiny instant later due a small finite net-baryon number in the system. (The K-N cross section is considerably larger than the K+N cross section.) Otherwise the K$`\pi `$ cross section, dominated by the K, yields almost identical freeze-out distributions for kaons and anti-kaons. Strangeness distillation  due to local large net-baryon numbers leading to emission asymetries that are observable in K+K- correlations  is not considered here.

Figure 2: Mean freeze-out time $`\tau _F`$ as a function of transverse momentum $`p_t`$ for kaons and anti-kaons in Au+Au collisions at RHIC. High $`p_t`$ kaons are correlated with early freeze-out times.

The coordinate system for the correlation analysis is defined by the long axis ($`z`$) (parallel to the beam axis), the out direction (parallel to the transverse momentum vector $`𝐊_𝐓=(𝐩_{\mathrm{𝐓𝟏}}+𝐩_{\mathrm{𝐓𝟐}})/2`$ of the pair), and the side direction (perpendicular to both). Due to the definition of the out and side direction, $`R_{\mathrm{out}}`$ probes the spatial and temporal extension of the source while $`R_{\mathrm{side}}`$ only probes the spatial extension. It has been suggested that the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ should increase strongly once the initial entropy density $`s_i`$ becomes substantially larger than that of the hadronic gas at $`T_c`$ . The Gaussian radius parameters can be obtained from a saddle-point integration over the classical phase space distribution of the hadrons at freeze-out (points of their last (strong) interaction) that is identified with the Wigner density of the source, $`S(x,K)`$ .

$`R_{\mathrm{side}}^2(𝐊_𝐓)`$ $`=`$ $`\stackrel{~}{y}^2(𝐊_𝐓),`$ (1)
$`R_{\mathrm{out}}^2(𝐊_𝐓)`$ $`=`$ $`(\stackrel{~}{x}\beta _t\stackrel{~}{t})^2(𝐊_𝐓)=\stackrel{~}{x}^2+\beta _t^2\stackrel{~}{t}^22\beta _t\stackrel{~}{x}\stackrel{~}{t}(𝐊_𝐓),`$ (2)
$`R_{\mathrm{long}}^2(𝐊_𝐓)`$ $`=`$ $`(\stackrel{~}{z}\beta _l\stackrel{~}{t})^2(𝐊_𝐓),`$ (3)
with

$$\stackrel{~}{x}^\mu (𝐊_𝐓)=x^\mu x^\mu (𝐊_𝐓)$$
(4)
being the space-time coordinates relative to the momentum dependent effective source centers. The average in (1)-(3) is taken over the emission function, i.e.,

$$f(K)=\frac{d^4xf(x)S(x,K)}{d^4xS(x,K)}$$
(5)
with $`K=(E_K,𝐊)`$. In the osl system $`\beta =(\beta _t,0,\beta _l)`$, where $`\beta =𝐊/E_K`$ and $`E_K=\sqrt{m^2+𝐊^2}`$ (on-shell approximation). For small $`\stackrel{~}{x}`$-$`\stackrel{~}{t}`$ correlations, i.e. in particular at small $`K_T`$, $`R_{\mathrm{out}}`$ is increased relative to $`R_{\mathrm{side}}`$ if the duration of emission $`\mathrm{\Delta }\tau =\sqrt{\stackrel{~}{t}^2}`$ is large . Strong (positive) $`\stackrel{~}{x}`$-$`\stackrel{~}{t}`$-correlations or large spatial anisotropies in out\- and side-direction ($`\stackrel{~}{y}^2>\stackrel{~}{x}^2`$) may, in principle, lead to $`R_{\mathrm{out}}R_{\mathrm{side}}`$. This is not seen within the model scenario discussed here.

From eqs. (1) and (2) one can easily calculate the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ for kaons. Kaons, compared to pions, are expected to be less contaminated by resonance decays . In addition the kaon phase space density is much smaller than the pion density . As a consequence, higher multiparticle correlation effects (that may be important for pions at RHIC energies ) are possibly under much better control for kaons. Fig. 3 shows the kaon $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio for SPS (initial conditions $`\tau _i=1`$fm/c (thermalization eigen-time) and $`s/\rho _B=45`$ (specific entropy density)) and RHIC ($`\tau _i=0.6`$, $`s/\rho _B=200`$. The bag parameter $`B`$ is varied from $`380`$MeV/fm3 to $`720`$MeV/fm3. This corresponds to a variation of the critical temperature from $`T_c160`$MeV to $`T_c200`$MeV.

Figure 3: $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ as obtained from eqs. (1) and (2) for kaons at RHIC (Au+Au, $`\sqrt{s}_{NN}=200`$GeV) (full symbols) and at SPS (Pb+Pb, $`\sqrt{s}_{NN}=17.4`$GeV) (open symbols), as a function of $`K_T`$ for critical temperatures $`T_c160`$MeV and $`T_c200`$MeV, respectively.

For small $`K_T`$, the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ increases slowly compared to the faster rise for pions . This is due to the different masses yielding different flow velocities at the same $`K_T`$ (see eq. (2)). Apparently, the sensitivity to the critical temperature $`T_c`$ and the specific entropy density increases strongly with $`K_T`$. For higher $`T_c`$ the hadronization is faster but the subsequent hadronic rescattering phase lasts longer. It was shown that this dissipative hadronic phase dominates the radii for pions . The ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ is always larger than unity and reaches values on the order of $`1.52`$ at large $`K_T`$. This value is considerably smaller than former expectations (e.g. ). On the other hand, first RHIC data for pion correlations show ratios that do not increase with $`K_T`$ and are even smaller than unity . This completely new behaviour has not been seen at SPS energies (see e.g. ). However, new data by the CERES collaboration show a similar trend . This observation would hint at a rather explosive scenario with very short emission times, not compatible with a picture of a thermalized quark-gluon plasma hadronizing via a first-order phase transition to an interacting hadron gas . Rather a shell-like emission as illustrated in Fig. 4 would be prefered. Thus, the further study of two-particle interferometry will provide extremly important information e.g. on the hadronization process or the question of thermalization in ultrarelativistic heavy ion collisions. Also, the comparison of pion and kaon interferometry data may further clarify the situation.

Figure 4: Illustration of a shell-like emission. The surface emission geometry corresponds to small values of the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$, indicated by the two dashed lines in the emission volume element relevant for an observer. The dashed line in the direction to the observer corresponds to the out direction, the orthogonal line is the side direction.

Kaons are less distorted by long-lived resonances and escape the opaque hadronic phase easier (compared to pions). About $`30\%`$ of the kaons at $`K_T1`$GeV/c are directly emitted from the phase-boundary (for $`T_c160`$MeV) as can be seen in Fig. 5. Complementary, we have seen already that large $`K_T`$ kaons and their $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio exhibit a strong sensitivity on the QCD equation of state. The fraction of resonance decays (K’s,$`\mathrm{}`$) is still quite large and decreases with $`K_T`$ from about 70$`\%`$ to 50$`\%`$. For the higher $`T_c`$ ($`200`$MeV), hadronization is earlier and the hadronic phase lasts longer, such that the system is more opaque for direct emission than in the lower $`T_c`$ case.

Figure 5: Fraction of kaons $`\mathrm{\Gamma }_{\mathrm{origin}}`$ that origin from a particular reaction channel prior to freeze-out. These are resonance decays (full circles), direct emission from the phase boundary (stars), elastic meson-meson (diamonds), or elastic meson-baryon (open circles) collisions. The upper and lower diagrams are for RHIC and SPS initial conditions for $`T_c160`$MeV (left) and $`T_c200`$MeV (right), respectively.

We now calculate the one-dimensional kaon kaon correlation functions $`C_2(Q_{\mathrm{inv}})`$ for various $`K_T`$-bins and determine the corresponding fit parameters. The corresponding results for the three-dimensional analysis in the out-side-long coordinate system are presented in . The parameters $`R_{\mathrm{inv}}`$ and $`\lambda _1`$ of the correlation functions are obtained by fitting a Gaussian as

$$C_2(Q_{\mathrm{inv}})=1+\lambda _1\mathrm{exp}(Q_{\mathrm{inv}}^2R_{\mathrm{inv}}^2).$$
(6)
The correlation functions themselves are calculated from the phase space distributions of kaons at freeze-out using the correlation after burner by Pratt . The correlation functions can be calculated as

$$C_2(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})1+\frac{\mathrm{d}^4xS(x,𝐊)\mathrm{d}^4yS(y,𝐊)\mathrm{exp}\left(2ik(xy)\right)}{|\mathrm{d}^4xS(x,𝐊)|^2}.$$
(7)
with $`2𝐊=𝐩_\mathrm{𝟏}+𝐩_\mathrm{𝟐}`$, $`2𝐤=𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}`$, and $`2k^0=E_{p_1}E_{p_2}`$.

Figure 6: Parameters $`R_{\mathrm{inv}}`$ and $`\lambda _1`$ (as obtained from a fit according to eq. (6)) as a function of $`K_T`$ for kaons in Au+Au collisions at RHIC for critical temperatures $`T_c160`$MeV (circles, squares) and $`T_c200`$MeV (triangles, diamonds), respectively. Calculations with and without taking momentum resolution effects into account are shown.

The results for RHIC are shown in Fig. 6. The $`\lambda _1`$ intercept parameters are multiplied by a factor 10. The correlation functions have been calculated for $`T_c160`$MeV (circles and squares) and $`T_c200`$MeV (triangles and diamonds). Moreover, calculations with and without taking finite momentum resolution (f.m.r.) effects have been performed. In the first case, the true particle momentum $`p`$ obtains an additional random component. This random component is assumed to be Gaussian with a width $`\delta p`$. A f.m.r. of $`2\%`$ of the center of each $`K_T`$ bin is considered, a value assumed for the STAR detector . The relative momenta of pairs are then calculated from these modified momenta. However, the correlator is calculated with the true relative momentum. $`R_{\mathrm{inv}}`$ is on the order of $`5.56.5`$fm without and $`4.56`$fm with taking f.m.r. into account. Hence, $`R_{\mathrm{inv}}`$ is reduced by the f.m.r.. This reduction grows with $`K_T`$. For the one-dimensional fit parameters, there is hardly a sensitivity to the critical temperature. However, the known trend (larger radii for higher $`T_c`$ due to a prolonged hadronic phase) is also indicated in these results. The sensitivity to the specific entropy density, another initial condition, is known to be rather weak , even for a change from SPS to RHIC initial conditions, due to the dominance of the late soft hadronic interactions. Hence, the differences in the correlation parameters between, for example, the RHIC 130 GeV and 200 GeV energies are minute; (the number of charged particles at midrapidity increases only by about $`15\%`$). Preliminary analysis of STAR KK data (not yet corrected for f.m.r.) at $`\sqrt{s}_{NN}=130`$GeV gives a value of $`R_{\mathrm{inv}}=4.5\pm 0.3`$fm (at $`150`$MeV$`<K_T<400`$MeV and rather central (11$`\%`$) collisions) ($`\lambda _1=0.92\pm 0.13`$, which is almost compatible to the calculations that take f.m.r. into account. The $`\lambda _1`$ intercept parameters are almost constant without f.m.r. effects but strongly decrease when taking them into account. For illustration, let us consider particle pairs with a small true relative momentum $`q`$ that will be redistributed by the random momentum smearing to, on the average, larger relative momenta $`\stackrel{~}{q}`$. Thus, the area around $`q0`$ is depleted from these true low momentum pairs that carry the full correlation strength. Complementary, some pairs with a larger true $`q`$ (and thus with a weaker correlation strength) are redistributed to the small $`\stackrel{~}{q}`$ area. As a consequence, the correlation strength is transported to larger $`q`$ values. Thus, the correlation function gets broader and the radius parameters become smaller. The strong decrease of $`\lambda _1`$ with $`K_T`$ is due to the absolute larger momentum smearing at high $`K_T`$. True correlated pairs at low $`q`$ are transported more likely to larger $`\stackrel{~}{q}`$ values. This reduced correlation strength in the first $`q`$ bins causes the small $`\lambda _1`$ values at larger $`K_T`$. Recent experimental data (corrected for f.m.r.!) for central Pb+Pb collisions at $`\sqrt{s}=17.4`$GeV yield $`R_{\mathrm{inv}}=6.2\pm 0.4\pm 0.9`$fm at low transverse momenta ($`p_t=250`$MeV) and $`R_{\mathrm{inv}}=3.3\pm 0.6`$fm at high transverse momenta ($`p_t=910`$MeV) . The corresponding $`\lambda _1`$ parameters are reported as $`\lambda _1=0.84\pm 0.11\pm 0.19`$ and $`\lambda _1=0.47\pm 0.12`$ for the low and high $`p_t`$ set, respectively.

## 3 Summary

The calculation of kaon correlation parameters for Au+Au collisions at RHIC energies, assuming a first-order phase transition from a thermalized QGP to a gas of hadrons shows

* an increasing $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio with $`K_T`$ (larger than unity),
* an increasing sensitivity of the $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio with $`K_T`$,
* a strong ($`30\%`$) direct emission component from the phase boundary at large $`K_T1`$GeV/c,
* no unusually large radius parameters $`R_{\mathrm{inv}}`$ that are reduced due to finite momentum resolution effects,
* $`\lambda _1`$ intercept parameters that get strongly reduced due to finite momentum resolution,
* the reduction of $`R_{\mathrm{inv}}`$ and $`\lambda _1`$ grows with $`K_T`$.

The kaon interferometry measurements at RHIC (at high $`K_T`$) will be, in combination with the pion data, an utmost important and valuable probe of the space-time dynamics (close to the phase boundary) and to the properties of the phase transition.

This work is supported by the Alexander von Humboldt-Foundation through a Feodor Lynen Fellowship (SS) and the Gesellschaft für Schwerionenforschung (Darmstadt). We thank A. Dumitru for many helpful comments, the UrQMD collaboration for permission to use the UrQMD transport model and S. Pratt for providing the correlation program CRAB. SS acknowledges support from DOE Grant No. DE-AC03-76SF00098 and SAB from DOE Grant No. DE-FG02-96ER40945 and DE-AC02-98CH10886.

## References


format: text

text:
# Global metallicity of globular cluster stars from colour-magnitude diagrams

## 1 Introduction

The colour-magnitude diagrams (CMDs) of globular clusters (GCs) present several features to be compared with the various constraints provided by stellar evolution theory. For this reason, they are often used to infer fundamental properties of ancient stellar populations and, in turn, of the early Universe.

It has been known for a long time that the absolute magnitude of the main-sequence turnoff (TO) can be calibrated in terms of the cluster age ($`t`$, in Gyr) and chemical composition parameters ($`Y`$, helium content and $`Z`$, global metallicity). In order to overcome the uncertainties of obtaining accurate GC distances, it is common to use the difference in magnitude between the TO and the zero-age-horizontal-branch (ZAHB) at the RR Lyrae instability strip (log$`T_e`$3.85), whose luminosity level depends on $`Y`$ and $`Z`$, with a negligible dependence on $`t`$. Besides this vertical method to get the cluster age, there is the so-called horizontal method, which is based on the determination of the TO colour relative to some fixed point on the red giant branch (RGB – see Stetson, Vandenberg & Bolte 1996). As a whole, both methods are affected by the empirical difficulties to get the precise position of the TO in the CMD, as well as by several still-open intrinsic uncertainties on the evolutionary models, most notably the dependence on $`Z`$ of the ZAHB luminosity, for the vertical method, and the efficiency of superadiabatic convection in the outer layers (i.e., the mixing length calibration) for the horizontal method, (see Rosenberg et al. 1999, hereafter R99, for details and references).

The determination of $`Z`$ for GC stars is by itself problematic. Under the classical assumption of scaled-solar metal distribution, where all of the abundance ratios of the various metals relative to iron have the solar value, the correlation between the global metallicity $`Z`$ and the measured iron-to-hydrogen ratio \[Fe/H\] is given by the classical relation log$`Z`$=\[Fe/H\]$``$1.70. One should however consider that current empirical scales for \[Fe/H\] (e.g., Zinn & West 1984, Zinn 1985, Rutledge, Hesser & Stetson 1997, and Carretta & Gratton 1997) may differ by up to $``$0.3 dex. In addition, it has to be noticed that the assumption of scaled-solar chemical compositions might be rather inadequate for Galactic GC stars which show a not negligible overabundance of the so-called $`\alpha `$elements ( O, Ne, Mg, Si, Ar, Ca, and Ti) with respect to iron-peak elements (\[$`\alpha `$/Fe\]$``$ 0.3, according to Carney 1996). It follows that the relationship between the measured \[Fe/H\] and the overall metallicity $`Z`$ should follow the relation

$$\mathrm{log}Z=[Fe/H]1.70+\mathrm{log}(0.638f+0.362)$$
$`(1)`$
where $`f`$ is the enhancement factor of $`\alpha `$elements with respect to iron (Salaris et al. 1993)111It is worth recalling that this relation holds strictly only in case of the scaled-solar metal mixture by Ross & Aller (1976). The coefficients of this relation has to be slightly changed when accounting for different solar metal distributions. However, these changes are very small (see Yi et al. 2001)..

The knowledge of the cluster age and global metallicity is fundamental for studying another important feature of CMDs, i.e. the luminosity level of the so-called RGB-bump. This is an intriguing feature of the RGB luminosity function which has been predicted for a long time, since the pioneering theoretical studies by Thomas (1967) and Iben (1968), but whose identification in GCs occurred years later (King et al. 1985). A first comparison (Fusi Pecci et al. 1990) of the observed difference in visual magnitude between RGB-bump and ZAHB level with theoretical models showed that the predicted dependence of $`\mathrm{\Delta }V`$(ZAHB-BUMP) on metallicity was nicely reproduced, but the zero point of the observed relation was $``$ 0.4 mag too faint.

An exhaustive investigation on the dependence of the parameter $`\mathrm{\Delta }V`$(ZAHB-BUMP) on the physical inputs adopted in stellar computations has been performed by Cassisi & Salaris (1997, hereinafter CS97). In the same work, it was clearly shown that a fine agreement between theory and observations does exist when updated RGB stellar models and observational data for GCs with spectroscopic measurements of both iron and $`\alpha `$element abundances, are considered. More recently, the same result was obtained by Zoccali et al. (1999) and Ferraro et al. (1999, hereinafter F99) by using a much larger sample of Galactic GCs than CS97. However, one has also to notice that Bergbusch & Vandenberg (2001) have recently found a discrepancy of the order of 0.25 mag, when comparing the $`\mathrm{\Delta }V`$(ZAHB-BUMP) values provided by their models with the observational data for 4 GCs at various metallicities.

In reality, the RGB-bump luminosity is also dependent on the age of the stellar population ($`\delta V`$(BUMP)/$`\delta `$log$`t`$ 1), while the above investigations do not take into account the age of individual GCs. It is well known the current debate whether the Galactic GCs formed all at the same time or there was a significantly protracted formation epoch. A recent analysis of an homogeneous photometric database containing 34 Galactic GCs (R99) shows that within the intermediate-metallicity clusters ($`1.2`$\[Fe/H\]$``$0.9) there is a clear evidence of age dispersion, with some clusters up to $``$30% younger than the oldest ones. As an example, we consider the globular cluster NGC 362 for which F99 give $`\mathrm{\Delta }V`$(ZAHB-BUMP)=0.10 mag$`\pm `$0.12 mag and a global metallicity222\[M/H\]=\[Fe/H\]+log(0.638$`f`$+0.362). \[M/H\]=$``$0.99.

From the relations given by these authors

$$M_{V,BUMP}=0.75+0.99\mathrm{log}t+1.58[M/H]+0.26[M/H]^2$$
$`(2)`$
$$M_{V,ZAHB}=1.00+0.35[M/H]+0.05[M/H]^2$$
$`(3)`$
one has that the agreement between theory and observations occurs for a cluster age of $``$14.9 Gyr, which is the value derived by R99 for the bulk of coeval GCs, using the Straniero et al. (1997) theoretical isochrones. However, R99 find that the actual age of NGC 362 is $``$ 11.5 Gyr, in which case one obtains a $`\mathrm{\Delta }V`$(ZAHB-BUMP) value by 0.21 mag, namely $``$0.1 mag larger than the observed one.

It should be noticed that R99 adopts the \[Fe/H\] scale by Carretta & Gratton (1997), without any correction due to $`\alpha `$enhancement and, in order to interpret the observed difference $`\mathrm{\Delta }V`$(TO-ZAHB), a linear relation for the absolute magnitude of the ZAHB as a function of the metal content ($`M_{V,ZAHB}`$=0.18(\[Fe/H\]+1.5)+0.65), instead of the quadratic relations used by F99 \[see Eq. (3)\] and CS97 ($`M_{V,ZAHB}`$=1.13+0.39\[M/H\]+0.06\[M/H\]2). In other words, the quoted extensive investigations on the RGB-bump luminosity and the age of Galactic GCs are not mutually consistent, since they use different metallicity scales and adopt, as theoretical values for the RGB-bump and TO luminosity, different calibrations of the ZAHB luminosity in terms of metallicity.

In order to make a complete and self-consistent study of these most prominent features in CMDs, we present in section 2 an homogeneous set of theoretical evolutionary results which predict the TO, RGB-bump and ZAHB luminosity for a wide range of ages (9$`t(Gyr)`$15) and metal abundances (0.0002$`Z`$0.004). In section 3 we introduce a new parameter $`A`$ defined as $`A=\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.566$`\mathrm{\Delta }M_V`$(TO-ZAHB), and we show that this depends only on the global metallicity $`Z`$, being unaffected by age. In order to test this result, the evolutionary models computed by Salaris & Weiss (1997, 1998) and by Vandenberg et al. (2000) are also considered. We show that the former ones are in close agreement with our computations, whereas in the latter case some significant differences are present, leading to a slightly different formulation of the $`\mathrm{\Delta }M_V`$(TO-BUMP) versus $`\mathrm{\Delta }M_V`$(TO-ZAHB) relationship. Nevertheless, the parameter $`A`$=$`\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.487$`\mathrm{\Delta }M_V`$(TO-ZAHB) inferred from the Vandenberg et al.’s (2000) computations is again strongly dependent on global metal abundance with no dependence on age.

On this basis, the TO, ZAHB and RGB-bump observations collected by F99 and R99 for a large sample of Galactic GCs are interpreted. The resulting global metal abundances are presented in section 4, together with some clues to the $`\alpha `$element enhancement with respect to iron in Galactic GCs. A summary is finally given in the last section.

## 2 Theoretical models

In order to understand the observed properties of CMDs and RGB luminosity functions of GCs, we need an homogeneous set of evolutionary models for both the H- and He-burning phases. In present work, we have adopted a set of theoretical models computed by using the FRANEC evolutionary code (Chieffi & Straniero 1989; CS97; Castellani et al. 1997). A subset of these stellar models and corresponding isochrones have been already presented by Cassisi et al. (1998, 1999). However, in order to improve the coverage of the wide metallicity range spanned by Galactic GCs, we have extended previous computations considering additional initial chemical compositions. The physical inputs adopted in the model computations have been already discussed by Cassisi et al. (1998). A short summary is given below.

We use the OPAL (Rogers et al. 1996) equation of state (EOS), supplemented, in the temperature-density regime not covered by the OPAL tables, by the EOS provided by Straniero (1988), plus a Saha EOS in the outer stellar layers. The radiative opacities provided by the OPAL group (Iglesias, Rogers & Wilson 1992) for temperatures larger than 10,000K are employed, together with the Alexander & Ferguson (1994) tables for lower temperatures. Electron conduction opacities come from Itoh et al. (1983). The boundaries of the convective regions are determined according to the canonical Schwarzschild criterion. Induced overshooting and semiconvection during the central He-burning phase have been accounted following the prescriptions given by Castellani et al. (1985). The standard mixing-length theory is used for evaluating the temperature gradient in the superadiabatic regions; the value of the mixing length parameter is fixed at 1.6 for all computations. This value allows to match the RGB effective temperatures of a selected sample of Galactic GCs provided by Frogel, Persson & Cohen (1983).

All the models include microscopic diffusion of both helium and heavy elements, using the formalism by Thoul et al. (1994). For a detailed discussion concerning the effects of atomic diffusion on both main-sequence stars and more evolved models, we refer the interested reader to Castellani et al. (1997) and Cassisi et al. (1998). It is well known that helioseismology provides clear evidence that helium and heavy element settling is occurring in the Sun (see Basu, Pinsonneault, & Bahcall 2000, and references therein). However, recent spectroscopic analyses of GC stars (Gratton et al. 2001; Castilho et al. 2000; Ramirez et al. 2001) show that the surface chemical patterns do not show the effects of atomic diffusion, at least for what heavy elements are concerned. This means that either chemical settling is less efficient in metal-poor stars, or some other physical mechanism at work in the stellar atmospheres counterbalances the effect of atomic diffusion. A firmer understanding of this topic has to wait for more accurate theoretical and observational investigations.

Evolutionary models with masses ranging from 0.6 to 1.3$`M_{}`$ have been computed from the Zero Age Main Sequence up to the tip of the RGB marking the triple-alpha ignition inside the degenerate He core. The explored range of metallicity is $`0.0002Z0.004`$, while for each metallicity the same initial He abundance $`Y`$=0.23 has been adopted. The He-core mass $`M_c`$ and the chemical stratification at the RGB He-flash provided by the different evolutionary models have been used to compute the ZAHB structures needed to estimate the ZAHB luminosity level at log$`T_e`$=3.85, the average effective temperature of the RR Lyrae instability strip. Isochrones have been then constructed in the age range 9 - 15 Gyr.

Finally, all the models and isochrones have been transposed from the theoretical H-R plane to the observational ($`V,BV`$) plane by using the bolometric corrections and colour-temperature relations provided by Castelli, Gratton & Kurucz (1997a,b). In Table 1, for the different chemical compositions adopted in present work, we report the TO, RGB-bump and ZAHB absolute visual magnitude, as a function of age.

## 3 Theoretical calibrations

The data listed in Table 1 show that, at a fixed age, TO, RGB-bump and ZAHB luminosities decrease with increasing the metal content. At a fixed metallicity the TO and RGB-bump become fainter for larger ages, whereas the ZAHB level becomes slightly brighter, due to the slight increase of the He-core mass $`M_c`$ with decreasing total mass of the progenitor on the RGB.

Figure 1 displays the two differential parameters $`\mathrm{\Delta }M_V`$(TO-ZAHB) and $`\mathrm{\Delta }M_V`$(TO-BUMP), as derived from the results in Table 1. The dashed lines show the theoretical trend at a constant age, in steps of 1 Gyr, from 9 Gyr (left) to 15 Gyr (right), while the solid lines refer to a constant $`Z`$ from 0.0002 to 0.004 (bottommost line). The difference $`\mathrm{\Delta }M_V`$(TO-BUMP) increases for larger ages, at fixed $`Z`$, and decreases for larger metallicities, at fixed $`t`$. As for the parameter $`\mathrm{\Delta }M_V`$(TO-ZAHB), it increases with increasing $`Z`$ up to log$`Z`$2.7, and then it decreases for larger metallicities, at least when $`Z`$0.004. On this subject, we recall that all our models are computed with an initial helium abundance $`Y`$=0.23. Within the quoted range of metal content, the expected He increase due to an helium-to-heavy element enrichment ratio $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$2.5 is only $`\mathrm{\Delta }Y`$=0.01, whereas for metallicities larger than $`Z`$=0.004 the original amount of helium is expected to become significantly larger. Furthermore, also the amount of extra-helium brought to the stellar surface by the first dredge up is expected to increase when moving towards larger metallicities. As a consequence, the above theoretical trends cannot be applied to very metal-rich clusters (\[Fe/H\]$`0.6`$).

The data in Fig. 1 disclose that, for ages between 9 and 15 Gyr, the variation of $`\mathrm{\Delta }M_V`$(TO-BUMP) with $`\mathrm{\Delta }M_V`$(TO-ZAHB) has a constant slope of 0.566 (solid lines), with the zero-point depending on $`Z`$. From the results listed in Table 1 we calculate then the mean value of the parameter $`A`$=$`\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.566$`\mathrm{\Delta }M_V`$(TO-ZAHB) for each given $`Z`$. The upper diagram in Fig. 2 shows that the results are well approximated by the following linear relations:

$$\mathrm{log}Z=+0.5181.933A$$
$`(4)`$
$$\mathrm{log}Z=1.1101.039A$$
$`(5)`$
with $`A1.82`$ and $`A<1.82`$, respectively. As it is shown in the lower panel of Fig. 2, the discrepancy between the evolutionary estimate of the global metallicity (hereafter $`Z_{ev}`$) based on the two differential parameters $`\mathrm{\Delta }M_V`$(TO-BUMP) and $`\mathrm{\Delta }M_V`$(TO-ZAHB), and the original value, is close to $`\pm `$0.03 dex.

Before proceeding further, it is worthwhile to compare our evolutionary predictions with similar computations provided by different authors. We considered the recent isochrones and ZAHB models computed by Salaris & Weiss (1997, 1998, hereinafter SW) and by Vandenberg et al. (2000, hereinafter Vdb, see also Bergbusch & Vandenberg 2001) with $`Z`$ ranging from 0.0002 to 0.004 and ages in the range of 10-14 Gyr and 10-16 Gyr, respectively. We wish to notice that both sets of models do not account for atomic diffusion and have been computed by adopting a slightly different $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$ ratio (3 and 2, respectively). Moreover, they are computed under different assumptions about the $`\alpha `$element enhancement. In our comparison we used the Vdb models based on a scaled-solar heavy element distribution.

As shown in the upper panel of Fig. 3, both the ZAHB (open circles) and RGB-bump (filled circles) absolute magnitudes inferred from SW models are in close agreement with our results, whereas the Vdb computations plotted in the lower panel suggest brighter RGB-bump and fainter ZAHB luminosities, with respect to the data in Table 1. The upper panel of Fig. 3 shows also the comparison between the RGB-bump absolute magnitudes provided by the recent Yale isochrones (Yi et al. 2001) and our results (asterisks), for the common values of metallicity ($`Z`$=0.001 and 0.004) and age (10 and 13 Gyrs). As a whole, one has that the Yale isochrones suggest brighter RGB-bumps than our models. Note that we were unable to carry out a similar comparison with the Yale ZAHB absolute visual magnitudes, as no data have been provided in the literature. Figure 4 deals with the TO luminosity provided by the quoted sets of computations, compared with the results in Table 1. We get that SW models are again in agreement with our predictions, whereas VdB and Yale models give brighter absolute magnitudes.

On the basis of the above discussion, we have considered the models by VdB for testing the effect of different input physics on the scenario suggested by our computations. Figure 5 presents $`\mathrm{\Delta }M_V`$(TO-BUMP) versus $`\mathrm{\Delta }M_V`$(TO-ZAHB), as derived from Vdb models. The dashed lines show the theoretical trend at a constant age, in steps of 2 Gyr, from 10 Gyr (left) to 16 Gyr (right), while the solid lines refer to a constant $`Z`$ from 0.0002 to 0.004 (bottommost line). The general trend in Fig. 1 is here confirmed, but the variation of $`\mathrm{\Delta }M_V`$(TO-BUMP) with $`\mathrm{\Delta }M_V`$(TO-ZAHB) has a slope of 0.487 (solid lines). The correlation between the mean value of the parameter $`A`$=$`\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.487$`\mathrm{\Delta }M_V`$(TO-ZAHB) and log$`Z`$ is presented in the upper diagram of Fig. 6, where the solid line depicts the relations

$$\mathrm{log}Z=+1.6172.119A$$
$`(6)`$
$$\mathrm{log}Z=0.1541.330A$$
$`(7)`$
with $`A2.24`$ and $`A<2.24`$, respectively. As shown in the lower panel of Fig. 6, also in this case the discrepancy between the evolutionary estimate of the global metallicity (hereafter $`Z_{ev}`$) and the original value is close to $`\pm `$0.03 dex.

## 4 Global metallicity for globular clusters in the Galaxy

In this section we provide accurate evolutionary estimates of the global metallicity of a large sample of Galactic GCs, by using the relationships discussed in the previous section.

The sample of Galactic GCs studied in this paper is listed in Table 2. The \[Fe/H\] values are listed both for the Zinn & West (1984) and the Carretta & Gratton (1997) metallicity scale, while the TO visual magnitudes are from R99, except for the following clusters: NGC 4833 and 5286 (Samus et al. 1995a,b), NGC 5694 (Ortolani & Gratton 1990), NGC 6584 (Sarajedini & Forrester 1995), NGC 6717 (Ortolani et al. 1999), NGC 6934 (Piotto et al. 1999) NGC 7006 (Buonanno et al. 1991) and NGC 7492 (Buonanno et al. 1987). The ZAHB and RGB-bump observations are from F99, except for the following additional measurements of the RGB-bump level: NGC 6352 and 6397 (Alves & Sarajedini 1999) and NGC 6362 (Piotto et al. 1999). We wish to notice that the two major sets of observations are mutually consistent as the ZAHB visual magnitudes listed by F99 are $``$ 0.07 mag fainter than the HB visual magnitudes in R99, except for NGC 6093 and NGC 6254. For these two clusters, the F99 values are brighter than those given by R99 (see later). Moreover, we have also checked that the ZAHB magnitudes given by F99 are consistent with the additional measurements of the TO and RGB-bump visual magnitudes.

As a first step, we determine the values of $`A`$=$`\mathrm{\Delta }V`$(TO-BUMP)-0.566$`\mathrm{\Delta }V`$(TO-ZAHB) from the data in Table 2, then we derive the global metallicity using Eq. (4) and Eq. (5). The results for $`Z_{ev}`$ are listed in column (4) of Table 3, together with the horizontal branch morphological type (HB-type) 333The HB-type is given by the ratio (B-R)/(B+V+R), where B and R denote the number of HB stars blueward and redward of the RR Lyrae instability strip, respectively, and V is the number of RR Lyrae stars. from Harris (1996). For the two clusters NGC 6093 and 6254, the last two lines in Table 3 refer to the HB level given by R99.

The correlation between log$`Z_{ev}`$ and the measured \[Fe/H\] parameter444Note that the measured \[Fe/H\] values are taken with an uncertainty of $`\pm `$0.10 dex. is presented in the upper panel of Fig. 7, where the solid line shows the canonical correlation log$`Z`$=\[Fe/H\]$``$1.7 for scaled-solar chemical compositions. One obtains that almost all the clusters present a mild overabundance of the global metallicity with respect to the measured iron-to-hydrogen content. In terms of $`\alpha `$element enhancement \[see Eq. (1)\], we show in the bottom panel in Fig. 7 that the clusters with \[Fe/H\]$``$1.7 show $`f`$ 1.5, whereas the four most metal-poor clusters show $`f`$ 3 (dashed lines).

Figure 8 displays the results based on Vdb computations. These models suggest slightly larger $`Z_{ev}`$ values (see data listed in column (5) of Table 3) in comparison with our models, with an average difference of 0.08$`\pm `$0.06 dex which is within the uncertainty of the derived global metallicity. Also the $`\alpha `$element enhancement is slightly increased, clusters with \[Fe/H\]$``$1.7 now showing $`f`$ 1.9, whereas the four most metal-poor clusters show $`f`$ 3.7 (dashed lines).

We remark here that the global metallicity given in Table 3 are derived from the differential parameters $`\mathrm{\Delta }V`$(TO-ZAHB) and $`\mathrm{\Delta }V`$(TO-BUMP) and are fully independent of the measured iron-to-hydrogen content. On the contrary, the f enhancement factors plotted in the lower panel of Figs. 7 and 8 obviously depend on the adopted empirical \[Fe/H\] scale. When we consider the Zinn & West (1984) scale the results are quite different. As shown in Fig. 9 and Fig. 10, the $`\alpha `$enhancement factor is $`f`$ 3 (our models) or $``$ 3.5 (Vdb models) with \[Fe/H\]$``$1.3, declining toward the solar value $`f`$1 at larger metallicities.

The abundances of $`\alpha `$elements in Population II stars provide fundamental constraints to the chemical evolution models of the Galaxy. Several spectroscopic studies have been devoted to this problem (see Carretta, Gratton & Sneden 2000 for discussion and references). As a whole, $`\alpha `$elements appear to be overabundant by a constant factor $``$2-3 in metal-poor stars, declining toward scaled solar ratios around \[Fe/H\]=$``$1.0. As the run of \[$`\alpha `$/Fe\]555\[$`\alpha `$/Fe\]=log$`f`$. with \[Fe/H\] is required to study the chronology of the Galactic halo formation (Matteucci & Francois 1992), we notice that mild variations of the empirical \[Fe/H\] scale may lead to quite different pictures for the dependence of \[$`\alpha `$/Fe\] on \[Fe/H\], and different \[Fe/H\] values for the transition toward scaled-solar ratios.

Both the distribution and absolute values of the $`\alpha `$element enhancement suggested by our method together with the Zinn & West (1984) \[Fe/H\] scale appear in fine agreement with the empirical results by Carney (1996). On the contrary, the results obtained using the Carretta & Gratton (1997) \[Fe/H\] scale do not appear fully supported by current spectroscopic measurements. Even though to investigate the accuracy of current empirical \[Fe/H\] scales is out of the aim of the present paper, we wish to emphasize that this is still a controversial issue. We also wish to notice that recently Vandenberg (2000) and Bergbusch & Vandenberg (2001) have provided some additional evidence, based on the comparison between their theoretical isochrones and CMDs of Galactic GCs, that the metallicities determined by Carretta & Gratton (1997), mainly for the intermediate-metallicity clusters, are too high, values close to the Zinn & West (1984) scale being more favourite. Present analysis provides independent support to their claim.

Before closing this section, it is worthwhile to emphasize that our method, being based only on well defined observational features of the CMD, appears quite suitable for deriving accurate estimates of the GC global metallicity. It uses only differential quantities, i.e. visual magnitude differences, and its results are therefore insensitive to uncertainties on the zero point of the CMD photometric calibration as well as on the cluster distance and reddening.

It is also worth noticing that the method relies on observables whose properties are not significantly affected by non-canonical processes occurring along the RGB, like deep-mixing processes, which are currently invoked to explain the chemical abundance patterns observed in RGB stars of many Galactic GCs (Sweigart & Mengel 1979, Langer, Hoffman & Sneden 1993). As suggested by theoretical predictions (Mestel 1957) and spectroscopic data (Charbonnel, Brown, & Wallerstein 1998, and references therein), molecular weight gradients, such as the H-discontinuity left over in the envelope of RGB stars by the first dredge up, should be able to inhibit these non-canonical mixing processes, which are therefore expected to be efficient only after the RGB bump.

## 5 Summary

The main results obtained in present work can be summarized as follows:

* We have developed an homogeneous theoretical scenario, based on updated stellar models including atomic diffusion, computed for a wide metallicity range; this allows us to predict the behaviour with the cluster metallicity and age of relevant observational features of the CM diagram, such as the Turn Off magnitude and the brightness of both the ZAHB at the level of the RR Lyrae instability strip, and the RGB-bump;
* we define a new observable based on the visual magnitude difference between the TO and the ZAHB $`\mathrm{\Delta }M_V`$(TO-ZAHB), and the TO and the RGB-bump $`\mathrm{\Delta }M_V`$(TO-BUMP), namely, $`A=\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.566$`\mathrm{\Delta }M_V`$(TO-ZAHB). We show that the parameter $`A`$ does not depend at all on the cluster age, whereas it does strongly correlate with the global cluster metallicity. The calibration of the parameter $`A`$ as a function of $`Z`$ is provided;
* by using this calibration, together with a large observational database of Galactic GCs, we estimate the global metallicity of a large sample of GCs. Our estimates are independent of empirical metallicity measurements, but strongly rely on the accuracy and reliability of the adopted theoretical models. By using our global metallicity determinations, together with presently available \[Fe/H\] determinations (Zinn & West 1984 or Carretta & Gratton 1997), and using the relation provided by Salaris et al. (1993) that connects global metallicity to \[Fe/H\] and $`\alpha `$element enhancement, we are able to provide an estimate of the \[$`\alpha `$/Fe\] ratio for the individual clusters in our sample;
* the trend of \[$`\alpha `$/Fe\] with \[Fe/H\] we predict when using the Zinn & West metallicity scale, appears in fine agreement with current spectroscopical measurements, which suggest a rather constant overabundance of $`\alpha `$elements in metal-poor stars, and a decline toward the solar ratios around \[Fe/H\]$``$1.0. This picture is not supported when we use the Carretta & Gratton \[Fe/H\] scale;
* the calibration of the $`AZ_{ev}`$ relation is model dependent, in the sense that the Vdb models, which are based on different input physics, yields a slightly different formulation of the parameter $`A`$, given by $`A`$=$`\mathrm{\Delta }M_V`$(TO-BUMP)$``$0.487$`\mathrm{\Delta }M_V`$(TO-ZAHB). However, the resulting $`Z_{ev}`$ values are only 0.08$`\pm `$0.06 dex larger than the results based on our computations;
* the good agreement between the global metallicity estimates provided by our models accounting for atomic diffusion, and the ones obtained by employing Vdb models which neglect atomic diffusion, shows that the theoretical calibration of the $`A`$ parameter is only marginally affected by the efficiency of this non-canonical process;
* the present analysis clearly demonstrates that the our newly defined observational parameter, independent of the uncertainty on the zero point of the photometric calibration as well as on the cluster distance and reddening, is a powerful indicator of globular cluster global metal abundances.

Acknowledgments:

We warmly thank D. Vandenberg, P. Bergbusch and M. Salaris for sending us data about their theoretical models. We gratefully also acknowledge M. Salaris for an accurate reading of an early version of this manuscript and for his helpful comments. We thank the anonymous referee for several helpful suggestions and comments. One of us (C.S.) wishes to acknowledge the Max Planck Institut für Astrophysik for the kind hospitality during the early phase of this work. Financial support for this work was provided by MIUR-COFIN2000, under the scientific project ”Stellar observables of cosmological relevance”.


text:
# Critical Behaviour of the Drossel-Schwabl Forest Fire Model

## I Introduction

Empirical analyses suggest that power laws with anomalous exponents are ubiquitous in nature, ranging from $`1/f`$ noise and earth quake distributions to fractal coast lines, species extinction rates, weather records, and the statistics of DNA mandelbrot ; bak ; buldyrev ; koscielny . On the other hand, it is well known that such scaling laws – most clearly seen by linear relationships in log-log plots – can be spurious. Log-log plots have a notorious tendency to suggest linear curves, even if there are no real power laws.

It would thus be extremely useful if these empirical observations could be backed by theoretical models where power laws can be either proven exactly or at least verified beyond doubt by high statistics simulations. Unfortunately, equilibrium systems in general show anomalous scaling laws only at critical points which are codimension one phenomena: One has to fine tune some parameter (e.g. temperature) to reach them, otherwise no power laws are obtained. Thus, they cannot be used to justify why such power laws should be seen in nature.

A possible solution of this puzzle was indicated in btw ; bak  where it was suggested that many non-equilibrium systems could be driven by their dynamics into a critical state. The main ingredients of this self organized criticality (SOC) is slow driving towards some instability and a mechanism to relax the tensions built up by the drive locally and partially. Since the tensions are not relaxed completely (in many models they are just redistributed), the state becomes marginally stable and apt to relaxation events on wider and wider length scales, which then lead to anomalous scaling laws. The paradigmatic model is the “sand pile” of btw  which does not describe real sand piles but which was proven exactly to show anomalous scaling laws at least for some of its observables dhar ; priezzhev ; glkp .

Another model which was proposed to show SOC is the forest fire model introduced independently by Henley henley  and by Drossel and Schwabl ds . In this lattice model with discrete time each site can be either occupied by a tree or empty. New trees are grown with small fixed rate $`p`$ on empty sites, and with a rate $`fp`$ sites are hit by lightning strokes which then burn the entire cluster of trees connected to this site. Burning happens infinitely fast, so the only relevant parameter is the ratio $`\theta =p/f`$ which also sets the scale for the average fire size (the number of trees burnt after one lightning). Criticality is observed in the limit $`\theta \mathrm{}`$.

The Drossel-Schwabl model (called DS model in the following) is different from other SOC models in two ways.

* It involves not only the separation between two time scales (the slow build-up of stress and the fast relaxation), but involves three time scales: The fast burning of connected clusters of trees, the slow re-growth, and the even slower rate of lightnings.
* The growth of trees does not lead to a state which is inherently unstable (as does the addition of sand grains to the top of a pile does), but only to a state susceptible to being burnt. Without lightning, the tree density in any patch of forest can go far beyond criticality. When the lightning strikes finally, the surviving trees have a density far above critical.

Indeed, it was observed in grass  that the stationary steady state of the model is not “critical” in most regions, in the sense that its tree density is not marginal for the spreading of fire. Rather, it is composed of large patches of roughly uniform density, most of which are either far below or far above the critical density for spreading. Nevertheless, power laws were observed for several observables, partly because these patches occur with all sizes, so that also the fires had a broad spectrum of sizes.

While normal scaling with mean field exponents had been seen in ds , all subsequent simulations grass ; henley ; christensen ; clar ; honecker  showed clear signs of anomalous scaling:

* The fire size distribution scaled, for small $`s`$ ($`s`$ is the number of trees burnt in one fire) and in the limit $`\theta \mathrm{}`$, as $`P(s)s^{1\tau }`$ with $`\tau 2.15`$;
* For finite $`\theta `$, $`P(s)`$ is essentially cut off at $`s_{\mathrm{max}}\theta ^\lambda `$ with $`\lambda 1.1`$;
* The average rms. radius $`R(\theta )=R^2^{1/2}`$ of all fires scaled as $`\theta ^\nu `$ with $`\nu 0.58`$. Here, $`R^2`$ is the Euclidean distance of a burning tree from the site of lightning, and the average is taken over all trees in all fires at fixed $`\theta `$; and
* The rms. radius $`R(s,\theta )`$ of fires of fixed size $`s`$ scaled as $`s^{1/D}`$ where the fractal dimension $`D`$ of fires is $`D1.96`$ clar  to $`2.0`$ grass .

Finally, the average density of trees was found to be $`\rho =0.408\mathrm{const}/\sqrt{\theta }`$.

There were however seen already at that time large corrections to this scaling law and deviations from conventional ansatzes. Thus,

* The determination of $`\tau `$ was subject to large systematic uncertainties grass ;
* The scaling $`Rs^{1/D}`$ was observed in grass  only for fires of typical size ($`s\theta `$). For large $`s`$, $`R`$ was significantly larger;
* For finite $`\theta `$, the fire distribution $`P(s)`$ did not follow the normal scaling ansatz $`P=s^{1\tau }\varphi (s/s_{\mathrm{max}})`$ grass ;
* There are (at least) two differently divergent length scales honecker : The correlation length evaluated from all pairs of sites scales differently with $`\theta `$ than $`R(\theta )`$ footnote1 ;
* Finite size behaviour is abnormal schenk0 .

In two recent publications, these problems were taken up again. In vespignani  it was claimed that they are due to non-leading corrections to scaling. In schenk  a more radical solution was proposed with two distinct classes of fires which both contribute to the scaling limit. In the latter paper also a connection to ordinary percolation was proposed, and the conclusion was backed by a “coarse grained” model which supposedly could mimic the DS model at extremely large $`\theta `$ and on extremely large lattices.

It is the purpose of this paper to present very large simulations which show quite unambiguously that none of these describe really the true critical behaviour of the DS model. While our simulations agree perfectly with previous ones for the lattice sizes and $`\theta `$ values used there, we shall see that the supposed scaling laws are just transients. Even the present simulations do not reach the true asymptotic regime, but some suggestive features of the true asymptotics do emerge.

We describe our simulations and their results in the next section. In the last section we draw our conclusions.

## II The Simulations

Our simulations are straightforward, with a few subtleties. They follow Ref.grass  in taking $`p0`$ by making no attempts to grow new trees while a cluster burns. As in grass  we made exactly $`\theta `$ growth attempts on randomly chosen sites between two successive lightnings, instead of letting this number fluctuate. In the large $`L`$ limit this should not make any difference. But in contrast to grass , where a depth first algorithm had been used, we used a breadth first algorithm to burn the cluster footnote2 .

In order to simulate very large lattices with minimal storage, we use multi-spin coding, i.e. we use only one bit to store the status of a site. In this way we could simulate lattices of size $`L\times L,L=65536`$ on computers with 1 GB main memory. Notice that we do not need to store for every tree whether it is burning or not, since the burning trees are stored in a separate list which is updated at each time step. Boundary conditions were helical, i.e. sites are indexed by one scalar index $`i`$, with neighbours having indices $`i\pm 1`$ and $`i\pm L`$, and with $`i+L^2i`$. The largest previous simulations vespignani  had used $`L=19000`$.

We were careful to discard sufficiently long transients (between $`1.6^6`$ and $`1.2\times 10^7`$ lightnings; this is up to one order of magnitude longer than those in honecker ; vespignani ) before taking statistics. This is needed since excessively large fires occur during these transients, and thus the tail of the distribution $`P(s)`$ is heavily influenced by them. We believe that previous analyses were affected by this problem which is easily overlooked since bulk properties (such as $`\rho `$ or $`P(s)`$ for typical $`s`$) show much faster convergence. The total number of fires in each run used for averaging was between $`10^9`$ (for $`\theta 250`$) and $`9.3\times 10^6`$ (for $`\theta =256000`$; see Table 1). Previous authors vespignani  went only to $`\theta =32768`$ with $`2.5\times 10^7`$ fires. Compared to that, our statistics is larger by roughly 1 order of magnitude. All simulations were done on fast Alpha workstations. The CPU times per run varied between 15h and 4 weeks. As random number generator we used Ziff’s four-tap generator with period $`2^{9689}1`$ ziff .

We tested for finite size effects by making for the same $`\theta `$ several runs on lattices of different sizes. Previously, finite size effects had been studied systematically in schenk0 , but only for much smaller lattices ($`L2000`$; the authors of vespignani  called their analysis a finite size analysis, but they actually made a conventional scaling analysis without checking the finite size behaviour). For $`\theta =64000`$, e.g., we made runs with $`L=2^{14},2^{15},`$ and $`2^{16}`$. We verified that distributions like $`P(s)`$, $`R(s,\theta )`$, or $`P(t)`$ ($`t`$ is the burning time of a fire) were independent of $`L`$ within the statistical accuracy, i.e. any systematic $`L`$-dependence was masked by statistical fluctuations.

Systematic dependencies were seen only for averaged quantities like $`s`$, $`\rho `$, or $`R(\theta )`$. Indeed, $`\rho `$ can be measured either immediately before or immediately after a lightning. Since each lightning burns in average $`s`$ trees (here $`s`$ is averaged over all lightnings, whether they hit a tree or an empty site), we have $`\rho _{\mathrm{before}}\rho _{\mathrm{after}}=s/L^2`$. We found that $`\rho _{\mathrm{before}}`$ decreased with $`L`$, while $`\rho _{\mathrm{after}}`$ increased with it (see Table 1). More precisely, we have

$$\rho _{\mathrm{after}}(L)\rho (\mathrm{})0.6s/L^2.$$
(1)
Notice that $`s`$ is given for $`L\mathrm{}`$ exactly by $`s=(1\rho )\theta `$ ds . For finite $`L`$, one has to replace just $`\rho `$ by $`(\rho _{\mathrm{before}}+\rho _{\mathrm{after}})/2`$, to leading order in $`1/L`$. This was verified in the simulations, but it tests just stationarity and the absence of gross mistakes.

The rms. radius of fires depended on $`L`$ and $`\theta `$ in a more complicated way, see Table 1. The data were less clear in that case, but they could be fitted by

$$R(\theta )_{L=\mathrm{}}R(\theta )_L0.004\theta ^{3/2}/L.$$
(2)
Results for the fire size distributions $`P(s)`$ are shown in Fig.1. We see the approximate power decay $`s^{1\tau }`$ for $`s<s_{\mathrm{max}}`$ with $`s_{\mathrm{max}}`$ roughly proportional to $`\theta `$. But we see also the strong deviations from this power law first observed in grass . Due to these deviations, $`P(s)`$ decreases with $`\theta `$ in the scaling region (making the effective value of $`\tau `$ increase with $`\theta `$), but has a growing bump at $`ss_{\mathrm{max}}`$. In schenk  it was conjectured that the asymptotic value of $`\tau `$, estimated from the scaling region $`1s\theta `$ in the limit $`\theta \mathrm{}`$, is $`\tau =2.45`$. This was based on heuristics and on simulations for small $`\theta `$ on very small lattices ($`L=1300`$). In order to test it, we plotted $`P(s)`$ in Fig.1b after multiplying it with a suitable power of $`s`$. It is seen that $`s^{1.19}P(s)`$ becomes flat in a wider and wider region of $`s`$ as $`\theta `$ increases. For very small fires ($`s<300`$) $`P(s)`$ decreases faster than $`s^{1.19}`$, but there is no indication that the slope increases with $`\theta `$ for $`\theta >8000`$. Based on this evidence we would thus conclude that $`\tau =2.19\pm 0.01`$, ruling thereby out the value of schenk  – provided, of course, that there is no change of behaviour as $`\theta `$ becomes even larger.

In grass  it was conjectured that the bump near $`ss_{\mathrm{max}}`$ is due to the cut-off. If the integrated distribution $`P_{\mathrm{int}}(s)=_{s^{}=s}^{\mathrm{}}P(s^{})`$ were just a power multiplied by a sharp cut-off, its derivative $`P(s)`$ would have a bump where the cut-off sets in. This bump would consist of those events which would have been in the tail which is cut off. If this were right and $`P_{\mathrm{int}}(s)`$ would indeed show normal scaling, we should expect the height of the bump in Fig.1b to be independent of $`\theta `$. But this is obviously not the case. Instead, it increases with $`\theta `$, suggesting a different scaling law $`s^{1\tau ^{}}`$ for the envelope of the curves in Fig.1a, with $`\tau ^{}=2.111\pm 0.006`$.

A similar conclusion is reached by looking directly at $`P_{\mathrm{int}}(s)`$. Indeed, log-log plots of $`P_{\mathrm{int}}(s)`$ versus $`s`$ are much more straight in the scaling region (Fig.2a). But again a blow-up after multiplication by a suitable power of $`s`$ shows that this is misleading (Fig.2b). Even $`P_{\mathrm{int}}(s)`$ is not convex for $`\theta >1000`$, and it develops an increasingly sharp shoulder near $`s=s_{\mathrm{max}}`$ as $`\theta \mathrm{}`$. We can try to obtain an alternative estimate of $`\tau `$ by fitting straight lines such that they touch both maxima in Fig.2b. Results of this are shown in Fig.3. In contrast to the previous estimate of $`\tau `$ they would indicate that the effective $`\tau `$ decreases with $`\theta `$ and is clearly less than 2.19 for $`\theta \mathrm{}`$ (also the value 2.159 of honecker  seems hardly compatible with the extrapolation of Fig.3 to $`\theta \mathrm{}`$). But a convergence to $`2.11`$ seems quite possible, as we would expect if the bumps in fig.1b have a width independent of $`\theta `$ (which is not excluded by the data).

This obviously means that these scaling violations are not corrections which disappear in the limit $`\theta \mathrm{}`$, as was claimed in vespignani . It rather indicates that the conventional scaling picture, and in particular the usual ansatz

$$P(s)=s^{1\tau }\varphi (s/s_{\mathrm{max}}),$$
(3)
is basically wrong, as was already conjectured in grass ; schenk .

Previous analyses indicated that $`s_{\mathrm{max}}`$ actually increased faster than $`\theta `$, roughly as

$$s_{\mathrm{max}}\theta ^\lambda $$
(4)
with $`\lambda 1.08`$. We verified this qualitatively, but verified also the finding of grass  that $`s_{\mathrm{max}}`$ and thus also $`\lambda `$ are not well defined since the sharpness of the cut-off increases with $`\theta `$. This is already obvious from Fig.2b, but it persists for larger values of $`s`$ not shown in this plot. It can also be seen by making copies of Fig.2a on transparencies and overlaying them.

Better scaling than for $`P(s)`$ was seen in grass  for $`R(\theta )`$. Our present data are fully compatible with those of grass ; clar ; honecker ; footnote1 , but involve much higher statistics and cover a wider range of $`\theta `$. Thus it might not be too surprising that we now do not see perfect scaling any more. But the observed deviations from a pure power law (see Fig.4) are much larger than expected from subasymptotic corrections footnote3 . They clearly show that the previously seen power law was spurious and does not describe the asymptotic behaviour. A power law fit through the last two points would give $`R(\theta )\theta ^{0.563}`$, but this is obviously not yet the asymptotic behaviour. Again, as for $`P(s)`$ and $`P_{\mathrm{int}}(s)`$, we cannot yet say what the correct asymptotic behaviour will be. In any case, the claim of honecker  that there are two different diverging correlation lengths becomes obsolete.

An indication of the origin of all these puzzles comes from looking at $`R(s,\theta )`$. This is defined analogously to $`R(\theta )`$, except that the averaging is now done over all fires of fixed size $`s`$. We expect $`R(s,\theta )s^{1/D}`$ if fires are fractal with dimension $`D`$. Previous analyses had given $`D2`$ (or slightly less). Thus we do not plot $`R(s,\theta )`$ versus $`s`$, but $`R(s,\theta )/s^{1/2}`$. From Fig.5 we see that there are 3 distinct regions which by and large coincide with the three regions (left bump, flat middle part, and right peak) in Fig.1b.

* Region I: These are very small clusters, of size $`s<300`$. For them $`R(s,\theta )`$ increases faster than $`\sqrt{s}`$, which might suggest $`D<2`$. But we rather interpret this as a finite-cluster artifact. Indeed, even compact clusters (with $`D=2`$) with a crumpled boundary will show an effective fractal dimension $`<2`$. These small clusters arise from lightnings which hit either very small regions with supercritical tree density which burn off completely, or regions of subcritical density on which the fires form subcritical percolation clusters. Both mechanisms give compact clusters with fuzzy boundary. In any case, region I becomes less and less important as $`\theta `$ increases.
* Region II: This is roughly equal to the scaling region in fig.1b. Here $`R(s)`$ is nearly proportional to $`\sqrt{s}`$, i.e. $`D`$ is very close to 2 – but not quite. Also, we see a clear decrease of the minimal slope in Fig.5 with $`\theta `$, but it seems not to be sufficient to give $`D=2`$ in the limit $`\theta \mathrm{}`$ (see Fig.6). According to Fig.6 the critical exponent $`\nu `$ defined by $`R(s,\theta )s^\nu `$ in region II seems to converge to 0.505 for $`\theta \mathrm{}`$, corresponding to $`D=1.98`$. This is half way between the best previous estimates $`D=1.96`$ clar  and $`D=2.0`$ grass . It indicates that clusters in Region II are fractal, but more compact than critical percolation clusters ($`D=1.89`$).
* Region III: In the region of very large fires, corresponding to the increasing bumps in Fig.1b, the apparent fractal dimension decreases again. Unfortunately, due to the strong curvatures of the curves in Fig.5 in this region, we cannot quote any definite dimension value. But by plotting $`R(s,\theta )/s^\nu `$ with suitable values of $`\nu `$ on log-log scales, we see that (i) the maximal values of $`\nu `$ increase slowly with $`\theta `$ (except for the very largest $`\theta `$ where statistics is poor); and (ii) for our largest $`\theta `$ and $`s`$, we have $`\nu 0.64`$. Thus the largest fires are definitely more fractal than critical percolation clusters!

In terms of the scenario with roughly homogeneous patches with constant tree density grass ; clar ; schenk , fires in region II correspond to single patches of typical size which are hit by lightning just when their density has reached about the critical percolation threshold. If lightning would always strike exactly at criticality, and if all trees in the patch would have burnt during the last fire so that all trees now are placed randomly, this would give $`D=1.89`$. Fires in region II are more compact mainly because it will take some time until a lightning strikes by chance the region, after it has reached the critical density.

The larger a patch is, the bigger will be its perimeter and thus also the chance that a fire “spills over” and burns also a neighbouring patch, and from this also a next patch, etc. Since this will happen only along parts of the perimeter, the resulting fires will be rather fuzzy, with effective dimension $`<2`$. We propose that such fires dominate in region III. Although they are rather few in number, they are very important since they burn large parts of the entire lattice, and they lead to rearrangements of the global pattern of patches. Notice that fires in region III burn only small parts of the entire region they cover, leaving behind more unburnt trees than fires of type II. Since they dominate more and more as $`\theta `$ increases, one might suspect that this leads to an increase of $`\rho `$ with $`\theta `$.

Indeed, a slight increase of $`\rho `$ with $`\theta `$ had been seen in all previous simulations. The best previous estimates were $`\rho =0.4075const/\sqrt{\theta }`$ grass  (unfortunately, the constant multiplying $`\theta ^{1/2}`$ was estimated wrongly in grass  due to a simple mistake) and $`\rho =(0.4084\pm 0.0005)const/\theta ^{0.47}`$ honecker ; vespignani . The results of our present simulations are shown in Fig.7. For $`\theta <10^4`$ we see a perfect agreement with previous results, but for larger $`\theta `$ there is dramatic disagreement: Our measured values are higher than predicted by extrapolation from small $`\theta `$, by up to 100 standard deviations. It is not clear why this was missed in vespignani . But we might mention that no data for $`\theta >10000`$ are shown in Fig.1 of vespignani , although the authors claimed to have made high statistics measurements up to $`\theta =32768`$.

When plotted against $`\mathrm{log}\theta `$, our values of $`\rho `$ follow roughly a straight line for $`\theta 4000`$. We should of course not take this as the asymptotic behaviour, since $`\rho `$ can never increase beyond $`p_c=0.5927\mathrm{}`$ which is the critical value for site percolation in 2 dimensions. But we can use it to obtain a very crude estimate of the order of magnitude when $`\rho p_c`$ should be reached. It is $`\theta 10^{40}`$.

Although any prediction based on such a large extrapolation should be taken with great care, we conjecture that indeed $`\rho `$ converges to $`p_c`$ for $`\theta \mathrm{}`$. The main reason is that we see no other plausible scenario compatible with our present numerics. It is not clear whether fires in this limit correspond to critical percolation clusters, since we must expect that weak correlations in the tree densities survive in this limit, sufficiently so to spoil any agreement with uncorrelated percolation on large scales.

A last hint in favour of our claim that large fires asymptotically are dominantly associated to regions of critical tree densities (i.e., $`\rho p_c`$), even though they might not form critical percolation clusters, is obtained by studying the mean sizes of fires which started in regions of given tree density. If a lightning hits a tree at site $`i`$, we define the local tree density $`\rho _{\mathrm{loc}}`$ at this site as the number of other trees in the surrounding square of size $`9\times 9`$ divided by 80. This size is of course rather arbitrary, but we obtained qualitatively similar results for squares of sizes $`7\times 7`$ and $`5\times 5`$. In Fig.8 we plot the average fire size $`s(\rho _{\mathrm{loc}})`$, divided by the overall mean $`s`$, against $`\rho _{\mathrm{loc}}`$. For small $`\theta `$ we see a monotonic increase which is easily understood: Since even the largest patches of uniform density are not much larger than the square, large fires can only result from regions with high local density. This is no longer true for large $`\theta `$. There, patches with very high density are probably small, otherwise they would already have been burnt down earlier. At values of $`\theta `$ reached in this work the largest $`s(\rho _{\mathrm{loc}})`$ are still for large $`\rho _{\mathrm{loc}}`$, but a pronounced peak at $`\rho p_c`$ develops where $`s(\rho _{\mathrm{loc}})`$ has a local maximum. Since also the number of sites with $`\rho _{\mathrm{loc}}p_c`$ is much larger than those with $`\rho _{\mathrm{loc}}p_c`$ (see next paragraph), we see that it is fires in regions with critical density which play an increasing important role as $`\theta \mathrm{}`$. In any case, the very strong dependence of $`s(\rho _{\mathrm{loc}})`$ on $`\theta `$ shows that we are still far from the asymptotic region where we expect this dependence to disappear.

Finally, we show in Fig.9 the distribution of local densities itself. It depends rather weakly on $`\theta `$ (less than $`s(\rho _{\mathrm{loc}})`$, at least), but the precise way it does depend on $`\theta `$ is rather surprising and not yet fully understood. First of all, it develops an increasingly sharp maximum which slowly shifts from $`\rho _{\mathrm{loc}}0.4`$ to $`\rho _{\mathrm{loc}}0.6`$ as $`\theta `$ increases. This is to be expected after our previous observations. What is unexpected and hard to explain is a shoulder at $`\rho _{\mathrm{loc}}0.8`$ which develops for the largest $`\theta `$ values. It is not very large but statistically highly significant (it was seen in all runs with $`\theta 32000`$). Presumably related to it is a shoulder at small $`\rho _{\mathrm{loc}}`$ ($`0.1`$) which also increases with $`\theta `$: If patches with density $`>0.7`$ burn, they leave behind extremely strongly depleted patches. One possible reason why such high density patches can survive at all is that many large fires are fractal and leave behind small disconnected regions of fairly high but subcritical density. These regions then are too small to have a chance to be hit by lightning until their density has grown far beyond $`p_c`$.

## III Discussion

The simulations reported in this paper leave little doubt that all scenarios proposed so far for the Drossel-Schwabl forest fire model do not describe the asymptotic behaviour for $`\theta \mathrm{}`$, where the model should show SOC according to the standard folklore. Indeed, we do not see much indications for any power laws in this model, as all proposed scaling laws seem to be just transient. There are a number of observables which do show straight lines on log-log plots (such as Fig.1a in the central region of $`s`$ or the envelope to Fig.1a), but it seems more likely that also these are spurious.

This situation is of course not altogether new. There are a number of models which were supposed to show anomalous scaling, but closer scrutiny proved otherwise. A good example is the Bak-Chen-Tang forest fire model bct  which at least in $`d=2`$ is not critical gk . Other examples include the Newman-Sneppen model newman  where one can prove semi-numerically that the observed power laws are transient grass-unpub  and maybe even the “classical” abelian sandpile model in $`d=2`$. While power laws for waves were proven rigorously in that model, it might well be that the observed deviations from finite size scaling demenech ; tebaldi  do not herald multifractality but just simply no scaling at all. One indication for the latter is the fact that some scaling laws show violations which do not seem do vanish for increasing system sizes glkp . Also, some other quantities in the sandpile model which involve superpositions of many waves depend qualitatively on the geometry of the lattice (square vs. strip) grass-unpub . For a system with true scaling one would not expect this.

The situation becomes even worse when going to real life phenomena. It does not seem unlikely that many of the observed scaling laws are just artifacts or transients. Problems of spurious scaling in models which can be simulated with very high precision such as the present one should be warnings that not every power law supposedly seen in nature is a real power law.

Acknowledgement: I am indebted to Walter Nadler for critically reading the manuscript.

Note added: After finishing this paper, I was informed that similar conclusions had been reached in a recent paper pruessner , where the authors studied $`P(s)`$ with comparable (or even larger) statistics but on somewhat smaller lattices and for somewhat smaller $`\theta `$ than in the present paper. Unfortunately, neither the spatial sizes of fires nor precise values of $`\rho (\theta )`$ were measured in pruessner . I am indebted to Gunnar Pruessner for sending me this paper, for very helpful correspondence and for pointing out some misprints in my manuscript.


text:
# 1 Introduction

## 1 Introduction

It has been shown that various observables of relativistic heavy-ion collisions can be well described by statistical-thermal or hydrodynamical models. In such a way, a selected subset of a large number of observables can be reproduced by a small number of characteristic parameters, such as temperature, density or flow velocity. It is the subject of the present note to pursue this idea and to analyze the centrality dependence of the thermal parameters describing hadron multiplicities, hadron spectra and dilepton spectra. This will provide further information about the effects of the size of the excited strongly interacting system and help in the systematic understanding of the experimental data.

## 2 Hadron Multiplicities

Hadron multiplicities can be reproduced  by the grand-canonical partition function $`𝒵_i(V,T,\mu ^\alpha )=\text{Tr}\left[\mathrm{exp}\{(\widehat{H}\mu ^\alpha Q_i^\alpha )/T\}\right]`$, where $`\widehat{H}`$ is the statistical operator of the system, $`T`$ denotes the temperature, and $`\mu ^\alpha `$ and $`Q_i^\alpha `$ represent the chemical potentials and corresponding conserved charges respectively. The net-zero strangeness and total electric charge constrain the components $`\mu ^\alpha `$. The particle numbers are accordingly

$$N_i^{\mathrm{prim}}=V(2J_i+1)\frac{d^3p}{(2\pi )^3}𝑑m_i\left[\gamma _s^{S_i}\text{e}^{\frac{E_i\mu ^\alpha Q_i^\alpha }{T}}\pm 1\right]^1\text{BW}(m_i),$$
(1)
where we include phenomenologically a strangeness saturation factor $`\gamma _s`$ (with $`S_i`$ the total number of strange quarks in hadron species $`i`$) to account for incomplete equilibration in this sector, $`E_i=\sqrt{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+m_i^2}`$; BW is the Breit-Wigner distribution (to be replaced by a $`\delta `$ function for stable hadrons). The final particle numbers are $`N_i=N_i^{\mathrm{prim}}+_j\text{Br}^{ji}N_j^{\mathrm{prim}}`$ due to decays of unstable particles with branching ratios $`\text{Br}^{ji}`$. Such a description can be justified for multiplicities measured over the whole phase-space, since many dynamical effects cancel out in ratios of hadron yields.

We have analyzed two data sets: (i) NA49 $`4\pi `$ multiplicities of $`\pi =\frac{1}{2}(\pi ^++\pi ^{})`$, $`K^\pm `$, $`\overline{p}`$, $`\varphi `$, and $`N_{\mathrm{part}}`$ (taken as the sum over all baryons) in 6 centrality bins in the reaction Pb(158 AGeV) + Pb  (it should be emphasized that protons are not included in our analysis  as they are not participants in non-central collisions), and (ii) PHENIX mid-rapidity densities of $`\pi ^\pm `$, $`K^\pm `$, and $`p^\pm `$ in the reaction Au + Au at $`\sqrt{s}=130`$ AGeV in 5 centrality bins . Results of our fits are displayed in fig. 1. Note that our model tends to underestimate the NA49 $`\varphi `$ yields for peripheral collisions.

A comparison of the individual thermal parameters of both data sets is displayed in fig. 2. Most remarkable is the drop of the baryo-chemical potential $`\mu _B`$ and the rise of the strangeness saturation factor $`\gamma _s`$ when going from $`\sqrt{s}=17`$ AGeV to 130 AGeV. The parameter $`\mu _B`$ is fairly independent of the centrality. $`T`$ seems to stay constant for $`\sqrt{s}=17`$ AGeV, while at 130 AGeV it rises with $`N_{\mathrm{part}}`$. The strangeness saturation factor and, of course, the volume-equivalent system size increase significantly with centrality. (Note the fiducial meaning of the volume for the mid-rapidity data of PHENIX.) Despite the rather limited set of analyzed hadron species, the extracted thermal parameters describe other hadron yields, which are at our disposal in central collisions, fairly well.

## 3 Hadron Spectra

For the reaction Pb(158 AGeV) + Pb, transverse momentum spectra at mid-rapidity are available also in 6 centrality bins . However, the raw data analysis in  needs refinements and cross checks222We are grateful to R. Stock, P. Seyboth, F. Sikler, and P. Jacobs for explaining us the status of the data analyses.. We, therefore, mention here only that our preliminary analysis along the lines in  seems to point to a larger transverse flow with increasing centrality.

## 4 Intermediate-Mass Dileptons

As pointed out in , the dilepton spectra can be analyzed in the same spirit as the hadron multiplicities and hadron momentum spectra, i.e., one discards any detail of the dynamics and asks only for a simple parameterization. As a result, one gets for the thermal dilepton spectrum

$$\frac{dN}{d^4Q}=\frac{5\alpha ^2}{36\pi ^4}N_{\mathrm{dil}}\mathrm{exp}\left\{\frac{M_{}\mathrm{cosh}(YY_{\mathrm{cms}})}{T_{\mathrm{dil}}}\right\},$$
(2)
where $`Q`$ is the lepton pair’s four-momentum, $`M_{}`$ its transverse mass and $`Y`$ its rapidity; $`Y_{\mathrm{cms}}`$ denotes the fire ball rapidity, and $`N_{\mathrm{dil}}`$ is a normalization factor characterizing the space-time volume of the fire ball; flow effects are negligible for invariant mass spectra. In  we have shown that the space-time averaged temperature parameter $`T_{\mathrm{dil}}170`$ MeV (i.e., a value coinciding with the chemical freeze-out temperature) provides a common description of the low-mass CERES data and the intermediate-mass NA50 data at $`\sqrt{s}=17`$ AGeV; also the WA98 photon data are consistent with this value. It was, therefore, a surprise to us that our analysis of the efficiency corrected and centrality binned NA50 data  separately give as optimum fit a temperature scale in the order of 250 MeV,333$`T_{\mathrm{dil}}`$ in bins 1 and 2 (very peripheral collisions) deviates from this value. The extracted values are not safe and depend sensitively on the normalization of the Drell-Yan yield. see figs. 3 and 4. At this temperature the spectral shapes of thermal dileptons and dileptons from correlated decays of open charm mesons are nearly identical within the NA50 acceptance.

We have also checked the temperature sensitivity of the parameterization eq. (2) by analyzing the new low-mass dilepton CERES data  for the reaction Pb(40 AGeV) + Au and find, when using as important ingredient the hadronic cocktail, $`T_{\mathrm{dil}}=145`$ MeV (which also coincides with the chemical freeze-out temperature deduced for this beam energy), see fig. 5.

## 5 Summary

The analyses of hadron multiplicities, hadron momentum spectra and intermediate-mass dilepton spectra point to a centrality independence of the chemical freeze-out temperature and baryo-chemical potential for $`\sqrt{s}=17`$ AGeV, while the strangeness saturation increases with centrality for both $`\sqrt{s}=17`$ AGeV and 130 AGeV. The very preliminary data on transverse momentum spectra of hadrons indicate a stronger flow in central events at $`\sqrt{s}=17`$ AGeV. The NA50 data indicate that the space-time averaged temperature determining the shape of dilepton spectra stays fairly independent of centrality.


text:
# Measurement of single electrons and implications for charm production in Au+Au collisions at √𝑠__{𝑁𝑁} = 130 GeV

## Abstract

Transverse momentum spectra of electrons from Au+Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV have been measured at midrapidity by the PHENIX experiment at RHIC. The spectra show an excess above the background from photon conversions and light hadron decays. The electron signal is consistent with that expected from semi-leptonic decays of charm. The yield of the electron signal $`dN_e/dy`$ for $`p_T>0.8`$ GeV/c is $`0.025\pm 0.004`$(stat.)$`\pm 0.010`$(sys.) in central collisions, and the corresponding charm cross section is $`380\pm 60`$(stat.)$`\pm 200`$(sys.) $`\mu `$b per binary nucleon-nucleon collision.

In this letter, we report the first measurement of single electron spectra, $`(e^++e^{})/2`$, in Au+Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV at the Relativistic Heavy Ion Collider (RHIC). The measurement of single leptons at high transverse momentum ($`p_T1`$ GeV/c) is a useful way to study heavy-quark production, an important probe of hot and dense matter created in high energy heavy ion collisions. Charm production is sensitive to the initial state gluon density. Nuclear and medium effects, such as shadowing and charm quark energy loss, can be studied by comparison of charm production in $`AA`$, $`pA`$, and $`pp`$ collisions. Measurement of charm is important for understanding J/$`\psi `$ suppression (a proposed signal of the deconfinement phase transition) and the dilepton mass distribution in $`1<M_{l^+l^{}}<3`$ GeV, where lepton pairs from charm make significant contributions. In $`pp`$ collisions at the ISR ($`\sqrt{s}=3063`$ GeV), production of single electrons was observed ($`e/\pi 10^4`$) for $`p_T>`$ 1 GeV/c, and interpreted as evidence of open charm production. In $`pp`$ collisions at RHIC energies, the signal level is expected to be higher, since charm production increases with $`\sqrt{s_{_{NN}}}`$ faster than pion production. We recently observed suppression of high $`p_T`$ pion production in Au+Au collisions at RHIC relative to binary nucleon-nucleon ($`NN`$) collision scaling. If charm production scales with $`NN`$ collisions, as expected in the absence of nuclear effects, the $`e/\pi `$ ratio will be even higher in Au+Au collisions at RHIC.

Data used for this analysis were recorded by the PHENIX west-arm spectrometer ($`\mathrm{\Delta }\varphi =90^o`$ in azimuth, $`|\eta |<0.35`$ in pseudo-rapidity), which consisted of a drift chamber (DC), a layer of pad chambers (PC1), a ring imaging Cerenkov detector (RICH), and a lead-scintillator electromagnetic calorimeter (EMCAL). The trigger was provided by beam-beam counters (BBC) and zero-degree calorimeters (ZDC). ZDC and BBC signals were combined to select centrality: central (0-10%), peripheral (60-80%), and minimum bias (0-92%).

The analysis uses 1.23 M minimum bias events with vertex position $`|z|<30`$ cm. Charged particle tracks are reconstructed by the DC and the PC1 with a momentum resolution $`\delta p/p0.6\%3.6\%p`$ (GeV/c). Tracks are confirmed by a matching hit in the EMCAL, which measures the energy $`E`$ deposited with a resolution of $`8.2\%/\sqrt{E(\mathrm{GeV})}1.9\%`$ for test beam electrons. Electron identification is performed using the RICH and the EMCAL. The RICH is filled with 1 atm CO2 and detects on average 10.8 photo-electrons per electron track, while a pion with $`p<4.7`$ GeV/c produces no signal. It is required that at least 3 RICH hits are associated with the track and that their hit pattern is consistent with that of an electron track. After these cuts, a clear electron signal is observed as a narrow peak centered at $`E/p`$ = 1.0. We select tracks in the peak as electron candidates. The $`E/p`$ cut reduces hadron background, and removes conversion electrons created far from the vertex. A hadron deposits only a fraction of its energy in the EMCAL, and the momentum of an off-vertex conversion electron is reconstructed incorrectly. The remaining background, about 10% of the electron candidates, is caused by accidental association of RICH hits with hadron tracks. The background level is measured statistically by an event mixing method, and is subtracted from the yield.

The electron acceptance ($`7.4`$ % of $`dN/dy`$) and efficiency ($`60`$ %) are determined using a detailed GEANT simulation, which satisfactorily reproduces the detector response. Additionally, a multiplicity dependent efficiency loss due to detector occupancy is evaluated by embedding simulated electrons into real events. This efficiency loss is $`27\pm 4`$% ($`4\pm 2`$%) for central (peripheral) collisions and has no significant $`p_T`$ dependence.

Fig. 1 shows the $`p_T`$ distributions of electrons in PHENIX for central, minimum bias, and peripheral collisions. Errors in the figure are statistical. The overall systematic uncertainty, which is the quadratic sum of several few percent effects, is about 11%. Expected sources of electrons are (1) Dalitz decays of $`\pi ^0`$, $`\eta `$, $`\eta ^{}`$, $`\omega `$, and $`\varphi `$, (2) di-electron decays of $`\rho `$, $`\omega `$, and $`\varphi `$, (3) photon conversions, (4) kaon decays ($`K^{0,\pm }\pi e\nu `$), (5) semi-leptonic decay of charm, and (6) other contributions such as bottom decays and thermal di-leptons. In this analysis, sources (1)-(4) are considered to be background.

We have calculated the contributions from Dalitz and di-electron decays with a hadron decay generator. PHENIX has measured the $`p_T`$ distributions of $`\pi ^\pm `$ in $`0.2<p_T<2.2`$ GeV/c and of $`\pi ^0`$ in $`1<p_T<4`$ GeV/c. Since the $`\pi ^\pm `$ and $`\pi ^0`$ data are consistent in the overlapping region, we fit a power law function to the combined data sets to determine the input $`\pi ^0`$ spectrum for the decay generator. The $`p_T`$ distribution of any other hadron $`h`$ is obtained from the $`\pi ^0`$ spectrum by replacing $`p_T`$ with $`\sqrt{p_T^2+m_h^2m_{\pi ^0}^2}`$. The shapes of the resulting $`p_T`$ spectra of $`K^\pm `$, $`p`$, and $`\overline{p}`$ agree with the PHENIX measurements within 20%. In this parameterization $`h/\pi ^0`$ ratios approach constants at high $`p_T`$. We assume the following asymptotic ratios to fix the relative normalizations: $`\eta /\pi ^0=0.55`$, $`\eta ^{}/\pi ^0=0.25`$, $`\rho /\pi ^0=\omega /\pi ^0=1.0`$, $`\varphi /\pi ^0=0.40`$. Except for the $`\varphi `$, these ratios are taken from proton beam data of SPS, FNAL, and ISR experiments. The $`\eta /\pi ^0`$ ratio is consistent with a measurement in Pb+Pb collisions at SPS. The $`\varphi /\pi ^0`$ ratio is based on the integrated ratio $`\varphi /h^{}0.02`$ in Au+Au collisions at $`\sqrt{s_{_{NN}}}`$ = 130 GeV. We assign to each ratio a conservative systematic uncertainty of 50%.

Photon conversions are evaluated using a combination of the GEANT simulation and the hadron decay generator. Since $`p_T`$ spectra of externally converted electrons are similar to those from Dalitz decay, the conversion spectra can be approximated by scaling the Dalitz decay spectra by an experiment specific factor, $`R_{conv}`$ = Conversion/Dalitz. $`R_{conv}`$ is evaluated using the GEANT simulation and is cross-checked by comparing the relative yield of reconstructed Dalitz and conversion pairs in the simulation and in data. The simulation shows that $`R_{conv}`$ has only a weak $`p_T`$ dependence, primarily due to the energy dependence of the pair creation cross section. $`R_{conv}`$ is parameterized as $`(1.9\pm 0.2)\times (10.0718\times p_T^{0.76})`$.

Background from kaon decays is also evaluated using the GEANT simulation, and is found to be negligible.

The upper panel of Fig. 2 shows the ratio of the measured electrons to the calculated background versus $`p_T`$ for minimum bias events. The shaded region is the quadratic sum of systematic errors in the electron measurement and in the background. The latter includes uncertainties in the normalization and the shape of the $`\pi ^0`$ spectrum, in the $`h/\pi ^0`$ ratios, and in $`R_{conv}`$. A significant electron excess above the background is observed for $`p_T>0.6`$ GeV/c. Central collisions show a similar excess. The peripheral collision data sample lacks sufficient statistics to reveal a signal in this analysis.

Fractional contributions to the background are shown in the lower panel of Fig. 2. More than 80% of the background is from $`\pi ^0`$ decay, directly from the Dalitz decay or indirectly from photon conversion. The $`\pi ^0`$ spectrum is well constrained by the PHENIX measurement. The next most important background source is $`\eta `$ decay. Given the assigned systematic error, the upper limit of the high $`p_T`$ asymptotic $`\eta /\pi ^0`$ ratio is 0.83. Since this ratio, corrected for feed-down, would imply that the primary $`\eta /\pi ^01`$, this provides a conservative limit on contributions from $`\eta `$’s. Contributions from all other hadrons combined are only a few percent of the total.

Background-subtracted electron spectra are shown in Fig. 3. The error bars on the data points represent the statistical errors, while the systematic error due to the background subtraction is indicated by brackets. The integrated yield of the electron signal $`dN_e/dy`$ for $`p_T>0.8`$ GeV/c is $`0.025\pm 0.004`$(stat.)$`\pm 0.010`$(sys.) for central collisions and is $`0.0079\pm 0.0006`$(stat.)$`\pm 0.0034`$(sys.) for minimum bias collisions.

Semi-leptonic decay of charmed hadrons is an expected source of the electron signal. We use the event generator PYTHIA to estimate electron spectra from charm decay. We tuned the parameters of PYTHIA such that charm production data at SPS and FNAL and single electron data at the ISR are well reproduced. The charm production cross section in $`pp`$ collisions from this PYTHIA calculation is $`\sigma _{c\overline{c}}=330\mu `$b at $`\sqrt{s}=`$130 GeV. The electron spectrum in Au+Au collisions is then calculated as $`EdN_e/dp^3=T_{AA}\times Ed\sigma _e/dp^3`$, where $`Ed\sigma _e/dp^3`$ is the electron spectrum from charm decay calculated with PYTHIA, and $`T_{AA}`$ (listed in Table I) is the nuclear overlap integral calculated from a Glauber model. The calculated electron spectra shown in Fig. 3 are in reasonable agreement with the data.

Before attributing the entire electron signal to open charm decays, it is necessary to quantify contributions from other possible sources. An analogous PYTHIA estimate of the bottom decay contribution is shown in Fig. 3. It becomes significant only above the measured $`p_T`$ range. Expected contributions from $`J/\mathrm{\Psi }`$ and Drell-Yan are negligible. In Pb+Pb collisions at SPS, direct photons and an enhanced yield of low mass di-leptons have been reported. If these are due to thermal radiation from hot matter, an even larger production is expected at RHIC energies, and can contribute to the electron signal. Since $`\rho e^+e^{}`$ contributes less than 1% to the calculated background as shown in Fig. 2, and since the dominant source of thermal di-leptons is $`\pi +\pi \rho e^+e^{}`$, a significant contribution from thermal di-leptons is unlikely. There are several predictions for direct photons at RHIC energies. The conversion electron spectrum calculated from a prediction in ref. is shown in Fig. 3 for central collisions. It could explain 10-20% of the signal, with large theoretical uncertainties.

Neglecting these other possible sources and assuming that all the electron signal is from charm, we derive the charm cross section corresponding to the electron data. We fit the charm electron spectrum from PYTHIA to the data for $`p_T>0.8`$ GeV/c, and obtain the rapidity density $`dN_{c\overline{c}}/dy|_{y=0}`$ and the total yield $`N_{c\overline{c}}`$ of open charm. They are then converted to cross sections per $`NN`$ collision: $`d\sigma _{c\overline{c}}/dy=(dN_{c\overline{c}}/dy)/T_{AA}`$ and $`\sigma _{c\overline{c}}=N_{c\overline{c}}/T_{AA}`$. Results are shown in Table I. The systematic error is a quadratic sum of many sources. For central collisions, they are background subtraction ($`\pm 44\%`$), uncertainties in the PYTHIA calculation ($`\pm 11\%`$ from $`<k_T>=1.5\pm 0.5`$, $`\pm 13\%`$ from $`D^+/D^0=0.65\pm 0.35`$, $`\pm 8\%`$ from PDFs), fit range ($`\pm 18\%`$), and $`T_{AA}`$ ($`\pm 7\%`$). Note that any finite contribution from neglected sources would reduce the derived charm cross section. Without nuclear or medium effects in charm production, $`\sigma _{c\overline{c}}`$ per $`NN`$ collision should be independent of centrality. Within uncertainties, our data are consistent with this expectation, in possible contrast to the attribution of increased charm production as the source of enhanced di-muon production reported in Pb+Pb collisions at SPS.

The single electron signal yield (divided by $`T_{AA}`$ to give the cross section per $`NN`$ collision) and the derived charm cross section are compared with single electron data of ISR experiments and charm data of fixed target experiments in Fig. 4. Cross section curves calculated with PYTHIA, which has been tuned to the charm data and the ISR electron data, and a charm cross section curve from a next-to-leading order (NLO) pQCD calculation are also shown in the figure. Our data are consistent with both of the calculations within large uncertainties.

In conclusion, we have observed single electrons above the expected background from decays of light hadrons and photon conversion in Au+Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV. The observed signal is consistent with semi-leptonic decay of charm. The forthcoming high statistics Au+Au data and $`pp`$ comparison data at full RHIC energy ($`\sqrt{s_{_{NN}}}=200`$ GeV) will be useful to clarify the nature of the single electron signal and to better determine heavy-quark production in Au+Au collisions at RHIC.

We thank the staff of the Collider-Accelerator and Physics Departments at BNL for their vital contributions. We acknowledge support from the Department of Energy and NSF (U.S.A.), MEXT and JSPS (Japan), RAS, RMAE, and RMS (Russia), BMBF, DAAD, and AvH (Germany), VR and KAW (Sweden), MIST and NSERC (Canada), CNPq and FAPESP (Brazil), IN2P3/CNRS (France), DAE and DST (India), KRF and CHEP (Korea), the U.S. CRDF for the FSU, and the US-Israel BSF.