Light-cone anisotropy in the 21 cm signal from the epoch of reionization

Abstract

INTRODUCTION

The epoch of reionization (EoR) extends from redshift z = 20–30 down to redshift 6. During this era, the cold and neutral intergalactic medium (IGM) is progressively ionized by the light of the first stars and galaxies. To this day we have very little observational information about the state of the IGM during this process. The Gunn–Peterson absorption in high-redshift quasars (QSOs; Gunn & Peterson 1965) apparently shows that by z ∼ 6 the Universe is more than 99.9 per cent reionized (Fan et al. 2006), although the current observational sample may not be statistically significant (see Mesinger, Furlanetto & Cen 2011). The optical depth of Thomson scattering on free electrons measured from cosmic microwave background (CMB) observations favours an extended EoR [the best fit from combined Wilkinson Microwave Anisotropy Probe (WMAP) and Planck data is τ = 0.089 ± 0.014; Planck Collaboration 2013]. Over the next few years, 21 cm observations of the neutral IGM are the most promising type of observations likely to improve our understanding of reionization.

A number of projects are currently under way to detect the 21 cm signal from reionization. The reionization project at the Giant Metre-wave Radio Telescope has published an upper limit of (248 mK)² for the power spectrum at wavenumber k = 0.5 h Mpc⁻¹ at z = 8.7 (Paciga et al. 2013). A similar upper limit, (300 mK)² for comoving wavenumber k = 0.046 Mpc⁻¹ at z = 9.5, was found by the Murchison Widefield Array (Dillon et al. 2014). A stronger upper limit was found by the Precision Array for Probing the Epoch of Re-ionization (PAPER): (52 mK)² for k = 0.11 h Mpc⁻¹ at z = 7.7 (Parsons et al. 2013). The PAPER constraint implies that a small amount of heating occurs at z > 7.7 in the neutral IGM (a few K), since otherwise the signal would have been detected. This is not surprising: the opposite would imply that almost no X-rays were produced prior to these redshifts. The Low-Frequency Array (LOFAR) for radio astronomy is in operation but has not yet published its results. The primary goal of all the Square Kilometre Array (SKA) pathfinders is to achieve a statistical detection of the signal in the form of its three-dimensional (3D) power spectrum.

In the process, these pathfinders will have to deal with foreground contamination, which is a thousand or more times brighter than the signal. The power spectrum is easier to detect than the flux density, as each single value is a statistic over a large number of Fourier modes. The level of noise affecting each single observed visibility is much reduced in the process, but a lot of information is lost. The full 3D imaging of the signal will mostly require SKA capabilities (see Mellema et al. 2013 for expected capabilities) although some low-resolution imaging should be possible with LOFAR (Zaroubi et al. 2012).

Numerical simulations are an invaluable tool to predict and interpret the upcoming observations. Since power spectrum measurements are expected to come first, the 3D spherically averaged power spectrum is the first observable quantity that simulations have focussed on. In the emission regime (probably covering most of the EoR), predictions are in broad agreement with one another (e.g. Iliev et al. 2008; Lidz et al. 2008). In the absorption regime, early during the EoR, predictions are also available (Santos et al. 2008; Baek et al. 2010), but are very sensitive to unknowns, such as the relative amounts of X-ray, ionizing and Lyman band photons, which depend on the nature of the sources and their formation history. Other statistical quantities able to detect non-Gaussianities, such as the pixel distribution function or the skewness of the brightness temperature distribution, have been studied (Ciardi, Ferrara & White 2003; Mellema et al. 2006; Harker et al. 2009; Baek et al. 2010; Ichikawa et al. 2010). How to exploit the information in the full 3D tomography, however, is a subject that has barely been mined yet (Vonlanthen et al. 2011; Datta et al. 2012a; Majumdar, Bharadwaj & Choudhury 2012), but will likely receive a lot of attention in the upcoming years.

Somewhere between the low-noise low-information spherically averaged power spectrum and the high-noise high-information 3D imaging lies the angular power spectrum or, equivalently, the angular correlation function. Indeed, several factors induce an anisotropy in the power spectrum in such a way that properties are different along the line of sight (LOS) and perpendicular to the LOS. The first effect comes from peculiar velocity gradients along the LOS that enhance or dim the 21 cm brightness temperature, also called the Kaiser effect (Kaiser 1987). Barkana & Loeb (2005) showed how this anisotropy, cosmological in nature, could be separated in the linear regime from astrophysical isotropic contributions to the brightness temperature fluctuations such as ionization patchiness and spin temperature fluctuations. Since then, their method has been tested for robustness in the non-linear regime and further refined (Lidz et al. 2007; Mao et al. 2012; Jensen et al. 2013; Shapiro et al. 2013). There is a second source of anisotropy that may interfere with the Kaiser effect in extracting cosmological information. Since we will observe a light cone, points farther away will be seen earlier in the history of the EoR when bubbles were smaller, the IGM less heated by X-rays and the Wouthuysen–Field coupling weaker. This again introduces an anisotropy where the LOS is a special direction.

Barkana & Loeb (2006, hereafter BL) first attempted to quantify this anisotropy using a simple theoretical model. They found that the anisotropy occurs during the later part of the reionization history, when the neutral fraction x_n is smaller than 0.5, on scales larger than 50 comoving Mpc. They quantified it by comparing the two-point correlation functions along and perpendicular to the LOS. In their fig. 4, they found that the peaks of the correlation functions computed along and perpendicular to the LOS occur at different times. They also found that the anisotropy is stronger in the case of an earlier reionization history (Pop III stars), and that the amplitude of the peaks is about 100 per cent higher in the case of Pop III stars than for the case of Pop II stars. Their model makes a number of assumptions whose robustness is difficult to evaluate. One of the key ingredients is an estimate of the average radius of ionized bubbles as a function of redshift, based mainly on the average overdensity inside the bubble. This ignores the wide scatter in ionized region sizes that occurs especially during the overlap phase. Chardin, Aubert & Ocvirk (2012) found that the linear size of the largest ionized region at x_n = 0.5 is ∼20 times larger than the average ionized region size. Moreover, BL found that the anisotropy occurs in late reionization, during the overlap, when collective effects are most important and the ionization field is most non-Gaussian, making simple analytic estimations of correlation functions more problematic. Their predictions have not been cross-checked with full numerical simulations before, mainly because they require very large simulation boxes, and thus are costly. However, note the work by Datta et al. (2012b) that examines the light-cone effect on the isotropic power spectrum in a 163 cMpc box. In this work, we run full radiative transfer simulations in a 400 h⁻¹ cMpc box¹, including both stars and X-ray sources, and Lyman α radiation for a fully self-consistent modelling of the early EoR. As we will see, even though we do find anisotropies on similar scales, our quantitative predictions differ from BL's.

DESCRIPTION OF THE SIMULATION SUITE

Hydrodynamic run

We run the hydrodynamical simulation using gadget2 (Springel 2005), with 2 × 512³ particles in a 400 h⁻¹ cMpc box. Snapshots are extracted using a fixed interval of the scale factor Δa = 0.001. This produces more than 100 snapshots by z = 6, with a varying time interval of the order of 10 Myr between snapshots. The cosmological parameters are chosen from WMAP7+BAO+H₀ data: Ω_m = 0.272, Ω_Λ = 0.728, Ω_b = 0.0455, h = 0.704 (Komatsu et al. 2011).

The mass of a dark matter particle in our simulation suite is ∼4 × 10¹⁰ M_⊙, and that of a baryonic particle is ∼8 × 10⁹ M_⊙. Consequently, only haloes with masses above ∼5 × 10¹¹ M_⊙ are resolved. It is believed, however, that 80 per cent of the photons contributing to reionization are emitted by galaxies with masses smaller than 10⁹ M_⊙ (e.g. Choudhury & Ferrrara 2007). This implies that in our simulations, lacking the contribution of small galaxies, the luminosity of massive galaxies is boosted to complete reionization by z = 6. This relocation of ionizing photons from small to massive galaxies changes the topology of the reionization field to some degree. Whether it affects the level of anisotropy in the 21 cm signal is difficult to assess, but is indeed a possibility. A definitive answer to this question would require running much larger simulations (8192³ particles in the same volume), both for the dynamics and for the radiative transfer.

Ionizing radiative transfer run

We next compute the radiative transfer of UV ionizing photons using the code licorice (Baek et al. 2009; Iliev et al. 2009). In the version of licorice used in Baek et al. (2010), a finite velocity (c = 3 × 10⁸ m s^{− 1}) was used only for X-ray photons. UV photon packets were propagated at an infinite speed until 99.99 per cent of their content was absorbed. In the version we use for this work, we implement the actual finite speed of light for both X-ray and UV photons. Between two snapshots, photon packets typically travel 25 cMpc. During the overlap phase, ionized regions can extend to more than 100 cMpc, which are fully sampled in our 400 h⁻¹ cMpc simulation volume. Since the luminosity of the sources (updated with every new snapshot) changes substantially over a few tens of Myr, using an infinite speed of light results in overestimating the speed of the ionization fronts. This issue is most sensitive if the simulation volume is large enough to produce very large ionized patches, which is our case. Another issue occurs when dealing with large ionized patches. The estimated mean free path of ionizing photons through the forest of Lyman limit systems (LLS) at z ∼ 6 is in the 50 cMpc range (e.g. Songaila & Cowie 2010). Such systems are definitely not resolved in large-volume simulations of the EoR. Consequently, the flux received by the ionization fronts of very large patches is somewhat overestimated. Simulations taking LLS into account through sub-grid recipes do not suffer as much from using an infinite speed of light. We do not include any treatment of LLS.

The star formation efficiency and escape fraction (f_esc = 6 per cent) are calibrated to reach complete reionization at z = 6. Sources are formed in regions where the gas is above a fixed density threshold using a standard Schmidt law. We use a simple model to account for the production of X-rays by QSOs, supernovae and other possible sources: each source produces 0.6 per cent of its luminosity in the form of X-rays. We use a 100 per cent escape fraction for X-rays. This level of X-ray emission is equivalent to f_X ∼ 5 (as defined by Furlanetto et al. 2006, for example), considering that we use the same source model as in Baek et al. (2010). The level of X-ray emission quoted in Baek et al. (2010) should be corrected by f_esc, since this factor was (incorrectly) applied to X-ray emission.

Lyman band radiative transfer

In order to compute the Lyman band radiative transfer part of the simulation pipeline, we first interpolate the output of the ionizing radiative transfer run on a 512³ grid. We then emit 1.6 × 10⁹ photons between each pair of snapshots and propagate them at the speed of light. The full radiative transfer in the first five Lyman resonances is computed, assuming a flat spectrum in that range, between z = 13.8 and 7.5. A complete description of the numerical scheme and the physics included in our simulations can be found in Semelin, Combes & Baek (2007) and Vonlanthen et al. (2011). At lower redshifts, we make the assumption that the spin temperature is fully coupled to the gas kinetic temperature in the neutral regions: T_s = T_K.

21 cm differential brightness temperature data cubes

To investigate the correlation function and the evolution of the power spectrum, we analyse the results of the numerical simulation suite, which consists of 77 cubes of differential brightness temperature δT_b in the range of redshifts from z = 13.84 to 6.06. We compute these cubes at different resolutions for the following cases:

• basic δT_b – excluding velocity gradients (δT_b);

• as for basic δT_b, but taking into account the velocity gradient term (δT_b+velocity gradient);

• basic δT_b degraded using SKA-like noise and resolution (δT_b+SKA);

• as for basic δT_b, but including both velocity gradients and SKA-like noise and resolution (δT_b+velocity gradient+SKA);

• basic δT_b degraded using LOFAR-like noise and a resolution of 15 arcmin (δT_b+LOFAR 15 arcmin).

Including the velocity gradient contribution

Local peculiar velocity gradients modify the intensity of the emitted 21 cm signal (Bharadwaj & Ali 2004; Barkana & Loeb 2005). Since their contribution is anisotropic, we must include them to check that they can be disentangled from the light-cone effect. While the velocity gradient term in equation (1) is straightforward to implement, it is a linear approximation for small gradients: if the value of a negative gradient reaches that of the local Hubble flow, the brightness temperature (unphysically) diverges. Mao et al. (2012) have proposed an improved quasi-linear scheme to compute this contribution in a way that avoids divergences. In this work we use the usual linear approximation and implement a cut-off to avoid divergences.

Indeed, without a cut-off, the minimum and maximum values of δT_b are much below −1000 K and above 4000 K, respectively. These rare spurious cells strongly influence the shape of the power spectrum and correlation function, and significantly distort them.

To avoid this problem we use a clipping method, i.e. we replace cells with δT_b greater than 100 mK or smaller than −300 mK by cells with δT_b exactly equal to 100 mK and −300 mK, respectively. We estimate the robustness of this method on data containing the basic δT_b (without a velocity gradient). We check that the power spectrum and correlation function for the basic δT_b remains unchanged when data are clipped in this range. We use a rather low cut-off since we want to be able to check that we do not alter the properties of the signal itself by using a cut-off, independently of the spurious divergent values. We also check that a much higher cut-off on data with velocity gradient works just as well to remove the effect of spurious cells, and that the result is independent from the value of the cut-off, as long as we cut all cells with δT_b below ∼−10⁵ mK and above a few 10⁴ mK.

Including SKA noise and resolution

In order to include SKA-like noise and resolution, we proceed in three steps. We first build cubes of pure noise, using the noise power spectrum given in Santos et al. (2011). The resolution of these noise cubes corresponds to the expected SKA resolution at the corresponding redshift. We assume a 5 km core, a total observation time of 1000 h and a maximal baseline of 10 km. In order to calibrate the noise rms of our realizations, we consider the value of 1 mK at z = 6 mentioned in the SKA Design Reference Mission, with 1 MHz frequency resolution and a 1.1 arcmin angular resolution [although Mellema et al. (2013) estimate that this level should be reached only at 5 arcmin resolution]. Assuming Gaussian statistics, we normalize that value to the size of our simulation cells to obtain an rms of about 9 mK at z = 6. We then deduce the necessary sensitivity to reach this rms, and assume that the effective collecting area A_eff evolves as the square of the observed wavelength for observed frequencies above 100 MHz. For lower observed frequencies, A_eff is constant and equal to its value at 100 MHz.

Then we add the noise cubes to the full resolution δT_b cubes. We smooth the resulting data by applying a 3D Gaussian smoothing with the full width at half-maximum equal to the resolution at the corresponding redshift.

Including LOFAR noise and resolution

To simulate LOFAR noise, we go through the same three steps as for the SKA case. Here, we normalize the noise rms to 56 mK at 150 MHz, as given in Zaroubi et al. (2012) for an angular resolution of 3 arcmin and a frequency resolution of 1 MHz. Following Zaroubi et al. (2012), we assume a constant resolution of 15 arcmin at all redshifts. We scale the rms value at 15 arcmin resolution from the value at 3 arcmin by assuming that the noise is Gaussian, although this is not strictly the case. The behaviour of the rms with wavelength is illustrated in Zaroubi et al. (2012, fig. 3) based on a detailed model of the instrument. As a rough fit we scale the rms as λ², where λ is the observed wavelength. A detailed model for the evaluation of the rms can be found in Labropoulos et al. (2009). We add the noise cubes to the full resolution data and smooth the result to the expected LOFAR resolution (15 arcmin).

Building light cones

First, we make a flat sky approximation. This introduces a distortion of ∼4 cMpc along the LOS at the edge of a 3° field of view at z = 13.6. This is close to the SKA resolution, which is acceptable for the purposes of this paper. Secondly, although we will refer to a light cone, we actually build a light cylinder, assuming that all the LOS are parallel. On the scales where the anisotropy is found (100 cMpc), a real light cone would involve ∼1 per cent change in the linear size for two cross-sections 100 cMpc apart along the LOS. The effect would be negligible on the correlation function amplitude. One non-negligible effect would be to increase the size of the sample of pairs of points between z = 6 and 13.6 by up to a factor of 2, and thus possibly induce a slight change in the variance of the correlation function. In view of other approximations made in our method, we feel that the light cylinder approximation is acceptable.

When building a light cone from simulations, a common trick is to translate and rotate the box every time the light cone enters a new periodical replica of the simulation box. This avoids encountering the same structure repeatedly along the LOS and other artefacts of using a simulation box that is smaller than the length of the light cone. While this works well in the context of building galaxy catalogues, which are discrete samples, it cannot be applied to 21 cm light cones. Indeed, large ionized regions straddling the box boundary would be truncated without any possible justification. Consequently, we must forgo these rotations and translations. To avoid the artefacts, we must use a large simulation box. With our 400 h⁻¹ cMpc box, we need only three replications along the light cone.

Figure 1.

Open in new tabDownload slide

The light cone which models the Universe.

To improve statistical significance, we build separate light cones along the X, Y and Z axes of the simulation boxes. We also rotate the simulation boxes 45° around the X, Y or Z axis and build the light cones along the Y and Z, or Z and X, or X and Y axes, respectively. In this way, we generate nine quasi-independent axes of space–time. Along each of these nine axes, we build four light cones, starting from different positions in the first snapshot, namely from cell index 1, 128, 256 and 384. We use these four light cones generated along the same axis to check whether the periodicity conditions affect the outcome of our analysis. If periodicity affected the correlation functions, it would create characteristic features in the curves that would be simply differently shifted between the four cases. We do not observe any such effect. Altogether, we create 36 quasi-independent light cones.

The light cone at low redshift is quite smooth, while at high redshift a banding is visible due to the slice structure (Fig. 2a). Extracting the snapshots with shorter time intervals is possible but costly in terms of CPU and storage. To eliminate these discontinuities, we linearly interpolate between |$\delta T_{\mathrm{b}}(\boldsymbol {x},z_i)$| and |$\delta T_{\mathrm{b}}(\boldsymbol {x},z_{i+1})$| as a function of redshift, between consecutive snapshots i and i + 1 at redshifts z_i and z_i+1 (Fig. 2b).

Figure 2.

Open in new tabDownload slide

A fragment of the cross-section of the light cone (512×512 cells) along the space–time axis before (a) and after interpolation (b). False colours describe the differential brightness temperature δT_b [mK]. Units on the x- and y-axes are in cMpc.

The interpolation smoothes the structure at high redshift and does not distort it at low redshift. Thus, all our calculations are based on the interpolated light cones. Examples of the cross-section of different light cones are depicted in Fig. 10.

Correlations

• ξ is parametrized as a function of Θ, the angle between the LOS and the line connecting the two points. Θ = 0° means two points along the LOS and Θ = 90° means two points whose separation vector is perpendicular to the LOS,

• ξ is a function of the comoving distance r between two points at two different redshifts (Θ = 0°) or at the same redshift (Θ = 90°). We calculate ξ for three distances: r = 100, 50 and 20 cMpc,

• in the Θ = 0° case, a single redshift z is taken at the mid-point between the pair of points in terms of comoving distance,

• the mean |$\delta \bar{T}_{\mathrm{b}}$| at the given redshifts is subtracted from δT_b at each cell because this is what upcoming interferometer data will look like: they will only contain fluctuations of δT_b.

Since variance will prove to be an issue, the choice of binning will be critical when presenting our results for the correlation function. Our resolution element is a cubic cell 1.1 cMpc in size. Thus, naturally, our bin size for r is Δr = 1.1 cMpc. The bin size for Θ is then ΔΘ = 0| $_{.}^{\circ}$|63 for r = 100 cMpc, ΔΘ = 1| $_{.}^{\circ}$|26 for r = 50 cMpc and ΔΘ = 3| $_{.}^{\circ}$|15 for r = 20 cMpc. Using a single cell for the bin size in z would yield very noisy results. Instead, we use bins of 40 cells, which correspond to Δz = 0.29 at z = 13.8 and Δz = 0.1 at z = 6. These bins overlap since we use steps of five cells between them. A much larger Δz bin size would smooth the features in the evolution as a function of z and thus weaken the discriminating power of ξ. Using larger bin sizes for r and Θ would not probe information independent of that found using the large Δz bin size. Finally, we average over the 36 light cones. We could have used a binning system closer to observations, using constant bin size in frequency and angular resolution instead of redshift and comoving position. However, this would involve a lot of interpolation and would only introduce small differences in the error bars.

RESULTS

Global reionization history

We have calibrated the product of the escape fraction and the star formation efficiency to reach complete reionization at z = 6, in agreement with observations of the Gunn–Peterson effect in the spectra of high-redshift QSOs (e.g. Fan et al. 2006). The resulting global ionization history is shown in Fig. 3, along with the evolution of the average spin and kinetic temperatures. The corresponding Thomson scattering optical depth is τ = 0.054. This value is 2.5σ below the value inferred from observation of the CMB (τ = 0.089 ± 0.014; Planck Collaboration 2013). This is common for numerical simulations of the EoR that end reionization at z = 6: it is difficult to match both observational constraints simultaneously. Although a higher resolution allows starting reionization earlier, it does not seem to be enough to bring τ in agreement with observation while completing reionization at z = 6. This fact suggests that either we misinterpret the observations or that something is missing in the simulations (e.g. Ahn et al. 2012).

Figure 3.

Open in new tabDownload slide

Top: the box-averaged spin, kinetic and CMB temperatures as a function of the redshift. At redshifts lower than 7.5, we make the assumption that, in neutral regions, the spin temperature is fully coupled to the gas kinetic temperature: T_s = T_K. Bottom: cosmic neutral fraction x_n and the ionization fraction x_H as a function of redshift during the process of cosmic reionization.

BL considered two different ionization histories, early reionization driven by Population III (Pop III) stars and late reionization driven by Pop II stars. Their early reionization scenario, however, is rather unlikely since it completes reionization at z ∼ 13, and is presented mainly as a point of comparison. Our reionization history is similar to their late reionization scenario.

Evolution of the power spectrum

Fig. 4 shows the evolution of the differential brightness temperature power spectrum as a function of redshift for three different wavevectors corresponding to wavelengths of 20, 50 and 100 cMpc. We present these plots mainly to allow comparison with previous works. The evolution is in agreement at a large scale with model S7 in Baek et al. (2010), which also assumes strong X-ray heating (their fig. 6). Our result is also in general agreement with Santos et al. (2008). The magenta curve in Fig. 4 shows the evolution of the power spectrum computed from 400 h⁻¹ cMpc boxes cut from the light cone, for wavevectors which correspond to distances of 100 cMpc. The bumps are smoothed in this case because within the boxes used to compute the power spectrum, one finds different stages in the history of reionization, and when P(k) is computed for a given k value, it is averaged over a whole era of reionization history. Indeed, 400 h⁻¹ cMpc is more than the (comoving) distance travelled by light between redshifts 8 and 6, and between these two redshifts, in our simulation, the ionization fraction rises from 0.1 to 1. In the snapshot version, bumps have maxima at specific z. In the light-cone version, one averages over a non-negligible Δz, so the bumps are smoothed. We observe the smoothing effect at all scales studied. This is similar to what Datta et al. (2012b) found with a 163 cMpc box size, showing that the 3D isotropic P(k) is not a good diagnostic when applied to light-cone data on large scales. Consequently, and to facilitate comparison with BL's results, we will focus our analysis on the correlation function.

Figure 4.

Open in new tabDownload slide

Evolution of the differential brightness temperature power spectrum with redshift at wavevectors which correspond to distances of 20, 50 and 100 cMpc. Magenta (dotted) curve shows the power spectrum at 100 cMpc computed from the light cone.

The effect of sample variance

Compared with the 3D isotropic correlation function, the correlation functions computed transversely, and even more so along the LOS, represent much smaller sample sizes, and thus are more noisy.

While we could compensate for this sample-size noise by initializing the high-z end of the light cone at different slices in our simulation snapshot, this will not be an option with observational data. Consequently, we need to estimate how large a survey size is needed to mitigate sample variance. The values that are usually quoted for the isotropic density field are of no help in our case.

Fig. 5 (top panel) shows the fluctuations in the mean δT_b calculated for the light cones created along three perpendicular directions. Each value in these curves is an average computed in 40-cell-thick slices (∼30 h⁻¹ cMpc), perpendicular to the LOS. While this is a volume equivalent to a ∼ (170 h⁻¹ cMpc)³ cube, we still find a non-negligible variance. When we consider the evolution of the correlation perpendicular to the LOS at ∼ 100 cMpc (Fig. 5, bottom panel), the variance is much worse. Obviously, a 3°20′ × 3°20′ survey size is insufficient to mitigate the variance on such directional diagnostics.

Figure 5.

Open in new tabDownload slide

Top: the mean differential brightness temperature δT_b of the 21 cm signal from simulation snapshots of basic δT_b (black smooth dotted line). Colour curves denote the mean δT_b along the light cones built along the X-axis (solid red), Y-axis (dotted cyan) or Z-axis (short dotted blue). Bottom: example of the correlations functions at r = 100 cMpc between points perpendicular to the LOS (Θ = 90°). The calculations were made on the basis of the light cones formed along the X- (solid red), Y- (dotted cyan) or Z-axis (short dotted blue). The large differences reveal a large sample variance.

To quantitatively estimate the variance in a 3°20′ × 3°20′ survey size, we build, as stated above, 36 nearly independent light cones. We rotate the LOS by ±45° from each of the three axes of the grid, resulting in nine different LOS and use four different starting positions for each LOS, separated from one another by 100 h⁻¹ cMpc in our highest redshift data cube. We estimate the sample variance for a 3°20′ × 3°20′ survey size by computing the sample standard deviation for each point in the curve over the 36 values obtained from the 36 light cones. Finally, for a rough estimate of the error in a 20° × 20° survey size, i.e. in a solid angle 36 times greater, we assumethat this consists of 36 independent samples with Gaussianerror distributions in each (i.e. we neglect cross-correlations), and thus divide the sample standard deviation by |$\sqrt{36}$|⁠. It turns out that using the average over such a survey size allows us to differentiate the correlation along and perpendicular to the LOS at the level of ∼2σ at the redshifts when the anisotropy is greatest.

Anisotropy in the raw signal

Averaging over the 36 light cones, we compute the correlation function of the fluctuations in the differential brightness temperature at three different scales: 20, 50 and 100 cMpc. The results are presented in Fig. 6. In these plots, the error bars include only the contribution of sample variance.

Figure 6.

Open in new tabDownload slide

Correlation function of δT_b simulations at different redshifts expressed as the neutral fraction x_n in the range (0–0.9), (0.9–1) and (0.99–1) from left to right. Red (solid) and green (dotted) curves correspond to 0° and 90°, respectively. The curves are averaged over 36 different light cones, and error bars are estimates of the standard deviation of the mean. The top panels show the mean δT_b from the original 77 snapshots.

In the emission regime (left-hand panel, x_n < 0.9), we find correlations whose amplitudes are in general agreement with the theoretical predictions in BL at all three scales. We also agree with them in the point that the relative amplitude of parallel to perpendicular correlations shows non-negligible anisotropy only on scales of ∼100 cMpc. However, BL found that perpendicular and parallel correlations peak at different stages of the reionization history, with the correlation along the LOS peaking first. We find a more complicated pattern, where the correlation along the LOS rises up to more than 2σ above the correlation perpendicular to the LOS, in two to three different peaks. The statistical significance of this behaviour is not strong enough that we should feel entitled to attempt a detailed interpretation of those two to three peaks. Only the general dominance of the correlation along the LOS seems robust at this point. The behaviour found by BL disagrees with our results by more than 3σ. However, the model of BL relies on bubbles of fixed sizes (evolving with redshift), without accounting for the total effect occurring during the overlap phase. This is exactly when they found an anisotropy. BL's model uses the semi-analytical method of Furlanetto, Zaldarriaga & Hernquist (2004) to compute the ionization field, and Santos et al. (2008) found that this method shows the largest difference with full simulations around x_n ∼ 0.2, again where the anisotropy is found. Consequently, we do not expect to find detailed quantitative agreement with BL.

BL's model is restricted to the emission regime when the spin temperature is assumed to be fully coupled to the kinetic temperature of the gas, itself much larger than the CMB temperature. Our simulation suite includes 3D radiative transfer for both the inhomogeneous heating by X-rays and the inhomogeneous coupling through Lyman α photons. So we can explore the light-cone anisotropies in the absorption regime during the early EoR. In the middle panel of Fig. 6, the 100 cMpc correlation along the LOS rises about 2σ in terms of sample variance above the correlations perpendicular to the LOS. This happens when the Universe is ∼2 per cent ionized, as the average differential brightness temperature rises again after the minimum in δT_b, i.e. the maximum in absorption. While the anisotropy in the emission regime, at x_n ∼ 0.2–0.5, was created by the time evolution of the ionization field, it is likely that the anisotropy at x_n ∼ 0.98 is created by the time evolution of the kinetic temperature fluctuations. Fig. 7 qualitatively confirms this interpretation. It shows a crosscut of the light cone for |$\delta T_{\mathrm{b}} - \delta \bar{T}_{\mathrm{b}}$| and T_K. Structures at x_n ∼ 0.2–0.5 and x_n ∼ 0.98 seem to be more extended along the LOS: growing bubbles of ionization ( x_n ∼ 0.2–0.5) and heated IGM regions (x_n ∼ 0.98) overlap preferentially along the time axis. In the right-hand panel of Fig. 6, anisotropy is again visible at the 2σ–3σ level in terms of sample variance at 100 cMpc scales, but also, to a lesser degree, in 50 cMpc correlations. It occurs in the x_n = 0.994–0.999 range and is probably created by the time evolution of coupling fluctuations: growing regions of coupled IGM. As we can see in Fig. 8, where the isocontours of the coupling coefficient x_α in a region around one of the first isolated sources are shown, the coupling fluctuations are clearly anisotropic, revealing a parabolic envelope on large scales. The axis of the paraboloid is the axis of the light cone. The analytic formula for this parabola is very simple to establish if one assumes a source of constant luminosity switching on at a given time, and that all the LOS in the region are parallel to each other (small angular size).

Figure 7.

Open in new tabDownload slide

Examples of the same cross-section of the light cone along a space–time axis from the basic δT_b case. Colours in the left-hand panel depict the difference between δT_b at given point and mean δT_b at a given redshift: |$T_{\mathrm{b}} = \delta T_{\mathrm{b}} - \bar{\delta T_{\mathrm{b}}}$| in mK; the right-hand panel shows kinetic temperature in K.

Figure 8.

Open in new tabDownload slide

Fragment of a cross-section of the light cone showing the first source. The colours depict the x_α factor on the logarithmic scale. Values below 10⁻² are set to 10⁻². The range of the x- and y-axes (in cMpc) corresponds to a linear size of 220 cMpc. Values of cells on the y-axis correspond to a range of neutral fraction x_n from 0.999 84 to 0.999 73.

Effect of peculiar velocities

Peculiar velocity gradients are a source of anisotropic fluctuations of δT_b (Barkana & Loeb 2005). Thus, we have to determine whether, on 100 cMpc scales, they can be disentangled from the light-cone anisotropy. In Fig. 9, we compare the correlation functions of δT_b with and without the contribution from peculiar velocity gradients, computed from a single light cone so as not to fold in variance effects. The velocity gradient term moderately increases the amplitude of the correlation but does not significantly change the overall shape of the curves. On the 100 cMpc scale, the net effect is generally smaller than the variance for a 20° × 20° light cone. We conclude that on the 100 cMpc scale, when light-cone anisotropy peaks, it dominates over peculiar velocity anisotropy. However, we do apply an approximate treatment of the effect of peculiar velocities, which is valid in the linear regime. A definitive answer would require a treatment such as that advocated by Mao et al. (2012) or Jensen et al. (2013).

Figure 9.

Open in new tabDownload slide

Example of the correlation functions obtained from a single light cone. Correlations along (0°) and perpendicular (90°) to the LOS are presented for two cases: basic δT_b and δT_b+velocity gradient. The left- and right-hand panels show different ranges for the neutral fraction x_n. Peculiar velocity increases the amplitude of the correlation by a moderate amount but does not significantly change the overall shape of the curves.

Anisotropies observed with LOFAR and SKA

Will it be possible to detect the light-cone anisotropy in LOFAR and SKA observations? We make the optimistic assumption that foregrounds have been subtracted with negligible residuals. We model only two instrumental effects: the noise resulting from limited sensitivity and integration time, and the limited resolution (see Section 2.4). Fig. 10 shows how the cosmological δT_b light cone is modified when observed with SKA at resolution ∼1.5 arcmin (depending on redshift) and with LOFAR at resolution 15 arcmin.

Figure 10.

Open in new tabDownload slide

Example of cross-sections of the light cones: basic δT_b (left), δT_b with SKA-like noise simulation (middle) and δT_b with LOFAR-like noise (right). False colours show the differential brightness temperature δT_b[mK]. For better visualization, 5 per cent of the brightest and 5 per cent of the darkest pixels are not shown. The y-axes in three pictures are the space–time axes and are denoted in three mutually dependent units: the neutral fraction of hydrogen x_n, redshift and cells.

Fig. 11 shows how the estimates of the correlation along and perpendicular to the LOS are affected. In this figure, the error bars combine contributions from sample variance and instrumental noise. We could have estimated the latter by generating many realizations of the noise and computing the rms for each point. This would have been costly in CPU time. Instead, we used an analytical estimate described in Appendix A. It is obvious that the instrumental noise from the SKA will not be an issue: it is small compared to the level of sample variance that can be expected for a 20° × 20° survey (36 times the size of our light cones). The size of the survey will be the critical parameter since the minimal value for the survey size considered in the current design is 5° × 5°. For the latter value, the anisotropic features of the correlation function may not be detected with statistical significance over the sample variance. In the case of LOFAR, the contribution of instrumental noise is similar to the sample variance for a 20° × 20° light cone. As a result, none of the anisotropy features are clearly detected. Reducing sample variance with a very large survey area might only result in a 2σ detection. Moreover, limited resolution simulations tend to delay the beginning of reionization and the absorption regime anisotropies may well occur at higher redshifts than in our simulations. Consequently, the anisotropies in the absorption regime will probably not fall in the LOFAR band anyway.

Figure 11.

Open in new tabDownload slide

Two-point correlation functions calculated at a comoving separation of 100 cMpc along (red solid) and perpendicular (green dotted) to the LOS for three different sets of data (from top to bottom): basic δT_b, δT_b+SKA and δT_b+LOFAR 15 arcmin. Left-hand panels show the range of the neutral fraction 0 ≤ x_n ≤ 0.9, and right-hand panels 0.9 ≤ x_n ≤ 1. Basic δT_b and δT_b+SKA are almost identical. Error bars for the LOFAR case come from sample variance and instrumental noise.

CONCLUSIONS

To evaluate the amplitude of light-cone anisotropies in the 21 cm signal from the EoR, we run a suite of simulations in a 400 h⁻¹ cMpc box with 2 × 512³ particles. Our simulations include dynamics (gravitation and hydrodynamics), source formation recipes, ionizing UV radiative transfer, X-ray radiative transfer and Lyman α radiative transfer: this allows us to model the signal both in the absorption and in the emission regime.

To quantify the anisotropy, we examine the two-point correlation function of the differential brightness temperature of the 21 cm signal perpendicular and parallel to the LOS, on scales of 20, 50 and 100 cMpc. Our first conclusion is that sample variance affects these correlation functions over much larger volumes than is usual for cosmological diagnostics. We find that a light cone >10° across (∼1.5 cGpc at z ∼ 10) is necessary to reach statistically significant results. This is due to the directional nature of our diagnostics. We are able to mitigate variance with our 400 h⁻¹ cMpc simulations by building several semi-independent light cones.

Using a sufficiently large sample size, we are able to find an anisotropy in the emission regime, but only on a 100 cMpc scale. This matches the theoretical prediction by BL. While the amplitudes of our correlation functions are also in broad agreement with their predictions, we do not agree on the detailed behaviour. We find that the parallel correlation generally dominates over the perpendicular correlation at neutral fractions in the range x_n = 0.2–0.5. We interpret this as ionized bubbles connecting preferentially along the LOS during the overlap phase.

Our analysis also includes the absorption regime. At the early stages of the EoR, the correlation function parallel to the LOS rises above the perpendicular correlation at x_n ≈ 0.98 and the opposite is true in the x_n ≈ 0.994–0.999 range (the exact values depend on the nature of the sources). The feature at x_n ≈ 0.98 corresponds to the onset of heated regions. We interpret the anisotropic feature in the x_n ≈ 0.994–0.999 range as the onset of coupled regions. These regions are elongated along the light cone.

We checked that, on 100 cMpc scales, the inclusion of the velocity gradient contribution to δT_b does not change the shape of the correlation functions much. The velocity gradient contribution only moderately enhances the amplitude, well below the amplitude difference of the perpendicular and parallel correlation when the light-cone anisotropy peaks. This is important since the velocity gradient contribution is also anisotropic and could confuse the light-cone anisotropy signal.

We assess the observability of these anisotropic features by SKA and LOFAR using standard survey parameters. We assume that the residuals from foreground subtraction are much smaller than the signal. We find that SKA can detect all the anisotropy peaks we have identified, provided that the size of the survey is larger than 10° × 10°. In the case of LOFAR, the emission regime anisotropy peak could only be marginally detected in a 15 arcmin resolution survey with a very large size.

There are several reasons why it will be interesting to search for light-cone anisotropy in the observations. First, the anisotropy peaks are connected to milestones in the global reionization history: the onset of coupling and heating, and the overlap of ionized bubbles. The global history can only be indirectly retraced by interferometers that do not give any information about the average signal. Secondly, the light-cone anisotropy is characteristic of the signal and should not show the same features in the residual of foreground subtraction. This could provide a test to check whether foregrounds have been correctly removed.

The authors wish to thank Garrelt Mellema and Boudewijn Roukema for a thorough critical reading of the manuscript and the anonymous referee for a useful report. Part of this work consists of research conducted in the scope of the HECOLS International Associated Laboratory. Some of this work was carried out in the context of the LIDAU project. The LIDAU project is financed by the ANR (Agence Nationale de la Recherche, France) grant ANR-09-BLAN-0030.

REFERENCES

APPENDIX A: CONTRIBUTION OF INSTRUMENTAL NOISE TO ERROR BARS OF THE TWO-POINT CORRELATION FUNCTION

Let us call δT_b(r_i, z) the differential brightness temperature of the 21 cm signal of a pixel in the light cone at redshift z, and consider position r_i in the plane perpendicular to the LOS. The index i fully covers the two-dimensional transverse maps. Then, let us write n(r_i, z) for the instrumental noise in the same pixel. In a given transverse map, n is a realization of a Gaussian random field with standard deviation σ_n(z).

Since |$n({r_j^1},z)$| and |$n({r_j^2},z)$| are independent, we can write |${\mathrm{Dev}}[{1 \over N} \sum \limits _{j=1}^N n({{r_j^1}}) n({{r_j^2}}) - \langle n\rangle ^2] = 2{ \sigma ^2_n(z) \over M^2} {.}$|

Finally, we get the following expression for the standard variation of the correlation function (error bar due to noise): |${\mathrm{Dev}}( \Delta \xi )= 2 { \sigma _n(z) \sigma _{T_{\mathrm{b}}}(z)\over M^2} + 2 { \sigma ^2_n(z) \over M^2}$| In observations it would not be possible to evaluate |$\sigma _{T_{\mathrm{b}}}(z)$|⁠. But in our case we can compute it from the simulated signal. The relative contribution of each term is plotted in Fig. A1.

Figure A1.

Open in new tabDownload slide

Contribution of instrumental noise to error bars of the two-point correlation function.