Inhomogeneous reionization of the intergalactic medium regulated by radiative and stellar feedbacks

Abstract

We study the inhomogeneous reionization in a critical density CDM universe resulting from stellar sources, including Population III objects. The spatial distribution of the sources is obtained from high-resolution numerical N-body simulations. We calculate the source properties, taking into account a self-consistent treatment of both radiative (i.e. ionizing and H₂-photodissociating photons) and stellar (i.e. SN explosions) feedbacks regulated by massive stars. This allows us to describe the topology of the ionized and dissociated regions at various cosmic epochs, and to derive the evolution of H, He and H₂ filling factors, soft UV background, cosmic star formation rate and the final fate of ionizing objects. The main results are: (i) galaxies reionize the intergalactic medium by z ≈ 10 (with some uncertainty related to the gas clumping factor), whereas H₂ is completely dissociated already by z ≈ 25; (ii) reionization is mostly caused by the relatively massive objects which collapse via H line cooling, while objects the formation of which relies on H₂ cooling alone are insufficient for this purpose; (iii) the diffuse soft UV background is the major source of radiative feedback effects for z ≤ 15; at higher z direct flux from neighbouring objects dominates; (iv) the match of the calculated cosmic star formation history with that observed at lower redshifts suggests that the conversion efficiency of baryons into stars is ≈1 per cent; (v) we find that a very large population of dark objects which failed to form stars is present by z ≈ 8. We discuss and compare our results with similar previous studies.

galaxies: formation, intergalactic medium, cosmology: theory

1 Introduction

At z≈1100 the intergalactic medium (IGM) is expected to recombine and remain neutral until the first sources of ionizing radiation form and reionize it. The application of the Gunn-Peterson (1965) test to QSOs absorption spectra suggests that the H i reionization is complete by z≈5. Recently, Lyα emitters have been seen up to a redshift of z± 6.68 (Chen, Lanzetta & Pascarelle 1999). This indicates that the IGM was already ionized by then, but the epoch of complete reionization has yet to be firmly established. Several authors (Shapiro, Giroux & Babul 1994; Madau, Haardt & Rees 1998, and references therein) have claimed that the known population of quasars [showing a pronounced cut-off at z≳3 (Warren, Hewett & Osmer 1994; Schmidt, Schneider & Gunn 1995; Shaver et al. 1996)] and galaxies provides ≈10 times fewer ionizing photons than are necessary to keep the observed IGM ionization level. Thus, additional sources of ionizing photons are required at high redshift. Although models have been proposed in which the UV photons responsible for the IGM reionization may be emitted by a cosmological distribution of decaying dark matter particles, such as neutrinos (see, e.g., Scott, Rees & Sciama 1991), the most promising sources are early galaxies and quasars.

Recent observational evidence suggests the existence of an early population of pregalactic objects which could have contributed to the reionization and metal enrichment of the IGM. Metals have been clearly detected in Lyα forest clouds with column densities low enough to be identified with a truly diffuse IGM (Cowie et al. 1995; Tytler et al. 1995; Lu et al. 1998; Cowie & Songaila 1998), although discrepant metallicities have been inferred. Naively, one can estimate that the amount of heavy elements associated with the number of photons required to ionize every baryon in the Universe a few times corresponds to an IGM metallicity Z≈3×10⁻³ Z_⊙. This figure is roughly consistent with the value obtained by Cowie & Songaila, but 10 times larger than the one derived by Lu et al. The reasons for this discrepancy are not yet clear. Whatever the correct level of metallicity is, these heavy elements must have been produced by objects hosting the sites of the earliest star formation activity in the Universe. The presence of a soft, stellar component of the UV background deduced from studies of the [Si/C] abundance ratios in low column density absorption systems (Giroux & Shull 1997; Savaglio et al. 1997), also supports this view. The question remains concerning the transport mechanism from the production regions, which are presumably associated with high peaks of the density fluctuation field, into the very diffuse medium probed by absorption line experiments. It is currently unclear if the mechanical energy input associated with the same massive stars producing the metals (i.e. supernovae) is sufficient to drive the metals far enough from their production site, or if additional mechanisms need be invoked. Although recently Gnedin (1998), using high-resolution cosmological simulations which include a coarse modelling of a two-phase interstellar medium, has suggested that merging is the dominant mechanism for transporting heavy elements from primeval galaxies into the IGM, direct ejection of the metal-enriched interstellar gas by supernovae or stellar winds may still play an important role.

A contribution to the ionizing flux at high redshift from an early population of quasars may provide the hard component required for He ii reionization. High-resolution spectra of a quasar at z≈3 have shown large fluctuations of the H i to He ii optical depth ratio (Reimers et al. 1997), interpreted as evidence for a patchy He ii ionization in the IGM. However, it is unclear if these fluctuations could rather be caused by statistical fluctuations of the IGM density or of the ionizing background flux (Miralda-Escudé 1998; Miralda-Escudé, Haehnelt & Rees 2000). Analytical results from Haiman & Loeb (1998) show that He ii reionization due to sources with quasar-like spectra should occur at approximately the same redshift as hydrogen reionization; if the predicted He ii reionization redshift is confirmed by future measurements, this would severely constrain the number of high-redshift quasars predicted by theory.

The study of the IGM reionization due to these primeval sources has been tackled by several authors, both via semi-analytical (Shapiro 1986; Fukugita & Kawasaki 1994; Tegmark, Silk & Blanchard 1994; Haiman & Loeb 1998; Miralda-Escudé et al. 2000; Valageas & Silk 1999, hereafter VS) and numerical (Gnedin & Ostriker 1997, hereafter GO; Norman, Paschos & Abel 1998; Razoumov & Scott 1999) approaches. One of the main difficulties of the problem is the treatment of cosmological radiative transfer. Only recently have numerical approaches been attempted (Norman et al. 1998; Abel, Norman & Madau 1999), and these still require severe approximations for the radiative transfer equation. Clearly, a numerical approach offers the advantage of a more detailed treatment of features such as the IGM clumpiness (GO; Shapiro, Raga & Mellema 1998) or the source spatial distribution. Conversely, semi-analytical approaches are more flexible and allow a more thorough exploration of the parameter space.

A proper treatment of the reionization problem must take into account the interplay of galaxy formation with the reionization process itself. Basically, two types of feedback, radiative and stellar, can be at work. In the standard cosmological hierarchical scenario for structure formation, the objects which form first are predicted to have masses corresponding to virial temperatures T_vir<10⁴ K. Once the gas has virialized in the potential wells of dark matter haloes, additional cooling is required to further collapse the gas and form stars. For a gas of primordial composition at such low temperatures the main coolant is molecular hydrogen (Peebles & Dicke 1968; Shapiro 1992; Haiman, Rees & Loeb 1996b; Abel et al. 1997; Tegmark et al. 1997; Ferrara 1998). We define Population (III) (Pop III) objects as those for which H₂ cooling is required for collapse. After a H₂ molecule has been rotationally or vibrationally excited through a collision with an H atom or another H₂ molecule, a radiative de-excitation leads to cooling of the gas. This mechanism has been studied in great detail by several authors (Lepp & Shull 1984; Hollenbach & McKee 1989; Martin, Schwarz & Mandy 1996; Galli & Palla 1998) and has been applied to the study of the formation of primordial low-mass objects (Haiman, Thoul & Loeb 1996a; Tegmark et al. 1997; Abel et al. 1998). Primordial H₂ forms with a fractional abundance of ≈10⁻⁷ at redshifts ≳400 via the formation channel. At redshifts ≲110, when the Cosmic Microwave Background (CMB) radiation intensity becomes weak enough to allow for significant formation of H⁻ ions, a primordial fraction of ƒ_H2≈2×10⁻⁶ (Shapiro 1992; Anninos & Norman 1996) is produced for model universes that satisfy the standard primordial nucleosynthesis constrain Ω_bh²=0.0125 (Copi, Schramm & Turner 1995), where Ω_b is the baryon density parameter and H₀=100 h km s⁻¹ Mpc⁻¹ is the Hubble constant. This primordial fraction is usually lower than that required for the formation of Pop III objects, but during the collapse phase the molecular hydrogen content can reach high enough values to trigger star formation. On the other hand, objects with virial temperatures (or masses) above that required for the hydrogen Lyα line cooling to be efficient do not rely on H₂ cooling to ignite internal star formation. As the first stars form, their photons in the energy range 11.26–13.6 eV are able to penetrate the gas and photodissociate H₂ molecules both in the IGM and in the nearest collapsing structures, if they can propagate that far from their source. Thus the existence of a UV flux below the Lyman limit due to primordial objects, capable of dissociating the H₂, could strongly influence subsequent small-structure formation. Haiman, Rees & Loeb (1997, hereafter HRL), for example, have argued that Pop III objects could depress the H₂ abundance in neighbour collapsing clouds, due to their UV photodissociating radiation, thus inhibiting subsequent formation of low-mass structures. On the other hand, Ciardi, Ferrara & Abel (2000, hereafter CFA) have shown that the 'soft-UV background’ (SUVB) produced by Pop III objects is well below the threshold required for negative feedback to be effective earlier than z≈20. In principle, the collapse of larger mass objects can also be influenced by an ionizing background, as gas in haloes with a circular velocity lower than the sound speed of ionized gas may be prevented from collapsing due to pressure support in the gravitational potential. Thoul & Weinberg (1996) (see also Babul & Rees 1992) have, however, shown that the collapse is only delayed by this process. We will refer to this complex network of processes as ‘radiative feedback’.

The other feedback mechanism is related to massive stars, and for this reason we will call it 'stellar feedback’. Once star formation begins in the central regions of collapsed objects, it may strongly influence their evolution via the effects of mass and energy deposition due to massive stars through winds and supernova explosions. These processes may induce two essentially different phenomena. Low-mass objects are characterized by shallow potential wells in which the baryons are only loosely bound and a relatively small energy injection may be sufficient to expel the entire gas content back into the IGM, i.e. a blowaway, thus quenching star formation. Larger objects may instead be able to at least partially retain their baryons, although a substantial fraction of the latter are lost in an outflow, i.e. a blowout (Ciardi & Ferrara 1997; Mac Low & Ferrara 1999, hereafter MF; Ferrara & Tolstoy 2000, hereafter FT). However, even in this case the outflow induces a decrease of the star formation rate due to the global heating and loss of the galactic ISM. Omukai & Nishi (1999) have pointed out that the ionizing radiation of the first stars formed in Pop III objects can also produce an abrupt interruption of the star formation by dissociating the internal H₂ content. The effect of this feedback is very similar to the blowaway, as both processes are regulated by massive stars.

Our aim here is to study and describe in detail the reionization process of the IGM. As the ionizing sources are not spatially homogeneously distributed, the pattern of the ionized regions is not uniform, i.e. inhomogeneous reionization takes place. Thus a crucial ingredient of such calculations is the realistic description of the ionized region topology in the Universe. This can be achieved only by numerical simulations which track the formation and merging of dark matter haloes. For our purposes, these simulations must reach the highest possible mass resolution in order to unambiguously identify the very small Pop III objects initiating the reionization process. The second fundamental ingredient is a proper treatment of radiative and stellar feedback processes on which a great deal of previous work, either by ourselves or other groups, has been accumulating, ranging from semi-analytical studies to numerical models and hydrodynamical simulations, particularly with regard to stellar feedback. Improving on such results, we have been able to effectively model and include in the computations a wide number of effects governed by ionizing/dissociating radiation and massive star energy injection which will be presented and discussed in detail. The delicate link between the properties of dark matter haloes, which can be reliably derived from N-body simulations, and those describing their baryonic counterparts is represented by the prescription used to model star formation. The advantages that the proposed approach brings along should become evident by reading the paper: it allows us to derive self-consistent inhomogeneous reionization models, and to infer the properties of the early galaxies producing it. In addition, its predictive power allow us to design specific tests that could validate the results.

The plan of the paper is the following. In Section 2 we present the N-body numerical simulations, and in Section 3 we describe our assumptions about star formation. Then we turn to the discussion of radiative (Section 4) and stellar (Section 5) feedbacks, and a brief sketch summarizing the possible galactic evolutionary tracks due to such processes is presented in Section 6. Results are given in Section 7, and further discussed in Section 9. Section 8 presents a comparison with previous work. A short summary concludes the paper.

2 Numerical simulations and halo selection

We have simulated structure formation within a periodic cube of comoving length L± 2.55 h⁻¹ Mpc for a critical density cold dark matter model (Ω₀=1, h± 0.5 with σ₈=0.6 at z± 0). Our choice for the normalization (σ₈) corresponds roughly to those inferred from the present-day cluster abundance (see, e.g., Eke, Cole & Frenk 1996 and Governato et al., in preparation). We remind the reader that once the present-day value of σ₈ for the SCDM run is selected, the redshift epoch of all other outputs with lower values of σ₈ is uniquely specified. For example, for the case where the present-day normalization is chosen to be σ₈=0.6, σ₈=0.3 output corresponds to the z± 1 epoch. It is useful to note that the rms linear overdensity of a sphere with a mass equal to the box is 2.1/(1+z). At early times when the amplitude of fluctuations are small for scales larger than the box size the simulation box should be large enough to be representative of the world model as a whole. However, as the scale of non-linear fluctuations grows, there will come a time when the simulation volume will not be large enough to provide a proper average. By redshift z± 8, for example, the rms linear overdensity of a sphere with the mass of the entire simulation is 0.23, and by then the mass spectrum of dark haloes in the simulation box may not match the real mass spectrum very well — particularly at the high-mass end.

The initial conditions for the simulations were set up at z± 100 and integrated in time to z± 8.3 (with intermediate outputs at z± 29.6, 25.3, 22.1, 19.8, 18.0, 16.5, 15.4, 14.3, 11.4 and 10.9). We used a transfer function for CDM calculated with cmbfast and a baryonic fraction Ω_b=0.06 (Copi et al. 1995). The simulation, which uses 256³ particles (about 17 million), was computed with AP3M (Pearce & Couchman 1997), a parallel P3M code. A cubic spline force softening of 0.05 h⁻¹ kpc Plummer equivalent was used so that the overall structure of the haloes could be resolved. The good mass (particle mass is ≈5×10⁵ M_⊙) and force resolution allows us to study in detail the evolution of all structures likely to host the formation of primordial stars.

In numerical simulations, haloes can be identified using a variety of schemes. Of these, we have chosen one that is available in the public domain: FOF¹ (Davis et al. 1985). In this scheme, all particle pairs separated by less than b times the mean interparticle separation are linked together. Sets of mutually linked particles form groups that are then identified as dark matter haloes. Other halo finders that are often used in literature to find virialized haloes are HOP (Eisenstein & Hut 1998), the 'spherical overdensity algorithm’ or SO, which finds spherically averaged haloes above a given overdensity (Lacey & Cole 1994), and the scheme recently developed by Gross et al. (1998). In the present study we adopted the linking length that Lacey & Cole deduced to identify virialized haloes with mean densities of ≈200 times the critical density at the epoch under consideration. Halo masses (and, ultimately, our findings) do not depend on the halo finder used (FOF, HOP or SO), as long as it selects haloes including all particles down to the same overdensity. Systematic deviations between different algorithms are of the order of a few per cent (see Lacey & Cole 1994 and Governato et al., in preparation). In our analysis, we consider only haloes consisting of 16 particles or more, corresponding to a mass of 8×10⁶ M_⊙.

We have compared the results of the numerical simulations for the dark matter halo distribution with the one obtained from the Press & Schechter (1974, hereafter PS) formalism, so that we can use this semi-analytical expression to calculate certain quantities, such as the soft-UV background (see Section 4.2). The power spectrum normalization is σ₈=0.6, consistent with the above numerical simulations. Although PS slightly overestimates the halo number density, the match is remarkably good over the entire range of redshifts considered.

Finally, the correlation function for galaxies at high redshift has been recently predicted (Governato et al. 1998) and then measured by Giavalisco et al. (1998), giving strong support to hierarchical models of galaxy formation. Bright galaxies at high redshift form on top of the highest peaks in the mass distribution, and so they are strongly biased with respect to the underlying dark matter distribution. In hierarchical models we expect a similar behaviour for the most massive haloes hosting Pop III stars. This measurement would give a strong constraint on the shape of the power spectrum at scales smaller than 1 Mpc, very relevant for all theories of galaxy formation. Indeed, numerical simulations have already shown a large discrepancy between the shape of the rotation curves of CDM haloes and those of real dark matter dominated galaxies (see Moore et al. 1999a) as CDM models predict dark matter haloes with steeper density profiles than those of real galaxies. We have preliminarily obtained the two-point correlation function of haloes significantly contributing to the photon production between redshifts ≈20 and 8.3. In our SCDM model the correlation function of haloes hosting Pop III objects has a (comoving) length-scale of about 0.15 h⁻¹ Mpc, slowly decreasing with redshift. This decrease is expected (Tegmark & Peebles 1998), as the number of haloes of mass sufficient to host Pop III objects increases with time. Feedback effects of the type discussed here are clearly strongly influencing the evolution of the correlation function. We plan to present the results concerning the correlation properties of these haloes in a future communication.

3 Forming the first stars

Once the gas, driven by gravitational instabilities, has been virialized in the potential well of the parent dark matter halo, further fragmentation of the gas and ignition of star formation is possible only if the gas can cool efficiently and lose pressure support. For a plasma of primordial composition at temperature T<10⁴ K, the typical virial temperature of the early bound structures, molecular hydrogen is the only efficient coolant. Thus a minimum H₂ fraction is required for a gas cloud to be able to cool in a Hubble time. As the intergalactic relic H₂ abundance falls short by at least two orders of magnitude with respect to the above value, the fate of a virialized lump depends crucially on its ability to rapidly increase its H₂ content during the collapse phase. Tegmark et al. (1997) have addressed this question in great detail by calculating the evolution of the H₂ abundance for different halo masses and initial conditions for a standard CDM cosmology. They conclude that if the prevailing conditions are such that a molecular hydrogen fraction of order of ƒ_H2≈5×10⁻⁴ is produced, then the lump will cool, fragment and eventually form stars. This criterion is met only by larger haloes, implying that for each virialization redshift there will exist some critical mass, M_crit, such that protogalaxies with total mass M>M_crit will be able to form stars and those with M<M_crit will fail (see their fig. 6 for the evolution of M_crit with the virialization redshift). In reality, even haloes with masses smaller than M_crit could eventually collapse at a later time (Haiman & Loeb 1997), with a delay increasing with decreasing baryonic mass for a fixed rms amplitude of the fluctuation. However, as we will see below, these structures will be strongly affected by various feedbacks which essentially erase their contribution to the reionization process; hence we neglect this effect in our calculations. In the absence of additional effects that could prevent or delay the collapse (such as the feedbacks discussed below), we can associate to each dark matter halo with M>M_crit a corresponding baryonic collapsed mass equal to M_b=Ω_bM. Throughout the paper we adopt a value of the baryon density parameter Ω_b=0.06 consistently with the numerical simulations; note that the set of cosmological parameters used here are the same as in Tegmark et al. (1997), thus allowing a direct use of their results.

As the gas collapses and the first stars form, stellar photons with energies in the Lyman-Werner (LW) band and above the Lyman limit, respectively, can photodissociate H₂ molecules and ionize H and He atoms in the surrounding IGM. By this process a photodissociated/ionized region will be produced around the collapsed object, whose size will depend on the source emission properties. The radiation spectrum, j(ν), adopted here is obtained from the recently revised version of the Bruzual & Charlot (1993, hereafter BC) spectrophotometric code, for a Salpeter initial mass function (IMF), a single burst mode of star formation, and a metallicity Z± 10⁻² Z_⊙. This choice is supported by recent observational results indicating that the IMF in nearby systems has a slope close the Salpeter value above a solar mass, while flattening at lower masses (see Scalo 1998). A time-independent IMF is assumed by most studies, but indirect evidence for an IMF biased towards massive stars at early times emerges (see Larson 1998, and references therein). Larson suggests that the IMF might be well represented by a universal power law at large masses, flattening below a characteristic mass whose value is changing with time. The rationale for this conclusion is that the mass-scale for star formation probably depends strongly on temperature and, since star-forming clouds were probably hotter at earlier cosmic times, the mass-scale should also have been correspondingly higher. This results in an increase of the relative number of high-mass stars formed, and consequently of the number of ionizing photons. The relevance of the Jeans mass-scale for the IMF is still subject of lively debate, and it is not clear to what extent it might be responsible for the mass distribution. In view of these uncertainties and for the sake of simplicity, we will adopt the Salpeter IMF throughout the paper.

Fig. 1 shows the spectral energy distribution of an object per solar mass of stars formed at four different burst evolutionary times t± (0,1,2,32) ×10⁷ yr. At earlier times the spectrum is rather flat between the 13.6 eV and the He ii edge, with a pronounced cut-off at higher energies; at later stages, it becomes softer with strong discontinuities at both 13.6 and 24.6 eV. The number of ionizing photons decreases considerably as the stars age.

Figure 1.

Evolution of the adopted emission spectrum as a function of photon energy, at four different evolutionary times t± (0,1,2,32) ×10⁷ yr. The specific luminosity corresponds to one solar mass of stars formed.

The total luminosity of an object is determined by the mass of stars formed. For a baryonic mass M_b=10⁵M_b,5 M_⊙, the corresponding stellar mass is

(1)

where ƒ_b=0.08 ƒ_b,8 is the fraction of virialized baryons that is able to cool and become available to form stars, and ƒ*=0.15 ƒ*_,15 is the star formation efficiency. With these assumptions the luminosity per unit frequency at the Lyman limit, j₀, as obtained from the adopted spectrum at early evolutionary times (Fig. 1), can be explicitly written in terms of M_b:

(2)

of which only a fraction j₀ƒ_esc≃9.6×10²²M_b,5ƒ_b,8ƒ*_,15ƒ_esc,20 erg s⁻¹ Hz⁻¹ (where ƒ_esc=0.2 ƒ_esc,20 is the photon escape fraction from the protogalaxy) is able to escape into the IGM. This accounts at least approximately for absorption occurring in the host galaxy.

As ƒ_b, ƒ* and ƒ_esc are the main parameters involved in the calculation, aside from the cosmological ones fixed by the simulations, it is useful to discuss them in more detail. Primordial stars can form only in the fraction of gas that has been able to cool, and the cooled-to-total baryonic mass ratio in the collapsed objects is only a few per cent. Numerical simulations (Abel et al. 1998; Bromm, Coppi & Larson 1999) have shown that ƒ_b can be ≈8 per cent, which we assume as our fiducial value. The star formation efficiency, ƒ*, is rather uncertain and dependent over the properties of the star formation environment. Even its definition is not completely unambiguous, as star formation usually takes place in the densest cores of molecular clouds and the efficiency in the cores is obviously different from the one averaged on the parent molecular cloud. We define ƒ* as the fraction of cooled (probably molecular) gas that has been converted into stars. The formation of a bound cluster system requires a star formation efficiency of nearly 50 per cent where the cloud disruption is sudden, and nearly 20 per cent where cloud disruption takes place on a longer timescale (Margulis & Lada 1983; Mathieu 1983). The highest efficiency estimated for a star-forming region is about 30 per cent for the Ophiuchi dark cloud (Wilking & Lada 1983; Lada & Wilking 1984). Lada, Evans II & Falgarone (1997) have studied the physical properties of molecular cloud cores in L1630, finding efficiencies ranging from 4 to 30 per cent for different cores. Pandey, Paliwal & Mahra (1990) have investigated the effect of variations in the IMF on the star formation efficiency in clouds of various masses. They conclude that the efficiency is lower if the most massive, and consequently most destructive, stars form earlier. On a larger scale, Planesas, Colina & Perez-Olea (1997) have studied the molecular gas properties and star formation in nearby nuclear starburst galaxies, showing the existence of giant molecular clouds (M_H2≈10⁸−10⁹ M_⊙) with associated H ii regions where the star formation process is characterized by being short-lived (<3×10⁷ yr) and with an overall gas to stars conversion less than 10 per cent of the gas mass. Giannakopoulou-Creighton, Fich & Wilson (1999) have recently determined the star formation efficiency in two M101 giant molecular clouds, finding the values 6 and >11 per cent. In spite of great efforts to derive such an important quantity, various authors obtain results uncomfortably different, even when studying the same star-forming complex (see the prototypical case of 30 Dor). Given this situation, we use the educated guess ƒ*=15 per cent. Finally, the value ƒ_esc=0.2 is an upper limit derived from observational (Leitherer et al. 1995; Hurwitz, Jelinsky & Dixon 1997) and theoretical (Dove & Shull 1994; Dove, Shull & Ferrara 1999) studies. The latter authors, in particular, have included the effects of superbubbles in the computation of ƒ_esc, finding a value ƒ_esc≈6 per cent for the case of coeval star formation history, which should be relevant to the present case. These studies concentrate on large disc galaxies like the Milky Way, and it is not clear to what extent they can be applied to the small objects populating the early Universe; also, the much lower dust-to-gas ratio could make the escaping probability higher. Again, ƒ_esc=0.2 can be regarded only as a reference value. It is worth noting that while for all processes considered here only the product ƒ_b׃* enters the model, ƒ_esc is instead an independent parameter; this reduces the number of effective free parameters to two.

The total ionizing photon rate, S_i(0), from the formed stellar cluster can be written as

(3)

where ν_l=13.6 eV and ν_u=150 eV. As seen from Fig. 1, the stellar spectrum drops off sharply above 100 eV even shortly after the burst; the choice of the upper integration limit follows from this remark. For the spectrum in Fig. 1 at t± 0, equation (3) becomes

(4)

analogously we can define the LW photon production rate, S_LW, which is typically S_LWβS_i(0), with β≈1.5 at early evolutionary times. The UV photons create a cosmological H ii region in the surrounding IGM, whose radius, R_i, can be approximated by the Strömgren (proper) radius, where n_H=8×10⁻⁶Ω_bh²(1+z)³ cm⁻³ is the IGM hydrogen number density, and α⁽²⁾ is the hydrogen recombination rate to levels ≥2. In general, R_s represents an upper limit for R_i, since the ionization front fills the time-varying Strömgren radius only at very high redshift, z≈100 (Shapiro & Giroux 1987). For our reference parameters it is

(5)

where S₄₇S_i(0)/(10⁴⁷ s⁻¹), and (1+z)₃₀=(1+z)/30. In our calculation we use the exact solution for R_i as numerically calculated in CFA, and whose analytical approximation is given by Shapiro & Giroux. Typically, R_i is about 1.5–2.0 times smaller than R_s, and cosmological expansion cannot be neglected. Clearly, these solutions assume that the ionizing object is irradiating a volume of the IGM which was previously neutral. If instead the object is located inside (or overlapping) a pre-existing ionized sphere, a detailed solution of the radiative transfer would be necessary. This occurrence can result in an underestimate of our calculated ionized volume when substantial sphere overlapping takes place, i.e. close to the reionization epoch.

Similarly to R_s, we can define the He ii Strömgren radius as R_s,He=[3S_i,He(0)/(4πn_He++n_eα_He++)]^1/3, where S_i,He(0) is the rate of ionizing photons with hν>54.4 eV, n_He++ and n_e are the He⁺⁺ and electron number density, respectively, and α_He++ is the He⁺⁺ recombination rate. For the adopted spectrum it is

(6)

We have derived n_He++ and n_e assuming that the helium and hydrogen neutral fraction is negligible. For our reference parameters and given equations (4) and (6), the radius of the sphere of doubly ionized helium is

(7)

thus R_s,He≃0.1R_s.

Given a point source that radiates S_LW photons per second in the LW bands, a good estimate of the maximum radius of the H₂ photodissociated sphere is the distance at which the (optically thin) photodissociation time becomes longer than the Hubble time:

(8)

where S_LW,47S_LW/(10⁴⁷ s⁻¹). CFA also considered the time-dependent evolution of both the ionization and dissociation fronts; here we use the exact value for R_i and the equilibrium values given by the equations above for both R_i,He and R_d. The latter values represent upper limits to the extent of the spheres as emphasized above: ionized regions at high redshift fill only partially the corresponding Strömgren spheres, while for the dissociated regions the radius R_d is completely filled, though on a relatively long time-scale. For example, an object born at z± 30 will have a fully developed surrounding dissociated sphere by z≈20.

To each dark matter halo in which stars can form we assign the ionized and dissociated spheres produced by its stellar cluster according to the above prescriptions. This allow us to derive the three-dimensional structure and topology of the ionized/dissociated regions as a function of cosmic time.

4 Radiative feedback

In addition to their local effects, the first objects will also produce UV radiation which could in principle introduce long-range feedbacks on nearby collapsing haloes. Particularly relevant is the SUVB in the LW bands, as by dissociating the H₂ it could influence the star formation history of other small objects preventing their cooling.

CFA have argued that the SUVB produced by Pop III objects is below the threshold required for the negative feedback on the subsequent galaxy formation to be effective before z≈20. At later times this feedback becomes important for objects with T_vir≲T_H=10 000 K, corresponding to a mass M_H=4.4×10⁹ M_⊙ (1+z_vir)^−1.5h⁻¹ (Padmanabhan 1993), where z_vir is the redshift of virialization, for which cooling via Lyα line radiation is possible. Protogalaxies with M>M_H are not affected by the negative feedback and are assumed to form. In principle, even for these larger objects a different type of feedback can be at work: haloes with a circular velocity lower than the sound speed of the gas in an ionized region may be prevented from collapsing due to gas pressure support in the gravitational potential. Thoul & Weinberg (1996) (see also Babul & Rees 1992) have shown that the collapse is only delayed by this process, and therefore we neglect this complication.

At higher redshift the radiative feedback can be induced by the direct dissociating flux from a nearby object. In practice, two different situations can occur: (i) the collapsing object is outside the dissociated spheres produced by pre-existent objects: then its formation could be affected only by the SUVB (J_LW,b), as by construction the direct flux (J_LW,d) can only dissociate molecular hydrogen on time-scales shorter than the Hubble time inside R_d; (ii) the collapsing object is located inside the dissociation sphere of a previously collapsed object: the actual dissociating flux in this case is essentially given by J_LW,max=(J_LW,b+J_LW,d). It is thus assumed that, given a forming Pop III object, if the incident dissociating flux (J_LW,b in the former case, J_LW,max in the latter) is higher than the minimum flux required for negative feedback (J_s), the collapse of the object is halted. This implies the existence of a population of ‘dark objects’ which were not able to produce stars and, hence, light.

4.1 Minimum flux for negative feedback

As already stated, in the absence of an external dissociating flux, an object can form if it satisfies the condition M>M_crit. In this subsection we investigate how this minimum mass is changed by the requirement that the object is able to self-shield from an external incident dissociating flux that could deplete the H₂ abundance and prevent the collapse of the gas. For this purpose, we consider the non-equilibrium multifrequency radiative transfer of an incident spectrum inside a homogeneous gas layer, and study the evolution of the following nine species: H, H⁻, H⁺, He, He⁺, He⁺⁺, H₂, and free electrons. The energy equation can be written as (Ferrara & Giallongo 1996)

(9)

where T, n and p are the gas temperature, number density and pressure, respectively. The gas is assumed of primordial composition with a helium abundance equal to n_He=0.1 n_H and total number density n. The symbol X± H, He, H₂ denotes the species considered, and i its state of ionization; obviously, i± 0 for hydrogen, and i± 0, 1 for helium; h_Xn_Xn is the relative abundance of the species X; with i_max=0, 1 for H and He, respectively, and x_Xi+1n_Xi+1n_X the fractional density of the ionization state i+1 of the element X;γ is the specific heat ratio, assumed to be constant and equal to 5/3. The function L represents the net cooling rate per unit volume and is given by the difference between heating and cooling:

(10)

The cooling term Λ (in erg cm⁻³ s⁻¹) includes collisional ionization and excitation of H, He, He⁺, recombination to H, He, He⁺, dielectronic recombination to He, free-free (Black 1981), Compton cooling (Peebles 1971) and H₂ cooling (Martin, Schwarz & Mandy 1996); H includes the heating terms due to photoionization of H, He, He⁺, H₂, photodissociation of and H₂, and H⁻ photodetachment and is given by

(11)

where ν_X indicates the ionization limit for each species, σ is the photoionization cross-section, and J_s is the incident flux, given below. The various cross-sections are given in Abel et al. (1997), apart from σ₂₁ (taken from Brown 1971) and σ₂₅ (Shapiro & Kang 1987), and they are numbered according to the nomenclature of Abel et al. (1998). We have assumed that He, H⁻ and are in equilibrium, as all reactions determining the H⁻ abundance occur on much shorter time-scales than those relevant to the H₂ chemistry; He and do not influence substantially the final results. The chemical network includes the 27 reactions listed in Abel et al. (1998). We have adopted the same rates of that paper except for k₁ (taken from Cen 1992), k₂ (Black 1981), k₁₁ and k₁₅ (Shapiro & Kang 1987). The photoionization rate γ_p is given by

(12)

where φ is the secondary ionization rate. We have adopted the ‘on the spot’ approximation in which the diffuse field photons are supposed to be absorbed close to the point where they have been generated.

To derive generally valid results for the negative feedback, we have taken the same incident spectrum presented in Fig. 1 but at a later evolutionary time (2×10⁷ yr after the burst) in order to average roughly over the different population evolutionary stages. An useful analytical approximation to the adopted spectrum is:

(13)

with β± 50 and Ehν. Here J_s,0 (units of erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹) is the parameter with respect to which we quantify the negative feedback at a given redshift.

Initially, we assume that the gas layer is homogeneous, with a mean density 〈ρ_b〉 18π² times higher than the background intergalactic medium. The fraction of free electrons is taken to have the relic value after recombination (see, e.g., Blumenthal et al. 1984); in principle, x could have changed during the virialization process: we neglect this possibility. We use equilibrium values for the He ionic abundances; H⁻ and abundances can be arbitrarily small, as their choice does not influence the results. The initial condition for the molecular abundance is slightly more delicate; we determine it as follows. For a given M_crit (the mass of an object able to collapse in the absence of radiation) we can determine the corresponding virial temperature T_vir. This represents the initial condition for the gas temperature. The initial H₂ abundance is set by the value required to satisfy the condition that the free-fall time of the object is equal to its cooling time at temperature T_vir. This assures that, if no feedback is imposed on the object, it would collapse and form stars.

Starting from these initial conditions, we have followed the chemical evolution up to a free-fall time. These conditions are such that in the absence of radiation the fraction of H₂ is high enough to allow for the collapse of the object on a time-scale comparable to the free-fall time. The effect of the radiation field is to decrease this abundance in the external regions, which therefore cannot cool. In the interior, on the other hand, where the H₂ abundance essentially remains equal to the initial one, the collapse can proceed unimpeded. For this reason the collapse conditions must be checked after a free-fall time. As an example we discuss the case of an object which virializes at z± 29.6 and to which an incident flux with spectrum given by equation (13) with J_s,0J₂₁10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ has been imposed. The critical mass for collapse, M_crit, corresponding to that redshift is M_crit=9×10⁵ M_⊙; this fixes the initial H₂ fraction as explained above. Figs 2– 4 show the corresponding hydrogen ionization fraction, x, the gas temperature, T, and the neutral molecular hydrogen fraction, ƒ_H2, after a free-fall time. For sake of comparison, we also show the cases relative to a power-law (PL) spectrum with α± 1.5 and a cut-off energy of 40 keV. Different values of the β parameter, describing the relative number of photodissociating and ionizing photons in the spectrum, are studied in order to cover a wide range of evolutionary phases of the stellar energy emission. As seen from Fig. 1, the value of β tends to increase as the stellar cluster ages. Our standard case is the BC spectrum with β± 50, appropriate to an age of 2×10⁷ yr.

Figure 2.

Evolution of the ionized atomic hydrogen fraction as a function of depth, for a Pop III at z± 29.6 and an incident flux J_s,0J₂₁×10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹. The curves are derived for the BC stellar spectrum calculated at t± 2×10⁷ yr and, for comparison, for a power-law (PL) spectrum with α± 1.5. Different values of the β parameter are shown. The horizontal line represents the initial condition.

Figure 3.

As Fig. 2 for the gas temperature.

Figure 4.

As Fig. 2 for the neutral molecular hydrogen fraction.

From Fig. 2 it is clear that the ionization fraction is mostly determined by photons with energies above the Lyman limit; however, in the cases with β± 1, as the final abundance of H₂ is increased with respect to the initial one (see Fig. 4), electrons are depleted by the H⁻ molecular hydrogen formation channel, with respect to the cases with β± 50. The high-energy photons present in the PL spectrum penetrate deeper inside the layer, and consequently the initial ionization is enhanced over a larger depth than for the BC spectrum, though with a lower ionization level; this feature is typical of hard ionizing spectra. For the same reason, a PL spectrum heats the cloud to higher temperatures and the plateau extends into deeper regions (Fig. 3). Lower temperatures and higher ƒ_H2 are obtained for smaller values of β (see Fig. 3), when the H₂ cooling is more efficient. It is clear from Fig. 3 that, although in the external regions of the layer heating is dominant due to photoelectric effect, in the inner regions, where the medium becomes optically thick to LW photons, after a free-fall time the temperature has decreased by about 30 per cent due to H₂ cooling. In Fig. 4 the ƒ_H2 profile is shown. In the external regions, the most important mechanisms for the H₂ destruction are direct photoionization by photons with energies above 15.4 eV and LW photons dissociation. A PL spectrum with the same value of β yields a larger H₂ destruction, while if the upper PL spectrum cut-off is increased the curves shift towards larger depths due to penetrating photons. Note that when β± 1, H₂ is formed rather than destroyed due to the larger abundance of free electrons and paucity of LW photons available. Finally, the ƒ_H2 smoothly approaches the initial value in the internal regions due to the loss of free electrons (see Fig. 2).

We would like to comment briefly on the differences between the above curves for the H₂ abundance and the analogous ones derived by HRL. These are mainly due to (i) a different incident flux, and (ii) their assumption of chemical and ionization equilibrium at a constant temperature, which could not be reached within a free-fall time. As they also note in a previous paper (Haiman et al. 1996b), this might lead to an overestimate of the H₂ abundance. In addition, they used the Lepp & Shull (1984) cooling function, which, at the low densities of interest here, is higher than the one we have adopted (for a comparison see Galli & Palla 1998).

As shown above, after a free-fall time the temperature in the inner regions of the protogalactic cloud has considerably decreased due to H₂ cooling. It is normally assumed that this process would eventually lead to the collapse, fragmentation and star formation. These regions are therefore not affected by the imposed external flux as they can efficiently self-shield. The gas in the outer regions, up to the shielding radius, R_sh(J_s,0), is instead globally heated by radiation and fails to collapse. On this basis, we can then define the minimum total mass required for an object to self-shield from an external flux of intensity J_s,0 at the Lyman limit, as where is the mean dark halo matter density. In practice, R_sh is defined as the point where the temperature profile has increased by 0.01 per cent from the inner flat curve (see Fig. 3). Values of M_sh for different values of J_s,0 have been obtained at various redshifts and will be then used in the calculations presented below; these curves are shown in Fig. 5. We point out that the results depend only very weakly on the detailed shape of the spectrum as long as this is produced by stellar sources with a relatively soft spectrum; differences up to a factor of a few can be found if instead a hard spectral component is present (Haiman et al. 1999). In this paper we restrict our analysis to the first case. Protogalaxies with masses above M_H, for which cooling is predominantly contributed by Lyα line cooling, will not be affected by the negative feedback studied here: these objects lie on the upper dashed portions of the curves in Fig. 5. The collapse of very small objects with mass <M_crit is, on the other hand, made impossible by the cooling time being longer than the Hubble time (lower dashed portion of the curves). Thus the only mass range in which negative feedback is important (solid portion) lies approximately in the range 10⁶−10⁸ M_⊙, depending on redshift. In order for the negative feedback to be effective, fluxes of the order of 10⁻²⁴−10⁻²³ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ are required. In turn, by using equation (4), we see that these fluxes are produced by a Pop III object with baryonic mass M_b,5 at distances closer than for the two above flux values, respectively, while the SUVB can reach an intensity in the above range only after z≈15 (see Section 7.2). This suggests that at high z negative feedback is driven primarily by the direct irradiation from neighbour objects in regions of intense clustering, while only for z≲15 does the SUVB becomes dominant.

Figure 5.

Minimum total mass for self-shielding from an external incident flux with intensity J_s,0 at the Lyman limit. The curves are for different redshift: from top to bottom z± 8.3, 11.4, 15.4, 18.0, 22.1 and 29.6. Circles show the values of M_H (open) and M_crit (filled). Radiative feedback works only in the solid portions of the curves.

4.2 Maximum incident flux

We define the maximum dissociating flux felt by a forming Pop III object as the sum J_LW,max=(J_LW,b+J_LW,d), where J_LW,b is the SUVB, and J_LW,d is the direct flux produced by a nearby luminous object. J_LW,d is defined as

(14)

r is the distance between the forming object and the source of the dissociating flux. The value of j₀ is that given in equation (2), and the value of β± 1.5 is adopted for the reasons explained in Section 3.

We now derive the intensity of the SUVB. As the LW range is very narrow (≈2 eV), we consider the average flux at the central frequency of the band, hν₀=12.45 eV. To properly calculate the intensity of the SUVB we must consider the intergalactic H₂ attenuation, including the effects of cosmological expansion, and treat in detail the radiative transfer through LW H₂ lines. LW lines are optically thin, with τ_iN_H2,iσ_i≈0.05, essentially for all the lines. This implies that as J_LW is redshifted due to cosmological expansion through LW lines; it is attenuated by each line by a factor e^−τi. Abgrall & Roueff (1989) have included in their study of classic H₂ PDRs more than 1000 LW lines. As in our case, H₂ formation, which leaves the molecule in excited roto/vibrational levels, is negligible, a smaller number of lines ≈70 — involving the ground state only — needs to be considered. Globally, J_LW is attenuated by a factor e^−τH2, where and i runs up to the 71 lines considered from the ground roto/vibrational states; we obtain τ_H2≲3, depending on the number of lines encountered, and thus on the photon energy. The intensity of the background is

(15)

(16)

where z_on=29.6 corresponds to the redshift at which the first objects start to form, and ν′=ν₀(1+z′)/(1+z); ε(ν′,z′) [erg cm⁻³ s⁻¹ Hz⁻¹] is the proper emissivity due to sources with M_min≤M≤M_max; A_H,H2 takes into account the H and H₂ line absorption in the LW band and the H absorption above the Lyman limit; the corresponding curves of growth are taken from Federman, Glassgold & Kwan (1979). j(ν′) is taken from the BC spectrum as described in Section 3 with intensity given by equation (2), and N(M, z′) is the source number density in the mass interval dM derived from the PS formalism. We use the PS distribution, which closely reproduces the mass distribution of our simulations (see Section 2), in order to account for haloes larger than those contained in our simulations; this correction is of the order of a few per cent. Also, it would be very difficult to obtain a reliable estimate of the SUVB directly from the simulations due to the uncertain treatment of photons leaking from the box. The upper cut-off is M_max=10¹⁴ M_⊙, the maximum statistically significant mass down to z± 8.3, the lowest redshift of the numerical simulation. M_min is the mass of the smallest halo that has been able to collapse at the previous redshift step; we are thus assuming that all objects with masses greater than M_min will contribute to the SUVB. This calculation of J_LW,b is correct as long as the filling factor of the ionized and dissociated regions is small. As the volume filling factor of these regions grows with time, their attenuation of the radiation emitted by the sources is diminished and eventually, when the dissociation/ionization process is complete, it vanishes. To take roughly into account these complications, we assume that only a fraction 1−ƒ_i of the IGM gas, where ƒ_i is the filling factor of either ionized or dissociated regions, contributes to the attenuation. We note that taking into account this effect, together with the determination of the lower limit M_min of the integral governed by the negative feedback, implies that equations (15) and (16) cannot be solved directly but must be computed iteratively as we derive the evolution of the reionization process.

5 Stellar feedback

Once the first stars have formed in the host protogalaxy, they can deeply influence the subsequent star formation process through the effects of mass and energy deposition due to winds and supernova explosions. While low-mass objects may experience a blowaway, expelling their entire gas content into the IGM and quenching star formation, larger objects may instead be able to at least partially retain their baryons. However, even in this case the blowout induces a decrease of the star formation rate due to the global heating and loss of the galactic ISM. These two regimes are separated by a critical mass, M_by, to be calculated in the following.

For the relatively small objects present during the reionization epoch 30≳z≳10 the importance of these stellar feedbacks can hardly be overlooked. To understand the role of stellar feedback, let us consider a collapsed object with mass lower than M_by. Then the star formation is suddenly halted as the entire gas content is removed. In this case, the ionizing photon production will last only for a time interval of order t_OB≈10⁷ yr, the mean lifetime of the massive stars produced initially; after this, the ionized gas around the source will start to recombine at a fast rate as a result of the highly efficient high-z Compton cooling, rapidly decreasing the temperature inside the ionized region. Due to the short lifetime and recombination time-scales, these object will only produce transient H ii regions which will rapidly disappear. In fact, the recombination time-scale, when the Compton cooling is taken into account, is of order t_rec≈(1–50) ×10⁶ yr at redshift z≈30–10, respectively, and therefore much shorter than the corresponding Hubble time. The dissociating photon production will nevertheless continue for a longer time, due to the important contribution of long-lived intermediate-mass stars formed in the same initial star formation burst. This, combined with the fact that there is no efficient mechanism available to re-form the destroyed H₂ in the IGM analogous to H recombination, implies that the dissociation is not impeded by blowaway and the contribution from these small objects should be still accounted for. However, this is not necessary. In fact, after blowaway, H₂ is efficiently formed in the shocked IGM gas, cooling under non-equilibrium conditions (Ferrara 1998). The final radius of the cooled shell behind which H₂ is formed by this process is

(17)

We notice that the formation process produces an amount of H₂ roughly similar to that destroyed by radiation (Ferrara 1998). As inside R_s H₂ is re-formed, the net effect on the destruction of intergalactic H₂ of blown-away objects is negligible (if not positive), and we can safely neglect them in the subsequent calculations.

In conclusion, low-mass objects with mass <M_by produce ionization regions which last only for a recombination time; larger objects can survive, but their star formation ability is impaired by ISM loss/heating. The rest of the section is devoted to quantifying these two cases.

5.1 Low-mass objects: M < M_by

Due to the small mass of these protogalaxies, massive star formation and, particularly, multisupernova explosions result in a blowaway. The blowaway is always preceded (unless the galaxy is perfectly spherical) by a blowout, as explained in FT (see also MF). The blowout and blowaway conditions can be derived by comparing the blowout velocity, v_b, and the escape velocity of the galaxy, v_e. From an analysis of the Kompaneets (1957) solution for the propagation of a shock produced by an explosion in a stratified medium, it follows that the shock velocity decreases down to a minimum occurring at about 3H, where H is the characteristic scaleheight, before being reaccelerated; v_b corresponds to such a minimum.

To calculate these two velocities, we define a protogalaxy as a two-component system made of a dark matter halo and a gaseous disc. We assume a modified isothermal halo density profile ρ_h(r)=ρ_c/[1+(rr₀)²] extending out to a radius r_h={3M_h/[4π(18π²)ρ_crit]}^1/3, defined as the characteristic radius within which the mean dark matter density is 18π² times the critical density g cm⁻³; M_h is the halo mass. For such a halo the escape velocity can be written as

(18)

where Φ(r_d) is the the gravitational potential of the halo (in this calculation we neglect the contribution of baryons to the potential) calculated at the radius of the disc, r_d; p± 1.65; the assumption r_d≫r₀ has been made.

The explicit expression for v_b has been obtained by FT:

(19)

where L is the mechanical luminosity, and 〈ρ_d〉 is the disc mean gas density. L at redshift z can be written as

(20)

where ε₀=10⁵¹ erg is the energy of a SN explosion; as previously stated, we assume a Salpeter IMF, according to which one supernova is produced for each 135 M_⊙ν⁻¹ of stars formed; the free-fall time is t_ff=(4πG〈ρ_h〉)^−1/2.

Due to the dissipative nature of the gas, the baryons which should be initially distributed approximately as the dark matter will lose pressure and collapse in the gravitational field of the dark matter. If the halo is rotating, the gas will collapse in a centrifugally supported disc and, assuming that all the baryons bound to the dark matter halo are able to collapse in the disc, The radius of the disc can be estimated by imposing that the specific angular momentum of the disc, j_d, is equal to that of the halo, j_h (Mo, Mao & White 1998; Weil, Eke & Efstathiou 1998):

(21)

where v_c is the circular velocity, ℓ_d is the disc scalelength, and λ the standard halo spin parameter. As shown by numerical simulations (Barnes & Efstathiou 1987; Steinmetz & Bartelmann 1995), λ is only very weakly dependent on M_h and on the density fluctuation spectrum; its distribution is approximately log-normal and peaks around λ± 0.04. From equation (21) we then obtain For an exponential disc, as implicitly assumed in deriving equation (21) above, the optical radius (i.e. the radius encompassing 83 per cent of the total integrated light) is 3.2 ℓ_d; then, the H i radius in galaxies is typically found to be ≈2 times larger than the optical radius (Salpeter & Hoffman 1996). Thus, for the radius of the gaseous disc we take

(22)

The scaleheight, H, is roughly given by

(23)

where c_s is the effective gas sound speed, including also a possible turbulent contribution; we take c_s≈1 km s⁻¹, corresponding to the minimum temperature allowed by molecular cooling. It is useful to note that i.e. it is independent of mass and redshift. We point out that these estimates are in excellent agreement with the properties of Pop III discs found in numerical simulations of Bromm et al. (1999).

In the following we derive the necessary condition for blowaway to occur. Following blowout, the pressure inside the hot gas bubble drops suddenly due to the inward propagation of a rarefaction wave. The lateral walls of the shell, moving along the galaxy major axis, will continue to expand unperturbed until they are overcome by the rarefaction wave at r_c, corresponding to a time t_c elapsed from the blowout. After that moment, the shell enters the momentum-conserving phase, since the driving pressure has been dissipated by the blowout. The requirement for the blowaway to take place is then that the momentum of the shell (of mass M_c at r_c) is larger than the momentum necessary to accelerate the gas outside r_c at a velocity larger than the escape velocity:

(24)

With the assumptions and algebra outlined in FT, one can write the condition for blowaway to occur:

(25)

where v_b is the blowout velocity, v_e is the escape velocity, a is a parameter equal to 2/3, and εr_dH is the ratio of the major axis to the minor axis of the protogalaxy. From equations (22) and (23), we find

(26)

Substituting equations (18), (19) and (26) into equation (25), we find the critical mass below which blowaway occurs. It is useful to recall that this prescription has been carefully tested against hydrodynamical simulations presented in MF.

In Fig. 6, M_by and M_crit are plotted for the usual cosmological parameters (Ω₀=1, h± 0.5 and Ω_b=0.06) and for the three runs A, B and C, whose parameters are defined in Table 1 and described in detail in Section 7. As M_by depends on the product ƒ_b׃*, all other runs give the same results as run A. As is seen from the figure, the blowaway mass depends very weakly on this product [approximately ∝( ƒ_bƒ*)^1/4], and it is always larger than M_crit, but the number of objects suffering blowaway decreases with redshift. Note as that the vast majority of the objects which virialized at redshifts ≳15 will be blown away due to their low mass predicted by hierarchical models. We also point out that at redshifts below ≈6, when objects with masses ≥10⁹ M_⊙ start to dominate the mass function, our estimate of M_by becomes less accurate. This is due to the assumption implicitly made that all SNe drive the same superbubble. Though valid for small objects, for larger galaxies this idealization becomes increasingly poorer, and one should consider the fact that in general SNe are distributed among OB associations with different values of L spread around the galactic disc. However, our calculations are not affected by this problem, as we stop the simulations at z≈8.

Figure 6.

Minimum total mass for collapse of an object as a function of redshift according to Tegmark et al. (1997) (M_crit, solid line) and maximum mass for blowaway as a function of redshift (M_by, dotted line). M_by is shown for the three runs A, B and C (see Table 1).

Table 1.

Parameters of the calculation: fraction of virialized baryons that are able to cool and be available for star formation, ƒ_b; star formation efficiency, ƒ*; photon escape fraction, ƒ_esc; object cooling mechanism; dark matter fluctuation power-spectrum normalization, σ₈.

5.2 Larger objects: M > M_by

Objects with M>M_by will not experience blowaway and, rather than explode, they continue to produce stars quiescently. Nevertheless, a blowout with consequent mass-loss and heating of the ISM by SN shock waves will very likely take place. MF found that galaxies with masses lower than a few ×10¹⁰ M_⊙ will suffer losses in an outflow, although the gas loss efficiency is not very high. Most important, perhaps, are the evaporation and destruction of molecular clouds as a result of an enhanced star formation activity which depresses the star formation rate. These processes are not yet fully understood in their complexity, and detailed hydrodynamical models implementing the physics regulating the multiphase structure of the ISM are required to make further progress. We parametrize these processes by an heuristic approach prescribing that the mass transformed into stars is lower than the value M* used for low-mass objects by a mass-dependent factor s_feed:

(27)

where γ is a free parameter usually taken to be equal to 2. We have checked that our results are almost insensitive to different values of γ. This parametrization of the stellar feedback is analogous to the one typically used in semi-analytical models of galaxy formation (White & Frenk 1991; Kauffmann 1995; Baugh et al. 1998; Guiderdoni et al. 1998). This factor allows for a maximum 50 per cent decrease of the stellar mass of a galaxy due to gas heating/loss; this stellar feedback becomes increasingly less important for larger objects.

6 Summary of evolutionary tracks

Given the large number of processes so far discussed and included in the model, it is probably worth summarizing them briefly before presenting the results. Fig. 7 illustrates all possible evolutionary tracks and final fates of primordial objects, together with the mass-scales determined by the various physical processes and feedbacks. We recall that there are four critical mass-scales in the problem: (i) M_crit, the minimum mass for an object to be able to cool in a Hubble time; (ii) M_H, the critical mass for which hydrogen Lyα line cooling is dominant; (iii) M_sh, the characteristic mass above which the object is self-shielded, and (iv) M_by, the characteristic mass for stellar feedback, below which blowaway cannot be avoided. Starting from a virialized dark matter halo, condition (i) produces the first branching, and objects failing to satisfy it will not collapse and form only a negligible amount of stars. In the following we will refer to these objects as dark objects. Protogalaxies with masses in the range M_crit<M<M_H are then subject to the effect of radiative feedback, which could either impede the collapse of those of them with mass M<M_sh, thus contributing to the class of dark objects, or allow the collapse of the remaining ones (M>M_sh) to join those with M>M_H in the class of luminous objects. This is the class of objects that convert a considerable fraction of their baryons in stars. Stellar feedback causes the final bifurcation by inducing a blowaway of the baryons contained in luminous objects with mass M<M_by; this separates the class in two subclasses, namely ‘normal’ galaxies (though of masses comparable to present-day dwarfs) that we dub gaseous galaxies, and tiny stellar aggregates with negligible traces (if any) of gas to which consequently we will refer to as naked stellar clusters. The role of these distinct populations for the reionization of the Universe and their density evolution will be clarified by the following results.

Figure 7.

Possible evolutionary tracks of objects as determined by the processes and feedbacks included in the model. Explanations are given in Section 6.

7 Results

In this section we present the results obtained by including all the effects discussed above in the calculation of the reionization history. As the cosmological model (CDM with Ω₀=1, h± 0.5, Ω_b=0.06) is fixed by the N-body simulations, there are three free parameters left, namely ƒ_b, ƒ* and ƒ_esc, defined in Section 3 (actually, as already pointed out, the number of effective free parameters is two, ƒ_b׃* and ƒ_esc). We have performed runs with different values for these parameters, labelled as runs A-D in Table 1; our reference case is run A ( ƒ_b=0.08, ƒ*=0.15 and ƒ_esc=0.2). We have made runs with the same parameters as in run A, but in which only objects with masses M>M_H are considered (run A1) and the normalization of the fluctuation spectrum has been varied (run A2). These cases are intended to isolate the effects of negative feedback and of a different evolution of the dark matter haloes, as explained below.

The present approach, based on a combination of very high-resolution N-body numerical simulations and detailed treatment of the radiative/stellar feedbacks, allows us to study the inhomogeneous ionization process and the nature of the sources producing it at a remarkably detailed level. The main output of our model consist of simulation boxes in which ionized/dissociated regions can be clearly identified, and for which we can derive the global three-dimensional structure as a function of cosmic time.

7.1 Filling factor

As a typical example, we show in Fig. 8 the actual outcome for run A. There, the spatial distribution of reionized regions is shown at two different redshifts (z± 22.1 and 19.8). The images shown in Fig. 8 give an immediate and qualitative view of the evolution of the volume occupied by ionized atomic hydrogen. A more quantitative measure of the filling factor of the ionized gas can easily be obtained.

Figure 8.

Global spatial structure of ionized regions at (a) z± 22.1 and (b) z± 19.8.

The filling factors of the dissociated H₂ and ionized H are defined as the box volume fraction occupied by those species. The results are shown in Figs 9(a) and (b) for different runs. The intergalactic relic molecular hydrogen is found to be completely dissociated at very high redshift (z≈25), independently of the parameters of the simulation. This follows from the fact that dissociation spheres are relatively large and overlap at early times. Ionization spheres, on the other hand, are always smaller than dissociation ones, and complete reionization occurs considerably later. Except for run C, when reionization occurs by z≈15, primordial galaxies are able to reionize the IGM at a redshift z≈10. The filling factor is approximately linearly proportional to the number of sources, but has a stronger (cubic) dependence on the ionization sphere radius. Thus, although the number of (relatively small) luminous objects present increases with redshift, the filling factor is only boosted by the appearance of larger, and therefore more luminous, objects which can ionize more efficiently. The subsequent flattening at lower redshifts is obviously due to the fact that when the volume fraction occupied by the ionized gas approaches unity, most photons are preferentially used to sustain the reached ionization level rather than to create new ionized regions.

Figure 9.

(a) Dissociated molecular hydrogen (top set of lines), ionized atomic hydrogen (middle set) and doubly ionized helium (bottom line) filling factor as a function of redshift for different runs: A (solid line), B (dotted), C (dot-dashed) and D (dashed). (b) Same as (a) for runs: A (solid line), A1 (dotted) and A2 (dot-dashed).

In principle, a higher photon injection in the IGM could result both in an increase (as larger H ii regions are produced) or a decrease (as the number of sources is reduced by the effect of radiative feedback) of the filling factor. Along the sequence run B, D, A, C the number of ionizing photons injected in the IGM is progressively increased. Then Fig. 9(a) allows us to conclude that the former effect is dominating, i.e. the radiative feedback is of minor importance.

As a comparison, we have calculated the evolution of the H ionization fraction in a simplified manner, by solving the time-dependent ionization-recombination equation, dx/dtγ_p(t)(1−x) −αx². This time we compute the photoionization rate using the total, volume-averaged ionizing photon rate S_i(0)^V, as deduced from our simulations. A reasonable fit to this quantity is S_i(0)^V=4×10⁵⁷ exp[− (1+z)^0.85]. Clearly, x must also be understood as the volume-averaged ionization fraction. For run A, the condition x± 1 (full ionization) is obtained at redshift z± 14.2, therefore slightly before that concluded above. This discrepancy is due to the fact that assuming a spatially constant (or average) photoionization rate neglects the fact that the IGM inside the simulation volume is optically thick on scales typical of the Strömgren spheres and there are ‘dark’ regions which are not affected by ionizing photons until a later epoch, when overlapping has become important and a true background originates.

For the specific case of run A (we do not discuss in further detail He reionization in this paper) we show the analogous filling factor for the doubly ionized He; the lowest curve in Fig. 9(a) shows its evolution. As R_s,He≃0.1R_s (see equations 5 and 7), the helium filling factor should be approximately 0.1 per cent of the hydrogen filling factor, but it increases when the effects of the ionized hydrogen sphere overlapping is included (see Fig. 9a). The stellar spectrum contributed by the reionizing galaxies is relatively soft, and the number of photons above the He⁺⁺ ionization threshold very limited. Hence He reionization does not occur in our model, and could possibly be achieved only through quasar-like sources with harder spectra.

In Fig. 9(b) we analyse the effect of negative feedback and the normalization of the dark matter fluctuation power spectrum, σ₈. By comparing run A with run A1, in which only objects with masses M>M_H are considered, we see that, apart from very high redshifts where few such objects are present and the bulk of photon production is due to objects with M<M_H, the main contribution to the IGM reionization comes from sources with M>M_H and Pop III objects (with M<M_H) alone would not be able the reionize the IGM. This statement remains valid even in the absence of radiative feedback effects on Pop III objects, as indicated from the results of a run A in which radiative feedback has not been allowed. In this case, although a slightly higher filling factor is obtained as the contribution from previously suppressed Pop III objects is now present, the difference with run A is never higher than 10 per cent.

Finally, in run A2, the normalization of the fluctuation spectrum has been decreased to σ₈=0.5. The evolution of the filling factor appears similar to that in run A, but shifted towards lower redshift. This is due to the delayed evolution of this model with respect to the reference one, which results in a lower number of luminous objects at a given redshift.

7.2 SUVB

The derived evolution of the SUVB is presented in Fig. 10(a) at the mean frequency of the LW band, hν₀=12.45 eV, for runs C, A, D and B (from top to bottom). The intensity scales with the product ƒ_b׃*׃_esc as deduced from equations (2), (15) and (16). Typically, a flux lower than 10⁻²⁴ erg cm⁻² s⁻¹ Hz⁻¹ sr⁻¹ is found at z≳15, except for case C for which it is found to be about one order of magnitude larger. As an intensity of at least few ×10⁻²⁴ erg cm⁻² s⁻¹ Hz⁻¹ sr⁻¹ at the Lyman limit is needed for an appreciable negative feedback at these redshifts (see Fig. 5), we conclude, in agreement with CFA, that the SUVB flux is too weak to affect low-mass structure formation at high redshift. At these early epochs, negative feedback is instead induced by the direct flux from pre-existing nearby objects. The SUVB becomes intense and dominates the direct flux for z<15 (earlier for run C), and consequently it governs the formation of late forming Pop III objects. The additional cases for runs A1 and A2 are shown in Fig. 10(b). Run A1 produces the same SUVB as run A at low redshift; at early times, where there is little contribution from objects with M>M_H, the SUVB intensity is lower. The case with a lower normalization (run A2) also produces a lower SUVB, as expected.

Figure 10.

(a) SUVB intensity as a function of redshift, at the central frequency of the Lyman-Werner band hν₀=12.45 eV for runs A-D. (b) Same as (a) for runs A, A1 and A2.

7.3 SFR

We have calculated the star formation rate per comoving volume (SFR) predicted by our model and we have compared it with the values derived by the most recent studies at lower redshift (Lilly et al. 1996; Connolly et al. 1997; Madau, Pozzetti & Dickinson 1998; Steidel et al. 1999). Inside each simulation box we have derived the SFR as where V is the comoving volume of the simulation box and sfr_i are the star formation rates of the individual luminous objects:

(28)

in units of M_⊙ yr⁻¹. Here M*_,i is given by equation (1), showing that the SFR depends crucially on the parameters ƒ_b and ƒ*; t_e is defined as the free-fall time for objects with mass M>M_by, while for low-mass objects (M<M_by), witnessing only a single burst of star formation, we assume it to be equal to the Hubble time:

(29)

where As we have already pointed out, the simulations, though resolving low-mass objects, might miss haloes larger than M_max,s; this value increases with decreasing redshift. To take into account this limitation and correct for it, we include also objects with masses M>M_max,s and model their contribution as follows. At each redshift, this contribution is given by

(30)

where M* is given in equation (1), while N(M) and M_max are defined in equation (16). We implicitly assume that these objects do not suffer from radiative feedback (this assumption is justified by the fact that these are usually objects with masses M>M_H), and that these haloes host luminous objects. We then add this contribution to the SFR derived from the simulations. The results are shown in Fig. 11(a) for runs A-D. Obviously, the higher the ƒ_b׃* product, the higher is the SFR obtained. Actually, although runs A and D use the same values for the above parameters, in run D more luminous objects are formed [as the SUVB is lower (see Section 7.2), more objects escape negative feedback] and this results in a slightly higher SFR. We compare the curves obtained with the most recent measurements of the SFR, namely Lilly et al. (1996) [circles], Madau et al. (1998) [triangles], Connolly et al. (1997) [squares] and Steidel et al. (1999) [crosses]. The recent lower limit to the SFR at z≈3 set by SCUBA (Hughes et al. 1998) coincides with the point corresponding to Steidel et al. (1999). The points from Madau et al. (1998) are derived from a study of the Lyman break galaxies in the Hubble Deep Field (HDF) and show a decline of the SFR at high redshift. The measurements by Steidel et al. do not show the same decline and favour a more or less constant SFR at high redshift. These determinations are based on a study of star-forming galaxies at z>3.8 in a field much wider than the HDF, though shallower. All points are corrected for dust extinction as in Steidel et al.'s paper. From Fig. 11(a) we see that runs B and C can be taken respectively as lower and upper limits to the SFR obtained from our model, thus constraining the product ƒ_b׃*. Runs A and D produce a trend which apparently matches the observations well, suggesting a likely value of the product ƒ_b׃* around ≈0.01. In Fig. 11(b) we study the effect of negative feedback and σ₈. From a comparison between runs A and A1, where only objects with M>M_H are considered, we see that, while at high redshift the main contribution to the SFR comes from low-mass objects (M<M_H), at lower redshift, their contribution becomes negligible. Finally, run A2, with σ₈=0.5, produces a lower SFR as less luminous objects are present at a given redshift, again the same point made in Section 7.1.

Figure 11.

(a) Comoving star formation rate as a function of redshift, for runs C, D, A and B, from top to bottom, respectively. Data points are from Lilly et al. (1996) [circles]; Connolly et al. (1997) [squares]; Madau, Pozzetti & Dickinson (1998) [triangles] and Steidel et al. (1999) [crosses]. (b) Same as (a) for runs A, A1 and A2.

Additionally, we have derived the evolution of the stellar density parameter, Ω*(z)=ρ*_,tot(z)/ρ_crit, produced by run A; ρ*_,tot(z) is the total mass density in stars at redshift z inside the comoving volume of the box V± 125 Mpc³, and is given by

(31)

where t_on is the Hubble time corresponding to z_on=29.6. We find that by the time the IGM is completely reionized, z≈10, only 2 per cent of the stars observed at z± 0 (Ω*(0) ≈4.9×10⁻³ (Fukugita, Hogan & Peebles 1998) has been produced. This apparently small number of stars is sufficient to produce the necessary ionizing flux. The number of ionizing photons produced by z≈10 can be derived as follows. At z≈10, M*_,tot≈3.8×10⁸ M_⊙, corresponding to an ionizing photon rate S_i(0) ≈2.62×10⁵⁴ s⁻¹ for the standard case A (see equation 4). The time at which the lowest mass OB stars expire is ≈40 Myr (Oey & Clarke 1997); thus the total number of ionizing photons escaping from luminous objects is N_i≈3.3×10⁶⁹. As the number of baryons in the box is N_bVΩ_bρ_crit/(μm_p) ≈4.7×10⁶⁸, where μ≈1.27 is the mean molecular weight, there are N_i/N_b≈7 ionizing photons per baryon available. At z≈10 the Hubble time is t_H≈3.59×10⁸ yr, while the mean hydrogen recombination time over the time interval z≈30-10 of interest here is 〈t_rec〉≈1.33×10⁸ yr. Thus, typically an atom has recombined about t_H/〈t_rec〉≈3 times. However, the ionizing photon budget provided by the formed stars is sufficient to balance this recombination and keep the gas ionized.

A rough estimate of the mean metallicity produced by these stars can be calculated. Given the adopted IMF and the total mass in stars at z≈10, M*_,tot, we can derive the number of SNe in the simulation box and, assuming that each SN produces ≈1 M_⊙ of heavy elements, we find that the total mass of heavy elements at z≈10 is M_met≈2.81×10⁶ M_⊙. As the total mass in baryons in the considered cosmic volume is M_b,tot≈5.06×10¹¹ M_⊙, we find that the mean metallicity at z≈10 is

(32)

7.4 Galaxy evolutionary tracks

In this section we discuss the final fates of primordial objects, and we show their relative numbers in Fig. 12. The straight lines represent, from top to bottom, the number of dark matter haloes, dark objects, naked stellar clusters and gaseous galaxies, respectively. The dotted curve represents the number of luminous objects with large enough mass (M>M_H) to make the H line cooling efficient and become insensitive to the negative feedback. We remind that the naked stellar clusters are the luminous objects with M<M_by, while the gaseous galaxies are the ones with M>M_by; thus the number of luminous objects present at a certain redshift is given by the sum of naked stellar clusters and gaseous galaxies. We first notice that the majority of the luminous objects that are able to form at high redshift will experience blowaway, becoming naked stellar clusters, while only a minor fraction, and only at z≲15, when larger objects start to form, will survive and become gaseous galaxies. An always increasing number of luminous objects is forming with decreasing redshift, until z≈15, where a flattening is seen. This is due to the fact that the dark matter halo mass function is still dominated by low-mass objects, but a large fraction of them cannot form due to the following combined effects: (i) towards lower redshift the critical mass for the collapse (M_crit) increases, and fewer objects satisfy the condition M>M_crit; (ii) the radiative feedback due to either the direct dissociating flux or the SUVB (see Section 4) increases at low redshift as the SUVB intensity reaches values significant for the negative feedback effect. When the number of luminous objects becomes dominated by objects with M>M_H, by z≈10 the population of luminous objects grows again, basically because their formation is now unaffected by negative feedback. A steadily increasing number of objects is prevented from forming stars and remains dark; this population is about ≈99 per cent of the total population of dark matter haloes at z≈8. This is also due to the combined effect of points (i) and (ii) mentioned above. This population of haloes which have failed to produce stars could be identified with the low-mass tail distribution of the dark galaxies that reveal their presence through gravitational lensing of quasars (Hawkins 1997; Jimenez et al. 1997). It has been argued that this population of dark galaxies outnumbers normal galaxies by a substantial amount, and Fig. 12 supports this view. At the same time, CDM models predict that a large number of satellites (a factor of 100 more than observed (Klypin et al. 1999; Moore et al. 1999b) should be present around normal galaxies. Many of them would form at redshifts higher than 5 and would survive merging and tidal stripping inside larger haloes to the present time. Their existence is then linked to that of low-mass primordial objects, and the natural question arises if these objects can be reconciled with the internal properties of haloes of present-day galaxies. Similarly to what found here, Abel & Mo (1998) have suggested that a large fraction of high-redshift Lyman limit absorption systems (LLS) seen in quasar spectra originate from gas trapped in minihaloes formed prior to reionization, which are able to survive the onset of the UV background at a later time and in which star formation is quenched by their inability to cool.

Figure 12.

(a) Number evolution of different objects in the simulation box for run A (see discussion in Section 7.4) (b) Same as (a) for run C.

A question that naturally arises is what is the fate of the naked stellar objects. The actual number of naked stellar objects should not be much higher than that at z≈8, for two reasons: (i) with decreasing redshift, a lower number of luminous objects is subject to blowaway, as increasingly larger objects form; (ii) as already pointed out, at redshifts below ≈6, our estimate of M_by becomes less accurate, actually larger than the real value (see FT); thus the number of luminous objects undergoing a blowaway is decreased also by this effect. With this caveat, approximately 2 naked stellar objects should end up in the Galactic halo, where, at the present time, only stars with masses in the range 0.1–1 M_⊙ are still shining. As the maximum total mass of a naked stellar object at z≈8 is ≈10⁹ M_⊙, given a Salpeter IMF, we find that the upper limit for the number of the above relic stars is ≈4×10⁶. As in the Galaxy, given the same IMF, ≈2×10¹¹ stars with masses in the range 0.1–1 M_⊙ are present, one out of ≈5×10⁴ stars comes from naked stellar objects. The upper limit on the metallicity of such stars is the IGM one at z≈10, i.e. Z≈3×10⁻⁴ Z_⊙, as derived in equation (32).

We have compared the results for runs A and C. The only significant difference is that in case C, where a higher SUVB is found, the redshift at which the M>M_H population peaks is higher, as from z≈15 there are very few objects that can avoid negative feedback, and the only surviving haloes are the ones for which H line cooling is efficient.

8 Comparison with previous work

In the last decade there have been a number of different works dedicated to reionization. The ‘first generation’ of studies used very simplified descriptions of the problem in terms of the treatment of the physical processes, cosmological evolution of the ionizing sources, or both. These pioneering works (Fukugita & Kawasaki 1994; Tegmark, Silk & Blanchard 1994; Giroux & Shapiro 1996) found results often in disagreement with each other, probably due to the different effects/cosmological models included or neglected. More recent works have tackled again the problem, improving on at least two crucial ingredients: a detailed cosmological hydrodynamic evolution of the IGM/galaxy formation, which gives a more detailed distribution of ionizing sources and of the gas, and the inclusion of newly discovered feedback effects, along the lines presented here. GO performed detailed simulations for a CDM+Λ cosmological model, with an approximate treatment of the radiative transfer and allowing for IGM clumping. The N-body simulations used here have a larger number of particles (256³ instead of their 64/128³), although the size of the box is practically the same. Their spatial resolution (≈1 kpc) can only marginally resolve the H ii spheres which, for Pop III objects, are often smaller than that size. Also, we treat in more detail the effects of radiative and stellar feedback, even if we cannot compete with the wealth of information on the IGM density distribution that can be extracted from their simulations. They conclude that the reionization redshift is ≈7, i.e. about 3 redshift units below the one found here. The discrepancy might depend both on the different cosmological model adopted and/or on the above differences. A low value for the reionization redshift is found also in the semi-analytical model presented by VS: for both an open and a critical density universe, they use a mass function alternative to the usually adopted PS formalism and, given a prescription for galaxy formation and photon production, they find that the reionization is complete by redshift 6.8 for the open universe and 5.6 for the critical-density one. As in their calculation they include the treatment of the gas clumpiness, this could again explain the discrepancy with the value found in the present work. Additionally, deriving the radius of the H ii region surrounding a source, they neglect the influence of the expansion of the universe, underestimating the real radius and obtaining a late reionization epoch. Haiman & Loeb (1997) used a semi-analytical model based on PS formalism which adopted a strong version of the radiative feedback (but not considering the stellar one), assuming the suppression of star formation inside objects with virial temperatures below 10⁴ K. They also neglect the effect of gas clumpiness. This approach has the advantage of being very flexible when exploring the dependence of reionization on various parameters, but clearly lacks the crucial information about the spatial distribution of ionizing sources. They find that for a large range of CDM models, reionization occurs at a redshift ≳10, consistent with the value we find in the present work. The basic agreement with their result derives from the fact that we also find that reionization is mostly driven by objects collapsed through H line cooling. Stellar feedback instead introduces a different prediction for the fate of ionizing objects as discussed above.

9 Discussion

In this paper we have studied the reionization of the Universe due to an inhomogeneous distribution of sources, including a number of relevant physical processes and feedback effects whose importance for this type of studies has only recently been recognized. Our approach provides a reliable picture of the actual process of IGM reionization.

May be the major feature not considered here is constituted by the possible IGM density inhomogeneities, as throughout the paper we have assumed an homogeneous gas distribution. The IGM clumpiness can have, in principle, different effects on the reionization of the Universe. Both GO and VS have studied the process of reionization, taking into account the effect of the inhomogeneous gas distribution through a clumping factor that, included in the calculation of the recombination process, increases the number of ionizing photons necessary to reionize the IGM by the amount C, thus delaying the IGM complete reionization with respect to the homogeneous case. The gas clumpiness influences also the shape of the ionized regions: indeed, the propagation of the ionizing front in an inhomogeneous medium, does not result in a spherically symmetric H ii region, with ionization proceeding in an inside-out fashion, but rather the reionization occurs outside-in, starting in voids and gradually penetrating in high-density regions, as pointed out by Miralda-Escudé et al. (2000). If this is the case, recombinations from overdense regions which are presumably very compact, would contribute to the average recombination rate only at later times, when most of the volume has already been ionized. Note that this effect is opposite to the first one, indicating that at a given redshift in the homogeneous case one might be underestimating the filling factor of the ionized gas. An even more refined approach to the problem should include the combined treatment of the complete radiative transfer equation in an inhomogeneous medium. The very first steps in this directions have been presented by Norman et al. (1998), who are trying to incorporate radiative transfer into 3D hydrodynamic cosmological simulations. As a preliminary study of the effect of the gas clumpiness on our results, we have made additional runs in which we include a clumping factor in the calculation of the recombination process. For C we have adopted the curves given in GO and VS. In Fig. 13 the H ii filling factor for run A is shown, together with the curves derived in the presence of clumping. The GO clumping factor (dashed curve) produces little difference, as it is close to unity until z≈10. On the other hand, the VS one (dotted line), reaching much higher values, results in a substantially lower filling factor and a delay in the IGM reionization, which is not complete by z≈8. The possibility exists that a non-negligible fraction of baryons might be in collapsed structures, while throughout the paper we assume that the baryons reside predominantly in the IGM. However, QSO absorption-line experiments show that at z≈3 a significant fraction of all baryons is in Lyα forest systems (perhaps as much as 50–100 per cent; Rauch et al. 1997; Weinberg et al. 1997); in particular, at z± 2.5 the baryon density parameter in the Lyα forest is greater than 0.05 (Pettini 1999), which implies that for our choice of Ω_b, more than 80 per cent of baryons are probably not in collapsed objects; this fraction very likely increases with redshift. Thus we consider our assumption at least as a reasonable one.

Figure 13.

Ionized atomic hydrogen filling factor as a function of redshift for different runs: A (solid line), A with the clumping factor curve of GO (dashed), and VS (dotted).

Another important output of the model is the star formation history. A direct comparison with most of the previous calculations of the SFR, both in numerical and semi-analytical approaches, can be only approximate, as those concentrated mainly on lower redshift ranges than the ones studied here (see, e.g., Guiderdoni et al. 1998, Baugh et al. 1999, Pei, Fall & Hauser 1999 and Nagamine, Cen & Ostriker 1999). VS, on the other hand, calculate the SFR up to z≈20, obtaining values slightly lower than both ours and those derived observationally.

Finally, the result that a high fraction of dark objects is present is intriguing, as the question of whether dark galaxies can exist is a long-standing one (Dekel & Silk 1986). Indeed, the fact that the luminosity function has a flatter faint-end slope than the associated halo mass function can be explained with an increasing fraction of dark galaxies towards lower masses. Some authors (Babul & Rees 1992; Metcalfe et al. 1997) claim that the initial burst of the star formation in low-mass objects may provide a population of disappearing dwarfs proposed for interpreting the faint galaxy counts. On the other hand, Jimenez et al. (1997) propose that a dark galaxy forms as a low-density disc in a dark halo of high spin parameter: such a galaxy may have a surface density too low to produce a significant star formation rate. Our model offers an alternative explanation to the nature of such objects.

A comparison of the predictions of our model is currently possible only in a indirect way, as done throughout the paper for several quantities. However, in the near future, thanks to the availability of powerful space and ground-based instrumentation, it will become possible to directly probe most of the activity occurring prior to complete reionization. For example, NGST is expected to detect a large number of high-z SNe (Marri & Ferrara 1998) that could be used to trace early star formation and stellar feedback in Pop III objects, and possibly the IGM, via absorption experiments. NGST will be also able to determine with great accuracy the reionization epoch from the spectra of high-redshift sources (Haiman & Loeb 1999a). The ionized regions can be mapped directly in free-free emission (Oh 1999) by the Square Kilometer Array (SKA). Reionization should leave a measurable imprint on the CMB, as secondary anisotropies are produced via the kinetic Sunyaev-Zel'dovich effect by inhomogeneous reionization (Gruzinov & Hu 1998; Knox, Scoccimarro & Dodelson 1998); this experiment is well into the capabilities of MAP and Planck missions. Submillimetre telescopes will be used to search for dust produced by primordial objects, and possibly molecules; the dust can be an important reservoir of the metals synthesized by the first SNe, and it might have important consequences for the exploration of the early Universe. Finally, radio surveys might reveal signatures of reionization and give insight into the thermal evolution of neutral hydrogen through its redshifted H i 21-cm line emission. These features are within reach of new generation radio telescope like SKA (Madau, Meiksin & Rees 1997; Tozzi et al. 2000). Clearly, exciting times are ahead of us.

10 Summary

We have studied inhomogeneous reionization in a critical density CDM model, including a detailed treatment of radiative and stellar feedback processes on galaxy formation. The main results discussed in this paper can be summarized as follows.

• (1)

  Galaxies are able to reionize the neutral atomic hydrogen by a redshift z≈10, while molecular hydrogen is completely dissociated at very high redshift (z≈25).

• (2)

  IGM reionization is basically driven by objects collapsed through H line cooling (M>M_H), while low-mass objects (M<M_H) play only a minor role and, even in the absence of a radiative negative feedback, they would not be able the reionize the IGM.

• (3)

  The soft-UV background intensity is too low to produce sensible negative feedback effects on low-mass galaxy formation at redshift z≳15, where the radiative feedback is dominated by the direct flux from pre-existing objects. At lower redshifts the SUVB is instead dominant and reaches interesting values to influence the subsequent galaxy formation process.

• (4)

  The evolution of the star formation rate obtained shows a trend consistent with the most recent measurements and the match constrains the baryon conversion efficiency into stars in a narrow range around ≈0.01. Only about 2 per cent of the stars observed at z± 0 are required to reionize the Universe. This corresponds to an average IGM metallicity at redshift z≈10 equal to 〈Z〉≈3×10⁻⁴ Z_⊙.

• (5)

  A consistent fraction of haloes is prevented from forming stars by either the condition M<M_crit or the effect of radiative feedback; this population of dark objects reaches ≈99 per cent of the dark matter halo population at z≈8.

Acknowledgments

We thank Z. Haiman, L. Pozzetti and M. Tegmark for providing useful data, V. Bromm, D. Galli, F. Haardt, F. Palla and M. Rees for stimulating discussions, and the referee, T. Abel, for useful comments. BC especially thanks A. Riccardi for his help with software. The work presented in this paper was carried out using data made available by the Virgo Supercomputing Consortium ( frazerp/virgo/virgo.html) and computers based at the Computing Centre of the Max-Planck Society in Garching and Edinburgh Parallel Computing Centre.

1

http://www-hpcc.astro.washington.edu/tools/FOF/