text:
# The chemical equilibration volume: measuring the degree of thermalization

## I INTRODUCTION

Experiments are presently underway at the relativistic heavy-ion collider (RHIC) at BNL and the SPS at CERN to collide large nuclei at ultra-relativistic energies. The aim is to create energy densities high enough for the production of a state of deconfined quarks and gluons: the quark gluon plasma (QGP). Inherent in this statement is the assumption that, in such a collision, we have created a system close to kinetic and possibly chemical equilibrium. The QGP is rather short lived and soon hadronizes, leading one to observe a plethora of mesons and baryons at the detectors. One has to surmise by indirect means as to whether a QGP or at least a thermalized system was created in the history of a given collision. On confronting experimental results, two possibilities of what may have occurred in such collisions emerge.

Even though most energetic partons are expected to undergo small angle (soft) scattering, it is possible that a sizable number of partons undergo large angle (hard) scattering followed by multiple rescattering in the central region. This may lead to the appearance of sufficient energy densities in the central region of such collisions that quarks and gluons that are usually bound in nuclei become deconfined over a region much larger than the size of a nucleon. One obtains a thermalized Quark Gluon Plasma (QGP) col75 ; shu80 . Such a state will undoubtably be short lived and is expected to quickly expand and hadronize into a hot plasma of mesons and baryons. The system will continue to expand and cool, finally interactions freeze out and the hadrons free stream to the detectors.

The above paragraph represents a rather optimistic possibility. The other extreme is to view the above collision as a superposition of a series of nucleon nucleon collisions. Each such collision produces jets of hard hadrons and some soft hadrons as is the case in a proton-proton (or proton-antiproton) collision. The hadrons mostly emerge without interaction from the central region. There is no multiple partonic rescattering, and as a result no widespread thermalization entails and no quark gluon plasma is formed. Another equivalent statement, would be that that the quarks and gluons that eventually lead to the formation of the hadrons are not deconfined over a region much larger than the volume of a nucleon.

It is reasonable to assume that actual heavy-ion experiments lie somewhere in between these two extremes. To obtain complete thermalization one needs an infinite amount of time, not available in a heavy-ion collision. On the other hand, it is rather unlikely that every hadron produced in each of the nucleon nucleon collisions will emerge untouched from a region with a density many times that of normal nuclear matter. Undoubtedly, there is some thermalization: we picture a scenario where an impinging nucleon has a probability distribution to undergo $`n`$ hard or semi-hard scatterings; each of these $`n`$ points becomes a source of particles (see e.g., Ref. hwa90 ). There is also a probability distribution of the number of particles produced at these $`n`$ points, greater particle production leads to a greater probability of final state multiple scattering. Multiple scattering leads to the onset of kinetic followed by chemical equilibrium. As a result, one obtains multiple small domains in the larger volume, each, very close to being in chemical and kinetic equilibrium. The size of the chemically equilibrated volume is assumed to be smaller than the kinetically equilibrated volume. The particles in these domains can undergo further rescattering with particles in surrounding domains leading to an increase in the domain size. If very little rescattering occurs then the size of these domains is closer to the size of a nucleon. Extensive rescattering leads to domain sizes of the order of the size of the system. We will concentrate on the size of a chemically equilibrated domain. Whether the size of such domains is of the order of the system size or the size of a nucleon will be the primary question addressed in this article.

The primary observables in such reactions are the populations of the various species of hadrons bra99 ; bra01 ; raf01 ; zsc02  and their momentum distributions bro01 ; bro02 ; fri03 . Interestingly, the broad characteristics of these observables fail to distinguish between the above two scenarios koc02 . The reason behind this strife is the observation that the abundances of various hadrons produced in proton-proton, or even $`e^+e^{}`$ collisions may be well reproduced by thermal models bec97 . To distinguish between these two possibilities one needs to look at finer features of the data (see e.g.,jeo00 ; ble00 ). One such feature is to look at the momentum spectra of the various hadrons, especially at high transverse momenta adl02 . Yet another feature is the abundances of hadrons which carry a conserved quantum number ada02 . In this article we will be concerned with the latter aspect.

There are many exactly conserved and almost exactly conserved quantum numbers in a heavy-ion collision. For the moment we will specifically concentrate on the mid-rapidity or central region (the full size of the region must be such that at freeze-out, the boundaries of the system are not beyond causal contact). We will discuss the effects on 4$`\pi `$ observables later. Net baryon number is exactly conserved. As the fireball lives for a time on the order of the strong interaction scale, weak interactions and hence flavor changing currents are considered absent. Thus net strangeness, and charm are almost exactly conserved. The application of thermal models to particle production involving conserved quantum numbers requires the use of the canonical ensemble for these numbers. This was first pointed out in proton-proton collisions in Ref. hag68 .

If the number of particles are large (or we are dealing with a rather large system), one approaches the grand canonical limit, where one may simply introduce a Lagrange multiplier for the conserved quantity and evaluate the much simpler grand canonical partition function. One may then fine tune the Lagrange multiplier to obtain the correct value of the conserved number on the average. If the number of particles are small, (or we are dealing with a smaller system), then canonical methods yield results in marked distinction from grand canonical estimates cle91 ; ham00 .

Using such models one may obtain a better fit to the data on particle abundances, both for those influenced by conservation laws and those free of such constraints. The constraint that we will be concerned with will be strangeness conservation. As the partons(hadrons) rescatter they will produce strange partons(hadrons). This is the onset of chemical equilibration. By the time freeze-out occurs, chemical equilibrium has taken place over a volume possibly smaller than the system size, this will be referred to as the domain or partition volume. The radius of such a domain will, in principle, be governed by the strangeness diffusion constant. We are assuming the full system to be made up of multiple such domains. Within each domain we will enforce exact strangeness conservation by calculating strange particle yields using the canonical ensemble. Particles not subject to this constraint will have their yields calculated via a grand canonical ensemble with a baryon chemical potential $`\mu `$. The chemical potential will be assumed to be a constant throughout the system; in this effort we will set it to zero. It should be pointed out that in a grand canonical calculation, the number of particles is strictly proportional to the volume. This implies that the yields of particles not constrained by strangeness is unaffected by the partitioning of the system into domains. The pion yield, for example, is directly proportional to the total freeze-out volume. In the following we will propose measurable observables that may be used to discern the domain volume (or the strangeness chemical equilibration volume).

Once again, we point out that a small domain volume does not necessarily indicate a small kinetically equilibrated volume. The entire system may have reached kinetic equilibrium and frozen out before having reached strangeness chemical equilibrium. The strange quark has a current mass more than ten times that of the up or down quark, while the kaons are about four times as heavy as the pions. Hence, if it turned out that the domain volume is of the order of the size of the system, then it would point towards the presence of an efficient strangeness equilibration mechanism in the history of the collision. The leading candidate for such a mechanism is the quark gluon plasma.

In the following, we begin by outlining the calculation of the partition function for the strange sector for a single domain in Sect. II. Here we will introduce the recursion relation technique which greatly simplifies the numerical calculation compared to the previous methods involving Bessel function techniques. In Sect. III we compute the probability to obtain $`A`$ distinguishable particles from the full freeze-out volume. Here a second recursion relation is introduced. With the calculation of the probability distributions, various statistical quantities may be computed as a function of the number of domains. This will be performed and the results discussed in Sect. IV. The reader not interested in the calculational details may skip directly to this section. Here, we will propose measurable observables that will allow one to narrow down the range of domain volumes. Our conclusions are summarized in Sect. V.

## II The partition function of a domain.

Let us assume that freeze-out occurs in a large volume of size $`V_f`$. If the system is not chemically equilibrated over this volume, but over a much smaller volume $`V`$, we divide up the system into $`p`$ partitions or domains, each with a volume $`V=V_f/p`$. In each domain we assume, we have full chemical and thermal equilibrium. Strangeness is exactly conserved in a domain, hence the partition function for net strangeness equal to $`N`$ is:

$`Z(N)`$ $`=`$ $`Z_0{\displaystyle \underset{n_1=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_1^{n_1}}{n_1!}}{\displaystyle \underset{n_1=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_1^{n_1}}{n_1!}}{\displaystyle \underset{n_2=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_2^{n_2}}{n_2!}}{\displaystyle \underset{n_2=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_2^{n_2}}{n_2!}}{\displaystyle \underset{n_3=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_3^{n_3}}{n_3!}}{\displaystyle \underset{n_3=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{S_3^{n_3}}{n_3!}}`$ (1)
$`\times \delta _{n_1n_1+2n_22n_2+3n_33n_3,N}.`$
In the above equation, $`Z_0`$ is the grand canonical partition function for non-strange particles:

$`Z_0={\displaystyle \underset{i}{}}\mathrm{exp}\left({\displaystyle \frac{d^3r_id^3p_i}{(2\pi \mathrm{})^3}e^{\beta (E_i\mu _i)}}\right),`$ (2)
where, $`i`$ runs over all the different kinds of non-strange particles and anti-particles. The $`\mu _i`$’s represent the respective chemical potentials. In Eq. (1) The $`S_k`$’s represent the single particle partition functions for particles with strangeness $`k`$. For instance for strangeness = -1,

$`S_1=S(K^{})+S(\overline{K}^0)+S(K^{})+S(K_1)+S(\mathrm{\Lambda })+S(\mathrm{\Sigma })+\mathrm{}`$ (3)
In Eq. (1), the Kroneker delta at the end ensures that the total strangeness is held fixed at $`N`$ (an integer value). The $`n!`$ terms in the denominators serve the dual purpose of symmetrization and the removal of double counting. In this paper we will deal with temperatures where quantum statistics will not play a role; in such a setting the $`n!`$ terms are introduced to avoid Gibbs’ paradox rei85

On calculating the partition function, one may evaluate its various derivatives to obtain various quantities of interest. The usual method used up to now is to take the Fourier transform of the Kroneker delta function, which allows one to express the strange partition function $`Z(N)/Z_0`$ in a closed form cle91 :

$`Z_S(N)`$ $`=`$ $`{\displaystyle \frac{Z(N)}{Z_0}}={\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }𝑑\varphi \mathrm{exp}\left(S_0+S_1e^{i\varphi }+S_1e^{i\varphi }+S_2e^{2i\varphi }+S_2e^{2i\varphi }+S_3e^{3i\varphi }S_3e^{3i\varphi }\right).`$ (4)
The solution of this equation however involves the use of series of sums of modified Bessel functions (see Eq. (24) of Ref. cle91 ). This turns out to be rather complicated to handle, both in analytic form and in numerical computation. If one is not interested in an analytic expression but rather in the numerical results, a simple and elegant means of calculating in the canonical ensemble exists. This involves the use of recursion relation techniques first introduced by Mekjian for intermediate energy heavy-ion reactions cha95 ; das98 ; das99 . The use of such techniques to solve the problem of exact strangeness conservation will be introduced in the following subsection. The method of recursion relations is very general and may be used to calculate a variety of problems with exact conservation conditions e.g., baryon number, charm, charge. The recursion relations can also be used to calculate partition functions with multiple conserved charges maj03 . A variant of this method will be used later to calculate strangeness production from all the domains together i.e., for the full system.

### II.1 The Recursion Relation Technique

We begin by expanding Eq. (1) for the simple and realistic case of $`N=0`$. Once again, we are using strangeness to illustrate the derivation of the recursion relations, however any conserved charge may be used. Making explicit use of the Kroneker delta function, we obtain,

$`Z_S={\displaystyle \frac{Z}{Z_0}}`$ $`=`$ $`1+S_1S_1+{\displaystyle \frac{S_1^2}{2!}}{\displaystyle \frac{S_1^2}{2!}}+{\displaystyle \frac{S_1^2}{2!}}S_2+S_2{\displaystyle \frac{S_1^2}{2!}}+\mathrm{}`$ (5)
$`=`$ $`1+S_1S_1+\left({\displaystyle \frac{S_1^2}{2!}}+S_2\right)\left({\displaystyle \frac{S_1^2}{2!}}+S_2\right)+\mathrm{}`$
Where we have rewritten $`Z(0)`$ as simply $`Z`$. In the above series, we identify the various terms in the brackets as,

$`Z_S={\displaystyle \frac{Z}{Z_0}}`$ $`=`$ $`1+Z_1Z_1+Z_2Z_2+\mathrm{}+Z_nZ_n+\mathrm{}`$ (6)
where $`Z_n`$ is the partition function of a box containing only strange particles, with net strangeness of the box equal to $`n`$. The term $`Z_n`$ represents a similar system with only antistrange particles with net strangeness $`n`$. One immediately notes that each term in the sum has net strangeness zero. The general term $`Z_n`$ may be expressed formally as

$`Z_n={\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=n}{}}{\displaystyle \underset{k}{}}{\displaystyle \frac{S_k^{n_k}}{n_k!}}.`$ (7)
We will eventually set $`S_k`$ to 0 for all $`k`$ greater than 3 i.e., there are no baryons with strangeness greater than 3; presently we keep them in for generality. The constraint on the numbers $`n_k`$, as the subscript on the sum indicates, is such that the total strangeness sums to $`n`$. This implies, that, with these value of $`n_k`$

$`{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}kn_k=1.`$ (8)
We now insert unity, as obtained in the above equation, into Eq. (7), to obtain,

$`Z_n`$ $`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=n}{}}{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}kn_k\left[{\displaystyle \underset{j}{}}{\displaystyle \frac{S_j^{n_j}}{n_j!}}\right]`$ (9)
$`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=n}{}}{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}kn_k{\displaystyle \frac{S_1^{n_1}}{n_1!}}\mathrm{}{\displaystyle \frac{S_k^{n_k}}{n_k!}}\mathrm{},`$
where we have expanded out the product. We now cancel the appearance of $`n_k`$ in numerator and denominator, and may reverse the order of the sums to obtain,

$`Z_n`$ $`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=n}{}}{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}k{\displaystyle \frac{S_1^{n_1}}{n_1!}}\mathrm{}{\displaystyle \frac{S_k^{n_k}}{(n_k1)!}}\mathrm{}`$ (10)
$`=`$ $`{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}S_k{\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=nk}{}}{\displaystyle \frac{S_1^{n_1}}{n_1!}}\mathrm{}{\displaystyle \frac{S_k^{n_k}}{(n_k)!}}`$
$`=`$ $`{\displaystyle \frac{1}{n}}{\displaystyle \underset{k}{}}S_kZ_{nk}.`$
Where in the last two lines, we have simply redefined $`n_k=n_k1`$, which changes the sum to $`kn_k=nk`$. Since $`k3`$, for hadrons, we obtain the simplified recursion relation:

$`Z_n`$ $`=`$ $`{\displaystyle \frac{1}{n}}\left[S_1Z_{n1}+2S_2Z_{n2}+3S_3Z_{n3}\right].`$ (11)
The actual method of computation involves simply calculating $`S_1,S_2,S_3`$, and using this to calculate $`Z_1,Z_2,Z_3`$ and so on in ascending order. A similar method may be used to compute the $`Z_n`$’s. The presence of the $`n!`$ in the denominator of $`Z_n`$ (see Eq. (1)) greatly reduces the magnitude of $`Z_n`$ for large $`n`$. As a result, one never needs to go beyond $`n=200`$ to obtain an accuracy of 12 significant digits. The use of the recursion relations, reduce such computations to a fraction of a second on a regular PC. Having calculated the partition function we may easily compute the strange free energy $`F_s`$, where

$`F=T\mathrm{log}(Z)=T\mathrm{log}\left(Z_0Z_S\right)=F_0+F_s.`$ (12)
This is plotted in Fig. (1) as a function of domain size and temperature.

### II.2 Moments of the distribution.

The recursion relations introduced in the previous subsection may also be suitably modified to calculate the moments of the distribution. Formally these correspond to various derivatives of the partition function for different numbers of particles and may be calculated numerically as such. However, they may be exactly evaluated without the use of any numerical differentiation techniques. In the following we simply illustrate the calculation of the first and second moments. The higher moments may be calculated in similar fashion. We will, only demonstrate the method for strange hadrons; the calculation for anti-strange hadrons is almost identical.

We focus on the general term in the series in Eq. (6). Expanding the strange part of this term, one obtains,

$`Z_0Z_nZ_n=Z_0\left[{\displaystyle \underset{n_1,n_2,\mathrm{}}{}}{\displaystyle \underset{k}{}}{\displaystyle \frac{S_k^{n_k}}{n_k!}}\right]Z_n.`$ (13)
Each term in the sum represents the partition function of the system with the restriction that there be $`n_1`$ particles with strangeness 1, $`n_2`$ particles with strangeness 2 and so on. Dividing this by the full partition function will give us the probability to obtain this partition (we represent the set of numbers $`n_1,n_2,\mathrm{}`$ by the vector $`\stackrel{}{n}`$),

$`P_\stackrel{}{n}={\displaystyle \frac{_k\frac{S_k^{n_k}}{n_k!}|_\stackrel{}{n}Z_n}{_nZ_nZ_n}}.`$ (14)
The probability to obtain a particular partition, represents a central quantity in our computations. One may use this to compute any statistical quantity, e.g., the mean number of particles of strangeness $`k`$ will thus be given by the general expression,

$`n_k=n_kP_{\stackrel{}{n}_1}+n_kP_{\stackrel{}{n}_2}+\mathrm{}`$ (15)
On the $`r.h.s.`$ of the above equation, $`n_k`$ is an operator acting on the partitions (or vectors) $`\stackrel{}{n}_1,\stackrel{}{n}_2`$ etc. The eigenvalue of this operator is the number of particles with strangeness $`k`$ in the vector $`\stackrel{}{n}_1,\stackrel{}{n}_2`$ etc. Substituting the value of $`P_\stackrel{}{n}`$, we obtain,

$`n_k`$ $`=`$ $`{\displaystyle \frac{_n_\stackrel{}{n}n_k_i\frac{S_i^{n_i}}{n_i!}|_{_iin_i=n}Z_n}{_nZ_nZ_n}}`$ (16)
$`=`$ $`{\displaystyle \frac{_n_\stackrel{}{n}|_{_{i=1}^{i_{max}}in_i=n}\frac{S_1^{n_1}}{n_1!}\mathrm{}\frac{S_k^{n_k}}{(n_k1)!}\mathrm{}\frac{S_{i_{max}}^{n_{i_{max}}}}{n_{i_{max}}!}Z_n}{_nZ_nZ_n}}`$
$`=`$ $`{\displaystyle \frac{_n_\stackrel{}{n}S_k_i\frac{S_i^{n_i}}{n_i!}|_{_iin_i=nk}Z_n}{_nZ_nZ_n}}.`$
Where, in the top line, the first sum represents the sum over different values of total strangeness, and the second sum represents the different partitions (or vectors) $`\stackrel{}{n}`$ that will lead to the same total strangeness $`n`$. In the second line we have, for illustration, expanded the product $`_i`$; here $`i_{max}`$ is equal to 3. In the last line we have, once again, shifted $`n_kn_k1`$, and hence the sum $`_iin_i`$ is equal to $`nk`$. Extracting $`S_k`$ from the sum, we note that the sum over all vectors $`\stackrel{}{n}`$, under the constraint $`_iin_i=nk`$, is simply the pure strange partition function with strangeness $`nk`$, i.e.,

$`n_k`$ $`=`$ $`{\displaystyle \frac{_{n=k}^{\mathrm{}}S_kZ_{nk}Z_n}{_{n=0}^{\mathrm{}}Z_nZ_n}}.`$ (17)
Using the same method as before to calculate the $`Z_n`$’s, we may simply calculate the mean number of particles with strangeness $`k`$. The ratio of the means for different strangeness to the number of pions is plotted in Fig. (2).

The second moment may also be easily evaluated. The natural quantity to evaluate in this recursion relation scenario is the term $`f_2n^2=n_k^2n_k=n_k(n_k1)`$, this is obtained from the general expression,

$`n_k^2n_k=n_k(n_k1)P_{\stackrel{}{n}_1}+n_k(n_k1)P_{\stackrel{}{n}_2}+\mathrm{}`$ (18)
Substituting the expressions for the probabilities, and using methods almost identical to the derivation for the first moment, we obtain

$`n_k(n_k1)={\displaystyle \frac{_{n=2k}^{\mathrm{}}S_k^2Z_{n2k}Z_n}{_{n=0}^{\mathrm{}}Z_nZ_n}}.`$ (19)
The method of evaluation of these moments, involves simply evaluating the various $`Z_n`$’s to the required degree of accuracy and then using those to calculate the sums in the above expressions. We also plot the ratio of the second moment to the square of the first moment in Fig. (2). In all these plots, we have plotted the simplest case of zero chemical potential. It should be pointed out that incorporating a finite chemical potential involves almost no increase in the complexity or running time of the computation.

## III The probability of $`N`$ particles from $`p`$ domains.

In the previous sections we have dealt with the partition function, and strangeness production from a single domain. In our picture, the central region in a heavy-ion collision has a full freeze-out volume $`V_f`$ which is divided into $`p`$ domains of volume $`V`$. In each domain, we have reached full thermal and chemical equilibrium in the canonical sense. As an aside, we point out here, as we did in the introduction that full chemical equilibrium may not have been reached, and the domain may in principle be chemically undersaturated or oversaturated. The usual method of dealing with such a scenario, is to introduce a fugacity factor $`\gamma `$ multiplicatively with the factors $`S_k`$. We leave this complication for a later effort. Yet another complication that will be ignored in this first effort, is a variation in the temperature and baryon chemical potential from domain to domain. All domains will be assumed to have the same temperature and chemical potential. In this section we will calculate the various probability distributions for the production of strange particles from the full system of volume $`V`$. The reason for concentrating on the probability distribution is obvious: it represents the most basic quantity that may be computed in a given situation. All other quantities such as the mean, variance, higher moments etc. may be computed using the probability distributions.

Of the variety of strange particles produced in a heavy-ion collision the Kaon is by far the most common and the easiest to identify. We thus compute the probability to obtain $`A`$ Kaons from the full system that consists of $`p`$ domains, i.e., $`P_p^A`$. The method employed will be almost identical for particles with strangeness -2 ($`\mathrm{\Xi }`$ particles) and -3 ($`\mathrm{\Omega }`$ particles). The $`A`$ Kaons are indistinguishable bosons: thus the probability to obtain $`A`$ indistinguishable bosons from $`p`$ domains may be written schematically as

$`P_p^A={\displaystyle \frac{D_p^A}{A!}},`$ (20)
where $`D_p^A`$ is the probability to obtain $`A`$ distinguishable Kaons from $`p`$ partitions. If there are $`n_1`$ Kaons emanating from domain 1, $`n_2`$ from domain 2, etc., then we define the quantity,

$`D_p^A(n_1,n_2\mathrm{}n_p)={\displaystyle \frac{A!}{n_1!n_2!\mathrm{}n_p!}}D_1^{n_1}D_1^{n_2}\mathrm{}D_1^{n_3}.`$ (21)
Using this we may express the probability to obtain $`A`$ distinguishable Kaons from $`p`$ domains as,

$$D_p^A=\underset{n_1,n_2\mathrm{}n_p|_in_i=A}{}D_p^A(n_1,n_2\mathrm{}n_p).$$

As a result we obtain,

$`P_p^A`$ $`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}n_p|_{k=1}^pn_k=A}{}}{\displaystyle \frac{D_1^{n_1}}{n_1!}}{\displaystyle \frac{D_1^{n_2}}{n_2!}}\mathrm{}{\displaystyle \frac{D_1^{n_p}}{n_p!}}`$ (22)
$`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}n_p|_{k=1}^pn_k=A}{}}P_1^{n_1}P_1^{n_2}\mathrm{}P_1^{n_p}.`$
Where $`P_1^{n_k}=\frac{D_1^{n_k}}{n_k!}`$ represents the probability to obtain $`n_k`$ indistinguishable Kaons from 1 domain. In Eq. (1) the $`S_k`$ represent the single particle partition functions for particles with strangeness $`k`$. There are many such particles. The quantity $`S_k`$ is composed of a sum of the single particle partition functions for each of these particles (see Eq. (3)). From Eq. (3) we extract $`S(K)`$, and define the quantity $`\overline{S}_1`$, as

$$S_1=S(K)+S(K^{})+S(K_1)+S(\mathrm{\Lambda })+S(\mathrm{\Sigma })+\mathrm{}=S(K)+\overline{S}_1.$$

With this, we define a partition function with total strangeness $`n`$, consisting of only strange particles but not containing any Kaons, as

$`\overline{Z}_n={\displaystyle \underset{n_1,n_2,\mathrm{}|{\scriptscriptstyle kn_k}=n}{}}{\displaystyle \frac{\overline{S}_1^{n_1}}{n_1!}}{\displaystyle \frac{S_2^{n_2}}{n_2!}}{\displaystyle \frac{S_3^{n_3}}{n_3!}}.`$ (23)
The above partition function, may also be easily evaluated with a similar recursion relation:

$$\overline{Z}_n=\frac{1}{n}\left[\overline{S}_1\overline{Z}_{n1}+2S_2\overline{Z}_{n2}+3S_3\overline{Z}_{n3}\right].$$

In the above $`S_2`$ and $`S_3`$ remain unaffected. Using the above two equations one may construct the probability $`P_1^n`$ of obtaining exactly $`n`$ Kaons from a single domain, with net strangeness zero. Using Eq. (14), this may be written as,

$`P_1^n={\displaystyle \frac{Z_0}{Z}}{\displaystyle \frac{S(K)^n}{n!}}{\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}\left[{\displaystyle \underset{n_1,n_2,n_3|{\scriptscriptstyle in_i}=N}{}}{\displaystyle \frac{\overline{S}_1^{n_1}}{n_1!}}{\displaystyle \frac{S_2^{n_2}}{n_2!}}{\displaystyle \frac{S_3^{n_3}}{n_3!}}\right]Z_{(N+n)}.`$ (24)
In the above equation, the term in square brackets represents the partition function for a box containing only strange particles, with net strangeness $`N`$, none of which are Kaons. Obviously, the term in square brackets is simply $`\overline{Z}_N`$. With this we may write down the full probability to obtain $`A`$ Kaons from the full system with $`p`$ partitions as

$`P_p^A`$ $`=`$ $`{\displaystyle \underset{n_1,n_2,\mathrm{}n_p|_kn_k=A}{}}\left({\displaystyle \frac{Z_0}{Z}}\right)^p{\displaystyle \underset{k=1}{\overset{p}{}}}\left({\displaystyle \frac{S(K)^{n_k}}{n_k!}}{\displaystyle \underset{t_k=0}{\overset{\mathrm{}}{}}}\overline{Z}_{t_k}Z_{(t_k+n_k)}\right)`$ (25)
$`=`$ $`\left({\displaystyle \frac{Z_0}{Z}}\right)^pS(K)^A{\displaystyle \underset{n_1,n_2,\mathrm{}n_p|{\scriptscriptstyle n_k}=A}{}}{\displaystyle \underset{k}{}}w_{n_k}.`$
In the above equation

$`w_{n_k}=\left({\displaystyle \underset{t_k=0}{\overset{\mathrm{}}{}}}\overline{Z}_{t_k}Z_{(t_k+n_k)}\right)/n_k!`$ (26)
Where $`w_{n_k}`$’s for $`0<n_k<A`$ may be easily calculated using the $`\overline{Z}_t`$ and the $`Z_t`$’s. The evaluation of the sum with the constraint in Eq. (25) seems to be extremely complicated. However, yet another recursion relation exists which allows this to be computed swiftly. Note that the sum

$$S_p^A=\underset{n_1,n_2,\mathrm{}n_p|{\scriptscriptstyle n_k}=A}{}\underset{k}{}w_{n_k},$$

requires us to distribute $`A`$ particles in $`p`$ distinguishable boxes. If we put $`q`$ particles in the first box, then we have to distribute $`Aq`$ particles in $`p1`$ boxes. Summing over $`q`$ gives us the recursion relation,

$`S_p^A={\displaystyle \underset{q=0}{\overset{A}{}}}w_qS_{p1}^{Aq}.`$ (27)
One may make the simple observation that $`S_1^k=w_k`$. This along with Eq. (27) allows to calculate the probability $`P_p^A`$ recursively. One essentially starts with $`S_1^k`$, uses these to build the $`S_2^k`$ and so on, until one obtains $`S_p^A`$ with the chosen values of $`p`$ and $`A`$. A plot of these probabilities is presented in Figs. (3,6,7).

## IV Results and discussions

We begin the analysis with the kaon distributions (Fig. (3)). In our calculations, the full freeze-out volume $`V_f`$ is taken to be about 400 times the volume of a nucleon. In the left panel of Fig. (3), various cases are presented where this $`V_f`$ is divided into 1,17,33…81 domains. In comparison, the right panel represents the probability distributions for one partition, but for different cases of a decreasing full freeze-out volume. The volumes range from $`V_f=400V_{nuc}`$ to $`V_f=350V_{nuc}`$. Note that, in the left panel, as the number of domains is increased, the Gaussian like curves slowly shift toward the left. A similar pattern is realised if instead of partitioning the volume, we simply reduced the full volume by 10%, constraining the system to only one domain. Although the probability distributions vary with the partitioning of the system as well as the reduction of the volume, a concise statement about the change, or even the difference between the two cases may not be made from Fig. (3). A clearer picture of the problem may be obtained from Fig. (4): here we plot the mean number of kaons $`n_k`$ as the number of domains is increased (left panel), and as the the full freeze-out volume is decreased (right panel). For an unpartitioned volume with $`V_f=400V_{nuc}`$ the mean $`n47`$ Kaons. Along with this, we also plot the probability to get $`n_k+27`$ and $`n_k27`$ kaons. These probabilities sample the edges of the probability distributions, which we believe will be most sensitive to a change of the macroscopic conditions. In the following, we will denote these as “edge” probabilities. For an unpartitioned volume with $`V_f=400V_{nuc}`$ these “edge” probabilities represent the probabilites to obtain 74 and 20 Kaons. The values of such probabilities are $`10^6`$ and are measurable in experiments. One may go even further towards the edges of the distribution to observe a larger effect of partitioning the system, however the values of such probabilities may be too small to be measurable in experiments. In the figure, all values are normalized to the values of the quantities obtained from the unpartitioned maximum volume ($`V_f=400V_{nuc}`$). Note that the mean hardly changes as we increase the number of domains, or reduce the volume. The “edge” probabilities however display substantial variation. In a real heavy-ion collision, one can never determine the volume to a great degree of accuracy (in the case of flow, the very meaning of the volume becomes complicated). Hence, the similarity of the two plots makes it difficult to distinguish between the two cases. If only the kaon probability distributions were measured, sole observation of these would not allow one to distinguish between these two cases. Measurements of other observables, not constrained by strangeness conservation (such as pion yields), will have to be invoked.

Even though the “edge” probabilities in the left and right panels of Fig. (4) look substantially different, one cannot simply superpose the two figures to make a comparison, the $`x`$-axis in both graphs represent very different quantities. One may use a limiting procedure by calculating the absolute number of particles not constrained by conservation laws, such as pions. These are sensitive to the total volume and can be used to set constraints on the total volume. This limits the range of the $`x`$axis of the right panel of Fig. (4). Such a constraint will set limits on the maximum spread between the two “edge” probabilities. If the spread between the “edge” probabilities in the left panel went beyond the limits set by the right panel, then one would be lead to the conclusion that a large thermalized system was not formed. Thus, to determine the amount of partitioning of the total volume, one would need to determine the pion or proton yields with an error somewhat smaller than the range of the $`x`$-axis of Fig.(4). Thus, discering the degree of thermalization will require accurate experimental results.

As mentioned before, a change in the total volume would lead to a noticeable change in the pion number. The kaon-to-pion ratio, however, would remain almost constant in the case of a single domain, as both are proportional to the volume. The behavior of this ratio versus the volume, for a single domain, is plotted in Fig. (2). In Fig. (5) (right panel), we plot a ratio of probabilities versus the $`K/\pi `$ ratio from both scenarios. The choice of which two probabilities to use is arbitrary. We chose the numerator near the peak of the distribution $`(P(47))`$, while the denominator is chosen near the edge $`(P(70))`$. In the following, this will be denoted as the “peak-edge” ratio. The motivation behind this choice is that as the probability distributions (in Fig. (3)) shift towards the left, the value near the peak changes marginally in comparison to the value at the edge, which drops rapidly. Experimentally, this will correspond to the case where the temperature $`T`$ has been measured from other ratios. This temperature $`T`$ is then used to plot the curves in Fig. (5). The vertical line represents the case of particle production from an unpartitioned volume. The $`K/\pi `$ ratio is almost a constant as expected. The “peak-edge” ratio simply grows as the volume is reduced. The slanted set of lines represents the case for multiple domains. The $`K/\pi `$ ratio drops as the number of domains increases. As the total volume is unchanged, the number of pions remains fixed. The different lines represent different cases where we have considered smaller freeze-out volumes $`V_f`$. The cases plotted are for $`V_f=400,\mathrm{\hspace{0.17em}390},\mathrm{\hspace{0.17em}380},\mathrm{\hspace{0.17em}370}V_{nuc}`$. Thus, the triangular region represents the “allowed region” where an experimental point may lie if such an observation were made. The location of the vertical line of the region depends solely on the temperature. The location of the actual experimental point allows one to decipher not only the total volume which contributes to the production of particles but also the size of the domain over which strangeness is conserved. However, as before, we note that, even though the “peak-edge” ratio varies over many orders of magnitude; the $`K/\pi `$ ratio shows only a 15 % variation. Thus, both the $`K/\pi `$ ratio and $`T`$ have to be determined with an error much less than 15 %. This reinstates our earlier result that accurate experimental results are required to quantitatively describe the degree of thermalization.

The importance of an accurate measurement of the pion yield may be further underlined if we plot the “peak-edge” ratio versus only the kaon yield, this is done in the left panel of Fig. (5). The solid line represents the variation of the “peak-edge” ratio of probabilities with the mean number of Kaons emanating from multiple domains in a full volume of $`400V_{nuc}`$. The dot-dashed line represents the same comparison but for a total volume descending from a maximum value of $`400V_{nuc}`$. Note that both curves lie almost on top of one another. The “peak-edge” ratio drops slightly more sharply for the partitioned case. The main reason for the two plots being similar is, not an ill choice of the probability ratio, but rather, the mean number of kaons, which varies almost linearly with change in volume for volumes in the range of $`350400V_{nuc}`$. This allows the kaon number to change by almost the same range in both plots. If we were to proceed to a much smaller number of kaons the two graphs would then separate from each other. This demonstrates, yet again, the similarity of the probability distributions in Fig. (3).

The primary reason for the requirement of accurate experimental data is the low sensitivity of the kaon yields to the domain volume. At temperatures in the vicinity of 170MeV a large number of kaons ($`50`$) are produced. This causes the kaon yields to approach the grand canonical limit, in other words the kaons yields do not experience sufficient canonical suppression. A similar trend has recently been observed in a partonic scenario at higher temperature pra03 , where very small volumes are required to observe canonical suppression. Naïvely, one would expect multiple strange particles to be subject to greater canonical strangeness conservation. We present the probability distributions for particles of strangeness -2 (i.e., $`\mathrm{\Xi }`$ particles) and strangeness -3 (i.e., $`\mathrm{\Omega }`$ particles). These are plotted in the left panels of Figs. (6,7). We note that, as expected, the difference between the distributions is more pronounced for the $`\mathrm{\Xi }`$ particles, and even more for the $`\mathrm{\Omega }`$ particles, as the number of domains is increased. The $`\mathrm{\Xi }`$ distributions are peaked near 1 to 3 particles. The mean number is about 2 $`\mathrm{\Xi }`$ particles. The $`\mathrm{\Omega }`$ distributions are peaked at zero and the mean number is a fraction. We note that there is an order of magnitude difference between the probability to obtain 10 $`\mathrm{\Xi }`$ particles from one domain compared to that from 51 domains. The same is true for the probability to obtain 5 $`\mathrm{\Omega }`$ particles. The minimum probabilities in this region are $`10^6`$ and hence, are in principle measurable at RHIC. It should be pointed out that, the effect of canonical strangeness conservation, for cases where we have greater than 20 partitions (or domain volumes smaller than $`20V_{nuc}`$), may be deduced by an accurate determination of the mean numbers of kaon, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ particles (see Fig. (8)). From Fig. (8) we note that a determination of the mean number of $`\mathrm{\Omega }`$ particles with an accuracy of 5% will allow one to estimate the domain volume up to a volume of $`V=20V_{nuc}`$ (this is for the case of full freeze-out volume $`V_f=400V_{nuc}`$ and a $`T=170`$ MeV ). For domain volumes greater than this, one will need to determine the mean numbers with very high accuracy (errors $`1`$ %). Alternatively, one may determine the probability to obtain 5 $`\mathrm{\Omega }`$ particles with an accuracy of 20%, or the probability to obtain 10 $`\mathrm{\Omega }`$ particles with an accuracy of 50%. Similar error bars may also be deduced for the cascades. Which observable is better suited to the determination of the domain size, now becomes a question of experiment. The behaviour of experimental data on the mean numbers of strange particles hoh03  (see also the canonical model analysis in Ref. ker02 ) indicates that the number of partitions, in central collisions, is in the range from 1 to 20 i.e., the region in Fig. (8) where the mean values reach a plateau. As a result of this, determination of the amount of partitioning of the system (i.e., resolving the domain size between $`V=20V_{nuc}`$ to domain size of the order of the freeze-out volume $`V=V_f`$) will require very accurate measurements of the means or the probability distributions.

In the above discussion, we have demonstrated that the “edges” of the probability distributions for $`\mathrm{\Xi }`$ particles and $`\mathrm{\Omega }`$ particles are sensitive to the partitioning of the system (we also pointed out that the mean yields also shift by about 20% as we go from 1 to 51 partitions). This variation in the probability distributions is an indication that these particles are experiencing canonical suppression. Recall the situation of pion production which is insensitive to the partitioning of the system, the pion yields are calculated in the grand canonical ensemble. The probability distributions from such an ensemble are Poissonian. A more quantitative measure of the departure of the $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ production from the grand canonical limit may be obtained from the ratio of the probability distributions to the corresponding Poisson distributions. Using the mean $`n`$ we may construct a Poisson distribution to obtain $`n`$ $`\mathrm{\Xi }`$ particles from $`p`$ domains as,

$$P_p^n=\frac{n^n}{n!}e^n.$$

The ratio of the actual, calculated, probability distribution to the Poisson probability distribution is plotted in the right panel of Fig. (6) for varying number of domains as indicated. We note that the probability distributions closely resemble a Poisson distribution up to the probability to obtain about eight $`\mathrm{\Xi }`$ particles. Beyond this the distributions begin to deviate from a Poisson distribution. The same ratio is also plotted for the $`\mathrm{\Omega }`$ particles in Fig. (7). The $`\mathrm{\Omega }`$ distribution remains Poissonian up to five $`\mathrm{\Omega }`$ particles.

An interesting pattern that may be noted both for the $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ distributions is that the probability distributions approach the Poisson distribution as the number of domains is increased compared to the number of particles (see the right panels of Figs. (6,7)). This effect is much more pronounced in the case of the $`\mathrm{\Omega }`$ distributions. This may be understood as follows: in the case of a single, small domain, the $`\mathrm{\Omega }`$ probability distribution achieves an almost binary configuration. The probability $`b`$ to produce zero $`\mathrm{\Omega }`$ particles is large (i.e., $`b1`$), followed by a small probability $`a`$ to produce one $`\mathrm{\Omega }`$ particle (i.e., $`a0`$). The probability to produce multiple $`\mathrm{\Omega }`$ particles is very small and hence may be ignored (i.e., $`a+b=1`$). If there are $`p`$ domains, the probability to obtain $`n`$ $`\mathrm{\Omega }`$ particles is given by the simple binary distribution,

$`P_p^n={\displaystyle \frac{p!}{n!(pn)!}}a^nb^{pn}.`$ (28)
As $`p`$ becomes much larger than $`n`$, the binary distribution turns into the Poisson distribution. This is why the ratio of probabilities in Figs. (6,7) (right panels) begins to deviate from 1 as we go to higher and higher numbers of particles. This corresponds to increasing $`n`$. Increasing the number of domains (keeping the volume of each domain constant) corresponds to increasing $`p`$. If instead, the total volume is held fixed and the number of domains is increased, then aside from the increase in $`p`$, the production of $`\mathrm{\Omega }`$ particles from each domain approaches the binary limit. As a result of these two effects, increasing the number of domains drives the probability distribution function towards the Poisson value. Similar effects in moderation may be ascribed to the $`\mathrm{\Xi }`$ probability distribution functions. Thus, when the number of domains is large compared to the number of particles produced, one obtains a Poisson like probability.

As pointed out in the previous paragraph, $`\mathrm{\Omega }`$ production from a single, small domain has a binary configuration, i.e., a very narrow distribution ($`P(n)0`$ only for $`n=0,\mathrm{\hspace{0.17em}1}`$ ). The distribution from many such domains (see left panel of Fig. (6)), is necessarily broader due to fluctuations. However, even for the case of no partitioning, the probability to produce two $`\mathrm{\Omega }`$ particles from the full freeze-out volume is suppressed by a factor of three compared to that for a single $`\mathrm{\Omega }`$, i.e., the distributions are still rather narrow (see left panel of Fig. (7)). An amusing consequence of this is that the variance of the $`\mathrm{\Omega }`$ distribution is equal to its mean. The variance is defined as

$$\sigma =n^2n^2=\underset{n}{}n^2P(n)\left(\underset{n}{}nP(n)\right)^2,$$

where $`P(n)`$ is the probability to obtain $`n`$ $`\mathrm{\Omega }`$ particles. This equality would indicate that the $`\mathrm{\Omega }`$ yields are tending towards the grand canonical limit, whereas the situation is quite the opposite. This is illustrated in Fig. (8), here we plot the mean and variance of the the Kaons, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ particles. For a grand canonical ensemble (no suppression), the distribution is Poissonian. As a result, the variance equals the mean. As we note from Fig. (8), the $`\mathrm{\Omega }`$ particles seem to exhibit a Poissonian variance, whereas the $`\mathrm{\Xi }`$ particles deviate slightly from this limit. The Kaons display a noticeably different variance from the Poissonian value. As pointed out previously, the reason behind this rather contradictory result is the fact that the $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ distributions descend rapidly from their peak values (see left panels of Figs. (6,7)). The decent is much more mild for the kaons (see left panel of Fig. (3)). All three types of particle distributions show a Poissonian behavior near the “peak” of the distribution, deviating from such behavior as one proceeds to the “edges” of the distribution. The Kaon distribution resembles a Poisson distribution over a much wider range than the $`\mathrm{\Xi }`$ or $`\mathrm{\Omega }`$ distributions (see Fig. (9)). While the $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ distributions predominantly sample the region in the vicinity of the “peak”, the variance of the kaons receives support from a wider region. This causes the $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ variance to approach that of a Poisson distribution, while the kaon variance falls below this limit. Hence, the appearance of a Poisson like variance is not indicative of a lack of canonical suppression, and vice-versa.

A study of Fig. (8) indicates that the Kaon variance does not tend towards its mean, even as the volume is increased. This may be understood by a simple analytical argument. Imagine a large freeze-out volume ($`V_f>400V_{nuc}`$) at a sufficiently high temperature ($`T170`$ MeV). If the only particles carrying a conserved strangeness quantum number are kaons and anti-kaons, one would expect a large number of these to be produced ( $`N_K=N_{\overline{K}}100`$). In fact, the mean number of kaons $`N_K`$, anti-kaons $`N_{\overline{K}}`$ and the total, mean, number of strangeness carrying particles $`N`$ can now be estimated by a grand canonical ensemble. Strangeness conservation, however, is imposed by the following relation

$`N_K+N_{\overline{K}}`$ $`=`$ $`N`$
$`N_KN_{\overline{K}}`$ $`=`$ $`0.`$ (29)
We may square the above equations, take the mean and sum them to obtain,

$`2N_K^2+2N_{\overline{K}}^2=N^2.`$ (30)
We may also take the mean and then square and sum to obtain,

$`2N_K^2+2N_{\overline{K}}^2=N^2.`$ (31)
We may now subtract Eq. (31) from Eq. (30) to obtain the variance. As the total number of particles carrying strangeness, $`N`$, is not a conserved quantity, its distribution may be estimated accurately by a grand canonical ensemble (i.e., a Poisson distribution). Hence, the variance of $`N`$ is equal to its mean, i.e.,

$`\sigma _N=N=N_K+N_{\overline{K}}.`$ (32)
From the second line of Eq. (29) we obtain

$`N_K`$ $`=`$ $`N_{\overline{K}}`$
$`\sigma _{N_K}`$ $`=`$ $`\sigma _{N_{\overline{K}}}.`$ (33)
From Eqs. (30-33) we obtain,

$`\sigma _{N_K}={\displaystyle \frac{N_K}{2}}.`$ (34)
Thus, the kaon variance is equal to half of that of a Poisson distribution, even for a large system. The presence of particles other than kaons in the partition sums calculated previously allow for more fluctuation of the conserved strangeness among different strangeness carriers. This increases the variance from $`\frac{N_K}{2}`$ to the value plotted in Fig. (8).

It should be pointed out that measurements of the variance of the $`K/\pi `$ ratio at the SPS indicate a Poisson value, at forward rapidities afa01 (central region in fixed target experiments). This has been shown to be consistent with a Poisson variance for the kaon number jeo99 . This is in contradistinction with Fig. (8). This indicates that there is a contamination of the kaon number from neighbouring rapidity intervals. This contamination is essentially random from event to event and thus leads to an increased fluctuation of the kaon numbers. This will cause the variance to increase and lead it to becoming Poissonian. If instead, it were possible to measure the kaon yields over a larger acceptance, the percentage of random fluctuation in the kaon variance would reduce and one should begin to approach the value in Fig. (8). This serves as an important test to ensure the observation of strangeness production from a closed volume; i.e., a region whose strange yields are not influenced (event by event) by strange yields from neighbouring regions.

## V Summary and Conclusions

The principal question addressed in this article has been the degree of thermalization achieved in heavy-ion collisions, and whether this question may be answered within the framework of statistical models (i.e., without invoking hard or penetrating probes). We began with the understanding that complete equilibrium and simple superposition of proton-proton collisions are extreme possibilities. We quantified the degree of thermalization as the size of the domain within which thermalization is achieved, independent of many other such domains which may be formed. The domain was defined as the region over which chemical equilibrium (in particular strangeness chemical equilibrium) was achieved. Given the fact that chemical equilibrium is achieved at a slower rate than kinetic equilibrium, estimating the size of such a domain would set a lower bound on the size of the kinetically equilibrated region.

We further restricted the definition of the domain as the volume within which strangeness in strictly conserved. A more complete calculation with baryon number also strictly conserved may also be performed; this has been left for a future effort maj03 . Mean baryon number is set by means of a chemical potential. In this preliminary effort, a study of the effect of varying chemical potential was not performed; $`\mu `$ was set to zero throughout. The full freeze-out volume was postulated to consist of many such domains. Throughout the calculation, the full freeze-out volume $`V_f`$ was set to $`400V_{nucleon}`$ and the temperature to $`170`$ MeV. These chosen values of the thermodynamic parameters $`T,V_f,\mu `$ are approximately in agreement with collisions at RHIC. As a result, the domain size $`V`$ became $`V_f/p`$, where $`p`$ is the number of domains. The estimation of the domain size was performed by noting the behavior of the total probability distributions of various strange particles as a function of $`p`$. These results are essentially contained in Figs. (3-8). These probability distributions are, in principle, measurable in heavy-ion collisions.

In the calculation of all the above quantities (e.g., mean, variance, probability distributions etc.), we have made express use of the recursion relation technique cha95 ; das98 ; das99 , grafted from intermediate energy heavy-ion physics. This technique greatly simplifies numerical calculations in the canonical ensemble with exact conservation conditions. Two different sets of recursion relations have been devised: the first one calculates the partition functions, and its derivatives in a given domain (Eqs. (11,17,19)), while the second convolutes the probability distributions from different domains (Eq. (27)). Without the use of recursion relation techniques, such probability distribution calculations may become prohibitively difficult.

It is clear from the behavior of the mean numbers of Kaons, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ particles in Fig. (8) and the numbers seen in experiments hoh03 , that domain sizes of the order of a nucleon volume are eliminated. The question reduces to resolving the domain size between $`V=20V_{nuc}`$ ($`p=20`$) to domain size of the order of the freeze-out volume $`V=V_f`$. A definite result for $`V`$ will involve an accurate determination of the probability distributions of the Kaons, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ particles, and comparison with Figs. (3-7). Alternatively, one may determine the yields with much greater accuracy than presently imposed in heavy ion collisions and compare with Fig. (8). Such a measurement will set a lower bound on the size of the region over which thermal equilibrium has been established. We also pointed out that measurements of two particle correlations, such as the variance, may lead to incorrect conclusions regarding the amount of canonical suppression faced by different strangeness carriers.

As a corollary, we note that partitioning the system introduces a new parameter $`p`$, which essentially changes the yields of strange particles keeping the pion yields fixed. Hence this provides a natural formalism to understand the difference between the $`K/\pi `$ ratio in central collisions of small systems compared to non-central collisions of large nuclei, with the same $`N_{part}`$ in both cases (see e.g., Ref. kam03 ). Fitting the experimental points involves introducing various refinements (e.g., Baryon number conservation, temperature and volume fluctuations among domains etc.) to the rather simple model presented above. These and other considerations will be addressed in a future effort maj03 .

## VI Acknowledgment

The authors wish to thank X. N. Wang for helpful discussions. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada and by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, and by the Office of Basic Energy Sciences, Division of Nuclear Sciences, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.


format: text

text:
# References

A correlation of the cosmic microwave sky with large scale structure

Stephen Boughn1 & Robert Crittenden2

1Haverford College, Haverford PA 19041 USA

2Institute of Cosmology and Gravitation, Portsmouth PO1 2EG UK

The fluctuations in the cosmic microwave background (CMB) have proved an invaluable tool for uncovering the nature of our universe. The recent dramatic data provided by the WMAP satellite  have confirmed previous indications that the expansion of the universe may be accelerating, driven by a cosmological constant or similar dark energy component. One consequence of dark energy is the suppression of the rate of gravitational collapse of matter at relatively recent times. This causes fluctuations in the CMB to be created as the photons pass through nearby large scale structures, a phenomenon known as the integrated Sachs-Wolfe (ISW) effect. The result is additional large scale fluctuations in the CMB which are correlated with the relatively nearby (i.e., at redshift $`z1`$) matter distribution . Here we report evidence of correlations between the WMAP data and two all sky probes of large scale structure, the hard X-ray background observed by the HEAO-1 satellite  and the NVSS survey of radio galaxies . Both observed correlations are consistent with an ISW origin, indicating that we are seeing the impact of dark energy on the growth of structure.

In the standard model of the origin of structure, most of the fluctuations were imprinted on the CMB at the epoch of last scattering, when the universe was 400,000 years old ($`z1100.`$) The ISW effect induces extra fluctuations only when matter domination ends and the dark energy becomes important dynamically ($`z1`$.) When this happens, the gravitational potentials of large, diffuse concentrations and rarefactions of matter begin to decay and the energy of photons passing through them changes by an amount that depends on the depth of the potentials. The amplitude of these ISW fluctuations tends to be small compared to the fluctuations originating at the epoch of last scattering except on very large scales. However, since ISW fluctuations were created more recently, it is expected that the CMB fluctuations should be partially correlated with tracers of the large scale matter distribution, e.g., with the distribution of distant galaxies.

Detecting the relatively weak correlation of the CMB with the distribution of galaxies requires nearly full sky surveys out to redshifts $`z1`$. Focus has thus has been on luminous active galaxies, which are believed to trace the mass distribution on large scales. While active galaxies emit at a wide range of frequencies, the most useful maps are in the hard X-rays (2-10 KeV), where they dominate the X-ray sky, and in the radio, where the number counts are dominated by sources at $`z1`$. The full sky map of the intensity of the hard X-ray background (XRB) made by the HEAO-1 satellite  and a map of the number density of radio sources provided by the NRAO VLA sky survey (NVSS)  are the two best maps to search for this effect. In order to predict the expected level of the ISW effect in these two surveys, it is essential to know both the inherent clustering of the sources and their distribution in redshift. The former can be determined from previous studies of the X-ray and radio auto-correlation functions  while the latter can be estimated from deep, pointed surveys. The detailed considerations for calculating the amplitude of the ISW effect can be found elsewhere.

Previous attempts at measuring the correlation by cross correlating the CMB with both the HEAO-1 and the NVSS maps have yielded only upper limits . Here we repeat these analyses with the WMAP satellite data. We use two CMB maps generated from the WMAP data: the “internal linear combination” (ILC) map ; and the “cleaned” map of Tegmark, de Oliveira-Costa and Hamilton . The dominant “noise” to the cross correlation signal is due to the fluctuations in the CMB itself and is well characterized, while the instrument noise is negligible on the angular scales relevant to the present analysis. To reduce possible residual Galaxy and nearby radio source contamination, these maps were masked with the most aggressive mask provided by the WMAP team. The X-ray map was similarly masked and, in addition, was corrected for several large scale systematics including a linear drift of the detectors . Finally, the NVSS catalog was corrected for a systematic variation of source counts with declination . The sky coverage was $`68\%`$ for the CMB maps, $`56\%`$ for the NVSS map, and $`33\%`$ for the X-ray map. The aggressive masking of the X-ray maps was chosen because of the low resolution of the map ($`3^{}`$) and the larger number of foreground sources. Significantly, however, our results are relatively insensitive to the corrections and level of masking of the maps.

The maps were expressed as 24576 $`1.3^{}\times 1.3^{}`$ pixels in an equatorial quadrilaterized spherical cube projection on the sky . The cross correlation function (CCF) of the X-ray and CMB maps was computed according to

$$CCF(\theta )=\frac{1}{N_\theta }\underset{i,j}{}(I_i\overline{I})(T_j\overline{T})$$
(1)
where the sum is over all pairs of pixels, $`i,j`$, with angular separation $`\theta `$, $`N_\theta `$ is the number of pairs of pixels separated by $`\theta `$, $`I_i`$ is the X-ray intensity of the $`i^{th}`$ pixel, $`T_j`$ is the temperature of the $`j^{th}`$ pixel of the the CMB map, and $`\overline{I}`$ and $`\overline{T}`$ are the mean values of these quantities. The CCF of the NVSS and CMB maps is similarly defined. Figures 1, 3 are plots of the CCF’s for angular separations from 0 to 30 degrees along with the expected ISW signal. The CCF’s using the ILC and “cleaned” CMB maps were entirely consistent with each other as expected and we plot the averages of the results for the two maps. Four hundred Monte Carlo simulations of the CMB sky were also cross correlated with the actual X-ray and radio maps. The $`rms`$ values of these trials are taken as the errors for each point. These errors are highly correlated due to the large scale structure in the maps. For comparison, the errors were also computed from the data themselves by rotating the maps with respect to each other so that for small angular separations, the pixel pairs were far removed from one another on the sky. The two error estimates were consistent. The significance of the detection of the ISW effect can be estimated by fitting the amplitude of the theoretical profile to the data. Minimum $`\chi ^2`$ fits indicate that the significance of the signal is roughly $`2.02.5\sigma `$ for the radio/CMB CCF and $`2.53.0\sigma `$ for the X-ray/CMB CCF, depending on the number of bins included in the fits. These values are consistent with number of times the CCF($`\theta `$) of the Monte Carlo trials exceed the values of the measured CCF($`\theta `$) for the first few data points where the ISW signal is largest.

The possibility that the signal is arising from unmasked foreground sources is unlikely since the observed correlations do not change when the Galactic cuts are varied nor when the aggressiveness of the masking is changed. It is significant that the amplitudes of the CCF’s are essentially unchanged when no masking at all is applied to the CMB maps. It is also unlikely that a small number of diffuse foreground sources is responsible for the signal. When the different hemispheres of the maps are analyzed separately the results are consistent. The fluctuations caused by inverse Compton scattering of CMB photons (as they travel through the ionized intergalactic medium of rich clusters of galaxies), the Sunyaev-Zeldovich effect, is expected to dominate only at much smaller angles and be anti-correlated with the matter distribution. Finally, it is possible that microwave emission from the radio/X-ray sources themselves could result in the positive correlation of the maps. However, extrapolations of the frequency spectra of the radio galaxies indicate that the microwave emission is much smaller than the observed signal . In addition, the clustering of radio sources is on a much smaller angular scale than the apparent signals, and any such contamination would appear only at $`\theta =0^{}`$ and not at larger angles.

While the detection of ISW type signals in the two CCF’s is only at the $`2`$ to $`3\sigma `$ level, we find it significant that both these signals are consistent with the theoretical predictions with no free parameters. However, it should be emphasized that these two measurements are by no means independent. Both of them use the same CMB maps and, furthermore, the X-ray background is highly correlated with the NVSS radio sources. If the distribution of sources of one of the maps by chance coincided with the fluctuations of the CMB, one would expect the other map also to be correlated with the CMB to some degree. The coincidence of the expected amplitudes of the two signals is encouraging but by no means definitive since the signal to noise is not large. On the other hand, the detection of the ISW type signal in both CCF’s gives a strong indication that they are not due to some unknown systematic effects in the maps. The radio and X-ray data were gathered by quite distinct methods and it would be surprising if unknown systematics in the two maps were correlated in any way. We will present more detailed analysis of all these issues elsewhere.

We tentatively conclude that we have observed the integrated Sachs-Wolfe effect. If so, these observations offer the first direct glimpse into the production of CMB fluctuations and provide important, independent confirmation of the standard cosmological model, dominated by dark energy. It is worth noting that the power spectrum of CMB fluctuations does not show evidence of increased power on large angular scales ($`\theta >20^{}`$) as predicted by the ISW effect, but rather indicates that there is a significant deficit of power at these angular scales. Our observed CCF’s also appear to drop off somewhat faster in $`\theta `$ than do the theoretical predictions indicating perhaps a deficit of large scale power in the CMB fluctuations produced locally. However, for $`\theta >8^{}`$ both the signal and theoretical prediction are less than $`1\sigma `$ so quantitative assessments are problematic. We also note that the WMAP-NVSS results have been independently analyzed by the WMAP team.

We are grateful to Mike Nolta, Lyman Page and the rest of the WMAP team, as well as Neil Turok, Bruce Partridge and Bruce Bassett, for useful conversations. RC acknowledges financial support from a PPARC fellowship.

Communicating author: Robert Crittenden, Robert.Crittenden@port.ac.uk


text:
# Bondi Accretion and the Problem of the Missing Isolated Neutron Stars

## 1 Introduction

Neutron stars (NSs) are a common endpoint of stellar evolution. Nucleosynthesis constraints on Galactic chemical evolution require that the total number of NSs in our Galaxy today be on the order of $`10^9`$ (Arnett, Schramm & Truran 1989). After they are born, NSs are visible for a short while either as X-ray sources from their cooling radiation or as radio sources from pulsar activity. These phases, however, only last $`10^6`$ and $`10^7`$ yr, respectively, which is a small fraction of the age $`10^{10}`$ yr of the Galaxy. After a few tens of millions of years both the thermal and the rotational energy supply will be exhausted, and from then on the NSs would just be cold, invisible objects, if it were not for accretion from the interstellar medium (ISM), which makes them shine again in the X rays, as it was realized some 30 years ago by Ostriker, Rees & Silk (1970). The problem received close attention in the early 1990s. According to the calculations of Treves & Colpi (1991) and Blaes & Madau (1993; BM in the following), the luminosity due to accretion can be significant, and as many as 10,000 NSs should be detectable with the ROSAT all sky survey.

The problem of the “Missing Neutron Stars” became evident when the observations returned only a handful of possible candidates. Subsequent revisions of the expected rates have somewhat alleviated the problem by accounting for other effects not included in the BM calculations, such as dynamical heating (Madau & Blaes 1994), coupling of the dynamical and magnetorotational evolution (Popov et al. 2000; Livio, Xu & Frank 1998; Colpi et al. 1998). Furthermore, in the last few years, more accurate estimates of pulsar velocities (Lorimer, Bailes & Harrison 1997; Hansen & Phinney 1997; Cordes & Chernoff 1998; Arzoumanian et al. 2002) seemed to point towards an increase in the mean velocity (and hence a decrease in the accretion rate) relative to the Narayan & Ostriker (1990) distribution used by BM. Whereas larger typical velocities imply smaller accretion luminosities, on the other hand, as we show here (see also Zane et al. 1995), the inclusion of more realistic spectra which go beyond the simple blackbody spectrum, tends to increase the predicted X-ray number counts111This is due to the fact that, for the same total luminosity, the spectrum from a light element atmosphere is significantly harder than a blackbody (see §2.3 for more details). and hence partly counterbalance the reduction due to other effects.

In this paper, we propose that a likely solution to the missing neutron stars problem lies in the use of the Bondi expression to calculate the accretion rate. Since the work of Fabian & Canizares (1988), it has been known that the accretion luminosities of supermassive black holes in most galactic nuclei are much too low to be consistent with radiatively efficient accretion at the Bondi rate from the surrounding ISM. The evidence has become compelling in recent years (Di Matteo et al. 2000, 2001; Loewenstein et al. 2001; Narayan 2002). Two solutions to this problem have been proposed (see e.g. Narayan 2002 for a review): (i) Radiatively inefficient accretion, as in an advection-dominated accretion flow (Ichimaru 1977; Rees et al. 1982; Narayan, Yi & Mahadevan 1995; Fabian & Rees 1995); (ii) Accretion at below the Bondi mass accretion rate (Blandford & Begelman 1999; Narayan, Igumenshchev & Abramowicz 2000; Quataert & Gruzinov 2000). The current consensus is that both effects operate, with radiatively inefficient accretion triggering a reduction in the mass accretion rate $`\dot{M}`$. The case for a reduced $`\dot{M}`$ has become especially strong with the recent detection of linear polarization at 230 GHz from Sgr A at the Galactic Center (Bower et al. 2003). This observation indicates that the gas density close to the black hole is much less than that predicted for Bondi accretion.

In the context of the missing neutron star problem, we note that radiatively inefficient accretion alone will not solve the problem since the advected energy will be radiated once the accreting gas hits the surface of the NS. A reduction in the mass accretion rate relative to the Bondi rate is essential. Such a reduction has been seen in recent MHD simulations of spherically symmetric accretion with weak, dynamically unimportant magnetic fields (Igumenshchev & Narayan 2002). The radial structure of the accretion flow is dramatically modified from the Bondi profile, and the mass accretion rate is decreased significantly. Other work on rotating flows also gives similar results (§2.3 below).

Backed by the above observational and theoretical results, we perform a new calculation of the count rates for isolated NSs that accounts for the recent estimate of the velocity distribution, more realistic atmospheric spectra, as well as the correction to the Bondi rates predicted by the numerical simulations. We find that the resulting numbers are consistent with null detection at the ROSAT sensitivity.

## 2 Dynamical evolution, emission properties and accretion rates

The computation of the number counts requires several ingredients, and we describe each of them in detail in the following.

### 2.1 Dynamical evolution and local space density of INSs

As discussed in §1, in the last few years several groups (Lorimer, Bailes & Harrison 1997; Hansen & Phinney 1997; Cordes & Chernoff 1998; Arzoumanian et al. 2002) have realized that the velocity distribution for pulsar velocities derived by Narayan & Ostriker (1990; NO in the following) overpredicts the fraction of slow moving stars and underpredicts the contribution from faster moving ones. Here, we adopt the estimates by Cordes & Chernoff (1998, CC in the following). Using a sample of 49 young pulsars, they found that the distribution $`p(v_{Ri},v_{\varphi i},v_{zi})`$ for the velocities of their progenitors can be best described by a two-component gaussian model with characteristic speeds of $`175`$ km/s and $`700`$ km/s, representing $`86\%`$ and $`14\%`$ of the population, respectively. The scale height for both components (assuming a gaussian distribution) is $`0.13`$ kpc. The radial birth distribution is assumed to follow that of the stellar disk, that is an exponential with scale length 4.5 kpc out to a maximum radius of 20 kpc.

For a given initial phase space distribution $`f_i(R_i,z_i,v_{Ri},v_{\varphi i},v_{zi})`$, the equilibrium distribution can be determined by Monte Carlo simulations of the orbits of the stars in the galactic potential (e.g. Paczynski 1990; Blaes & Rajagopal 1991). However, this is a very time-consuming method and, as BM demonstrated, an excellent approximation to the numerical simulations is provided by the use of Liouville’s theorem under the “thin-disk approximation”, where the motions along $`R`$ and $`z`$ can be decoupled. In this case, the equilibrium phase space distribution is related to the initial one through the relation

$$𝑑R𝑑z𝑑v_R𝑑v_\varphi 𝑑v_zf_{\mathrm{eq}}(L_z,E_R,E_z)=𝑑R_i𝑑z_i𝑑v_{Ri}𝑑v_{\varphi i}𝑑v_{zi}f_i(R_i,z_i,v_{Ri},v_{\varphi i},v_{zi}),$$
(1)
where the integrals are taken over the regions of phase space around given values of $`L_z`$, $`E_R`$, and $`E_z`$.

Figure 1 shows a comparison between the equilibrium speed distribution derived from the NO initial distribution, and the one derived from the CC birth distribution. It can be seen that there is a significantly larger fraction of high velocity stars. When using the birth parameters derived by CC, the local NS density $`n_{\mathrm{NS}}(R=8.5\mathrm{kpc},z=0)n_{\mathrm{NS}}^0`$ is reduced from the value $`n_{\mathrm{NS}}^07\times 10^{13}N_{\mathrm{tot}}`$ $`\mathrm{pc}^3`$ found by BM for the NO population, to $`n_{\mathrm{NS}}^04\times 10^{13}N_{\mathrm{tot}}`$ $`\mathrm{pc}^3`$.

It should be noted that the relation in equation (1) does not account for the fraction of high velocity stars that escape the Galaxy. As a matter of fact, this has very little effect on the statistics of the high-$`\dot{M}`$ tail of the accretion rate distribution, which is heavily dominated by the low-$`v`$ tail of the velocity distribution. What could be affected more significantly is the local density $`n_{\mathrm{NS}}^0`$ of stars, although this too is mostly determined by the slowest local stars. In fact, since most NSs are born in the vicinity of the Galactic plane and move away from it, the slowest stars, having smaller orbits, will pass more frequently through the detectable volume. A proper computation of the local density corrected for the fraction of escaped stars, would require a knowledge of the birth rate distribution during the lifetime of the Galaxy (as the oldest stars are more likely to have escaped than the younger ones). However, it is possible to place an upper limit to the error made in the local number density by neglecting the fraction of escaped NSs. This can be done by assuming that all NSs were born at $`t=0`$ (Galaxy birth time), and therefore all the ones with velocity larger than the escape velocity have abandoned the Galaxy by the present time. Using again equation (1) with a cutoff in the initial velocity distribution for $`v>v_{\mathrm{esc}}`$, and for an escape velocity $`500`$ km/s (Cordes & Chernoff 1998), we find that the correction to the local number density is within 10%. This upper limit is well within the uncertainties due to other effects, and in particular the details of the low-$`v`$ tail of the velocity distribution which, according to the estimates made by other groups (e.g. Arzoumanian et al. 2002) can be more substantial222The Arzoumanian et al. distribution has also a more substantial high velocity tail, and they find that about 50% of all pulsars will escape the Galaxy..

For the total number $`N_{\mathrm{tot}}`$ of neutron stars populating the Galaxy at the present time, we use the estimate $`N_{\mathrm{tot}}10^9`$ based on nucleosynthesis constraints (Arnett et al. 1989), and in particular on the condition that the number of Type II SNe in the past be sufficient to account for the total mass fraction of Oxygen observed in the Galactic disk today. Note that this implies that either the Type II SN rate was much higher in the past than it is today333The currently observed Type II supernova rate from external Sb and Sc galaxies is about one every 50 years (e.g. Vandenberg, McClure & Evans 1987)., or that the observed rate today is only a small fraction of the actual core collapse events444Arnett et al. argue that this could be the case if many core collapse events were no more luminous than SN 1987A, or if a large fraction of them occurs in dense, obscuring regions, as in the case of Cas A..

### 2.2 Emission properties of accreting NSs

Previous investigations of the expected source counts from accreting isolated NSs have assumed a blackbody spectrum, although it has been recognized, since the early work of Zeldovich & Shakura (1969), that the spectrum of accreting NSs will depart from a purely Planckian shape. Zampieri et al. (1995) worked out the case of pure H atmospheres with negligible fields at low luminosities, and found that the emerging spectrum is harder than a blackbody at $`T_{\mathrm{eff}}`$, and the hardening increases as the temperature decreases. The spectral shape, as well as this behaviour with $`T_{\mathrm{eff}}`$, is very similar to that found in unmagnetized, H atmospheres around cooling NSs (e.g. Romani 1987; Zavlin et al. 1996), as the hardening is not related to the heat source but to the frequency dependence of the free-free opacity. The detailed spectra resulting from accretion at low rates onto a strongly magnetized NS have been computed by Zane et al. (2000). They found that the hard tail present in nonmagnetic models with comparable luminosity is suppressed, and the X ray spectrum is much closer to a blackbody.

For a given accretion luminosity, the geometry of the incoming flow is strongly dependent on the strength of the magnetic field of the star. In the simplest model (Davidson & Ostriker 1973), one can assume that the radius of the accretion column is delineated by the field lines that pass through the Alfven radius. This leads to an area for the accretion region $`10^8B_{12}^{4/7}\dot{M}_{10}^{2/7}`$ cm2, although it does not account for the fact that the Rayleigh-Taylor instability can lead to penetration of the accretion material through the field lines, resulting in much larger emitting areas (Arons & Lea 1980). For definiteness, in our calculations we assume a polar cap of area $`1`$ km2.

Given the considerations above, for the purpose of our analysis we use spectra of unmagnetized, hydrogen atmospheres when considering a model in which the emission occurs isotropically over the whole star (as would be typical of low magnetic fields), and blackbody spectra in a model where accretion occurs through the polar caps (at high magnetic fields). Note that, if accretion occurs through a fraction $`f`$ of the area of the star, then the corresponding temperature at the star surface will be larger by a factor $`f^{1/4}`$. For completeness, we also consider another spectral model which is a combination of isotropic emission and a power-law. This is empirically motivated, as many accreting systems display evidence for both a power law and a thermal component in their spectra (White, Stella & Parmar 1988), although note that the formation of standing shocks during the process of accretion can also give rise to significant power law tails (Zeldovich & Shakura 1969). For definiteness and in agreement with observations (Campana & Stella 2000), when considering this model, we assume that about half of the accretion energy gets reprocessed into a thermal component, and half comes out in a non-thermal, powerlaw component. For the latter, we adopt a photon index of 1.5, up to an energy of $`10`$ keV.

### 2.3 Accretion rates

The standard Bondi formula for the accretion rate onto the surface of a compact star of mass $`M`$ is given by the expression

$$\dot{M}_{\mathrm{Bondi}}\frac{4\pi G^2M^2\rho }{(v^2+c_s^2)^{3/2}},$$
(2)
where $`c_s`$ and $`\rho `$ are the sound speed and density of the ambient ISM, and $`v`$ is the velocity of the star with respect to the medium. Previous investigations of the statistics of accreting NSs have been based on the Bondi formula to compute accretion rates555See however Rutledge (2001) for a discussion of a subpopulation of highly magnetic $`(B10^{14}G)`$, accreting NSs whose accretion rate scales as $`v^{1/3}`$ rater than $`v^3`$ as in Bondi.. However, as discussed in §1, a measurement of the gas density in the vicinity of the black hole in Sgr A indicates that the Bondi expression highly overestimates the true accretion rate in that object.

The observational result on Sgr A is supported by three-dimensional MHD simulations by Igumenshchev & Narayan (2002) and Proga & Begelman (2003) which show that, in the presence of a weak, dynamically unimportant, large-scale magnetic field in the external medium666In their simulations, Igumenshchev & Narayan assumed that the initial interstellar field had magnetic pressure on the order of 1 percent of the gas pressure. If the field is stronger, then the effects are only more dramatic., the structure of a spherically-symmetric radiatively inefficient accretion flow is dramatically different from the Bondi prediction. As the gas flows in, the magnetic field is amplified. When the field reaches equipartition, it begins to reconnect, heating the gas, and setting up turbulent motions. These motions, plus the tension in the amplified field, strongly suppress the accretion rate onto the black hole. The resulting accretion rate scales roughly as

$$\dot{M}\left(\frac{R_{\mathrm{in}}}{R_a}\right)^p\dot{M}_{\mathrm{Bondi}},$$
(3)
where $`R_aGM/c_s^2`$ is the “accretion radius” at which the ambient material begins to feel the gravitational influence of the accreting star, $`R_{\mathrm{in}}`$ is an effective inner radius which is of order a few to few tens of Schwarzschild radii, and the index $`p`$ is uncertain but expected to be less than or of order unity.

The work described above refers to the infall of non-rotating or slowly-rotating magnetized gas, and is most relevant for the problem at hand. However, the main result from the simulations, namely that the accretion rate is strongly suppressed, is more general and has been noted also in many studies of strongly rotating flows. The evidence has come from a variety of hydrodynamic simulations — Stone, Pringle & Begelman (1999; two-dimensional simulations with viscosity), Igumenshchev & Abramowicz (2000; 2D simulations with $`\alpha `$ viscosity), Igumenshchev et al. (2000, 3D $`\alpha `$ viscosity), Proga & Begelman (2002; 2D inviscid) — and MHD simulations — Stone & Pringle (2001; 2D simulations), Machida, Matsumoto & Mineshige (2001; 2D and 3D), Hawley, Balbus & Stone (2001; 3D), Igumenshchev, Narayan & Abramowicz (2003; 3D). All of these studies indicate a reduced $`\dot{M}`$, as in equation (3), with $`p`$ in the range $`0.51`$.

Two different explanations have been advanced for the reduction in $`\dot{M}`$ in radiatively inefficient accretion flows: (i) winds and outflows (Blandford & Begelman 1999; Narayan & Yi 1994, 1995), (ii) convection (Narayan et al. 2000; Quataert & Gruzinov 2000). In the former explanation, gas accretes at the Bondi rate at the accretion radius $`R_a`$. However, as the gas flows in towards smaller radii, much of it is lost to an outflow so that only a tiny fraction finally reaches the central star. In the latter explanation, convection sets up a barrier to accretion, such that the net accretion rate is reduced at all radii out to $`R_a`$. Note that, within the latter scenario, a stationary solution can be obtained as convection changes the conditions at $`R_a`$ by increasing the pressure and preventing the gas from falling in supersonically. The physical picture is somewhat similar to that envisaged by Ostriker et al. (1976), where preheating of the infalling gas by the emergent X rays can suppress accretion (see in this regard also Blaes, Madau & Warren 1995). Detailed MHD simulations provide support for both the outflow picture (Stone & Pringle 2001; Hawley et al. 2001) and the convection idea (Machida et al. 20001; Igumenshchev et al. 2003). It is possible that both effects operate, e.g., with convection driving an outflow.

All the simulations described above refer to accretion onto a black hole. Most of the results are expected to carry over to the NS case. If the star is unmagnetized, the estimate of $`\dot{M}`$ given in equation (3) should be valid with essentially no change. If the star is strongly magnetized, then equation (3) will only be valid under the conditions for which the star is able to accrete at a steady pace. Indeed, in order for a rotating, magnetic star to be able to accrete steadily, three conditions need to be satisfied (BM, Popov et al. 2000, Treves et al. 2000, Ikhsanov 2003). Firstly, the Alfven radius $`r_A`$ (where the magnetic energy density is comparable to the kinetic energy of the accreting gas) must not exceed the accretion radius $`R_a`$, otherwise magnetic pressure would dominate everywhere over the gravitational pull, and the star would be in the so-called georotator state777 This regime was also investigated by Toropina et al. (1991), who called it “magnetic plow”. Some accretion may occur in this regime due to 3D MHD instabilities (Arons & Lea 1976, 1980). (Liuponov 1992). This condition is typically fulfilled for the stars in the low-$`v`$ tail of the velocity distribution, which dominate the accretion statistics. Once the constraint $`r_A<R_a`$ is satisfied, there are still two barriers that need to be overcome for accretion to proceed steadily, that is: (i) at the accretion radius, the ram pressure of the accreting material must be larger than the radiation pressure from the pulsar (otherwise the star would be in the ejector phase), and (ii) the velocity of the magnetosphere at the Alfven radius must be smaller than the Keplerian velocity at that same radius. If this latter condition were not satisfied, the star would act as a propeller and prevent material from penetrating inside $`r_A`$ (Illiaronov & Sunyaev 1975).

The relative importance of the various phases discussed above during the life of the NS, depends on a combination of its velocity, spin period, magnetic field, and accretion rate. In particular, the strength and the evolution of the magnetic field in the star play an especially important role. Livio et al. (1998) and Colpi et al (1998) first pointed out that significant field decay can cause the NS to linger in the ejector or propeller state for a time longer than the lifetime of the Galaxy. The fraction of stars which can be in the accreting phase at the present time is therefore sensitive to both the initial value of the magnetic field, and to its decay time. For an initial field strength on the order of $`10^{12}`$ G, and an intermediate value of the decay time $`\tau _d10^9`$ yr, they found that about half of the stars are allowed to accrete at the present time. As a matter of fact, the question of whether the magnetic field of the stars decays is, by itself, a very controversial issue, both on theoretical and on observational grounds. Theoretical predictions range from an exponential/power-law decay of the field (e.g. Ostriker & Gunn 1969; Goldreich & Reisenegger 1992; Urpin & Muslimov 1992) to little or no decay at all within the lifetime of the Galaxy (e.g. Romani 1990; Srinivasan et al. 1990). Observationally, low magnetic fields ($`10^9`$ G or so), possibly the result of a decay, have been inferred in low-mass X-ray binaries and millisecond pulsars, but no firm conclusions have been established yet for isolated objects. However, if evolution occurs but still leaves a residual field $`B10^810^9`$ G, a substantial fraction of the population should still be able to accrete at the present time (BM, Colpi et al. 1998). Keeping these uncertainties in mind, here we assume that all the NSs have reached the accretion phase at the present time and are accreting at the rate given by equation (3). As a matter of fact, it is sufficient that this condition holds only for a substantial fraction of the the low-$`v`$ tail of the velocity distribution, which is what dominates the accretion statistics.

For the calculations at hand, we need an estimate of the factor by which the mass accretion rate $`\dot{M}`$ onto the NS is suppressed relative to the Bondi rate. This requires an estimate of $`R_{\mathrm{in}}`$ and $`p`$ in equation (3). Abramowicz et al. (2002) found that $`R_{\mathrm{in}}`$ is of order 50 Schwarzschild radii for a rotating convection-dominated flow; $`R_{\mathrm{in}}`$ is probably somewhat smaller, by a factor of a few, for a non-rotating flow. A more serious uncertainty is in the value of $`p`$. The outflow model does not make any specific prediction other than that $`p`$ should lie somewhere in the range 0 to 1 (Blandford & Begelman 1999). The idealized convection-dominated accretion flow model, on the other hand, predicts a precise value: $`p=1`$ (Narayan et al. 2000; Quataert & Gruzinov 2000). Numerical simulations give a range of values for $`p`$, between about 0.5 and 1; the deviation from unity may be explained within the convection model as the result of “vertical” rather than “radial” convection (Igumenshchev et al. 2003). Finally, we note that a recent analysis of Sgr A (Yuan et al. 2003) indicates that an accretion model with $`p0.3`$ fits the observed spectral data and polarization constraints fairly well.

Given the uncertainties, especially in the value of $`p`$, it is not possible to derive a precise estimate of $`\dot{M}`$ for the accreting isolated NS problem. The only statement one can make with confidence, based on the evidence at hand, is that a radiatively inefficient accretion flow will not have $`\dot{M}`$ equal to the Bondi rate; the accretion rate is always suppressed, but it is not clear by precisely how much. We thus try a range of suppression factors in this paper: $`\dot{M}=(10^2,10^3,10^4)\dot{M}_{\mathrm{Bondi}}`$, with $`10^3\dot{M}_{\mathrm{Bondi}}`$ as our best-guess estimate. The selected range corresponds roughly to $`p0.51`$.

## 3 Source counts: theory confronts observations

To derive the probability distribution for the accretion rate, $`f(\dot{M},r)`$, at distance $`r`$ from us, we also need an estimate for the distribution $`n_{\mathrm{HI}}(r)`$ of hydrogen in the ambient medium surrounding the solar region. Measurements of interstellar absorption (Paresce 1984; Warwick et al. 1993) suggest that the local interstellar medium consists of smoothly distributed gas with density $`0.07`$ cm-3 up to a distance of about 100 pc, and denser gas, $`0.51`$ cm-3 beyond that. As a representative line of sight to compute source counts per unit solid angle, we will consider a line of sight with $`n_{\mathrm{HI}}(r)=0.07`$ cm-3 for $`r100`$ pc and $`n_{\mathrm{HI}}(r)=1`$ cm-3 for $`r>100`$ pc. Whereas this is a good approximation as an average in the local bubble, as a matter of fact, the ISM is multiphase and comprised of cold (neutral) gas, as well as warm and hot (ionized) gas. The hotter component is more tenuous, and as such the accretion luminosity for stars embedded within it is further reduced. On the other hand, hot gas does not contribute much to photoelectric absorption, and the resulting larger sampling depth along ionized lines of sight tends to counterbalance the reduction in luminosity. Here we have assumed for simplicity all the material to be neutral. A detailed modelling of the various ISM phases, which would vary depending on the line of sight under consideration, is beyond the scope of this work.

Let now $`C`$ be the detector count rate. The number of sources along a given line of sight which are detectable above the threshold $`C`$ is

$$\frac{dN}{d\mathrm{\Omega }}(>C)=n_{\mathrm{NS}}^0_0^{\mathrm{}}𝑑\dot{M}_0^{d(\dot{M},C)}𝑑rr^2f(\dot{M},r),$$
(4)
where the maximum distance $`d(\dot{M},C)`$, up to which sources above the threshold $`C`$ can be detected, is obtained by solving the implicit equation

$$C=\frac{1}{4\pi d^2(\dot{M},C)}_0^{\mathrm{}}L(\nu ,\dot{M})A_{\mathrm{eff}}(\nu )e^{\tau (\nu )}\frac{d\nu }{h\nu }.$$
(5)
Here $`L(\nu ,\dot{M})`$ is the luminosity of the star, $`A_{\mathrm{eff}}(\nu )`$ is the energy-dependent effective area of the detector under consideration, and $`\tau (\nu )`$ is the photoelectric absorption to the source, which we compute using the fit to the cross sections by Morrison & McCammon (1983). For the ROSAT PSPC response, we have taken the version pspcb-gain1-256888 ftp://legacy.gsfc.nasa.gov/caldb/data/rosat/pspc/cpf/matrices/ and used it in Eq.(5) by performing an interpolation to the data.

As a first step towards estimating the contribution of the various effects described above to the source statistics, in Figure 2 we compare the number counts obtained under the same assumptions made by BM (i.e. NO pulsar velocity distribution, Bondi accretion rate, and blackbody spectra) with those obtained assuming the CC velocity distribution and atmospheric spectra but still with the Bondi accretion rate. As the figure shows, despite the reduction in the number counts due to the higher mean velocity, the predicted rates are still well above the limits set by the ROSAT All Sky Survey. In particular, in the case of isotropic emission with atmospheric spectra, the reduction in the counts due to the higher mean velocities is considerably offset by the increase due to the atmospheric rather than BB spectra. In fact, as explained above, for the same value of the star temperature, the atmospheric spectra are harder than BB, hence producing more emission in the X-ray band.

In Figure 3, we show the predicted source counts once the correction to the Bondi rate discussed in §2.3 is introduced. Note that, while in Figure 2 we showed the counts per unit solid angle, here we have considered the numbers over the entire sky. This is because, in this case, the distances over which the sources are detectable are comparable to the size of the local bubble, over which the density distribution can be taken as isotropic, and therefore the line of sight that we have used can be considered as representative of the whole sky. It can be seen from Figure 3, that the predicted counts are consistent with a null detection at the ROSAT PSPC sensitivity.

As discussed in §2.3, a rough estimate for the reduction in the Bondi rate is $`10^3`$; however, for completeness, we explore the effects on the source counts for the range of $`\dot{M}10^210^4\dot{M}_{\mathrm{Bondi}}`$. This is shown in Figure 4, for the case of polar cap emission and the CC velocity distribution. Whereas for $`\dot{M}10^3\dot{M}_{\mathrm{Bondi}}`$ the counts are consistent with no detection at the ROSAT sensitivity, a rate $`\dot{M}10^2\dot{M}_{\mathrm{Bondi}}`$ would still require that a good fraction of the sources is not in the accretion phase for most of the time.

## 4 Summary

We have revisited the problem of the statistics of old and isolated neutron stars, which ought to be detectable in X-rays as a result of accretion from the ISM. We have used the recent estimates for the velocity distribution made by Cordes & Chernoff (1998), we have gone beyond the blackbody assumption for the spectra, using atmospheric models, and we have used a “revised Bondi” accretion rate, according to the results of MHD simulations and supported by direct observations in the case of accretion around supermassive black holes. Whereas the first two effects together do not result in a significant reduction of the predicted number counts, the modification to the Bondi formula dramatically reduces the number of observable sources.

As discussed in §2.3, other effects can also influence the detection statistics. Particularly important is the role played by the magnetic field of the star. This determines the spin evolution of the star itself, and the consequent fraction of the lifetime during which the star can be in a regime of steady accretion, rather than in other phases such as ejector or propeller. The importance of these various phases during the evolution of the NS life strongly depends on the evolution of the magnetic field within the star, which is by itself a very controversial issue. The source counts that we have computed here should therefore be considered as upper limits in that they are representative of a population that is accreting for most of its lifetime.

The number of observable sources that we have derived, with the inclusion of the revised velocity distribution, atmospheric spectra and reduced Bondi rate, is consistent with the detection of only a few candidates in the ROSAT all sky survey. As of now, there are only two objects which have been confirmed to be isolated NSs (see Rutledge et al. (2003) for a comprehensive review of the observations so far); moreover it is hard to firmly establish that they are actually powered by accretion and not by cooling. In the X-ray range, the spectrum from an accreting star looks quite similar to that of a cooling NS; however, Zane et al. (2000) noted that the two scenarios can be distinguished from optical observations, as accretion spectra have been found to display an excess emission over the best fitting blackbody X-ray distribution, which, on the other hand, does not appear in cooling models with similar temperatures.

We thank Omer Blaes, Bob Rutledge, Roberto Turolla and Silvia Zane for very useful discussions related to various aspects of this work. We particularly thank Omer Blaes for providing us with a routine for the space density distribution. We also thank the referee for helpful comments on the manuscript. RP aknowledges support by the Harvard Society of Fellows and a research grant by the Harvard-Smithsonian Center for Astrophysics. RN was supported in part by NASA grant NAG5-10780.