The Impact of Neutral Intergalactic Gas on Lyα Intensity Mapping during Reionization

When and how the Epoch of Reionization (EoR) occurred is one of the biggest unsolved problems in cosmology (for a recent review see McQuinn 2016). A well-established observable of this epoch is damping wing absorption in the spectrum of galaxies and quasars, which is sensitive to neutral intergalactic gas that lies tens of megaparsecs in the foreground (Miralda-Escude 1998). Such absorption is the leading explanation for the rapid decline in the number of Lyα-emitting galaxies above z ≈ 6.6 (Iye et al. 2006; Ono et al. 2012; Schenker et al. 2012). Indeed, many models find that this requires fast evolution in the neutral fraction over this redshift interval (Choudhury et al. 2015; Mesinger et al. 2015).

There should also be a diffuse H i Lyα component that owes to re-emissions following intergalactic medium (IGM) absorptions; the feasibility of mapping this diffuse intensity is our focus. Previous work suggests that detecting the total diffuse Lyα flux is possible with the intensity mapping satellites SPHEREx and the Cosmic Dawn Intensity Mapper (CDIM; Silva et al. 2013; Doré et al. 2014; Pullen et al. 2014; Comaschi & Ferrara 2016; Cooray et al. 2016; Heneka et al. 2017). However, previous work has generally only considered the direct diffuse emission from the sum total of all galaxies, ignoring scattering from a neutral IGM that should imprint the structure of reionization on this signal.

This Letter presents calculations of the diffuse Lyα intensity field from the EoR. In Sections 2 and 3 we describe these calculations. Section 4 presents simulated Lyα intensity maps. Section 5 concludes with a brief discussion of our results and the major takeaways. Throughout, we assume a ΛCDM cosmology with parameters consistent with those reported by the Planck Collaboration (2014): Ω_m = 0.32, Ω_Λ = 0.68, Ω_b = 0.049, h = 0.67, σ₈ = 0.83, and n_s = 0.96. All cosmological distances are given in comoving units.

We use the 21cmFAST code to generate 200 Mpc cubic realizations of reionization at z = 7 with a resolution of 256³ voxels, outputting the IGM density, the IGM neutral fraction, and the halo mass density field (Mesinger et al. 2011). The density and velocity fields in 21cmFAST are computed with the Zeldovich approximation; the halo positions and ionization fields are calculated using excursion set models. Excursion set algorithms for reionization successfully reproduce the morphology of reionization in full radiative transfer simulations (Zahn et al. 2007). There are three main astrophysical parameters that determine the properties of reionization in 21cmFAST: the minimum virial temperature of halos hosting ionizing sources, T_min, the maximum radius of individual H ii regions, R_min, and the ionizing efficiency, ζ, which encodes quantities such as the star formation efficiency and the escape fraction of ionizing photons (for details see Mesinger et al. 2011). For all runs, we set T_min = 5 × 10⁴ K (corresponding to a minimum halo mass of M_halo,min = 1.5 × 10⁹ M_⊙) and R_max = 15 Mpc. While sources down to M_halo,min ≈ 10⁸ M_⊙ can host substantial star formation, feedback likely makes these halos less efficient (Furlanetto et al. 2017); our choice of M_halo,min accounts for this inefficiency within the parametrization of 21cmFAST. Furthermore, R_max is set by the mean free path of ionizing photons, which current models suggest is less than 15 Mpc at z = 7 (Worseck et al. 2014; D'Aloisio et al. 2018).³ Smaller R_max or M_halo,min would increase the sizes of the effects that we investigate. With these parameters fixed, we simulate three different ionizing efficiencies, ζ = 10, 15, and 20, in order to compare three different ionized fractions.

21cmFAST outputs a neutral fraction field that is either equal to 0 or 1. This binarity does not account for the small residual hydrogen within ionized regions, which can scatter Lyα photons that redshift through line center. We set the neutral fraction inside ionized regions to x_{H i} = 10⁻⁴ by hand, with the value chosen to yield a significant line center optical depth (mimicking the saturation seen in the highest redshift Lyα forest). However, varying this value a factor of a few has no impact on our results. We do not use 21cmFAST to compute temperatures of the gas and instead for simplicity assume that neutral regions have a uniform temperature of T = 10³ K and that ionized regions have a temperature of T = 10⁴ K. These temperatures are what could be expected due to X-ray heating in the neutral IGM (e.g., Pritchard & Furlanetto 2007) and heating associated with reionization in the ionized IGM. We find that setting the entire box to T = 10³ K does not change the Lyα power spectrum by more that 10 percent on the spatial scales examined in this Letter.

We assume that Lyα photons are produced by hydrogen recombinations in the interstellar medium (ISM) of galaxies in halos larger than T_min. The total Lyα luminosity of a galaxy is related to the star formation rate (SFR), , by

where f_esc is the fraction of ionizing photons that escape into the IGM.⁴ This equation assumes no dust absorption so that every ionization results in 0.6 Lyα photons, and a Salpeter initial mass function (IMF; Schaerer 2003) over a mass range of 1–100 M_⊙ with metallicity Z = 0.04; other empirically motivated PopII IMFs result in factor of ∼2 differences. Our model assumes SFR proportional to halo mass for all halos with T > T_min, where the global star formation rate density is normalized to ρ_⋆ = 0.015 M_⊙ yr⁻¹ Mpc⁻³. This value is similar to ρ_⋆ = 0.02 M_⊙ yr⁻¹ Mpc⁻³ measured at z ≈ 6.8 in Bouwens et al. (2015). Our predicted intensities scale linearly with this overall normalization. Because we always assume the same SFR prescription, modifying ζ to change the mean neutral fraction is equivalent to changing f_esc. However, we use f_esc = 0.1 in Equation (1) to fix the total luminosity (instead of allowing it to change by 5% across our simulations), permitting a more straightforward comparison. Unless otherwise noted, we assume that the emergent Lyα spectrum escaping from the ISM of each star-forming galaxy is a Gaussian with a 1 − σ width of 100 km s⁻¹. Lastly, to simplify the calculation, we take 10% of halos to be active Lyα emitters at a given time, making active emitters 10 times more luminous to preserve the global emissivity; we comment when this assumption affects our results.

To compute the impact of Lyα scattering on simulated intensity maps, we have developed a parallelized 3D Monte Carlo Lyα radiative transfer code, similar to the one described in Faucher-Giguère et al. (2010; see also Zheng & Miralda-Escudé 2002; Cantalupo et al. 2005; Dijkstra et al. 2006; Laursen & Sommer-Larsen 2007). For a detailed description see Appendix C in Faucher-Giguère et al. (2010). To generate a distant-observer image, our calculations follow the paths of Lyα photon packets as they scatter throughout the simulated neutral gas field. We follow these packets for spatial steps that are much smaller than our simulation resolution, with specifications chosen for convergence. We also utilize the "accelerated scheme" described in Appendix C of Faucher-Giguère et al. (2010) with x_crit = 5.0, which greatly speeds up our calculation without significantly impacting the results. We performed a series of tests on our code including the static sphere, redistribution function, and inflowing/outflowing spheres, described in Appendix C of Faucher-Giguère et al. (2010), and find good agreement with the established solutions. With regard to the calculations in this Letter, convergence is achieved following 3 × 10⁶ photon packets and for the spatial resolution of our 21cmFAST realizations.

Here we present our simulated Lyα intensity maps. These maps all have the same underlying dark matter density distribution, but different reionization histories (created by varying the ionizing efficiency parameter in 21cmFAST). In Figure 1, we plot projections of the neutral fraction and Lyα intensity from the z = 7 snapshot of two of our simulations that have mean neutral fractions of (top row) and 0.36 (bottom row). We separate the Lyα intensity into components with frequencies that are close (x ≤ 4) and far (x > 4) from line center at the final scatter, where and ν_α is the frequency of Lyα. The component close to line center includes photons that do not scatter at all or that scatter while redshifting through line center within ionized regions. The component far from line center comes mainly from photons that scatter off of highly neutral gas in the IGM and need to redshift into the wing of the Lyα line to escape. The map (top row) has a significantly larger neutral IGM scattered component than the . The enhanced scattering in the former occurs because the ionized bubbles in the IGM are smaller, leading to Lyα photons being less redshifted by the time they reach the edge of the bubbles. Note that for , 0.63, and 0.36, a fraction 0.94, 0.86, and 0.76 of the Lyα photons scatter at least one time, respectively. For the case, a fraction of 0.34 (0.87) photons scatter in the highly neutral IGM with x > 4.

Zoom In Zoom Out Reset image size

Figure 2 shows the distribution of distances from the source at which the observed photon last scatters, . For the case, has FWHM that spans a decade and peaks at r ≈ 5 Mpc, slightly smaller than the ≈8 Mpc ionized bubble in a neutral IGM for which a Lyα photon emitted at line (and bubble) center has optical depth unity. (Another useful limit to note is that of a fully neutral IGM around sources, which yields the ∼1 Mpc in extent; Loeb & Rybicki 1999, Lyα "halos.") For the case of , there is again a peak at r ≈ 5 Mpc if the residual neutral fraction is not included inside of the ionized bubbles. However, when this residual H i is (realistically) included, there is an additional peak at r ≈ 1 Mpc that obscures the r ≈ 5 Mpc peak, which primarily arises from photons emitted on the blue side and that last scatter as they redshift through line center within ionized regions.

Zoom In Zoom Out Reset image size

The IGM scattering damps the small-scale fluctuations in the maps. Figure 3 illustrates this damping in the intensity power spectrum. The model curves' normalizations scale quadratically with the global star formation rate density, and so differences in the overall normalization likely cannot be used to detect reionization. However, differences in the slope of P_Lyα on the linear scales that are shown likely cannot be mimicked by uncertainties in the galaxy formation model. For our fiducial case between k = 0.1–1 Mpc⁻¹, the power spectrum scales roughly as for , for , for , and for . The most significant slope differences require relatively high neutral fractions, which likely occur at higher redshifts than z ≈ 7.

Zoom In Zoom Out Reset image size

One of the main challenges for detecting this signal is contamination due to line emission from foreground galaxies (Gong et al. 2014; Comaschi et al. 2016). The most problematic for z ∼ 7 Lyα is z ∼ 1 Hα, which can dwarf the Lyα signal. One way to mitigate this contamination is to cross correlate with a different probe of the same underlying density field. Here we explore the cross correlation of Lyα intensity maps and a future galaxy redshift survey. We compute the Lyα-galaxy cross-power spectrum and the corresponding sensitivity for two hypothetical space-based intensity mapping telescopes, one roughly modeled after CDIM and the other after SPHEREx. Unfortunately, a suitable high-redshift galaxy sample does not yet exist. We assume that the galaxies in our survey have a volume density of n_gal = 10⁻⁴ Mpc⁻³ and a mean bias of b_gal = 3 (similar to the Lyα emitters measured with Subaru (Ouchi et al. 2018), although better redshift determinations would be needed for this population). The dimensions of our survey are assumed to be 3° × 3° on the sky and Δz = 0.8 in redshift along the line of sight. The signal-to-noise of the cross-power spectrum scales as , where V is the survey volume.

In Figure 4, we plot the cross-power spectrum and associated 1 − σ sensitivity for our hypothetical surveys. The cross-power spectrum is estimated by cross-correlating our Lyα intensity maps with the density field from 21cmFAST multiplied by b_gal. The error on the power spectrum is given by , where P_N is the noise power from our hypothetical Lyα intensity-mapping telescopes, P_Hα is the power from interloping Hα sources, and N_k is the number of independent modes per k-bin. The noise power is calculated using Equations (16) and (17) in Comaschi et al. (2016) assuming a telescope diameter of 1.5 m for the CDIM-like instrument and 0.2 m for the SPHEREx-like instrument, an integration time of 10⁵ s, a zodiacal light background intensity of νI_ν = 500 nW m⁻² sr⁻¹, and an observational efficiency of detecting a photon accounting for losses in the instrument of  = 0.5. We assume spectral resolutions of R = 300 and R = 40 for our CDIM-like and SPHEREx-like instruments, respectively.

Zoom In Zoom Out Reset image size

The power from foreground Hα remaining after masking out the brightest sources is taken to be P_Hα = 0.04 ×(Mpc⁻¹/k) nW m⁻² sr⁻¹ Mpc³. This is roughly the same as computed by Pullen et al. (2014) for z = 6 and a limiting flux above which Hα sources can be removed of 10⁻¹⁷ erg s⁻¹ cm⁻² (see their Figure 13). This flux cut corresponds to an r-band AB magnitude of m_r ≈ 26.5, which will be observable over large areas with telescopes such as the Hyper Suprime-Cam (Pullen et al. 2014). We note that the errors in cross power are dominated by the Hα interlopers on large scales for the SPHEREx-like instrument, and on all scales shown for the CDIM-like instrument. Reducing the size of the CDIM-like instrument to 0.5 m has very little impact on the results (for a survey that is this small). We also note that if the flux limit to remove Hα sources were an order of magnitude higher (corresponding to m_r ≈ 24), the unmasked P_Hα would increase by roughly a factor of 30 (Pullen et al. 2014). For the CDIM-like instrument this increases the errors on the cross-power spectrum by ∼5.5 for all wavenumbers shown in Figure 4 (as this is the dominant source of error). For the SPHEREx-like instrument, this increases the error by a factor of ∼5 at k ∼ 0.1 Mpc⁻¹ and by a factor of ∼2 at k ∼ 1.0 Mpc⁻¹. This pessimistic masking criterion would make it challenging to observe the shape of the power spectrum with either instrument. As in the case of the Lyα power spectrum, the shape of the Lyα-galaxy cross-power spectrum also depends sensitively on the mean neutral fractions of our reionization simulations. At k = 1 Mpc⁻¹, this results in the cross-power spectrum being 1.4, 2.4, and 6.7 times lower for , 0.63, and 0.86, respectively, than for .

We have investigated two potentially smoking gun signatures for diffuse IGM emission: an imaginary cross power, and non-isotropy in k. There is ∼10 Mpc (comoving) spatial offset between observed galaxies and Lyα intensity maps along the line of sight due to Lyα photons needing to be redshifted approximately ∼1000 km s⁻¹ in order to escape to the observer. This leads to a non-zero imaginary component of the cross-power spectrum, , which is shown along with the real component in Figure 4. However, its low amplitude will make it more challenging to observe. Additionally, because large ionized structures oriented along the line of sight enhance transmission, modes oriented along the line of sight (which lead to smaller ionized structures in this direction) have suppressed power relative to modes perpendicular to the line of sight (which lead to larger ones). We find that the anisotropy in P_Lyα(k) is only a large effect for high mean neutral fraction, and leave a detailed analysis for future work.

Our choice of duty cycle, _duty = 0.1, has no impact on the cross-power spectrum in our calculation (but it does impact high-k in the auto, increasing P_Lyα(k) by a factor of ∼2 at k ∼ 1 Mpc⁻¹). We note that in reality there is likely to be some cross-shot power between the observed galaxies and the Lyα intensity map that will enhance the signal at high wavenumbers (with the neglected term scaling as ∼k³ in ).

We have simulated the Lyα intensity field from the EoR, including for the first time scattering by neutral intergalactic gas. We find that the scattering of Lyα photons in the IGM produces a scale-dependent suppression of spatial fluctuations on linear scales that likely cannot be mimicked by other processes, with the magnitude of this suppression dependent on the mean IGM neutral fraction. We find that using a Lyα intensity map supplied by an instrument similar to SPHEREx and especially CDIM, this signal at z ∼ 7 could be detected in cross correlation with a future galaxy redshift survey. In future work, we intend to investigate the dependencies of this signal on different reionization scenarios, as well as the effect of additional processes that could enhance the EoR signal (such as Lyα cooling in the ionizing fronts of expanding H ii regions and from galaxy continuum photons being redshifted into Lyman-series lines).

The Flatiron Institute is supported by the Simons Foundation. The simulations were carried out on the Flatiron Institute supercomputer, Rusty. M.M. was supported by NASA ATP award NNX17AH68G and by the Alfred P. Sloan Foundation.