A Halo Model Approach to the $\rm{21\, cm}$ and $\rm{Ly}\alpha$ Cross-correlation

TL;DR

This paper develops a halo-model approach to predict the cross-correlation between 21 cm HI fluctuations and Lyα emitters during the epoch of reionization, accounting for foreground contamination and forecasting detectability.

Contribution

It introduces a novel halo-model-based method to calculate the 21 cm-Lyα cross-correlation, including foreground mitigation and parameter forecasting.

Findings

Cross-correlation transitions from negative to positive during reionization.

Masking high flux Lyα sources reduces foreground contamination.

Forecasts indicate potential detectability with upcoming experiments.

Abstract

We present a halo-model-based approach to calculate the cross-correlation between $\rm{21\,cm}$ HI intensity fluctuations and $\rm{Ly}\alpha$ emitters (LAE) during the epoch of reionization (EoR). Ionizing radiation around dark matter halos are modeled as bubbles with the size and growth determined based on the reionization photon production, among other physical parameters. The cross-correlation shows a clear negative-to-positive transition, associated with transition from ionized to neutral hydrogen in the intergalactic medium during EoR. The cross-correlation is subject to several foreground contaminants, including foreground radio point sources important for $\rm{21\,cm}$ experiments and low-$z$ interloper emission lines, such as $\rm{H}\alpha$, OIII, and OII, for $\rm{Ly}\alpha$ experiments. Our calculations show that by masking out high fluxes in the $\rm{Ly}\alpha$ measurement, the correlated foreground contamination on the $\rm{21\,cm}$-$\rm{Ly}\alpha$ cross-correlation can be dramatically reduced. We forecast the detectability of $\rm{21\,cm}$-$\rm{Ly}\alpha$ cross-correlation at different redshifts and adopt a Fisher matrix approach to estimate uncertainties on the key EoR parameters that have not been well constrained by other observations of reionization. This halo-model-based approach enables us to explore the EoR parameter space rapidly for different $\rm{21\,cm}$ and $\rm{Ly}\alpha$ experiments.

Tables

Table 1. Table 1: 21 ​ cm 21 cm \rm{21\,cm} Experiment

Exp.	${\lambda }^{\mathrm{rest}}\left(m\right)$	$D_{\max }\left(m\right)$	$T_{\mathrm{sys}}\left(K\right)$	$t_0\left(\mathrm{hr}\right)$	$A_{\mathrm{tot}}\left(m^2\right)$	FOV (${\mathrm{deg}}^2$)	BW	$\delta \nu$	$N_a$	$z$
SKA1-LOW	0.21	1000	400	1000	925	13	18MHz	3.9kHz	433	8

Table 2. Table 2: Ly ​ α Ly 𝛼 \rm{Ly}\alpha Experiment

Exp.	${\lambda }^{\mathrm{rest}}$ (Å)	${\Omega }_{\mathrm{pix}}{(}^{\prime \prime })$	$A_s\left({\mathrm{deg}}^2\right)$	${\sigma }_{\mathrm{pix}}\left(\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}{\mathrm{sr}}^{-1}\right)$	$\lambda /\delta \lambda$	BW ($\mu$m)
CDIM	$1216$	1	300	$3\times {10}^{-6}$	500	0.7-8.0

Equations

$$\bar{x}_{e}(z)=\frac{1}{2}\Big{[}1+\tanh\Big{(}\frac{(1+z_{\rm{re}})^{3/2}-(1+z)^{3/2}}{\Delta y}\Big{)}\Big{]},$$

$$\tau ={\int }_0^{{\chi }_{\mathrm{re}}}d\ell n_e({\chi }^{\prime }){\sigma }_T,$$

$$T_{0}(z)=27\Big{(}\frac{1-Y_{\rm p}}{1-0.248}\Big{)}\Big{(}\frac{\Omega_{b}}{0.044}\Big{)}\Big{[}\Big{(}\frac{0.27}{\Omega_{m}}\Big{)}\Big{(}\frac{1+z}{10}\Big{)}\Big{]}^{1/2}\,\,(\rm{mK}),$$

$$\psi (\text{x},z)={\bar{x}}_H(1+{\delta }_x)(1+\delta )={\bar{x}}_H(1+{\delta }_x+\delta +{\delta }_x\delta ).$$

$${\xi }_{\psi \psi }={\xi }_{xx}(1+{\xi }_{\delta \delta })+{\bar{x}}_H^2{\xi }_{\delta \delta }+{\xi }_{x\delta }(2{\bar{x}}_H+{\xi }_{x\delta })+{\xi }_{x\delta x\delta }.$$

$$P_{\psi \psi }^{(2h)}$$

$$\frac{T_{\rm{vir}}}{10^{4}\rm{K}}=1.1\Big{(}\frac{\Omega_{m}h^{2}}{0.15}\Big{)}^{1/3}\Big{(}\frac{1+z}{10}\Big{)}\Big{(}\frac{M_{\rm{th}}}{10^{8}M_{\odot}}\Big{)}^{2/3}.$$

$${\bar{V}}_b=\int dRP(R)V_b(R)=\frac{4\pi {\bar{R}}^3}{3}{e}^{9{\sigma }_{\ln R}^2/2},$$

$$P(R)=\frac{1}{R}\frac{1}{\sqrt{2\pi\sigma^{2}_{\ln R}}}e^{-\frac{1}{2}\Big{(}\frac{\ln(R/\bar{R})}{\sigma_{\ln R}}\Big{)}^{2}}.$$

$$⟨x_e(\text{x})⟩=1-e^{-n_b(\text{x}){\bar{V}}_b},$$

$$n_b(\text{x})={\bar{n}}_b(1+b{\delta }_L(\text{x})),$$

$${\bar{W}}_R^n(k)=\frac{1}{{\bar{V}}_b^n}{\int }_0^{\infty }dRP(R)[V_b(R)W(kR){]}^n.$$

$$X_l^{(x)}(k,M,z)={\bar{x}}_H\ln {\bar{x}}_H{b}_{\mathrm{bubble}}{\bar{W}}_R(k)u_1,$$

$$b_{\mathrm{bubble}}=\frac{1}{{\bar{n}}_b}{\int }_{M_{\mathrm{th}}}^{\infty }b(M,z)\frac{dn}{d\ln M}\frac{dM}{M}.$$

$$u(k,M,z)=\frac{1}{M}{\int }_0^{r_{\mathrm{vir}}}dr4\pi r^2\frac{\sin kr}{kr}{\rho }_{\mathrm{NFW}},$$

$$\rho_{\rm{NFW}}=\rho_{s}\Big{(}\frac{r}{r_{s}}\Big{)}^{-1}\Big{(}1+\frac{r}{r_{s}}\Big{)}^{-2}.$$

$$P_{\psi \psi }^{(1h)}$$

$${\xi }_{\psi \eta }={\bar{x}}_H({\xi }_{{\delta }_{\eta }\delta }+{\xi }_{{\delta }_{\eta }x\delta })+{\xi }_{x{\delta }_{\eta }},$$

$$P_{\psi \eta }={\bar{x}}_H(P_{{\delta }_{\eta }\delta }+P_{{\delta }_{\eta }x\delta })+P_{x{\delta }_{\eta }}.$$

$$P_{\psi \eta }^{(2h)}$$

$$P_{\psi \eta }^{(1h)}$$

$$P_{\psi \eta }\simeq P_{x{\delta }_{\eta }}^{(2h)}+{\bar{x}}_H{P}_{{\delta }_{\eta }\delta }^{(2h)}.$$

$$P^{1h,XY}(k,z)=\int dM\frac{dn}{dM}{X}_l(k,M,z)Y_l(k,M,z)$$

$$P^{2h,XY}(k,z)$$

$$\frac{\mathrm{SFR}(M,z)}{M_{\odot }{\mathrm{yr}}^{-1}}$$

$$L_{\mathrm{rec}}^{\mathrm{GAL}}(M,z)=1.55\times 10^{42}(1-f_{\mathrm{esc}})f_{\mathrm{Lya}}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}(\mathrm{erg}{\mathrm{s}}^{-1}).$$

$$L_{\mathrm{exc}}^{\mathrm{GAL}}(M,z)=4.03\times 10^{41}(1-f_{\mathrm{esc}})f_{\mathrm{Lya}}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}(\mathrm{erg}{\mathrm{s}}^{-1}),$$

$$L_{\mathrm{cooling}}^{\mathrm{GAL}}(M,z)$$

$$L_{\mathrm{stellar}}^{\mathrm{GAL}}(M,z)=5.12\times 10^{40}{f}_{\mathrm{Lya}}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}(\mathrm{erg}{\mathrm{s}}^{-1}).$$

$$X_l(k,M,z)=\frac{L(M,z)}{4\pi D_L^2}y{D}_A^{2}u(k,M,z).$$

Full text

I Introduction

The early universe, initially filled with hot plasma, became neutral as hydrogen ions captured electrons that were decoupled from cosmic microwave background (CMB) photons at a redshift of 1100. A cosmic “dark age” subsequently ensued in the universe until the linear density fluctuations seeded by inflation were amplified, forming the first stars and galaxies Loeb and Barkana (2001). The X-rays from mini quasars and ultraviolet radiation from the massive stars in first-light galaxies heated and ionized the neutral hydrogen, and the universe gradually transformed from completely neutral to fully ionized during the epoch of reionization (EoR). Today the EoR still remains largely unexplored, as signatures imprinted on the intergalactic medium (IGM) in the early universe are too faint to be detected.

The physical processes present during the EoR are of extreme importance to our understanding of the universe and the structure that formed in it. The clustering of neutral hydrogen (HI) down to the Jeans length scale contains a wealth of information about certain fundamental physics, including dark matter. The HI tomography is not subject to small-scale physical effects such as photon diffusion damping present in the CMB power spectrum. The timing and duration of the EoR can help interpret other cosmological measurements, such as the kinetic Sunyaev-Zel’dovich (kSZ) effect McQuinn et al. (2005). Moreover, some exotic physics such as primordial magnetic fields Schleicher et al. (2009) and decaying dark matter Furlanetto et al. (2006a) could be probed during the EoR. To date, the neutral fraction during the EoR was measured from quasar absorption spectra Fan et al. (2006) and $\rm{Ly}\alpha$-emitting galaxy luminosity functions Konno et al. (2014); Malhotra and Rhoads (2004, 2006) around $z\sim 6$. Another important quantity of the EoR, the Thomson scattering optical depth, is constrained to $\tau$ = 0.088 $\pm$ 0.014 by WMAP Komatsu et al. (2011) and $\tau$ = 0.058 $\pm$ 0.012 by Planck satellites Planck Collaboration et al. (2016).

The best way to measure the HI content prior to and during reionization is through the $\rm{21\,cm}$ HI fine-structure spin-flip transition. A number of experiments have been targeting the $\rm{21\,cm}$ emission, such as the Low Frequency Array (LOFAR) van Haarlem et al. (2013), the Murchison Widefield Array (MWA) Tingay et al. (2013), the Precision Array for Probing the Epoch of Reionization (PAPER) Parsons et al. (2010), the Hydrogen Epoch of Reionization Array (HERA) DeBoer et al. (2016) and the Square Kilometer Array (SKA) Koopmans et al. (2015). The redshifted $\rm{21\,cm}$ emission is contaminated by both galactic and extragalactic foregrounds that consist of galactic synchrotron, supernovae remnants, free-free emission, and radio point sources Furlanetto et al. (2006b). The Galactic synchrotron emission is the dominant contribution, as it is three to four orders of magnitude stronger than the background brightness temperature fluctuations. By performing a component separation or subtracting the $\rm{21\,cm}$ foreground, the $\rm{21\,cm}$ brightness fluctuations could be measured Bonaldi and Brown (2015). This, however, relies on the reliability of the foreground estimation. The radio point sources are also thought to be another foreground issue for $\rm{21\,cm}$ experiments, but this signal is very likely to be a subdominant contamination Liu et al. (2009a). The expected $\rm{21\,cm}$ signal is at the level of 10 $\rm{mK}^{2}$ at $k=0.3$ $\rm{Mpc}^{-1}$ Mesinger et al. (2011), while recent measurements from PAPER set a 2$\sigma$ upper limit as ($22.4\,\rm{mK})^{2}$ in the range $0.15<k<0.5$ $h\,\rm{Mpc}^{-1}$ at $z=8.4$ Ali et al. (2015).

During the EoR the ultraviolet $\rm{Ly}\alpha$ emission was created by the first stars and galaxies. The $\rm{Ly}\alpha$ background traces the underlying dark matter distribution and also affects the spin-temperature distribution. By directly measuring the $\rm{Ly}\alpha$ emissions, we get an additional observable on EoR physics as well Jensen et al. (2013). However, the $\rm{Ly}\alpha$ background is contaminated by low-$z$ foregrounds, such as $\rm{H}\alpha$ at $z$ = 0.5, OIII at $z$ = 0.9, and OII at $z$ = 1.6. These low-$z$ components are much brighter than $\rm{Ly}\alpha$, precluding a clean detection. On the other hand, such low-$z$ foregrounds can be easily masked out since they are very bright Pullen et al. (2014); Gong et al. (2014). Therefore, a simple masking procedure would recover the genuine $\rm{Ly}\alpha$ background from experiments.

The $\rm{21\,cm}$ and $\rm{Ly}\alpha$ emission is anti-correlated at large angular scales because they originate from IGM and galaxies, respectively, and ionized bubbles around $\rm{Ly}\alpha$ galaxies are devoid of HI that is seen with $\rm{21\,cm}$ experiments. The transition in the cross-correlation from negative to positive indicates a characteristic size for the average of HII regions around halos. Therefore, the cross-correlation between $\rm{21\,cm}$ and $\rm{Ly}\alpha$ can be viewed as a complementary probe of EoR physics. The cross-correlation could be more advantageous in terms of foreground removal as the two sets of aforementioned foregrounds would be largely uncorrelated, potentially allowing a higher signal-to-noise detection and an easy confirmation of the EoR signature. Previously, the cross-correlation between $\rm{21\,cm}$ experiments and galaxies was studied for $\rm{21\,cm}$ experiments such as MWA and LOFAR, using both analytical and numerical calculations Furlanetto and Lidz (2007); Lidz et al. (2009); Wyithe and Loeb (2007), as well as for LOFAR and Subaru’s Hyper Suprime-cam Vrbanec et al. (2016). The cross-correlation between $\rm{21\,cm}$ and CO/kSZ also shows a similar transition in the correlation sign Lidz et al. (2011); Jelić et al. (2010).

So far, different approaches have been used to model reionization. The large-scale $N$-body and radiative transfer simulations, while desirable, are challenging because it is computationally intensive due to the large dynamic range Iliev et al. (2006, 2015). Another approach involves semi-analytical/semi-numerical models by taking a halo catalog generated from $N$-body simulations and generating a reionization field by smoothly filtering the halo field Zahn et al. (2007, 2011). A more simplified idea of this semi-numerical simulation is to make the density field from Gaussian random variables instead of relying on the $N$-body simulations. The simulation can be done efficiently within a small box for EoR Mesinger and Furlanetto (2007); Mesinger et al. (2011). However, this numerical solution becomes ineffective when the box is too large or the simulated epoch is far beyond the EoR when the CMB temperature $T_{\rm{cmb}}$ is coupled to the spin temperature $T_{\rm{s}}$ and the assumption $T_{\rm{s}}\gg T_{\rm{cmb}}$ breaks down. An upgraded version of this implementation uses a very similar algorithm to extend to large boxes Santos et al. (2010). Here, we apply a very simple ionizing bubble model Furlanetto et al. (2004a) to the calculations of $\rm{21\,cm}$ brightness temperature anisotropy and its cross-correlation with $\rm{Ly}\alpha$ analytically, so we can quickly forecast the detectability of the signal for different combinations of $\rm{21\,cm}$ and $\rm{Ly}\alpha$ experiments, and explore the EoR parameter space without significant computational cost. This approach would be very beneficial when the cross-correlation measurements with different experiments and major foreground or instrumental issues need to be identified in the early stage of the development.

This paper is organized as follows. In Section II, we introduce the halo model for the ionizing bubble as well as the cross-correlation. In Section III, the $\rm{Ly}\alpha$ luminosity is discussed. Then we focus on the low-$z$ foregrounds for both $\rm{21\,cm}$ and $\rm{Ly}\alpha$ measurements in Section IV and estimate signal-to-noise for the detectability of different experiments, as well as the uncertainties on the EoR parameters in Section V. We conclude in Section VI. We use the Planck cosmological parameters: $\Omega_{b}h^{2}=0.02230$, $\Omega_{c}h^{2}=0.1188$, $H_{0}=67.74\,\rm{km}/\rm{s}/\rm{Mpc}$, $Y_{\rm p}=0.249$, $\ln(10^{10}A_{s})=3.064$ at $k_{\ast}=0.05\,\rm{Mpc}^{-1}$, $n_{s}=0.9667$, and $\tau=0.058$.

II Theoretical Model of the Cross-correlation

Here we describe the basic ingredients of our halo model. Since the mean ionizing fraction is not precisely constrained by current observations, we use the CAMB’s reionization model Lewis (2008); i.e.,

[TABLE]

where the redshift $z_{\rm{re}}$ is derived from the optical depth $\tau$ today; i.e.,

[TABLE]

and $\Delta y$ = 1.5\sqrt{1+z_{\rm{re}}}$$\Delta_{z}. Here, $\sigma_{T}$ is the Thomson cross-section, the electron density is $n_{e}=(1-3/4Y_{\rm p})\rho_{b,0}/m_{\rm{H}}a^{-3}x_{e}$, the comoving length $d\ell=ad\chi^{\prime}$, the Helium fraction is $Y_{\rm p}$, the proton mass is $m_{\rm{H}}$, and the mean neutral hydrogen fraction is $\bar{x}_{H}=1-\bar{x}_{e}$. We assume that helium is singly ionized along with hydrogen, while the double ionization of helium is neglected.

The $\rm{21\,cm}$ brightness temperature can be split into two components, $T_{0}(z)\psi(\textbf{x},z)$, in which the isotropic background temperature is

[TABLE]

and spatial fluctuation $\psi$ is Meerburg et al. (2013)

[TABLE]

Here $\delta_{x}$ is the density contrast of ionizing field ($x$) and we neglect perturbations introduced by spin-temperature fluctuations and peculiar velocities. For the $\rm{21\,cm}$ field, we only consider the signals from IGM as galaxy contributions are $10^{-4}$ times smaller Gong et al. (2011), and model the ionizing field with “bubbles” Meerburg et al. (2013). From Eq. (4), the two-point correlation functions for ionizing and matter density contrasts are $\langle\delta_{x}\delta_{x}\rangle=\xi_{xx}/\bar{x}^{2}_{H}$, $\langle\delta_{x}\delta\rangle=\xi_{x\delta}/\bar{x}_{H}$, and $\langle\delta\delta\rangle=\xi_{\delta\delta}$.

The auto-correlation function of the isotropic $\rm{21\,cm}$ spatial fluctuation field is Zaldarriaga et al. (2004); Furlanetto et al. (2004b)

[TABLE]

We should note that this $\rm{21\,cm}$ auto-correlation function is only an approximation in that the three-point correlation terms neglected are generally substantial, as Lidz et al. (2007) pointed out. However, we mainly focused on the cross-correlation calculations, for which we did consider all of the higher order terms. Although the assumption that the density field is Gaussian is a reasonable approximation on most of the scales of interest, we create the ionizing field from a Poisson process, as we will discuss later, to account for its non-Gaussianity. We only use power spectrum to do the statistics so that the non-Gaussianity of the field is not captured Zaldarriaga et al. (2004). We neglect the redshift distortions and make use of the fact that the spin temperature is significantly higher than CMB at $z<$ 10 Thomas and Zaroubi (2011). We Fourier transform the correlation function, assuming that the quadratic terms are negligibly small. The $\rm{21\,cm}$ power spectrum is

[TABLE]

The power spectrum can be calculated from a halo model by describing HII regions as bubbles. The two-halo term of $P_{x\delta x\delta}$ is higher order and negligible Furlanetto et al. (2004b).

The viral temperature of halo $T_{\rm{vir}}=5\times 10^{4}\,\rm{K}$ (suggested by Meerburg et al. (2013)) sets the minimum halo mass

[TABLE]

With this threshold mass, we can calculate the mean number density of bubble $\bar{n}_{b}=\int_{M_{\rm{th}}}^{\infty}dM\frac{dn}{dM}$ and the average bubble size from $\bar{n}_{b}=-(\ln\bar{x}_{H})/\bar{V}_{b}$. When it is compared to the predicted value

[TABLE]

the bubble radius is constrained. The bubble radius $R$ is assumed to satisfy a logarithmic distribution Lidz et al. (2009) as

[TABLE]

In Figure 1, we show the distribution of the bubble size at different redshifts.

Given the average volume and number density, the ionizing field is generated from a Poisson process

[TABLE]

and its number density is

[TABLE]

where the density contrast $\delta_{L}$ is the matter density $\delta$ smoothed by a top-hat window of radius $R$. The top-hat window in Fourier space is

[TABLE]

The shape factor for the ionizing field $x$ is defined as

[TABLE]

with the bubble bias given by

[TABLE]

Here, $u_{1}$ is the Fourier transform of the NFW profile Navarro et al. (1996); i.e., $u_{1}=M/\bar{\rho}_{m}u$. The NFW Fourier transform is

[TABLE]

which is derived from a standard NFW profile

[TABLE]

The detailed discussions of $\rho_{s}$, $r_{s}$ and $r_{\rm{vir}}$ can be found in Ref. Feng et al. (2016).

The 1-halo term of the $\rm{21\,cm}$ field Mortonson and Hu (2007); Wang and Hu (2006) is

[TABLE]

where $P_{xx}^{(\rm{1h})}$ = $(\bar{x}_{e}-\bar{x}_{e}^{2})V_{b}\bar{W}_{R}^{2}$ and $P_{x\delta x\delta}^{(\rm{1h})}$ = $(\bar{x}_{e}-\bar{x}_{e}^{2})\tilde{P}_{\delta\delta}$. The term $P_{x\delta}^{(\rm{1h})}$ is zero as we assume the bubble is completely ionized. Here $\sigma_{R}^{2}$ = $\int dkk^{2}/2\pi^{2}\bar{W}_{R}^{2}(k)P_{\delta\delta}(k)$ and $\tilde{P}_{\delta\delta}$ = $P_{\delta\delta}V_{b}\sigma_{R}^{2}/\sqrt{P_{\delta\delta}^{2}+(V_{b}\sigma_{R}^{2})^{2}}$. The two-halo term can be easily calculated with the shape factor in Eq. (13).

On the other hand, the intensity mapping (IM) of $\rm{Ly}\alpha$ emitters (LAEs) is a biased tracer of the same dark matter distribution, i.e., $\sim(1+\delta_{\eta})$. For simplicity, we use $\delta_{\rm{Ly}\alpha}=\delta_{\eta}$. The cross-correlation between $\rm{21\,cm}$ and LAEs is

[TABLE]

and the 3D power spectrum is

[TABLE]

Here, the two-halo and one-halo terms are given by

[TABLE]

and

[TABLE]

respectively. The subscript $\delta_{\eta}x\delta$ essentially refers to $\delta_{\eta},x\delta$ and the “,” is omitted for simplicity.

On small scales $P^{1h}_{\delta_{\eta}x\delta}$ cancels out the term $P^{1h}_{\delta_{\eta}\delta}$, so the summation is almost zero Lidz et al. (2009). Also, the large-scale information of $P^{2h}_{\delta_{\eta}x\delta}$ should be very negligible. With all of these approximations, the final power spectrum of the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-correlation is

[TABLE]

The halo-model approach, i.e.,

[TABLE]

and

[TABLE]

can be used to calculate each power spectrum in Eq. (22). In these equations, $dn/dM$ is the mass function and $b(M,z)$ is the bias. The linear matter power spectrum is $P_{\rm{lin}}$. We will work out the shape factors $X_{l}(k,M,z)$ and $Y_{l}(k,M,z)$ (or $\tilde{X}_{l}(k,M,z)$ and $\tilde{Y}_{l}(k,M,z)$) for $\rm{Ly}\alpha$ in the next section.

III $\rm{Ly}\alpha$ emission

The UV radiation emitted from massive and short-lived stars can ionize the neutral hydrogen in the interstellar medium (ISM) in galaxies and the number of ionizing photons closely depends on the star formation rate (SFR). In this work we consider an SFR model that is consistent with numerical simulations. The fitted SFR Silva et al. (2013) is

[TABLE]

where $a=2.8$, $b=-0.94$, $d=-1.7$, $M_{1}=10^{9}M_{\odot}$, and $M_{2}=7\times 10^{10}M_{\odot}$.

The ionizing photons could escape the galaxies with a fraction $f_{\rm{esc}}(M,z)=e^{-\alpha(M/M_{\odot})^{\beta}}$, but the remains will ionize the hydrogen and 66% of the ionization will result in a recombination process that produces $\rm{Ly}\alpha$ photons. The dust in the ISM can also absorb the $\rm{Ly}\alpha$ emissions, and the remaining fraction that survives the dust extinction is $f_{\rm{Lya}}(z)$. The luminosity due to the recombination is then calculated as

[TABLE]

The ionizing radiation can heat the gas so that the process of hydrogen excitation and cooling produces $\rm{Ly}\alpha$ emission as well. The luminosities due to excitation and cooling are

[TABLE]

and

[TABLE]

respectively.

Besides these line emissions, the continuum produces $\rm{Ly}\alpha$ photons through stellar radiation, free-free (ff), free-bound (fb), and two-photon (2$\gamma$) processes. Among these contributions, the stellar emission with a blackbody spectrum below the Lyman limit is dominant and its luminosity is

[TABLE]

Our calculation takes all of these continuum lines into account, and the detailed line luminosity can be found in Ref. Silva et al. (2013).

The total $\rm{Ly}\alpha$ luminosity $L(M,z)$ from a galaxy is a summation of all of the above components and the shape factor for $\rm{Ly}\alpha$ field is

[TABLE]

Here the conversion factor from frequency to comoving distance is $y$ = $d\chi/d\nu$ = $\lambda/(ca^{2})d\chi/dz$, $\lambda$ is the line rest frame wavelength, and $D_{L}$ and $D_{A}$ are luminosity and angular comoving distances, respectively.

This construction of shape factor is only a mathematical definition that facilitates the halo-model calculations in Eqs. (23) and (24). We should note that for individual point-like $\rm{Ly}\alpha$ emitters, it would appear to be extended due to spatial diffusion of $\rm{Ly}\alpha$ photons Zheng et al. (2010, 2011). Also, the scattering of photons on the red-side of $\rm{Ly}\alpha$ in the IGM would further damp the $\rm{Ly}\alpha$ flux along the line of sight Miralda-Escudé (1998). We first consider that the $\rm{Ly}\alpha$ emission is a biased tracer of the underlying dark matter distribution and phenomenologically account for the extended structure by the mass- and redshift-dependent quantity $b(M,z)$ in Eq. (24) and the luminosity function $L(M,z)$. Next we will discuss some dominating effects of IGM on the $\rm{Ly}\alpha$ emissions, but effects such as the damping wing of $\rm{Ly}\alpha$ have to rely on a numerical simulation.

The mean $\rm{Ly}\alpha$ intensity varies at different redshifts as

[TABLE]

where $M_{\rm{min}}=10^{8}M_{\odot}$ and $M_{\rm{max}}=10^{13}M_{\odot}$.

The escaped photons from galaxies can ionize the IGM, which can also emit $\rm{Ly}\alpha$ photons due to the recombination process. The recombination rate is

[TABLE]

where $n_{e}=x_{e}n_{b}$, $n_{\rm{HII}}=x_{e}n_{b}C$, $C=(1-Y_{\rm p})/(1-3Y_{\rm p}/4)$, and $\alpha$ is a case A comoving recombination coefficient. The luminosity function of IGM is $L^{\rm{IGM}}_{\rm{rec}}=f_{\rm{rec}}\dot{n}_{\rm{rec}}E_{\rm{Ly}\alpha}$ and the fraction $f_{\rm{rec}}$ is spin-temperature dependent. We show the contribution of IGM in Figure 2 and it is seen that the IGM contribution is negligible for both auto- and cross-power spectra.

Another IGM contribution to $\rm{Ly}\alpha$ emission comes from the scattering of Ly$n$ photons escaping from galaxies. From the previous calculations Silva et al. (2013, 2016); Pullen et al. (2014), it is found that the diffuse IGM contribution is a few orders of magnitude smaller than the galaxies. Therefore, we ignore this contribution to the overall $\rm{Ly}\alpha$ signal.

IV Low-z Foregrounds

The $\rm{Ly}\alpha$ emission at the EoR can be significantly contaminated by low-$z$ foregrounds. The foreground at $z_{f}$ projected onto the source plane $z_{s}$ becomes anisotropic as the wave vector of the foreground power spectrum in Fourier space becomes $\textbf{k}_{f\rightarrow s}=(\chi_{s}/\chi_{f}k_{\perp},\,y_{s}/y_{f}k_{\parallel})$, which is not radially symmetric. The low-$z$ foregrounds are identified as $\rm{H}\alpha$ [$6563\,\AA$, $z$ = 0.5], OIII [$5007\,\AA$, $z$ = 0.9], and OII [$3727\,\AA$, $z$ = 1.6] with luminosities $L_{H\alpha}=1.3\times 10^{41}\frac{\rm{SFR}}{M_{\odot}\rm{yr}^{-1}}$, $L_{OII}=7.1\times 10^{40}\frac{\rm{SFR}}{M_{\odot}\rm{yr}^{-1}}$, and $L_{\rm{OIII}}=1.3\times 10^{41}\frac{\rm{SFR}}{M_{\odot}\rm{yr}^{-1}}$, respectively. The low-$z$ SFR is exclusively modeled as

[TABLE]

for the foreground line emissions. This SFR model is fitted to the numerical simulations below $z$ = 2 and the parameters are constrained as $a$ = $-9.097$, $b$ = 0.484, $c$ = 2.7, $d$ = $-4.0$, $M_{1}$ = $10^{8}M_{\odot}$, and $M_{2}$ = $8\times 10^{11}M_{\odot}$ Gong et al. (2014).

The projected power spectrum of the foreground is then expressed as

[TABLE]

and $k_{f}$ = $\sqrt{(\chi_{s}/\chi_{f})^{2}k^{2}_{\perp}+(y_{s}/y_{f})^{2}k^{2}_{\parallel}}$. In Figure 3, we show the unprojected and projected power spectra for $\rm{H}\alpha$. OIII and OII show similar contours so they are neglected. In Figure 4, the blue and red curves are radially averaged from the anisotropic power spectra.

In Figures 4 and 5, we show the power spectra for $\rm{Ly}\alpha$ and the foreground lines $\rm{H}\alpha$, OIII, and OII with projection and with flux masking at redshifts $z$ = 7 and 9. The projected foreground emissions are much higher than the $\rm{Ly}\alpha$ lines. By selecting the brightest sources at the flux detection threshold and forming a mask, we can effectively remove those “hot” pixels which only account for a very tiny fraction of the sky coverage Pullen et al. (2014); Gong et al. (2014). As can be seen in Figure 6, we show the percentage of the removed pixels as the threshold flux changes. We find that a flux cut at $10^{-18}\rm{W/m^{2}}$ can significantly lower amplitudes of the foreground power spectra while only removing less than 0.1% of the pixels. Therefore, the flux masking procedure makes the low-$z$ foregrounds negligible.

The synchrotron radiation dominates the $\rm{21\,cm}$ signals but its smooth spectral feature can easily be used to isolate this component in frequency domain, also the galactic foregrounds are not correlated with extragalactic line emissions at low-$z$. So we do not expect any noticeable cross-correlations between galactic synchrotron and $\rm{Ly}\alpha$ foregrounds. However, the radio point sources that are too faint to be resolved are indeed correlated with the low-$z$ foregrounds within $0.5<z<1.6$, so this component would be picked up in the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-correlation, making it not as systematic-free as expected. To estimate its contribution, we use the model in Refs. Gleser et al. (2008); Singal et al. (2010); Serra et al. (2008). The model is described as

[TABLE]

based on the fact that the radio point source is a tracer of underlying density field. Here we have defined $I_{\nu}^{\rm{radio}}$ = $\int_{0}^{S_{\rm{cut}}}dSSdN/dS$ and assume a flux limit $S_{\rm{cut}}$ = 1$\,\rm{mJy}$, above which the radio point sources are bright enough to be resolved. The flux distribution is a simple power-law, i.e., $dN/dS\,=\,A(S/S_{0})^{\alpha}$, where $A$ = 4$\,\rm{mJy}^{-1}\,\rm{sr}^{-1}$, $S_{0}$ = 880$\,\rm{mJy}$, and $\alpha\,=\,-1.75$ Liu et al. (2009b). Also, $\partial B_{\nu}/\partial T$ = 99.27\,\rm{Jy}\,\rm{sr}^{-1}/(\mu K)$$x^{4}e^{x}$$/(e^{x}-$1)^{2}$ and $x$ = $h\nu/k_{b}T_{\rm{CMB}}$ = $\nu/56.84\,(\rm{GHz})$. In Fourier space, the shape functions of the halo model for the radio sources are described as

[TABLE]

and

[TABLE]

which are directly inserted into Eqs. (23) and (24) to obtain the one-halo and two-halo terms for the point-source clustering. Both $X_{l}(k,M,z)$ and $\tilde{X}_{l}(k,M,z)$ describe the Fourier-space profile of a point source with mass $M$ located at redshift $z$. Here the central and satellite galaxy numbers are

[TABLE]

and

[TABLE]

The mean galaxy number density is

[TABLE]

where $N_{g}(M)=N_{c}(M)+N_{s}(M)$. The parameters determined from luminosity and color dependence of galaxy clustering in the SDSS DR7 main galaxy sample are $M_{\rm{min}}$ = $10^{9}M_{\odot}$, $\sigma_{M}$ = 0.2, $M_{s}$ = $5\times 10^{10}M_{\odot}$, and $\alpha_{s}$ = 1 Zehavi et al. (2011).

We estimate that the radio foreground contributions at the $\rm{Ly}\alpha$ foreground redshifts and the raw power spectra are a few orders of magnitude higher than the $\rm{21\,cm}$ signal as revealed by Liu et al. (2009b); Alonso et al. (2014). Therefore, the foreground suppression is very crucial; the spectral fitting procedure studied in Liu et al. (2009b) demonstrated that the radio foregrounds can be reduced by six orders of magnitude in map space, and it has been validated that this is true from flux cut 0.1–100 mJy. Consequently, the radio foreground contamination becomes negligible, and we show all of the power spectra in Figures 7 and 8 at redshifts $z$ = 7 and 9. We see that the resulting radio point sources have very negligible contaminating power on the $\rm{21\,cm}$ measurements. Finally, we show the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-power spectra at redshifts $z$ = 7 and 9 in Figures 9 and 10 with both foreground separation schemes incorporated. In Figure 11, we show both the evolutions and cross-correlation coefficients of the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-power spectrum as a function of $k$ from $z$ = 6 to $z$ = 10.

Despite the fact that all of the foreground cross-correlations are small, this has to rely on the assumption that we have very efficient foreground mitigation strategies for both $\rm{21\,cm}$ and $\rm{Ly}\alpha$ measurements. Non-negligible foreground residuals on $\rm{21\,cm}$ and $\rm{Ly}\alpha$ would make a great impact on the power spectrum uncertainties, even if they are uncorrelated.

V Forecast for the Experiments

In this section, we consider two experiments, SKA and Cosmic Dawn Intensity Mapper (CDIM) Pritchard et al. (2015); Cooray et al. (2016), and investigate the detectability of the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-correlation with both instrumental noises and foregrounds. We list the experimental specifications in Tables 1 and 2.

The instrumental noise of the $\rm{21\,cm}$ experiment is entirely determined by some key factors such as integration time $t_{0}$, system temperature $T_{\rm{sys}}$, maximum baseline $D_{\rm{max}}$, collecting area $A_{\rm{tot}}$, antenna number $N_{a}$, and frequency resolution $d\nu$. The noise is then given by

[TABLE]

For the $\rm{Ly}\alpha$ experiments, the noise is

[TABLE]

where the comoving volume subtended by the detector pixel

[TABLE]

depends on the pixel area $A_{\rm{pix}}$ and frequency resolution $\delta_{\nu}$. The general Knox formula Bowman et al. (2006) for the measured signal $\rm{21\,cm}$ or $\rm{Ly}\alpha$ is

[TABLE]

and the number of modes in the bin $k$ is $N_{m}\,=\,2\pi k^{2}\Delta k\frac{V_{s}}{(2\pi)^{3}}$. Here $X,\,Y$ = $\{\rm{21\,cm},\,\rm{Ly}\alpha\}$, the survey volume is $V_{s}\,=\,\chi^{2}A_{s}yB_{\nu}$, $A_{s}$ is the survey area, and $B_{\nu}$ is the bandwidth (BW). For $\rm{Ly}\alpha$ experiments, the minimum and maximum scales are determined by the survey and pixel areas. For $\rm{21\,cm}$, we normally consider modes at scales below $k=10\,\rm{Mpc}^{-1}$ and get the minimum $k$ from the total survey area. For the cross-correlation, the common $k$ range is chosen from two experiments and the minimum volume between the $\rm{21\,cm}$ and $\rm{Ly}\alpha$ experiments is taken to calculate the number of modes. The noise- and foreground-included power spectra are formed as $\tilde{P}_{\rm{21\,cm}}\,=\,P_{\rm{21\,cm}}+P^{N}_{\rm{21\,cm}}+P^{\rm{radio}}_{\rm{21\,cm}}$, $\tilde{P}_{\rm{Ly}\alpha}\,=\,P_{\rm{Ly}\alpha}+P^{N}_{\rm{Ly}\alpha}+P^{\rm{low}\hbox{-}z}_{\rm{Ly}\alpha}$, and $\tilde{P}_{\rm{21\,cm}\mbox{--}\rm{Ly}\alpha}\,=\,P_{\rm{21\,cm}\mbox{--}\rm{Ly}\alpha}+P^{\rm{radio}\mbox{--}\rm{low}\hbox{-}z}_{\rm{21\,cm}\mbox{--}\rm{Ly}\alpha}$. The $k$ region is binned in log-space and we calculate the errors using Eq. (V) for each $k$-band.

Following Ref. Dillon et al. (2014), we also consider the foreground wedge in Fourier space and its impact on the power spectrum sensitivity. For SKA, the characteristic angle $\theta_{0}$ is chosen to be $\sqrt{\rm{FOV}}\sim 4^{\circ}$. The wedge cut reduces the effective modes $N_{m}$ in Fourier space and decreases the overall signal-to-noise by 8%, which could be much larger if a bigger angle $\theta_{0}$ is assumed. In Figure 12, we show the reduced region by the foreground wedge using a 2D $\rm{21\,cm}$ power spectrum. In addition to the modes cut by the wedge, low $k_{\parallel}$ modes are strongly contaminated by foreground and should be excluded as well. A large $k_{\rm max}$ can compensate for the loss of modes due to the horizontal cut. For example, at $k_{\rm max}\,=\,10\,h\rm{Mpc}^{-1}$, a horizontal cut between $0.05\,h{\rm Mpc}^{-1}<k_{\parallel}<0.1\,h{\rm Mpc}^{-1}$ introduces negligible changes to the overall signal-to-noise. But for a small $k_{\rm max}$, such as $0.5\,h\rm{Mpc}^{-1}$, the total wedge could reduce the signal-to-noise by 11% with a horizontal cut at $k_{\parallel}$ = 0.1 $h\rm{Mpc}^{-1}$.

In Figures 13, 14 and 15, we show all of the power spectra and their band errors for $\rm{21\,cm}$ and $\rm{Ly}\alpha$ at $z$ = 7 and 9. As can be seen from Figure 15, the anti-correlations between neutral hydrogen and galaxies can be probed at very high signal-to-noise ratios.

From the forecasted power spectra in Figure 15, we can further try to constrain the EoR parameters defined as P = {$\tau$, $\Delta y$, $\sigma_{\ln R}$} = {0.058, 6.0, 1.0}. The Fisher matrix Pober et al. (2014) is

[TABLE]

Here $\Delta\eta(k,z)$ is the error on the cross-power spectrum $\Delta^{2}$ and $p_{i}$ refers to any parameters in the set P and $\Delta^{2}$ = $k^{3}/(2\pi^{2})P_{\rm{21\,cm}\mbox{--}\rm{Ly}\alpha}$. All of the 1$\sigma$ confidence levels, as well as the likelihood functions in Figure 16, are calculated from the $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-power spectra at $z$ = 7 and 9. As can be seen in the figure, the bubble size can be constrained from the cross-correlation while the errors on the optical depth and duration of the reionization transition are large. This is due to the fact that the cross-correlation is proportional to $x_{e}(z)$ and not $x^{2}_{e}(z)$, to which the $\rm{21\,cm}$–$\rm{21\,cm}$ and $\rm{21\,cm}$–$\tau$ Meerburg et al. (2013) are proportional. Therefore, the cross-power spectrum at a single redshift is less sensitive to the reionization history. However, the cross-correlations at different redshifts might be able to break the degeneracies among the parameters, and so are useful as complementary probes to cosmological and astrophysical problems. In Figure 17, we reduce the bandwidth by a factor of 5 and combine the cross-power spectra measured at $z$ = 6, 7, 8, 9, and 10. It is seen that the Fisher matrix error bars on all of the parameters are significantly reduced, and that multi-redshift measurements within narrow bins can effectively break the parameter degeneracies.

VI Conclusion

In this work, we applied a bubble model to the computation of $\rm{21\,cm}$ and $\rm{Ly}\alpha$ cross-correlation at the EoR. Making use of the empirical relation between $\rm{Ly}\alpha$ luminosity and mass for the line emissions, we also calculated the power spectra for $\rm{Ly}\alpha$. The $\rm{21\,cm}$–$\rm{Ly}\alpha$ cross-power spectrum in this fast approach can reproduce the key features of the one made by detailed numerical simulations, and we can use it to quickly assess the overall performance of the future EoR experiments.

The cross-correlation is contaminated by the low-$z$ foregrounds for both $\rm{21\,cm}$ and $\rm{Ly}\alpha$. We studied the radio galaxies for the $\rm{21\,cm}$ experiments and H$\alpha$, OIII, and OII line emissions for the $\rm{Ly}\alpha$ experiments. All of these foregrounds could be a few orders of magnitude higher than the signals we are probing if the foreground mitigation is not incorporated. The map-space spectral fitting can effectively remove the radio point-source contaminations, while a flux masking for the intensity mapping experiments have been shown to be a good and easy foreground-removal method.

We take advantage of this efficient algorithm and estimate the errors on the EoR parameters $\tau$, $\Delta y$ and $\sigma_{\ln R}$, based on the Fisher matrix formalism. For other physical processes during the EoR, such as X-ray heating, supernovae explosion, and shock heating, numerical simulations with these effects or an extension to this work should be devised. We will discuss these in the future works.

VII Acknowledgements

A.C. and C.F. acknowledge support from NASA grants NASA NNX16AJ69G, NASA NNX16AF39G and Ax Foundation for Cosmology at UC San Diego.