text:
# Alternatives to Quintessence Model Building

## I Introduction

Currently available observations, especially from high-$`z`$ type Ia supernovae combined with cosmic microwave background (CMB) results Jaffe et al. (2001), suggest that about one third of the critical energy density of the Universe is in the form of ordinary matter (including classical dark matter), while the remaining two thirds are in an unclustered form which is commonly called dark energy. Among other effects, this unknown component produces a recent accelerated expansion, a behaviour which standard decelerating Friedmann models are unable to match—though see Meszaros (2002). The cosmological constant $`\mathrm{\Lambda }`$ is arguably the simplest candidate for this dark energy, although it is well known that theoretical predictions for its value are many orders of magnitude off from observationally acceptable values.

Noteworthy among the many proposed alternatives are time varying scalar fields, dubbed quintessence Caldwell et al. (1998); Wang et al. (2000). Quintessence models typically involve a scalar field function and in some cases more than one. However, quintessence models often suffer from a major problem that also afflicts the cosmological constant: fine-tuning. This is often referred to as the ‘why now?’ problem: why is the cosmological constant (or a quintessence field) so small, and why does it become dominant over the matter content of the Universe right about the present day? There are so-called ‘tracking’ models where one obtains that quintessence energy density is reasonably independent of initial conditions, but on the other hand one does have to tweak some parameters in the scalar field potential in order to obtain the desired behaviour, so this can’t really be claimed as an satisfactory solution to the problem.

On the other hand, given that one has yet to see a scalar field in action, it is clear that all such toy models are not much better justified than a classical cosmological constant (despite some claims to the contrary). This is further compounded by the fact that, given some time dependence for the scale factor and energy density, one is always able to *construct* a potential for a quintessence-type model which reproduces them (see for example Padmanabhan (2002)). One is therefore reminded of Occam’s razor and can legitimately ask if observational data provide any strong justification for them, as compared to the conceptually simpler cosmological constant.

Here we make a contribution to this ongoing discussion, by studying two particularly illuminating such toy models. Cardassian expansion Freese and Lewis (2002); Freese and Kinney (2002) has recently been suggested as a model for an accelerating flat Universe without any use of a cosmological constant or vacuum energy whatsoever but solely depending on a purely matter driven acceleration. This has been accomplished by the use of a modified version of the Friedmann equation where an additional *empirical* term has been added. Unfortunately, as we show here, the Cardassian model doesn’t bring anything particularly new since, for most practical purposes, it reduces to a class of quintessence models. An example of an alternative to quintessence is the so-called Chaplygin gas Kamenshchik et al. (2001), which obeys a physically exotic equation of state (although it can be motivated, at some level, within the context of higher dimensions brane theory). In this case we will show that, as expected, current observations restrict the model parameters to a gravitational behaviour that is very similar to that of a cosmological constant.

## II Quintessence Models

Quintessence is a generic name given to a fluid component, other than ordinary pressureless matter and radiation, parameterized by an equation of state of the form $`p_i=\omega _i\rho _i`$ (in fundamental units) where $`\omega _i`$ is related to a scalar field function. Given that our chances of distinguishing between several scalar field models yielding various different time varying $`\omega _i`$ seem rather bleak Kujat et al. (2001); Maor et al. (2001), we will for the time being take $`\omega _i`$ to be constant. We shall start by reviewing some basic properties of this class of models.

### II.1 Basic Dynamical Properties

The total energy density appearing in the Friedmann equation can be expanded into a sum of quintessential fluid components, therefore turning

$$\mathrm{H}^2=\frac{8\pi G}{3}\rho \frac{k}{a^2},$$
(1)
where $`a`$ is the scale factor, $`\mathrm{H}=\dot{a}/a`$, the dot stands for time derivative and $`k`$, the curvature term, into

$$\mathrm{H}^2=\underset{i}{}\mathrm{\Omega }_i\frac{\rho _i}{\rho _{i0}}+\left(1\underset{i}{}\mathrm{\Omega }_i\right)a^2.$$
(2)
Here $`\mathrm{\Omega }_i=\rho _{i0}/\rho _{c0}`$ is the $`i`$-th density parameter ($`\rho _c`$ being the critical density) and the subscript 0 refers to present time. We have also taken $`a_0=\mathrm{H}_0=1`$.

Now, an adiabatic flow is usually assumed for each individual fluid component, so that $`d\left(\rho _ia^3\right)+p_id\left(a^3\right)=0`$. This in turn implies, for each fluid

$$\rho _i=\rho _{i0}a^{3(1+\omega _i)}=\rho _{i0}\left(1+z\right)^{3(1+\omega _i)},$$
(3)
where $`z=1/a1`$ is the usual redshift parameter. Of course, for $`\omega _i=0`$ and $`\omega _i=1/3`$ we retrieve the usual relations for dust and radiation. Inserting this adiabatic condition into Eq. (2) we finally obtain

$$\mathrm{H}^2=\underset{i}{}\mathrm{\Omega }_ia^{3(1+\omega _i)}+\left(1\underset{i}{}\mathrm{\Omega }_i\right)a^2.$$
(4)
Notice also how for $`\omega _i=1`$ a constant $`\mathrm{\Omega }_i`$ term appears, formally acting as a $`\mathrm{\Lambda }`$ cosmological constant.

A physical restriction on $`\omega _i`$ stems from the square of the sound speed velocity of each fluid component, given by $`p_i/\rho _i`$ for constant entropy. The speed at which information is carried by the fluid must necessarily be smaller or equal to unity, so that $`\omega _i1`$ in the present case. We emphasize that this constraint is only valid for constant $`\omega _i`$ since in general $`p_i/\rho _i`$ will involve an additional term due to a possible dependence of $`\omega _i`$ on $`\rho _i`$.

The dynamics of a collection of two or more different fluids can easily be studied if we interpret the Friedmann equation as an energy integral of motion of a one-dimensional fictitious particle moving with an $`a`$ coordinate. Consider Eq. (2) rewritten in the following form

$$\dot{a}^2\underset{i}{}\mathrm{\Omega }_ia^{\left(1+3\omega _i\right)}=\left(1\underset{i}{}\mathrm{\Omega }_i\right),$$
(5)
and compare it to the standard equation of motion, $`E=K+V`$, of the fictitious particle. This allows us to relate the curvature term to the mechanical energy of the system, $`\dot{a}^2`$ to its kinetic energy and the remaining term to the potential felt by the particle. We can push forward this Newtonian analogy by calculating the force $`dV/da`$ felt by the particle and see the Raychaudhuri equation emerge. Each fluid then contributes with a partial force of

$$dV_i/da=(1+3\omega _i)\mathrm{\Omega }_ia^{(2+3\omega _i)},$$
(6)
so that fluids with $`1/3<\omega <1`$ decelerate expansion while fluids with $`\omega <1/3`$ accelerate it.

### II.2 The Cardassian Model in a Quintessence Framework

Consider the Friedmann equation for a flat matter dominated Universe, i.e., $`\mathrm{H}^2=A\rho _m`$, where $`A=8\pi G/3`$. Since $`\rho _ma^3`$ we have $`\dot{a}^2=A\rho _{m0}a^1`$, hence a decreasing function of the scale factor. Now, the easiest way of making this decreasing function of the scale factor a growing one is by empirically adding some other function that counters the $`a^1`$ decreasing term. Arguably, the simplest function one can add is a power law of the scale factor, $`a^m`$ with $`m>0`$; let us choose $`m=23n`$. So now we have the following form for the time derivative of the scale factor

$$\dot{a}^2=A\rho _{m0}a^1+Ba^{23n},$$
(7)
which we can easily show to be re-writable as

$$\mathrm{H}^2=A\rho _m+B^{}\rho _m^n,$$
(8)
where $`B^{}=B/\rho _{m0}^n`$. This is precisely the Cardassian model which consists of a modification to the Friedmann equation motivated by theories with extra dimensions.

However, since $`\rho _m/\rho _{m0}=a^3`$ we can immediately rewrite Eq. (8) in the form

$$\mathrm{H}^2=\mathrm{\Omega }_ma^3+\mathrm{\Omega }_\omega a^{3n},$$
(9)
where we have defined $`A=\mathrm{\Omega }_m/\rho _{m0}`$, $`B^{}=\mathrm{\Omega }_\omega /\rho _{m0}^n`$ and $`n=1+\omega `$ (note that $`\mathrm{\Omega }_m+\mathrm{\Omega }_\omega =1`$).

Notice the equivalence between equations (8) and (9). We can interpret the Cardassian empirical term in the modified Friedmann equation as the superposition of a quintessential fluid with $`\omega =n1`$, to a background of dust. Note that since $`m>0`$ then $`n<2/3`$, implying $`\omega <1/3`$.

We thus see that any Cardassian model can readily be expressed as a quintessence model with constant $`\omega `$ and is therefore, for most practical purposes, indistinguishable from it. We say ‘for most practical purposes’ due to the fact that Cardassian models do not specify the behaviour of cosmological density fluctuations on scales larger than the horizon, although it is assumed that Newtonian gravity holds on small scales. Cardassian models are therefore incomplete, effective toy models which describe the average universe, which is one needs for most practical purposes. On the other hand, if one wanted to calculate, e.g. to calculate the CMB anisotropy on large angular scales, one would need to go beyond this simplified procedure. it is also worth emphasizing that one can’t meaningfully claim that one model is much better justified than the other Freese and Lewis (2002), at least at this stage, since both are no more than toy models.

## III The Chaplygin Gas Model

There are various ways in which one can explain the transition from a matter dominated Universe to one with an accelerated expansion using fluid components with more exotic state equations. Examples include brane-inspired models Avelino and Martins (2002) and vacuum metamorphosis Bassett et al. (2002); another is the generalized Chaplygin gas Kamenshchik et al. (2001) obeying an equation of state of the form

$$p=\frac{A}{\rho ^\alpha },$$
(10)
where $`A`$ is a positive constant and $`\alpha 0`$. The original form of the Chaplygin gas Kamenshchik et al. (2001); Bilic et al. (2002) where $`\alpha =1`$ can be motivated by considering a $`d`$-brane in a space-time of $`d+2`$ dimensions; no such motivation exists for $`\alpha 1`$. In this case its Nambu-Goto action can be seen as describing a Newtonian fluid with an equation of state given by Eq. (10). We should note that as of now the Chaplygin gas is the only fluid known to admit a super-symmetric generalization. At a toy model level, one can trivially generalize this to a different dependence on the density, as described in Kamenshchik et al. (2001); Bento et al. (2002). We emphasize that, at this simple toy-model level, this should be seen as a generalization simply for the purpose of phenomenological analysis.

As discussed in the above references, the Chaplygin gas can be described in terms of a scalar field with a given potential. As for the quintessence case, we’re able to retrieve an integral of motion from energy-momentum conservation,

$$d\left[a^{3(1+\alpha )}\left(\rho ^{1+\alpha }A\right)\right]=0,$$
(11)
and from it obtain the equivalent expression to Eq. (3), namely

$$\rho =\rho _0\left[\overline{A}+(1\overline{A})a^{3(1+\alpha )}\right]^{1/1+\alpha },$$
(12)
where $`\overline{A}=A/\rho _0^{1+\alpha }`$. This quantity can be related to the Chaplygin gas sound speed which we can easily show to be $`v_s^2=\alpha A/\rho ^{1+\alpha }`$ so that, at present, we have $`v_{s0}^2=\alpha \overline{A}`$.

Again, restrictions must be imposed on the values of $`\alpha `$ since we expect this velocity to be limited by the speed of light. In Fig. (1) we find the sound speed evolution with the scale factor for a generalized Chaplygin gas. We note that it is bounded by $`\alpha `$ when $`a\mathrm{}`$; hence $`\alpha \overline{A}\alpha 1`$ implying $`\alpha 1`$ and $`\overline{A}1`$. If we assume the generalized Chaplygin gas to be a good approximation of the Universe for all times then these later constraints apply. On the other hand, if we chose to make no assumptions about the future behaviour of the Universe and treat the Chaplygin gas only as an *effective* model in the sense of effectively reproducing the Universe today without assuming anything beyond, then only $`\alpha \overline{A}1`$ applies. We will take this approach and use supernovae results in order to constrain the values of $`\alpha `$, in the phenomenological spirit discussed above.

One of the most interesting characteristics of the Chaplygin gas resides in its ability to mimic the cosmological constant when $`\overline{A}=1`$ and ordinary matter when $`\overline{A}=0`$. We note that the accelerating regime begins after $`a_{}`$ given by

$$a_{}=\left[(1\overline{A})/2\overline{A}\right]^{1/3(1+\alpha )}.$$
(13)
Thus, for $`aa_{}`$ we’re able to expand the energy density in Eq. (12) into $`\rho B+Ca^{3(1+\alpha )}`$ where $`B`$ and $`C`$ are constants. What this tells us is that the Chaplygin gas in the accelerating regime describes a kind of mixture between a cosmological constant and a type of matter known as stiff matter obeying the equation of state $`p=\alpha \rho `$. Note that this is again very similar to the quintessence results, though not quite identical: a Chaplygin cosmology can be interpreted as an interpolation from a dust to a de Sitter Universe. This we can see most clearly in Fig. (1) where we have also drawn the evolution of $`\omega =p/\rho `$ for a generalized Chaplygin gas which resembles a phase transition studied in Bassett et al. (2002).

## IV Constraints from Observational Results

We have shown the Cardassian model to be formally equivalent to a quintessence model if we take $`\omega =n1`$. Therefore, it is possible to use recent constraints on $`\omega `$ to place limits on the Cardassian model parameter $`n`$. Several papers using supernova data have constrained the $`\omega `$ and the dark energy density, namely Perlmutter et al. (1999a); Weller and Albrecht (2001); Huterer and Turner (2001); Wang and Garnavich (2001). Some even introduced a spatial (or temporal) dependency on $`\omega \omega (z)`$ and tried to ascertain the first expansion coefficients Goliath et al. (2001). However, this work is still severely impaired by the small data set available today.

For a flat Universe, these analysis determined with a $`95\%`$ confidence level that $`\omega <0.45`$ which includes the $`\mathrm{\Lambda }`$ cosmological constant. On the other hand, CMB restrictions Bean and Melchiorri (2002) place $`\omega `$ much closer to a $`\mathrm{\Lambda }`$ scenario, in fact, with a $`68\%`$ confidence level $`\omega `$ should be smaller than $`0.85`$. Without further observational data, no significantly different analysis from the work by previous authors can be done. Therefore, we concentrate upon the addition of a Chaplygin gas to a standard model.

We should also point out that a conceptually very different possibility for a redshift-dependent equation of state is that of an inhomogeneous universe, for example one permeated by a network of light domain walls, which divides it into regions with different cosmological parameters Avelino et al. (2001a). Since these regions tend to become more and more different as the universe evolves, it turns out to be quite difficult to exclude this type of scenario, or even distinguish it from the standard one. In particular, the presently available supernovae data still allow for significant variations in the matter and vacuum densities Avelino et al. (2001b), whereas the CMB constraints are somewhat more restrictive Avelino et al. (2002a). Note that in this type of model the equation of state changes abruptly as the line of sight goes across a domain wall separating two such regions. In this sense, this model is similar, for observational purposes, to vacuum metamorphosis Bassett et al. (2002).

### IV.1 Type Ia Supernovae

Following the release of the results from the Supernova Cosmology Project (SCP) Perlmutter et al. (1999b) and the High-$`z`$ Supernova Search Team (H$`z`$ST) Riess et al. (1998) there has been a surging interest in the study of the energy content of the Universe using Type Ia Supernovae. Here we’ll use the combined observational results from both groups, using a procedure first described in Wang (2000), to produce constraints on the parameters of the Chaplygin gas—much simpler estimates were done in Fabris et al. (2002).

As usual, the parameter fit is based upon the luminosity distance $`d_L`$ defined through $`=/4\pi d_L^2`$ where $``$ stands for the intrinsic luminosity of the source and $``$ for the measured flux. From the Friedmann metric Hogg (1999) it follows that the luminosity distance, for a flat geometry, as a function of redshift is given by

$$d_L=d_\mathrm{H}(1+z)_0^z\frac{dz^{}}{\mathrm{E}(z^{})},$$
(14)
where $`d_\mathrm{H}`$ is the Hubble distance ($`c/\mathrm{H}_0=1`$ in geometrized units) and

$$\mathrm{E}^2(z)=\mathrm{\Omega }_mx^3+\mathrm{\Omega }_{cg}\left(\overline{A}+\frac{(1\overline{A})}{x^{3(1+\alpha )}}\right)^{1/1+\alpha }$$
(15)
where $`x=1+z`$, for the case of a mixture of two fluids: one of ordinary matter and the other of generalized Chaplygin gas. The apparent magnitude $`m`$ of a supernova (a parameter more often used than the measured flux $``$ to which it is related) at a given magnitude is then given by

$$m=M+5\mathrm{log}d_L+25,$$
(16)
$`M`$ being the absolute magnitude of the supernova (related to its intrinsic luminosity $``$). Following Wang Wang (2000), we use results from both the SCP and the H$`z`$ST even though their published data sets differ in their presentation. We define the distance modulus to be

$$\mu _0=5\mathrm{log}d_L+25,$$
(17)
as presented in the H$`z`$ST comprising $`50`$ supernovae. Comparatively, SCP published its measured effective rest-frame $`B`$-magnitude $`m_b^{\mathrm{eff}}`$ for $`60`$ supernovae which relate to the H$`z`$ST results through

$$m_b^{\mathrm{eff}}=M_b+\mu _0,$$
(18)
where $`M_b`$ is the peak $`B`$-band absolute magnitude of a standard supernova. The published results of the SCP and the H$`z`$ST groups have $`18`$ common supernovae, $`16`$ of which are from the Calán-Tololo Survey Phillips et al. (1999). If we calculate $`M_b`$ by comparing results from these $`18`$ supernovae (using the results from the H$`z`$ST estimated by means of the MLCS method), we’re able to get

$$M_b=m_b^{\mathrm{eff}}\mu _0=19.33\pm 0.25.$$
(19)
Hence, assuming the value $`M_b=19.33`$ for the absolute magnitude of the supernovae, we can convert SCP results to distance modulus through Eq. (18). We then add $`42`$ of these supernovae to the data set from H$`z`$ST leaving out the $`18`$ already present, thus making a total of $`92`$ supernovae used.

Our dataset is displayed in Fig. 2, together with the ranges which the models under discussion can take within the diagram. Note that for each region a larger $`\overline{A}`$ corresponds to a larger value of $`mM`$ at high redshift. The same is true for large values of $`\alpha `$, though this effect is subdominant relative to that of $`\overline{A}`$. This reflects the fact that the luminosity distance is fairly insensitive to time variations of $`w`$ Maor et al. (2001).

Results for the constraints imposed on the values of the Chaplygin gas density and its $`\overline{A}`$ parameter, for $`\alpha =0.5,1,2`$ are shown in Fig. 3 where we have performed a $`\chi ^2`$ analysis using the combined SCP and H$`z`$ST data set of the $`92`$ measured supernovae, assuming a flat universe filled only by ordinary matter and a Chaplygin gas: $`\mathrm{\Omega }_m+\mathrm{\Omega }_{cg}=1`$. We should point out, however, that there are potential caveats to the chi-squared analysis: see for example Avelino et al. (2002b). Here we have also integrated the Hubble parameter so as to remove its uncertainties from the results.

In analogy with quintessence models, where a large amount of dark energy with a less negative $`\omega `$ plays the same role as a smaller amount of energy with a more negative $`\omega `$ Perlmutter et al. (1999b), it is perceptible that the Chaplygin gas models with a higher value of $`\alpha `$ need a greater $`\overline{A}`$ to reproduce models with an inferior $`\alpha `$, even though all $`\alpha `$ values indicate roughly the same Chaplygin gas density.

The $`\overline{A}`$ parameter is restricted with a $`95\%`$ confidence level to be in the region $`0.66<\overline{A}<1`$ in the case of a pure Chaplygin gas form, to $`0.62<\overline{A}<1`$ for $`\alpha =0.5`$ and to $`0.75<\overline{A}<1`$ when $`\alpha =2`$. Notice, however, that in the latter case the limitation $`\alpha \overline{A}<1`$ implies sound velocities for $`\overline{A}>0.5`$ greater than the speed of light. Hence the case $`\alpha =2`$ is in fact completely ruled out. Similar arguments can be used to obtain sound velocity restrictions on the values of $`\overline{A}`$ and $`\alpha `$ (more specifically, their product): it is often the case that the values given a higher probability by the $`\chi ^2`$ analysis are physically unreasonable.

We should point out, however, that since one of the strongest claims of the Chaplygin gas is that of a unified explanation for the dark matter and dark energy, one might expect that the only components of the universe would be the Chaplygin gas and standard baryonic matter. In this case, since $`\mathrm{\Omega }_b=1\mathrm{\Omega }_{cg}0.04`$, we see that the supernovae data strongly exclude a $`\mathrm{\Lambda }`$-like behaviour, that is the case $`\overline{A}=1`$.

Therefore, we lift the flatness restriction (modifying the luminosity distance definition (14) accordingly) and assume the energy content to be comprised of a pure Chaplygin gas and baryonic matter with a present day density of $`\mathrm{\Omega }_b=0.04`$. The results of this analysis are depicted in Fig. 4. A large degeneracy is clearly evident. However, the plot does seem to indicate that $`\overline{A}`$ should be around $`0.8`$ and we can say with a $`95\%`$ confidence level that $`\overline{A}>0.69`$. Note that for lower values of the Chaplygin gas density we’re forced to approach the cosmological constant $`\overline{A}1`$. Note also that for the case of a plane geometry the acceptable values for $`\overline{A}`$ already exclude a $`\mathrm{\Lambda }`$ scenario.

Finally we have tried to ascertain how future SNAP results will improve our results. Assuming a flat Universe filled with ordinary matter and a Chaplygin gas portrayed by $`\mathrm{\Omega }_m=1\mathrm{\Omega }_{cg}=0.3`$ we repeat the $`\chi ^2`$ analysis as detailed above, and show the corresponding results in Fig. 5.

Obviously the constraints are now much tighter, though it is also clear that there exists a fundamental degeneracy in the $`(\mathrm{\Omega }_{cg},\overline{A})`$ plane. Again we notice that for the case $`\mathrm{\Omega }_m=0.3`$ the model is restricted to be very similar to a standard cosmological constant $`\overline{A}1`$, while for $`\mathrm{\Omega }_m=\mathrm{\Omega }_b0.04`$ the $`\mathrm{\Lambda }`$-like case is excluded outright.

### IV.2 The matter power spectrum

The length scale of the co-moving horizon at the epoch of equality between matter and radiation is directly related to the shape of the power spectrum of matter perturbations. It is straightforward to show that this scale is proportional to the scale factor at that epoch

$$\frac{\mathrm{H}_{eq}^1}{a_{eq}}a_{eq}.$$
(20)
In the standard model it is a simple matter to show that $`a_{eq}(\mathrm{\Omega }_mh^2)^1`$. The introduction of a generalized Chaplygin gas naturally implies a change in this expression which becomes

$$a_{eq}\left(\left(\mathrm{\Omega }_m+\left(1\mathrm{\Omega }_m\right)(1\overline{A})^{1/1+\alpha }\right)h^2\right)^1,$$
(21)
where we have only assumed a flat universe. Note that since the Chaplygin gas behaves as CDM except near the present time (see Fig. 1) the shape of the power spectrum goes unmodified by the evolution of density perturbations deep in the matter era, except for an overall amplitude change. Indeed, the shape of the power spectrum remains the same as for CDM except for the slight modification that Eq. (21) introduces. It is common to express the wave number $`k`$ in units of $`\mathrm{\Omega }_mh^2`$ or $`\mathrm{\Gamma }h`$, where $`\mathrm{\Gamma }`$ is known as the shape parameter (here we ignore the small correction due to the dependence on the the baryon density). In the case where a generalized Chaplygin gas is also present, we can show by using Eqs (20-21) that the shape parameter reduces to

$$\mathrm{\Gamma }=\left(\mathrm{\Omega }_m+\left(1\mathrm{\Omega }_m\right)(1\overline{A})^{1/1+\alpha }\right)h.$$
(22)
Recent work by Percival et al. (2001) using data from the 2dFGRS galaxy survey has constrained this parameter to be of the order of $`0.2\pm 0.03`$, in agreement with preliminary SDSS results Dodelson et al. (2001). Therefore, this allows us to obtain simple bounds on $`\overline{A}`$ as a function of $`\alpha `$, $`h`$ and $`\mathrm{\Omega }_m`$. Fig. 6 shows two such bounds, for the ‘extreme’ cases $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_m=0.04`$; in both cases we have assumed a Hubble parameter $`h=0.65`$. We see that $`\overline{A}0.65`$ for all generalized Chaplygin models, a result indicative of the preference for $`\mathrm{\Lambda }`$-type scenarios. In the best-motivated case $`\alpha =1`$, the constraint is much tighter, $`\overline{A}0.85`$.

## V Results and Comments

We have studied some cosmological implications of two dark energy models presented as an alternative to the now standard cosmological constant scenario, and discussed how these toy models relate to the more familiar quintessence paradigm. We have shown that the Cardassian model is, for most cosmological purposes, a standard quintessence model, whereas the Chaplygin gas does in principle have distinguishing characteristics.

On the other hand, by using supernova and density perturbation growth constraints, we have shown that any Chaplygin gas type component must have a behaviour that is very close to that of a ‘standard’ cosmological constant $`\mathrm{\Lambda }`$, assuming that the density of ‘normal’ (clustered) matter is $`\mathrm{\Omega }_m0.3`$. Of course, this result is also known to apply for the standard quintessence models Bean and Melchiorri (2002). Indeed, if by an independent method we were able to determine the total matter content of the Universe (including dark matter) to be around $`0.3`$, then in the context of this model we would in fact require a cosmological constant so as to account for the current observational results (see Fig. 3).

Conversely, the (arguably best-motivated) case where the matter content is entirely baryonic ($`\mathrm{\Omega }_b0.04`$), so that the Chaplygin gas provides both the dark matter and the dark energy, is the one where the differences with respect to the standard case would be maximal. In this case, the $`\mathrm{\Lambda }`$-like limit of this model is already strongly disfavoured by observations. It should also be noted that future observations may be able to provide fairly tight constraints on the exponent $`\alpha `$, in particular if the degeneracies discussed above are broken.

We thus conclude that the potential relevance of the Chaplygin model, in the sense of yielding observational consequences that are significantly different from those of the simpler cosmological constant, is strongly dependent on the total matter content of the universe.

###### Acknowledgements.

We are grateful to Irit Maor and Alessandro Melchiorri for enlightening discussions of this topic. C.M. is funded by FCT (Portugal), under grant no. FMRH/BPD/1600/2000 . Additional support for this project came from grant CERN/FIS/43737/2001.


format: text

text:
# 1 Schematic picture of the leakage phenomenon. Between the Lorentz contracted remnants “B” and “A” of the nuclei which have collided charmonia move in the created parton/hadron plasma. Those ψ’s with large transverse velocities vₜ (case (1)) may leak out and escape suppression, while low vₜ particles may remain (case(2)) leading to an effective increase of ⟨p²ₜ⟩ for the surviving (observed) ψ’s.

Time structure of anomalous $`J\mathbf{/}\psi `$ and $`\psi ^{\mathbf{}}`$ suppression in nuclear collisions
Jörg Hüfnera, Pengfei Zhuangb

a Institute for Theoretical Physics, University of Heidelberg, D-69120 Heidelberg, Germany.
b Physics Department, Tsinghua University, Beijing 100084, PR China

## Abstract

The data for the mean squared transverse momentum $`p_t^2(E_t)`$ as function of transverse energy $`E_t`$ of $`J/\psi `$ and $`\psi ^{}`$ produced in Pb-Pb collisions at the CERN-SPS are analyzed and it is claimed that they contain information about the time structure of anomalous suppression. A transport equation which describes transverse motion of $`J/\psi `$ and $`\psi ^{}`$ in the absorptive medium is proposed and solved for a QGP and a comover scenario of suppression. While the comover approach accounts for the data fairly well without adjusting any parameter, the fit to the data within the QGP scenario requires to assume anomalous suppression to become effective rather late, $`34`$ fm/c after the nuclear overlap.

The discovery in 1996 of anomalous $`J/\psi `$ suppression in Pb-Pb collisions at the SPS has been one of the highlights of the research with ultrarelativistic heavy ions at CERN . Does it point to the discovery of the predicted quark gluon plasma (QGP)? Six years later the situation is still confused, since several models - with and without the assumption of a QGP - describe the observed suppression, after at least one parameter is adjusted. The data on the mean squared transverse momentum $`p_t^2(E_t)`$  for the $`\psi `$ (this symbol stands for $`J/\psi `$ and $`\psi ^{}`$) in the regime of anomalous suppression and as a function of transverse energy $`E_t`$ have received less attention - for no good reason. We claim: $`p_t^2(E_t)`$ contains additional information about the nature of anomalous suppression and may help to distinguish between different scenarios. In this paper we investigate how the time structure of anomalous suppression influences the values $`p_t^2(E_t)`$. This idea has already been considered more than 10 years ago  (c.f. also more recent works ) and is based on the following observation (Fig.1): Anomalous suppression is not an instantaneous process, but takes a certain time depending on the mechanism. During this time $`\psi `$’s produced with high transverse momenta may leak out of the parton/hadron plasma and escape suppression. As a consequence, low $`p_t`$ $`\psi `$’s are absorbed preferentially and the $`p_t^2`$ of the surviving (observed) $`\psi `$’s show an increase $`\delta p_t^2`$, which grows monotonically with the mean time $`t_A`$, when anomalous absorption acts . In this letter we propose a general formalism of how to incorporate the effect of leakage into the various models, which have been proposed to describe anomalous suppression and we extract information about the time $`t_A`$ from a comparison with experiment.

It has become customary to distinguish between normal and anomalous values of suppression $`S(E_t)=\sigma ^\psi (E_t)/\sigma ^{DY}(E_t)`$ for $`\psi `$’s produced in nuclear collisions, c.f. reviews . Here, $`\sigma ^\psi `$ and $`\sigma ^{DY}`$ are the production cross section for a $`\psi `$ and a Drell-Yan pair in an AB collision, respectively. By definition, $`\psi `$’s produced in $`pA`$ collisions show normal suppression via inelastic $`\psi N`$ collisions in the final state and normal increase of $`p_t^2`$ (above $`p_t^2_{NN}`$ in $`NN`$ collisions) via gluon rescattering in the initial state. These normal effects are also present in nucleus-nucleus collisions and happen, while projectile and target nuclei overlap. Anomalous values of $`S`$ and $`p_t^2`$ are attributed to the action on the $`\psi `$ by the mostly baryon free phase of partons and/or hadrons (we call it parton/hadron plasma) which is formed after the nuclear overlap. It may lead to deconfinement of the $`\psi `$ via colour screening in the QGP, dissociation via gluon absorption or inelastic collisions by the comoving hadrons during the later period of the plasma evolution. In this letter we describe anomalous $`\psi `$ suppression within a transport theory and apply it to two rather different scenarios: (I) Absorption involving a threshold in the energy density like in a QGP scenario, (II) continuous absorption via comovers.

We denote by $`d\sigma ^\psi /d\stackrel{}{p_t}(\stackrel{}{p}_t,E_t)`$ the cross section for the production of a $`\psi `$ with given $`p_t`$ and in an event with fixed transverse energy $`E_t`$. It can be related to the phase space density $`f^\psi `$ via

$$\frac{d\sigma ^\psi }{d\stackrel{}{p}_t}(\stackrel{}{p}_t,E_t)=\underset{t\mathrm{}}{lim}𝑑\stackrel{}{b}P(E_t;b)𝑑\stackrel{}{s}f^\psi (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b}).$$
(1)
Here, $`P(E_t;b)`$ describes the distribution of transverse energy $`E_t`$ in events with a given impact parameter $`\stackrel{}{b}`$ between projectile A and target B. We follow ref.  in notation for $`P(E_t;\stackrel{}{b})`$ and the values of the numerical constants. The function $`f^\psi (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})`$ is the distribution of $`\psi `$’s in the transverse phase space $`(\stackrel{}{s},\stackrel{}{p}_t)`$ at time $`t`$ for given $`\stackrel{}{b}`$.

We define $`t=0`$ as the time, when the process of normal suppression and normal generation of $`p_t^2`$ has ceased and denote by $`f_N^\psi (\stackrel{}{s},\stackrel{}{p}_t;\stackrel{}{b})`$ the distribution of $`\psi `$’s at this time. $`f_N^\psi `$ is taken as initial condition for the motion and absorption of the $`\psi `$’s during the action of anomalous interactions. The evolution of the $`\psi `$ is described by a transport equation

$$\frac{}{t}f^\psi +\stackrel{}{v}_t\stackrel{}{}_sf^\psi =\alpha f^\psi .$$
(2)
The time dependence arises from the free streaming of the $`\psi `$ with transverse velocity $`\stackrel{}{v}_t=\stackrel{}{p}_t/\sqrt{m_\psi ^2+p_t^2}`$ (l.h.s.) and an absorptive term on the r.h.s., where the function $`\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})`$ contains all details about the surrounding matter and the absorption process. We have left out effects from a mean field, because the elastic $`\psi N`$ cross section is very small and have also neglected a gain term on the r.h.s., because recombination processes $`c+\overline{c}+N\psi +X`$ seem unimportant at SPS energies where at most one $`c\overline{c}`$ is created per event.

Eq. (2) can be solved analytically with the result

$$f^\psi (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})=\mathrm{exp}\left(_0^t𝑑t^{}\alpha (\stackrel{}{s}\stackrel{}{v}_t(tt^{}),\stackrel{}{p}_t,t^{};\stackrel{}{b})\right)f_N^\psi (\stackrel{}{s}\stackrel{}{v}_tt,\stackrel{}{p}_t;\stackrel{}{b}),$$
(3)
which for $`t=0`$ reduces to $`f^\psi =f_N^\psi `$. If we denote by $`t_f`$ the time when anomalous suppression has ceased, $`\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})=0`$ for $`t>t_f`$, the limit $`t\mathrm{}`$ in eq. (1) can be replaced by setting $`t=t_f`$, since the distribution in $`p_t`$ does not change for larger $`t`$’s.

There is little controversy about $`\psi `$ production and suppression in the normal phase: The gluons, which fuse to the $`c\overline{c}`$, collide with nucleons before fusion and gain additional $`p_t^2`$. The $`\psi `$ on its way out is suppressed by inelastic $`\psi N`$ collisions without any change in $`p_t^2`$. Neglecting effects of formation time, one has

$`f_N^\psi (\stackrel{}{s},\stackrel{}{p}_t;\stackrel{}{b})`$ $`=`$ $`\sigma _{NN}^\psi {\displaystyle }dz_Adz_B\rho _A(\stackrel{}{s},z_A)\rho _B(\stackrel{}{b}\stackrel{}{s},z_B)`$ (4)
$`\mathrm{exp}\left(\sigma _{abs}^\psi [T_A(\stackrel{}{s},z_A,\mathrm{})+T_B(\stackrel{}{b}\stackrel{}{s},\mathrm{},z_B)]\right)p_t^2_N^1\mathrm{exp}\left(p_t^2/p_t^2_N\right),`$
where

$$p_t^2_N(\stackrel{}{b},\stackrel{}{s},z_A,z_B)=p_t^2_{NN}^\psi +a_{gN}\rho _0^1[T_A(\stackrel{}{s},\mathrm{},z_A)+T_B(\stackrel{}{b}\stackrel{}{s},z_B,+\mathrm{})]$$
(5)
with the thickness function $`T(\stackrel{}{s},z_1,z_2)=_{z_1}^{z_2}𝑑z\rho (\stackrel{}{s},z)`$. All densities $`\rho _A,\rho _B`$ are normalized to the number of nucleons ($`\rho _0`$ is the nuclear matter density). We shortly explain eqs. (4) and (5): For given values $`\stackrel{}{b}`$ and $`\stackrel{}{s}`$ in the transverse plane the $`\psi `$ is produced at coordinates $`z_A`$ and $`z_B`$ in nuclei $`A`$ and $`B`$, respectively. On its way out, the $`\psi `$ experiences the thickness $`T_A(\stackrel{}{s},z_A,\mathrm{})`$ and $`T_B(\stackrel{}{b}\stackrel{}{s},\mathrm{},z_B)`$ in nuclei $`A`$ and $`B`$, respectively and is suppressed with an effective absorption cross section $`\sigma _{abs}^\psi `$. The two gluons which fuse carry transverse momentum from two sources: (i) Intrinsic $`p_t`$, because they had been confined to a nucleon. The intrinsic part is observable in $`NN\psi `$ collisions and leads to $`p_t^2_{NN}^\psi `$ in eq. (5). (ii) Collisional contribution to $`p_t`$, because in a nuclear collision, the gluons traverse thicknesses $`T_A(\stackrel{}{s},\mathrm{},z_A)`$ and $`T_B(\stackrel{}{b}\stackrel{}{s},z_B,+\mathrm{})`$ of nuclear matter in A and B, respectively, and acquire additional transverse momentum via $`gN`$ collisions. This is the origin of the second term in eq. (5).

The constants $`\sigma _{abs}^{J/\psi },\sigma _{abs}^\psi ^{}`$ and $`a_{gN}`$ are usually adjusted to the data from $`pA`$ collisions, before one investigates anomalous suppression. Fig. 2 shows (dashed curves) the results for normal suppression $`S^\psi (E_t)`$ and $`p_t^2^\psi (E_t)`$ calculated with $`f_N^\psi `$ eq. (4) in eq. (1). While the difference between calculation and data is enormous for the suppression, it is rather small for $`p_t^2^\psi (E_t)`$.

Since the physical origin of anomalous suppression is not yet settled, we investigate suppression $`S^\psi (E_t)`$ and $`p_t^2^\psi (E_t)`$ for two models, which have rather contradictory assumptions.

1. Threshold (QGP-) scenario: $`\psi `$’s are totally and rapidly destroyed, when they are in a medium with energy density above a critical one, and nothing happens elsewhere. As a representative model we use the approach by Blaizot et al. .
2. Comover scenario: The plasma of comoving partons and/or hadrons leads to a continuous absorption of long duration due to inelastic collisions with the comoving particles. As a representative model, we use the approaches by Capella et al.  and Kharzeev et al. .

In this letter we study the effect of leakage on the observed values of $`p_t^2`$ within two well established scenarios, using their assumptions and parameters. We do not introduce any modifications like (i) the $`p_t`$ dependence of the absorption process (i.e. $`p_t`$ dependence of $`\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;b)`$ in eq. (2)) and (ii) expansion of the plasma during absorption. Both effects may contribute to $`p_t^2`$ (we will present a qualitative discussion at the end), but both require a detailed study by themselves.

We begin with model I: In their schematic approach Blaizot et al.  include anomalous suppression via

$$f^\psi (\stackrel{}{s},\stackrel{}{p}_t;\stackrel{}{b})=\theta (n_cn_p(\stackrel{}{b},\stackrel{}{s}))f_N^\psi (\stackrel{}{s},\stackrel{}{p}_t;\stackrel{}{b}).$$
(6)
Here $`n_c`$ is a critical density and $`n_p(\stackrel{}{b},\stackrel{}{s})`$ is the density of participant nucleons

$$n_p(\stackrel{}{b},\stackrel{}{s})=T_A(\stackrel{}{s},\mathrm{},+\mathrm{})[1\mathrm{exp}(\sigma _{in}^{NN}T_B(\stackrel{}{b}\stackrel{}{s},\mathrm{},+\mathrm{}))]+(AB).$$
(7)
According to eq. (6) all $`\psi `$’s are destroyed if the energy density (which is directly proportional to the participant density) at the location $`\stackrel{}{s}`$ of the $`\psi `$ is larger than the critical density. All other $`\psi `$’s survive. While the prescription eq.(6) successfully describes the data in the full $`E_t`$ range of anomalous suppression after the only one free parameter, $`n_c`$, is adjusted, the predictions for $`p_t^2(E_t)`$ are significantly below the data, especially at large $`E_t`$ (see below).

The expression eq. (6) for the phase space distribution $`f^\psi `$ including anomalous suppression within the threshold model is recovered within our transport approach eq. (3) by setting

$$\alpha (\stackrel{}{s},\stackrel{}{p},t;\stackrel{}{b})=\alpha _0\theta (n_p(\stackrel{}{b},\stackrel{}{s})n_c)\delta (t)$$
(8)
and taking the limit $`\alpha _0\mathrm{}`$. The delta function $`\delta (t)`$ has to be included to recover the expression eq. (6) and may be understood by the physical picture that the energy density is highest at $`t=0`$ and anomalous absorption most effective then.

There are various ways to introduce another time structure into the absorption term. We will try two options, one being

$$\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})=\alpha _0\theta (n_p(\stackrel{}{b},\stackrel{}{s})n_c)\delta (tt_A).$$
(9)
The idea of a threshold density is kept, but suppression does not act at $`t=0`$ but at a later time $`t_A`$, which time is then determined from a comparison with the data. For times $`t>t_A`$, one finds from the general solution eq. (3) (and $`\alpha _0\mathrm{}`$)

$$f^\psi (\stackrel{}{s},\stackrel{}{p}_t,t)=\theta (n_cn_p(\stackrel{}{b},\stackrel{}{s}))f_N^\psi (\stackrel{}{s}\stackrel{}{v}_tt_A,\stackrel{}{p}_t),t>t_A$$
(10)
which differs from the expression (6), by the motion in phase space of $`f_N^\psi `$ during the time $`0tt_A`$. For $`t>t_A`$ the momentum distribution derived from $`f^\psi `$ does not change any more.

The suppression $`S^\psi (E_t)`$ and $`p_t^2\psi (E_t)`$ calculated with the distribution $`f^\psi `$ from eq. (10) depends on the parameter $`t_A`$. We use the parameters of ref.  where available, i.e. for $`J/\psi `$: $`\sigma _{abs}^{J/\psi }=6.4`$ mb, $`n_c=3.75`$ fm-2. The parameters for the generation of $`p_t^2`$ by gluon rescattering are taken from a fit to the $`pA`$ data : $`p_t^2_{NN}=1.11`$ (GeV$`/c)^2`$ and $`a_{gN}=0.081`$ (GeV$`/c^2)`$ fm-1. We also account for the transverse energy fluctuations  which have been shown to be significant for the explanation of the sharp drop of $`J/\psi `$ suppression in the domain of very large $`E_t`$ values, by replacing $`n_p`$ by $`\frac{E_t}{E_t}n_p`$ where $`E_t`$ is the mean transverse energy at given $`b`$. We then calculate $`\sigma ^{J/\psi }/\sigma ^{DY}`$ as a function of $`t_A`$. Since the critical density for $`J/\psi `$ is quite high, the leakage affects only the very high momentum $`\psi `$’s. Since their number is small, the calculated results for the suppression $`S^{J/\psi }(E_t)`$ depend only little on $`t_A`$ (Fig. 2).

We turn to a discussion of $`p_t^2^{J/\psi }(E_t)`$. Fig. 2 shows calculated curves for values of $`t_A=0`$ to 4 fm/$`c`$. The dotted line ($`t_A=0`$) is the result of the original threshold model with immediate anomalous suppression (and has been predicted in ). It fails badly at large values of $`E_t`$. Also no other curve with a given $`t_A`$ describes the data for all values of $`E_t`$. We have to conclude that $`t_A`$ depends on $`E_t`$: $`t_A(E_t)`$. The larger the values $`E_t`$ the later anomalous suppression acts. From a comparison with data we have $`t_A<\mathrm{\hspace{0.25em}2}`$ fm$`/c`$ for $`E_t<80`$ GeV, $`t_A2.5`$ fm$`/c`$ for $`80E_t100`$ GeV, and $`t_A3.5`$ fm/c for $`E_t>100`$ GeV.

A similar analysis is performed for the $`\psi ^{}`$. The data for the suppression are taken from , those for $`p_t^2^\psi ^{}`$ from . Since the $`\psi ^{}`$ has not been treated by Blaizot et al., we fit the suppression data and find values for the absorption parameters, $`\sigma _{abs}^\psi ^{}=7`$ mb, $`n_c^\psi ^{}=2.3`$ fm-2, leaving the parameters $`p_t^2_{NN}`$ and $`a_{gN}`$ unchanged. Since the critical density for $`\psi ^{}`$ is smaller than that for $`J/\psi `$, more $`\psi ^{}`$s leak out of the anomalous suppression region, the change in suppression due to the increase of $`t_A`$ is noticeable. Fig. 2 shows the results with a good fit to the suppression data.

The data for $`p_t^2^\psi ^{}`$ have rather large error bars. The calculated curves show a strong dependence on $`t_A`$, again with a trend that $`p_t^2`$ at larger values of $`E_t`$ require larger values of $`t_A`$. However, the numerical values $`t_A(E_t)`$ for $`\psi ^{}`$ are above those for the $`J/\psi `$ by about $`12fm/c`$. This result is strange, because we expect the $`\psi ^{}`$ to be destroyed more easily and therefore more rapidly. We will come back to this point in the conclusions.

The time structure introduced via eq. (9) is certainly oversimplified. Rather than having it act at one particular time $`t_A`$, it is more reasonable to assume that it acts during a time interval. Therefore we have also investigated the following form of the absorptive term

$$\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})=\frac{\alpha _0}{t_1t_0}\theta (n_p(\stackrel{}{b},\stackrel{}{s})n_c)\theta (t_1t)\theta (tt_0).$$
(11)
Anomalous suppression acts for times $`t`$ between $`t_0`$ and $`t_1`$ (which may be a function of $`E_t`$). Especially $`t_0`$ is an interesting quantity, because it gives the starting time for anomalous suppression.

We calculate the suppression $`S^{J/\psi }(E_t)`$ and $`p_t^2^{J/\psi }(E_t)`$ as a function of $`[t_0,t_1]`$ and for $`\alpha _0=5`$. All other parameters are unchanged. The suppression is well described for all intervals $`[t_0,t_1]`$, but the calculations for $`p_t^2^{J/\psi }(E_t)`$ depend strongly on the choice of this interval. Fig.3 shows some representative examples: The four windows in Fig.3 display events with $`t_0=0,1,2,3fm/c`$, respectively, and several values for $`t_1`$. It is obvious that there is not one curve, which describes all the data. Rather we find that the best curve for $`E_t<60`$ GeV corresponds to $`[0,1]`$, for $`60<E_t<80`$ GeV to $`[1,3]`$, for $`80<E_t<100`$ GeV to $`[2,5]`$, and for $`E_t>100`$ GeV to $`[3,5]`$. The length $`t_1t_0`$ of the interval remains approximately constant, while the beginning time $`t_0`$ increases with $`E_t`$. We will discuss the significance of these results after we have treated leakage within the comover model.

We proceed to model II, the scenario of comovers: Partons and/or hadrons which move with the $`\psi `$ may destroy the $`\psi `$ with a cross section $`\sigma _{co}^\psi `$. The comover density $`n_{co}(t)`$ depends on time, for which the Bjorken scenario of longitudinal expansion predicts $`t^1`$. The absorptive term $`\alpha `$ in eq. (2) then takes the form

$$\alpha (\stackrel{}{s},\stackrel{}{p}_t,t;\stackrel{}{b})=\sigma _{co}^\psi \frac{n_c(\stackrel{}{b},\stackrel{}{s})}{t}\theta (\frac{n_c(\stackrel{}{b},\stackrel{}{s})}{n_f}t_0t)\theta (tt_0),$$
(12)
i.e. absorption by comovers starts at $`t=t_0`$ and ends at $`t_1`$, when the comover density $`n_c(\stackrel{}{b},\stackrel{}{s})t_0/t_1`$ has reached a value $`n_f`$ independent of $`\stackrel{}{b}`$ and $`\stackrel{}{s}`$. The comover approach contains a definite time structure for anomalous suppression and we have not changed it. We also account for the transverse energy fluctuations, by replacing $`n_c`$ by $`\frac{E_t}{<E_t>}n_c`$, and the transverse energy loss induced by the $`J/\psi `$ trigger, by rescaling $`<E_t>`$ by a factor $`\frac{n_p2}{n_p}`$, which have been shown to be significant for the explanation of the sharp decrease of $`S^{J/\psi }(E_t)`$ at $`E_t>100`$ GeV. In the choice of parameters we have followed : $`n_{co}(\stackrel{}{b},\stackrel{}{s})=1.5n_p(\stackrel{}{b},\stackrel{}{s})`$ with the participant density from eq. (7), $`t_0=1`$ fm/c, $`n_f=1`$ fm-2$`\sigma _{abs}^{J/\psi }=4.5`$ mb, $`\sigma _{abs}^\psi ^{}=6`$ mb, $`\sigma _{co}^{J/\psi }=1`$ mb and $`\sigma _{co}^\psi ^{}=3`$ mb, no additional parameter has been introduced.

The calculated $`E_t`$ dependence of the suppression shown in Fig.4 fits the data acceptably well like model I. For $`p_t^2`$, we compare the results of two calculations with the data. The dashed lines in Fig. 4 show the calculation leaving out leakage. Formally this limit is obtained from eq. (3) by setting $`\stackrel{}{v}_t=0`$ in the exponent and in $`f_N^\psi `$. Due to the introduction of the $`E_t`$ loss which is necessary to recover the $`J/\psi `$ suppression for large values of $`E_t`$, the case without leakage does not fit the data $`p_t^2`$ even in the domain of low $`E_t`$ values. Only when the leakage effect is taken into account, the calculation agrees well with the data. We stress: the calculation of $`p_t^2(E_t)`$ in the comover model is a true prediction in the sense that no parameter is adjusted above those which are fitted to the suppression $`S^\psi (E_t)`$. We have also calculated the mean time $`t_A^\psi `$ for comover action by studying the suppression $`S^\psi (E_t)`$ as a function of time and taking the mean of $`t`$ with the weight $`dS^\psi (E_t)/dt`$ and find

$`t_A^{J/\psi }`$ $`=`$ $`3.5\mathrm{fm}/c`$
$`t_A^\psi ^{}`$ $`=`$ $`3.0\mathrm{fm}/c,`$ (13)
which values include the time $`t_0`$, eq. (12), between the end of normal suppression and the beginning of comover action. The values eq. (S0.Ex2) are found to be rather independent of $`E_t`$.

We summarize: In this letter we have investigated the influence of leakage on the calculation of $`p_t^2(E_t)`$ for $`J/\psi `$ and $`\psi ^{}`$ produced in Pb-Pb collisions at SPS energies. This effect is closely related to time structure of anomalous suppression. The evolution of the $`\psi `$ during anomalous suppression including leakage is described within a general transport equation. The formalism is applied to two models with rather contradictory underling physical assumptions, the threshold (QGP) and the comover models with the following results:

1. Calculations within the original models, where leakage is left out, do not describe the data for $`p_t^2^\psi (E_t)`$, the discrepancy being particularly strong at high values of $`E_t`$.
2. Including leakage into the comover model, without changing its structure nor its parameters leads to a good agreement with the data for $`p_t^2(E_t)`$ for $`J/\psi `$ over the full range of values $`E_t`$.
3. The assumption in the threshold model that anomalous suppression acts instantaneously at $`t_A=0`$, i.e. right after normal suppression is not supported by experiment. Rather for central collisions, the data of $`p_t^2`$ are described best, if one assumes that anomalous suppression acts at a time $`t_A=34`$ fm$`/c`$ after the nuclear overlap.
4. The situation for the $`\psi ^{}`$ is less clear. While the data for suppression can be fitted by properly adjusting the parameters, both models underpredict the data for $`p_t^2^\psi ^{}(E_t)`$. This is evident for the comover model. Within the threshold model the values $`t_A`$ required to fit the $`\psi ^{}`$ data are larger than those for the $`J/\psi `$, which seems unreasonable to us. It could be that the error bars on the experiments are too small.

We conclude: In this letter we have investigated how leakage (escape of high $`p_t\psi ^{}`$s) when introduced into exciting models of anomalous suppression influences the calculated values of $`p_t^2(E_t)`$. As mentioned already above, there could also be other effects which could influence $`p_t^2`$. We discuss them briefly: $`p_t`$ dependence of the mechanism responsible for anomalous suppression and transverse expansion of the plasma. Both effects can be treated within the transport approach eqs. (2),(3) by introducing an explicit dependence on $`p_t`$ and a modified dependence on $`t`$ into the function $`\alpha (s,p_t,t;b)`$. A first and schematic investigation on the $`p_t`$ has been made by Dingh . One sees without calculation that if $`\alpha /p_t<0`$, i.e. high $`p_t\psi ^{}`$s are absorbed less (whatever the mechanism may be), the calculated values of $`p_t^2(E_t)`$ increase, thus having the same effect as an increase of the time $`t_A`$ eq. (9), when anomalous suppression acts. As for the transverse expansion of the plasma, its effect is qualitatively clear: It reduces the effect of leakage since it makes it harder for the $`\psi ^{}`$s to escape. Although the qualitative changes on $`p_t^2`$ of the two effects are clear, their quantitative treatment needs a rather careful study of the underlying physics, necessitates to modify existing models, to change their input parameters and to introduce new parameters. This is beyond the scope of this letter.

Acknowledgments: We are grateful to our friends and colleagues A. Gal, J. Dolejšı, Yu. Ivanov, B.Z. Kopeliovich, H.J. Pirner and C.Volpe for help and valuable comments. One of the authors (P.Z.) thanks for the hospitality at the Institute of Theoretical Physics. This work has been supported by the grant 06HD954 of the German Federal Ministry of Science and Research and by the Chinese National Science Foundation.


text:
# STATISTICS OF NEUTRAL REGIONS DURING HYDROGEN RE-IONIZATION

## 1 Introduction

Early hydrogen reionization is a particularly interesting process in the high redshift universe and is inevitably linked to the appearance of the first star forming objects, at least those that served as sources of the ionizing radiation. If it occurred early enough, reionization imprints distinct features in maps of the cosmic microwave background (CMB) on arcminute angular scales (Vishniac, 1987; Bruscoli et al., 2000; Benson et al., 2001). Several aspects of hydrogen reionization remain uncertain despite the rapidly accumulating data on the high redshift universe. For example, it is unclear what objects produce most of the ionizing radiation, although high redshift galaxies are very strong candidates (Couchman & Rees, 1986; Haiman & Loeb, 1996; Ciardi et al., 2000). It is also unclear how the ionized regions develop in space (Miralda-Escudé et al., 2000). The ionizing sources are likely to lie in high density regions, but those regions do not necessarily ionize first; the ionizing photons may tunnel into less dense regions and ionize those first. Further, the duration of reionization is unknown and only a lower limit on the redshift marking the end of that epoch exists (Becker et al., 2001; Gunn & Peterson, 1965).

In previous papers (Benson et al. 2001, Liu et al. 2001) we examined how the CMB is affected by the reionization process and how future CMB maps can be used to extract information on that process. However, hydrogen reionization is currently best probed by spectra of high redshift QSOs. Unlike maps of the CMB which are sensitive to line of sight integrals over the density and velocity of ionized gas, QSO spectra contain direct information on the local distribution of neutral hydrogen.

Recently, Becker et al. (2001) analyzed spectra of a sample of QSOs with redshifts between $`z=5.82`$ and $`z=6.28`$. In the spectrum of their highest redshift QSO ($`z=6.28`$), the transmitted flux in the Ly$`\alpha `$ and Ly$`\beta `$ forest in the redshift stretch $`5.95<z<6.16`$ is consistent with zero, with a lower limit of 20 on the Ly$`\alpha `$ optical depth. This long stretch in redshift corresponds to a comoving distance of $`60\mathrm{h}^1\mathrm{Mpc}`$ in a universe with a cosmological constant of $`\mathrm{\Lambda }_0=0.7`$ and matter density of $`\mathrm{\Omega }_0=0.3`$. Becker et al. suggest that this long neutral region is a detection of the end of the hydrogen reionization era. To increase the signal-to-noise ratio, Becker et al. binned their spectra in $`4`$Å pixels. This prevented them from detecting small scale dark windows that are likely to appear in the spectra as left overs from the reionization epoch. A high resolution spectrum for one of the quasars observed by Becker et al. was obtained by Djorgovski et al. (2001). The spectrum of this quasar ($`z=5.73`$) was thorough analyzed by Djorgovski et al. (2001) and was found to contain several small dark windows signifying the detection of the trailing edge of the reionization epoch.

Motivated by the results of Becker et al. (2001) and Djorgovski et al. (2001), we examine here how the key ingredients in the reionization process can be probed by QSO spectra and future high redshift galaxy surveys. The methodology of the present paper has been previously developed in Benson et al. (2001, hereafter BNSL). As in BNSL, we obtain the distribution of ionized gas in an N-body simulation. Using a semi-analytic model for galaxy formation (Kauffmann et al., 1993; Somerville & Primack, 1999; Cole et al., 2000) we identify mock galaxies in the simulation and estimate the ionizing flux emitted by stars at high redshift. We then use several schemes to follow the development of ionized regions in the simulation. Using the output of this procedure we obtain predictions for two statistics. The first is the expected number of neutral segments longer than a given length, and the second is the cross-correlation function between the galaxies and neutral regions. Both of these measures rely on observations of QSO spectra, and the latter also on a sample of high redshift candidates for the sources of ionizing radiation.

## 2 Modeling the development of ionized regions

BNSL employed a semi-analytical model for galaxy formation in a high resolution N-body simulation of dark matter to estimate the amount of ionizing radiation produced by stars in high redshift galaxies. Here we follow a similar procedure using the latest version of the galform galaxy formation model (Cole et al., 2000), and the same $`\mathrm{\Lambda }`$CDM simulation as BNSL (c.f. Jenkins et al. 1998). This simulation has $`\mathrm{\Omega }_0=0.3`$, a cosmological constant $`\mathrm{\Lambda }_0=0.7`$, a Hubble constant of $`h=0.7`$ in units of $`100\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$, and is normalized to produce the observed abundance of rich clusters at $`z0`$ according to Eke, Cole & Frenk (1996). The simulation has a box of length 141.3 $`\mathrm{h}^1\mathrm{Mpc}`$ and contains $`256^3`$ dark matter particles.

Most of the ionizing photons emitted by stars are likely to be absorbed by gas and dust inside galaxies and only a small fraction, $`f_{\mathrm{esc}}`$, escapes and becomes available for hydrogen ionization in the intergalactic medium (IGM) (Leitherer et al., 1995). Assuming a value111BNSL considered several models for the variation of $`f_{\mathrm{esc}}`$ from galaxy to galaxy. Here we adopt the simplest model, in which $`f_{\mathrm{esc}}`$ is constant for all galaxies. for $`f_{\mathrm{esc}}`$, BNSL used the following models to follow the propagation of ionized regions in the simulation (see BSNL for details).

Model A (Growing front model) Ionize a spherical volume around each source (halo) with a radius equal to the ionization front radius for that halo assuming a large-scale uniform distribution of neutral hydrogen. Since the neutral hydrogen in the simulation is not uniformly distributed, and also because some spheres will overlap, the ionized volume will not contain the correct total mass of hydrogen. We therefore scale the radius of each sphere by a constant factor and keep repeating the procedure until the correct total mass has been ionized.

Model B (High density model) We simply rank the cells in the simulation volume by their density. We then completely ionize the gas in the densest cell. If this has not ionized enough hydrogen we ionize the second densest cell. This process is repeated until the correct total mass of hydrogen has been ionized.

Model C (Low density model) As model B, but we begin by ionizing the least dense cell, and work our way up to cells of greater and greater density (Miralda-Escudé et al., 2000).

Model D (Random spheres model) As Model A but the spheres are placed in the simulation entirely at random rather than on the dark matter halos.

Model E (Boundary model) Ionize a spherical region around each halo with a radius equal to the ionization front radius for that halo. This may ionize too much or not enough neutral hydrogen depending on the density of gas around each source. We therefore begin adding or removing cells at random from the boundaries of the already ionized regions until the required mass is ionized.

Guided by the observations of Becker et al. (2001) we will compute the number of neutral segments and the cross correlations at three output redshifts, $`z=6.67`$, $`6.22`$, and $`5.80`$. Table 1 lists the volume filling factors (ratio of volume of ionized regions to total volume of the simulation box) in each of our five models for $`f_{\mathrm{esc}}=0.01`$ (column 2 in table 1) and $`f_{\mathrm{esc}}=0.05`$ (column 3). We will see later that the results of Becker et al. imply that the amount of ionizing radiation increases significantly between $`z6.2`$ and $`5.8`$. Therefore we also show results for a variable escape fraction,$`f_{\mathrm{esc}}^{\mathrm{var}}`$ (column 4) which equals $`0.01`$ before $`z=6.22`$ and increases linearly with time to $`0.1`$ at $`z=5.8`$. As expected, model C (low density model) has the highest filling factor for a given $`f_{\mathrm{esc}}`$. For $`f_{\mathrm{esc}}=0.05`$ the simulation box is fully ionized at $`z=5.8`$ in all models.

## 3 The statistical measures

### 3.1 The number of neutral segments

We are now in a position to compute the proposed statistics. We begin with $`N(>L)`$, the mean number of neutral (unionized) segments of length greater than $`L`$ in a given redshift range in a line of sight (see Barkana 2002, for a similar statistic). We have the spatial distribution of ionized and neutral regions in the simulation as a function of redshift, for each of our five reionization scenarios (see Table 1). To compute $`N(>L)`$ we choose several random “lines of sight” in the output of the simulation at a given redshift. In each line of sight we identify the neutral segments and tabulate their lengths, $`L`$, in comoving $`\mathrm{h}^1\mathrm{Mpc}`$. A line of sight is obtained by starting from a grid point at the boundary of the simulation and going around the boundary of a rectangular slice of perimeter $`4\times 141\mathrm{h}^1\mathrm{Mpc}`$ (comoving) until we return to the starting point. (The shape of the path used to extract a line of sight makes no difference to our results. Using a rectangular path helps reduce the chance of pattern repetition.) This yields a total of $`256`$ lines of sight each spanning a redshift range corresponding to $`564\mathrm{h}^1\mathrm{Mpc}`$ (comoving). We then compute the mean $`N(>L)`$ from these lines of sight. For convenience we normalize $`N(>L)`$ to a redshift span of $`1\mathrm{h}^1\mathrm{Gpc}`$ by multiplying the direct result obtained from the simulation by $`1000/564`$. The $`N(>L)`$ (normalized to $`1\mathrm{h}^1\mathrm{Gpc}`$) is shown in figure (1) for $`z=6.66`$ (top), 6.22 (middle), and 5.8 (bottom). The panels to the left show $`N(>L)`$ computed for $`f_{\mathrm{esc}}=0.01`$, at these three redshifts. To the right we show curves computed with $`f_{\mathrm{esc}}=0.05`$ at $`z=6.66`$ (top) and $`6.22`$ (middle), and a variable $`f_{\mathrm{esc}}^{\mathrm{var}}`$ at $`z=5.8`$ (bottom). We have also computed $`N(>L)`$ from lines of sight each of length $`141\mathrm{h}^1\mathrm{Mpc}`$ passing through in the simulation box at random positions and found very similar results to those shown in the figure.

In observed spectra ionized regions with large optical depths can be confused with completely unionized regions, and vice versa. By measuring Ly$`\beta `$ absorption lower limits of about 20 on the Ly$`\alpha `$ optical depth can be obtained. In CDM-like models ionized regions with optical depth larger than this lower limit are small. A careful analysis and modeling of spectra can therefore help reduce this confusion. Nevertheless to estimate the degree of the confusion, we have computed $`N(>L)`$ assuming that regions with optical depth larger than 20 are identified as unionized222The neutral hydrogen fraction at each point was estimated assuming photoionization equilibrium (Theuns et al. 1998) and a hydrogen ionization rate of $`\mathrm{\Gamma }=4.3\times 10^{13}`$s-1. The Ly$`\alpha `$ line-width and effects of peculiar velocities were also taken into account.. We found that in this more detailed calculation $`N(>L)`$ is shifted to smaller $`L`$ by less than a factor of 2 at large $`L`$. At small $`L`$, $`N(>L)`$ is reduced by a factor of less than 2, although this could likely be reduced with a more careful analysis as suggested above.

The cells of our computational grid are approximately $`0.55h^1`$Mpc (comoving) in extent. According to Benson et al. (2002) the characteristic smoothing length in the IGM at this redshift is approximately $`0.2h^1`$Mpc (comoving). Thus, our current simulation does not fully resolve the structure of gas in the IGM. To assess the consequences of this limitation we repeated our calculations using a lower resolution computational grid ($`128^3`$ cells instead of $`256^3`$). We find that, at large $`L`$, the lengths of neutral regions are approximately 30% smaller when the higher resolution grid is used. Doubling the grid resolution to $`512^3`$ cells (and therefore almost fully resolving the smoothing length) should make little difference to our results. We also checked that noise due to the finite number of particles is unimportant for the results presented here.

### 3.2 The cross correlation

In addition to the distribution of neutral hydrogen, the cross-correlation, $`\xi `$, requires the spatial distribution of (potential) sources of the ionizing photons. In real observations the distribution of neutral hydrogen is obtained from QSO spectra and the positions of sources from a high redshift galaxy or QSO survey. Our goal here is simply to demonstrate that $`\xi `$ can distinguish between various models for the propagation of the ionized regions. We will therefore simply compute $`\xi `$ from the three dimensional distribution of neutral hydrogen in the simulation, i.e. we will not address the question of how well $`\xi `$ can be estimated from realistic mock observations where the neutral regions are identified in lines of sight.

The simulation box is divided into a $`256^3`$ cubic grid. At each grid point $`𝐱`$ we define a quantity $`h`$ to be zero if that point has been ionized and unity otherwise. We also use the cloud-in-cell (CIC) method to derive the number density $`n(𝐱)`$ of galaxies from the galaxy positions in the simulation. We select galaxies on the basis of their ionizing luminosity and include only those with luminosities sufficiently high to ensure the population is fully resolved in the N-body simulation. Because ionizing luminosity correlates only weakly with halo mass—since it depends so strongly on the star formation rate—this means we select only the most luminous sources—$`6`$, $`6`$, and $`10\times 10^{54}h^2`$photons/s at $`z=5.8`$, $`6.2`$ and $`6.7`$ respectively. These sources are rare and contribute only a small fraction to the total ionizing luminosity density of the universe (about 15% and 18% at z=6.222 and 5.80, respectively). Denoting the average values of $`n`$ and $`h`$ by $`\overline{n}`$ and $`\overline{h}`$ respectively, we define the cross-correlation, $`\xi `$, as

$$\xi (𝐫)=<[n(𝐱)/\overline{n}1][h(𝐱+𝐫)/\overline{h}1]>$$
(1)
where the symbol $`<.>`$ implies averaging over all grid points $`𝐱`$. In practice the calculation of $`\xi `$ is done using the technique of Fast Fourier Transforms. The results are shown in Figure 2. This figure demonstrates that $`\xi `$ is sensitive to the ionization model. It is positive for model C (low density), almost vanishes for model D (random spheres), and is negative for models A, B, and E, which ionize dense regions first. However, $`\xi `$ has a similar shape for all the high density models (A, B and E) and we expect that it will be difficult to distinguish between them in observational data. Nevertheless they all are significantly different from either model C or D. So $`\xi `$ should successfully discriminate among low density, random ionization, and high density models.

## 4 Discussion

The statistic $`N(>L)`$ giving the number of neutral segments of length greater than $`L`$ for a given total length of a QSO spectrum is sensitive to the filling factor and the way in which ionization proceeds. The cross correlation between candidate ionizing sources and neutral regions is less sensitive to the filling factor but is a more direct and robust probe of the propagation of the ionization fronts. In addition to QSO spectra, the cross-correlation requires a sample of candidates (galaxies and QSO) for the ionizing radiation. Catalogs of galaxies and QSO at high redshift are rapidly accumulating, making it possible to compute the QSO-flux and galaxy-flux correlations. Comparison between galaxy-flux and QSO-flux correlation functions will tell us whether galaxies or QSO contributed most of the ionizing radiation.

Current observations do not allow a robust determination of $`N(>L)`$, The number of spectra needed to determine $`N(>L)`$ to within a given accuracy can be estimated by noting that the relative error on this function is $`1/\sqrt{MN}`$, where $`M`$ is the number of observed QSO spectra covering the same redshift range. Nevertheless we still can make general conclusions based on the Becker et. al. (2001) result, assuming that the long Gunn-Peterson trough they observe is indeed a signature of reionization. Let us take the length of a spectrum in the Ly$`\alpha `$ forest at $`z6`$ to be $`250\mathrm{h}^1\mathrm{Mpc}`$ corresponding to the comoving distance between Ly$`\alpha `$ and Ly$`\beta `$ emission lines. An inspection of Figure 1 shows that: $`1)`$ A completely neutral stretch of a comoving length of $`60\mathrm{h}^1\mathrm{Mpc}`$ at $`z6.2`$ is inconsistent with a large filling factor. 2) The observations indicate that the chances of finding long neutral regions at $`z<5.94`$ are tiny while they are significant at higher redshift. This behavior seems inconsistent with our theoretical $`N(>L)`$ computed with constant $`f_{\mathrm{esc}}`$. If $`f_{\mathrm{esc}}=0.01`$ then there are similar probabilities for finding long segments at $`z=6.22`$ and $`z=5.80`$. If $`f_{\mathrm{esc}}=0.05`$ then the box is fully ionized at $`5.80`$, while at $`z=6.22`$ long segments are very rare. Therefore the data favor models in which there is a significant increase in the amount of ionizing radiation in the IGM, either due to an increasing escape fraction over this redshift range, or a much stronger evolution in the galaxy/QSO population than is predicted by our model. And, $`3)`$ a model in which ionization proceeds from low to high density regions seems to be inconsistent with a $`60\mathrm{h}^1\mathrm{Mpc}`$ neutral region even for escape fractions as low as $`f_{\mathrm{esc}}=0.01`$.

## Acknowledgments

This research is supported by the EC RTN network “The Physics of the Intergalactic medium”. A.N. is supported by the Israeli Academy of Science, the German Israeli Foundation for Scientific Research and Development, and the Technion V. P. R. Fund and Henri Gutwirth Promotion of Research. N.S. is supported by the Alexander von Humboldt Foundation and a Japanese Grant-in-Aid for Science Research Fund of the Ministry of Education, No. 14540290. CGL is supported by PPARC. A.N. and A.J.B. wish to thank the National Astronomical Observatory of Japan in Mitaka for its hospitality and support.