Lecture notes for the SAAS-Fee Winter School, April 2006 (to be published by Springer Verlag)
astro-ph/0603360

For a PDF version of the article, click here.

Abstract. The first dwarf galaxies, which constitute the building blocks of the collapsed objects we find today in the Universe, had formed hundreds of millions of years after the big bang. This pedagogical review describes the early growth of their small-amplitude seed fluctuations from the epoch of inflation through dark matter decoupling and matter-radiation equality, to the final collapse and fragmentation of the dark matter on all mass scales above ~ 10^-4 M. The condensation of baryons into halos in the mass range of ~ 10⁵ - 10¹⁰ M led to the formation of the first stars and the re-ionization of the cold hydrogen gas, left over from the big bang. The production of heavy elements by the first stars started the metal enrichment process that eventually led to the formation of rocky planets and life.

A wide variety of instruments currently under design [including large-aperture infrared telescopes on the ground or in space (JWST), and low-frequency arrays for the detection of redshifted 21cm radiation], will establish better understanding of the first sources of light during an epoch in cosmic history that was largely unexplored so far. Numerical simulations of reionization are computationally challenging, as they require radiative transfer across large cosmological volumes as well as sufficently high resolution to identify the sources of the ionizing radiation. The technological challenges for observations and the computational challenges for numerical simulations, will motivate intense work in this field over the coming decade.

Disclaimer: This review was written as an introductory text for a series of lectures at the SAAS-FEE 2006 winter school, and so it includes a limited sample of references on each subject. It does not intend to provide a comprehensive list of all up-to-date references on the topics under discussion, but rather to raise the interest of beginning graduate students in the related literature.

Table of Contents

OPENING REMARKS

EXCAVATING THE UNIVERSE FOR CLUES ABOUT ITS HISTORY

BAKGROUND COSMOLOGICAL MODEL

The Expanding Universe

Composition of the Universe

Linear Gravitational Growth

The Smallest-Scale Power Spectrum of Cold Dark Matter

Structure of the Baryons

Formation of Nonlinear Objects

Spherical Collapse

Halo Properties

NONLINEAR GROWTH

The Abundance of Dark Matter Halos

The Excursion-Set (Extended Press-Schechter) Formalism

Response of Baryons to Nonlinear Dark Matter Potentials

FRAGMENTATION OF THE FIRST GASEOUS OBJECTS TO STARS

Star Formation

The Mass Function of Stars

Gamma-ray Bursts: Probing the First Stars One Star at a Time

Emission Spectrum of Metal-Free Stars

Emission of Recombination Lines from the First Galaxies

SUPERMASSIVE BLACK HOLES

The Principle of Self-Regulation

Feedback on Large Intergalactic Scales

What seeded the growth of the supermassive black holes?

RADIATIVE FEEDBACK FROM THE FIRST SOURCES OF LIGHT

Escape of Ionizing Radiation from Galaxies

Propagation of Ionization Fronts in the IGM

Reionization of Hydrogen

Photo-evaporation of Gaseous Halos After Reionization

Suppression of the Formation of Low Mass Galaxies

FEEDBACK FROM GALACTIC OUTFLOWS

Propagation of Supernova Outflows in the IGM

Effect of Outflows on Dwarf Galaxies and on the IGM

THE FRONTIER OF 21CM COSMOLOGY

Mapping Hydrogen Before Reionization

The Characteristic Observed Size of Ionized Bubbles at the End of Reionization

Separating the "Physics" from the "Astrophysics" of the Reionization Epoch with 21cm Fluctuations

MAJOR CHALLENGE FOR FUTURE THEORETICAL RESEARCH: radiative transfer during reionization requires a large dynamic range, challenging the capabilities of existing simulation codes

REFERENCES

INDEX

1. OPENING REMARKS

When I open the daily newspaper as part of my morning routine, I often see lengthy descriptions of conflicts between people on borders, properties, or liberties. Today's news is often forgotten a few days later. But when one opens ancient texts that have appealed to a broad audience over a longer period of time, such as the Bible, what does one often find in the opening chapter?... a discussion of how the constituents of the Universe (including light, stars and life) were created. Although humans are often occupied with mundane problems, they are curious about the big picture. As citizens of the Universe, we cannot help but wonder how the first sources of light formed, how life came to existence, and whether we are alone as intelligent beings in this vast space. As astronomers in the twenty first century, we are uniquely positioned to answer these big questions with scientific instruments and a quantitative methodology. In this pedagogical review, intended for students preparing to specialize in cosmology, I will describe current ideas about one of these topics: the appearance of the first sources of light and their influence on the surrounding Universe. This topic is one of the most active frontiers in present-day cosmology. As such it is an excellent area for a PhD thesis of a graduate student interested in cosmology. I will therefore highlight the unsolved questions in this field as much as the bits we understand.

2. EXCAVATING THE UNIVERSE FOR CLUES ABOUT ITS HISTORY

When we look at our image reflected off a mirror at a distance of 1 meter, we see the way we looked 6 nano-seconds ago, the light travel time to the mirror and back. If the mirror is spaced 10¹⁹ cm = 3pc away, we will see the way we looked twenty one years ago. Light propagates at a finite speed, and so by observing distant regions, we are able to see how the Universe looked like in the past, a light travel time ago. The statistical homogeneity of the Universe on large scales guarantees that what we see far away is a fair statistical representation of the conditions that were present in in our region of the Universe a long time ago.

Figure 1. Cosmology is like archeology. The deeper one looks, the older is the layer that one is revealing, owing to the finite propagation speed of light.

We can in principle image the Universe only if it is transparent. Earlier than 0.4 million years after the big bang, the cosmic plasma was ionized and the Universe was opaque to Thomson scattering by the dense gas of free electrons that filled it. Thus, telescopes cannot be used to image the infant Universe at earlier times (or redshifts > 10³). The earliest possible image of the Universe was recorded by COBE and WMAP (see Fig. 2).

Figure 2. Images of the Universe shortly after it became transparent, taken by the COBE and WMAP satellites (see http://map.gsfc.nasa.gov/ for details). The slight density inhomogeneties in the otherwise uniform Universe, imprinted hot and cold brightness map of the cosmic microwave background. The existence of these anisotropies was predicted three decades before the technology for taking this image became available in a number of theoretical papers, including [355, 308, 297, 338, 282].

Figure 3. The optical depth of the Universe to electron scattering (upper panel) and the ionization fraction (lower panel) as a function of redshift before reionization. Observatories of electromagnetic radiation cannot image the opaque Universe beyond a redshift of z ~ 1100.

3.1. The Expanding Universe

The modern physical description of the Universe as a whole can be traced back to Einstein, who argued theoretically for the so-called "cosmological principle": that the distribution of matter and energy must be homogeneous and isotropic on the largest scales. Today isotropy is well established (see the review by Wu, Lahav, & Rees 1999 [389]) for the distribution of faint radio sources, optically-selected galaxies, the X-ray background, and most importantly the cosmic microwave background (hereafter, CMB; see, e.g., Bennett et al. 1996 [36]). The constraints on homogeneity are less strict, but a cosmological model in which the Universe is isotropic but significantly inhomogeneous in spherical shells around our special location, is also excluded [155].

In General Relativity, the metric for a space which is spatially homogeneous and isotropic is the Friedman-Robertson-Walker metric, which can be written in the form

(1)

where a(t) is the cosmic scale factor which describes expansion in time, and (R, , ) are spherical comoving coordinates. The constant k determines the geometry of the metric; it is positive in a closed Universe, zero in a flat Universe, and negative in an open Universe. Observers at rest remain at rest, at fixed (R, , ), with their physical separation increasing with time in proportion to a(t). A given observer sees a nearby observer at physical distance D receding at the Hubble velocity H(t)D, where the Hubble constant at time t is H(t) = d a(t) / dt. Light emitted by a source at time t is observed at t = 0 with a redshift z = 1/a(t) - 1, where we set a(t = 0) 1 for convenience (but note that old textbooks may use a different convention).

The Einstein field equations of General Relativity yield the Friedmann equation (e.g., Weinberg 1972 [376]; Kolb & Turner 1990 [205])

(2)

which relates the expansion of the Universe to its matter-energy content. For each component of the energy density , with an equation of state p = p(), the density varies with a(t) according to the equation of energy conservation

(3)

With the critical density

(4)

defined as the density needed for k = 0, we define the ratio of the total density to the critical density as

(5)

With _m, , and _r denoting the present contributions to from matter (including cold dark matter as well as a contribution _b from baryons), vacuum density (cosmological constant), and radiation, respectively, the Friedmann equation becomes

(6)

where we define H₀ and ₀ = _m + + _r to be the present values of H and , respectively, and we let

(7)

In the particularly simple Einstein-de Sitter model (_m = 1, = _r = _k = 0), the scale factor varies as a(t) t^2/3. Even models with non-zero or _k approach the Einstein-de Sitter behavior at high redshift, i.e. when (1 + z) >> |_m^-1 - 1| (as long as _r can be neglected). In this high-z regime the age of the Universe is,

(8)

The Friedmann equation implies that models with _k = 0 converge to the Einstein-de Sitter limit faster than do open models.

In the standard hot Big Bang model, the Universe is initially hot and the energy density is dominated by radiation. The transition to matter domination occurs at z ~ 3500, but the Universe remains hot enough that the gas is ionized, and electron-photon scattering effectively couples the matter and radiation. At z ~ 1100 the temperature drops below ~ 3000K and protons and electrons recombine to form neutral hydrogen. The photons then decouple and travel freely until the present, when they are observed as the CMB [348].

3.2. Composition of the Universe

According to the standard cosmological model, the Universe started at the big bang about 14 billion years ago. During an early epoch of accelerated superluminal expansion, called inflation, a region of microscopic size was stretched to a scale much bigger than the visible Universe and our local geometry became flat. At the same time, primordial density fluctuations were generated out of quantum mechanical fluctuations of the vacuum. These inhomogeneities seeded the formation of present-day structure through the process of gravitational instability. The mass density of ordinary (baryonic) matter makes up only a fifth of the matter that led to the emergence of structure and the rest is the form of an unknown dark matter component. Recently, the Universe entered a new phase of accelerated expansion due to the dominance of some dark vacuum energy density over the ever rarefying matter density.

The basic question that cosmology attempts to answer is:

What are the ingredients (composition and initial conditions) of the Universe and what processes generated the observed structures in it?

In detail, we would like to know:

(a) Did inflation occur and when? If so, what drove it and how did it end?

(b) What is the nature of of the dark energy and how does it change over time and space?

(c) What is the nature of the dark matter and how did it regulate the evolution of structure in the Universe?

Before hydrogen recombined, the Universe was opaque to electromagnetic radiation, precluding any possibility for direct imaging of its evolution. The only way to probe inflation is through the fossil record that it left behind in the form of density perturbations and gravitational waves. Following inflation, the Universe went through several other milestones which left a detectable record. These include: baryogenesis (which resulted in the observed asymmetry between matter and anti-matter), the electroweak phase transition (during which the symmetry between electromagnetic and weak interactions was broken), the QCD phase transition (during which protons and neutrons were assembled out of quarks and gluons), the dark matter freeze-out epoch (during which the dark matter decoupled from the cosmic plasma), neutrino decoupling, electron-positron annihilation, and light-element nucleosynthesis (during which helium, deuterium and lithium were synthesized). The signatures that these processes left in the Universe can be used to constrain its parameters and answer the above questions.

Half a million years after the big bang, hydrogen recombined and the Universe became transparent. The ultimate goal of observational cosmology is to image the entire history of the Universe since then. Currently, we have a snapshot of the Universe at recombination from the CMB, and detailed images of its evolution starting from an age of a billion years until the present time. The evolution between a million and a billion years has not been imaged as of yet.

Within the next decade, NASA plans to launch an infrared space telescope (JWST) that will image the very first sources of light (stars and black holes) in the Universe, which are predicted theoretically to have formed in the first hundreds of millions of years. In parallel, there are several initiatives to construct large-aperture infrared telescopes on the ground with the same goal in mind ^{1, 2, 3}. The neutral hydrogen, relic from cosmological recombination, can be mapped in three-dimensions through its 21cm line even before the first galaxies formed [226]. Several groups are currently constructing low-frequency radio arrays in an attempt to map the initial inhomogeneities as well as the process by which the hydrogen was re-ionized by the first galaxies.

Figure 4. A sketch of the current design for the James Webb Space Telescope, the successor to the Hubble Space Telescope to be launched in 2011 (http://www.jwst.nasa.gov/). The current design includes a primary mirror made of beryllium which is 6.5 meter in diameter as well as an instrument sensitivity that spans the full range of infrared wavelengths of 0.6 - 28 µm and will allow detection of the first galaxies in the infant Universe. The telescope will orbit 1.5 million km from Earth at the Lagrange L2 point.

The next generation of ground-based telescopes will have a diameter of twenty to thirty meter. Together with JWST (that will not be affected by the atmospheric backgound) they will be able to image the first galaxies. Given that these galaxies also created the ionized bubbles around them, the same galaxy locations should correlate with bubbles in the neutral hydrogen (created by their UV emission). Within a decade it would be possible to explore the environmental influence of individual galaxies by using the two sets of instruments in concert [390].

Figure 5. Artist conception of the design for one of the future giant telescopes that could probe the first generation of galaxies from the ground. The Giant Magellan Telescope (GMT) will contain seven mirrors (each 8.4 meter in diameter) and will have the resolving power equivalent to a 24.5 meter (80 foot) primary mirror. For more details see http://www.gmto.org/

3.3. Linear Gravitational Growth

Observations of the CMB (e.g., Bennett et al. 1996 [36]) show that the Universe at recombination was extremely uniform, but with spatial fluctuations in the energy density and gravitational potential of roughly one part in 10⁵. Such small fluctuations, generated in the early Universe, grow over time due to gravitational instability, and eventually lead to the formation of galaxies and the large-scale structure observed in the present Universe.

As before, we distinguish between fixed and comoving coordinates. Using vector notation, the fixed coordinate r corresponds to a comoving position x = r / a. In a homogeneous Universe with density , we describe the cosmological expansion in terms of an ideal pressureless fluid of particles each of which is at fixed x, expanding with the Hubble flow v = H(t) r where v = dr / dt. Onto this uniform expansion we impose small perturbations, given by a relative density perturbation

(9)

where the mean fluid density is , with a corresponding peculiar velocity u v - H r. Then the fluid is described by the continuity and Euler equations in comoving coordinates [283, 284]:

(10) (11)

The potential is given by the Poisson equation, in terms of the density perturbation:

(12)

This fluid description is valid for describing the evolution of collisionless cold dark matter particles until different particle streams cross. This "shell-crossing" typically occurs only after perturbations have grown to become non-linear, and at that point the individual particle trajectories must in general be followed. Similarly, baryons can be described as a pressureless fluid as long as their temperature is negligibly small, but non-linear collapse leads to the formation of shocks in the gas.

For small perturbations << 1, the fluid equations can be linearized and combined to yield

(13)

This linear equation has in general two independent solutions, only one of which grows with time. Starting with random initial conditions, this "growing mode" comes to dominate the density evolution. Thus, until it becomes non-linear, the density perturbation maintains its shape in comoving coordinates and grows in proportion to a growth factor D(t). The growth factor in the matter-dominated era is given by [283]

(14)

where we neglect _r when considering halos forming in the matter-dominated regime at z << 10⁴. In the Einstein-de Sitter model (or, at high redshift, in other models as well) the growth factor is simply proportional to a(t).

The spatial form of the initial density fluctuations can be described in Fourier space, in terms of Fourier components

(15)

Here we use the comoving wavevector k, whose magnitude k is the comoving wavenumber which is equal to 2 divided by the wavelength. The Fourier description is particularly simple for fluctuations generated by inflation (e.g., Kolb & Turner 1990 [205]). Inflation generates perturbations given by a Gaussian random field, in which different k-modes are statistically independent, each with a random phase. The statistical properties of the fluctuations are determined by the variance of the different k-modes, and the variance is described in terms of the power spectrum P(k) as follows:

(16)

where ⁽³⁾ is the three-dimensional Dirac delta function. The gravitational potential fluctuations are sourced by the density fluctuations through Poisson's equation.

In standard models, inflation produces a primordial power-law spectrum P(k) kⁿ with n ~ 1. Perturbation growth in the radiation-dominated and then matter-dominated Universe results in a modified final power spectrum, characterized by a turnover at a scale of order the horizon cH^-1 at matter-radiation equality, and a small-scale asymptotic shape of P(k) k^n-4. The overall amplitude of the power spectrum is not specified by current models of inflation, and it is usually set by comparing to the observed CMB temperature fluctuations or to local measures of large-scale structure.

Since density fluctuations may exist on all scales, in order to determine the formation of objects of a given size or mass it is useful to consider the statistical distribution of the smoothed density field. Using a window function W(r) normalized so that d³ r W(r) = 1, the smoothed density perturbation field, d³ r (x) W(r), itself follows a Gaussian distribution with zero mean. For the particular choice of a spherical top-hat, in which W = 1 in a sphere of radius R and is zero outside, the smoothed perturbation field measures the fluctuations in the mass in spheres of radius R. The normalization of the present power spectrum is often specified by the value of ₈ (R = 8 h^-1 Mpc). For the top-hat, the smoothed perturbation field is denoted _r or _m, where the mass M is related to the comoving radius R by M = 4 _m R³/3, in terms of the current mean density of matter _m. The variance < _m >² is

(17)

where j₁(x) = (sinx - x cosx) / x². The function (M) plays a crucial role in estimates of the abundance of collapsed objects, as we describe later.

Species that decouple from the cosmic plasma (like the dark matter or the baryons) would show fossil evidence for acoustic oscillations in their power spectrum of inhomogeneities due to sound waves in the radiation fluid to which they were coupled at early times. This phenomenon can be understood as follows. Imagine a localized point-like perturbation from inflation at t = 0. The small perturbation in density or pressure will send out a sound wave that will reach the sound horizon c_s t at any later time t. The perturbation will therefore correlate with its surroundings up to the sound horizon and all k-modes with wavelengths equal to this scale or its harmonics will be correlated. The scales of the perturbations that grow to become the first collapsed objects at z < 100 cross the horizon in the radiation dominated era after the dark matter decouples from the cosmic plasma. Next we consider the imprint of this decoupling on the smallest-scale structure of the dark matter.

3.4. The Smallest-Scale Power Spectrum of Cold Dark Matter

A broad range of observational data involving the dynamics of galaxies, the growth of large-scale structure, and the dynamics and nucleosynthesis of the Universe as a whole, indicate the existence of dark matter with a mean cosmic mass density that is ~ 5 times larger than the density of the baryonic matter [189, 348]. The data is consistent with a dark matter composed of weakly-interacting, massive particles, that decoupled early and adiabatically cooled to an extremely low temperature by the present time [189]. The Cold Dark Matter (CDM) has not been observed directly as of yet, although laboratory searches for particles from the dark halo of our own Milky-Way galaxy have been able to restrict the allowed parameter space for these particles. Since an alternative more-radical interpretation of the dark matter phenomenology involves a modification of gravity [253], it is of prime importance to find direct fingerprints of the CDM particles. One such fingerprint involves the small-scale structure in the Universe [158], on which we focus in this section.

The most popular candidate for the CDM particle is a Weakly Interacting Massive Particle (WIMP). The lightest supersymmetric particle (LSP) could be a WIMP (for a review see [189]). The CDM particle mass depends on free parameters in the particle physics model but typical values cover a range around M ~ 100 GeV (up to values close to a TeV). In many cases the LSP hypothesis will be tested at the Large Hadron Collider (e.g. [33]) or in direct detection experiments (e.g. [16]).

The properties of the CDM particles affect their response to the small-scale primordial inhomogeneities produced during cosmic inflation. The particle cross-section for scattering off standard model fermions sets the epoch of their thermal and kinematic decoupling from the cosmic plasma (which is significantly later than the time when their abundance freezes-out at a temperature T ~ M). Thermal decoupling is defined as the time when the temperature of the CDM stops following that of the cosmic plasma while kinematic decoupling is defined as the time when the bulk motion of the two species start to differ. For CDM the epochs of thermal and kinetic decoupling coincide. They occur when the time it takes for collisions to change the momentum of the CDM particles equals the Hubble time. The particle mass determines the thermal spread in the speeds of CDM particles, which tends to smooth-out fluctuations on very small scales due to the free-streaming of particles after kinematic decoupling [158, 159]. Viscosity has a similar effect before the CDM fluid decouples from the cosmic radiation fluid [182]. An important effect involves the memory the CDM fluid has of the acoustic oscillations of the cosmic radiation fluid out of which it decoupled. Here we consider the imprint of these acoustic oscillations on the small-scale power spectrum of density fluctuations in the Universe. Analogous imprints of acoustic oscillations of the baryons were identified recently in maps of the CMB [348], and the distribution of nearby galaxies [119]; these signatures appear on much larger scales, since the baryons decouple much later when the scale of the horizon is larger. The discussion in this section follows Loeb & Zaldarriaga (2005) [228].

Formalism

Kinematic decoupling of CDM occurs during the radiation-dominated era. For example, if the CDM is made of neutralinos with a particle mass of ~ 100 GeV, then kinematic decoupling occurs at a cosmic temperature of T_d ~ 10 MeV [182, 87]. As long as T_d << 100 MeV, we may ignore the imprint of the QCD phase transition (which transformed the cosmic quark-gluon soup into protons and neutrons) on the CDM power spectrum [321]. Over a short period of time during this transition, the pressure does not depend on density and the sound speed of the plasma vanishes, resulting in a significant growth for perturbations with periods shorter than the length of time over which the sound speed vanishes. The transition occurs when the temperature of the cosmic plasma is ~ 100-200 MeV and lasts for a small fraction of the Hubble time. As a result, the induced modifications are on scales smaller than those we are considering here and the imprint of the QCD phase transition is washed-out by the effects we calculate.

At early times the contribution of the dark matter to the energy density is negligible. Only at relatively late times when the cosmic temperature drops to values as low as ~ 1 eV, matter and radiation have comparable energy densities. As a result, the dynamics of the plasma at earlier times is virtually unaffected by the presence of the dark matter particles. In this limit, the dynamics of the radiation determines the gravitational potential and the dark matter just responds to that potential. We will use this simplification to obtain analytic estimates for the behavior of the dark matter transfer function.

The primordial inflationary fluctuations lead to acoustic modes in the radiation fluid during this era. The interaction rate of the particles in the plasma is so high that we can consider the plasma as a perfect fluid down to a comoving scale,

(18)

where _d = ₀^t_d dt / a(t) is the conformal time (i.e. the comoving size of the horizon) at the time of CDM decoupling, t_d; is the scattering cross section and n is the relevant particle density. (Throughout this section we set the speed of light and Planck's constant to unity for brevity.) The damping scale depends on the species being considered and its contribution to the energy density, and is the largest for neutrinos which are only coupled through weak interactions. In that case N ~ (T / T_d)³ where T_d ~ 1 MeV is the temperature of neutrino decoupling. At the time of CDM decoupling N ~ M / T_d ~ 10⁴ for the rest of the plasma, where M is the mass of the CDM particle. Here we will consider modes of wavelength larger than _f, and so we neglect the effect of radiation diffusion damping and treat the plasma (without the CDM) as a perfect fluid.

The equations of motion for a perfect fluid during the radiation era can be solved analytically. We will use that solution here, following the notation of Dodelson [109]. As usual we Fourier decompose fluctuations and study the behavior of each Fourier component separately. For a mode of comoving wavenumber k in Newtonian gauge, the gravitational potential fluctuations are given by:

(19)

where = k / 3^1/2 is the frequency of a mode and _p is its primordial amplitude in the limit 0. In this section we use conformal time = dt / a(t) with a(t) t^1/2 during the radiation-dominated era. Expanding the temperature anisotropy in multipole moments and using the Boltzmann equation to describe their evolution, the monopole ₀ and dipole ₁ of the photon distribution can be written in terms of the gravitational potential as [109]:

(20)

where x k and a prime denotes a derivative with respect to x.

The solutions in equations (19) and (20) assume that both the sound speed and the number of relativistic degrees of freedom are constant over time. As a result of the QCD phase transition and of various particles becoming non-relativistic, both of these assumptions are not strictly correct. The vanishing sound speed during the QCD phase transition provides the most dramatic effect, but its imprint is on scales smaller than the ones we consider here because the transition occurs at a significantly higher temperature and only lasts for a fraction of the Hubble time [321].

Before the dark matter decouples kinematically, we will treat it as a fluid which can exchange momentum with the plasma through particle collisions. At early times, the CDM fluid follows the motion of the plasma and is involved in its acoustic oscillations. The continuity and momentum equations for the CDM can be written as:

(21)

where a dot denotes an -derivative, _c is the dark matter density perturbation, _c is the divergence of the dark matter velocity field and _c denotes the anisotropic stress. In writing these equations we have followed Ref. [230]. The term _c^-1 (₁ - _c) encodes the transfer of momentun between the radiation and CDM fluids and _c^-1 provides the collisional rate of momentum transfer,

(22)

with n being the number density of particles with which the dark matter is interacting, (T) the average cross section for interaction and M the mass of the dark matter particle. The relevant scattering partners are the standard model leptons which have thermal abundances. For detailed expressions of the cross section in the case of supersymmetric (SUSY) dark matter, see Refs. [87, 159]. For our purpose, it is sufficient to specify the typical size of the cross section and its scaling with cosmic time,

(23)

where the coupling mass M is of the order of the weak-interaction scale (~ 100 GeV) for SUSY dark matter. This equation should be taken as the definition of M, as it encodes all the uncertainties in the details of the particle physics model into a single parameter. The temperature dependance of the averaged cross section is a result of the available phase space. Our results are quite insensitive to the details other than through the decoupling time. Equating _c^-1 / a to the Hubble expansion rate gives the temperature of kinematic decoupling:

(24)

The term k² c_s² _c in Eq. (21) results from the pressure gradient force and c_s is the dark matter sound speed. In the tight coupling limit, _c << H^-1 we find that c_s² f_c T / M and that the shear term is k² _c f_v c_s² _{c c}. Here f_v and f_c are constant factors of order unity. We will find that both these terms make a small difference on the scales of interest, so their precise value is unimportant.

By combining both equations in (21) into a single equation for _c we get

(25)

where x_d = k _d and _d denotes the time of kinematic decoupling which can be expressed in terms of the decoupling temperature as,

(26)

with T₀ = 2.7K being the present-day CMB temperature and z_d being the redshift at kinematic decoupling. We have also introduced the source function,

(27)

For x << x_d, the dark matter sound speed is given by

(28)

where c_s²(x_d) is the dark matter sound speed at kinematic decoupling (in units of the speed of light),

(29)

In writing (28) we have assumed that prior to decoupling the temperature of the dark matter follows that of the plasma. For the viscosity term we have,

(30)

Free streaming after kinematic decoupling

In the limit of the collision rate being much slower than the Hubble expansion, the CDM is decoupled and the evolution of its perturbations is obtained by solving a Boltzman equation:

(31)

where f(r, q, ) is the distribution function which depends on position, comoving momentum q, and time. The comoving momentum 3-components are dx_i / d = q_i / a. We use the Boltzman equation to find the evolution of modes that are well inside the horizon with x >> 1. In the radiation era, the gravitational potential decays after horizon crossing (see Eq. 19). In this limit the comoving momentum remains constant, dq_i / d =0 and the Boltzman equation becomes,

(32)

We consider a single Fourier mode and write f as,

(33)

where f₀(q) is the unperturbed distribution,

(34)

where n_CDM and T_CDM are the present-day density and temperature of the dark matter.

Our approach is to solve the Boltzman equation with initial conditions given by the fluid solution at a time _* (which will depend on k). The simplified Boltzman equation can be easily solved to give _f(q, ) as a function of the initial conditions _f(q, _*),

(35)

The CDM overdensity _c can then be expressed in terms of the perturbation in the distribution function as,

(36)

We can use (35) to obtain the evolution of _c in terms of its value at _*,

(37)

where k_f^-2 = [(T_d / M)]^1/2 _d. The exponential term is responsible for the damping of perturbations as a result of free streaming and the dispersion of the CDM particles after they decouple from the plasma. The above expression is only valid during the radiation era. The free streaming scale is simply given by dt(v / a) dta^-2 which grows logarithmically during the radiation era as in equation (37) but stops growing in the matter era when a t^2/3.

Equation (37) can be used to show that even during the free streaming epoch, _c satisfies equation (25) but with a modified sound speed and viscous term. For x >> x_d one should use,

(38)

The differences between the above scalings and those during the tight coupling regime are a result of the fact that the dark matter temperature stops following the plasma temperature but rather scales as a^-2 after thermal decoupling, which coincides with the kinematic decoupling. We ignore the effects of heat transfer during the fluid stage of the CDM because its temperature is controlled by the much larger heat reservoir of the radiation-dominated plasma at that stage.

To obtain the transfer function we solve the dark matter fluid equation until decoupling and then evolve the overdensity using equation (37) up to the time of matter - radiation equality. In practice, we use the fluid equations up to x_* = 10 max(x_d, 10) so as to switch into the free streaming solution well after the gravitational potential has decayed. In the fluid equations, we smoothly match the sound speed and viscosity terms at x = x_d. As mentioned earlier, because c_s(x_d) is so small and we are interested in modes that are comparable to the size of the horizon at decoupling, i.e. x_d ~ few, both the dark matter sound speed and the associated viscosity play only a minor role, and our simplified treatment is adequate.

In Figure 6 we illustrate the time evolution of modes during decoupling for a variety of k values. The situation is clear. Modes that enter the horizon before kinematic decoupling oscillate with the radiation fluid. This behavior has two important effects. In the absence of the coupling, modes receive a "kick" by the source term S(x) as they cross the horizon. After that they grow logarithmically. In our case, modes that entered the horizon before kinematic decoupling follow the plasma oscillations and thus miss out on both the horizon "kick" and the beginning of the logarithmic growth. Second, the decoupling from the radiation fluid is not instantaneous and this acts to further damp the amplitude of modes with x_d >> 1. This effect can be understood as follows. Once the oscillation frequency of the mode becomes high compared to the scattering rate, the coupling to the plasma effectively damps the mode. In that limit one can replace the forcing term ₀^' by its average value, which is close to zero. Thus in this regime, the scattering is forcing the amplitude of the dark matter oscillations to zero. After kinematic decoupling the modes again grow logarithmically but from a very reduced amplitude. The coupling with the plasma induces both oscillations and damping of modes that entered the horizon before kinematic decoupling. This damping is different from the free streaming damping that occurs after kinematic decoupling.

Figure 6. The normalized amplitude of CDM fluctuations / p for a variety of modes with comoving wavenumbers log(kd) = (0,1/3,2/3,1,4/3,5/3,2) as a function of x k, where = 0t dt / a(t) is the conformal time coordinate. The dashed line shows the temperature monopole 30 and the uppermost (dotted) curve shows the evolution of a mode that is uncoupled to the cosmic plasma.

Figure 7. Transfer function of the CDM density perturbation amplitude (normalized by the primordial amplitude from inflation). We show two cases: (i) Td / M = 10-4 and Td / Teq = 107; (ii) Td / M = 10-5 and Td / Teq = 107. In each case the oscillatory curve is our result and the other curve is the free-streaming only result that was derived previously in the literature [4,7,8].

In Figure 7 we show the resulting transfer function of the CDM overdensity. The transfer function is defined as the ratio between the CDM density perturbation amplitude _c when the effect of the coupling to the plasma is included and the same quantity in a model where the CDM is a perfect fluid down to arbitrarily small scales (thus, the power spectrum is obtained by multiplying the standard result by the square of the transfer function). This function shows both the oscillations and the damping signature mentioned above. The peaks occur at multipoles of the horizon scale at decoupling,

(39)

This same scale determines the "oscillation" damping. The free streaming damping scale is,

(40)

where T_eq is the temperature at matter radiation equality, T_eq 1 eV. The free streaming scale is parametrically different from the "oscillation" damping scale. However for our fiducial choice of parameters for the CDM particle they roughly coincide.

The CDM damping scale is significantly smaller than the scales observed directly in the Cosmic Microwave Background or through large scale structure surveys. For example, the ratio of the damping scale to the scale that entered the horizon at Matter-radiation equality is _d / _eq ~ T_eq / T_d ~ 10^-7 and to our present horizon _d / ₀ ~ (T_eq T₀)^1/2 / T_d ~ 10^-9. In the context of inflation, these scales were created 16 and 20 e - folds apart. Given the large extrapolation, one could certainly imagine that a change in the spectrum could alter the shape of the power spectrum around the damping scale. However, for smooth inflaton potentials with small departures from scale invariance this is not likely to be the case. On scales much smaller than the horizon at matter radiation equality, the spectrum of perturbations density before the effects of the damping are included is approximately,

(41)

where the first term encodes the shape of the primordial spectrum and the second the transfer function. Primordial departures from scale invariance are encoded in the slope n and its running . The effective slope at scale k is then,

(42)

For typical values of (n - 1) ~ 1/60 and ~ 1/60² the slope is still positive at k ~ _d^-1, so the cut-off in the power will come from the effects we calculate rather than from the shape of the primordial spectrum. However given the large extrapolation in scale, one should keep in mind the possibility of significant effects resulting from the mechanisms that generates the density perturbations.

Implications We have found that acoustic oscillations, a relic from the epoch when the dark matter coupled to the cosmic radiation fluid, truncate the CDM power spectrum on a comoving scale larger than effects considered before, such as free-streaming and viscosity [158, 159, 182]. For SUSY dark matter, the minimum mass of dark matter clumps that form in the Universe is therefore increased by more than an order of magnitude to a value of ⁴

(43)

where _crit = (H₀² / 8 G) = 9 × 10^-30 g cm^-3 is the critical density today, and _m is the matter density for the concordance cosmological model [348]. We define the cut-off wavenumber k_cut as the point where the transfer function first drops to a fraction 1 / e of its value at k 0. This corresponds to k_cut 3.3 _d^-1.

Figure 8. A slice through a numerical simulation of the first dark matter condensations to form in the Universe. Colors represent the dark matter density at z = 26. The simulated volume is 60 comoving pc on a side, simulated with 64 million particles each weighing 1.2 × 10-10 M (!). (from Diemand, Moore, & Stadel 2005 [105]).

Recent numerical simulations [105, 146] of the earliest and smallest objects to have formed in the Universe (see Fig. 3.4), need to be redone for the modified power spectrum that we calculated in this section. Although it is difficult to forecast the effects of the acoustic oscillations through the standard Press-Schechter formalism [291], it is likely that the results of such simulations will be qualitatively the same as before except that the smallest clumps would have a mass larger than before (as given by Eq. 43).

Potentially, there are several observational signatures of the smallest CDM clumps. As pointed out in the literature [105, 353], the smallest CDM clumps could produce -rays through dark-matter annihilation in their inner density cusps, with a flux in excess of that from nearby dwarf galaxies. If a substantial fraction of the Milky Way halo is composed of CDM clumps with a mass ~ 10^-4 M, the nearest clump is expected to be at a distance of ~ 4 × 10¹⁷ cm. Given that the characteristic speed of such clumps is a few hundred km s^-1, the -ray flux would therefore show temporal variations on the relatively long timescale of a thousand years. Passage of clumps through the solar system should also induce fluctuations in the detection rate of CDM particles in direct search experiments. Other observational effects have rather limited prospects for detectability. Because of their relatively low-mass and large size (~ 10¹⁷ cm), the CDM clumps are too diffuse to produce any gravitational lensing signatures (including femto-lensing [161]), even at cosmological distances.

The smallest CDM clumps should not affect the intergalactic baryons which have a much larger Jeans mass. However, once objects above ~ 10⁶ M start to collapse at redshifts z < 30, the baryons would be able to cool inside of them via molecular hydrogen transitions and the interior baryonic Jeans mass would drop. The existence of dark matter clumps could then seed the formation of the first stars inside these objects [66].

3.5. Structure of the Baryons

Early Evolution of Baryonic Perturbations on Large Scales

The baryons are coupled through Thomson scattering to the radiation fluid until they become neutral and decouple. After cosmic recombination, they start to fall into the potential wells of the dark matter and their early evolution was derived by Barkana & Loeb (2005) [29].

On large scales, the dark matter (dm) and the baryons (b) are affected only by their combined gravity and gas pressure can be ignored. The evolution of sub-horizon linear perturbations is described in the matter-dominated regime by two coupled second-order differential equations [284]:

(44)

where _dm(t) and _b(t) are the perturbations in the dark matter and baryons, respectively, the derivatives are with respect to cosmic time t, H(t) = / a is the Hubble constant with a = (1 + z)^-1, and we assume that the mean mass density _m(t) is made up of respective mass fractions f_dm and f_b = 1 - f_dm. Since these linear equations contain no spatial gradients, they can be solved spatially for _dm(x, t) and _b(x, t) or in Fourier space for _dm(k, t) and _b(k, t).

Defining _tot f_{b b} + f_{dm dm} and _{b- b} - _tot , we find

(45)

Each of these equations has two independent solutions. The equation for _tot has the usual growing and decaying solutions, which we denote D₁(t) and D₄(t), respectively, while the _b- equation has solutions D₂(t) and D₃(t); we number the solutions in order of declining growth rate (or increasing decay rate). We assume an Einstein-de Sitter, matter-dominated Universe in the redshift range z = 20 - 150, since the radiation contributes less than a few percent at z < 150, while the cosmological constant and the curvature contribute to the energy density less than a few percent at z > 3. In this regime a t^2/3 and the solutions are D₁(t) = a / a_i and D₄(t) = (a / a_i)^-3/2 for _tot, and D₂(t) = 1 and D₃(t) = (a / a_i)^-1/2 for _b-, where we have normalized each solution to unity at the starting scale factor a_i, which we set at a redshift z_i = 150. The observable baryon perturbation can then be written as

(46)

and similarly for the dark matter perturbation,

(47)

where C_i = D_i for i = 1,4 and C_i = -(f_b / f_dm)D_i for i = 2,3. We may establish the values of _m(k) by inverting the 4 × 4 matrix A that relates the 4-vector (₁, ₂, ₃, ₄) to the 4-vector that represents the initial conditions (_b, _dm, _b, _dm) at the initial time.

Figure 9. Redshift evolution of the amplitudes of the independent modes of the density perturbations (upper panel) and the temperature perturbations (lower panel), starting at redshift 150 (from Barkana & Loeb 2005 [29]). We show m = 1 (solid curves), m = 2 (short-dashed curves), m = 3 (long-dashed curves), m = 4 (dotted curves), and m = 0 (dot-dashed curve). Note that the lower panel shows one plus the mode amplitude.

Figure 10. Power spectra and initial perturbation amplitudes versus wavenumber (from [29]). The upper panel shows b (solid curves) and dm (dashed curves) at z = 150 and 20 (from bottom to top). The lower panel shows the initial (z = 150) amplitudes of 1 (solid curve), 2 (short-dashed curve), 3 (long-dashed curve), 4 (dotted curve), and T(ti) (dot-dashed curve). Note that if 1 is positive then so are 3 and T(ti), while 2 is negative at all k, and 4 is negative at the lowest k but is positive at k > 0.017 Mpc-1.

Next we describe the fluctuations in the sound speed of the cosmic gas caused by Compton heating of the gas, which is due to scattering of the residual electrons with the CMB photons. The evolution of the temperature T of a gas element of density _b is given by the first law of thermodynamics:

(48)

where dQ is the heating rate per particle. Before the first galaxies formed,

(49)

where _T is the Thomson cross-section, x_e(t) is the electron fraction out of the total number density of gas particles, and is the CMB energy density at a temperature T. In the redshift range of interest, we assume that the photon temperature (T= T⁰ / a) is spatially uniform, since the high sound speed of the photons (i.e., c / 3^1/2) suppresses fluctuations on the sub-horizon scales that we consider, and the horizon-scale ~ 10^-5 fluctuations imprinted at cosmic recombination are also negligible compared to the smallWe establish the values of _m(k) by inverting the 4 × 4 matrix A that relates the 4-vector (₁, ₂, ₃, ₄) to the 4-vector that represents the initial conditions (_b, _dm, _b, _dm) at the initial time. er-scale fluctuations in the gas density and temperature. Fluctuations in the residual electron fraction x_e(t) are even smaller. Thus,

(50)

where t^{-1 0} (8_T c / 3 m_e) = 8.55 × 10^-13 yr^-1. After cosmic recombination, x_e(t) changes due to the slow recombination rate of the residual ions:

(51)

where _B(T) is the case-B recombination coefficient of hydrogen, _H is the mean number density of hydrogen at time t, and y = 0.079 is the helium to hydrogen number density ratio. This yields the evolution of the mean temperature, d / dt = - 2 H + x_e(t) t^-1 (T - ) a^-4. In prior analyses [284, 230] a spatially uniform speed of sound was assumed for the gas at each redshift. Note that we refer to p / as the square of the sound speed of the fluid, where p is the pressure perturbation, although we are analyzing perturbations driven by gravity rather than sound waves driven by pressure gradients.

Instead of assuming a uniform sound speed, we find the first-order perturbation equation,

(52)

where we defined the fractional temperature perturbation _T. Like the density perturbation equations, this equation can be solved separately at each x or at each k. Furthermore, the solution _T (t) is a linear functional of _b(t) [for a fixed function x_e(t)]. Thus, if we choose an initial time t_i then using Eq. (46) we can write the solution in Fourier space as

(53)

where D_m^T(t) is the solution of Eq. (52) with _T = 0 at t_i and with the perturbation mode D_m(t) substituted for _b(t), while D₀^T(t) is the solution with no perturbation _b(t) and with _T = 1 at t_i. By modifying the CMBFAST code (http://www.cmbfast.org/), we can numerically solve Eq. (52) along with the density perturbation equations for each k down to z_i = 150, and then match the solution to the form of Eq. (53).

Figure 11. Relative sensitivity of perturbation amplitudes at z = 150 to cosmological parameters (from [29]). For variations in a parameter x, we show d log [P(k)]1/2 / d log(x). We consider variations in dm h2 (upper panel), in b h2 (middle panel), and in the Hubble constant h (lower panel). When we vary each parameter we fix the other two, and the variations are all carried out in a flat total = 1 universe. We show the sensitivity of 1 (solid curves), 2 (short-dashed curves), 3 (long-dashed curves), and T(ti) (dot-dashed curves).

Figure 9 shows the time evolution of the various independent modes that make up the perturbations of density and temperature, starting at the time t_i corresponding to z_i = 150. D₂^T(t) is identically zero since D₂(t) = 1 is constant, while D₃^T(t) and D₄^T(t) are negative. Figure 10 shows the amplitudes of the various components of the initial perturbations. We consider comoving wavevectors k in the range 0.01 - 40 Mpc^-1, where the lower limit is set by considering sub-horizon scales at z = 150 for which photon perturbations are negligible compared to _dm and _b, and the upper limit is set by requiring baryonic pressure to be negligible compared to gravity. ₂ and ₃ clearly show a strong signature of the large-scale baryonic oscillations, left over from the era of the photon-baryon fluid before recombination, while ₁, ₄, and _T carry only a weak sign of the oscillations. For each quantity, the plot shows [k³ P(k) / (2 ²)]^1/2, where P(k) is the corresponding power spectrum of fluctuations. ₄ is already a very small correction at z = 150 and declines quickly at lower redshift, but the other three modes all contribute significantly to _b, and the _T(t_i) term remains significant in _T(t) even at z < 100. Note that at z = 150 the temperature perturbation _T has a different shape with respect to k than the baryon perturbation _b, showing that their ratio cannot be described by a scale-independent speed of sound.

The power spectra of the various perturbation modes and of _T(t_i) depend on the initial power spectrum of density fluctuations from inflation and on the values of the fundamental cosmological parameters (_dm, _b, , and h). If these independent power spectra can be measured through 21cm fluctuations, this will probe the basic cosmological parameters through multiple combinations, allowing consistency checks that can be used to verify the adiabatic nature and the expected history of the perturbations. Figure 11 illustrates the relative sensitivity of [P(k)]^1/2 to variations in _dm h², _b h², and h, for the quantities ₁, ₂, ₃, and _T(t_i). Not shown is ₄, which although it is more sensitive (changing by order unity due to 10% variations in the parameters), its magnitude always remains much smaller than the other modes, making it much harder to detect. Note that although the angular scale of the baryon oscillations constrains also the history of dark energy through the angular diameter distance, we have focused here on other cosmological parameters, since the contribution of dark energy relative to matter becomes negligible at high redshift.

Cosmological Jeans Mass

The Jeans length _J was originally defined (Jeans 1928 [187]) in Newtonian gravity as the critical wavelength that separates oscillatory and exponentially-growing density perturbations in an infinite, uniform, and stationary distribution of gas. On scales smaller than _J, the sound crossing time, / c_s is shorter than the gravitational free-fall time, (G )^-1/2, allowing the build-up of a pressure force that counteracts gravity. On larger scales, the pressure gradient force is too slow to react to a build-up of the attractive gravitational force. The Jeans mass is defined as the mass within a sphere of radius _J / 2, M_J = (4 / 3) (_J / 2)³. In a perturbation with a mass greater than M_J, the self-gravity cannot be supported by the pressure gradient, and so the gas is unstable to gravitational collapse. The Newtonian derivation of the Jeans instability suffers from a conceptual inconsistency, as the unperturbed gravitational force of the uniform background must induce bulk motions (compare to Binney & Tremaine 1987 [43]). However, this inconsistency is remedied when the analysis is done in an expanding Universe.

The perturbative derivation of the Jeans instability criterion can be carried out in a cosmological setting by considering a sinusoidal perturbation superposed on a uniformly expanding background. Here, as in the Newtonian limit, there is a critical wavelength _J that separates oscillatory and growing modes. Although the expansion of the background slows down the exponential growth of the amplitude to a power-law growth, the fundamental concept of a minimum mass that can collapse at any given time remains the same (see, e.g. Kolb & Turner 1990 [205]; Peebles 1993 [284]).

We consider a mixture of dark matter and baryons with density parameters _dm^z = _dm / _c and _b^z = _b / _c, where _dm is the average dark matter density, _b is the average baryonic density, _c is the critical density, and _dm^z + _b^z = _m^z is given by equation(83). We also assume spatial fluctuations in the gas and dark matter densities with the form of a single spherical Fourier mode on a scale much smaller than the horizon,

(54) (55)

where _dm(t) and _b(t) are the background densities of the dark matter and baryons, _dm(t) and _b(t) are the dark matter and baryon overdensity amplitudes, r is the comoving radial coordinate, and k is the comoving perturbation wavenumber. We adopt an ideal gas equation-of-state for the baryons with a specific heat ratio = 5/3. Initially, at time t = t_i, the gas temperature is uniform T_b(r, t_i) = T_i, and the perturbation amplitudes are small _dm,i, _b,i << 1. We define the region inside the first zero of sin(kr) / (kr), namely 0 < kr < , as the collapsing "object".

The evolution of the temperature of the baryons T_b(r, t) in the linear regime is determined by the coupling of their free electrons to the CMB through Compton scattering, and by the adiabatic expansion of the gas. Hence, T_b(r, t) is generally somewhere between the CMB temperature, T (1 + z)^-1 and the adiabatically-scaled temperature T_ad (1 + z)^-2. In the limit of tight coupling to T, the gas temperature remains uniform. On the other hand, in the adiabatic limit, the temperature develops a gradient according to the relation

(56)

The evolution of a cold dark matter overdensity, _dm(t), in the linear regime is described by the equation (44),

(57)

whereas the evolution of the overdensity of the baryons, _b(t), with the inclusion of their pressure force is described by (see Section 9.3.2 of [205]),

(58)

Here, H(t) = / a is the Hubble parameter at a cosmological time t, and µ = 1.22 is the mean molecular weight of the neutral primordial gas in atomic units. The parameter distinguishes between the two limits for the evolution of the gas temperature. In the adiabatic limit = 1, and when the baryon temperature is uniform and locked to the background radiation, = 0. The last term on the right hand side (in square brackets) takes into account the extra pressure gradient force in (_b T) = (T _b + _b T), arising from the temperature gradient which develops in the adiabatic limit. The Jeans wavelength _J = 2 / k_J is obtained by setting the right-hand side of equation (58) to zero, and solving for the critical wavenumber k_J. As can be seen from equation (58), the critical wavelength _J (and therefore the mass M_J) is in general time-dependent. We infer from equation (58) that as time proceeds, perturbations with increasingly smaller initial wavelengths stop oscillating and start to grow.

To estimate the Jeans wavelength, we equate the right-hand-side of equation (58) to zero. We further approximate _b ~ _dm, and consider sufficiently high redshifts at which the Universe is matter dominated and flat, (1 + z) >> max[(1 - _m - ) / _m, ( / _m)^1/3]. In this regime, _b << _m 1, H 2 / (3t), and a = (1 + z)^-1 (3H₀ (_m)^1/2 / 2)^2/3 t^2/3, where _m = _dm + _b is the total matter density parameter. Following cosmological recombination at z 10³, the residual ionization of the cosmic gas keeps its temperature locked to the CMB temperature (via Compton scattering) down to a redshift of [284]

(59)

In the redshift range between recombination and z_t, = 0 and

(60)

so that the Jeans mass is therefore redshift independent and obtains the value (for the total mass of baryons and dark matter)

(61)

Based on the similarity of M_J to the mass of a globular cluster, Peebles & Dicke (1968) [281] suggested that globular clusters form as the first generation of baryonic objects shortly after cosmological recombination. Peebles & Dicke assumed a baryonic Universe, with a nonlinear fluctuation amplitude on small scales at z ~ 10³, a model which has by now been ruled out. The lack of a dominant mass of dark matter inside globular clusters makes it unlikely that they formed through direct cosmological collapse, and more likely that they resulted from fragmentation during the process of galaxy formation.

Figure 12. Thermal history of the baryons, left over from the big bang, before the first galaxies formed. The residual fraction of free electrons couple the gas temperture Tgas to the cosmic microwave background temperature [T (1 + z)] until a redshift z ~ 200. Subsequently the gas temperature cools adiabatically at a faster rate [Tgas (1 + z)2]. Also shown is the spin temperature of the 21cm transition of hydrogen Ts which interpolates between the gas and radiation temperature and will be discussed in detail later in this review.

At z < z_t, the gas temperature declines adiabatically as [(1 + z) / ( 1 + z_t)]² (i.e., = 1) and the total Jeans mass obtains the value,

(62)

It is not clear how the value of the Jeans mass derived above relates to the mass of collapsed, bound objects. The above analysis is perturbative (Eqs. 57 and 58 are valid only as long as _b and _dm are much smaller than unity), and thus can only describe the initial phase of the collapse. As _b and _dm grow and become larger than unity, the density profiles start to evolve and dark matter shells may cross baryonic shells [167] due to their different dynamics. Hence the amount of mass enclosed within a given baryonic shell may increase with time, until eventually the dark matter pulls the baryons with it and causes their collapse even for objects below the Jeans mass.

Even within linear theory, the Jeans mass is related only to the evolution of perturbations at a given time. When the Jeans mass itself varies with time, the overall suppression of the growth of perturbations depends on a time-weighted Jeans mass. Gnedin & Hui (1998) [150] showed that the correct time-weighted mass is the filtering mass M_f = (4 /3) (2 a / k_f)³, in terms of the comoving wavenumber k_f associated with the "filtering scale" (note the change in convention from / k_J to 2 / k_f). The wavenumber k_f is related to the Jeans wavenumber k_J by

(63)

where D(t) is the linear growth factor. At high redshift (where _m^z 1), this relation simplifies to [153]

(64)

Then the relationship between the linear overdensity of the dark matter _dm and the linear overdensity of the baryons _b, in the limit of small k, can be written as [150]

(65)

Linear theory specifies whether an initial perturbation, characterized by the parameters k, _dm,i, _b,i and t_i, begins to grow. To determine the minimum mass of nonlinear baryonic objects resulting from the shell-crossing and virialization of the dark matter, we must use a different model which examines the response of the gas to the gravitational potential of a virialized dark matter halo.

3.6. Formation of Nonlinear Objects

3.7. Spherical Collapse

Let us consider a spherically symmetric density or velocity perturbation of the smooth cosmological background, and examine the dynamics of a test particle at a radius r relative to the center of symmetry. Birkhoff's (1923) [44] theorem implies that we may ignore the mass outside this radius in computing the motion of our particle. We further find that the relativistic equations of motion describing the system reduce to the usual Friedmann equation for the evolution of the scale factor of a homogeneous Universe, but with a density parameter that now takes account of the additional mass or peculiar velocity. In particular, despite the arbitrary density and velocity profiles given to the perturbation, only the total mass interior to the particle's radius and the peculiar velocity at the particle's radius contribute to the effective value of . We thus find a solution to the particle's motion which describes its departure from the background Hubble flow and its subsequent collapse or expansion. This solution holds until our particle crosses paths with one from a different radius, which happens rather late for most initial profiles.

As with the Friedmann equation for a smooth Universe, it is possible to reinterpret the problem into a Newtonian form. Here we work in an inertial (i.e. non-comoving) coordinate system and consider the force on the particle as that resulting from a point mass at the origin (ignoring the possible presence of a vacuum energy density):

(66)

where G is Newton's constant, r is the distance of the particle from the center of the spherical perturbation, and M is the total mass within that radius. As long as the radial shells do not cross each other, the mass M is constant in time. The initial density profile determines M, while the initial velocity profile determines dr / dt at the initial time. As is well-known, there are three branches of solutions: one in which the particle turns around and collapses, another in which it reaches an infinite radius with some asymptotically positive velocity, and a third intermediate case in which it reachs an infinite radius but with a velocity that approaches zero. These cases may be written as [164]:

(67)

(68)

(69)

where A³ = GMB² applies in all cases. All three solutions have r³ = 9GMt² / 2 as t goes to zero, which matches the linear theory expectation that the perturbation amplitude get smaller as one goes back in time. In the closed case, the shell turns around at time B and radius 2A and collapses to zero radius at time 2 B.

We are now faced with the problem of relating the spherical collapse parameters A, B, and M to the linear theory density perturbation [283]. We do this by returning to the equation of motion. Consider that at an early epoch (i.e. scale factor a_i << 1), we are given a spherical patch of uniform overdensity _i (the so-called `top-hat' perturbation). If is essentially unity at this time and if the perturbation is pure growing mode, then the initial velocity is radially inward with magnitude _i H(t_i)r / 3, where H(t_i) is the Hubble constant at the initial time and r is the radius from the center of the sphere. This can be easily seen from the continuity equation in spherical coordinates. The equation of motion (in noncomoving coordinates) for a particle beginning at radius r_i is simply

(70)

where M = (4 / 3) r_i³ _i (1 + _i) and _i is the background density of the Universe at time t_i. We next define the dimensionless radius x = ra_i / r_i and rewrite equation (70) as

(71)

Our initial conditions for the integration of this orbit are

(72)

(73)

where H(t₁) = H₀[_m / a³(t₁) + (1 - _m)]^1/2 is the Hubble parameter for a flat Universe at a a cosmic time t₁. Integrating equation (71) yields

(74)

where K is a constant of integration. Evaluating this at the initial time and dropping terms of O(a_i) (but _i ~ a_i, so we keep ratios of order unity), we find

(75)

If K is sufficiently negative, the particle will turn-around and the sphere will collapse at a time

(76)

where a_max is the value of a which sets the denominator of the integral to zero.

For the case of = 0, we can determine the spherical collapse parameters A and B. K > 0 (K < 0) produces an open (closed) model. Comparing coefficients in the energy equations [eq. (74) and the integration of (66)], one finds

(77)

(78)

where _k = 1 - _m. In particular, in an = 1 Universe, where 1 + z = (3H₀ t / 2)^-2/3, we find that a shell collapses at redshift 1 + z_c = 0.5929_i / a_i, or in other words a shell collapsing at redshift z_c had a linear overdensity extrapolated to the present day of ₀ = 1.686(1 + z_c).

While this derivation has been for spheres of constant density, we may treat a general spherical density profile _i(r) up until shell crossing [164]. A particular radial shell evolves according to the mass interior to it; therefore, we define the average overdensity

(79)

so that we may use in place of _i in the above formulae. If is not monotonically decreasing with R, then the spherical top-hat evolution of two different radii will predict that they cross each other at some late time; this is known as shell crossing and signals the breakdown of the solution. Even well-behaved profiles will produce shell crossing if shells are allowed to collapse to r = 0 and then reexpand, since these expanding shells will cross infalling shells. In such a case, first-time infalling shells will never be affected prior to their turn-around; the more complicated behavior after turn-around is a manifestation of virialization. While the end state for general initial conditions cannot be predicted, various results are known for a self-similar collapse, in which (r) is a power-law [132, 40], as well as for the case of secondary infall models [156, 165, 181].

3.8. Halo Properties

The small density fluctuations evidenced in the CMB grow over time as described in the previous subsection, until the perturbation becomes of order unity, and the full non-linear gravitational problem must be considered. The dynamical collapse of a dark matter halo can be solved analytically only in cases of particular symmetry. If we consider a region which is much smaller than the horizon cH^-1, then the formation of a halo can be formulated as a problem in Newtonian gravity, in some cases with minor corrections coming from General Relativity. The simplest case is that of spherical symmetry, with an initial (t = t_i << t₀) top-hat of uniform overdensity _i inside a sphere of radius R. Although this model is restricted in its direct applicability, the results of spherical collapse have turned out to be surprisingly useful in understanding the properties and distribution of halos in models based on cold dark matter.

The collapse of a spherical top-hat perturbation is described by the Newtonian equation (with a correction for the cosmological constant)

(80)

where r is the radius in a fixed (not comoving) coordinate frame, H₀ is the present-day Hubble constant, M is the total mass enclosed within radius r, and the initial velocity field is given by the Hubble flow dr / dt = H(t) r. The enclosed grows initially as _L = _i D(t) / D(t_i), in accordance with linear theory, but eventually grows above _L. If the mass shell at radius r is bound (i.e., if its total Newtonian energy is negative) then it reaches a radius of maximum expansion and subsequently collapses. As demonstrated in the previous section, at the moment when the top-hat collapses to a point, the overdensity predicted by linear theory is _L = 1.686 in the Einstein-de Sitter model, with only a weak dependence on _m and . Thus a tophat collapses at redshift z if its linear overdensity extrapolated to the present day (also termed the critical density of collapse) is

(81)

where we set D(z = 0) = 1.

Even a slight violation of the exact symmetry of the initial perturbation can prevent the tophat from collapsing to a point. Instead, the halo reaches a state of virial equilibrium by violent relaxation (phase mixing). Using the virial theorem U = -2K to relate the potential energy U to the kinetic energy K in the final state (implying that the virial radius is half the turnaround radius - where the kinetic energy vanishes), the final overdensity relative to the critical density at the collapse redshift is _c = 18² 178 in the Einstein-de Sitter model, modified in a Universe with _m + = 1 to the fitting formula (Bryan & Norman 1998 [71])

(82)

where d _m^z-1 is evaluated at the collapse redshift, so that

(83)

A halo of mass M collapsing at redshift z thus has a virial radius

(84)

and a corresponding circular velocity,

(85)

In these expressions we have assumed a present Hubble constant written in the form H₀ = 100 h km s^-1Mpc^-1. We may also define a virial temperature

(86)

where µ is the mean molecular weight and m_p is the proton mass. Note that the value of µ depends on the ionization fraction of the gas; for a fully ionized primordial gas µ = 0.59, while a gas with ionized hydrogen but only singly-ionized helium has µ = 0.61. The binding energy of the halo is approximately ⁵

(87)

Note that the binding energy of the baryons is smaller by a factor equal to the baryon fraction _b / _m.

Although spherical collapse captures some of the physics governing the formation of halos, structure formation in cold dark matter models procedes hierarchically. At early times, most of the dark matter is in low-mass halos, and these halos continuously accrete and merge to form high-mass halos. Numerical simulations of hierarchical halo formation indicate a roughly universal spherically-averaged density profile for the resulting halos (Navarro, Frenk, & White 1997, hereafter NFW [266]), though with considerable scatter among different halos (e.g., [72]). The NFW profile has the form

(88)

where x = r / r_vir, and the characteristic density _c is related to the concentration parameter c_N by

(89)

The concentration parameter itself depends on the halo mass M, at a given redshift z [377].

More recent N-body simulations indicate deviations from the original NFW profile; for details and refined fitting formula see [268].

4. NONLINEAR GROWTH

4.1. The Abundance of Dark Matter Halos

In addition to characterizing the properties of individual halos, a critical prediction of any theory of structure formation is the abundance of halos, i.e. the number density of halos as a function of mass, at any redshift. This prediction is an important step toward inferring the abundances of galaxies and galaxy clusters. While the number density of halos can be measured for particular cosmologies in numerical simulations, an analytic model helps us gain physical understanding and can be used to explore the dependence of abundances on all the cosmological parameters.

A simple analytic model which successfully matches most of the numerical simulations was developed by Press & Schechter (1974) [291]. The model is based on the ideas of a Gaussian random field of density perturbations, linear gravitational growth, and spherical collapse. To determine the abundance of halos at a redshift z, we use _m, the density field smoothed on a mass scale M, as defined in Section 3.3. Since _m is distributed as a Gaussian variable with zero mean and standard deviation (M) [which depends only on the present linear power spectrum, see equation (17)], the probability that _m is greater than some equals

(90)

The fundamental ansatz is to identify this probability with the fraction of dark matter particles which are part of collapsed halos of mass greater than M, at redshift z. There are two additional ingredients: First, the value used for is _crit(z) given in equation (81), which is the critical density of collapse found for a spherical top-hat (extrapolated to the present since (M) is calculated using the present power spectrum); and second, the fraction of dark matter in halos above M is multiplied by an additional factor of 2 in order to ensure that every particle ends up as part of some halo with M > 0. Thus, the final formula for the mass fraction in halos above M at redshift z is

(91)

This ad-hoc factor of 2 is necessary, since otherwise only positive fluctuations of _m would be included. Bond et al. (1991) [52] found an alternate derivation of this correction factor, using a different ansatz. In their derivation, the factor of 2 has a more satisfactory origin, namely the so-called "cloud-in-cloud" problem: For a given mass M, even if _m is smaller than _crit(z), it is possible that the corresponding region lies inside a region of some larger mass M_L > M, with _ML > _crit(z). In this case the original region should be counted as belonging to a halo of mass M_L. Thus, the fraction of particles which are part of collapsed halos of mass greater than M is larger than the expression given in equation (90). Bond et al. showed that, under certain assumptions, the additional contribution results precisely in a factor of 2 correction.

Figure 13. Fraction of baryons that assembled into dark matter halos with a virial temperature of Tvir > 104K as a function of redshift. These baryons are above the temperature threshold for gas cooling and fragmentation via atomic transitions. After reionization the temperature barrier for star formation in galaxies is raised because the photo-ionized intergalactic medium is already heated to ~ 104K and it can condense only into halos with Tvir > 105K.

Differentiating the fraction of dark matter in halos above M yields the mass distribution. Letting dn be the comoving number density of halos of mass between M and M + dM, we have

(92)

where _c = _crit(z) / (M) is the number of standard deviations which the critical collapse overdensity represents on mass scale M. Thus, the abundance of halos depends on the two functions (M) and _crit (z), each of which depends on the energy content of the Universe and the values of the other cosmological parameters. Since recent observations confine the standard set of parameters to a relatively narrow range, we illustrate the abundance of halos and other results for a single set of parameters: _m = 0.3, = 0.7, _b = 0.045, ₈ = 0.9, a primordial power spectrum index n = 1 and a Hubble constant h = 0.7.

Figure 14 shows (M) and _crit(z), with the input power spectrum computed from Eisenstein & Hu (1999) [118]. The solid line is (M) for the cold dark matter model with the parameters specified above. The horizontal dotted lines show the value of _crit(z) at z = 0, 2, 5, 10, 20 and 30, as indicated in the figure. From the intersection of these horizontal lines with the solid line we infer, e.g., that at z = 5 a 1- fluctuation on a mass scale of 2 × 10⁷ M will collapse. On the other hand, at z = 5 collapsing halos require a 2- fluctuation on a mass scale of 3 × 10¹⁰ M, since (M) on this mass scale equals about half of _crit(z = 5). Since at each redshift a fixed fraction (31.7%) of the total dark matter mass lies in halos above the 1- mass, Figure 14 shows that most of the mass is in small halos at high redshift, but it continuously shifts toward higher characteristic halo masses at lower redshift. Note also that (M) flattens at low masses because of the changing shape of the power spectrum. Since as M 0, in the cold dark matter model all the dark matter is tied up in halos at all redshifts, if sufficiently low-mass halos are considered.

Figure 14. Mass fluctuations and collapse thresholds in cold dark matter models. The horizontal dotted lines show the value of the extrapolated collapse overdensity crit(z) at the indicated redshifts. Also shown is the value of (M) for the cosmological parameters given in the text (solid curve), as well as (M) for a power spectrum with a cutoff below a mass M = 1.7 × 108 M (short-dashed curve), or M = 1.7 × 1011 M (long-dashed curve). The intersection of the horizontal lines with the other curves indicate, at each redshift z, the mass scale (for each model) at which a 1- fluctuation is just collapsing at z (see the discussion in the text).

Also shown in Figure 14 is the effect of cutting off the power spectrum on small scales. The short-dashed curve corresponds to the case where the power spectrum is set to zero above a comoving wavenumber k = 10 Mpc^-1, which corresponds to a mass M = 1.7 × 10⁸ M. The long-dashed curve corresponds to a more radical cutoff above k = 1 Mpc^-1, or below M = 1.7 × 10¹¹ M. A cutoff severely reduces the abundance of low-mass halos, and the finite value of (M = 0) implies that at all redshifts some fraction of the dark matter does not fall into halos. At high redshifts where _crit(z) >> (M = 0), all halos are rare and only a small fraction of the dark matter lies in halos. In particular, this can affect the abundance of halos at the time of reionization, and thus the observed limits on reionization constrain scenarios which include a small-scale cutoff in the power spectrum [21].

In figures 15 - 18 we show explicitly the properties of collapsing halos which represent 1-, 2-, and 3- fluctuations (corresponding in all cases to the curves in order from bottom to top), as a function of redshift. No cutoff is applied to the power spectrum. Figure 15 shows the halo mass, Figure 16 the virial radius, Figure 17 the virial temperature (with µ in equation (86) set equal to 0.6, although low temperature halos contain neutral gas) as well as circular velocity, and Figure 18 shows the total binding energy of these halos. In figures 15 and 17, the dotted curves indicate the minimum virial temperature required for efficient cooling with primordial atomic species only (upper curve) or with the addition of molecular hydrogen (lower curve). Figure 18 shows the binding energy of dark matter halos. The binding energy of the baryons is a factor ~ _b / _m ~ 15% smaller, if they follow the dark matter. Except for this constant factor, the figure shows the minimum amount of energy that needs to be deposited into the gas in order to unbind it from the potential well of the dark matter. For example, the hydrodynamic energy released by a single supernovae, ~ 10⁵¹ erg, is sufficient to unbind the gas in all 1- halos at z > 5 and in all 2- halos at z > 12.

Figure 15. Characteristic properties of collapsing halos: Halo mass. The solid curves show the mass of collapsing halos which correspond to 1-, 2-, and 3- fluctuations (in order from bottom to top). The dotted curves show the mass corresponding to the minimum temperature required for efficient cooling with primordial atomic species only (upper curve) or with the addition of molecular hydrogen (lower curve).

Figure 16. Characteristic properties of collapsing halos: Halo virial radius. The curves show the virial radius of collapsing halos which correspond to 1-, 2-, and 3- fluctuations (in order from bottom to top).

Figure 17. Characteristic properties of collapsing halos: Halo virial temperature and circular velocity. The solid curves show the virial temperature (or, equivalently, the circular velocity) of collapsing halos which correspond to 1-, 2-, and 3- fluctuations (in order from bottom to top). The dotted curves show the minimum temperature required for efficient cooling with primordial atomic species only (upper curve) or with the addition of molecular hydrogen (lower curve).

Figure 18. Characteristic properties of collapsing halos: Halo binding energy. The curves show the total binding energy of collapsing halos which correspond to 1-, 2-, and 3- fluctuations (in order from bottom to top).

In Figure 19 we show the halo mass function dn / dln(M) at several different redshifts: z = 0 (solid curve), z = 5 (dotted curve), z = 10 (short-dashed curve), z = 20 (long-dashed curve), and z = 30 (dot-dashed curve). Note that the mass function does not decrease monotonically with redshift at all masses. At the lowest masses, the abundance of halos is higher at z > 0 than at z = 0.

Figure 19. Halo mass function at several redshifts: z = 0 (solid curve), z = 5 (dotted curve), z = 10 (short-dashed curve), z = 20 (long-dashed curve), and z = 30 (dot-dashed curve).

4.2. The Excursion-Set (Extended Press-Schechter) Formalism

The usual Press-Schechter formalism makes no attempt to deal with the correlations between halos or between different mass scales. In particular, this means that while it can generate a distribution of halos at two different epochs, it says nothing about how particular halos in one epoch are related to those in the second. We therefore would like some method to predict, at least statistically, the growth of individual halos via accretion and mergers. Even restricting ourselves to spherical collapse, such a model must utilize the full spherically-averaged density profile around a particular point. The potential correlations between the mean overdensities at different radii make the statistical description substantially more difficult.

The excursion set formalism (Bond et al. 1991 [52]) seeks to describe the statistics of halos by considering the statistical properties of (R), the average overdensity within some spherical window of characteristic radius R, as a function of R. While the Press-Schechter model depends only on the Gaussian distribution of for one particular R, the excursion set considers all R. Again the connection between a value of the linear regime and the final state is made via the spherical collapse solution, so that there is a critical value _c (z) of which is required for collapse at a redshift z.

For most choices of window function, the functions (R) are correlated from one R to another such that it is prohibitively difficult to calculate the desired statistics directly [although Monte Carlo realizations are possible [52]]. However, for one particular choice of a window function, the correlations between different R greatly simplify and many interesting quantities may be calculated [52, 212]. The key is to use a k-space top-hat window function, namely W_k = 1 for all k less than some critical k_c and W_k = 0 for k > k_c. This filter has a spatial form of W(r) j₁ (k_c r) / k_c r, which implies a volume 6² / k_c³ or mass 6² _b / k_c³. The characteristic radius of the filter is ~ k_c^-1, as expected. Note that in real space, this window function converges very slowly, due only to a sinusoidal oscillation, so the region under study is rather poorly localized.

The great advantage of the sharp k-space filter is that the difference at a given point between on one mass scale and that on another mass scale is statistically independent from the value on the larger mass scale. With a Gaussian random field, each _k is Gaussian distributed independently from the others. For this filter,

(93)

meaning that the overdensity on a particular scale is simply the sum of the random variables _k interior to the chosen k_c. Consequently, the difference between the (M) on two mass scales is just the sum of the _k in the spherical shell between the two k_c, which is independent from the sum of the _k interior to the smaller k_c. Meanwhile, the distribution of (M) given no prior information is still a Gaussian of mean zero and variance

(94)

If we now consider as a function of scale k_c, we see that we begin from = 0 at k_c = 0 (M = ) and then add independently random pieces as k_c increases. This generates a random walk, albeit one whose stepsize varies with k_c. We then assume that, at redshift z, a given function (k_c ) represents a collapsed mass M corresponding to the k_c where the function first crosses the critical value _c (z). With this assumption, we may use the properties of random walks to calculate the evolution of the mass as a function of redshift.

It is now easy to rederive the Press-Schechter mass function, including the previously unexplained factor of 2 [52, 212, 382]. The fraction of mass elements included in halos of mass less than M is just the probability that a random walk remains below _c (z) for all k_c less than K_c, the filter cutoff appropriate to M. This probability must be the complement of the sum of the probabilities that: (a) (K_c) > _c(z); or that (b) (K_c) < _c(z) but (k'_c) > _c(z) for some k'_c < K_c. But these two cases in fact have equal probability; any random walk belonging to class (a) may be reflected around its first upcrossing of _c(z) to produce a walk of class (b), and vice versa. Since the distribution of (K_c) is simply Gaussian with variance ²(M), the fraction of random walks falling into class (a) is simply (1 / (2 ²)^1/2) _{c (z)} d exp{- ² / 2² (M)}. Hence, the fraction of mass elements included in halos of mass less than M at redshift z is simply

(95)

which may be differentiated to yield the Press-Schechter mass function. We may now go further and consider how halos at one redshift are related to those at another redshift. If we are given that a halo of mass M₂ exists at redshift z₂, then we know that the random function (k_c) for each mass element within the halo first crosses (z₂) at k_c2 corresponding to M₂. Given this constraint, we may study the distribution of k_c where the function (k_c) crosses other thresholds. It is particularly easy to construct the probability distribution for when trajectories first cross some _c (z₁) > _c (z₂) (implying z₁ > z₂); clearly this occurs at some k_c1 > k_c2. This problem reduces to the previous one if we translate the origin of the random walks from (k_c, ) = (0,0) to (k_c2, _c(z₂)). We therefore find the distribution of halo masses M₁ that a mass element finds itself in at redshift z₁ given that it is part of a larger halo of mass M₂ at a later redshift z₂ is [52, 55])

(96)

This may be rewritten as saying that the quantity

(97)

is distributed as the positive half of a Gaussian with unit variance; equation (97) may be inverted to find M₁().

We seek to interpret the statistics of these random walks as those of merging and accreting halos. For a single halo, we may imagine that as we look back in time, the object breaks into ever smaller pieces, similar to the branching of a tree. Equation (96) is the distribution of the sizes of these branches at some given earlier time. However, using this description of the ensemble distribution to generate Monte Carlo realizations of single merger trees has proven to be difficult. In all cases, one recursively steps back in time, at each step breaking the final object into two or more pieces. An elaborate scheme (Kauffmann & White 1993 [195]) picks a large number of progenitors from the ensemble distribution and then randomly groups them into sets with the correct total mass. This generates many (hundreds) possible branching schemes of equal likelihood. A simpler scheme (Lacey & Cole 1993 [212]) assumes that at each time step, the object breaks into two pieces. One value from the distribution (96) then determines the mass ratio of the two branchs.

One may also use the distribution of the ensemble to derive some additional analytic results. A useful example is the distribution of the epoch at which an object that has mass M₂ at redshift z₂ has accumulated half of its mass [212]. The probability that the formation time is earlier than z₁ is equal to the probability that at redshift z₁ a progenitor whose mass exceeds M₂/2 exists:

(98)

where dP / dM is given in equation (96). The factor of M₂ / M corrects the counting from mass weighted to number weighted; each halo of mass M₂ can have only one progenitor of mass greater than M₂ / 2. Differentiating equation (98) with respect to time gives the distribution of formation times. This analytic form is an excellent match to scale-free N-body simulations [213]. On the other hand, simple Monte Carlo implementations of equation (96) produce formation redshifts about 40% higher [212]. As there may be correlations between the various branches, there is no unique Monte Carlo scheme.

Numerical tests of the excursion set formalism are quite encouraging. Its predictions for merger rates are in very good agreement with those measured in scale-free N-body simulations for mass scales down to around 10% of the nonlinear mass scale (that scale at which _m = 1), and distributions of formation times closely match the analytic predictions [213]. The model appears to be a promising method for tracking the merging of halos, with many applications to cluster and galaxy formation modeling. In particular, one may use the formalism as the foundation of semi-analytic galaxy formation models [196]. The excursion set formalism may also be used to derive the correlations of halos in the nonlinear regime [258].

4.3. Response of Baryons to Nonlinear Dark Matter Potentials

The dark matter is assumed to be cold and to dominate gravity, and so its collapse and virialization proceeds unimpeded by pressure effects. In order to estimate the minimum mass of baryonic objects, we must go beyond linear perturbation theory and examine the baryonic mass that can accrete into the final gravitational potential well of the dark matter.

For this purpose, we assume that the dark matter had already virialized and produced a gravitational potential (r) at a redshift z_vir (with 0 at large distances, and < 0 inside the object) and calculate the resulting overdensity in the gas distribution, ignoring cooling (an assumption justified by spherical collapse simulations which indicate that cooling becomes important only after virialization; see Haiman et al. 1996 [168]).

After the gas settles into the dark matter potential well, it satisfies the hydrostatic equilibrium equation,

(99)

where p_b and _b are the pressure and mass density of the gas. At z < 100 the gas temperature is decoupled from the CMB, and its pressure evolves adiabatically (ignoring atomic or molecular cooling),

(100)

where a bar denotes the background conditions. We substitute equation (100) into (99) and get the solution,

(101)

where = _b µ m_p / (k _b) is the background gas temperature. If we define T_vir= -1/3m_p / k as the virial temperature for a potential depth -, then the overdensity of the baryons at the virialization redshift is

(102)

This solution is approximate for two reasons: (i) we assumed that the gas is stationary throughout the entire region and ignored the transitions to infall and the Hubble expansion at the interface between the collapsed object and the background intergalactic medium (henceforth IGM), and (ii) we ignored entropy production at the virialization shock surrounding the object. Nevertheless, the result should provide a better estimate for the minimum mass of collapsed baryonic objects than the Jeans mass does, since it incorporates the nonlinear potential of the dark matter.

We may define the threshold for the collapse of baryons by the criterion that their mean overdensity, _b, exceeds a value of 100, amounting to > 50% of the baryons that would assemble in the absence of gas pressure, according to the spherical top-hat collapse model. Equation (102) then implies that T_vir > 17.2 .

As mentioned before, the gas temperature evolves at z < 160 according to the relation 170 [(1 + z) / 100]² K. This implies that baryons are overdense by _b > 100 only inside halos with a virial temperature T_vir > 2.9 × 10³ [(1 + z) / 100]² K. Based on the top-hat model, this implies a minimum halo mass for baryonic objects of

(103)

where we consider sufficiently high redshifts so that _m^z 1. This minimum mass is coincidentally almost identical to the naive Jeans mass calculation of linear theory in equation (62) despite the fact that it incorporates shell crossing by the dark matter, which is not accounted for by linear theory. Unlike the Jeans mass, the minimum mass depends on the choice for an overdensity threshold [taken arbitrarily as _b > 100 in equation (103)]. To estimate the minimum halo mass which produces any significant accretion we set, e.g., _b = 5, and get a mass which is lower than M_min by a factor of 27.

Of course, once the first stars and quasars form they heat the surrounding IGM by either outflows or radiation. As a result, the Jeans mass which is relevant for the formation of new objects changes [148, 152]). The most dramatic change occurs when the IGM is photo-ionized and is consequently heated to a temperature of ~ (1 - 2) × 10⁴ K.

5. FRAGMENTATION OF THE FIRST GASEOUS OBJECTS TO STARS

5.1. Star Formation

As mentioned in the preface, the fragmentation of the first gaseous objects is a well-posed physics problem with well specified initial conditions, for a given power-spectrum of primordial density fluctuations. This problem is ideally suited for three-dimensional computer simulations, since it cannot be reliably addressed in idealized 1D or 2D geometries.

Recently, two groups have attempted detailed 3D simulations of the formation process of the first stars in a halo of ~ 10⁶ M by following the dynamics of both the dark matter and the gas components, including H₂ chemistry and cooling. Bromm, Coppi, & Larson (1999) [57] have used a Smooth Particle Hydrodynamics (SPH) code to simulate the collapse of a top-hat overdensity with a prescribed solid-body rotation (corresponding to a spin parameter = 5%) and additional small perturbations with P(k) k^-3 added to the top-hat profile. Abel et al. (2002) [5] isolated a high-density filament out of a larger simulated cosmological volume and followed the evolution of its density maximum with exceedingly high resolution using an Adaptive Mesh Refinement (AMR) algorithm.

Figure 20. Cooling rates as a function of temperature for a primordial gas composed of atomic hydrogen and helium, as well as molecular hydrogen, in the absence of any external radiation. We assume a hydrogen number density nH = 0.045 cm-3, corresponding to the mean density of virialized halos at z = 10. The plotted quantity / nH2 is roughly independent of density (unless nH > 10 cm-3), where is the volume cooling rate (in erg/sec/cm3). The solid line shows the cooling curve for an atomic gas, with the characteristic peaks due to collisional excitation of H1 and He2. The dashed line shows the additional contribution of molecular cooling, assuming a molecular abundance equal to 1% of nH.

The generic results of Bromm et al. (1999 [57]; see also Bromm 2000 [58]) are illustrated in Figure 21. The collapsing region forms a disk which fragments into many clumps. The clumps have a typical mass ~ 10² - 10³ M. This mass scale corresponds to the Jeans mass for a temperature of ~ 500K and the density ~ 10⁴ cm^-3 where the gas lingers because its cooling time is longer than its collapse time at that point (see Fig. 22). Each clump accretes mass slowly until it exceeds the Jeans mass and collapses at a roughly constant temperature (isothermally) due to H₂ cooling that brings the gas to a fixed temperature floor. The clump formation efficiency is high in this simulation due to the synchronized collapse of the overall top-hat perturbation.

Figure 21. Numerical results from Bromm et al. (1999) [57], showing gas properties at z = 31.2 for a collapsing slightly inhomogeneous top-hat region with a prescribed solid-body rotation. (a) Free electron fraction (by number) vs. hydrogen number density (in cm-3). At densities exceeding n ~ 103 cm-3, recombination is very efficient, and the gas becomes almost completely neutral. (b) Molecular hydrogen fraction vs. number density. After a quick initial rise, the H2 fraction approaches the asymptotic value of f ~ 10-3, due to the H- channel. (c) Gas temperature vs. number density. At densities below ~ 1 cm-3, the gas temperature rises because of adiabatic compression until it reaches the virial value of Tvir 5000 K. At higher densities, cooling due to H2 drives the temperature down again, until the gas settles into a quasi-hydrostatic state at T ~ 500 K and n ~ 104 cm-3. Upon further compression due to accretion and the onset of gravitational collapse, the gas shows a further modest rise in temperature. (d) Jeans mass (in M) vs. number density. The Jeans mass reaches a value of MJ ~ 103 M for the quasi-hydrostatic gas in the center of the potential well, and reaches the resolution limit of the simulation, Mres 200 M, for densities close to n = 108 cm-3.

Figure 22. Gas and clump morphology at z = 28.9 in the simulation of Bromm et al. (1999) [57]. Top row: The remaining gas in the diffuse phase. Bottom row: Distribution of clumps. The numbers next to the dots denote clump mass in units of M. Left panels: Face-on view. Right panels: Edge-on view. The length of the box is 30 pc. The gas has settled into a flattened configuration with two dominant clumps of mass close to 20,000 M. During the subsequent evolution, the clumps survive without merging, and grow in mass only slightly by accretion of surrounding gas.

Bromm (2000) [58] has simulated the collapse of one of the above-mentioned clumps with ~ 1000 M and demonstrated that it does not tend to fragment into sub-components. Rather, the clump core of ~ 100M free-falls towards the center leaving an extended envelope behind with a roughly isothermal density profile. At very high gas densities, three-body reactions become important in the chemistry of H₂. Omukai & Nishi (1998) [274] have included these reactions as well as radiative transfer and followed the collapse in spherical symmetry up to stellar densities. Radiation pressure from nuclear burning at the center is unlikely to reverse the infall as the stellar mass builds up. These calculations indicate that each clump may end as a single massive star; however, it is conceivable that angular momentum may eventually halt the collapsing cloud and lead to the formation of a binary stellar system instead.

The Jeans mass, which is defined based on small fluctuations in a background of uniform density, does not strictly apply in the context of collapsing gas cores. We can instead use a slightly modified critical mass known as the Bonnor-Ebert mass [53, 114]. For baryons in a background of uniform density _b, perturbations are unstable to gravitational collapse in a region more massive than the Jeans mass. Instead of a uniform background, we consider a spherical, non-singular, isothermal, self-gravitating gas in hydrostatic equilibrium, i.e., a centrally-concentrated object which more closely resembles the gas cores found in the above-mentioned simulations. In this case, small fluctuations are unstable and lead to collapse if the sphere is more massive than the Bonnor-Ebert mass M_BE, given by the same expression the Jeans Mass but with a different coefficient (1.2 instead of 2.9) and with _b denoting in this case the gas (volume) density at the surface of the sphere,

(104)

In their simulation, Abel et al. (2000) [4] adopted the actual cosmological density perturbations as initial conditions. The simulation focused on the density peak of a filament within the IGM, and evolved it to very high densities (Fig. 23). Following the initial collapse of the filament, a clump core formed with ~ 200 M, amounting to only ~ 1% of the virialized mass. Subsequently due to slow cooling, the clump collapsed subsonically in a state close to hydrostatic equilibrium (see Fig. 24). Unlike the idealized top-hat simulation of Bromm et al. (2001) [59], the collapse of the different clumps within the filament is not synchronized. Once the first star forms at the center of the first collapsing clump, it is likely to affect the formation of other stars in its vicinity.

As soon as nuclear burning sets in the core of the proto-star, the radiation emitted by the star starts to affect the infall of the surrounding gas towards it. The radiative feedback involves photo-dissociation of H₂, Ly radiation pressure, and photo-evaporation of the accretion disk. Tan & McKee [357] studied these effects by extrapolating analytically the infall of gas from the final snapshot of the above resolution-limited simulations to the scale of a proto-star; they concluded that nuclear burning (and hence the feedback) starts when the proton-star accretes ~ 30 M and accretion is likely to be terminated when the star reaches ~ 200 M.

Figure 23. Zooming in on the core of a star forming region with the Adaptive Mesh Refinement simulation of Abel et al. (2000) [4]. The panels show different length scales, decreasing clockwise by an order of magnitude between adjacent panels. Note the large dynamic range of scales which are being resolved, from 6 kpc (top left panel) down to 10,000 AU (bottom left panel).

Figure 24. Gas profiles from the simulation of Abel et al. (2000) [4]. The cell size on the finest grid corresponds to 0.024 pc, while the simulation box size corresponds to 6.4 kpc. Shown are spherically-averaged mass-weighted profiles around the baryon density peak shortly before a well defined fragment forms (z = 19.1). Panel (a) shows the baryonic number density, enclosed gas mass in solar mass, and the local Bonnor-Ebert mass MBE (see text). Panel (b) plots the molecular hydrogen fraction (by number) fH2 and the free electron fraction x. The H2 cooling time, tH2, the time it takes a sound wave to travel to the center, tcross, and the free - fall time tff = [3 / (32G )]1/2 are given in panel (c). Panel (d) gives the temperature in K as a function of radius. The bottom panel gives the local sound speed, cs (solid line with circles), the rms radial velocities of the dark matter (dashed line) and the gas (dashed line with asterisks) as well as the rms gas velocity (solid line with square symbols). The vertical dotted line indicates the radius (~ 5 pc) at which the gas has reached its minimum temperature allowed by H2 cooling. The virial radius of the 5.6 × 106 M halo is 106 pc.

If the clumps in the above simulations end up forming individual very massive stars, then these stars will likely radiate copious amounts of ionizing radiation [50, 370, 59] and expel strong winds. Hence, the stars will have a large effect on their interstellar environment, and feedback is likely to control the overall star formation efficiency. This efficiency is likely to be small in galactic potential wells which have a virial temperature lower than the temperature of photoionized gas, ~ 10⁴K. In such potential wells, the gas may go through only a single generation of star formation, leading to a "suicidal" population of massive stars.

The final state in the evolution of these stars is uncertain; but if their mass loss is not too extensive, then they are likely to end up as black holes [50, 137]. The remnants may provide the seeds of quasar black holes [215]. Some of the massive stars may end their lives by producing gamma-ray bursts. If so then the broad-band afterglows of these bursts could provide a powerful tool for probing the epoch of reionization [214, 94]). There is no better way to end the dark ages than with -ray burst fireworks.

Where are the first stars or their remnants located today? The very first stars formed in rare high- peaks and hence are likely to populate the cores of present-day galaxies [380]. However, the bulk of the stars which formed in low-mass systems at later times are expected to behave similarly to the collisionless dark matter particles and populate galaxy halos [221].

5.2. The Mass Function of Stars

Currently, we do not have direct observational constraints on how the first stars, the so-called Population III stars, formed at the end of the cosmic dark ages. It is, therefore, instructive to briefly summarize what we have learned about star formation in the present-day Universe, where theoretical reasoning is guided by a wealth of observational data (see [293] for a recent review).

Population I stars form out of cold, dense molecular gas that is structured in a complex, highly inhomogeneous way. The molecular clouds are supported against gravity by turbulent velocity fields and pervaded on large scales by magnetic fields. Stars tend to form in clusters, ranging from a few hundred up to ~ 10⁶ stars. It appears likely that the clustered nature of star formation leads to complicated dynamics and tidal interactions that transport angular momentum, thus allowing the collapsing gas to overcome the classical centrifugal barrier [216]. The initial mass function (IMF) of Pop I stars is observed to have the approximate Salpeter form (e.g., [208])

(105)

where

(106)

The lower cutoff in mass corresponds roughly to the opacity limit for fragmentation. This limit reflects the minimum fragment mass, set when the rate at which gravitational energy is released during the collapse exceeds the rate at which the gas can cool (e.g., [298]). The most important feature of the observed IMF is that ~ 1 M is the characteristic mass scale of Pop I star formation, in the sense that most of the mass goes into stars with masses close to this value. In Figure 25, we show the result from a recent hydrodynamical simulation of the collapse and fragmentation of a molecular cloud core [31, 32]. This simulation illustrates the highly dynamic and chaotic nature of the star formation process ⁶.

Figure 25. A hydrodynamic simulation of the collapse and fragmentation of a turbulent molecular cloud in the present-day Universe (from [32]). The cloud has a mass of 50 M. The panels show the column density through the cloud, and span a scale of 0.4 pc across. Left: The initial phase of the collapse. The turbulence organizes the gas into a network of filaments, and decays thereafter through shocks. Right: A snapshot taken near the end of the simulation, after 1.4 initial free-fall times of 2 × 105yr. Fragmentation has resulted in ~ 50 stars and brown dwarfs. The star formation efficiency is ~ 10% on the scale of the overall cloud, but can be much larger in the dense sub-condensations. This result is in good agreement with what is observed in local star-forming regions.

Bromm & Loeb (2004) [67] carried out idealized simulations of the protostellar accretion problem and estimated the final mass of a Population III star. Using the smoothed particle hydrodynamics (SPH) method, they included the chemistry and cooling physics relevant for the evolution of metal-free gas (see [62] for details). Improving on earlier work [57, 62] by initializing the simulations according to the CDM model, they focused on an isolated overdense region that corresponds to a 3-peak [67]: a halo containing a total mass of 10⁶ M, and collapsing at a redshift z_vir 20. In these runs, one high-density clump has formed at the center of the minihalo, possessing a gas mass of a few hundred solar masses. Soon after its formation, the clump becomes gravitationally unstable and undergoes runaway collapse. Once the gas clump has exceeded a threshold density of 10⁷ cm^-3, it is replaced by a sink particle which is a collisionless point-like particle that is inserted into the simulation. This choice for the density threshold ensures that the local Jeans mass is resolved throughout the simulation. The clump (i.e., sink particle) has an initial mass of M_Cl 200M, and grows subsequently by ongoing accretion of surrounding gas. High-density clumps with such masses result from the chemistry and cooling rate of molecular hydrogen, H₂, which imprint characteristic values of temperature, T ~ 200 K, and density, n ~ 10⁴ cm^-3, into the metal-free gas [62]. Evaluating the Jeans mass for these characteristic values results in M_J ~ a few × 10² M, which is close to the initial clump masses found in the simulations.

Figure 26. Collapse and fragmentation of a primordial cloud (from [67]). Shown is the projected gas density at a redshift z 21.5, briefly after gravitational runaway collapse has commenced in the center of the cloud. Left: The coarse-grained morphology in a box with linear physical size of 23.5 pc. At this time in the unrefined simulation, a high-density clump (sink particle) has formed with an initial mass of ~ 103 M. Right: The refined morphology in a box with linear physical size of 0.5 pc. The central density peak, vigorously gaining mass by accretion, is accompanied by a secondary clump.

Figure 27. Accretion onto a primordial protostar (from [67]). The morphology of this accretion flow is shown in Fig. 26. Left: Accretion rate (in M yr-1) vs. time (in yr) since molecular core formation. Right: Mass of the central core (in M) vs. time. Solid line: Accretion history approximated as: M* t0.45. Using this analytical approximation, we extrapolate that the protostellar mass has grown to ~ 150 M after ~ 105 yr, and to ~ 700 M after ~ 3 × 106 yr, the total lifetime of a very massive star.

The high-density clumps are clearly not stars yet. To probe the subsequent fate of a clump, Bromm & Loeb (2004) [67] have re-simulated the evolution of the central clump with sufficient resolution to follow the collapse to higher densities. Figure 26 (right panel) shows the gas density on a scale of 0.5 pc, which is two orders of magnitude smaller than before. Several features are evident in this plot. First, the central clump does not undergo further sub-fragmentation, and is likely to form a single Population III star. Second, a companion clump is visible at a distance of ~ 0.25 pc. If negative feedback from the first-forming star is ignored, this companion clump would undergo runaway collapse on its own approximately ~ 3 Myr later. This timescale is comparable to the lifetime of a very massive star (VMS) [59]. If the second clump was able to survive the intense radiative heating from its neighbor, it could become a star before the first one explodes as a supernova (SN). Whether more than one star can form in a low-mass halo thus crucially depends on the degree of synchronization of clump formation. Finally, the non-axisymmetric disturbance induced by the companion clump, as well as the angular momentum stored in the orbital motion of the binary system, allow the system to overcome the angular momentum barrier for the collapse of the central clump (see [216]).

The recent discovery of stars like HE0107-5240 with a mass of 0.8 M and an iron abundance of [Fe/H] = -5.3 [90] shows that at least some low mass stars could have formed out of extremely low-metallicity gas. The above simulations show that although the majority of clumps are very massive, a few of them, like the secondary clump in Fig. 26, are significantly less massive. Alternatively, low-mass fragments could form in the dense, shock-compressed shells that surround the first hypernovae [234].

How massive were the first stars? Star formation typically proceeds from the `inside-out', through the accretion of gas onto a central hydrostatic core. Whereas the initial mass of the hydrostatic core is very similar for primordial and present-day star formation [274], the accretion process - ultimately responsible for setting the final stellar mass, is expected to be rather different. On dimensional grounds, the accretion rate is simply related to the sound speed cubed over Newton's constant (or equivalently given by the ratio of the Jeans mass and the free-fall time): M_acc ~ c_s³ / G T^3/2. A simple comparison of the temperatures in present-day star forming regions (T ~ 10 K) with those in primordial ones (T ~ 200-300 K) already indicates a difference in the accretion rate of more than two orders of magnitude.

The above refined simulation enables one to study the three-dimensional accretion flow around the protostar (see also [276, 304, 357]). The gas may now reach densities of 10¹² cm^-3 before being incorporated into a central sink particle. At these high densities, three-body reactions [280] convert the gas into a fully molecular form. Figure 27 shows how the molecular core grows in mass over the first ~ 10⁴ yr after its formation. The accretion rate (left panel) is initially very high, M_acc ~ 0.1 M yr^-1, and subsequently declines according to a power law, with a possible break at ~ 5000 yr. The mass of the molecular core (right panel), taken as an estimator of the proto-stellar mass, grows approximately as: M_* ~ _acc dt t^0.45. A rough upper limit for the final mass of the star is then: M_*(t = 3 × 10⁶ yr) ~ 700 M. In deriving this upper bound, we have conservatively assumed that accretion cannot go on for longer than the total lifetime of a massive star.

Can a Population III star ever reach this asymptotic mass limit? The answer to this question is not yet known with any certainty, and it depends on whether the accretion from a dust-free envelope is eventually terminated by feedback from the star (e.g., [276, 304, 357, 277]). The standard mechanism by which accretion may be terminated in metal-rich gas, namely radiation pressure on dust grains [386], is evidently not effective for gas with a primordial composition. Recently, it has been speculated that accretion could instead be turned off through the formation of an H II region [277], or through the photo-evaporation of the accretion disk [357]. The termination of the accretion process defines the current unsolved frontier in studies of Population III star formation. Current simulations indicate that the first stars were predominantly very massive (> 30 M), and consequently rather different from present-day stellar populations. The crucial question then arises: How and when did the transition take place from the early formation of massive stars to that of low-mass stars at later times? We address this problem next.

The very first stars, marking the cosmic Renaissance of structure formation, formed under conditions that were much simpler than the highly complex environment in present-day molecular clouds. Subsequently, however, the situation rapidly became more complicated again due to the feedback from the first stars on the IGM. Supernova explosions dispersed the nucleosynthetic products from the first generation of stars into the surrounding gas (e.g., [241, 261, 361]), including also dust grains produced in the explosion itself [222, 364]. Atomic and molecular cooling became much more efficient after the addition of these metals. Moreover, the presence of ionizing cosmic rays, as well as of UV and X-ray background photons, modified the thermal and chemical behavior of the gas in important ways (e.g., [232, 233]).

Early metal enrichment was likely the dominant effect that brought about the transition from Population III to Population II star formation. Recent numerical simulations of collapsing primordial objects with overall masses of ~ 10⁶ M, have shown that the gas has to be enriched with heavy elements to a minimum level of Z_crit 10^-3.5 Z, in order to have any effect on the dynamics and fragmentation properties of the system [275, 60, 64]. Normal, low-mass (Population II) stars are hypothesized to only form out of gas with metallicity Z Z_crit. Thus, the characteristic mass scale for star formation is expected to be a function of metallicity, with a discontinuity at Z_crit where the mass scale changes by ~ two orders of magnitude. The redshift where this transition occurs has important implications for the early growth of cosmic structure, and the resulting observational signature (e.g., [392, 141, 234, 322]) include the extended nature of reionization [144].

For additional detailes about the properties of the first stars, see the comprehensive review by Bromm & Larson (2004) [66].

5.3. Gamma-ray Bursts: Probing the First Stars One Star at a Time

Gamma-Ray Bursts (GRBs) are believed to originate in compact remnants (neutron stars or black holes) of massive stars. Their high luminosities make them detectable out to the edge of the visible Universe [94, 214]. GRBs offer the opportunity to detect the most distant (and hence earliest) population of massive stars, the so-called Population III (or Pop III), one star at a time. In the hierarchical assembly process of halos which are dominated by cold dark matter (CDM), the first galaxies should have had lower masses (and lower stellar luminosities) than their low-redshift counterparts. Consequently, the characteristic luminosity of galaxies or quasars is expected to decline with increasing redshift. GRB afterglows, which already produce a peak flux comparable to that of quasars or starburst galaxies at z ~ 1-2, are therefore expected to outshine any competing source at the highest redshifts, when the first dwarf galaxies have formed in the Universe.

Figure 28. Illustration of a long-duration gamma-ray burst in the popular "collapsar" model. The collapse of the core of a massive star (which lost its hydrogen envelope) to a black hole generates two opposite jets moving out at a speed close to the speed of light. The jets drill a hole in the star and shine brightly towards an observer who happened to be located within with the collimation cones of the jets. The jets emenating from a single massive star are so bright that they can be seen across the Universe out to the epoch when the first stars have formed. Upcoming observations by the Swift satellite will have the sensitivity to reveal whether the first stars served as progenitors of gamma-ray bursts (for updates see http://swift.gsfc.nasa.gov/).

The first-year polarization data from the Wilkinson Microwave Anisotropy Probe (WMAP) indicates an optical depth to electron scattering of ~ 17 ± 4% after cosmological recombination [203, 348]. This implies that the first stars must have formed at a redshift z ~ 10 - 20, and reionized a substantial fraction of the intergalactic hydrogen around that time [83, 93, 345, 394, 403]. Early reionization can be achieved with plausible star formation parameters in the standard CDM cosmology; in fact, the required optical depth can be achieved in a variety of very different ionization histories since WMAP places only an integral constraint on these histories [176]. One would like to probe the full history of reionization in order to disentangle the properties and formation history of the stars that are responsible for it. GRB afterglows offer the opportunity to detect stars as well as to probe the metal enrichment level [141] of the intervening IGM.

GRBs, the electromagnetically-brightest explosions in the Universe, should be detectable out to redshifts z > 10 [94, 214]. High-redshift GRBs can be identified through infrared photometry, based on the Ly break induced by absorption of their spectrum at wavelengths below 1.216 µm [(1 + z) / 10]. Follow-up spectroscopy of high-redshift candidates can then be performed on a 10-meter-class telescope. Recently, the ongoing Swift mission [147] has detected a GRB originating at z 6.3 (e.g., [179]), thus demonstrating the viability of GRBs as probes of the early Universe.

There are four main advantages of GRBs relative to traditional cosmic sources such as quasars:

(i) The GRB afterglow flux at a given observed time lag after the -ray trigger is not expected to fade significantly with increasing redshift, since higher redshifts translate to earlier times in the source frame, during which the afterglow is intrinsically brighter [94]. For standard afterglow lightcurves and spectra, the increase in the luminosity distance with redshift is compensated by this cosmological time-stretching effect.

Figure 29. GRB afterglow flux as a function of time since the -ray trigger in the observer frame (taken from [67]). The flux (solid curves) is calculated at the redshifted Ly wavelength. The dotted curves show the planned detection threshold for the James Webb Space Telescope (JWST), assuming a spectral resolution R = 5000 with the near infrared spectrometer, a signal to noise ratio of 5 per spectral resolution element, and an exposure time equal to 20% of the time since the GRB explosion (see http://www.ngst.stsci.edu/nms/main/). Each set of curves shows a sequence of redshifts, namely z = 5, 7, 9, 11, 13, and 15, respectively, from top to bottom.

(ii) As already mentioned, in the standard CDM cosmology, galaxies form hierarchically, starting from small masses and increasing their average mass with cosmic time. Hence, the characteristic mass of quasar black holes and the total stellar mass of a galaxy were smaller at higher redshifts, making these sources intrinsically fainter [391]. However, GRBs are believed to originate from a stellar mass progenitor and so the intrinsic luminosity of their engine should not depend on the mass of their host galaxy. GRB afterglows are therefore expected to outshine their host galaxies by a factor that gets larger with increasing redshift.

(iii) Since the progenitors of GRBs are believed to be stellar, they likely originate in the most common star-forming galaxies at a given redshift rather than in the most massive host galaxies, as is the case for bright quasars [26]. Low-mass host galaxies induce only a weak ionization effect on the surrounding IGM and do not greatly perturb the Hubble flow around them. Hence, the Ly damping wing should be closer to the idealized unperturbed IGM case and its detailed spectral shape should be easier to interpret. Note also that unlike the case of a quasar, a GRB afterglow can itself ionize at most ~ 4 × 10⁴ E₅₁ M of hydrogen if its UV energy is E₅₁ in units of 10⁵¹ ergs (based on the available number of ionizing photons), and so it should have a negligible cosmic effect on the surrounding IGM.

(iv) GRB afterglows have smooth (broken power-law) continuum spectra unlike quasars which show strong spectral features (such as broad emission lines or the so-called "blue bump") that complicate the extraction of IGM absorption features. In particular, the continuum extrapolation into the Ly damping wing (the so-called Gunn-Peterson absorption trough) during the epoch of reionization is much more straightforward for the smooth UV spectra of GRB afterglows than for quasars with an underlying broad Ly emission line [26].

The optical depth of the uniform IGM to Ly absorption is given by (Gunn & Peterson 1965 [163]),

(107)

where H 100h km s^-1 Mpc^-1 _m^1/2(1 + z_s)^3/2 is the Hubble parameter at the source redshift z_s >> 1, f = 0.4162 and = 1216 Å are the oscillator strength and the wavelength of the Ly transition; n_HI(z_s) is the neutral hydrogen density at the source redshift (assuming primordial abundances); _m and _b are the present-day density parameters of all matter and of baryons, respectively; and x_HI is the average fraction of neutral hydrogen. In the second equality we have implicitly considered high-redshifts, (1 + z) >> max[(1 - _m - ) / _m, ( / _m)^1/3], at which the vacuum energy density is negligible relative to matter ( << _m) and the Universe is nearly flat; for _m = 0.3, = 0.7 this corresponds to the condition z >> 1.3 which is well satisfied by the reionization redshift.

At wavelengths longer than Ly at the source, the optical depth obtains a small value; these photons redshift away from the line center along its red wing and never resonate with the line core on their way to the observer. The red damping wing of the Gunn-Peterson trough (Miralda-Escudé 1998 [254])

(108)

where _s is given in equation (107), also we define

(109)

and

(110)

Figure 30. Expected spectral shape of the Ly absorption trough due to intergalactic absorption in GRB afterglows (taken from [67]). The spectrum is presented in terms of the flux density F versus relative observed wavelength , for a source redshift z = 7 (assumed to be prior to the final reionization phase) and the typical halo mass M = 4 × 108 M expected for GRB host galaxies that cool via atomic transitions. Top panel: Two examples for the predicted spectrum including IGM HI absorption (both resonant and damping wing), for host galaxies with (i) an age ts = 107 yr, a UV escape fraction fesc = 10% and a Scalo initial mass function (IMF) in solid curves, or (ii) ts = 108 yr, fesc = 90% and massive (> 100 M) Pop III stars in dashed curves. The observed time after the -ray trigger is one hour, one day, and ten days, from top to bottom, respectively. Bottom panel: Predicted spectra one day after a GRB for a host galaxy with ts = 107 yr, fesc = 10% and a Scalo IMF. Shown is the unabsorbed GRB afterglow (dot-short dashed curve), the afterglow with resonant IGM absorption only (dot-long dashed curve), and the afterglow with full (resonant and damping wing) IGM absorption (solid curve). Also shown, with 1.7 magnitudes of extinction, are the afterglow with full IGM absorption (dotted curve), and attempts to reproduce this profile with a damped Ly absorption system in the host galaxy (dashed curves). (Note, however, that damped absorption of this type could be suppressed by the ionizing effect of the afterglow UV radiation on the surrounding interstellar medium of its host galaxy [289].) Most importantly, the overall spectral shape of the Ly trough carries precious information about the neutral fraction of the IGM at the source redshift; averaging over an ensemble of sources with similar redshifts can reduce ambiguities in the interpretation of each case due to particular local effects.

Although the nature of the central engine that powers the relativistic jets of GRBs is still unknown, recent evidence indicates that long-duration GRBs trace the formation of massive stars (e.g., [365, 383, 45, 211, 47, 264]) and in particular that long-duration GRBs are associated with Type Ib/c supernovae [351]. Since the first stars in the Universe are predicted to be predominantly massive [5, 62, 66], their death might give rise to large numbers of GRBs at high redshifts. In contrast to quasars of comparable brightness, GRB afterglows are short-lived and release ~ 10 orders of magnitude less energy into the surrounding IGM. Beyond the scale of their host galaxy, they have a negligible effect on their cosmological environment ⁷. Consequently, they are ideal probes of the IGM during the reionization epoch. Their rest-frame UV spectra can be used to probe the ionization state of the IGM through the spectral shape of the Gunn-Peterson (Ly) absorption trough, or its metal enrichment history through the intersection of enriched bubbles of supernova (SN) ejecta from early galaxies [141]. Afterglows that are unusually bright (> 10 mJy) at radio frequencies should also show a detectable forest of 21 cm absorption lines due to enhanced HI column densities in sheets, filaments, and collapsed minihalos within the IGM [76, 140].

Another advantage of GRB afterglows is that once they fade away, one may search for their host galaxies. Hence, GRBs may serve as signposts of the earliest dwarf galaxies that are otherwise too faint or rare on their own for a dedicated search to find them. Detection of metal absorption lines from the host galaxy in the afterglow spectrum, offers an unusual opportunity to study the physical conditions (temperature, metallicity, ionization state, and kinematics) in the interstellar medium of these high-redshift galaxies. As Figure 30 indicates, damped Ly absorption within the host galaxy may mask the clear signature of the Gunn-Peterson trough in some galaxies [67]. A small fraction ( ~ 10) of the GRB afterglows are expected to originate at redshifts z > 5 [61, 68]. This subset of afterglows can be selected photometrically using a small telescope, based on the Ly break at a wavelength of 1.216 µm< [(1 + z) / 10], caused by intergalactic HI absorption. The challenge in the upcoming years will be to follow-up on these candidates spectroscopically, using a large (10-meter class) telescope. GRB afterglows are likely to revolutionize observational cosmology and replace traditional sources like quasars, as probes of the IGM at z > 5. The near future promises to be exciting for GRB astronomy as well as for studies of the high-redshift Universe.

It is of great importance to constrain the Pop III star formation mode, and in particular to determine down to which redshift it continues to be prominent. The extent of the Pop III star formation will affect models of the initial stages of reionization (e.g., [394, 93, 343, 403, 12]) and metal enrichment (e.g., [234, 141, 144, 320, 340]), and will determine whether planned surveys will be able to effectively probe Pop III stars (e.g., [319]). The constraints on Pop III star formation will also determine whether the first stars could have contributed a significant fraction to the cosmic near-IR background (e.g., [311, 310, 193, 242, 113]). To constrain high-redshift star formation from GRB observations, one has to address two major questions:

(1) What is the signature of GRBs that originate in metal-free, Pop III progenitors? Simply knowing that a given GRB came from a high redshift is not sufficient to reach a definite conclusion as to the nature of the progenitor. Pregalactic metal enrichment was likely inhomogeneous, and we expect normal Pop I and II stars to exist in galaxies that were already metal-enriched at these high redshifts [68]. Pop III and Pop I/II star formation is thus predicted to have occurred concurrently at z > 5. How is the predicted high mass-scale for Pop III stars reflected in the observational signature of the resulting GRBs? Preliminary results from numerical simulations of Pop III star formation indicate that circumburst densities are systematically higher in Pop III environments. GRB afterglows will then be much brighter than for conventional GRBs. In addition, due to the systematically increased progenitor masses, the Pop III distribution may be biased toward long-duration events.

(2) The modelling of Pop III cosmic star formation histories has a number of free parameters, such as the star formation efficiency and the strength of the chemical feedback. The latter refers to the timescale for, and spatial extent of, the distribution of the first heavy elements that were produced inside of Pop III stars, and subsequently dispersed into the IGM by supernova blast waves. Comparing with theoretical GRB redshift distributions one can use the GRB redshift distribution observed by Swift to calibrate the free model parameters. In particular, one can use this strategy to measure the redshift where Pop III star formation terminates.

Figure 31. Theoretical prediction for the comoving star formation rate (SFR) in units of M yr-1 Mpc-3, as a function of redshift (from [68]). We assume that cooling in primordial gas is due to atomic hydrogen only, a star formation efficiency of * = 10%, and reionization beginning at zreion 17. Solid line: Total comoving SFR. Dotted lines: Contribution to the total SFR from Pop I/II and Pop III for the case of weak chemical feedback. Dashed lines: Contribution to the total SFR from Pop I/II and Pop III for the case of strong chemical feedback. Pop III star formation is restricted to high redshifts, but extends over a significant range, z ~ 10-15.

Figures 31 and 32 illustrate these issues (based on [68]). Figure 32 leads to the robust expectation that ~ 10% of all Swift bursts should originate at z > 5. This prediction is based on the contribution from Population I/II stars which are known to exist even at these high redshifts. Additional GRBs could be triggered by Pop III stars, with a highly uncertain efficiency. Assuming that long-duration GRBs are produced by the collapsar mechanism, a Pop III star with a close binary companion provides a plausible GRB progenitor. The Pop III GRB efficiency, reflecting the probability of forming sufficiently close and massive binary systems, to lie between zero (if tight Pop III binaries do not exist) and ~ 10 times the empirically inferred value for Population I/II (due to the increased fraction of black hole forming progenitors among the massive Pop III stars).

A key ingredient in determining the underlying star formation history from the observed GRB redshift distribution is the GRB luminosity function, which is only poorly constrained at present. The improved statistics provided by Swift will enable the construction of an empirical luminosity function. With an improved luminosity function it would be possible to re-calibrate the theoretical prediction in Figure 2 more reliably.

Figure 32. Predicted GRB rate to be observed by Swift (from [68]). Shown is the observed number of bursts per year, dNGRBobs / dln(1 + z), as a function of redshift. All rates are calculated with a constant GRB efficiency, GRB 2 × 10-9 bursts M-1, using the cosmic SFRs from Fig. 31. Dotted lines: Contribution to the observed GRB rate from Pop I/II and Pop III for the case of weak chemical feedback. Dashed lines: Contribution to the GRB rate from Pop I/II and Pop III for the case of strong chemical feedback. Filled circle: GRB rate from Pop III stars if these were responsible for reionizing the Universe at z ~ 17.

In order to predict the observational signature of high-redshift GRBs, it is important to know the properties of the GRB host systems. Within variants of the popular CDM model for structure formation, where small objects form first and subsequently merge to build up more massive ones, the first stars are predicted to form at z ~ 20 - 30 in minihalos of total mass (dark matter plus gas) ~ 10⁶ M [359, 23, 403]. These objects are the sites for the formation of the first stars, and thus are the potential hosts of the highest-redshift GRBs. What is the environment in which the earliest GRBs and their afterglows did occur? This problem breaks down into two related questions: (i) what type of stars (in terms of mass, metallicity, and clustering properties) will form in each minihalo?, and (ii) how will the ionizing radiation from each star modify the density structure of the surrounding gas? These two questions are fundamentally intertwined. The ionizing photon production strongly depends on the stellar mass, which in turn is determined by how the accretion flow onto the growing protostar proceeds under the influence of this radiation field. In other words, the assembly of the Population III stars and the development of an HII region around them proceed simultaneously, and affect each other. As a preliminary illustration, Figure 27 describes the photo-evaporation as a self-similar champagne flow [337] with parameters appropriate for the Pop III case.

Figure 33. Effect of photoheating from a Population III star on the density profile in a high-redshift minihalo (from [69]). The curves, labeled by the time after the onset of the central point source, are calculated according to a self-similar model for the expansion of an HII region. Numerical simulations closely conform to this analytical behavior. Notice that the central density is significantly reduced by the end of the life of a massive star, and that a central core has developed where the density is constant.

Notice that the central density is significantly reduced by the end of the life of a massive star, and that a central core has developed where the density is nearly constant. Such a flat density profile is markedly different from that created by stellar winds ( r^-2). Winds, and consequently mass-loss, may not be important for massive Population III stars [18, 210], and such a flat density profile may be characteristic of GRBs that originate from metal-free Population III progenitors.

Figure 34. Supernova explosion in the high-redshift Universe (from [65]). The snapshot is taken ~ 106 yr after the explosion with total energy ESN 1053 ergs. We show the projected gas density within a box of linear size 1 kpc. The SN bubble has expanded to a radius of ~ 200 pc, having evacuated most of the gas in the minihalo. Inset: Distribution of metals. The stellar ejecta (gray dots) trace the metals and are embedded in pristine metal-poor gas (black dots).

The first galaxies may be surrounded by a shell of highly enriched material that was carried out in a SN-driven wind (see Fig. 34). A GRB in that galaxy may show strong absorption lines at a velocity separation associated with the wind velocity. Simulating these winds and calculating the absorption profile in the featureless spectrum of a GRB afterglow, will allow us to use the observed spectra of high-z GRBs and directly probe the degree of metal enrichment in the vicinity of the first star forming regions (see [141] for a semi-analytic treatment).

As the early afterglow radiation propagates through the interstellar environment of the GRB, it will likely modify the gas properties close to the source; these changes could in turn be noticed as time-dependent spectral features in the spectrum of the afterglow and used to derive the properties of the gas cloud (density, metal abundance, and size). The UV afterglow radiation can induce detectable changes to the interstellar absorption features of the host galaxy [289]; dust destruction could have occurred due to the GRB X-rays [375, 136], and molecules could have been destroyed near the GRB source [112]. Quantitatively, all of the effects mentioned above strongly depend on the exact properties of the gas in the host system.

Most studies to date have assumed a constant efficiency of forming GRBs per unit mass of stars. This simplifying assumption could lead, under different circumstances, to an overestimation or an underestimation of the frequency of GRBs. Metal-free stars are thought to be massive [5, 62] and their extended envelopes may suppress the emergence of relativistic jets out of their surface (even if such jets are produced through the collapse of their core to a spinning black hole). On the other hand, low-metallicity stars are expected to have weak winds with little angular momentum loss during their evolution, and so they may preferentially yield rotating central configurations that make GRB jets after core collapse.

What kind of metal-free, Pop III progenitor stars may lead to GRBs? Long-duration GRBs appear to be associated with Type Ib/c supernovae [351], namely progenitor massive stars that have lost their hydrogen envelope. This requirement is explained theoretically in the collapsar model, in which the relativistic jets produced by core collapse to a black hole are unable to emerge relativistically out of the stellar surface if the hydrogen envelope is retained [231]. The question then arises as to whether the lack of metal line-opacity that is essential for radiation-driven winds in metal-rich stars, would make a Pop III star retain its hydrogen envelope, thus quenching any relativistic jets and GRBs.

Aside from mass transfer in a binary system, individual PopIII stars could lose their hydrogen envelope due to either: (i) violent pulsations, particularly in the mass range 100 - 140 M, or (ii) a wind driven by helium lines. The outer stellar layers are in a state where gravity only marginally exceeds radiation pressure due to electron-scattering (Thomson) opacity. Adding the small, but still non-negligible contribution from the bound-free opacity provided by singly-ionized helium, may be able to unbind the atmospheric gas. Therefore, mass-loss might occur even in the absence of dust or any heavy elements.

5.4. Emission Spectrum of Metal-Free Stars

The evolution of metal-free (Population III) stars is qualitatively different from that of enriched (Population I and II) stars. In the absence of the catalysts necessary for the operation of the CNO cycle, nuclear burning does not proceed in the standard way. At first, hydrogen burning can only occur via the inefficient PP chain. To provide the necessary luminosity, the star has to reach very high central temperatures (T_c 10^8.1 K). These temperatures are high enough for the spontaneous turn-on of helium burning via the triple- process. After a brief initial period of triple- burning, a trace amount of heavy elements forms. Subsequently, the star follows the CNO cycle. In constructing main-sequence models, it is customary to assume that a trace mass fraction of metals (Z ~ 10^-9) is already present in the star (El Eid 1983 [115]; Castellani et al. 1983 [78]).

Figures 35 and 36 show the luminosity L vs. effective temperature T for zero-age main sequence stars in the mass ranges of 2 - 90 M (Fig. 35) and 100 - 1000 M (Fig. 36). Note that above ~ 100M the effective temperature is roughly constant, T_eff ~ 10⁵ K, implying that the spectrum is independent of the mass distribution of the stars in this regime (Bromm, Kudritzky, & Loeb 2001 [59]). As is evident from these figures (see also Tumlinson & Shull 2000 [370]), both the effective temperature and the ionizing power of metal-free (Pop III) stars are substantially larger than those of metal-rich (Pop I) stars. Metal-free stars with masses > 20M emit between 10⁴⁷ and 10⁴⁸ H I and He I ionizing photons per second per solar mass of stars, where the lower value applies to stars of ~ 20 M and the upper value applies to stars of > 100 M (see Tumlinson & Shull 2000 [370] and Bromm et al. 2001 [59] for more details). Over a lifetime of ~ 3 × 10⁶ years these massive stars produce 10⁴ - 10⁵ ionizing photons per stellar baryon. However, this powerful UV emission is suppressed as soon as the interstellar medium out of which new stars form is enriched by trace amounts of metals. Even though the collapsed fraction of baryons is small at the epoch of reionization, it is likely that most of the stars responsible for the reionization of the Universe formed out of enriched gas.

Figure 35. Luminosity vs. effective temperature for zero-age main sequences stars in the mass range of 2 - 90M (from Tumlinson & Shull 2000 [370]). The curves show Pop I (Z = 0.02) and Pop III stars of mass 2, 5, 8, 10, 15, 20, 25, 30, 35, 40, 50, 60, 70, 80, and 90 M. The diamonds mark decades in metallicity in the approach to Z = 0 from 10-2 down to 10-5 at 2 M, down to 10-10 at 15 M, and down to 10-13 at 90 M. The dashed line along the Pop III ZAMS assumes pure H-He composition, while the solid line (on the left) marks the upper MS with Zc = 10-10 for the M 15 M models. Squares mark the points corresponding to pre-enriched evolutionary models from El Eid et al. (1983) [115] at 80 M and from Castellani et al. (1983) [78] for 25 M.

Figure 36. Same as Figure 35 but for very massive stars above 100M (from Bromm, Kudritzki, & Loeb 2001 [59]). Left solid line: Pop III zero-age main sequence (ZAMS). Right solid line: Pop I ZAMS. In each case, stellar luminosity (in L) is plotted vs. effective temperature (in K). Diamond-shaped symbols: Stellar masses along the sequence, from 100 M (bottom) to 1000 M (top) in increments of 100 M. The Pop III ZAMS is systematically shifted to higher effective temperature, with a value of ~ 105 K which is approximately independent of mass. The luminosities, on the other hand, are almost identical in the two cases.

Will it be possible to infer the initial mass function (IMF) of the first stars from spectroscopic observations of the first galaxies? Figure 37 compares the observed spectrum from a Salpeter IMF (dN_* / dM M^-2.35) and a heavy IMF (with all stars more massive than 100 M) for a galaxy at z_s = 10. The latter case follows from the assumption that each of the dense clumps in the simulations described in the previous section ends up as a single star with no significant fragmentation or mass loss. The difference between the plotted spectra cannot be confused with simple reddening due to normal dust. Another distinguishing feature of the IMF is the expected flux in the hydrogen and helium recombination lines, such as Ly and He II 1640 Å, from the interstellar medium surrounding these stars. We discuss this next.

Figure 37. Comparison of the predicted flux from a Pop III star cluster at zs = 10 for a Salpeter IMF (Tumlinson & Shull 2000 [370]) and a massive IMF (Bromm et al. 2001 [59]). Plotted is the observed flux (in nJy per 106 M of stars) vs. observed wavelength (in µm) for a flat Universe with = 0.7 and h = 0.65. Solid line: The case of a heavy IMF. Dotted line: The fiducial case of a standard Salpeter IMF. The cutoff below obs = 1216 Å (1 + zs) = 1.34 µm is due to complete Gunn-Peterson absorption (which is artificially assumed to be sharp). Clearly, for the same total stellar mass, the observable flux is larger by an order of magnitude for stars which are biased towards having masses > 100 M.

5.5. Emission of Recombination Lines from the First Galaxies

The hard UV emission from a star cluster or a quasar at high redshift is likely reprocessed by the surrounding interstellar medium, producing very strong recombination lines of hydrogen and helium (Oh 1999 [270]; Tumlinson & Shull 2000 [370]; see also Baltz, Gnedin & Silk 1998 [17]). We define _ion to be the production rate of ionizing photons by the source. The emitted luminosity L_line^em per unit stellar mass in a particular recombination line is then estimated to be

(111)

where p_line^em is the probability that a recombination leads to the emission of a photon in the corresponding line, is the frequency of the line and p_cont^esc and p_line^esc are the escape probabilities for the ionizing photons and the line photons, respectively. It is natural to assume that the stellar cluster is surrounded by a finite H II region, and hence that p_cont^esc is close to zero [387, 302]. In addition, p_line^esc is likely close to unity in the H II region, due to the lack of dust in the ambient metal-free gas. Although the emitted line photons may be scattered by neutral gas, they diffuse out to the observer and in the end survive if the gas is dust free. Thus, for simplicity, we adopt a value of unity for p_line^esc.

As a particular example we consider case B recombination which yields p_line^em of about 0.65 and 0.47 for the Ly and He II 1640 Å lines, respectively. These numbers correspond to an electron temperature of ~ 3 × 10⁴K and an electron density of ~ 10² - 10³ cm^-3 inside the H II region [354]. For example, we consider the extreme and most favorable case of metal-free stars all of which are more massive than ~ 100 M. In this case L_line^em = 1.7 × 10³⁷ and 2.2 × 10³⁶ erg s^-1 M^-1 for the recombination luminosities of Ly and He II 1640 Åper stellar mass [59]. A cluster of 10⁶ M in such stars would then produce 4.4 and 0.6 × 10⁹ L in the Ly and He II 1640 Å lines. Comparably-high luminosities would be produced in other recombination lines at longer wavelengths, such as He II 4686 Å and H [270, 271].

The rest - frame equivalent width of the above emission lines measured against the stellar continuum of the embedded star cluster at the line wavelengths is given by

(112)

where L is the spectral luminosity per unit wavelength of the stars at the line resonance. The extreme case of metal-free stars which are more massive than 100 M yields a spectral luminosity per unit frequency L = 2.7 × 10²¹ and 1.8 × 10²¹ erg s^-1 Hz^-1 M^-1 at the corresponding wavelengths [59]. Converting to L, this yields rest-frame equivalent widths of W = 3100 Å and 1100 Å for Ly and He II 1640 Å, respectively. These extreme emission equivalent widths are more than an order of magnitude larger than the expectation for a normal cluster of hot metal-free stars with the same total mass and a Salpeter IMF under the same assumptions concerning the escape probabilities and recombination [209]. The equivalent widths are, of course, larger by a factor of (1 + z_s) in the observer frame. Extremely strong recombination lines, such as Ly and He II 1640 Å, are therefore expected to be an additional spectral signature that is unique to very massive stars in the early Universe. The strong recombination lines from the first luminous objects are potentially detectable with JWST [271].

6. SUPERMASSIVE BLACK HOLES

6.1. The Principle of Self-Regulation

The fossil record in the present-day Universe indicates that every bulged galaxy hosts a supermassive black hole (BH) at its center [206]. This conclusion is derived from a variety of techniques which probe the dynamics of stars and gas in galactic nuclei. The inferred BHs are dormant or faint most of the time, but ocassionally flash in a short burst of radiation that lasts for a small fraction of the Hubble time. The short duty cycle acounts for the fact that bright quasars are much less abundant than their host galaxies, but it begs the more fundamental question: why is the quasar activity so brief? A natural explanation is that quasars are suicidal, namely the energy output from the BHs regulates their own growth.

Supermassive BHs make up a small fraction, < 10^-3, of the total mass in their host galaxies, and so their direct dynamical impact is limited to the central star distribution where their gravitational influence dominates. Dynamical friction on the background stars keeps the BH close to the center. Random fluctuations in the distribution of stars induces a Brownian motion of the BH. This motion can be decribed by the same Langevin equation that captures the motion of a massive dust particle as it responds to random kicks from the much lighter molecules of air around it [86]. The characteristic speed by which the BH wanders around the center is small, ~ (m_* / M_BH)^1/2 _*, where m_* and M_BH are the masses of a single star and the BH, respectively, and _* is the stellar velocity dispersion. Since the random force fluctuates on a dynamical time, the BH wanders across a region that is smaller by a factor of ~ (m_* / M_BH)^1/2 than the region traversed by the stars inducing the fluctuating force on it.

The dynamical insignificance of the BH on the global galactic scale is misleading. The gravitational binding energy per rest-mass energy of galaxies is of order ~ (_* / c)² < 10^-6. Since BH are relativistic objects, the gravitational binding energy of material that feeds them amounts to a substantial fraction its rest mass energy. Even if the BH mass occupies a fraction as small as ~ 10^-4 of the baryonic mass in a galaxy, and only a percent of the accreted rest-mass energy leaks into the gaseous environment of the BH, this slight leakage can unbind the entire gas reservoir of the host galaxy! This order-of-magnitude estimate explains why quasars are short lived. As soon as the central BH accretes large quantities of gas so as to significantly increase its mass, it releases large amounts of energy that would suppress further accretion onto it. In short, the BH growth is self-regulated.

The principle of self-regulation naturally leads to a correlation between the final BH mass, M_bh, and the depth of the gravitational potential well to which the surrounding gas is confined which can be characterized by the velocity dispersion of the associated stars, ~ _*². Indeed such a correlation is observed in the present-day Universe [368]. The observed power-law relation between M_bh and _* can be generalized to a correlation between the BH mass and the circular velocity of the host halo, v_c [130], which in turn can be related to the halo mass, M_halo, and redshift, z [394]

(113)

where _o 10^-5.7 is a constant, and as before [(_m / _m^z)(_c / 18²)], _m^z [1 + ( / _m)(1 + z)^-3]^-1, _c = 18² + 82d - 39d², and d = _m^z - 1. If quasars shine near their Eddington limit as suggested by observations of low and high-redshift quasars [134, 384], then the above value of _o implies that a fraction of ~ 5 - 10% of the energy released by the quasar over a galactic dynamical time needs to be captured in the surrounding galactic gas in order for the BH growth to be self-regulated [394].

Figure 38. Comparison of the observed and model luminosity functions (from [394]). The data points at z < 4 are summarized in Ref. [286], while the light lines show the double power-law fit to the 2dF quasar luminosity function [56]. At z ~ 4.3 and z ~ 6.0 the data is from Refs. [125]. The grey regions show the 1 - range of logarithmic slope ([-2.25,-3.75] at z ~ 4.3 and [-1.6,-3.1] at z ~ 6), and the vertical bars show the uncertainty in the normalization. The open circles show data points converted from the X-ray luminosity function [20] of low luminosity quasars using the median quasar spectral energy distribution. In each panel the vertical dashed lines correspond to the Eddington luminosities of BHs bracketing the observed range of the Mbh - vc relation, and the vertical dotted line corresponds to a BH in a 1013.5 M galaxy.

Various prescriptions for self-regulation were sketched by Silk & Rees [339]. These involve either energy or momentum-driven winds, where the latter type is a factor of ~ v_c / c less efficient [35, 199, 262]. Wyithe & Loeb [394] demonstrated that a particularly simple prescription for an energy-driven wind can reproduce the luminosity function of quasars out to highest measured redshift, z ~ 6 (see Figs. 38 and 40), as well as the observed clustering properties of quasars at z ~ 3 [398] (see Fig. 41). The prescription postulates that: (i) self-regulation leads to the growth of M_bh up the redshift-independent limit as a function of v_c in Eq. (113), for all galaxies throughout their evolution; and (ii) the growth of M_bh to the limiting mass in Eq. (113) occurs through halo merger episodes during which the BH shines at its Eddington luminosity (with the median quasar spectrum) over the dynamical time of its host galaxy, t_dyn. This model has only one adjustable parameter, namely the fraction of the released quasar energy that couples to the surrounding gas in the host galaxy. This parameter can be fixed based on the M_bh - _* relation in the local Universe [130]. It is remarkable that the combination of the above simple prescription and the standard CDM cosmology for the evolution and merger rate of galaxy halos, lead to a satisfactory agreement with the rich data set on quasar evolution over cosmic history.

The cooling time of the heated gas is typically longer than its dynamical time and so the gas should expand into the galactic halo and escape the galaxy if its initial temperature exceeds the virial temperature of the galaxy [394]. The quasar remains active during the dynamical time of the initial gas reservoir, ~ 10⁷ years, and fades afterwards due to the dilution of this reservoir. Accretion is halted as soon as the quasar supplies the galactic gas with more than its binding energy. The BH growth may resume if the cold gas reservoir is replenished through a new merger.

Following the early analytic work, extensive numerical simulations by Springel, Hernquist, & Di Matteo (2005) [350] (see also Di Matteo et al. 2005 [108]) demonstrated that galaxy mergers do produce the observed correlations between black hole mass and spheroid properties when a similar energy feedback is incorporated. Because of the limited resolution near the galaxy nucleus, these simulations adopt a simple prescription for the accretion flow that feeds the black hole. The actual feedback in reality may depend crucially on the geometry of this flow and the physical mechanism that couples the energy or momentum output of the quasar to the surrounding gas.

Figure 39. Simulation images of a merger of galaxies resulting in quasar activity that eventually shuts-off the accretion of gas onto the black hole (from Di Matteo et al. 2005 [108]). The upper (lower) panels show a sequence of snapshots of the gas distribution during a merger with (without) feedback from a central black hole. The temperature of the gas is color coded.

Figure 40. The comoving density of supermassive BHs per unit BH mass (from [394]). The grey region shows the estimate based on the observed velocity distribution function of galaxies in Ref. [336] and the Mbh - vc relation in Eq. (113). The lower bound corresponds to the lower limit in density for the observed velocity function while the grey lines show the extrapolation to lower densities. We also show the mass function computed at z = 1, 3 and 6 from the Press-Schechter [292] halo mass function and Eq. (113), as well as the mass function at z ~ 2.35 and z ~ 3 implied by the observed density of quasars and a quasar lifetime of order the dynamical time of the host galactic disk, tdyn (dot-dashed lines).

Agreement between the predicted and observed correlation function of quasars (Fig. 41) is obtained only if the BH mass scales with redshift as in Eq. (113) and the quasar lifetime is of the order of the dynamical time of the host galactic disk [398],

(114)

This characterizes the timescale it takes low angular momentum gas to settle inwards and feed the black hole from across the galaxy before feedback sets in and suppresses additional infall. It also characterizes the timescale for establishing an outflow at the escape speed from the host spheroid.

The inflow of cold gas towards galaxy centers during the growth phase of the BH would naturally be accompanied by a burst of star formation. The fraction of gas that is not consumed by stars or ejected by supernovae, will continue to feed the BH. It is therefore not surprising that quasar and starburst activities co-exist in Ultra Luminous Infrared Galaxies [356], and that all quasars show broad metal lines indicating a super-solar metallicity of the surrounding gas [106]. Applying a similar self-regulation principle to the stars, leads to the expectation [394, 197] that the ratio between the mass of the BH and the mass in stars is independent of halo mass (as observed locally [243]) but increases with redshift as (z)^1/2(1 + z)^3/2. A consistent trend has indeed been inferred in an observed sample of gravitationally-lensed quasars [305].

Figure 41. Predicted correlation function of quasars at various redshifts in comparison to the 2dF data [101] (from [398]). The dark lines show the correlation function predictions for quasars of various apparent B-band magnitudes. The 2dF limit is B ~ 20.85. The lower right panel shows data from entire 2dF sample in comparison to the theoretical prediction at the mean quasar redshift of <z> = 1.5. The B = 20.85 prediction at this redshift is also shown by thick gray lines in the other panels to guide the eye. The predictions are based on the scaling Mbh vc5 in Eq. (113).

The upper mass of galaxies may also be regulated by the energy output from quasar activity. This would account for the fact that cooling flows are suppressed in present-day X-ray clusters [123, 91, 273], and that massive BHs and stars in galactic bulges were already formed at z ~ 2. The quasars discovered by the Sloan Digital Sky Survey (SDSS) at z ~ 6 mark the early growth of the most massive BHs and galactic spheroids. The present-day abundance of galaxies capable of hosting BHs of mass ~ 10⁹ M (based on Eq. 113) already existed at z ~ 6 [225]. At some epoch, the quasar energy output may have led to the extinction of cold gas in these galaxies and the suppression of further star formation in them, leading to an apparent "anti-hierarchical" mode of galaxy formation where massive spheroids formed early and did not make new stars at late times. In the course of subsequent merger events, the cores of the most massive spheroids acquired an envelope of collisionless matter in the form of already-formed stars or dark matter [225], without the proportional accretion of cold gas into the central BH. The upper limit on the mass of the central BH and the mass of the spheroid is caused by the lack of cold gas and cooling flows in their X-ray halos. In the cores of cooling X-ray clusters, there is often an active central BH that supplies sufficient energy to compensate for the cooling of the gas [91, 123, 35]. The primary physical process by which this energy couples to the gas is still unknown.

Figure 42. The global influence of magnetized quasar outflows on the intergalactic medium (from [139]). Upper Panel: Predicted volume filling fraction of magnetized quasar bubbles F(z), as a function of redshift. Lower Panel: Ratio of normalized magnetic energy density, b / -1, to the fiducial thermal energy density of the intergalactic medium ufid = 3 n(z) k TIGM, where TIGM = 104 K, as a function of redshift (see [139] for more details). In each panel, the solid curves assume that the blast wave created by quasar ouflows is nearly (80%) adiabatic, and that the minimum halo mass of galaxies, Mh,min, is determined by atomic cooling before reionization and by suppression due to galactic infall afterwards (top curve), Mh,min = 109 M (middle curve), and Mh,min = 1010 M (bottom curve). The dashed curve assumes a fully-radiative blast wave and fixes Mh,min by the thresholds for atomic cooling and infall suppression. The vertical dotted line indicates the assumed redshift of complete reionization, zr = 7.

6.2. Feedback on Large Intergalactic Scales

Aside from affecting their host galaxy, quasars disturb their large-scale cosmological environment. Powerful quasar outflows are observed in the form of radio jets [34] or broad-absorption-line winds [160]. The amount of energy carried by these outflows is largely unknown, but could be comparable to the radiative output from the same quasars. Furlanetto & Loeb [139] have calculated the intergalactic volume filled by such outflows as a function of cosmic time (see Fig. 42). This volume is likely to contain magnetic fields and metals, providing a natural source for the observed magnetization of the metal-rich gas in X-ray clusters [207] and in galaxies [103]. The injection of energy by quasar outflows may also explain the deficit of Ly absorption in the vicinity of Lyman-break galaxies [7, 100] and the required pre-heating in X-ray clusters [54, 91].

Beyond the reach of their outflows, the brightest SDSS quasars at z > 6 are inferred to have ionized exceedingly large regions of gas (tens of comoving Mpc) around them prior to global reionization (see Fig. 43 and Refs. [381, 400]). Thus, quasars must have suppressed the faint-end of the galaxy luminosity function in these regions before the same occurred throughout the Universe. The recombination time is comparable to the Hubble time for the mean gas density at z ~ 7 and so ionized regions persist [272] on these large scales where inhomogeneities are small. The minimum galaxy mass is increased by at least an order of magnitude to a virial temperature of ~ 10⁵ K in these ionized regions [23]. It would be particularly interesting to examine whether the faint end (_* < 30 km s^-1) of the luminosity function of dwarf galaxies shows any moduluation on large-scales around rare massive BHs, such as M87.

To find the volume filling fraction of relic regions from z ~ 6, we consider a BH of mass M_bh ~ 3 × 10⁹ M. We can estimate the comoving density of BHs directly from the observed quasar luminosity function and our estimate of quasar lifetime. At z ~ 6, quasars powered by M_bh ~ 3 × 10⁹ M BHs had a comoving density of ~ 0.5G pc^-3 [394]. However, the Hubble time exceeds t_dyn by a factor of ~ 2 × 10² (reflecting the square root of the density contrast of cores of galaxies relative to the mean density of the Universe), so that the comoving density of the bubbles created by the z ~ 6 BHs is ~ 10²Gpc^-3 (see Fig. 40). The density implies that the volume filling fraction of relic z ~ 6 regions is small, < 10%, and that the nearest BH that had M_bh ~ 3 × 10⁹ M at z ~ 6 (and could have been detected as an SDSS quasar then) should be at a distance d_bh ~ (4 / 3 × 10²)^1/3 Gpc ~ 140 Mpc which is almost an order-of-magnitude larger than the distance of M87, a galaxy known to possess a BH of this mass [135].

Figure 43. Quasars serve as probes of the end of reionization. The measured size of the HII regions around SDSS quasars can be used [396, 251] to demonstrate that a significant fraction of the intergalactic hydrogen was neutral at z ~ 6.3 or else the inferred size of the quasar HII regions would have been much larger than observed (assuming typical quasar lifetimes [248]). Also, quasars can be used to measure the redshift at which the intergalactic medium started to transmit Ly photons [381, 400]. The upper panel illustrates how the line-of-sight towards a quasar intersects this transition redshift. The resulting Ly transmission of the intrinsic quasar spectrum is shown schematically in the lower panel.

6.3. What seeded the growth of the supermassive black holes?

The BHs powering the bright SDSS quasars possess a mass of a few × 10⁹ M, and reside in galaxies with a velocity dispersion of ~ 500 km s^-1 [24]. A quasar radiating at its Eddington limiting luminosity, L_e = 1.4 × 10⁴⁶ erg s^-1(M_bh / 10⁸ M), with a radiative efficiency, _rad = L_e / c² would grow exponentially in mass as a function of time t, M_bh = M_seed exp{t / t_e} on a time scale, t_e = 4.1 × 10⁷ yr (_rad / 0.1). Thus, the required growth time in units of the Hubble time t_hubble = 9 × 10⁸ yr[(1 + z) / 7]^-3/2 is

(115)

The age of the Universe at z ~ 6 provides just sufficient time to grow an SDSS BH with M_bh ~ 10⁹ M out of a stellar mass seed with _rad = 10% [175]. The growth time is shorter for smaller radiative efficiencies, as expected if the seed originates from the optically-thick collapse of a supermassive star (in which case M_seed in the logarithmic factor is also larger).

Figure 44. SPH simulation of the collapse of an early dwarf galaxy with a virial temperature just above the cooling threshold of atomic hydrogen and no H2 (from [63]). The image shows a snapshot of the gas density distribution at z 10, indicating the formation of two compact objects near the center of the galaxy with masses of 2.2 × 106 M and 3.1 × 106 M, respectively, and radii < 1 pc. Sub-fragmentation into lower mass clumps is inhibited as long as molecular hydrogen is dissociated by a background UV flux. These circumstances lead to the formation of supermassive stars [220] that inevitably collapse and trigger the birth of supermassive black holes [220, 309]. The box size is 200 pc.

What was the mass of the initial BH seeds? Were they planted in early dwarf galaxies through the collapse of massive, metal free (Pop-III) stars (leading to M_seed of hundreds of solar masses) or through the collapse of even more massive, i.e. supermassive, stars [ 220] ? Bromm & Loeb [63] have shown through a hydrodynamical simulation (see Fig. 44) that supermassive stars were likely to form in early galaxies at z ~ 10 in which the virial temperature was close to the cooling threshold of atomic hydrogen, ~ 10⁴ K. The gas in these galaxies condensed into massive ~ 10⁶ M clumps (the progenitors of supermassive stars), rather than fragmenting into many small clumps (the progenitors of stars), as it does in environments that are much hotter than the cooling threshold. This formation channel requires that a galaxy be close to its cooling threshold and immersed in a UV background that dissociates molecular hydrogen in it. These requirements should make this channel sufficiently rare, so as not to overproduce the cosmic mass density of supermassive BH.

The minimum seed BH mass can be identified observationally through the detection of gravitational waves from BH binaries with Advanced LIGO [395] or with LISA [393]. Most of the mHz binary coalescence events originate at z > 7 if the earliest galaxies included BHs that obey the M_bh - v_c relation in Eq. (113). The number of LISA sources per unit redshift per year should drop substantially after reionization, when the minimum mass of galaxies increased due to photo-ionization heating of the intergalactic medium. Studies of the highest redshift sources among the few hundred detectable events per year, will provide unique information about the physics and history of BH growth in galaxies [327].

The early BH progenitors can also be detected as unresolved point sources, using the future James Webb Space Telescope (JWST). Unfortunately, the spectrum of metal-free massive and supermassive stars is the same, since their surface temperature ~ 10⁵ K is independent of mass [59]. Hence, an unresolved cluster of massive early stars would show the same spectrum as a supermassive star of the same total mass.

It is difficult to ignore the possible environmental impact of quasars on anthropic selection. One may wonder whether it is not a coincidence that our Milky-Way Galaxy has a relatively modest BH mass of only a few million solar masses in that the energy output from a much more massive (e.g. ~ 10⁹ M) black hole would have disrupted the evolution of life on our planet. A proper calculation remains to be done (as in the context of nearby Gamma-Ray Bursts [316]) in order to demonstrate any such link.

7. RADIATIVE FEEDBACK FROM THE FIRST SOURCES OF LIGHT

7.1. Escape of Ionizing Radiation from Galaxies

The intergalactic ionizing radiation field, a key ingredient in the development of reionization, is determined by the amount of ionizing radiation escaping from the host galaxies of stars and quasars. The value of the escape fraction as a function of redshift and galaxy mass remains a major uncertainty in all current studies, and could affect the cumulative radiation intensity by orders of magnitude at any given redshift. Gas within halos is far denser than the typical density of the IGM, and in general each halo is itself embedded within an overdense region, so the transfer of the ionizing radiation must be followed in the densest regions in the Universe. Reionization simulations are limited in resolution and often treat the sources of ionizing radiation and their immediate surroundings as unresolved point sources within the large-scale intergalactic medium (see, e.g., Gnedin 2000 [152]). The escape fraction is highly sensitive to the three-dimensional distribution of the UV sources relative to the geometry of the absorbing gas within the host galaxy (which may allow escape routes for photons along particular directions but not others).

The escape of ionizing radiation (h > 13.6 eV, < 912 Å) from the disks of present-day galaxies has been studied in recent years in the context of explaining the extensive diffuse ionized gas layers observed above the disk in the Milky Way [300] and other galaxies [295, 183]. Theoretical models predict that of order 3 - 14% of the ionizing luminosity from O and B stars escapes the Milky Way disk [111, 110]. A similar escape fraction of f_esc = 6% was determined by Bland-Hawthorn & Maloney (1999) [46] based on H measurements of the Magellanic Stream. From Hopkins Ultraviolet Telescope observations of four nearby starburst galaxies (Leitherer et al. 1995 [217]; Hurwitz, Jelinsky, & Dixon 1997 [185]), the escape fraction was estimated to be in the range 3% < f_esc < 57%. If similar escape fractions characterize high redshift galaxies, then stars could have provided a major fraction of the background radiation that reionized the IGM [236, 238]. However, the escape fraction from high-redshift galaxies, which formed when the Universe was much denser ( (1 + z)³), may be significantly lower than that predicted by models ment to describe present-day galaxies. Current reionization calculations assume that galaxies are isotropic point sources of ionizing radiation and adopt escape fractions in the range 5% < f_esc < 60% [152].

Clumping is known to have a significant effect on the penetration and escape of radiation from an inhomogeneous medium [49, 385, 269, 173, 42]. The inclusion of clumpiness introduces several unknown parameters into the calculation, such as the number and overdensity of the clumps, and the spatial correlation between the clumps and the ionizing sources. An additional complication may arise from hydrodynamic feedback, whereby part of the gas mass is expelled from the disk by stellar winds and supernovae (Section 8).

Wood & Loeb (2000) [387] used a three-dimensional radiation transfer code to calculate the steady-state escape fraction of ionizing photons from disk galaxies as a function of redshift and galaxy mass. The gaseous disks were assumed to be isothermal, with a sound speed c_s ~ 10 km s^-1, and radially exponential, with a scale-length based on the characteristic spin parameter and virial radius of their host halos. The corresponding temperature of ~ 10⁴ K is typical for a gas which is continuousely heated by photo-ionization from stars. The sources of radiation were taken to be either stars embedded in the disk, or a central quasar. For stellar sources, the predicted increase in the disk density with redshift resulted in a strong decline of the escape fraction with increasing redshift. The situation is different for a central quasar. Due to its higher luminosity and central location, the quasar tends to produce an ionization channel in the surrounding disk through which much of its ionizing radiation escapes from the host. In a steady state, only recombinations in this ionization channel must be balanced by ionizations, while for stars there are many ionization channels produced by individual star-forming regions and the total recombination rate in these channels is very high. Escape fractions > 10% were achieved for stars at z ~ 10 only if ~ 90% of the gas was expelled from the disks or if dense clumps removed the gas from the vast majority (> 80%) of the disk volume (see Fig. 45). This analysis applies only to halos with virial temperatures > 10⁴ K. Ricotti & Shull (2000) [302] reached similar conclusions but for a quasi-spherical configuration of stars and gas. They demonstrated that the escape fraction is substantially higher in low-mass halos with a virial temperature < 10⁴ K. However, the formation of stars in such halos depends on their uncertain ability to cool via the efficient production of molecular hydrogen.

Figure 45. Escape fractions of stellar ionizing photons from a gaseous disk embedded within a 1010 M halo which have formed at z = 10 (from Wood & Loeb 2000 [387]). The curves show three different cases of clumpiness within the disk. The volume filling factor refers to either the ionizing emissivity, the gas clumps, or both, depending on the case. The escape fraction is substantial (> 1%) only if the gas distribution is highly clumped.

The main uncertainty in the above predictions involves the distribution of the gas inside the host galaxy, as the gas is exposed to the radiation released by stars and the mechanical energy deposited by supernovae. Given the fundamental role played by the escape fraction, it is desirable to calibrate its value observationally. Steidel, Pettini, & Adelberger [352] reported a detection of significant Lyman continuum flux in the composite spectrum of 29 Lyman break galaxies (LBG) with redshifts in the range z = 3.40 ± 0.09. They co-added the spectra of these galaxies in order to be able to measure the low flux. Another difficulty in the measurement comes from the need to separate the Lyman-limit break caused by the interstellar medium from that already produced in the stellar atmospheres. After correcting for intergalactic absorption, Steidel et al. [352] inferred a ratio between the emergent flux density at 1500 Å and 900 Å (rest frame) of 4.6 ± 1.0. Taking into account the fact that the stellar spectrum should already have an intrinsic Lyman discontinuity of a factor of ~ 3 - 5, but that only ~ 15 - 20% of the 1500 Å photons escape from typical LBGs without being absorbed by dust (Pettini et al. 1998 [290]; Adelberger et al. 2000 [6]), the inferred 900 Å escape fraction is f_esc ~ 10 - 20%. Although the galaxies in this sample were drawn from the bluest quartile of the LBG spectral energy distributions, the measurement implies that this quartile may itself dominate the hydrogen-ionizing background relative to quasars at z ~ 3.

7.2. Propagation of Ionization Fronts in the IGM

The radiation output from the first stars ionizes hydrogen in a growing volume, eventually encompassing almost the entire IGM within a single H II bubble. In the early stages of this process, each galaxy produces a distinct H II region, and only when the overall H II filling factor becomes significant do neighboring bubbles begin to overlap in large numbers, ushering in the "overlap phase" of reionization. Thus, the first goal of a model of reionization is to describe the initial stage, when each source produces an isolated expanding H II region.

We assume a spherical ionized volume V, separated from the surrounding neutral gas by a sharp ionization front. Indeed, in the case of a stellar ionizing spectrum, most ionizing photons are just above the hydrogen ionization threshold of 13.6 eV, where the absorption cross-section is high and a very thin layer of neutral hydrogen is sufficient to absorb all the ionizing photons. On the other hand, an ionizing source such as a quasar produces significant numbers of higher energy photons and results in a thicker transition region.

In the absence of recombinations, each hydrogen atom in the IGM would only have to be ionized once, and the ionized proper volume V_p would simply be determined by

(116)

where _H is the mean number density of hydrogen and N is the total number of ionizing photons produced by the source. However, the increased density of the IGM at high redshift implies that recombinations cannot be neglected. Indeed, in the case of a steady ionizing source (and neglecting the cosmological expansion), a steady-state volume would be reached corresponding to the Strömgren sphere, with recombinations balancing ionizations:

(117)

where the recombination rate depends on the square of the density and on the case B recombination coefficient _b = 2.6 × 10^-13 cm³ s^-1 for hydrogen at T = 10⁴ K. The exact evolution for an expanding H II region, including a non-steady ionizing source, recombinations, and cosmological expansion, is given by (Shapiro & Giroux 1987 [329])

(118)

In this equation, the mean density n_H varies with time as 1 / a³(t). A critical physical ingredient is the dependence of recombination on the square of the density. This means that if the IGM is not uniform, but instead the gas which is being ionized is mostly distributed in high-density clumps, then the recombination time is very short. This is often dealt with by introducing a volume-averaged clumping factor C (in general time-dependent), defined by ⁸

(119)

If the ionized volume is large compared to the typical scale of clumping, so that many clumps are averaged over, then equation (118) can be solved by supplementing it with equation (119) and specifying C. Switching to the comoving volume V, the resulting equation is

(120)

where the present number density of hydrogen is

(121)

This number density is lower than the total number density of baryons _b⁰ by a factor of ~ 0.76, corresponding to the primordial mass fraction of hydrogen. The solution for V(t) (generalized from Shapiro & Giroux 1987 [329]) around a source which turns on at t = t_i is

(122)

where

(123)

At high redshift (when (1 + z) >> |_m^-1 - 1|), the scale factor varies as

(124)

and with the additional assumption of a constant C the function F simplifies as follows. Defining

(125)

we derive

(126)

where the last equality assumes C = 10 and our standard choice of cosmological parameters: _m = 0.3, = 0.7, and _b = 0.045. Although this expression for F(t', t) is in general an accurate approximation at high redshift, in the particular case of the CDM model (where _m + = 1) we get the exact result by replacing equation (125) with

(127)

The size of the resulting H II region depends on the halo which produces it. Consider a halo of total mass M and baryon fraction _b / _m. To derive a rough estimate, we assume that baryons are incorporated into stars with an efficiency of f_star = 10%, and that the escape fraction for the resulting ionizing radiation is also f_esc = 10%. If the stellar IMF is similar to the one measured locally [315], then N 4000 ionizing photons are produced per baryon in stars (for a metallicity equal to 1/20 of the solar value). We define a parameter which gives the overall number of ionizations per baryon,

(128)

If we neglect recombinations then we obtain the maximum comoving radius of the region which the halo of mass M can ionize,

(129)

for our standard set of parameters. The actual radius never reaches this size if the recombination time is shorter than the lifetime of the ionizing source. For an instantaneous starburst with the Scalo (1998) [315] IMF, the production rate of ionizing photons can be approximated as

(130)

where N = 4000, = 4.5, and the most massive stars fade away with the characteristic timescale t_s = 3 × 10⁶ yr. In figure 46 we show the time evolution of the volume ionized by such a source, with the volume shown in units of the maximum volume V_max which corresponds to r_max in equation (129). We consider a source turning on at z = 10 (solid curves) or z = 15 (dashed curves), with three cases for each: no recombinations, C = 1, and C = 10, in order from top to bottom (Note that the result is independent of redshift in the case of no recombinations). When recombinations are included, the volume rises and reaches close to V_max before dropping after the source turns off. At large t recombinations stop due to the dropping density, and the volume approaches a constant value (although V << V_max at large t if C = 10).

Figure 46. Expanding H II region around an isolated ionizing source. The comoving ionized volume V is expressed in units of the maximum possible volume, Vmax = 4 rmax3 / 3 [with rmax given in equation (129)], and the time is measured after an instantaneous starburst which produces ionizing photons according to equation (130). We consider a source turning on at z = 10 (solid curves) or z = 15 (dashed curves), with three cases for each: no recombinations, C = 1, and C = 10, in order from top to bottom. The no-recombination curve is identical for the different source redshifts.

We obtain a similar result for the size of the H II region around a galaxy if we consider a mini-quasar rather than stars. For the typical quasar spectrum (Elvis et al. 1994 [122]), roughly 11,000 ionizing photons are produced per baryon incorporated into the black hole, assuming a radiative efficiency of ~ 6%. The efficiency of incorporating the baryons in a galaxy into a central black hole is low (< 0.6% in the local Universe, e.g. Magorrian et al. 1998 [243]), but the escape fraction for quasars is likely to be close to unity, i.e., an order of magnitude higher than for stars (see previous sub-section). Thus, for every baryon in galaxies, up to ~ 65 ionizing photons may be produced by a central black hole and ~ 40 by stars, although both of these numbers for N_ion are highly uncertain. These numbers suggest that in either case the typical size of H II regions before reionization may be < 1 Mpc or ~ 10 Mpc, depending on whether 10⁸ M halos or 10¹² M halos dominate.

The ionization front around a bright transient source like a quasar expands at early times at nearly the speed of light. This occurs when the HII region is sufficiently small so that the production rate of ionizing photons by the central source exceeds their consumption rate by hydrogen atoms within this volume. It is straightforward to do the accounting for these rates (including recombinations) taking the light propagation delay into account. This was done by Wyithe & Loeb [396] [see also White et al. (2003) [381]] who derived the general equation for the relativistic expansion of the comoving radius [r = (1 + z)r_p] of the quasar H II region in an IGM with a neutral filling fraction x_HI (fixed by other ionizing sources) as,

(131)

where c is the speed of light, C is the clumping factor, _b = 2.6 × 10^-13 cm³ s^-1 is the case-B recombination coefficient at the characteristic temperature of 10⁴ K, and is the rate of ionizing photons crossing a shell at the radius of the HII region at time t. Indeed, for the propagation speed of the proper radius of the HII region r_p = r / (1 + z) approaches the speed of light in the above expression, (dr_p / dt) c. The actual size of the HII region along the line-of-sight to a quasar can be inferred from the extent of the spectral gap between the quasar's rest-frame Ly wavelength and the start of Ly absorption by the IGM in the observed spectrum. Existing data from the SDSS quasars [396, 251, 401] provide typical values of r_p ~ 5 Mpc and indicate for plausible choices of the quasar lifetimes that x_HI > 0.1 at z > 6. These ionized bubbles could be imaged directly by future 21cm maps of the regions around the highest-redshift quasars [367, 397, 390].

The profile of the Ly emission line of galaxies has also been suggested as a probe of the ionization state of the IGM [223, 314, 81, 177, 240, 227, 246]. If the IGM is neutral, then the damping wing of the Gunn-Peterson trough in equation (108) is modified since Ly absorption starts only from the near edge of the ionized region along the line-of-sight to the source [81, 240]. Rhoads & Malhotra [246] showed that the observed abundance of galaxies with Ly emission at z ~ 6.5 indicates that a substantial fraction (tens of percent) of the IGM must be ionized in order to allow transmission of the observed Ly photons. However, if these galaxies reside in groups, then galaxies with peculiar velocities away from the observer will preferentially Doppler-shift the emitted Ly photons to the red wing of the Ly resonance and reduce the depression of the line profile [227, 85]. Additional uncertainties in the intrinsic line profile based on the geometry and the stellar or gaseous contents of the source galaxy [227, 314], as well as the clustering of galaxies which ionize their immediate environment in groups [400, 145], limits this method from reaching robust conclusions. Imaging of the expected halos of scattered Ly radiation around galaxies embedded in a neutral IGM [223, 307] provide a more definitive test of the neutrality of the IGM, but is more challenging observationally.

7.3. Reionization of Hydrogen

In this section we summarize recent progress, both analytic and numerical, made toward elucidating the basic physics of reionization and the way in which the characteristics of reionization depend on the nature of the ionizing sources and on other input parameters of cosmological models.

The process of the reionization of hydrogen involves several distinct stages. The initial, "pre-overlap" stage (using the terminology of Gnedin [152]) consists of individual ionizing sources turning on and ionizing their surroundings. The first galaxies form in the most massive halos at high redshift, and these halos are biased and are preferentially located in the highest-density regions. Thus the ionizing photons which escape from the galaxy itself (see Section 7.1) must then make their way through the surrounding high-density regions, which are characterized by a high recombination rate. Once they emerge, the ionization fronts propagate more easily into the low-density voids, leaving behind pockets of neutral, high-density gas. During this period the IGM is a two-phase medium characterized by highly ionized regions separated from neutral regions by ionization fronts. Furthermore, the ionizing intensity is very inhomogeneous even within the ionized regions, with the intensity determined by the distance from the nearest source and by the ionizing luminosity of this source.

The central, relatively rapid "overlap" phase of reionization begins when neighboring H II regions begin to overlap. Whenever two ionized bubbles are joined, each point inside their common boundary becomes exposed to ionizing photons from both sources. Therefore, the ionizing intensity inside H II regions rises rapidly, allowing those regions to expand into high-density gas which had previously recombined fast enough to remain neutral when the ionizing intensity had been low. Since each bubble coalescence accelerates the process of reionization, the overlap phase has the character of a phase transition and is expected to occur rapidly, over less than a Hubble time at the overlap redshift. By the end of this stage most regions in the IGM are able to see several unobscured sources, and therefore the ionizing intensity is much higher than before overlap and it is also much more homogeneous. An additional ingredient in the rapid overlap phase results from the fact that hierarchical structure formation models predict a galaxy formation rate that rises rapidly with time at the relevant redshift range. This process leads to a state in which the low-density IGM has been highly ionized and ionizing radiation reaches everywhere except for gas located inside self-shielded, high-density clouds. This marks the end of the overlap phase, and this important landmark is most often referred to as the 'moment of reionization'.

Some neutral gas does, however, remain in high-density structures which correspond to Lyman Limit systems and damped Ly systems seen in absorption at lower redshifts. The high-density regions are gradually ionized as galaxy formation proceeds, and the mean ionizing intensity also grows with time. The ionizing intensity continues to grow and to become more uniform as an increasing number of ionizing sources is visible to every point in the IGM. This "post-overlap" phase continues indefinitely, since collapsed objects retain neutral gas even in the present Universe. The IGM does, however, reach another milestone at z ~ 1.6, the breakthrough redshift [239]. Below this redshift, all ionizing sources are visible to each other, while above this redshift absorption by the Ly forest implies that only sources in a small redshift range are visible to a typical point in the IGM.

Semi-analytic models of the pre-overlap stage focus on the evolution of the H II filling factor, i.e., the fraction of the volume of the Universe which is filled by H II regions. We distinguish between the naive filling factor F_{H II} and the actual filling factor or porosity Q_{H II}. The naive filling factor equals the number density of bubbles times the average volume of each, and it may exceed unity since when bubbles begin to overlap the overlapping volume is counted multiple times. However, as explained below, in the case of reionization the linearity of the physics means that F_{H II} is a very good approximation to Q_{H II} up to the end of the overlap phase of reionization.

The model of individual H II regions presented in the previous section can be used to understand the development of the total filling factor. Starting with equation (120), if we assume a common clumping factor C for all H II regions then we can sum each term of the equation over all bubbles in a given large volume of the Universe, and then divide by this volume. Then V is replaced by the filling factor and N by the total number of ionizing photons produced up to some time t, per unit volume. The latter quantity equals the mean number of ionizing photons per baryon times the mean density of baryons n_b. Following the arguments leading to equation (129), we find that if we include only stars then

(132)

where the collapse fraction F_col is the fraction of all the baryons in the Universe which are in galaxies, i.e., the fraction of gas which settles into halos and cools efficiently inside them. In writing equation (132) we are assuming instantaneous production of photons, i.e., that the timescale for the formation and evolution of the massive stars in a galaxy is short compared to the Hubble time at the formation redshift of the galaxy. In a model based on equation (120), the near-equality between F_{H II} and Q_{H II} results from the linearity of this equation. First, the total number of ionizations equals the total number of ionizing photons produced by stars, i.e., all ionizing photons contribute regardless of the spatial distribution of sources; and second, the total recombination rate is proportional to the total ionized volume, regardless of its topology. Thus, even if two or more bubbles overlap the model remains an accurate approximation for Q_{H II} (at least until Q_{H II} becomes nearly equal to 1). Note, however, that there still are a number of important simplifications in the model, including the assumption of a homogeneous (though possibly time-dependent) clumping factor, and the neglect of feedback whereby the formation of one galaxy may suppress further galaxy formation in neighboring regions. These complications are discussed in detail below and in Section 7.5 and Section 8.

Under these assumptions we convert equation (120), which describes individual H II regions, to an equation which statistically describes the transition from a neutral Universe to a fully ionized one (compare to Madau et al. 1999 [239] and Haiman & Loeb 1997 [171]):

(133)

where we assumed a primordial mass fraction of hydrogen of 0.76. The solution (in analogy with equation (122)) is

(134)

where F(t', t) is determined by equations (123) - (127).

A simple estimate of the collapse fraction at high redshift is the mass fraction (given by equation (91) in the Press-Schechter model) in halos above the cooling threshold, which is the minimum mass of halos in which gas can cool efficiently. Assuming that only atomic cooling is effective during the redshift range of reionization, the minimum mass corresponds roughly to a halo of virial temperature T_vir = 10⁴ K, which can be converted to a mass using equation (86). With this prescription we derive (for N_ion = 40) the reionization history shown in Fig. 47 for the case of a constant clumping factor C. The solid curves show Q_{H II} as a function of redshift for a clumping factor C = 0 (no recombinations), C = 1, C = 10, and C = 30, in order from left to right. Note that if C ~ 1 then recombinations are unimportant, but if C > 10 then recombinations significantly delay the reionization redshift (for a fixed star-formation history). The dashed curve shows the collapse fraction F_col in this model. For comparison, the vertical dotted line shows the z = 5.8 observational lower limit (Fan et al. 2000 [124]) on the reionization redshift.

Figure 47. Semi-analytic calculation of the reionization of the IGM (for Nion = 40), showing the redshift evolution of the filling factor QH II. Solid curves show QH II for a clumping factor C = 0 (no recombinations), C = 1, C = 10, and C = 30, in order from left to right. The dashed curve shows the collapse fraction Fcol, and the vertical dotted line shows the z = 5.8 observational lower limit (Fan et al. 2000 [124]) on the reionization redshift.

Clearly, star-forming galaxies in CDM hierarchical models are capable of ionizing the Universe at z ~ 6 - 15 with reasonable parameter choices. This has been shown by a large number of theoretical, semi-analytic calculations [138, 330, 171, 373, 89, 92, 392, 83, 371] as well as numerical simulations [79, 148, 152, 2, 296, 95, 342, 204, 186]. Similarly, if a small fraction (< 1%) of the gas in each galaxy accretes onto a central black hole, then the resulting mini-quasars are also able to reionize the Universe, as has also been shown using semi-analytic models [138, 172, 373, 392].

Although many models yield a reionization redshift around 7 - 12, the exact value depends on a number of uncertain parameters affecting both the source term and the recombination term in equation (133). The source parameters include the formation efficiency of stars and quasars and the escape fraction of ionizing photons produced by these sources. The formation efficiency of low mass galaxies may also be reduced by feedback from galactic outflows. These parameters affecting the sources are discussed elsewhere in this review (see Section 7.1, and 8). Even when the clumping is inhomogeneous, the recombination term in equation (133) is generally valid if C is defined as in equation (119), where we take a global volume average of the square of the density inside ionized regions (since neutral regions do not contribute to the recombination rate). The resulting mean clumping factor depends on the density and clustering of sources, and on the distribution and topology of density fluctuations in the IGM. Furthermore, the source halos should tend to form in overdense regions, and the clumping factor is affected by this cross-correlation between the sources and the IGM density.

Miralda-Escudé, Haehnelt, & Rees (2000) [256] presented a simple model for the distribution of density fluctuations, and more generally they discussed the implications of inhomogeneous clumping during reionization. They noted that as ionized regions grow, they more easily extend into low-density regions, and they tend to leave behind high-density concentrations, with these neutral islands being ionized only at a later stage. They therefore argued that, since at high-redshift the collapse fraction is low, most of the high-density regions, which would dominate the clumping factor if they were ionized, will in fact remain neutral and occupy only a tiny fraction of the total volume. Thus, the development of reionization through the end of the overlap phase should occur almost exclusively in the low-density IGM, and the effective clumping factor during this time should be ~ 1, making recombinations relatively unimportant (see Fig. 47). Only in the post-reionization phase, Miralda-Escudé et al. (2000) [256] argued, do the high density clouds and filaments become gradually ionized as the mean ionizing intensity further increases.

The complexity of the process of reionization is illustrated by the numerical simulation of Gnedin [152] of stellar reionization (in CDM with _m = 0.3). This simulation uses a formulation of radiative transfer which relies on several rough approximations; although it does not include the effect of shadowing behind optically-thick clumps, it does include for each point in the IGM the effects of an estimated local optical depth around that point, plus a local optical depth around each ionizing source. This simulation helps to understand the advantages of the various theoretical approaches, while pointing to the complications which are not included in the simple models. Figures 48 and 49, taken from Figure 3 in [152], show the state of the simulated Universe just before and just after the overlap phase, respectively. They show a thin (15 h^-1 comoving kpc) slice through the box, which is 4 h^-1 Mpc on a side, achieves a spatial resolution of 1 h^-1 kpc, and uses 128³ each of dark matter particles and baryonic particles (with each baryonic particle having a mass of 5 × 10⁵ M). The figures show the redshift evolution of the mean ionizing intensity J₂₁ (upper right panel), and visually the logarithm of the neutral hydrogen fraction (upper left panel), the gas density (lower left panel), and the gas temperature (lower right panel). Note the obvious features resulting from the periodic boundary conditions assumed in the simulation. Also note that the intensity J₂₁ is defined as the intensity at the Lyman limit, expressed in units of 10^-21 erg cm^-2 s^-1 sr ^-1Hz^-1. For a given source emission, the intensity inside H II regions depends on absorption and radiative transfer through the IGM (e.g., Haardt & Madau 1996 [166]; Abel & Haehnelt 1999 [1])

Figure 48. Visualization at z = 7.7 of a numerical simulation of reionization, adopted from Figure 3c of [152]. The panels display the logarithm of the neutral hydrogen fraction (upper left), the gas density (lower left), and the gas temperature (lower right). Also shown is the redshift evolution of the logarithm of the mean ionizing intensity (upper right). Note the periodic boundary conditions.

Figure 49. Visualization at z = 6.7 of a numerical simulation of reionization, adopted from Figure 3e of [152]. The panels display the logarithm of the neutral hydrogen fraction (upper left), the gas density (lower left), and the gas temperature (lower right). Also shown is the redshift evolution of the logarithm of the mean ionizing intensity (upper right). Note the periodic boundary conditions.

Figure 48 shows the two-phase IGM at z = 7.7, with ionized bubbles emanating from one main concentration of sources (located at the right edge of the image, vertically near the center; note the periodic boundary conditions). The bubbles are shown expanding into low density regions and beginning to overlap at the center of the image. The topology of ionized regions is clearly complex: While the ionized regions are analogous to islands in an ocean of neutral hydrogen, the islands themselves contain small lakes of dense neutral gas. One aspect which has not been included in theoretical models of clumping is clear from the figure. The sources themselves are located in the highest density regions (these being the sites where the earliest galaxies form) and must therefore ionize the gas in their immediate vicinity before the radiation can escape into the low density IGM. For this reason, the effective clumping factor is of order 100 in the simulation and also, by the overlap redshift, roughly ten ionizing photons have been produced per baryon. Figure 49 shows that by z = 6.7 the low density regions have all become highly ionized along with a rapid increase in the ionizing intensity. The only neutral islands left are the highest density regions (compare the two panels on the left). However, we emphasize that the quantitative results of this simulation must be considered preliminary, since the effects of increased resolution and a more accurate treatment of radiative transfer are yet to be explored. Methods are being developed for incorporating a more complete treatment of radiative transfer into three dimensional cosmological simulations (e.g., [2, 296, 95, 342, 204, 186]).

Gnedin, Ferrara, & Zweibel (2000) [151] investigated an additional effect of reionization. They showed that the Biermann battery in cosmological ionization fronts inevitably generates coherent magnetic fields of an amplitude ~ 10^-19 Gauss. These fields form as a result of the breakout of the ionization fronts from galaxies and their propagation through the H I filaments in the IGM. Although the fields are too small to directly affect galaxy formation, they could be the seeds for the magnetic fields observed in galaxies and X-ray clusters today.

If quasars contribute substantially to the ionizing intensity during reionization then several aspects of reionization are modified compared to the case of pure stellar reionization. First, the ionizing radiation emanates from a single, bright point-source inside each host galaxy, and can establish an escape route (H II funnel) more easily than in the case of stars which are smoothly distributed throughout the galaxy (Section 7.1). Second, the hard photons produced by a quasar penetrate deeper into the surrounding neutral gas, yielding a thicker ionization front. Finally, the quasar X-rays catalyze the formation of H₂ molecules and allow stars to keep forming in very small halos.

Oh (1999) [270] showed that star-forming regions may also produce significant X-rays at high redshift. The emission is due to inverse Compton scattering of CMB photons off relativistic electrons in the ejecta, as well as thermal emission by the hot supernova remnant. The spectrum expected from this process is even harder than for typical quasars, and the hard photons photoionize the IGM efficiently by repeated secondary ionizations. The radiation, characterized by roughly equal energy per logarithmic frequency interval, would produce a uniform ionizing intensity and lead to gradual ionization and heating of the entire IGM. Thus, if this source of emission is indeed effective at high redshift, it may have a crucial impact in changing the topology of reionization. Even if stars dominate the emission, the hardness of the ionizing spectrum depends on the initial mass function. At high redshift it may be biased toward massive, efficiently ionizing stars, but this remains very much uncertain.

Semi-analytic as well as numerical models of reionization depend on an extrapolation of hierarchical models to higher redshifts and lower-mass halos than the regime where the models have been compared to observations (see e.g. [392, 83, 371]). These models have the advantage that they are based on the current CDM paradigm which is supported by a variety of observations of large-scale structure, galaxy clustering, and the CMB. The disadvantage is that the properties of high-redshift galaxies are derived from those of their host halos by prescriptions which are based on low redshift observations, and these prescriptions will only be tested once abundant data is available on galaxies which formed during the reionization era (see [392] for the sensitivity of the results to model parameters). An alternative approach to analyzing the possible ionizing sources which brought about reionization is to extrapolate from the observed populations of galaxies and quasars at currently accessible redshifts. This has been attempted, e.g., by Madau et al. (1999) [239] and Miralda-Escudé et al. (2000) [256]. The general conclusion is that a high-redshift source population similar to the one observed at z = 3 - 4 would produce roughly the needed ionizing intensity for reionization. However, Dijkstra, Haiman, & Loeb (2004) [107] constrained the role of quasars in reionizing the Universe based on the unresolved flux of the X-ray background. At any event, a precise conclusion remains elusive because of the same kinds of uncertainties as those found in the models based on CDM: The typical escape fraction, and the faint end of the luminosity function, are both not well determined even at z = 3 - 4, and in addition the clumping factor at high redshift must be known in order to determine the importance of recombinations. Future direct observations of the source population at redshifts approaching reionization may help resolve some of these questions.

7.4. Photo-evaporation of Gaseous Halos After Reionization

The end of the reionization phase transition resulted in the emergence of an intense UV background that filled the Universe and heated the IGM to temperatures of ~ 1 - 2 × 10⁴ K (see the previous section). After ionizing the rarefied IGM in the voids and filaments on large scales, the cosmic UV background penetrated the denser regions associated with the virialized gaseous halos of the first generation of objects. A major fraction of the collapsed gas had been incorporated by that time into halos with a virial temperature < 10⁴ K, where the lack of atomic cooling prevented the formation of galactic disks and stars or quasars. Photoionization heating by the cosmic UV background could then evaporate much of this gas back into the IGM. The photo-evaporating halos, as well as those halos which did retain their gas, may have had a number of important consequences just after reionization as well as at lower redshifts.

In this section we focus on the process by which gas that had already settled into virialized halos by the time of reionization was evaporated back into the IGM due to the cosmic UV background. This process was investigated by Barkana & Loeb (1999) [22] using semi-analytic methods and idealized numerical calculations. They first considered an isolated spherical, centrally-concentrated dark matter halo containing gas. Since most of the photo-evaporation occurs at the end of overlap, when the ionizing intensity builds up almost instantaneously, a sudden illumination by an external ionizing background may be assumed. Self-shielding of the gas implies that the halo interior sees a reduced intensity and a harder spectrum, since the outer gas layers preferentially block photons with energies just above the Lyman limit. It is useful to parameterize the external radiation field by a specific intensity per unit frequency, ,

(135)

where _L is the Lyman limit frequency, and J₂₁ is the intensity at _L expressed in units of 10^-21 erg cm^-2 s^-1 sr ^-1Hz^-1. The intensity is normalized to an expected post - reionization value of around unity for the ratio of ionizing photon density to the baryon density. Different power laws can be used to represent either quasar spectra ( ~ 1.8) or stellar spectra ( ~ 5).

Once the gas is heated throughout the halo, some fraction of it acquires a sufficiently high temperature that it becomes unbound. This gas expands due to the resulting pressure gradient and eventually evaporates back into the IGM. The pressure gradient force (per unit volume) k (T / µm_p) competes with the gravitational force of G M / r². Due to the density gradient, the ratio between the pressure force and the gravitational force is roughly equal to the ratio between the thermal energy ~ k T and the gravitational binding energy ~ µm_p G M / r (which is ~ k T_vir at the virial radius r_vir) per particle. Thus, if the kinetic energy exceeds the potential energy (or roughly if T> T_vir), the repulsive pressure gradient force exceeds the attractive gravitational force and expels the gas on a dynamical time (or faster for halos with T >> T_vir).

The left panel of Figure 50 (adopted from Fig. 3 of Barkana & Loeb 1999 [22]) shows the fraction of gas within the virial radius which becomes unbound after reionization, as a function of the total halo circular velocity, with halo masses at z = 8 indicated at the top. The two pairs of curves correspond to spectral index = 5 (solid) or = 1.8 (dashed). In each pair, a calculation which assumes an optically-thin halo leads to the upper curve, but including radiative transfer and self-shielding modifies the result to the one shown by the lower curve. In each case self-shielding lowers the unbound fraction, but it mostly affects only a neutral core containing ~ 30% of the gas. Since high energy photons above the Lyman limit penetrate deep into the halo and heat the gas efficiently, a flattening of the spectral slope from = 5 to = 1.8 raises the unbound gas fraction. This figure is essentially independent of redshift if plotted in terms of circular velocity, but the conversion to a corresponding mass does vary with redshift. The characteristic circular velocity where most of the gas is lost is ~ 10 - 15 km s^-1, but clearly the effect of photo-evaporation is gradual, going from total gas removal down to no effect over a range of a factor of ~ 100 in halo mass.

Figure 50. Effect of photo-evaporation on individual halos and on the overall halo population. The left panel shows the unbound gas fraction (within the virial radius) versus total halo velocity dispersion or mass, adopted from Figure 3 of Barkana & Loeb (1999) [22]. The two pairs of curves correspond to spectral index = 5 (solid) or = 1.8 (dashed), in each case at z = 8. In each pair, assuming an optically-thin halo leads to the upper curve, while the lower curve shows the result of including radiative transfer and self shielding. The right panel shows the total fraction of gas in the Universe which evaporates from halos at reionization, versus the reionization redshift, adopted from Figure 7 of Barkana & Loeb (1999) [22]. The solid line assumes a spectral index = 1.8, and the dotted line assumes = 5.

Given the values of the unbound gas fraction in halos of different masses, the Press-Schechter mass function (Section 4.1) can be used to calculate the total fraction of the IGM which goes through the process of accreting onto a halo and then being recycled into the IGM at reionization. The low-mass cutoff in this sum over halos is given by the lowest mass halo in which gas has assembled by the reionization redshift. This mass can be estimated by the linear Jeans mass M_J in equation (62). The Jeans mass does not in general precisely equal the limiting mass for accretion (see the discussion in the next section). Indeed, at a given redshift some gas can continue to fall into halos of lower mass than the Jeans mass at that redshift. On the other hand, the larger Jeans mass at higher redshifts means that a time-averaged Jeans mass may be more appropriate, as indicated by the filtering mass. In practice, the Jeans mass is sufficiently accurate since at z ~ 10 - 20 it agrees well with the values found in the numerical spherical collapse calculations of Haiman, Thoul, & Loeb (1996) [168].

The right panel of Figure 50 (adopted from Fig. 7 of Barkana & Loeb 1999 [22]) shows the total fraction of gas in the Universe which evaporates from halos at reionization, versus the reionization redshift. The solid line assumes a spectral index = 1.8, and the dotted line assumes = 5, showing that the result is insensitive to the spectrum. Even at high redshift, the amount of gas which participates in photo-evaporation is significant, which suggests a number of possible implications as discussed below. The gas fraction shown in the figure represents most (~ 60 - 80% depending on the redshift) of the collapsed fraction before reionization, although some gas does remain in more massive halos.

The photo-evaporation of gas out of large numbers of halos may have interesting implications. First, gas which falls into halos and is expelled at reionization attains a different entropy than if it had stayed in the low-density IGM. The resulting overall reduction in the entropy is expected to be small - the same as would be produced by reducing the temperature of the entire IGM by a factor of ~ 1.5 - but localized effects near photo-evaporating halos may be more significant. Furthermore, the resulting ~ 20 km s^-1 outflows induce small-scale fluctuations in peculiar velocity and temperature. These outflows are usually well below the resolution limit of most numerical simulations, but some outflows were resolved in the simulation of Bryan et al. (1998) [70]. The evaporating halos may consume a significant number of ionizing photons in the post-overlap stage of reionization [174, 186], but a definitive determination requires detailed simulations which include the three-dimensional geometry of source halos and sink halos.

Although gas is quickly expelled out of the smallest halos, photo-evaporation occurs more gradually in larger halos which retain some of their gas. These surviving halos initially expand but they continue to accrete dark matter and to merge with other halos. These evaporating gas halos could contribute to the high column density end of the Ly forest [51]. Abel & Mo (1998) [3] suggested that, based on the expected number of surviving halos, a large fraction of the Lyman limit systems at z ~ 3 may correspond to mini-halos that survived reionization. Surviving halos may even have identifiable remnants in the present Universe. These ideas thus offer the possibility that a population of halos which originally formed prior to reionization may correspond almost directly to several populations that are observed much later in the history of the Universe. However, the detailed dynamics of photo-evaporating halos are complex, and detailed simulations are required to confirm these ideas. Photo-evaporation of a gas cloud has been followed in a two dimensional simulation with radiative transfer, by Shapiro & Raga (2000) [331]. They found that an evaporating halo would indeed appear in absorption as a damped Ly system initially, and as a weaker absorption system subsequently. Future simulations [186] will clarify the contribution to quasar absorption lines of the entire population of photo-evaporating halos.

7.5. Suppression of the Formation of Low Mass Galaxies

At the end of overlap, the cosmic ionizing background increased sharply, and the IGM was heated by the ionizing radiation to a temperature > 10⁴ K. Due to the substantial increase in the IGM temperature, the intergalactic Jeans mass increased dramatically, changing the minimum mass of forming galaxies [299, 117, 148, 255].

Gas infall depends sensitively on the Jeans mass. When a halo more massive than the Jeans mass begins to form, the gravity of its dark matter overcomes the gas pressure. Even in halos below the Jeans mass, although the gas is initially held up by pressure, once the dark matter collapses its increased gravity pulls in some gas [168]. Thus, the Jeans mass is generally higher than the actual limiting mass for accretion. Before reionization, the IGM is cold and neutral, and the Jeans mass plays a secondary role in limiting galaxy formation compared to cooling. After reionization, the Jeans mass is increased by several orders of magnitude due to the photoionization heating of the IGM, and hence begins to play a dominant role in limiting the formation of stars. Gas infall in a reionized and heated Universe has been investigated in a number of numerical simulations. Thoul & Weinberg (1996) [363] inferred, based on a spherically-symmetric collapse simulation, a reduction of ~ 50% in the collapsed gas mass due to heating, for a halo of circular velocity V_c ~ 50 km s^-1 at z = 2, and a complete suppression of infall below V_c ~ 30 km s^-1. Kitayama & Ikeuchi (2000) [201] also performed spherically-symmetric simulations but included self-shielding of the gas, and found that it lowers the circular velocity thresholds by ~ 5 km s^-1. Three dimensional numerical simulations [294, 378, 267] found a significant suppression of gas infall in even larger halos (V_c ~ 75 km s^-1), but this was mostly due to a suppression of late infall at z < 2.

When a volume of the IGM is ionized by stars, the gas is heated to a temperature T_IGM ~ 10⁴ K. If quasars dominate the UV background at reionization, their harder photon spectrum leads to T_IGM > 2 × 10⁴ K. Including the effects of dark matter, a given temperature results in a linear Jeans mass corresponding to a halo circular velocity of

(136)

where we used equation (85) and assumed µ = 0.6. In halos with V_c > V_J, the gas fraction in infalling gas equals the universal mean of _b / _m, but gas infall is suppressed in smaller halos. Even for a small dark matter halo, once it collapses to a virial overdensity of _c / _m ^z relative to the mean, it can pull in additional gas. A simple estimate of the limiting circular velocity, below which halos have essentially no gas infall, is obtained by substituting the virial overdensity for the mean density in the definition of the Jeans mass. The resulting estimate is

(137)

This value is in rough agreement with the numerical simulations mentioned before. A more recent study by Dijkstra et al. (2004) [107] indicates that at the high redshifts of z > 10 gas could nevertheless assemble into halos with circular velocities as low as v_c ~ 10 km s^-1, even in the presence of a UV background.

Although the Jeans mass is closely related to the rate of gas infall at a given time, it does not directly yield the total gas residing in halos at a given time. The latter quantity depends on the entire history of gas accretion onto halos, as well as on the merger histories of halos, and an accurate description must involve a time-averaged Jeans mass. Gnedin [153] showed that the gas content of halos in simulations is well fit by an expression which depends on the filtering mass, a particular time-averaged Jeans mass (Gnedin & Hui 1998 [150]). Gnedin [153] calculated the Jeans and filtering masses using the mean temperature in the simulation to define the sound speed, and found the following fit to the simulation results:

(138)

where _g is the average gas mass of all objects with a total mass M, f_b = _b / _m is the universal baryon fraction, and the characteristic mass M_c is the total mass of objects which on average retain 50% of their gas mass. The characteristic mass was well fit by the filtering mass at a range of redshifts from z = 4 up to z ~ 15.

The reionization process was not perfectly synchronized throughout the Universe. Large-scale regions with a higher density than the mean tend to form galaxies first and reionize earlier than underdense regions (see detailed discussion in §168). The suppression of low-mass galaxies by reionization will therefore be modulated by the fluctuations in the timing of reionization. Babich & Loeb (2005) [14] considered the effect of inhomogeneous reionization on the power-spectrum of low-mass galaxies. They showed that the shape of the high redshift galaxy power spectrum on small scales in a manner which depends on the details of epoch of reionization. This effect is significantly larger than changes in the galaxy power spectrum due to the current uncertainty in the inflationary parameters, such as the tilt of the scalar power spectrum n and the running of the tilt . Therefore, future high redshift galaxies surveys hoping to constrain inflationary parameters must properly model the effects of reionization, but conversely they will also be sensitive to the thermal history of the high redshift intergalactic medium.

8. FEEDBACK FROM GALACTIC OUTFLOWS

8.1. Propagation of Supernova Outflows in the IGM

Star formation is accompanied by the violent death of massive stars in supernova explosions. In general, if each halo has a fixed baryon fraction and a fixed fraction of the baryons turns into massive stars, then the total energy in supernovae outflows is proportional to the halo mass. The binding energy of the gas in the halo is proportional to the halo mass squared. Thus, outflows are expected to escape more easily out of low-mass galaxies, and to expel a greater fraction of the gas from dwarf galaxies. At high redshifts, most galaxies form in relatively low-mass halos, and the high halo merger rate leads to vigorous star formation. Thus, outflows may have had a great impact on the earliest generations of galaxies, with consequences that may include metal enrichment of the IGM and the disruption of dwarf galaxies. In this subsection we present a simple model for the propagation of individual supernova shock fronts in the IGM. We discuss some implications of this model, but we defer to the following subsection the brunt of the discussion of the cosmological consequences of outflows.

For a galaxy forming in a given halo, the supernova rate is related to the star formation rate. In particular, for a Scalo (1998) [315] initial stellar mass function, if we assume that a supernova is produced by each M > 8 M star, then on average one supernova explodes for every 126 M of star formation, expelling an ejecta mass of ~ 3 M including ~ 1 M of heavy elements. We assume that the individual supernovae produce expanding hot bubbles which merge into a single overall region delineated by an outwardly moving shock front. We assume that most of the baryons in the outflow lie in a thin shell, while most of the thermal energy is carried by the hot interior. The total ejected mass equals a fraction f_gas of the total halo gas which is lifted out of the halo by the outflow. This gas mass includes a fraction f_eject of the mass of the supernova ejecta itself (with f_eject 1 since some metals may be deposited in the disk and not ejected). Since at high redshift most of the halo gas is likely to have cooled onto a disk, we assume that the mass carried by the outflow remains constant until the shock front reaches the halo virial radius. We assume an average supernova energy of 10⁵¹ E₅₁ erg, a fraction f_wind of which remains in the outflow after it escapes from the disk. The outflow must overcome the gravitational potential of the halo, which we assume to have a Navarro, Frenk, & White (1997) [266] density profile [NFW; see equation (88)]. Since the entire shell mass must be lifted out of the halo, we include the total shell mass as well as the total injected energy at the outset. This assumption is consistent with the fact that the burst of star formation in a halo is typically short compared to the total time for which the corresponding outflow expands.

The escape of an outflow from an NFW halo depends on the concentration parameter c_N of the halo. Simulations by Bullock et al. (2000) [72] indicate that the concentration parameter decreases with redshift, and their results may be extrapolated to our regime of interest (i.e., to smaller halo masses and higher redshifts) by assuming that

(139)

Although we calculate below the dynamics of each outflow in detail, it is also useful to estimate which halos can generate large-scale outflows by comparing the kinetic energy of the outflow to the potential energy needed to completely escape (i.e., to infinite distance) from an NFW halo. We thus find that the outflow can escape from its originating halo if the circular velocity is below a critical value given by

(140)

where the efficiency is the fraction of baryons incorporated in stars, and

(141)

Note that the contribution to f_gas of the supernova ejecta itself is 0.024 f_eject, so the ejecta mass is usually negligible unless f_gas < 1%. Equation (140) can also be used to yield the maximum gas fraction f_gas which can be ejected from halos, as a function of their circular velocity. Although this equation is most general, if we assume that the parameters f_gas and f_wind are independent of M and z then we can normalize them based on low-redshift observations. If we specify c_N ~ 10 (with g(10) = 6.1) at z = 0, then setting E₅₁ = 1 and = 10% yields the required energy efficiency as a function of the ejected halo gas fraction:

(142)

A value of V_crit ~ 100 km s^-1 is suggested by several theoretical and observational arguments which are discussed in the next subsection. However, these arguments are not conclusive, and V_crit may differ from this value by a large factor, especially at high redshift (where outflows are observationally unconstrained at present). Note the degeneracy between f_gas and f_wind which remains even if V_crit is specified. Thus, if V_crit ~ 100 km s^-1 then a high efficiency f_wind ~ 1 is required to eject most of the gas from all halos with V_c < V_crit, but only f_wind ~ 10% is required to eject 5 - 10% of the gas. The evolution of the outflow does depend on the value of f_wind and not just the ratio f_wind / f_gas, since the shell accumulates material from the IGM which eventually dominates over the initial mass carried by the outflow.

We solve numerically for the spherical expansion of a galactic outflow, elaborating on the basic approach of Tegmark, Silk, & Evrard (1993) [358]. We assume that most of the mass m carried along by the outflow lies in a thin, dense, relatively cool shell of proper radius R. The interior volume, while containing only a fraction f_int << 1 of the mass m, carries most of the thermal energy in a hot, isothermal plasma of pressure p_int and temperature T. We assume a uniform exterior gas, at the mean density of the Universe (at each redshift), which may be neutral or ionized, and may exert a pressure p_ext as indicated below. We also assume that the dark matter distribution follows the NFW profile out to the virial radius, and is at the mean density of the Universe outside the halo virial radius. Note that in reality an overdense distribution of gas as well as dark matter may surround each halo due to secondary infall.

The shell radius R in general evolves as follows:

(143)

where the right-hand-side includes forces due to pressure, sweeping up of additional mass, gravity, and a cosmological constant, respectively. The shell is accelerated by internal pressure and decelerated by external pressure, i.e., p = p_int - p_ext. In the gravitational force, M(R) is the total enclosed mass, not including matter in the shell, and 1/2m is the effective contribution of the shell mass in the thin-shell approximation [279]. The interior pressure is determined by energy conservation, and evolves according to [358]:

(144)

where the luminosity L incorporates heating and cooling terms. We include in L the supernova luminosity L_sn (during a brief initial period of energy injection), cooling terms L_cool, ionization L_ion, and dissipation L_diss. For simplicity, we assume ionization equilibrium for the interior plasma, and a primordial abundance of hydrogen and helium. We include in L_cool all relevant atomic cooling processes in hydrogen and helium, i.e., collisional processes, Bremsstrahlung emission, and Compton cooling off the CMB. Compton scattering is the dominant cooling process for high-redshift outflows. We include in L_ion only the power required to ionize the incoming hydrogen upstream, at the energy cost of 13.6 eV per hydrogen atom. The interaction between the expanding shell and the swept-up mass dissipates kinetic energy. The fraction f_d of this energy which is re-injected into the interior depends on complex processes occurring near the shock front, including turbulence, non-equilibrium ionization and cooling, and so (following Tegmark et al. 1993 [358]) we let

(145)

where we set f_d = 1 and compare below to the other extreme of f_d = 0.

In an expanding Universe, it is preferable to describe the propagation of outflows in terms of comoving coordinates since, e.g., the critical result is the maximum comoving size of each outflow, since this size yields directly the total IGM mass which is displaced by the outflow and injected with metals. Specifically, we apply the following transformation [328, 374]:

(146)

For = 0, Voit (1996) [374] obtained (with the time origin = 0 at redshift z₁):

(147)

while for _m + = 1 there is no simple analytic expression. We set = / _vir, in terms of the virial radius r_vir [equation (84)] of the source halo. We define _s¹ as the ratio of the shell mass m to 4/3 _{b vir}³, where _b = _b(z = 0) is the mean baryon density of the Universe at z = 0. More generally, we define

(148)

Here we assumed, as noted above, that the shell mass is constant until the halo virial radius is reached, at which point the outflow begins to sweep up material from the IGM. We thus derive the following equations:

(149)

along with

(150)

In the evolution equation for , for < 1 we assume for simplicity that the baryons are distributed in the same way as the dark matter, since in any case the dark matter halo dominates the gravitational force. For > 1, however, we correct (via the last term on the right-hand side) for the presence of mass in the shell, since at >> 1 this term may become important. The > 1 equation also includes the braking force due to the swept-up IGM mass. The enclosed mean overdensity for the NFW profile [Eq. (88)] surrounded by matter at the mean density is

(151)

The physics of supernova shells is discussed in Ostriker & McKee (1988) [279] along with a number of analytical solutions. The propagation of cosmological blast waves has also been computed by Ostriker & Cowie (1981) [278], Bertschinger (1985) [40] and Carr & Ikeuchi (1985) [74]. Voit (1996) [374] derived an exact analytic solution to the fluid equations which, although of limited validity, is nonetheless useful for understanding roughly how the outflow size depends on several of the parameters. The solution requires an idealized case of an outflow which at all times expands into a homogeneous IGM. Peculiar gravitational forces, and the energy lost in escaping from the host halo, are neglected, cooling and ionization losses are also assumed to be negligible, and the external pressure is not included. The dissipated energy is assumed to be retained, i.e., f_d is set equal to unity. Under these conditions, the standard Sedov self-similar solution [324, 325] generalizes to the cosmological case as follows [374]:

(152)

where = 2.026 and ₀ = E₀ / (1 + z₁)² in terms of the initial (i.e., at t = = 0 and z = z₁) energy E₀. Numerically, the comoving radius is

(153)

In solving the equations described above, we assume that the shock front expands into a pre-ionized region which then recombines after a time determined by the recombination rate. Thus, the external pressure is included initially, it is turned off after the pre-ionized region recombines, and it is then switched back on at a lower redshift when the Universe is reionized. When the ambient IGM is neutral and the pressure is off, the shock loses energy to ionization. In practice we find that the external pressure is unimportant during the initial expansion, although it is generally important after reionization. Also, at high redshift ionization losses are much smaller than losses due to Compton cooling. In the results shown below, we assume an instantaneous reionization at z = 9.

Figure 51 shows the results for a starting redshift z = 15, for a halo of mass 5.4 × 10⁷ M, stellar mass 8.0 × 10⁵ M, comoving r_vir = 12 kpc, and circular velocity V_c = 20 km/s. We show the shell comoving radius in units of the virial radius of the source halo (top panel), and the physical peculiar velocity of the shock front (bottom panel). Results are shown (solid curve) for the standard set of parameters f_int = 0.1, f_d = 1, f_wind = 75%, and f_gas = 50%. For comparison, we show several cases which adopt the standard parameters except for no cooling (dotted curve), no reionization (short-dashed curve), f_d = 0 (long-dashed curve), or f_wind = 15% and f_gas = 10% (dot-short dashed curve). When reionization is included, the external pressure halts the expanding bubble. We freeze the radius at the point of maximum expansion (where d / d = 0), since in reality the shell will at that point begin to spread and fill out the interior volume due to small-scale velocities in the IGM. For the chosen parameters, the bubble easily escapes from the halo, but when f_wind and f_gas are decreased the accumulated IGM mass slows down the outflow more effectively. In all cases the outflow reaches a size of 10 - 20 times r_vir, i.e., 100 - 200 comoving kpc. If all the metals are ejected (i.e., f_eject = 1), then this translates to an average metallicity in the shell of ~ 1 - 5 × 10^-3 in units of the solar metallicity (which is 2% by mass). The asymptotic size of the outflow varies roughly as f_wind^1/5, as predicted by the simple solution in equation (152), but the asymptotic size is rather insensitive to f_gas (at a fixed f_wind) since the outflow mass becomes dominated by the swept-up IGM mass once > 4_vir. With the standard parameter values (i.e., those corresponding to the solid curve), Figure 51 also shows (dot-long dashed curve) the Voit (1996) [374] solution of equation (152). The Voit solution behaves similarly to the no-reionization curve at low redshift, although it overestimates the shock radius by ~ 30%, and the overestimate is greater compared to the more realistic case which does include reionization.

Figure 51. Evolution of a supernova outflow from a z = 15 halo of circular velocity Vc = 20 km/s. Plotted are the shell comoving radius in units of the virial radius of the source halo (top panel), and the physical peculiar velocity of the shock front (bottom panel). Results are shown for the standard parameters fint = 0.1, fd = 1, fwind = 75%, and fgas = 50% (solid curve). Also shown for comparison are the cases of no cooling (dotted curve), no reionization (short-dashed curve), fd = 0 (long-dashed curve), or fwind = 15% and fgas = 10% (dot-short dashed curve), as well as the simple Voit (1996) [374] solution of equation (152) for the standard parameter set (dot-long dashed curve). In cases where the outflow halts, we freeze the radius at the point of maximum expansion.

Figure 52 shows different curves than Figure 51 but on an identical layout. A single curve starting at z = 15 (solid curve) is repeated from Figure 51, and it is compared here to outflows with the same parameters but starting at z = 20 (dotted curve), z = 10 (short-dashed curve), and z = 5 (long-dashed curve). A V_c = 20 km/s halo, with a stellar mass equal to 1.5% of the total halo mass, is chosen at the three higher redshifts, but at z = 5 a V_c = 42 km/s halo is assumed. Because of the suppression of gas infall after reionization, we assume that the z = 5 outflow is produced by supernovae from a stellar mass equal to only 0.3% of the total halo mass (with a similarly reduced initial shell mass), thus leading to a relatively small final shell radius. The main conclusion from both figures is the following: In all cases, the outflow undergoes a rapid initial expansion over a fractional redshift interval z / z ~ 0.2, at which point the shell has slowed down to ~ 10 km/s from an initial 300 km/s. The rapid deceleration is due to the accumulating IGM mass. External pressure from the reionized IGM completely halts all high-redshift outflows, and even without this effect most outflows would only move at ~ 10 km/s after the brief initial expansion. Thus, it may be possible for high-redshift outflows to pollute the Lyman alpha forest with metals without affecting the forest hydrodynamically at z < 4. While the bulk velocities of these outflows may dissipate quickly, the outflows do sweep away the IGM and create empty bubbles. The resulting effects on observations of the Lyman alpha forest should be studied in detail (some observational signatures of feedback have been suggested recently by Theuns, Mo, & Schaye 2000 [362]).

Figure 52. Evolution of supernova outflows at different redshifts. The top and bottom panels are arranged similarly to Figure 51. The z = 15 outflow (solid curve) is repeated from Figure 51, and it is compared here to outflows with the same parameters but starting at z = 20 (dotted curve), z = 10 (short-dashed curve), and z = 5 (long-dashed curve). A Vc = 20 km/s halo is assumed except for z = 5, in which case a Vc = 42 km/s halo is assumed to produce the outflow (see text).

Furlanetto & Loeb (2003) [141] derived the evolution of the characteristic scale and filling fraction of supernova-driven bubbles based on a refinement of this formalism (see also their 2001 paper for quasar-driven outflows). The role of metal-rich outflows in smearing the transition epoch between Pop-III (metal-free) and Pop II (metal-enriched) stars, was also analysed by Furlanetto & Loeb (2005) [144], who concluded that a double-reionization history in which the ionization fraction goes through two (or more) peaks is unlikely.

8.2. Effect of Outflows on Dwarf Galaxies and on the IGM

Galactic outflows represent a complex feedback process which affects the evolution of cosmic gas through a variety of phenomena. Outflows inject hydrodynamic energy into the interstellar medium of their host galaxy. As shown in the previous subsection, even a small fraction of this energy suffices to eject most of the gas from a dwarf galaxy, perhaps quenching further star formation after the initial burst. At the same time, the enriched gas in outflows can mix with the interstellar medium and with the surrounding IGM, allowing later generations of stars to form more easily because of metal-enhanced cooling. On the other hand, the expanding shock waves may also strip gas in surrounding galaxies and suppress star formation.

Dekel & Silk (1986) [104] attempted to explain the different properties of diffuse dwarf galaxies in terms of the effect of galactic outflows. They noted the observed trends whereby lower-mass dwarf galaxies have a lower surface brightness and metallicity, but a higher mass-to-light ratio, than higher mass galaxies. They argued that these trends are most naturally explained by substantial gas removal from an underlying dark matter potential. Galaxies lying in small halos can eject their remaining gas after only a tiny fraction of the gas has turned into stars, while larger galaxies require more substantial star formation before the resulting outflows can expel the rest of the gas. Assuming a wind efficiency f_wind ~ 100%, Dekel & Silk showed that outflows in halos below a circular velocity threshold of V_crit ~ 100 km/s have sufficient energy to expel most of the halo gas. Furthermore, cooling is very efficient for the characteristic gas temperatures associated with V_crit < 100 km/s halos, but it becomes less efficient in more massive halos. As a result, this critical velocity is expected to signify a dividing line between bright galaxies and diffuse dwarf galaxies. Although these simple considerations may explain a number of observed trends, many details are still not conclusively determined. For instance, even in galaxies with sufficient energy to expel the gas, it is possible that this energy gets deposited in only a small fraction of the gas, leaving the rest almost unaffected.

Since supernova explosions in an inhomogeneous interstellar medium lead to complicated hydrodynamics, in principle the best way to determine the basic parameters discussed in the previous subsection (f_wind, f_gas, and f_eject) is through detailed numerical simulations of individual galaxies. Mac Low & Ferrara (1999) [235] simulated a gas disk within a z = 0 dark matter halo. The disk was assumed to be azimuthally symmetric and initially smooth. They represented supernovae by a central source of energy and mass, assuming a constant luminosity which is maintained for 50 million years. They found that the hot, metal-enriched ejecta can in general escape from the halo much more easily than the colder gas within the disk, since the hot gas is ejected in a tube perpendicular to the disk without displacing most of the gas in the disk. In particular, most of the metals were expelled except for the case with the most massive halo considered (with 10⁹ M in gas) and the lowest luminosity (10³⁷ erg/s, or a total injection of 2 × 10⁵² erg). On the other hand, only a small fraction of the total gas mass was ejected except for the least massive halo (with 10⁶ M in gas), where a luminosity of 10³⁸ erg/s or more expelled most of the gas. We note that beyond the standard issues of numerical resolution and convergence, there are several difficulties in applying these results to high-redshift dwarf galaxies. Clumping within the expanding shells or the ambient interstellar medium may strongly affect both the cooling and the hydrodynamics. Also, the effect of distributing the star formation throughout the disk is unclear since in that case several characteristics of the problem will change; many small explosions will distribute the same energy over a larger gas volume than a single large explosion [as in the Sedov (1959) [324] solution; see, e.g., equation (152)], and the geometry will be different as each bubble tries to dig its own escape route through the disk. Also, high-redshift disks should be denser by orders of magnitude than z = 0 disks, due to the higher mean density of the Universe at early times. Thus, further numerical simulations of this process are required in order to assess its significance during the reionization epoch.

Some input on these issues also comes from observations. Martin (1999) [247] showed that the hottest extended X-ray emission in galaxies is characterized by a temperature of ~ 10^6.7 K. This hot gas, which is lifted out of the disk at a rate comparable to the rate at which gas goes into new stars, could escape from galaxies with rotation speeds of < 130 km/s. However, these results are based on a small sample which includes only the most vigorous star-forming local galaxies, and the mass-loss rate depends on assumptions about the poorly understood transfer of mass and energy among the various phases of the interstellar medium.

Many authors have attempted to estimate the overall cosmological effects of outflows by combining simple models of individual outflows with the formation rate of galaxies, obtained via semi-analytic methods [98, 358, 374, 265, 129, 317] or numerical simulations [148, 149, 80, 9]. The main goal of these calculations is to explain the characteristic metallicities of different environments as a function of redshift. For example, the IGM is observed to be enriched with metals at redshifts z < 5. Identification of C IV, Si IV an O VI absorption lines which correspond to Ly absorption lines in the spectra of high-redshift quasars has revealed that the low-density IGM has been enriched to a metal abundance (by mass) of Z_IGM ~ 10^{-2.5 ( ± 0.5)} Z, where Z = 0.019 is the solar metallicity [252, 372, 347, 229, 99, 346, 121]. The metal enrichment has been clearly identified down to H I column densities of ~ 10^14.5 cm^-2. The detailed comparison of cosmological hydrodynamic simulations with quasar absorption spectra has established that the forest of Ly absorption lines is caused by the smoothly-fluctuating density of the neutral component of the IGM [84, 405, 180]. The simulations show a strong correlation between the H I column density and the gas overdensity _gas [102], implying that metals were dispersed into regions with an overdensity as low as _gas ~ 3 or possibly even lower.

In general, dwarf galaxies are expected to dominate metal enrichment at high-redshift for several reasons. As noted above and in the previous subsection, outflows can escape more easily out of the potential wells of dwarfs. Also, at high redshift, massive halos are rare and dwarf halos are much more common. Finally, as already noted, the Sedov (1959) [324] solution [or equation (152)] implies that for a given total energy and expansion time, multiple small outflows fill large volumes more effectively than would a smaller number of large outflows. Note, however, that the strong effect of feedback in dwarf galaxies may also quench star formation rapidly and reduce the efficiency of star formation in dwarfs below that found in more massive galaxies.

Cen & Ostriker (1999) [80] showed via numerical simulation that metals produced by supernovae do not mix uniformly over cosmological volumes. Instead, at each epoch the highest density regions have much higher metallicity than the low-density IGM. They noted that early star formation occurs in the most overdense regions, which therefore reach a high metallicity (of order a tenth of the solar value) by z ~ 3, when the IGM metallicity is lower by 1 - 2 orders of magnitude. At later times, the formation of high-temperature clusters in the highest-density regions suppresses star formation there, while lower-density regions continue to increase their metallicity. Note, however, that the spatial resolution of the hydrodynamic code of Cen & Ostriker is a few hundred kpc, and anything occurring on smaller scales is inserted directly via simple parametrized models. Scannapieco & Broadhurst (2000) [317] implemented expanding outflows within a numerical scheme which, while not a full gravitational simulation, did include spatial correlations among halos. They showed that winds from low-mass galaxies may also strip gas from nearby galaxies (see also Scannapieco, Ferrara, & Broadhurst 2000 [318]), thus suppressing star formation in a local neighborhood and substantially reducing the overall abundance of galaxies in halos below a mass of ~ 10¹⁰ M. Although quasars do not produce metals, they may also affect galaxy formation in their vicinity via energetic outflows [116, 15, 339, 263].

Gnedin & Ostriker (1997) [148] and Gnedin (1998) [149] identified another mixing mechanism which, they argued, may be dominant at high redshift (z > 4). In a collision between two protogalaxies, the gas components collide in a shock and the resulting pressure force can eject a few percent of the gas out of the merger remnant. This is the merger mechanism, which is based on gravity and hydrodynamics rather than direct stellar feedback. Even if supernovae inject most of their metals in a local region, larger-scale mixing can occur via mergers. Note, however, that Gnedin's (1998) [149] simulation assumed a comoving star formation rate at z > 5 of ~ 1 M per year per comoving Mpc³, which is 5 - 10 times larger than the observed rate at redshift 3 - 4. Aguirre et al. [9] used outflows implemented in simulations to conclude that winds of ~ 300 km/s at z < 6 can produce the mean metallicity observed at z ~ 3 in the Ly forest. In a separate paper Aguirre et al. [10] explored another process, where metals in the form of dust grains are driven to large distances by radiation pressure, thus producing large-scale mixing without displacing or heating large volumes of IGM gas. The success of this mechanism depends on detailed microphysics such as dust grain destruction and the effect of magnetic fields. The scenario, though, may be directly testable because it leads to significant ejection only of elements which solidify as grains.

Feedback from galactic outflows encompasses a large variety of processes and influences. The large range of scales involved, from stars or quasars embedded in the interstellar medium up to the enriched IGM on cosmological scales, make possible a multitude of different, complementary approaches, promising to keep galactic feedback an active field of research.

9. THE FRONTIER OF 21CM COSMOLOGY

9.1. Mapping Hydrogen Before Reionization

The small residual fraction of free electrons after cosmological recombination coupled the temperature of the cosmic gas to that of the cosmic microwave background (CMB) down to a redshift, z ~ 200 [284]. Subsequently, the gas temperature dropped adiabatically as T_gas (1 + z)² below the CMB temperature T (1 + z). The gas heated up again after being exposed to the photo-ionizing ultraviolet light emitted by the first stars during the reionization epoch at z < 20. Prior to the formation of the first stars, the cosmic neutral hydrogen must have resonantly absorbed the CMB flux through its spin-flip 21cm transition [131, 323, 367, 404]. The linear density fluctuations at that time should have imprinted anisotropies on the CMB sky at an observed wavelength of = 21.12[(1 + z) / 100] meters. We discuss these early 21cm fluctuations mainly for pedagogical purposes. Detection of the earliest 21cm signal will be particularly challenging because the foreground sky brightness rises as ^2.5 at long wavelengths in addition to the standard ^1/2 scaling of the detector noise temperature for a given integration time and fractional bandwidth. The discussion in this section follows Loeb & Zaldarriaga (2004) [226].

Figure 53. The 21cm transition of hydrogen. The higher energy level the spin of the electron (e-) is aligned with that of the proton (p+). A spin flip results in the emission of a photon with a wavelength of 21cm (or a frequency of 1420MHz).

We start by calculating the history of the spin temperature, T_s, defined through the ratio between the number densities of hydrogen atoms in the excited and ground state levels, n₁ / n₀ = (g₁/ g₀)exp{-T_* / T_s},

(154)

where subscripts 1 and 0 correspond to the excited and ground state levels of the 21cm transition, (g₁ / g₀) = 3 is the ratio of the spin degeneracy factors of the levels, n_H = (n₀ + n₁) (1 + z)³ is the total hydrogen density, and T_* = 0.068K is the temperature corresponding to the energy difference between the levels. The time evolution of the density of atoms in the ground state is given by,

(155)

where a(t) = (1 + z)^-1 is the cosmic scale factor, A's and B's are the Einstein rate coefficients, C's are the collisional rate coefficients, and I is the blackbody intensity in the Rayleigh-Jeans tail of the CMB, namely I = 2kT / ² with = 21 cm [306]. Here a dot denotes a time-derivative. The 0 1 transition rates can be related to the 1 0 transition rates by the requirement that in thermal equilibrium with T_s = T = T_gas, the right-hand-side of Eq. (155) should vanish with the collisional terms balancing each other separately from the radiative terms. The Einstein coefficients are A₁₀ = 2.85 × 10^-15 s^-1, B₁₀ = (³ / 2hc) A₁₀ and B₀₁ = (g₁ / g₀)B₁₀ [131, 306]. The collisional de-excitation rates can be written as C₁₀ = 4/3 (1 - 0) n_H, where (1 - 0) is tabulated as a function of T_gas [11, 406].

Equation (155) can be simplified to the form,

(156)

where n₀ / n_H, H H₀ (_m)^1/2(1 + z)^3/2 is the Hubble parameter at high redshifts (with a present-day value of H₀), and _m is the density parameter of matter. The upper panel of Fig. 54 shows the results of integrating Eq. (156). Both the spin temperature and the kinetic temperature of the gas track the CMB temperature down to z ~ 200. Collisions are efficient at coupling T_s and T_gas down to z ~ 70 and so the spin temperature follows the kinetic temperature around that redshift. At much lower redshifts, the Hubble expansion makes the collision rate subdominant relative the radiative coupling rate to the CMB, and so T_s tracks T again. Consequently, there is a redshift window between 30 < z < 200, during which the cosmic hydrogen absorbs the CMB flux at its resonant 21cm transition. Coincidentally, this redshift interval precedes the appearance of collapsed objects [23] and so its signatures are not contaminated by nonlinear density structures or by radiative or hydrodynamic feedback effects from stars and quasars, as is the case at lower redshifts [404].

Figure 54. Upper panel: Evolution of the gas, CMB and spin temperatures with redshift [4]. Lower panel: dTb / dH as function of redshift. The separate contributions from fluctuations in the density and the spin temperature are depicted. We also show dTb / dH a dTb / dH × H, with an arbitrary normalization.

During the period when the spin temperature is smaller than the CMB temperature, neutral hydrogen atoms absorb CMB photons. The resonant 21cm absorption reduces the brightness temperature of the CMB by,

(157)

where the optical depth for resonant 21cm absorption is,

(158)

Small inhomogeneities in the hydrogen density _H (n_H - _H) / _H result in fluctuations of the 21cm absorption through two separate effects. An excess of neutral hydrogen directly increases the optical depth and also alters the evolution of the spin temperature. For now, we ignore the additional effects of peculiar velocities (Bharadwaj & Ali 2004 [41]; Barkana & Loeb 2004 [27]) as well as fluctuations in the gas kinetic temperature due to the adiabatic compression (rarefaction) in overdense (underdense) regions [29]. Under these approximations, we can write an equation for the resulting evolution of fluctuations,

(159)

leading to spin temperature fluctuations,

(160)

The resulting brightness temperature fluctuations can be related to the derivative,

(161)

The spin temperature fluctuations T_s / T_s are proportional to the density fluctuations and so we define,

(162)

through T_b = (d T_b / d _H) _H. We ignore fluctuations in C_ij due to fluctuations in T_gas which are very small [11]. Figure 54 shows dT_b / d_H as a function of redshift, including the two contributions to dT_b / d_H, one originating directly from density fluctuations and the second from the associated changes in the spin temperature [323]. Both contributions have the same sign, because an increase in density raises the collision rate and lowers the spin temperature and so it allows T_s to better track T_gas. Since _H grows with time as _H a, the signal peaks at z ~ 50, a slightly lower redshift than the peak of dT_b / d_H.

Next we calculate the angular power spectrum of the brightness temperature on the sky, resulting from density perturbations with a power spectrum P(k),

(163)

where _H(k) is the Fourier tansform of the hydrogen density field, k is the comoving wavevector, and < … > denotes an ensemble average (following the formalism described in [404]). The 21cm brightness temperature observed at a frequency corresponding to a distance r along the line of sight, is given by

(164)

where n denotes the direction of observation, W(r) is a narrow function of r that peaks at the distance corresponding to . The details of this function depend on the characteristics of the experiment. The brightness fluctuations in Eq. 164 can be expanded in spherical harmonics with expansion coefficients a_lm(). The angular power spectrum of map C_l() = < |a_lm()|² > can be expressed in terms of the 3D power spectrum of fluctuations in the density P(k),

(165)

Our calculation ignores inhomogeneities in the hydrogen ionization fraction, since they freeze at the earlier recombination epoch (z ~ 10³) and so their amplitude is more than an order of magnitude smaller than _H at z < 100. The gravitational potential perturbations induce a redshift distortion effect that is of order ~ (H / ck)² smaller than _H for the high - l modes of interest here.

Figure 55. Angular power spectrum of 21cm anisotropies on the sky at various redshifts. From top to bottom, z = 55,40,80,30,120,25,170.

Figure 55 shows the angular power spectrum at various redshifts. The signal peaks around z ~ 50 but maintains a substantial amplitude over the full range of 30 < z < 100. The ability to probe the small scale power of density fluctuations is only limited by the Jeans scale, below which the dark matter inhomogeneities are washed out by the finite pressure of the gas. Interestingly, the cosmological Jeans mass reaches its minimum value, ~ 3 × 10⁴ M, within the redshift interval of interest here which corresponds to modes of angular scale ~ arcsecond on the sky. During the epoch of reionization, photoionization heating raises the Jeans mass by several orders of magnitude and broadens spectral features, thus limiting the ability of other probes of the intergalactic medium, such as the Ly forest, from accessing the same very low mass scales. The 21cm tomography has the additional advantage of probing the majority of the cosmic gas, instead of the trace amount (~ 10^-5) of neutral hydrogen probed by the Ly forest after reionization. Similarly to the primary CMB anisotropies, the 21cm signal is simply shaped by gravity, adiabatic cosmic expansion, and well-known atomic physics, and is not contaminated by complex astrophysical processes that affect the intergalactic medium at z < 30.

Characterizing the initial fluctuations is one of the primary goals of observational cosmology, as it offers a window into the physics of the very early Universe, namely the epoch of inflation during which the fluctuations are believed to have been produced. In most models of inflation, the evolution of the Hubble parameter during inflation leads to departures from a scale-invariant spectrum that are of order 1/N_efold with N_efold ~ 60 being the number of e - folds between the time when the scale of our horizon was of order the horizon during inflation and the end of inflation [218]. Hints that the standard CDM model may have too much power on galactic scales have inspired several proposals for suppressing the power on small scales. Examples include the possibility that the dark matter is warm and it decoupled while being relativistic so that its free streaming erased small-scale power [48], or direct modifications of inflation that produce a cut-off in the power on small scales [192]. An unavoidable collisionless component of the cosmic mass budget beyond CDM, is provided by massive neutrinos (see [198] for a review). Particle physics experiments established the mass splittings among different species which translate into a lower limit on the fraction of the dark matter accounted for by neutrinos of f > 0.3 %, while current constraints based on galaxies as tracers of the small scale power imply f < 12 % [360].

Figure 56 shows the 21cm power spectrum for various models that differ in their level of small scale power. It is clear that a precise measurement of the 21cm power spectrum will dramatically improve current constraints on alternatives to the standard CDM spectrum.

Figure 56. Upper panel: Power spectrum of 21cm anisotropies at z = 55 for a CDM scale-invariant power spectrum, a model with n = 0.98, a model with n = 0.98 and r 1/2 (d2lnP / dln k2) = -0.07, a model of warm dark matter particles with a mass of 1 keV, and a model in which f = 10% of the matter density is in three species of massive neutrinos with a mass of 0.4 eV each. Lower panel: Ratios between the different power spectra and the scale-invariant spectrum.

The sensitivity of an experiment depends strongly on its particular design, involving the number and distribution of the antennae for an interferometer. Crudely speaking, the uncertainty in the measurement of [l(l + 1)C_l / 2]^1/2 is dominated by noise, N, which is controlled by the sky brightness I at the observed frequency [404],

(166)

where l_min is the minimum observable l as determined by the field of view of the instruments, l_max is the maximum observable l as determined by the maximum separation of the antennae, f_cover is the fraction of the array area thats is covered by telescopes, t₀ is the observation time and is the frequency range over which the signal can be detected. Note that the assumed sky temperature of 0.7 × 10⁴ K at = 50 MHz (corresponding to z ~ 30) is more than six orders of magnitude larger than the signal. We have already included the fact that several independent maps can be produced by varying the observed frequency. The numbers adopted above are appropriate for the inner core of the LOFAR array (http://www.lofar.org), planned for initial operation in 2006. The predicted signal is ~ 1 mK, and so a year of integration or an increase in the covering fraction are required to observe it with LOFAR. Other experiments whose goal is to detect 21cm fluctuations from the subsequent epoch of reionization at z ~ 6-12 (when ionized bubbles exist and the fluctuations are larger) include the Mileura Wide-Field Array (MWA; http://web.haystack.mit.edu/arrays/MWA/), the Primeval Structure Telescope (PAST; http://arxiv.org/abs/astro-ph/0502029), and in the more distant future the Square Kilometer Array (SKA; http://www.skatelescope.org). The main challenge in detecting the predicted signal from higher redshifts involves its appearance at low frequencies where the sky noise is high. Proposed space-based instruments [194] avoid the terrestrial radio noise and the increasing atmospheric opacity at < 20 MHz (corresponding to z > 70).

Figure 57. Prototype of the tile design for the Mileura Wide-Field Array (MWA) in western Australia, aimed at detecting redshifted 21cm from the epoch of reionization. Each 4m×4m tile contains 16 dipole antennas operating in the frequency range of 80 - 300MHz. Altogether the initial phase of MWA (the so-called "Low-Frequency Demostrator") will include 500 antenna tiles with a total collecting area of 8000 m2 at 150MHz, scattered across a 1.5 km region and providing an angular resolution of a few arcminutes.

The 21cm absorption is replaced by 21cm emission from neutral hydrogen as soon as the intergalactic medium is heated above the CMB temperature by X-ray sources during the epoch of reionization [88]. This occurs long before reionization since the required heating requires only a modest amount of energy, ~ 10^-2 eV[(1 + z) / 30], which is three orders of magnitude smaller than the amount necessary to ionize the Universe. As demonstrated by Chen & Miralda-Escude (2004) [88], heating due the recoil of atoms as they absorb Ly photons [237] is not effective; the Ly color temperature reaches equilibrium with the gas kinetic temperature and suppresses subsequent heating before the level of heating becomes substantial. Once most of the cosmic hydrogen is reionized at z_reion, the 21cm signal is diminished. The optical depth for free-free absorption after reionization, ~ 0.1 [(1 + z_reion) / 20]^5/2, modifies only slightly the expected 21cm anisotropies. Gravitational lensing should modify the power spectrum [287] at high l, but can be separated as in standard CMB studies (see [326] and references therein). The 21cm signal should be simpler to clean as it includes the same lensing foreground in independent maps obtained at different frequencies.

Figure 58. Schematic sketch of the evolution of the kinetic temperature (Tk) and spin temperature (Ts) of cosmic hydrogen. Following cosmological recombination at z ~ 103, the gas temperature (orange curve) tracks the CMB temperature (blue line; T (1 + z)) down to z ~ 200 and then declines below it (Tk (1 + z)2) until the first X-ray sources (accreting black holes or exploding supernovae) heat it up well above the CMB temperature. The spin temperature of the 21cm transition (red curve) interpolates between the gas and CMB temperatures. Initially it tracks the gas temperature through collisional coupling; then it tracks the CMB through radiative coupling; and eventually it tracks the gas temperature once again after the production of a cosmic background of UV photons between the Ly and the Lyman-limit frequencies that redshift or cascade into the Ly resonance (through the Wouthuysen-Field effect [Wouthuysen 1952 [388]; Field 1959 [131]]). Parts of the curve are exaggerated for pedagogic purposes. The exact shape depends on astrophysical details about the first galaxies, such as their production of X-ray binaries, supernovae, nuclear accreting black holes, and their generation of relativistic electrons in collisionless shocks which produce UV and X-ray photons through inverse-Compton scattering of CMB photons.

The large number of independent modes probed by the 21cm signal would provide a measure of non-Gaussian deviations to a level of ~ N_{21 cm}^-1/2, constituting a test of the inflationary origin of the primordial inhomogeneities which are expected to possess deviations > 10^-6 [245].

9.2. The Characteristic Observed Size of Ionized Bubbles at the End of Reionization

Figure 59. 21cm imaging of ionized bubbles during the epoch of reionization is analogous to slicing swiss cheese. The technique of slicing at intervals separated by the typical dimension of a bubble is optimal for revealing different pattens in each slice.

The first galaxies to appear in the Universe at redshifts z > 20 created ionized bubbles in the intergalactic medium (IGM) of neutral hydrogen (H I) left over from the Big-Bang. It is thought that the ionized bubbles grew with time, surrounded clusters of dwarf galaxies [67, 143] and eventually overlapped quickly throughout the Universe over a narrow redshift interval near z ~ 6. This event signaled the end of the reionization epoch when the Universe was a billion years old. Measuring the unknown size distribution of the bubbles at their final overlap phase is a focus of forthcoming observational programs aimed at highly redshifted 21cm emission from atomic hydrogen. In this sub-section we follow Wyithe & Loeb (2004) [399] and show that the combined constraints of cosmic variance and causality imply an observed bubble size at the end of the overlap epoch of ~ 10 physical Mpc, and a scatter in the observed redshift of overlap along different lines-of-sight of ~ 0.15. This scatter is consistent with observational constraints from recent spectroscopic data on the farthest known quasars. This result implies that future radio experiments should be tuned to a characteristic angular scale of ~ 0.5° and have a minimum frequency band-width of ~ 8 MHz for an optimal detection of 21cm flux fluctuations near the end of reionization.

During the reionization epoch, the characteristic bubble size (defined here as the spherically averaged mean radius of the H II regions that contain most of the ionized volume [143]) increased with time as smaller bubbles combined until their overlap completed and the diffuse IGM was reionized. However the largest size of isolated bubbles (fully surrounded by H I boundaries) that can be observed is finite, because of the combined phenomena of cosmic variance and causality. Figure 61 presents a schematic illustration of the geometry. There is a surface on the sky corresponding to the time along different lines-of-sight when the diffuse (uncollapsed) IGM was most recently neutral. We refer to it as the Surface of Bubble Overlap (SBO). There are two competing sources for fluctuations in the SBO, each of which is dependent on the characteristic size, R_SBO, of the ionized regions just before the final overlap. First, the finite speed of light implies that 21cm photons observed from different points along the curved boundary of an H II region must have been emitted at different times during the history of the Universe. Second, bubbles on a comoving scale R achieve reionization over a spread of redshifts due to cosmic variance in the initial conditions of the density field smoothed on that scale. The characteristic scale of H II bubbles grows with time, leading to a decline in the spread of their formation redshifts [67] as the cosmic variance is averaged over an increasing spatial volume. However the 21cm light-travel time across a bubble rises concurrently. Suppose a signal 21cm photon which encodes the presence of neutral gas, is emitted from the far edge of the ionizing bubble. If the adjacent region along the line-of-sight has not become ionized by the time this photon reaches the near side of the bubble, then the photon will encounter diffuse neutral gas. Other photons emitted at this lower redshift will therefore also encode the presence of diffuse neutral gas, implying that the first photon was emitted prior to overlap, and not from the SBO. Hence the largest observable scale of H II regions when their overlap completes, corresponds to the first epoch at which the light crossing time becomes larger than the spread in formation times of ionized regions. Only then will the signal photon leaving the far side of the HII region have the lowest redshift of any signal photon along that line-of-sight.

Figure 60. Spectra of 19 quasars with redshifts 5.74 < z < 6.42 from the Sloan Digital Sky Survey [128]. For some of the highest-redshift quasars, the spectrum shows no transmitted flux shortward of the Ly wavelength at the quasar redshift (the so-called "Gunn-Peterson trough"), indicating a non-negligible neutral fraction in the IGM (see the analysis of Fan et al. [128] for details).

The observed spectra of some quasars beyond z ~ 6.1 show a Gunn-Peterson trough [163, 127] (Fan et al. 2005 [128]), a blank spectral region at wavelengths shorter than Ly at the quasar redshift, implying the presence of H I in the diffuse IGM. The detection of Gunn-Peterson troughs indicates a rapid change [126, 288, 381] in the neutral content of the IGM at z ~ 6, and hence a rapid change in the intensity of the background ionizing flux. This rapid change implies that overlap, and hence the reionization epoch, concluded near z ~ 6. The most promising observational probe [404, 259] of the reionization epoch is redshifted 21cm emission from intergalactic H I. Future observations using low frequency radio arrays (e.g. LOFAR, MWA, and PAST) will allow a direct determination of the topology and duration of the phase of bubble overlap. In this section we determine the expected angular scale and redshift width of the 21cm fluctuations at the SBO theoretically, and show that this determination is consistent with current observational constraints.

We start by quantifying the constraints of causality and cosmic variance. First suppose we have an H II region with a physical radius R / (1 + < z >). For a 21cm photon, the light crossing time of this radius is