## Abstract

We present the first hydrodynamic *N*-body simulations of primordial gas clouds responsible for the reionization process in dark energy cosmologies. We compare the cosmological constant scenario with a SUGRA quintessence model with marked dynamics in order to highlight effects due to the different acceleration histories imposed by the dark energy. We show that both the number density of gas clouds and their clumpiness keep a record of the expansion rate during evolution, similar to the non-linear dark matter profile at virialization, as was recently demonstrated by Dolag et al. Varying the shape of the primordial power spectrum, we show how this effect is mitigated by a running spectral index decreasing the power at small scales. Our results demonstrate that, in order to constrain the dark energy from large-scale structures, one must track its effects down to the distribution of luminous matter.

## 1 INTRODUCTION

Different observational data, like those from high-redshift supernovae (Riess et al. 2004; Astier et al. 2006), the cosmic microwave background (CMB) (Spergel et al. 2003, 2006), and large-scale structure (Tegmark et al. 2004; Cole et al. 2005), are now giving a consistent picture of our Universe, which can be described by the so-called ‘concordance’Λ cold dark matter (ΛCDM) model: a spatially flat universe whose expansion accelerates in the present epoch because of a dominant dark energy component. Historically, the first and simplest candidate for that has been the cosmological constant Λ, corresponding to matter with an equation of state *p*=*w*ρ*c*^{2}, with constant *w*=−1. However, observations suggest a value for Λ which is more than a hundred orders of magnitude smaller than the energy scales expected to be responsible in the very early Universe.

To avoid inelegant solutions based on parameter fine-tuning, the general idea of dark energy is extended to encompass the so-called quintessence, possibly corresponding to a suitably self-interacting scalar field, whose pressure and energy density evolve during cosmic history. The currently available observational data sets do not allow yet to place strong constraints on *w*. We know that it has to come close to −1 with a precision of about 10 per cent at the present epoch (see e.g. Riess et al. 2004; Spergel et al. 2006), but its redshift evolution is still poorly constrained and can hopefully be determined only with new-generation data; see Seljak et al. (2005) for one of the very first attempts to measure the high-redshift dark energy equation of state.

One consequence of *w* > −1 earlier is that the formation of cosmic structures sets in earlier compared to the ΛCDM model. At a given redshift, this corresponds to a higher abundance of more concentrated haloes (see e.g. Dolag et al. 2004). This may open alternative ways to study quintessence models based on galaxy-cluster counts (Wang & Steinhardt 1998; Haiman, Mohr & Holder 2001; Battye & Weller 2003; Majumdar & Mohr 2003) and strong gravitational lensing (Bartelmann et al. 2003; Meneghetti et al. 2005).

The different history of structure formation must also affect the reionization epoch. To produce the required ionizing radiation, a sufficiently high number of early stellar sources are needed, such as supermassive Population III stars or more ‘standard’, massive Population II stars in protogalaxies. In particular, Population III stars are supposed to be composed mainly of hydrogen and helium with their primordial abundances, to have masses ≳100 M_{⊙} much larger than the standard Population I and II stars, and to form in dark matter (DM) haloes with typical masses of 10^{5…6} M_{⊙} (the so-called minihaloes).

The recent analysis of the temperature–polarization cross-correlation and the polarization auto-correlation measured by the *Wilkinson Microwave Anisotropy Probe* (*WMAP*) satellite now allows the estimation of the epoch when reionization occurs. The Thomson optical depth of τ∼ 0.17 ± 0.04 extracted from the first-year *WMAP* data had been interpreted as signalling the high redshift of *z*= 17 ± 5 for the primordial epoch of global reionization (Kogut et al. 2003). The recently released three-year *WMAP* data allowed an improved control in particular of the polarized foreground emission, yielding τ= 0.09 ± 0.03. This corresponds to *z*= 10^{+2.7}_{−2.3} for the completion of the reionization process, even if some level of parameter degeneracy still remains (see e.g. fig. 3 in Spergel et al. 2006).1 Different authors (see e.g. Ciardi, Ferrara & White 2003; Sokasian et al. 2003; Wyithe & Loeb 2003; Sokasian et al. 2004) suggested that these reionization data can place new complementary constraints on the cosmological parameters, and in particular on the nature of dark energy.

So far, studies of the structure formation process in dark energy cosmologies were focused on modifications of DM structures due to the different expansion histories (see e.g. Klypin et al. 2003; Dolag et al. 2004). They revealed that virialized objects keep a record of the expansion rate at the time of their formation. The general picture is that in cosmologies with *w* > −1, the increase in the dark energy density with redshift enhances structure growth, and virialized haloes become more concentrated because of the denser environment they form in (Dolag et al. 2004). This effect is stronger in tracking quintessence models, which we describe in the next section, in which the increase in the dark energy density with redshift is enhanced. Complementary hydrodynamical simulations still need to be carried out. This work is a first step in this direction, beginning with the numerical study of the dependence of the primeval gas clouds responsible for the reionization process on the cosmic expansion rate at their formation time.

Due to the higher concentration of structures in dark energy scenarios, which has so far been verified only for the DM, it is mandatory for our study to control the shape of the primordial power spectrum. When combined with the Lyα forest data and with the analysis of the 2dF Galaxy Redshift Survey, the first-year *WMAP* results supported a ‘running’ spectral index which tilts the spectrum slightly towards small scales, starting at *k*≳ 1 Mpc^{−1}. The three-year *WMAP* data are compatible with this (Spergel et al. 2006), albeit with less emphasis than the earlier results. A running index would cause a slower growth and evolution of small-scale cosmic structures compared to models with a constant index *n*. A reduction in power on small scales may alleviate several potential discrepancies in the ΛCDM model, such as the abundance of substructures in galaxies and the high central concentration of galactic haloes. However, as shown by analytic and numerical work (see e.g. Somerville, Bullock & Livio 2003; Yoshida et al. 2003b), models with a running spectral index (RSI) may have severe problems in producing enough objects to allow a global reionization at high redshift. Since this might be balanced by the enhanced structure growth due to the dark energy, it is important to jointly study these two aspects.

In this paper, we study the high-redshift structure formation in quintessence models based on the results of high-resolution cosmological *N*-body simulations combined with hydrodynamics. Specifically, we consider four different cases combining two flat cosmological models, that is, the concordance ΛCDM model and a quintessence model with a SUGRA potential, with two types of the primordial power spectrum, one with a constant spectral index *n*= 1 and an RSI model assuming the best-fitting relation found by Spergel et al. (2003). More detail will be given in the next section. This paper is organized as follows. In Section 2, we describe the general characteristics of the hydrodynamical simulations used below and introduce the quintessence models chosen. The techniques adopted for identifying DM haloes and the corresponding results are presented in Section 3. The abundances of gas clouds are presented in Section 4 together with the clumping factors and recombination times as computed from the simulation output. The implications of the previous results in terms of reionization are discussed in Section 5. The final discussion and our main conclusions are drawn in Section 6.

## 2 THE SIMULATIONS

We first review the cosmologies chosen, and then the computational and hydrodynamical features of our numerical simulations. As mentioned in the Introduction section, this work is a first step towards investigating the dependence of hydrodynamic *N*-body simulated structures on the global cosmic expansion rate, and it is appropriate to study the connection with the shape of the primordial power spectrum. Thus, we need dark energy cosmologies with different dynamics as well as constant and running index in the primordial power spectrum.

The standard, concordance ΛCDM model is our reference case. It is spatially flat, has the present matter density Ω_{0m}= 0.3, the baryon density Ω_{0b}= 0.04, and a dominating cosmological constant, Ω_{0Λ}= 0.7. The Hubble constant is assumed to be *h*= 0.7 in units of 100 km s^{−1} Mpc^{−1}. The second model is chosen so as to emphasize the dynamical effects of the dark energy. The spatial curvature is also vanishing, with Ω_{0m}= 0.3 and 70 per cent of the critical density contributed by a quintessence scalar field Φ evolving in time under the effect of a potential energy *V*(Φ). We review here the main characteristics of this model, referring to Peebles & Ratra (2003) and references therein for further detail.

The field Φ must satisfy the Klein–Gordon equation

The equation-of-state parameter*w*≡

*p*/(ρ

*c*

^{2}) is Evidently,

*w*evolves with redshift, but the behaviour of the cosmological constant is reproduced as

*w*→−1 in the static limit .

The attractive feature of such models is their attractor solutions which exist for exponential and inverse power-law potentials or combinations of those. Suitably tuning the potential amplitude, they can reach the present level of dark energy starting from a large range of initial conditions. The evolution of the equation of state may be more or less pronounced along these trajectories. We wish to have pronounced dark energy dynamics in order to clearly identify its effects. Thus, we choose a SUGRA potential

where*G*is the Newton constant, as suggested by supergravity corrections (Brax & Martin 1999, 2000). By exploiting the properties of tracking trajectories, we fix

*M*by requiring that the dark energy makes up 70 per cent of the critical density today, and α= 6.5 so that

*w*→

*w*

_{0}≡

*w*(

*z*= 0) =−0.85 today. At higher redshift,

*w*>

*w*

_{0}(Brax & Martin 2000), which meets our purpose of emphasizing differences with the ΛCDM model. It should be noted that this model is already at variance with the observational constraints (see e.g. Knop et al. 2003; Spergel et al. 2003; Riess et al. 2004).

The upper panel of Fig. 1 shows the redshift evolution of *w* for this model. The differences with the standard ΛCDM model (constant *w*=−1) are large at all redshifts. In particular, during the epoch of reionization which we are mostly interested in, the SUGRA model has *w*≈−0.35. Note that *w* strongly affects the formation history of cosmic structures through its modified expansion history. This is evident from the lower panel in Fig. 1, which shows the linear growth factor (normalized to coincide with the ΛCDM model at *z*= 0) for the model shown above. In quintessence models where the dark energy is more abundant in the past than it is now, structures tend to form earlier if the perturbation power spectrum is normalized at the present epoch. This is due to the earlier dark energy dominance, slowing down the perturbation growth. In our SUGRA model with *w*_{0}=−0.85, fluctuations at *z*= 20 are expected to be ≈15 per cent higher than that in the ΛCDM model. Consequently, we expect that non-linearity to be reached earlier, with obvious consequences for the reionization epoch.

Dynamical dark energy generally increases the concentration of virialized structures (Dolag et al. 2004), which can also be caused by a modified primordial power spectrum. Therefore, we vary the simulation set-up using two different shapes of the primordial power spectrum. Specifically, for both cosmological models (ΛCDM and SUGRA), we construct initial conditions corresponding to two scenarios, namely a power spectrum with the standard constant *n*= 1 normalized to σ_{8}= 0.9. The second model also has σ_{8}= 0.9, but allows for a reduced primordial power on small scales as suggested by the combined analysis of the first-year *WMAP* data, galaxy surveys and Lyα data (see e.g. Spergel et al. 2003), within the range compatible with the three-year *WMAP* data (Spergel et al. 2006).

The case for an RSI is interestingly controversial. Re-analyses of the Lyα data (see e.g. Viel, Weller & Haehnelt 2004; Seljak et al. 2005; Viel & Haehnelt 2006; Zaroubi et al. 2006) lend much less support for an RSI than that which was reportedbased on the first-year *WMAP* data and favour a constant primordial spectral index very close to unity (see, however, Desjacques & Nusser 2005). On the other hand, an RSI is still viable and fits the three-year *WMAP* data well (Spergel et al. 2006), while the support of a pure power law from large-scale structures persists (see e.g. Seljak, Slosar & McDonald 2006; Viel, Haehnelt & Lewis 2006). As mentioned before, we choose our cosmological models so as to cover a broad range of phenomenologies rather than to fit the existing data well. Thus, we choose to adopt a primordial power spectrum given by *P*(*k*) ∝*k*^{n(k)}, with

*k*

_{0}= 0.05 Mpc

^{−1},

*n*(

*k*

_{0}) = 0.93 and This choice conveniently allows a direct comparison with the similar analysis by Yoshida et al. (2003b), who adopted the same primordial power spectrum, but only considered the concordance ΛCDM model. It is worth recalling that this represents the best fit to the first-year

*WMAP*data and remains compatible with the three-year data, d

*n*/d ln

*k*≲−0.06 (Spergel et al. 2003, 2006). Fig. 2 shows the two primordial power spectra considered here. Evidently, they differ at large

*k*where the RSI spectrum is steeper. Less power on small scales will delay the structure formation (see also the corresponding discussion in Yoshida et al. 2003b).

Combining two cosmological models (ΛCDM and SUGRA) with two primordial power spectra (constant *n*= 1 and RSI), we ran four different simulations: ΛCDM and SUGRA with *n*= 1, and ΛCDM-RSI and SUGRA-RSI with an RSI. The initial conditions were generated at *z*≈ 120 using a comoving periodic box of 1 Mpc side length. We followed Dolag et al. (2004) in adapting the initial conditions for ΛCDM to those for the SUGRA model. Thus, by construction, only the pairs of models sharing the same primordial power spectrum have the same random phases, allowing a direct comparison of the structures formed. The models with different primordial power spectra have different random phases and can thus not be directly compared (see Fig. 3).

Since we are interested in structures during cosmic reionization, the box size used here can be considered sufficiently large to suppress cosmic variance (see the discussion on effects of finite box sizes in Yoshida et al. 2003b; Bagla & Prasad 2006). Moreover, to prevent the non-linear regime from reaching the fundamental fluctuation mode in the simulations, we stopped the simulation at the latest at *z*∼ 15.

When the gas starts cooling significantly within the first forming structure, the simulation time-step can become very small (since we neglect feedback in these simulations), effectively stopping the simulation. Accordingly, we had to stop some simulations at a higher redshift, for example, the SUGRA simulation, which could only be evolved until *z*∼ 19 because of its more rapid evolution, as we will see in a later analysis. None the less, all simulations are followed long enough for them to form the sufficient amount of structure we need for our analysis.

The simulations were performed using the gadget2 code (Springel, Yoshida & White 2001; Springel 2005) on the IBM-SP4/5 at CINECA (Bologna). gadget2 is based on the combination of a tree-particle-mesh algorithm to solve the gravitational forces, and the smoothed particle hydrodynamics (SPH) scheme to describe the hydrodynamics of gas particles. The code follows the non-equilibrium reactions of different chemical species (e^{−}, H, H^{+}, H^{−}, He, He^{+}, He^{++}, H _{2}, H^{+}_{2}), using the reaction coefficients computed by Abel et al. (1997) (for more detail, see also Anninos et al. 1997; Yoshida et al. 2003a), and adopts the cooling rate due to molecular hydrogen estimated by Galli & Palla (1998). The density field was sampled with 324^{3} DM particles and an equal number of gas particles, having a mass of ∼1040 and 160 M_{⊙}, respectively. The comoving Plummer-equivalent gravitational softening length was fixed to ε= 0.143 kpc.

In Fig. 3, we show the projected gas and (mass-weighted) temperature distributions for all four models at *z*= 20.6. As expected, the RSI simulations appear delayed, exhibiting smoother density and temperature distributions. The SUGRA models reveal their expected behaviour with a more evolved density field than their ΛCDM counterparts. Since the models with and without an RSI have different random phases in their initial conditions, they are not supposed to reproduce the same structures.

## 3 THE ABUNDANCE OF DARK MATTER HALOES

The abundance of virialized DM haloes is usually expressed by the mass function *N*(*M*). From a theoretical point of view, different relations are available to predict the redshift evolution of *N*(*M*), once the cosmological model is fixed. Some of these were obtained from a specific model for the collapse of structures (Press & Schechter 1974; Sheth & Tormen 2002); in other cases, they represent best fits to *N*-body simulations with very high resolution (Sheth & Tormen 1999; Jenkins et al. 2001; Warren et al. 2006). However, the theoretical mass functions have not yet been extensively tested in the regime of very low masses and very high redshifts which we are mainly interested in, where their predictions differ by orders of magnitude. This is due to the difficulty of combining sufficient resolution in space and mass with regions large enough to be considered a fair sample of our Universe.

The existing studies yield contradictory results. Simulations covering large volumes, such as those presented in Reed et al. (2003) (432^{3} particles in a 50 *h*^{−1} Mpc box) and in Springel et al. (2005) (2160^{3} particles in a 500 *h*^{−1} Mpc box), found that the Sheth & Tormen (1999) (hereafter ST99) formula fairly reproduces the numerical results, but they can safely identify haloes with masses ≳10^{10} M_{⊙} at *z*≲ 10. At *z*= 15, Reed et al. (2003) found a slight tendency of the ST99 model to overpredict the halo counts.

Numerical work based on much smaller boxes, such as Jang-Condell & Hernquist (2001) (128^{3} particles in a 1 *h*^{−1} Mpc box) or Yoshida et al. (2003b,c) (2 × 324^{3} particles in a 1 *h*^{−1} Mpc box), found very good agreement with the predictions of the Press & Schechter (1974) (hereafter PS74) model. Note that these simulations are more affected by their finite box size (see the following discussion). Further support for the PS74 formalism at high redshifts is provided by the analysis of Gao et al. (2005), who showed that the extended Press–Schechter theory (Bond et al. 1991) correctly predicts the structure growth and the dependence of halo abundances on the density of the surrounding region in the redshift range between 50 and 30. One of the best compromises between mass resolution and volume coverage was achieved by Heitmann et al. (2006), who used an extended set of simulations with 256^{3} particles in boxes sized between (4…126) *h*^{−1} Mpc. Their mass function at high redshift turns out to deviate significantly from the PS74 prediction for *M* > 10^{7} *h*^{−1} M_{⊙}, while the agreement with the Warren et al. (2006) (hereafter W06) relation is better. Note that the W06 formula predicts a substantial suppression of the halo abundance in the high-mass regime compared to the ST99 model. Finally, using an array of high-resolution *N*-body simulations with box sizes in the range 1–3000 *h*^{−1} Mpc, Reed et al. (2006) recently determined the mass function of DM haloes at redshifts between 10 and 30. They found that the PS74 model gives a poor fit at all epochs, while the ST99 model gives a better match to the simulations, but still overpredicts the abundance of most massive objects by up to 50 per cent.

We identify DM haloes in our simulations by running a friends-of-friends algorithm on the DM particles only and setting the linking length to 20 per cent of the mean interparticle separation. The resulting mass of each halo is then corrected by the factor (1040 + 160)/1040 ≈ 1.15 to account for the additional contribution by the gas particles. Fig. 4 displays the redshift evolution of the number of haloes with masses exceeding 7 × 10^{5} M_{⊙} for all four simulations. The figure compares the results with the predictions of the PS74, ST99 and W06 models, indicated by the solid, dashed and dot–dashed lines, respectively. Note that we do not consider the relation found by Jenkins et al. (2001) because it cannot extrapolate to the mass and redshift ranges considered here.

At first sight, the dark energy dynamics as well as the modified initial power spectrum generate significant differences. First, as expected, the number of massive haloes is always lower in the RSI models than in their scale-invariant counterparts. While the first objects with *M* > 7 × 10^{5} M_{⊙} start appearing near *z*≃ 30 in the SUGRA and the ΛCDM simulations, this is delayed to *z*≃ 25 and ≃20 for SUGRA-RSI and ΛCDM-RSI, respectively.

Secondly, the different expansion histories leave clear records in the structure population. Considering models with the same primordial power spectrum, the halo abundance is always larger in the SUGRA than in the ΛCDM models by a factor of ≃2, confirming the earlier structure formation in dynamical quintessence models expected from the modified linear growth factor in Fig. 1.

We now consider the difference of the theoretical expectations from the PS74, ST99 and W06 mass functions. All analytic relations were modified so as to account for the different expansion histories and initial power spectra in the different cosmological models. For the two simulations with a cosmological constant (ΛCDM and ΛCDM-RSI), the results indicate that the W06 formula agrees better with the simulations than the ST99 and PS74 formulae. Some small deviations are seen at very high redshifts where, however, counts are low and the statistical uncertainty is high. The ST99 mass function tends to slightly overestimate the simulation results, again mainly at high redshift, where the discrepancies with W06 are more evident. In contrast, the PS74 mass function always underestimates the halo abundances, with differences reaching up to a factor of ∼3. Thus, our results for the cosmological-constant models agree with the previous analyses by Reed et al. (2003), Springel et al. (2005) and Heitmann et al. (2006). However, this conclusion is in conflict with other studies (Jang-Condell & Hernquist 2001; Yoshida et al. 2003b,c).

The situation is different for both SUGRA models. The mass function in the simulation with constant *n*= 1 is always larger than all theoretical predictions. In particular, the ST99 and W06 relations give very similar estimates which fall below the simulation results by 30 per cent, while the difference with PS74 is always a factor of 3–4. The tendency of PS74 to underestimate the numerical mass function is also confirmed for SUGRA-RSI (a factor of 3–5 everywhere), while the agreement for W06 and mostly for ST99 is substantially better. The deviations reach 20 per cent only at low redshift. These results show that caution must be applied when applying the standard theoretical mass functions to very high redshift in quintessence models.

We must finally remark that, owing to the much higher mass resolution, the mass range covered in our simulations is quite different compared to earlier work. At the same time, our results could be affected by the small box size. In fact, as discussed in Yoshida et al. (2003b), the ‘effective’ variance, that is, the amount of power on the scales covered by the simulation, can be significantly smaller than that obtained by integrating the power over all scales. This difference can systematically delay structure formation and possibly cause the discrepancies between the simulation results from different box sizes. To quantify the effect of the finite box, we apply the analytic expressions proposed by Bagla & Prasad (2006), which allow an evaluation of the required correction to the mass function. For example, assuming the PS74 approach and the redshift range covered by the data shown in Fig. 4, we find that the correction for *M*= 7 × 10^{5} M_{⊙} leads to a slight underestimate of 10–15 per cent, which does not significantly affect our previous results. Similar conclusions can be reached on the ST99 relation.

## 4 THE PROPERTIES OF THE GAS DISTRIBUTION

We now show our results on the statistics and the physical properties of the gas distribution. We concentrate on two aspects, namely the formation time of the gas clouds and their clumping factor, focusing on their dependence on the underlying dark energy cosmology and primordial power.

### 4.1 The formation of gas clouds

We first study the formation of the gas clouds which are responsible for reionization. For this purpose, we run a friends-of-friends algorithm to group cold and dense gas particles, assuming a linking length of 1/20th of the mean interparticle separation to locate the densest gas concentrations. Following Yoshida et al. (2003b,c), we only consider clouds exceeding 3000 M_{⊙}, which corresponds to groups with at least 19 gas particles in our simulations.

The results are shown in Fig. 5 in terms of the redshift evolution of the number of gas clouds, *n*_{clouds}, identified in the four simulations. It is not surprising that the trend between the models is similar to that for the DM haloes in Fig. 4. Gas clouds start forming earlier and are more abundant in the SUGRA than in the ΛCDM model with the same primordial power spectrum.

The suppression of power on small scales is, however, the dominant effect. In fact, the two RSI models have their structure formation delayed and have the smallest number of gas clouds in the final output. At *z*≃ 19, there are order-of-magnitude differences between the models. We find nearly 300 clouds for SUGRA, nearly 50 for ΛCDM, seven for SUGRA-RSI and no clouds for ΛCDM-RSI. Also, the formation of the first gas clouds strongly depends on the cosmological parameters. It happens much before *z*= 30 for SUGRA, while it is *z*≃ 18 for ΛCDM-RSI. The reason for this pronounced difference is the integral nature of the observable *n*_{clouds} considered here. The production rate of clouds per unit time is the key quantity, which differs between the cosmologies considered because of the different expansion histories and primordial power spectra. The SUGRA model produces clouds earliest, followed by the ΛCDM, while the same models with an RSI stay in the same order but the structure formation process is delayed because of the modified primordial power. The integrated cloud count, *n*_{clouds}, accumulates these differences during the whole epoch of cloud production.

Finally, we have to note that in our study we neglected the feedback effects which might influence the structure formation history of primordial clumps. For example, radiative feedbacks, such as photoionization or molecule photodissociation, could stop cooling suppressing the formation of small objects; mechanical feedbacks and shocks can determine gas removal from the structure, on the one hand, or induce further shell fragmentation, on the other hand. Many studies on this subject are present in the literature (see the corresponding discussion in Ciardi & Ferrara 2005, and references therein), but up to now comprehensive conclusions have not been reached: in fact the variety of approaches, approximations and parametrizations adopted by different authors lead to results which are often discordant.

### 4.2 Clumping factor and recombination time

Following Hernquist & Springel (2003), we adopt a simple statistical description of the reionization process. Noting that each ionizing photon emitted is absorbed either by a newly ionized or by a recombining hydrogen atom, the filling factor *Q*(*t*) of regions of ionized hydrogen at any given time *t* follows from subtracting the total number of recombinations per atom from the total number of ionizing photons per hydrogen atom emitted before *t*. In this way, it is possible to derive a simple differential equation describing the transition from the neutral to the completely reionized universe (Madau, Haardt & Rees 1999):

*t*

_{rec}, is where α

_{B}is the recombination coefficient and χ is the relative abundance of helium with respect to hydrogen. The dimensionless clumping factor

*C*of hydrogen is given by It is evident from equation (7) that the recombination time is a function of redshift. Assuming a gas temperature of 10

^{4}K above which line cooling by atomic hydrogen becomes efficient, the previous relation can be written as (Madau et al. 1999).

The left-hand panel in Fig. 6 shows the redshift evolution of the clumping factor *C*, which has been estimated in our SPH simulations to be

*C*due to their high density, the sum extends over all gas particles whose density ρ

_{i}is smaller than a given threshold. In order to bracket realistic values for the threshold, we show the results (by shaded regions) assuming a maximum overdensity in the range between 100 and 500. We tested that the following analysis depends only slightly on the exact threshold.

In general, the growth of *C* with time is quite regular for all models. For this reason, we decided to fit its behaviour assuming a parabolic relation, *C*=*a*_{0}+*a*_{1}*z*+*a*_{2}*z*^{2}. The resulting best-fitting parameters *a _{i}*, obtained using

*C*the central values of the shaded regions with the corresponding dispersion, are reported in Table 1. Note that the relation for the ΛCDM-RSI model is valid only for

*z*< 22, because of the strong flattening at higher redshifts. In line with our previous results, we find that gas clumping decreases from SUGRA to ΛCDM, SUGRA-RSI and ΛCDM-RSI. The dominant effect is caused by the shape of primordial power spectrum. In fact, the differences between the two RSI models are very small at all redshifts, where

*C*is always between 1 (corresponding to an almost perfectly smooth gas distribution) and 3. For the ΛCDM model instead, the clumping factor starts to significantly deviate from 1 at

*z*≃ 25, reaching 4 at

*z*∼ 18. As mentioned before, the gas distribution appears much more clumpy in SUGRA, where

*C*starts to grow already from

*z*≃ 30, reaching

*C*∼ 6 at

*z*≃ 19. Note that our values for

*C*in ΛCDM are larger by 50 per cent compared to those obtained for the same model by Springel & Hernquist (2003) (see the open squares and the dashed line in the plot). This discrepancy is caused both by the higher resolution of our simulation and by the different treatment of the gas component in the numerical code. Even if star formation is not included like in Springel & Hernquist (2003), the fluid quantities (density, temperature and chemical abundances) are accurately and self-consistently evolved in our simulations, in particular also the cooling by molecular hydrogen. Thus, they are expected to produce a more realistic gas distribution in the diffuse intergalactic medium.

Model | a_{0} | a_{1} | a_{2} |

SUGRA | 30.379 | −1.858 | 0.030 |

ΛCDM | 17.295 | −1.093 | 0.019 |

SUGRA-RSI | 9.592 | −0.585 | 0.010 |

ΛCDM-RSI | 20.086 | −1.655 | 0.037 |

Model | a_{0} | a_{1} | a_{2} |

SUGRA | 30.379 | −1.858 | 0.030 |

ΛCDM | 17.295 | −1.093 | 0.019 |

SUGRA-RSI | 9.592 | −0.585 | 0.010 |

ΛCDM-RSI | 20.086 | −1.655 | 0.037 |

Knowing the clumpiness factor *C*, it is easy to estimate the recombination time *t*_{rec} from equation (9). The corresponding results are plotted in the right-hand panel of Fig. 6. Of course, the maximum is obtained for the model with the lowest *C*, that is, ΛCDM-RSI, in which the recombination process requires the long time of almost 40 Myr at *z*≃ 18. The values for *t*_{rec} slightly decrease to about 30 Myr for SUGRA-RSI and to 20 Myr for the ΛCDM at the same redshift, and to less than 10 Myr at *z*≃ 19 for the SUGRA model. The differences in both *C* and *t*_{rec} are again spread over almost an order of magnitude because they are integrated observables like *n*_{clouds}.

## 5 IMPLICATIONS FOR REIONIZATION

It is interesting to draw consequences from our results regarding the phenomenology of the overall reionization process, checking, in particular, how these scenarios might be constrained by data on the reionized optical depth. Even if more accurate and quantitative estimates would require detailed simulations including radiative transfer, we can reliably discuss the problem of reionization based on preceding work.

We first need a relation between the reionization epoch and the total Thomson optical depth τ, which can be derived from the CMB data. To compare the *WMAP* results to the predictions for the cosmological model considered here, we compute the redshift evolution of τ, adopting a simple model which assumes that complete reionization occurs instantaneously at some redshift *z*. The differential Thomson optical depth for complete ionization is given by

*c*is the speed of light and

*H*(

*z*) is the Hubble parameter at

*z*. The mean electron number density is given in terms of the ionization fraction

*x*

_{e}(

*z*) by

*n*(

*z*) =

*n*(

*z*= 0)(1 +

*z*)

^{3}

*x*

_{e}(

*z*). Note that in the previous equations the only dependence on the cosmological model appears in

*H*(

*z*) through the Hubble constant and the expansion rates. For this reason, the optical depth is completely independent of the spectral parameters (such as the power spectrum index

*n*), and the RSI predictions agree exactly with the corresponding

*n*= 1 model.

The integrated optical depth, computed assuming *x*_{e}= 1 always, is shown in Fig. 7 for both he ΛCDM and the SUGRA with *w*_{0}=−0.85. The horizontal dotted lines represent the observational estimate (with its 1σ error bars) obtained from the first-year *WMAP* data (Kogut et al. 2003), while the shaded region shows the more recent values from the three-year *WMAP* data (Page et al. 2006).2 When compared to standard ΛCDM, SUGRA has a lower optical depth, with a difference of about Δτ= 0.02 at *z*≃ 17 and Δτ= 0.01 at *z*≃ 10. As expected, this implies an earlier reionization epoch, causing tension between the present model and the measured reionization optical depth (Spergel et al. 2006), in addition to that regarding other cosmological probes (see e.g. Knop et al. 2003; Spergel et al. 2003; Riess et al. 2004). Although a more accurate estimate would require a realistic description of the reionization history, which in turn strongly depends on the free parameters of the adopted modelling, we can then assume that the reionization in both models (ΛCDM and SUGRA) must occur and conclude at early epoch, ranging between *z*≃ 15 and ≃10, to agree with the measured optical depth within 1σ. In general, these results demonstrate that the records of the high-redshift expansion rate in the physics of the primordial gas clouds may be used to constrain the dark energy from reionization data. Mainini, Colombo & Bonometto (2003) came to similar conclusions without considering modifications to the physics of gas clouds.

We can now study in more detail how the reionization process is going on in the redshift range covered by our simulations. First, we compute the total number of ionizing photons *n*_{ion}. For this goal, we can follow the usual ‘one star per halo’ assumption for minihaloes (see e.g. Yoshida et al. 2003a,b) and safely use the cloud number as a good approximation for the number of very massive stars: in fact strong radiative feedback disfavours the formation of multiple stars inside each primordial gas cloud. We then set the mass of a Population III star to be 500 M_{⊙} and we use the tables in Schaerer (2002) to estimate its lifetime (∼1.9 Myr) and the mean ionizing flux along the evolutionary track (∼6.8 × 10^{50} photons s^{−1}). We also assume an optimistic constant photon escape fraction of unity, which, however, can be considered realistic when the gas distribution is reasonably smooth (Oh et al. 2001). Detailed calculations by Kitayama et al. (2004) indeed find that such a large escape fraction is plausible for small-mass haloes. In order to complete the reionization process, the resulting values for *n*_{ion} must be at least as large as the total number of hydrogen atoms in the simulation volume, which in our case corresponds to *n*_{tot}∼ 4 × 10^{66}. In Fig. 8, the dashed lines show the redshift evolution for *n*_{ion} for the different models. We note that only the SUGRA model is able to produce a sufficient number of ionizing photons in the redshift range considered by our simulations, while for all remaining models the ratio *n*_{ion}/*n*_{tot} is smaller than unity. In particular, the two RSI models at the final simulation output reach only the 10–15 per cent of the required photons.

Secondly, we have to consider that, as discussed before, recombination is counteracting the reionization process. Using the values for the clumping factor and recombination time accurately derived from the simulations (and shown in Fig. 6), we can easily compute the cumulative number of recombined atoms, numerically solving equation (6). The results are also displayed in Fig. 8, where the shaded regions correspond to the uncertainties for *C* and *t*_{rec} discussed in the previous section. The plot is suggesting that the recombination process is rapid enough to significantly affect the fraction of ionized atoms. At the final simulation output, we find that the ratio between recombined and ionized atoms is about 75 per cent for the SUGRA and ΛCDM, and about 25 for the two RSI models.

Our analysis suggests that all models considered here, with the exception of SUGRA, are producing an insufficient number of collapsed objects to achieve complete reionization in the redshift interval covered by the simulations. This is true not only for both RSI models, whose suppression of small-scale power prevents the reionization, but also for the standard ΛCDM model. These conclusions are in qualitative agreement with the previous analysis performed by Yoshida et al. (2003b), who, considering a simple model based on the formation of first stars in minihaloes by molecular cooling, demonstrated that complete reionization is reached when a density of 50–100 very massive stars per comoving Mpc^{3} is turned on per ∼1 recombination time (see also Yoshida et al. 2003a,c). We note that for these models the results are not conflicting with the three-year *WMAP* data which suggest a later reionization epoch compared to the first-year analysis. Allowing the simulations to evolve up to *z*≃ 10 would certainly increase both the number of gas clouds and the clumpiness, completing the reionization process.

Finally, we note that the formation and evolution of cosmic structures in SUGRA are so fast that the corresponding number of gas clouds becomes large and *t*_{rec} becomes low that a global reionization is expected at redshifts even higher than *z*= 20. As mentioned already, this anticipation of the reionization epoch could be in conflict with observations (see e.g. the contour levels in fig. 3 of Spergel et al. 2006), indicating the potential constraining power of the dark energy records in the formation of primordial gas clouds which we pointed out in this work.

## 6 CONCLUSIONS

We presented hydrodynamic *N*-body simulations of the formation of primordial gas clouds on scales of tens of kpc in a variety of cosmological models, characterized by different dynamics in the dark energy component. Our main results are that the records of the modified expansion rate are well evident in the population and the clumpiness of such clouds. Cosmological models with the same power spectrum normalization at present show earlier cloud formation if the dynamics of the dark energy is enhanced, represented by an equation-of-state parameter *w* > −1 as in quintessence models.

Within dark energy models compatible with the present data on the CMB and large-scale structure, the difference in the integral population of clouds may vary by up to an order of magnitude, as a consequence of the different differential efficiency for structure formation. This is consistent with earlier results indicating a higher concentration in DM haloes under similar conditions (Dolag et al. 2004).

Since abundance and clumpiness of structures are directly related to the amount of primordial power on the corresponding scales, we varied the shape of the primordial power spectrum by an RSI reducing power on small scales within the confidence level of the three-year *WMAP* data (Spergel et al. 2006). As expected, we find that the extra population and clumpiness of clouds produced by a higher dark energy abundance compared with its level today might be mitigated if the primordial spectral index is running, decreasing the power on small scales.

Adopting a simple picture for the reionization process, we derived consequences for the reionization itself, leading to an earlier beginning of the reionization process in models where cloud formation starts earlier. Based on these results, we are able to identify possible tension between the *WMAP* data on the reionization optical depth and cosmological models whose dark energy is as dynamical as in SUGRA quintessence models.

Our results demonstrate that the effects on cosmological structure formation from a modified expansion history through different dark energy models must be traced back to the formation of the first clouds. In turn, this means that constraints on the dark energy, and in particular its abundance at high redshifts, may be obtained by forthcoming experiments aiming at measuring the abundance and the clumpiness of primordial gas clouds. In particular, the Atacama Large Millimeter Array (ALMA3) and the Mileura Widefield Array (MWA4) will probe the early stages of structure formation, say, between *z*= 6 and 10, where the reionization in progress should keep a record of the population of reionizing primordial gas clouds.

## Footnotes

*WMAP*data are suggesting a lower value not only for τ, but also for the power spectrum normalization σ

_{8}: these two effects nearly cancel in terms of early structure formation, that is, the inferred small τ does not really allow slow reionization for a given σ

_{8}.

*WMAP*: since

*x*

_{e}vanishes, the curves should indeed display a horizontal floor.

Computations were performed on the IBM-SP4/5 at CINECA, Bologna, with CPU time assigned under an INAF-CINECA grant, and the IBM-SP4 at the computer centre of the Max Planck Society with CPU time assigned to the Max-Planck-Institut für Astrophysik. We are grateful to C. Gheller for his assistance. We acknowledge useful discussions with Benedetta Ciardi, Bepi Tormen and Licia Verde.

## REFERENCES

*et al.*,

*et al.*,