© ESO, 2011

1. Introduction

The period of formation of the very first stars is one of the least understood epochs in the history of the Universe. The 21 cm spin-flip transition of neutral hydrogen at high redshifts has the potential to open a new observational window to study this early period, where the most distant, first galaxies reside, and even beyond that. Moreover, 21 cm observations may provide the only means to detect signatures of the first galaxies at z > 15 in the foreseeable future. While the James Webb Space Telescope (JWST¹) will be able to detect galaxies at z ≤ 10 (e.g., Stiavelli et al. 2004), it will not image all of the sources responsible for the reionization of the Universe. If current estimates in the literature are correct, even deep fields with JWST will only detect  ~30% of the sources responsible for reionization or maintaining reionization at z ~ 12 (Salvaterra et al. 2010). At z ~ 20, direct emission from stars will be below the threshold of JWST observations and new techniques will be required to probe the onset of reionization. One possibility is Lyman-α emitter surveys over the next decade, but even these surveys will only be able to detect the brightest galaxies for selected redshift ranges at z ≤ 10 (e.g., Ouchi et al. 2008; Nilsson et al. 2007; Cuby et al. 2007; Stark et al. 2007; Willis et al. 2008; McMahon et al. 2008; Hibon et al. 2010), but not the fainter, first galaxies at redshift z ≥ 20. The next generation of large ground based telescopes, like the European Extremely Large Telescope (E-ELT²), Thirty Meter Telescope (TMT³), and Giant Magellan Telescope (GMT⁴), will only be capable of detecting at most the very brightest galaxies at redshifts z ≈ 20. Most of the star formation however will be in fainter galaxies making it impossible for these telescopes alone to survey the full population. In combination with 21 cm observations, which are sensitive to the integrated Lyα emission of the galaxies, it should be possible to form a complete census of the first sources.

With a wide frequency coverage, observations of the 21 cm signal will provide a 3-dimensional tomographic view of the inter-galactic medium (IGM), before and during the cosmological reionization (for a review see Furlanetto et al. 2006). As soon as the first galaxies appear, at a redshift of z ~ 20−30 in the standard cold dark matter cosmological model (e.g., Komatsu et al. 2011), photons emitted at frequencies between the Lyα and Lyman limit quickly couple the spin temperature of the hyperfine level population to the IGM gas temperature through the Wouthuysen-Field effect (Wouthuysen 1952; Field 1959). Provided the IGM has been cooling adiabatically (X-ray heating is still negligible during this early stage of galaxy formation), the 21 cm signal will be strong and observed in absorption against the cosmic microwave background (CMB). Because the Lyα coupling depends on the local Lyα radiation, which is correlated with galaxies, the 21 cm signal fluctuations will trace the location, mass, emission spectra, luminosity and redshift evolution of high redshift galaxies.

Since the first galaxies represent rare high-σ peaks in the cosmic density field (e.g., Trac & Cen 2007), their spatial distribution have large fluctuations and we expect the Lyα contribution to the overall 21 cm signal to be appreciable, considerably improving the prospects for detection of the Lyα fluctuations through 21 cm observations. Although the high redshift-end of the reionization process will be inaccessible to the first generation experiments, such as the LOw Frequency ARray (LOFAR⁵) and the Murchison Widefield Array (MWA⁶), the second generation experiments such as the Square Kilometre Array (SKA⁷) should have enough collecting area to statistically probe the signal.

On the theoretical front, analytical models have been useful in predicting the probable evolution of the high redshift signal (Barkana & Loeb 2005b; Pritchard & Furlanetto 2007; Pritchard & Loeb 2008), but rely on linear approximations and require further comparison with simulations. Numerical simulations, on the other hand, can provide a self-consistent treatment of the Lyα radiative transfer, taking into account effects such as the scattering in the Lyα line wings (Semelin et al. 2007; Baek et al. 2009), although they are typically slow to run and are limited to small volumes (≤100 Mpc/h). Recent developments using semi-numerical algorithms have allowed the rapid generation of the high redshift 21 cm signal during the pre-reionization epoch (Santos et al. 2008, 2010; Mesinger et al. 2010) in large boxes, while maintaining the 3-d structure of the signal on large scales as seen in the full numerical simulations.

In this paper we will generate fast simulations as described in Santos et al. (2010) to analyze the dependence of the overall 21 cm signal on several parameters related to the first galaxies in the Universe and consider the ability of an experiment like SKA to constrain these parameters. The layout of this paper is as follows. We start in Section 2 by describing the 21 cm signal and the free parameters of the model that affect the Lyα fluctuations. In Sect. 3 we present the experimental setup used for SKA, calculating the expected error on the 3-d power spectrum. In Sect. 4 we study the possibility of constraining the signal in a model independent way, using the redshift-space distortions and the corresponding measurements on the 3-d power spectrum as a function of the angle with respect to the line of sight. We end with our conclusions and a discussion of the prospects for SKA in Sect. 5.

Throughout this paper where cosmological parameters are required we use the standard set of values Ω_m = 0.28, Ω_Λ = 0.72, Ω_b = 0.046, H = 100   h   km   s^-1   Mpc^-1 (with h = 0.70), n_S = 0.96, and σ₈ = 0.82, consistent with the latest measurements (Komatsu et al. 2011).

2. The 21 cm signal: Lyα fluctuations

The 21 cm brightness temperature corresponds to the change in the intensity of the CMB radiation due to absorption or emission when it travels through a patch of neutral hydrogen. It is given, at an observed frequency ν in the direction , by (see e.g. Santos et al. 2008) $\begin{array}{ccc}\mathit{\delta }{\mathit{T}}_{\mathrm{b}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}& \mathrm{\approx }& \mathrm{23}{\mathit{x}}_{\mathit{HI}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{)}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathit{\gamma }}}{{\mathit{T}}_{\mathrm{S}}}\right){\left(\frac{\mathit{h}}{\mathrm{0.7}}\right)}^{-1}\left(\frac{{\mathrm{\Omega }}_{\mathrm{b}}{\mathit{h}}^{\mathrm{2}}}{\mathrm{0.02}}\right)\\ & & \mathrm{\times }{\left[\left(\frac{\mathrm{0.15}}{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{h}}^{\mathrm{2}}}\right)\left(\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{10}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\mathit{H}\hspace{0.17em}\mathrm{d}{\mathit{v}}_{\mathit{r}}\mathit{/}\mathrm{d}\mathit{r}}\right)\mathrm{mK}\mathit{,}\end{array}$(1)where x_HI is the fraction of neutral hydrogen (mass weighted), dv_r/dr is the comoving gradient of the line of sight component of the comoving velocity and we use δ_a for the fractional value of the quantity a (${\mathit{\delta }}_{\mathit{a}}\mathrm{\equiv }\frac{\mathit{a}\mathrm{-}\mathrm{⟨}\mathit{a}\mathrm{⟩}}{\mathrm{⟨}\mathit{a}\mathrm{⟩}}$) with δ for the fluctuation in the matter density.

The spin temperature (T_S) is coupled to the hydrogen gas temperature (T_K) through the spin-flip transition, which can be excited by collisions or by the absorption of Lyα photons (Wouthuysen-Field effect) and we can write: $\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathit{\gamma }}}{{\mathit{T}}_{\mathrm{S}}}\mathrm{=}\frac{{\mathit{x}}_{\mathrm{tot}}}{\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathrm{tot}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathit{\gamma }}}{{\mathit{T}}_{\mathit{K}}}\right)\mathit{,}$(2)where x_tot = x_α + x_c is the sum of the radiative and collisional coupling parameters and we are already assuming that the color temperature of the Lyα radiation field at the Lyα frequency is equal to T_K. Note that although we talk here about Lyα radiation, in effect we consider the contribution from all the photons up to the Lyman limit frequency, since photons redshifting into Lyman series resonances can produce Lyα photons as a result of atomic cascades (Hirata 2006; Pritchard & Furlanetto 2006). When the coupling to the gas temperature is negligible (e.g. x_tot ~ 0), T_S ~ T_γ and there is no signal. On the other hand, for large x_tot, T_S simply follows T_K.

2.1. Simulations

We use the code SimFast21⁸ described in Santos et al. (2010) to calculate the spin temperature and the corresponding 21 cm signal. The code starts by calculating the linear density field in a box followed by the corresponding halo and velocity fields. Then, it computes the nonlinear density and halo fields, the collapsed mass distribution, the corresponding star formation rate, the IGM gas temperature, the Lyα coupling, the collisional coupling and finally the 21 cm brightness temperature accounting for all the contributions including the correction due to redshift space distortions using the nonlinear velocity field. In order to probe the full range of k space, we generated two end to end simulations using the SimFast21 code, one large simulation with 1Gpc in size and 1.6 Mpc resolution (S1) and another with higher resolution (0.186 Mpc) but smaller 143 Mpc in size (simulation S2). This code has the advantage of being modular, allowing to replace some of the steps by simulation boxes generated by other techniques. Therefore, for comparison, we also applied our Lyα calculation to the output (star formation rate, matter density) from the Trac et al. (2008) simulation, thus providing a check of the approximation used to generate the halo mass function. We call this simulation of the 21 cm signal, S3, with a size of 143 Mpc and a resolution of 0.186 Mpc. Table 1 shows a summary of the simulations used. All simulations use halos with masses down to 10⁸    M_⊙ corresponding to a minimum virial temperature of T ~ 10⁴ K. Simulation S1 cannot resolve halos down to these mass scales and we populate the remaining halos in each (empty) cell using a Poisson sampling biased with the underlying density field, giving a mass function consistent with N-body simulations once non-linear corrections are applied, as described in Santos et al. (2010). With the new version of the code, we also allow a hybrid approach instead of the “Poisson” method, where, after the halo finding algorithm, we calculate the unresolved collapsed mass in each empty cell using the collapsed fraction f_coll from the extended Press-Schechter formula (see e.g. Zahn et al. 2007). This is fine as long as we do not want to resolve halos within a cell. Note that this is used at the cell level, allowing to also apply the non-linear corrections to each cell. Both approaches give similar results in terms of the 21 cm signal, although the later method can be faster as we decrease the halo minimum mass used in the simulation.

Simulation S3 was run using a hybrid simulation code in modeling cosmic reionization (Trac et al. 2008), incorporating N-body, hydrodynamic, and RT algorithms to solve the coupled evolution of the dark matter, baryons, and radiation (Trac & Pen 2004, 2006; Trac & Cen 2007). It covers a large dynamic range and satisfies the requirements of having sufficiently high resolution to capture small-scale structure, and simultaneously a sufficiently large volume to reduce sample variance (Barkana & Loeb 2004). The simulation is run in two steps. The first step involves running a high-resolution N-body simulation with 3072³ dark matter particles on an effective mesh with 11   520³ cells in a comoving box, 143 Mpc on a side. We identified collapsed dark matter halos on the fly using a friends-of-friends algorithm, with a linking length b = 0.2 times the mean inter-particle spacing, in order to model radiation sources and sinks. With a particle mass resolution of 2.68 × 10⁶   h^-1    M_⊙, we can reliably locate all dark matter halos with virial temperatures above the atomic cooling limit (T ~ 10⁴ K) with a minimum of  ~40 particles (Heitmann et al. 2008), and half of this collapsed mass budget is resolved with  >400 particles per halo. Our halo mass functions are in very good agreement with other recent work (e.g. Reed et al. 2007; Lukic et al. 2007; Cohn & White 2008). The second step produces a hydrodynamic + RT simulation with moderate resolution, but incorporating subgrid physics modeled using the high-resolution information from the large N-body simulation. Radiation sources are prescribed and star formation rates calculated using the halo model described in Trac & Cen (2007). We consider only Population II stars from starbursts (Schaerer 2003a) as contributing to the ionizing photon budget. The hydro+RT simulation utilizes equal numbers (N = 1536³) of dark matter particles, gas cells, and adaptive rays, where for the latter, we track 5 frequencies above the hydrogen ionizing threshold of 13.6 eV. The photo-ionization and photo-heating rates for each cell are calculated from the incident radiation flux and used in the non-equilibrium solvers for the ionization and energy equations. The initial conditions are generated with a common white noise field and a linear transfer function calculated with CAMB (Lewis et al. 2000).

In order to calculate the Lyα coupling in a given cell, the algorithm assumes that the scattering rate is proportional to the total Lyman series flux arriving in that cell (see again Santos et al. 2010). This Lyman series flux is obtained through a 3-dimensional integration of the comoving photon emissivity, ϵ(x,ν,z) (defined as the number of photons emitted at position x, redshift z and frequency ν per comoving volume, per proper time and frequency), which is assumed proportional to the star formation rate: ${\mathit{ϵ}}_{\mathit{\alpha }}\mathrm{\left(}{x}\mathit{,\nu ,z}\mathrm{\right)}\mathrm{=}\mathrm{SFRD}\mathrm{\left(}{x}\mathit{,z}\mathrm{\right)}{\mathit{ϵ}}_{\mathrm{b}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{,}$(3)where SFRD(x,z) is the star formation rate density from the simulation (in terms of the number of baryons in stars per comoving volume and proper time). As stated above, in the simulation S3 we used the star formation rate obtained from the Trac et al. (2008) simulation. ϵ_b(ν) is the spectral distribution function of the sources (defined as the number of photons per unit frequency emitted at ν per baryon in stars). We assumed a power law model for ϵ_b(ν): ${\mathit{ϵ}}_{\mathrm{b}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathit{A}{\mathit{\nu }}^{\mathit{\alpha }}\mathit{,}$(4)between ν_α(10.2 eV) and the Lyman limit frequency (13.6 eV). This will allow us to calculate how sensitive the 21 cm power spectrum is to changes in the parameters of the emission model.

We take α = −0.9 and A such that the integration between the Lyα and Lyman limit frequencies gives a total emission of 20 000 photons per baryon. These numbers are based on the expected PopII spectra from Schaerer (2003b). Note that the parameter A will be totally degenerate with the star formation rate efficiency in our model. Since the initial mass function (IMF) is likely to be dominated by massive stars, this single power law model is likely to be a good description of the emission spectra (see e.g. Leitherer et al. 1999), although a broken power law can sometimes provide a slightly better fit (Pritchard & Furlanetto 2006). Moreover, as we shall see later, results are very insensitive to α, essentially depending on the total number of photons emitted per baryon. We also note that the above calculation of Lyα radiation does not take into account full radiative transfer effects via their scattering off of neutral hydrogen atoms (scattering in the wings of the Lyα line, Chuzhoy & Zheng 2007; Semelin et al. 2007, could increase the amplitude of the fluctuations on small scales). As a result, we expect that some details of the signal on small scales due to Lyα coupling to be partly inaccurate.

Since the details of the galaxy emissivity depend upon the assumed initial mass function (IMF), which is poorly known (Barkana & Loeb 2001), these values are subject to considerable uncertainty. Constraints would therefore provide useful information about the properties of the first galaxies. These uncertainties make the mapping between the mean value of the Lyα coupling  ⟨ x_α ⟩  and redshift uncertain. We will therefore describe the shape of the Lyα power spectrum for a given  ⟨ x_α ⟩  and caution that the redshift at which that actually occurs may be very different from our fiducial model.

2.2. Comparison

The 21 cm signal occurs over a wide range of scales, with Lyα fluctuations contributing at wave-numbers 0.01 ≲ k/(h/ Mpc) ≲ 10. Since none of our simulations can cover this broad set of scales individually, we exploit the modularity of the code to piece together two simulations to cover the full range. To check that this gives self-consistent results, in Fig. 1 we plot the full 21 cm brightness temperature for the three simulations at a few values of the mean Lyα coupling  ⟨ x_α ⟩ . In each case, we take into account fluctuations from the density, velocity and Lyα terms as well as the collisional coupling, ionization fraction (although at high redshifts this is negligible) and gas temperature (note again that for simulation S3 we calculate the temperature using our semi-numerical code based on the simulation SFRD). We see a smooth transition from the S1 to the S2 simulation, showing for the first time, the full range of scales from k ~ 0.01 h/Mpc to k ~ 20 h/Mpc with an amplitude that is roughly between 100 mK² and 1000 mK² on the range of interest, which should help the detection of the signal. The increase in the amplitude as the coupling increases is due primarily to the cooling of the IGM gas in this redshift range, which leads to a stronger absorption signal. Note that, as explained in Santos et al. (2010), when integrating the Lyα flux, we assume the star formation rate is homogeneous above a certain large scale in order to make the calculation faster. In this paper this scale was set to 143 Mpc which is the origin of the slight break we see in Fig. 1 for x_α ≳ 1. The black dash-dotted line in this figure shows the correction for x_α = 0.4 if instead we did the inhomogeneous integration up to 300 Mpc. Since this is a small correction, we decided to keep the limit to 143 Mpc.

For comparison, the dashed line shows the power spectrum from simulation S3, which agrees quite well with results using the semi-numerical code validating our halo/star formation prescription (although in S3 there is a decrease of the signal on very small scales, which is related to a smoothing of the velocity term on these scales). An extensive comparison made by Baek et al. (2010) using radiative transfer simulations (although for smaller box size  ~100 Mpc/h and lower mass resolution  ≳ 10¹⁰ M_⊙) shows a similar spectrum for models in the same parameter range.

Since it will be the focus of the present paper, we now examine the Lyα fluctuations in more detail. In Fig. 2 we show the power spectrum of the quantity ${\mathit{A}}_{\mathit{\alpha }}\mathrm{\equiv }\frac{{\mathit{x}}_{\mathit{\alpha }}}{\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{\alpha }}}$ for the three simulations at four different values of  ⟨ x_α ⟩ . A_α represents the contribution of the Lyα coupling to the 21 cm signal and varies between 0 and 1. Note that the x_α field tends to be dominated by the large values close to individual sources, which obscures the details of the coupling since in these regions A_α saturates. The amplitude of the power spectrum principally depends on  ⟨ x_α ⟩  with the fluctuations peaking at  ⟨ x_α ⟩  ~ 0.5, suggesting that this is the regime to best extract information related to the Wouthuysen-Field effect (z ~ 20 in our model), and we will concentrate on these redshifts from now on. Figure 2 shows again the correction to the calculation if we do not impose a limit to the Lyα inhomogeneous integration (top, cyan dashed line, for x_α = 0.9). This correction will be more evident as x_α increases (see the “break” for x_α = 8.4) but at the same time the contribution to the overall 21 cm signal will become smaller so that it should be safe to neglect this correction.

Comparing to Fig. 14 in Baek et al. (2009) we see that the fluctuation level is similar, although we do not see their small scale increase in power due to point sources for k > 1 Mpc/h. The difference on small scales (k > 1 Mpc/h) between our results and those of (Baek et al. 2009) is, at least in part, likely due to lack of Lyα radiative transfer treatment in ours (although differences in size and mass resolution of the simulations should also be considered). As will become clear later, it is the 21 cm observations at large scales that offer the best means to extract Lyα physics, so the inaccuracies at small scales do not pose a significant problem.

The full N-body simulation (S3) agrees reasonable well with our calculation and again we see a smooth transition from the very large scale simulation to the high resolution one, giving us a complete view over all scales of the expected signal. Although it is clear that there is room for more detailed numerical work on the 21 cm signal, these comparisons give us confidence in the semi-analytic simulation and we now turn to exploring the consequences for observations.

2.3. Identifying the Lyα epoch

In principle, there are several contributions to the 21 cm brightness temperature and full modeling of the signal would be required to disentangle the part due to Lyα fluctuations and obtain parameter constraints on the first galaxies. However, as we can see from Figs. 3 and 4, both the fluctuations in the collisional coupling and in the gas temperature should be subdominant at the redshifts where the Lyα fluctuations are strongest (when  ⟨ x_α ⟩  < 1 and k < 1 h/ Mpc) making the analysis of the signal more straightforward (the ionization is very small at these redshifts, so that their contribution to the fluctuation power is completely negligible). The fluctuations in the gas temperature shown in Fig. 4 only include the heating due to x-rays, which can potentially be harder to separate from the Lyα fluctuations. On top of this we also have fluctuations in the gas temperature due to the adiabatic cooling process which are essentially proportional to the matter density perturbations. Both contributions are intrinsically included in the simulation, although, as we see, the effect from x-ray heating is negligible.

The question then is if whether we can actually identify this epoch from observations rather than theoretical guesses. In many models, the global (average) 21 cm signal clearly shows the transition from Lyα coupling to X-ray heating and finally to emission once the gas is heated above the CMB temperature (Pritchard & Loeb 2010). This average signal will not be accessible to interferometers, but as was shown in Santos et al. (2008), Fig. 9, the large scale 21 cm power spectrum can also depict the evolution of the signal.

The rms of the signal (⟨ (δT_b −  ⟨ δT_b ⟩ )² ⟩ ) can also provide that information and should be easier to measure with interferometers. As we can see in Fig. 5, at z > 25 the signal is initially small since coupling to the gas temperature is low. As the Lyα coupling increases, the spin temperature approaches the gas temperature and the rms of the signal increases, not only due to fluctuations on the Lyα field but also because the factor A_α(1 − T_γ/T_K) is increasing (in absolute terms) since the gas is still cooling adiabatically. Once the gas temperature begins to increase due to X-ray heating the signal decreases and we see the turning point at z ~ 16 (the transition from absorption to emission occurs later at z ~ 13. Finally, the gas is heated far above the CMB temperature and the signal plateaus until reionization finally causes it to die away at z ~ 7.

Note that as X-rays start heating the IGM, the signal does not drop immediately since this heating is inhomogeneous and there will be some region still cooling adiabatically as well as temperature fluctuations contributing to the rms. Therefore we see a small plateau at the peak of the signal starting at z ~ 19. The conclusion is that it should be safe to neglect the contributions of the X-ray heating to the 21 cm signal between the points when the rms (or the average) starts rising at z ≳ 24 and before we reach the plateau at the maximum z ≲ 19 (x_α ~ 0.9), at least on large scales where the Lyα contribution dominates over the temperature as we can see in Fig. 4.

In this case, the equation simplifies and we can write the full 21 cm signal as: $\mathit{\delta }{\mathit{T}}_{\mathrm{b}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathit{C}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\delta }}_{{\mathit{A}}_{\mathit{\alpha }}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\delta }}_{{\mathit{A}}_{\mathit{v}}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta \delta }\mathrm{)}\mathit{,}$(5)where (note the different definition from Barkana & Loeb 2005b), corresponding to the adiabatic cooling process (T_K ∝ ρ^2/3) and $\mathit{T̅}{}_{\mathit{K}}\mathrm{\approx }\mathrm{180}{{}^{\mathrm{\right(}}\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{100}}^{\mathrm{\left)}}}^{\mathrm{2}}$ K. Also, $\begin{array}{ccc}\mathit{C}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{\approx }& \mathrm{23}\left(\frac{\mathrm{0.7}}{\mathit{h}}\right)\left(\frac{{\mathrm{\Omega }}_{\mathrm{b}}{\mathit{h}}^{\mathrm{2}}}{\mathrm{0.02}}\right){\left[\left(\frac{\mathrm{0.15}}{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{h}}^{\mathrm{2}}}\right)\left(\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{10}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }\mathrm{⟨}\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathit{\gamma }}}{{\mathit{T}}_{\mathit{K}}}\mathrm{⟩}\mathrm{⟨}{\mathit{A}}_{\mathit{\alpha }}\mathrm{⟩}\mathrm{⟨}{\mathit{A}}_{\mathit{v}}\mathrm{⟩}\mathrm{mK}\end{array}$(6)with ${\mathit{A}}_{\mathit{v}}\mathrm{=}{}^{\mathrm{\right(}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\mathit{H}\hspace{0.17em}\mathrm{d}{\mathit{v}}_{\mathit{r}}\mathit{/}\mathrm{d}\mathit{r}}^{\mathrm{\left)}}$ (⟨ A_v ⟩  ~ 1). Figure 6 shows the contribution to the full 21 cm signal from each of the terms considered on Eq. (5) using simulations S1 and S2. Again we see that the Lyα term dominates over large scales k < 1 h/Mpc, which makes it all important to be able to generate large volume simulations (>100 Mpc/h) and give us confidence that it should be possible to extract the Lyα signal over the redshift range proposed without confusion from other contributions. Moreover, the small contribution at these scales from the other terms in Eq. (5) can be modeled with reasonable accuracy if we know the matter power spectrum.

3. Measurements with the SKA

In this section, we outline details related to the experimental measurement of the 21 cm signal at the very high redshifts during the epoch where Lyα dominates.

3.1. Noise power spectrum

Before we go into the details of the experimental setup, we quickly review the noise calculation that goes into our analysis. Following Bowman et al. (2006); Mao et al. (2008); McQuinn et al. (2006), the expected error ΔP(k,θ) on the measurement of the 3-d power spectrum of the signal P_S(k,θ) with noise P_N(k,θ) is given by $\mathrm{\Delta }\mathit{P}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{{\mathit{N}}_{\mathrm{m}}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}}}\left[{\mathit{P}}_{\mathrm{S}}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathit{N}}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}\right]\mathit{,}$(7)where we are assuming that the power spectrum depends only on the moduli (k) of the vector k and on the angle θ between k and the line of sight. N_m(k,θ) is the total number of modes in k space contributing to the measurement (note that the sum is only done over half the sphere). In order to obtain the above expression, we “grid” the k space into pixels of size d${\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\mathit{d}{\mathit{k}}_{\mathrm{\parallel }}$, where  ⊥  denotes the component perpendicular to the line of sight and  ∥  the parallel one and assume that measurements in different pixels are uncorrelated. The size of these pixels is just determined by the real space volume of the observed region: $\mathrm{Vol}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathit{y}{\mathit{A}}_{\mathrm{FoV}}\mathit{B,}$(8)where A_FoV is the field of view of the experiment (in radians), B is the bandwidth for this particular measurement, $\mathit{r}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{z}}\mathit{c}{\mathit{H}}^{-1}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}$ is the comoving distance to redshift z and y ≡ λ₂₁(1 + z)²/H(z) is a conversion factor between frequency intervals and comoving distances (H(z) is the Hubble expansion rate and λ₂₁ ≈ 21.1 cm the rest frame wavelength of the 21 cm line). We can relate dk _⊥  to the resolution in the u − v (visibility) space, du (the Fourier dual of the angular coordinates on the sky) by dk _⊥  = 2πdu/r (using our Fourier conventions). Note that $\mathrm{d}\mathit{u}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{{\mathit{A}}_{\mathrm{FoV}}}$(9)and is related to the effective collecting area of one element of the interferometer, A_e = λ²du², where λ is the wavelength of observation (but care must be taken in this relation due to multi-beaming). The resolution along the line of sight is just dk _∥  = 2π/(yB) (visibilities are Fourier transformed along the frequency coordinate).

Finally, the noise power spectrum is given by ${\mathit{P}}_{\mathit{N}}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathit{y}\frac{{\mathit{T}}_{\mathrm{sys}}^{\mathrm{2}}}{\mathrm{d}{\mathit{u}}^{\mathrm{2}}\mathit{t}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}}\mathit{,}$(10)where t(k,θ) is the time spent observing a given k pixel which is just the time spent observing the corresponding pixel of resolution du² on the u − v plane, u = rk _⊥ /(2π). For a total observation time, t₀, this will be related to the baseline density distribution, n(|u|), through t(k,θ) = t₀du²n(|u|) = t₀du²n(rksin(θ)/(2π)) assuming the baseline density is rotationally invariant (which is a good approximation due to Earth rotation). If we make the further assumption that the baseline density distribution is constant on the u − v plane up to a maximum baseline D_max (which is not the same as assuming an uniform distribution of antennas), we have $\mathit{n}\mathrm{\left(}\mathrm{\right|}{u}\mathrm{\left|}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\lambda }}^{\mathrm{2}}{\mathit{N}}_{\mathit{a}}^{\mathrm{2}}}{\mathit{\pi }{\mathit{D}}_{\max }^{\mathrm{2}}}\mathit{,}$(11)noting that the integration over half the plane should give approximately ${\mathit{N}}_{\mathit{a}}^{\mathrm{2}}\mathit{/}\mathrm{2}$ (N_a is the number of elements in the interferometer). The noise power spectrum then further simplifies to: ${\mathit{P}}_{\mathit{N}}\mathrm{\left(}\mathit{k,\theta }\mathrm{\right)}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathit{y}\frac{\mathit{\pi }{\mathit{\lambda }}^{\mathrm{2}}{\mathit{D}}_{\max }^{\mathrm{2}}{\mathit{T}}_{\mathrm{sys}}^{\mathrm{2}}}{{\mathit{t}}_{\mathrm{0}}{\mathit{A}}_{\mathrm{tot}}^{\mathrm{2}}}\mathit{,}$(12)where A_tot is the total collecting area of the telescope. The system temperature was modeled by ${\mathit{T}}_{\mathrm{sys}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{50}\mathrm{+}\mathrm{60}{\left(\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{4.73}}\right)}^{\mathrm{2.55}}\mathrm{K}$(13)where the first term is the receiver noise temperature (assumed to be 50 K) and the second term is the much larger sky temperature which is dominated by the Galactic Synchrotron at the redshifts of interest. Note in particular that neither the bandwidth or the frequency resolution show up in the noise calculation (they can be used instead to increase the number of modes N_m available to reduce the total error). Moreover, looking at the last expression we see that P_N decrease as $\mathrm{1}\mathit{/}{\mathit{A}}_{\mathrm{tot}}^{\mathrm{2}}$ once we fix t₀ and the maximum baseline. Once we achieve P_S ~ P_N the only way to decrease the total error is through N_m which does not depend on the total collecting area or the baseline distribution.

3.2. Experimental setup

The current generation of radio-interferometers now in construction such as the MWA and LOFAR will not have the collecting area necessary to probe the high redshift universe (z > 15). This will require a second generation experiment such as the SKA, with an observational window between 70 MHz and 10 GHz and planned to be completed by 2020. The SKA will probably be made up of three different instruments (Schilizzi et al. 2007; Faulkner 2010): SKA-low, tailored for the low frequency signal between 70MHz to 450MHz, probably made up of sparse Aperture Arrays so that the collecting area scales as λ²; SKA-mid between 400 MHz and 1.4 GHz, using dense Aperture Arrays; and SKA-high, using dishes for frequencies above 1.2 GHz (note the small overlap between ranges). In this paper, we will assume a low frequency “SKA type” experiment with a setup capable of probing the high redshift, pre-reionization epoch.

The analysis in the previous Section showed that in order to measure the imprint of the first galaxies on the 21 cm signal, we need to observe at least up to z = 20, corresponding to a frequency of ν ~ 68 MHz. Allowing for some flexibility, we consider an experiment that would be able to measure frequencies down to 60 MHz. Taking into account the design reference for SKA, we assumed a sensitivity of 4000 m²/K at 70 MHz (which in fact requires a total collecting area of  ~14 km², well above the 1 km² used to baptize the telescope) and scaling as λ² around these frequencies. Note that the planned collecting area has been subject to some change along the years and the current design for what is called SKA Phase 1 suggests a sensitivity of 2000 m²/K instead (Garrett et al. 2010). In order to allow for designs with lower sensitivities, we also consider instruments with 20% and 10% of the total collecting area of the design reference SKA. For the distribution of antennas, the SKADS⁹ reference design assumes they would be collected in 250 stations of 180 m diameter each, with 66% in a 5 km diameter core and the rest along 5 spiral arms out to a 180 km radius.

The signal we are looking for dominates on scales k ≲ 1 h/Mpc, requiring baselines up to  ~10 km for proper sampling at the relevant frequencies, e.g. a resolution of 0.76 arcmin at z = 20 (k_max ⊥  ~ 1.8 h/Mpc – higher resolution can be achieved along the frequency direction although the error would be much larger since we would only have modes along the line of sight). Taking a reasonable setup, we assumed that 70% of the total collecting area would be concentrated within this core of 10 km in diameter. We do not use the rest of the collecting area in the noise calculation, assuming instead this would be used for point source removal and calibration. The same applies when the total collecting area is only 10% and 20% of the reference design.

For the field of view (FoV) we used 10 × 10 deg², which, taking z = 20 as reference, sets a resolution of roughly dk _⊥  ~ 4.5 × 10^-3 h/Mpc. Note that the FoV can be traded with the number of beams (in case we want to observe separate fields) and the total instantaneous bandwidth (in case we want to probe a larger redshift range). Assuming that, in the core, correlations within stations are also possible, we considered a minimum baseline of 50 m giving k_min ⊥  ~ 9 × 10^-3 h/Mpc which is enough to probe the scales of interest. The frequency interval for the analysis should be chosen carefully at these frequencies since for a fixed frequency interval, dz increases with redshift. For instance, an interval of 8 MHz used in previous analysis (e.g. Bowman et al. 2007) corresponds to integrating the signal over dz ~ 2.5 which will average out important features during the period when Lyα fluctuations are important. On the other hand, using a small interval means we will have fewer modes along the line of sight. Taking a compromise, we assumed that each measurement was done with an interval of 4 MHz and a resolution of 0.01 MHz, giving dk _∥  = 0.08 h/Mpc and k_max ∥  = 16.6 h/Mpc. Note that the total instantaneous bandwidth should be much larger than this (~50 MHz), so we can measure several redshifts at the same time. Finally, the total integration time was taken to be 1000 h (see Table 2 for a summary of the experimental setup).

We followed the error calculation described above for the power spectrum assuming a distribution for the antenna/stations of the telescope such that a constant baseline density is obtained. Although this is difficult to achieve in practice since the density of stations normally decreases from the center, it turns out that this configuration gives the best signal to noise in terms of the brightness temperature maps and it should therefore give a good indication of the telescope capabilities to probe the 21 cm signal power spectrum. A final note on foregrounds: they can be the most damaging factor at these low frequencies (together with the calibration issues raised by the ionosphere) although it is expected that if the foregrounds are smooth enough along the frequency direction, they can be simply removed by fitting out a smooth function (Jelić et al. 2008; Morales et al. 2006; Santos et al. 2005). In this analysis we assume that foreground cleaning was applied on a region of 32 MHz (so as to avoid edge effects on our 4 MHz interval) and was successful on scales smaller than this 32 MHz, e.g. we ignore scales with k < 0.01 h/Mpc (see Harker et al. 2010). Calibration issues raised by the ionosphere can also be quite damaging in particular due to the large field of view of this low frequency interferometers (Cohen & Röttgering 2009; Matejek & Morales 2009; Liu et al. 2010) and the situation should become more clear once we have results from the first generation of EoR experiments. In this paper we concentrate on analysing the required minimum setup for an SKA type experiment in order to make a statistical detection of the signal, neglecting for now the challenges posed by the ionosphere.

3.3. Power spectrum constraints

Figure 7 shows the expected error on the 3 dimensional power spectrum, P(k) (integrated over θ), assuming the above setup. Even if foreground removal affects larger k-modes than assumed here we still have a huge range of measurements available up to k ~ 7 h/Mpc and in particular, measurements are quite good on the interval where the Lyα contribution is more important. On very large scales it should be possible to measure the power spectrum even with only 10% of the collecting area. Sample variance dominates on large scales while noise (dashed green line) is dominating on small scales (the number of modes on small scales is limited due to the size of the maximum baselines). Note that we could redistribute the collecting area so that the noise power spectrum followed more closely the signal. The error on large scales is kept small because there are still a large number of measurements N_m due to the high resolution on the u − v space from the large field of view (see Eq. (7)).

With the above measurements of the 3d power spectra, we should be able to constrain some of the characteristics of the first galaxies such as their emission spectra. Although, as we have seen, the Lyα signal dominates on larger scales over other contributions, these constraints will be particularly strong if we fix the cosmology, which should be a reasonable assumption given the expected level of constraints from the Planck satellite.

Up to now we have been assuming a type II emission spectrum with a total of 20   000 photons per baryon emitted between the Lyα and Lyman limit frequency. If instead we assume that the stars responsible for the Lyα coupling are Pop III, we should have around 5000 photons emitted per baryon used in stars with a spectral index of α = 0.29. Figure 8 shows the results of changing the emission model with the dashed curves corresponding to the PopIII case (note that due to the lower number of photons, these lines also correspond to a lower  ⟨ x_α ⟩  of 0.1. The dotted curve corresponds to changing the spectral index to α = 0.29 while maintaining the total emission to 5000 photons, which will be difficult to distinguish under the current noise expectations (the total number of emitted photons per baryon is the most relevant parameter for the Lyα flux). Another possibility is to assume that Pop III stars are also being formed in halos with mass less than 10⁸    M_⊙. The green dot-dashed line in Fig. 8 shows the effect of considering halos down to 10⁶    M_⊙ from simulation S1 at the same redshift. Note that this has a much higher coupling and if we consider higher redshifts so that  ⟨ x_α ⟩  ~ 0.1 the line will be closer to the dashed one.

These calculations show that it should be possible for SKA to measure  ⟨ x_α ⟩  from observations of the 21 cm power spectrum at z ≈ 20. Given a detailed model and small error bars this could be measured directly from the shape of the power spectrum. More generally, if measured over only a small range of wave-numbers, the redshift dependence of  ⟨ x_α ⟩  could be extracted by looking at the redshift evolution of the amplitude and slope in a manner analogous to that suggested by Lidz et al. (2008) for reionization. Note that, although during Lyα domination different models have similar power spectra for the same  ⟨ x_α ⟩  , its evolution with redshift could be used to distinguish these models.

Measurements of  ⟨ x_α ⟩  in a number of different redshift bins would be instrumental in constructing a “21 cm Madau plot” (Madau et al. 1996). The Madau plot shows the evolution of the star formation rate as a function of redshift and producing such a plot in the very early Universe would be a major accomplishment of SKA. Measurements of  ⟨ x_α ⟩  do not give this directly, since in our model  ⟨ x_α ⟩  is determined by the product of the star-formation rate and the emissivity of the galaxies. Breaking this degeneracy observationally will be difficult and would lead to a plot whose absolute normalisation was uncertain, but whose shape tracked the star-formation rate (assuming no evolution of the galaxies). However, for many models of early galaxy formation (Leitherer et al. 1999; Bromm & Larson 2004) the variation in the number of Lyα photons produced per baryon is only a factor of two or so. Thus with a reasonable theoretical prior a useful measure of the star-formation history could be established.

4. Redshift space distortions

The signal we are trying to probe is actually anisotropic due to the redshift space distortions set by the peculiar velocity of the HI gas, which is already being taken into account in Eq. (5). Using a linear approximation, ${\mathit{A}}_{\mathit{\alpha }}\mathrm{\approx }\frac{\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{}^{\mathrm{(}}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{\mathit{\delta }}_{{\mathit{x}}_{\mathit{\alpha }}}{}^{\mathrm{)}}$ and ${\mathit{A}}_{\mathit{\nu }}\mathrm{\approx }\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathit{H}}\frac{\mathrm{d}{\mathit{v}}_{\mathit{r}}}{\mathrm{d}\mathit{r}}$ so to first order the fluctuation in T_b becomes $\mathit{\delta }{\mathit{T}}_{\mathrm{b}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{\approx }\mathit{C}\mathrm{\left(}\mathit{z}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{\delta }\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{\mathit{\delta }}_{{\mathit{x}}_{\mathit{\alpha }}}\mathrm{-}\frac{\mathrm{1}}{\mathit{H}}\frac{\mathrm{d}{\mathit{v}}_{\mathit{r}}}{\mathrm{d}\mathit{r}}\right)\mathrm{·}$(14)The peculiar velocity field originates from inhomogeneities in the density field, so that in Fourier space we can use Kaiser’s approximation $\frac{\mathrm{1}}{\mathit{H}}\frac{\mathrm{d}{\mathit{v}}_{\mathit{r}}}{\mathrm{d}\mathit{r}}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathit{\delta }\mathrm{\left(}{k}\mathrm{\right)}$ (Kaiser 1987) to write the velocity as a function of its dependence in μ = cos(θ). This linear approximation in the velocity contribution allows the brightness temperature power spectra to be separated in only three powers of μ: ${\mathit{P}}_{{\mathit{T}}_{\mathrm{b}}}\mathrm{\left(}\mathit{k,\mu }\mathrm{\right)}\mathrm{=}{\mathit{\mu }}^{\mathrm{4}}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{4}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{\mu }}^{\mathrm{2}}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{0}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(15)with ${\mathit{P}}_{{\mathit{\mu }}^{\mathrm{4}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{C}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{P}}_{\mathit{\delta }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(16)${\mathit{P}}_{{\mathit{\mu }}^{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathit{C}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\left[\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}{\mathit{P}}_{\mathit{\delta }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{\mathit{P}}_{\mathit{\delta }{\mathit{\delta }}_{{\mathit{x}}_{\mathit{\alpha }}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right]\mathit{,}$(17)and $\begin{array}{ccc}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{0}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& {\mathit{C}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}[\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }{\mathrm{)}}^{\mathrm{2}}{\mathit{P}}_{\mathit{\delta }}\mathrm{(}\mathit{k}\mathrm{)}\\ & & \end{array}$(18)where the power spectra P_ab(k) for generic quantities a and b is defined as: $\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}_{\mathit{D}}^{\mathrm{3}}\mathrm{\right(}{k}\mathrm{-}{{k}}^{\mathrm{\prime }}\mathrm{\left)}{\mathit{P}}_{\mathit{ab}}\mathrm{\right(}\mathit{k}\mathrm{)}\mathrm{\equiv }\frac{\mathrm{⟨}\mathit{a}\mathrm{\left(}{k}\mathrm{\right)}\mathit{b}\mathrm{(}{{k}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\mathit{b}\mathrm{(}{k}\mathrm{\left)}\mathit{a}\mathrm{\right(}{{k}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{⟩}}{\mathrm{2}}\mathrm{·}$(19)This decomposition can be particularly useful at the high redshifts we are considering, where, as can be seen from Fig. 6, the contribution from the velocity compression can be relatively strong. By fitting the observed signal to the above polynomial (Eq. (15)) at each k, we can in principle separate the cosmological and astrophysical contributions to the brightness temperature fluctuations (Barkana & Loeb 2005a,b). One can use P_μ² and P_μ⁰ to learn about the first sources of radiation while the P_μ⁴ term could be used to measure the matter power spectrum directly, with no interference from any other sources of fluctuations.

4.1. Assumptions on linearity

The decomposition discussed in the previous section relies on assumptions of linearity in both the velocity field and the Lyα fluctuations. On small scales, both of these assumptions may break down and we consider both in turn. As we go to small scales, fluctuations on the density field increase and the non-linearities in the velocity contribution can become important, so that the first order approximation to the full velocity term in Eq. (1) breaks down. Figure 9 compares the brightness temperature power spectrum from the full calculation with the case where we assume the first order approximation for the velocity. We see that for k ≳ 1.0 h/Mpc the non-linearities can become important. This result implies that for smaller scales the correct brightness temperature power spectra has a much more complicated dependence in μ than the one given by Eq. (15) which will complicate the model independent extraction of the astrophysical contributions (Shaw & Lewis 2008).

Fluctuations in the Lyα background too can leave the linear regime. The information from the Lyα background and consequently from the physical processes that originated it, is contained in δ_xα and it can only be directly obtained from the power spectrum angular separation in powers of μ for scales where the linear approximation of A_α is valid. To show the limits where the first order approximation can be used, we plotted in Fig. 10 the power spectra of A_α compared to its linear expression, using the values from simulation S1. The linear approximation should be valid for k ≲ 0.2 h/Mpc.

Looking at the histogram of δ_xα in Fig. 11 we see that there are few points where δ_xα > 1. Although this number is small compared to the overall distribution (most cells in the simulation have actually x_α ~ 0 at these high redshifts), they are enough to introduce non-linearities on the power spectrum calculation. This basically implies that for small scales the fluctuations in x_α will have a large dependence on the random distribution of the first galaxies so it will only be possible to do a linear approximation on A_α for scales large enough that a few small peaks of large intensity do not dominate the power spectra of the fluctuation. Some of the strong coupling regions coincide with ionized regions which may diminish the effect of strong fluctuations in the brightness temperature power spectra, however according with simulations S1, S2 and S3, at the high redshifts we are interested (z > 20) the ionized regions are so small that including this quantity has no meaningful effect in the linear approximation used.

From this study we can conclude that the brightness temperature power spectrum angular separation in Eq. (15) should be reasonable for k < 1 h/Mpc and the linear approximations in Eqs. (17) and (18) should be valid for larger scales, k ≲ 0.2 h/Mpc, allowing the direct extraction of information about the Lyα field (Barkana & Loeb 2005b). This is fortunate since it corresponds to the scales where the Lyα fluctuations dominate the contributions to the brightness temperature (Fig. 6) as well as the scales where the SKA should perform better (Fig. 7). For 0.2 ≲ k ≲ 1.0 h/Mpc we can still use the angular decomposition but Eqs. (17) and (18) will include higher order terms. Given the scales we are interested in, from now on we will concentrate on simulation S1 (1000 Mpc) which can be used to calculate the power spectrum up to k ~ 2.3 h/Mpc.

4.2. Constraints on P (k,μ)

As we have seen in the previous section, the redshift space distortions introduce an anisotropy in the 21 cm power spectrum that, in the linear regime, depends on k and even powers of μ. More generally non-linearities in the velocity field will introduce terms with more complicated dependence (Shaw & Lewis 2008). Measurements of P(k,μ) might provide more information on the Lyα signal (Barkana & Loeb 2005a) than the spherical average P(k) we considered so far. Figure 12 shows the expected errors in the measurements of the 3d power spectrum as a function of μ, for an SKA type experiment (Table 2). The calculation follows Eq. (7) where the number of modes in each bin is now a function of k and μ. We used bins with a logarithmic spacing in k and linear in μ (with 0 ≤ μ ≤ 1 since the power spectrum only depends on μ²). The number of modes in each bin is accounted for by gridding the k space with a resolution along the line of sight set by the interval in frequency used for the analysis and a resolution perpendicular to the line of sight set by the field of view (see Sect. 3.2). The minimum size of the bins is given by this resolution. Note that the number of independent μ bins increase with k, although, as can be seen from Fig. 12, the errors will also increase due to the experimental noise. Moreover, the power spectrum is quite flat for small k which will make it harder to extract the polynomial dependence (this is because the P_μ⁰ term is dominating in Eq. (15)).

Although difficult as we have seen, we can try to obtain P_μ⁰, P_μ² and P_μ⁴ by fitting Eq. (15) to the measured P(k,μ) for k ≲ 1.0 h/Mpc. This in turn should allow us to obtain direct constraints on a combination of P_δxαδxα and P_δδxα. For small k there are fewer bins so the uncertainty in P_μⁿ will be high even for those scales where the errors in P(k,μ) are small. The expected errors on the P_μⁿ term can be obtained using a Fisher matrix approach (Fisher 1935). For parameters a,b =  { P_μ⁰,P_μ²,P_μ⁴ } , the matrix is calculated through ${\mathit{F}}_{\mathit{ab}}\mathrm{\left(}{\mathit{k}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}}\frac{\mathrm{1}}{\mathrm{\Delta }\mathit{P}\mathrm{(}{\mathit{k}}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}}\frac{\mathit{\partial P}\mathrm{\left(}{\mathit{k}}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathit{j}}\mathrm{\right)}}{\mathit{\partial a}}\frac{\mathit{\partial P}\mathrm{\left(}{\mathit{k}}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathit{j}}\mathrm{\right)}}{\mathit{\partial b}}$(20)where again P(k,μ) is given by Eq. (15) and ΔP(k_i,μ_j) is the error for each bin. The covariance of the P_μⁿ terms is just the inverse of the Fisher matrix and Fig. 13 shows the expected 1-sigma errors ($\sqrt{\mathrm{(}{\mathit{F}}^{-1}{\mathrm{)}}_{\mathit{aa}}}$). We see that it should be possible to measure with reasonable accuracy both the P_μ⁰ and P_μ² terms which can already give us interesting constraints on the Lyα field. In fact, if we fix the cosmology (e.g. P_δ), then we should be able to extract information from P_δδxα and P_δxαδxα separately. On the other hand, the P_μ⁴ will be hard to measure, which will make the extraction of cosmological information more difficult, simply because it is an order of magnitude smaller than the other terms on the scales we are focusing here.

In order to improve this last constraint we would have to increase the collecting area of the instrument and the frequency interval used (4 MHz in this case) so to have more modes along the line of sight (note however that this will lead to cosmological evolution of the signal along the frequency bin, as already discussed). The SKA should be able to measure modes up to k ~ 10 h/Mpc which could have higher values of P_μ⁴, however for k > 1 h/Mpc the angular decomposition is no longer valid and the separate measurements of the P(k,μ) terms will be of little use (except maybe for foreground removal).

In this case, it might be better to just try to fit the parameters of the model directly to the averaged power spectrum P(k) using the simulations (as shown in Fig. 8). If on the other hand we assume we already know P_μ⁴ with reasonable accuracy then the errors on P_μ² will improve considerably (Fig. 13) and any degeneracies between P_μ⁰ and P_μ² will be broken.

4.3. Extracting the Poisson term

Finally, by measuring the P(k,μ) terms, we can in principle use an estimator to probe the properties of the first sources of radiation directly through the Lyα field, without having to go through a full, model dependent, parameter fit to the data. This is based on the decomposition proposed in Barkana & Loeb (2005a), valid in the linear regime. In that case, δ_xα can be expressed as a function of its dependence on effects correlated and uncorrelated with δ as: ${\mathit{\delta }}_{{\mathit{x}}_{\mathit{\alpha }}}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathit{W}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{\delta }\mathrm{\left(}{k}\mathrm{\right)}\mathrm{+}{\mathit{\delta }}_{\mathrm{p}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{.}$(21)Note that this decomposition is actually quite general and the main assumption here is the linearity on δ. δ_P is the Poisson contribution that arises from the statistical fluctuations in the number density of the rare first galaxies (galaxies being discrete objects and therefore not tracking the continuous density field perfectly) and which to first order is uncorrelated with the density fluctuations.

The window function, W(k), contains several effects by which δ_xα is related to the underlying matter density fluctuations, δ (see Barkana & Loeb 2005b for details). The strongest effect between δ and δ_xα is the biasing of the number density of galaxies with respect to the density fluctuations by a factor b(z) (the average halo bias) so that an overdense region will have a factor (1 + b(z)δ) more sources (Mo & White 1996). This is the only term that is implicitly included in our simulations through the way we do the 3d integration of the Lyα flux. Note however that this term completely dominates over the other terms in Barkana & Loeb (2005b) since b(z) ~ 15 at z ~ 20. In that case, we can express W(k) as: $\mathit{W}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}{\mathrm{\int }}_{\mathit{z}}^{{\mathit{z}}_{\max \mathrm{\left(}\mathrm{2}\mathrm{\right)}}}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\frac{\mathit{D}\mathrm{\left(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{D}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\frac{\mathrm{d}{\mathit{x}}_{\mathit{\alpha }}}{\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}}\mathit{b}\mathrm{\left(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\frac{\sin \mathrm{\left(}\mathit{kr}\mathrm{\right)}}{\mathit{kr}}\mathit{,}$(22)where D(z) is the linear growth function and we are just considering photons emitted below the Lyman beta limit for comparison, so that z_max(2) = (32/27)(1 + z) − 1 (Barkana & Loeb 2005b; Pritchard & Furlanetto 2007).

Using the parametrization of δ_xα we can write P_xα(k) = W²(k)P_δ(k) + P_p(k) and P_δδxα(k) = W(k)P_δ(k), meaning that we can isolate W(k) from the simulations (W(k) = P_δδxα(k)/P_δ(k)). In Fig. 14 we show the window function obtained from Eq. (22), with values from simulation S1, which has the same behavior as the one computed by Pritchard & Furlanetto (2007). Figure 14 also shows that the window function obtained from Eq. (22) has a higher amplitude for large scales and small scales than the one obtained from simulation S1 using P_δ and P_δδxα. On large scales this is probably a consequence of our box still not being quite large enough to encompass the rarest, most highly biased objects leading to an underestimate of the power on large scales. On small scales, this seems to be a consequence of the onset of non-linearities in the density and Lyα fields. On scales close enough to clusters of sources that the Lyα coupling saturates, we might expect to see a steeper fall off in the 21 cm power as seen here. Nevertheless, W(k) obtained from the simulation follows quite closely the expected theoretical values, which supports the decomposition given in Eq. (21) and validates the signal extraction we discuss next.

Also, in Fig. 15 we show P_δxα, W²P_δ and P_p, built with P_δ, P_δxα and P_δδxα from simulation S1. This shows the behavior expected from the analytic model of Barkana & Loeb (2005b) and gives an amplitude of P_p that is broadly consistent with the values that we calculate from their model. The Poisson contribution dominates for k > 0.2 h/Mpc so that there is just a small interval where we can expect to probe this term before non-linear effects correlated to the density field become relevant. Note that the Poisson contribution flattens on small scales k ≳ 1 (h/Mpc) slightly more than predicted by the analytic model. This is similar to the drop off in power on small scales seen in W(k) and probably occurs for similar reasons.

For large enough scales, we can use the observational data to obtain the window function through $\frac{{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{4}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}\mathrm{2}\mathrm{\times }\left[\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}\mathit{W}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right]\mathit{,}$(23)and the Poisson contribution uncorrelated with δ: ${\mathit{P}}_{\mathit{un}\mathrm{-}\mathit{\delta }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{0}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\frac{{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{2}}}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{4}{\mathit{P}}_{{\mathit{\mu }}^{\mathrm{4}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}{\mathit{C}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\left(\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{⟨}{\mathit{x}}_{\mathit{\alpha }}\mathrm{⟩}}\right)}^{\mathrm{2}}{\mathit{P}}_{\mathrm{p}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(24)The use of Eqs. (23) and (24) allows us to characterize the Lyα background and the number and distribution of galaxies at these early redshifts. In Fig. 16, we used simulation S1 to show P_un − δ and P_μ², which are consistent with the results obtained by Barkana & Loeb (2005b). Figure 17 shows the power spectra of the Poisson contribution (Eq. (24)) and the expected errors associated with this measurement. As expected, due to the high error in P_μ⁴, the measurement of this Poisson term will be quite difficult even with the assumed experimental setup (dashed line). If we assume that the matter power spectrum is known with reasonable accuracy from CMB data and galaxy surveys, then we can combine the measurement of the P_μ⁴ term at several k to constrain C²(z) (Eq. (16)), e.g. the mean brightness temperature on the sky. Unfortunately the constraint is still large, C²(z) ≈ (6.42 ± 5.60) × 10³mK², which translates into large errors for the Poisson term (dotted lines). Nevertheless, a detection will still be possible, with a signal to noise of 2.6 (dashed line) and 3.8 (dotted line).

If we further assume that we combine this result with the measurements from a global experiment in order to measure the mean brightness temperature with high accuracy, then we can extract all the relevant power spectra reasonably well (dot-dashed line).

5. Conclusions

In this paper, we have made use of a recently developed fast semi-numerical code to explore the possibility for 21 cm observations of the time of the first galaxies with SKA. We demonstrated that this code allows the simulation of the full range of scales likely to be accessible to observations. Here we have focused on the impact of 21 cm fluctuations from large fluctuations in the Lyα coupling, which arise due to the clustering of the rare first galaxies and considered different emission models.

Measurement of these fluctuations would provide insight into the formation of the very first galaxies and would be highly complementary with the next generation of large optical/IR telescopes. We have shown that SKA-pathfinders with  ~10% of the full collecting area should be capable of making a statistical detection of the 21 cm power spectrum at redshifts z ≲ 20. With the full SKA sensitivity this detection would become a measurement allowing astrophysical properties of the first galaxies to be determined. Note that it is crucial that these experiments are able to go below 70 MHz in order to probe the signal.

Observations with an SKA-like instrument would enable the determination of  ⟨ x_α ⟩  as a function of redshift from the amplitude and shape of the 21 cm fluctuations. Although different models with the same  ⟨ x_α ⟩  give similar power spectra, this evolution with redshift should allow to distinguish them. Even though there is a strong degeneracy between the star formation rate and the UV spectral properties of the galaxies themselves, these observations would enable constraints to be placed on the star formation rate in the earliest generations of galaxies. Although crude, the resulting “Madau plot” would be very useful for understanding the formation processes of the first galaxies.

Moving beyond the angle-averaged power spectrum, we have investigated the use of redshift-space distortions to separate out different components of the power spectrum as suggested by Barkana & Loeb (2005b). These measurements are difficult due to instrumental limitations and are further complicated by non-linearities on small scales. For the first time, we have investigated the non-linearity of Lyα fluctuations showing that they should be important on scales k ≳ 0.2 (h/Mpc). Nonetheless, on larger scales the detection of these effects is possible with SKA level sensitivity and should allow to put direct constraints on the Lyα power spectrum.

Redshift space distortions further offer the possibility of extracting components of the power spectrum that do not correlate with the density field. This Poissonian component contains direct information about the number density of sources and so gives complementary information to the other sources of fluctuations. We have shown that a detection is possible and if strong assumptions about the underlying cosmology are possible and combined with information about the mean 21 cm signal, then these fluctuations may be picked out clearly by SKA.

These results illustrate the potential of 21 cm observations to shed new light on the astrophysics during the pre-Reionization epoch. In the future, our improved knowledge of cosmological parameters will provide a firm foundation to pick out the details of galaxy formation in the early Universe. SKA and pathfinders capable of observing at frequencies ν ≲ 100   MHz will begin to access this interesting period and transform our understanding of the cosmic dawn.

---

Acknowledgments

This work was partially supported by FCT-Portugal under grants PTDC/FIS/66825/2006 and PTDC/FIS/100170/2008. J.R.P. is supported by NASA through Hubble Fellowship grant HST-HF-01211.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. A.C. acknowledges support of NSF grant CAREER AST-0645427. M.G.S. was a visitor at UCI when this work was concluded.

References

All Tables

All Figures