PROBING REIONIZATION WITH THE CROSS-POWER SPECTRUM OF 21 cm AND NEAR-INFRARED RADIATION BACKGROUNDS

Ignoring peculiar velocities, the 21 cm brightness temperature relative to the cosmic microwave background (CMB) at an observed frequency ν (i.e., redshift z) can be written as (e.g., Field 1958, 1959; Zaldarriaga et al. 2004)

$$\begin{eqnarray} T_{21}(\hat{\mathbf {n}},z)&\approx& 26 \,\bar{x}_{{{\rm H}\,\scriptsize{I}}}(1+\delta _{\rm \rho })(1+\delta _{\rm x}) \left(\frac{T_s-T_{\rm CMB}}{T_s}\right) \left(\frac{\Omega _{\rm b}h^2}{0.022}\right)\nonumber\\ &&\times \left(\frac{1+z}{10}\,\frac{0.15}{\Omega _{\rm M}h^2}\right)^{1/2}\, {\rm mK}. \end{eqnarray} \tag{ 1 }$$

In this equation, $\bar{x}_{{{\rm H}\,\scriptsize{I}}}$ is the volume-averaged neutral hydrogen fraction. The field δ_ρ is the fractional gas density fluctuation, while $\delta _{\rm x}=(x_{{{\rm H}\,\scriptsize{I}}}-\bar{x}_{{{\rm H}\,\scriptsize{I}}})/\bar{x}_{{{\rm H}\,\scriptsize{I}}}$ is the perturbation in the neutral hydrogen fraction. T_s is the spin temperature of the 21 cm transition and T_CMB is the CMB temperature. The other symbols have their usual meanings. Here, we adopt the simplifying assumption that T_s ≫ T_CMB globally, which should be a good approximation during most of the reionization epoch (Mesinger & Furlanetto 2007). We can now write the observed 21 cm brightness temperature on the sky as

$$\begin{equation} T_{21}(\hat{\mathbf {n}},z)=T_0(z)\,\bar{x}_{{{\rm H}\,\scriptsize{I}}}\,\int dr\,W(r)\,\psi (\hat{\mathbf {n}},r), \end{equation} \tag{ 2 }$$

where we define

$$\begin{equation} T_0=26\,\left(\frac{\Omega _{\rm b}h^2}{0.022}\right) \left(\frac{1+z}{10}\,\frac{0.15}{\Omega _{\rm M}h^2}\right)^{1/2}\, {\rm mK} \end{equation} \tag{ 3 }$$

as the mean 21 cm brightness temperature of the IGM, r is the radial distance around the target redshift z, and $\psi (\hat{\mathbf {n}},r)= (1+\delta _{\rm \rho })(1+\delta _{\rm x})$. Note that W(r) is the window function representing the frequency resolution of the telescope.

We first decompose the 21 cm temperature field into spherical harmonic multiple moments,

$$\begin{eqnarray} a^{21}_{\ell m}&=&\int d \hat{\mathbf {n}}\,T_{21}(\hat{\mathbf {n}},z)Y^{\ast }_{\ell m}(\hat{\mathbf {n}}) \nonumber \\ &=&4\pi i^{\ell }\int \frac{d^3 k}{(2\pi)^3}\,\psi _{k}(\mathbf {k})\,I^{21}_{\ell }(k,z)\,Y^{\ast }_{\ell m}(\mathbf {k})\protect, \protect \end{eqnarray} \tag{ 4 }$$

in which

$$\begin{equation} I^{21}_{\ell }(k,z)=\int dr\,W^{21}(r,z)j_{\ell }(kr), \end{equation} \tag{ 5 }$$

with the weighting function

$$\begin{equation} W^{21}(r,z)=T_0(z)\,\bar{x}_{{{\rm H}\,\scriptsize{I}}}\,W(r). \end{equation} \tag{ 6 }$$

Next, we calculate the angular power spectrum of the 21 cm fluctuations, which is defined by the two-point correlation function

$$\begin{equation} \left\langle a^{21}_{\ell m}a^{21}_{\ell ^{\prime } m^{\prime }}\right\rangle =\delta _{\ell \ell ^{\prime }}\delta _{mm^{\prime }}C^{\psi \psi }_{\ell }. \end{equation} \tag{ 7 }$$

Using the weighting function in radial space, we can finally construct the angular power spectrum as a line-of-sight projection (Song et al. 2003)

$$\begin{equation} C^{\psi \psi }_{\ell }(z)=\int \frac{dr}{r^2}[W^{21}(r,z)]^2\,P_{21}\left(k=\frac{\ell }{r},z\right). \end{equation} \tag{ 8 }$$

We have used the Limber approximation (Limber 1954) here. Furthermore, the 21 cm power spectrum P₂₁(k, z) can be decomposed into the sum of several terms (generalizing the formula from Furlanetto et al. 2006a)

$$\begin{eqnarray} P_{21}(k,z)&=&P_{\rm {\rho,\rho }}(k,z)+P_{\rm {x,x}}(k,z)+2P_{\rm {\rho,x}}(k,z)+2P_{\rm {x\rho,x}}(k,z)\nonumber\\ &&+2P_{\rm {x\rho,\rho }}(k,z)+P_{\rm {x\rho,x\rho }}(k,z). \end{eqnarray} \tag{ 9 }$$

The quantity P_{a, b} indicates the dimensional power spectrum of two random fields, a and b. As shown in Lidz et al. (2007), the last three terms, referred to as the higher-order terms, contribute significantly to the 21 cm power spectrum only at small scales k > 1 Mpc⁻¹, and hence are ignored in our calculations. The first three terms represent the matter power spectrum, the power spectrum of neutral hydrogen fluctuations, and the cross-correlation power, respectively. We have neglected the cross-correlation P_{ρ, x}(k, z) because it is small on the scales of interest (Santos et al. 2003).

Under the halo model, we can write the matter power spectrum as (Seljak 2000; Cooray & Sheth 2002)

$$\begin{eqnarray} P_{\rm {\rho,\rho }}(k,z)&=&P^{1h}_{\rm {\rho,\rho }}(k,z)+P^{2h}_{\rm {\rho,\rho }}(k,z) \nonumber \\ P^{1h}_{\rm {\rho,\rho }}(k,z)&=&\int dM \frac{dn(z,M)}{dM}\,\left |\frac{\rho _h(z,M,k)}{\bar{\rho }(z)}\right |^2 \nonumber \\ P^{2h}_{\rm {\rho,\rho }}(k,z)&=&\int dM_1\,\frac{dn(z,M_1)}{dM_1}\,\frac{\rho _h(z,M_1,k)}{\bar{\rho }(z)}\,b_h(z,M_1) \nonumber \\ &&\times \int dM_2\,\frac{dn(z,M_2)}{dM_2}\,\frac{\rho _h(z,M_2,k)}{\bar{\rho }(z)}\,b_h(z,M_2) \nonumber \\ &&\times P_{\rm {lin}}(k,z). \end{eqnarray} \tag{ 10 }$$

Here, M is the halo mass, dn/dM(z, M) is the halo mass function (Sheth & Tormen 1999), ρ_h(z, M, k) is the Fourier transform of the Navarro–Frenk–White halo density profile (Navarro et al. 1997), $\bar{\rho }(z)$ is the cosmic mean matter density, $b_h(z,M_)$ is the halo bias, and P_lin(k, z) denotes the linear matter power spectrum. For precise calculations, we have included the baryon content in the computation of the matter power spectrum.

Considering the case where $\bar{x}_{{{\rm H}\,\scriptsize{I}}}\,\delta _{\rm x}=-\bar{x}_{\rm i}\,\delta _{\rm i}$, P_{x, x} can be expressed in terms of the ionization fraction equation:

$$\begin{equation} {\bar{x}_{{{\rm H}\,\scriptsize{I}}}}^2\,P_{\rm {x,x}}(k,z)={\bar{x}_{\rm i}}^2\,P_{\rm {i,i}}(k,z). \end{equation} \tag{ 11 }$$

Here, $\delta _{\rm i}=(x_{\rm i}-\bar{x}_{\rm i})/\bar{x}_{\rm i}$ is the perturbation in the ionization fraction and P_{i, i}(k, z) is the power spectrum from δ_i. Under the assumption that the ionizing radiation was coming from stars formed from gas clouds that cooled in dark matter halos, we are thus able to write the auto correlation of the ionization fraction field as a multiple of the dark matter correlation (Santos et al. 2003, 2005):

$$\begin{equation} P_{\rm {i,i}}(k,z)={b^2_{\rm {eff}}}(z)\,P_{\rm {\rho,\rho }}(k,z)\,e^{-k^2R^2}, \end{equation} \tag{ 12 }$$

where b_eff(z) is the mean bias weighted by the different halo properties and R is the mean radius of the H ii patches. We follow the same approach as Santos et al. (2003) and Santos et al. (2005) for these two variables.

We now need to calculate the hydrogen ionization fraction $\bar{x}_{\rm i}$ as a function of redshift. Following Cooray et al. (2012), $\bar{x}_{\rm i}$ can be estimated as

$$\begin{equation} \frac{d \bar{x}_{\rm i}}{d t}=\frac{f_{\rm {esc}}\,\varphi (z)q(z)}{\bar{n}_{\rm H}(z)}-\frac{\bar{x}_{\rm i}}{\bar{t}_{\rm {rec}}}\protect, \protect \end{equation} \tag{ 13 }$$

where f_esc is the escape fraction of the ionization photons and $\bar{n}_{\rm H}(z)=1.905\times 10^{-7}(1+z)^3\,\rm {cm}^{-3}$ is the mean hydrogen number density (Shull et al. 2012). Here, φ(z) is the comoving star formation rate density (SFRD) given by the ongoing star formation model (Santos et al. 2002). The function q(z) is defined as $q(z)\equiv (\bar{Q}_{{{\rm H}\,\scriptsize{I}}}/\langle M_{\ast } \rangle)\,\langle \tau _{\ast } \rangle$, where $\bar{Q}_{{{\rm H}\,\scriptsize{I}}}$ is the time-averaged hydrogen photoionization rate, 〈M_*〉 is the average star mass, and 〈τ_*〉 is the average stellar lifetime (Lejeune & Schaerer 2001; Schaerer 2002). The $\bar{t}_{\rm {rec}}$ denotes the volume-averaged recombination time (Madau et al. 1999).

Figure 1 shows the evolution of the hydrogen ionization fraction with redshift. Hereafter, we set the star formation efficiency as f_* = 0.04, the escape fraction of ionization photons as f_esc = 0.5, and a gas temperature of T = 3 × 10⁴ K. In addition, M_min = 10⁶ M_☉ and M_max = 10¹⁴ M_☉ denote the minimum and maximum halo mass to host galaxies, respectively. With these parameters, we find that reionization ends around a redshift of six, which is consistent with current studies. Furthermore, we follow the approach of Cooray et al. (2012) to estimate the electron-scattering optical depth with the ionization fraction. We find that the derived optical depth is τ = 0.084, which is close to the result of the WMAP seven-year data with τ = 0.088 ± 0.014 (Komatsu et al. 2011).

Zoom In Zoom Out Reset image size

Similar to Knox et al. (2001), we write the antenna temperature of the NIRB at a given frequency ν and toward a direction $\hat{\mathbf {n}}$ as a product of the mean NIRB emissivity and its fluctuation:

$$\begin{eqnarray} &&T_{\rm NIR}(\hat{\mathbf {n}},\nu)=\int dr\,a(r)\,\bar{j}(\nu,r)\,\left[1+\frac{\delta j(r\hat{\mathbf {n}},\nu,r)}{\bar{j}(\nu,r)}\right]. \end{eqnarray} \tag{ 14 }$$

Here, r is the radial distance from us to a redshift of z, a = (1 + z)⁻¹ is the scale factor, and $\bar{j}(\nu,z)$ is the mean emissivity per comoving unit volume at frequency ν and redshift z. To model the spatial fluctuations related to the emissivity, we assume that the sources of this emissivity are galaxies, such that $\delta j(r\hat{\mathbf {n}},\nu,z)/\bar{j}(\nu,z)=\delta _{\rm gal}$, where $n_{\rm gal}=\bar{n}_{\rm gal}(1+\delta _{\rm gal})$ is the comoving number density of galaxies. Note that the NIRB would probe all of the galaxies as a whole, not just the brightest ones that have been detected individually in galaxy surveys. For a specific NIR experiment, we integrate $\bar{j}(\nu,r)$ over the observed bandpass and further define the observed NIRB temperature as

$$\begin{eqnarray} &&T_{\rm NIR}(\hat{\mathbf {n}})=\int dr\,a^2(r)\,\bar{j}(r)\,\phi (\hat{\mathbf {n}},r), \end{eqnarray} \tag{ 15 }$$

where $\bar{j}(z)$ is the band-averaged mean emissivity and $\phi (\hat{\mathbf {n}},r)=1+\delta _{\rm gal}$. Note that here we have a factor a² instead of a in Equation (14) since we are no longer estimating the temperature at only one frequency. This dependence has been explained in the Appendix of Fernandez et al. (2010).

In this work, we consider the emission from stars in the first-light galaxies present during reionization and neglect the emission from the IGM. There are two stellar populations in our calculation. The first, referred to as Pop II stars, are metal-poor stars with a metallicity of Z = 1/50 Z_☉, and the second, Pop III stars, are metal-free stars with Z = 0. For Pop II stars, we adopt the stellar initial mass functions (IMF) given by Salpeter (1995) with a mass range from 3 to 150 M_☉. For Pop III stars, we make use of the IMF obtained by Larson (1999), and the mass range is from 3 to 500 M_☉. Collecting all of the emission from these two stellar populations, we obtain the band-averaged mean emissivity as

$$\begin{equation} \bar{j}(z)=f_{\rm P}\,\bar{j}^{\rm {PopIII}}(z)+(1-f_{\rm P})\,\bar{j}^{\rm {PopII}}(z). \end{equation} \tag{ 16 }$$

Here, f_P(z) is the population fraction describing the relative fraction of the Pop II and Pop III stars at different redshifts (Cooray et al. 2012). The functions $\bar{j}^{\rm {PopII}}(z)$ and $\bar{j}^{\rm {PopIII}}(z)$ denote the band-averaged mean emissivities for Pop II and Pop III stars, respectively (Cooray et al. 2012):

$$\begin{equation} \bar{j}^{i}(z)=\frac{1}{4\pi }\,\bar{l}^{i}\,\langle \tau ^{i}_{\ast }\rangle \,\varphi (z), \end{equation} \tag{ 17 }$$

where 〈τ_*〉 is the mean stellar lifetime of each of the stellar type and φ(z) is the comoving SFRD. We integrate the luminosity mass density l_ν over the band of observed frequencies ν₁ to ν₂ to obtain the band-averaged $\bar{l}$ as (Fernandez et al. 2010)

$$\begin{equation} \bar{l}(z)=\int _{\nu _1(1+z)}^{\nu _2(1+z)} d\nu \,l_{\nu }(z). \end{equation} \tag{ 18 }$$

Following previous studies in the literature (Fernandez & Komatsu 2006; Fernandez et al. 2010; Cooray et al. 2012), we examine the luminosity mass density l_ν(z) for various emission processes that contribute to the NIRB, such as stellar emission, Lyman-α line, and free–free, free-bound, and two-photon emission. Here, we do not present the detailed theoretical calculations of the luminosity mass density for each radiative process and the basic formulae can be found in Fernandez & Komatsu (2006). We finally derive the total luminosity mass density from the stellar nebulae

$$\begin{equation} l_{\nu }=l^{\ast }_{\nu }+(1-f_{\rm {esc}})\left (l^{\rm {Ly\alpha }}_{\nu }+l^{\rm {ff}}_{\nu }+l^{\rm {fb}}_{\nu }+l^{\rm {2ph}}_{\nu }\right), \end{equation} \tag{ 19 }$$

where f_esc = 0.5 is the escape fraction of the ionization photons and $l^{i}_{\nu }$ denotes the luminosity mass density for a certain emission process. The relevant stellar parameters, such as the intrinsic bolometric luminosity, the effective temperature, the main-sequence lifetime, and the time-averaged hydrogen photoionization rate, are calculated using the fitting results from Lejeune & Schaerer (2001) and Schaerer (2002). Moreover, the NIR bandpasses ν₁ to ν₂, which are taken from Sterken & Manfroid (1992), are defined in Table 1.

In Figure 2, we show the luminosity mass density l_ν as a function of the rest-frame wavelength λ for Pop II and Pop III stars at z = 10. To specify the various emission processes contributing to the NIRB, we plot the l_ν of stellar emission (dot dashed line), the Lyman-α line (dotted line), and the free–free (thin solid line), free–bound (triple-dot dashed line), and two photon (dashed line) emission. For Pop II stars (left panel), we can see that the stellar emission is dominant over the wavelength range λ ∼ 0.1–1 μm. However, the reprocessed light, such as free–bound and two photon, is comparable to or even larger than the stellar spectrum for Pop III stars (right panel). As illustrated in these two panels, the total l_ν (thick solid line) from Pop II stars is similar to that from Pop III stars. Similar results were measured in Fernandez & Komatsu (2006).

Zoom In Zoom Out Reset image size

If we expand $\phi (\hat{\mathbf {n}},r)$ in a Fourier series and use the spherical harmonic decomposition in exact analogy with Equations (4)–(8), we can define the NIRB angular power spectrum given by

$$\begin{eqnarray} \!\!\!\!\!\!\!C^{\phi \phi }_{\ell }&=&\int \frac{dr}{r^2}[W^{\rm NIR}(z)]^2\,P_{\rm {gal,gal}}\left(k=\frac{\ell }{r},z\right) \nonumber \\ \!\!\!\!\!\!\!&=&\int _{z_{\rm {min}}}^{z_{\rm {max}}} dz\!\frac{dr}{dz}\,\frac{1}{r^2} [W^{\rm {NIR}}(z) ]^2\,P_{\rm {gal,gal}}\left(\!k=\frac{\ell }{r},z\!\right). \end{eqnarray} \tag{ 20 }$$

Here,

$$\begin{equation} W^{\rm NIR}(z)=a^2(z)\,\bar{j}(z) \end{equation} \tag{ 21 }$$

is the weighting function. Here, we examine the high-redshift component of the NIRB, coming from the first galaxies during reionization. For NIR observation in the H band, we calculate the angular power spectrum by integrating Equation (20) over the redshift range from z_min = 6 to z_max = 17. We do not consider the emission coming from sources at z > 17 since at the center wavelength 1.65 μm, photons should be more energetic than hν = 13.6 eV in the rest frame beyond z = 17. For NIR observations in the L band, the redshift interval of the NIRB integral is (z_min, z_max) = (6, 30). The function P_{gal, gal}(k, z) is the power spectrum of first galaxies. The spatial distribution of these galaxies is described through the relation of the tracer field to the dark matter halo distribution. Besides halo model, the other important ingredient is how these galaxies occupy dark matter halos. We use the HOD model for the first generation of galaxies. Then, the galaxy power spectrum can be written as

$$\begin{equation} P_{\rm {gal,gal}}(k,z)=P^{1h}_{\rm {gal,gal}}(k,z)+P^{2h}_{\rm {gal,gal}}(k,z). \end{equation} \tag{ 22 }$$

Here, $P^{1h}_{\rm {gal,gal}}$ and $P^{2h}_{\rm {gal,gal}}$ denote the power spectra contributed by galaxies in a single dark matter halo and galaxies in two different dark matter halos, respectively, which can be given by (Cooray & Sheth 2002)

$$\begin{eqnarray} P^{1h}_{\rm {gal,gal}}(k,z)&=&\int dM \frac{dn(z,M)}{dM}\,\frac{\langle N_{\rm {gal}}(M)[N_{\rm {gal}}(M)-1]\rangle }{\bar{n}^2_{\rm {gal}}(z)}\nonumber\\ &&\times \left [\frac{\rho _h(z,M,k)}{M}\right ]^p \nonumber \\ P^{2h}_{\rm {gal,gal}}(k,z)&=&\left [\int dM\,\frac{dn(z,M)}{dM}\, \frac{\langle N_{\rm {gal}}(M)\rangle }{\bar{n}_{\rm {gal}}(z)}\,\frac{\rho _h(z,M,k)}{M}\right. \nonumber\\ &&\left. \times b_h(z,M)\right ]^2 P_{\rm {lin}}(k,z). \end{eqnarray} \tag{ 23 }$$

We use N_gal(M) to indicate the average number of galaxies in each dark matter halo of mass M. The $\bar{n}_{\rm {gal}}(z)$ is the mean number density of the galaxies. The parameter p = 1 when 〈N_gal(N_gal − 1)〉 ⩽ 1, and p = 2 otherwise (Cooray & Sheth 2002).

In Figure 3, we show the 21 cm (top) and NIRB (bottom) anisotropy angular power spectra predicted by the theoretical model. Here, Δ² = ℓ(ℓ + 1)C_ℓ/2π denotes the dimensionless power spectrum. For the 21 cm experiment, we take the width of one frequency bin Δν = 0.1 MHz. Note that the amplitudes of Δ_NIR shown here are nearly 10 times smaller than those shown in Figure 12 of Cooray et al. (2012). This is because we represent the band-averaged power spectrum rather than that at only one wavelength. As mentioned above, the difference between our Equation (15) and the calculation in Cooray et al. (2012) is a factor of a = 1/(1 + z), which is only ∼0.1 during the EoR. Finally, the shot noise is ignored in our computation since it is much smaller than the expected NIRB signal at l < 10⁴ (Fernandez et al. 2010; Cooray et al. 2012).

Zoom In Zoom Out Reset image size

Based on the derived auto-correlations of 21 cm and NIR backgrounds, we examine their cross-correlation following Song et al. (2003)

$$\begin{equation} C^{\psi \phi }_{\ell }(z)=\int \frac{dr}{r^2}\,W^{21}(r,z)\,W^{\rm NIR}(z)\,P_{\rm {21,gal}}\left(k=\frac{\ell }{r},z\right). \end{equation} \tag{ 24 }$$

In this equation, P_{21, gal}(k, z) denotes the cross-power spectrum between the 21 cm fluctuation and the galaxy overdensity, which can be further decomposed into the sum of several contributing terms (Lidz et al. 2009):

$$\begin{equation} P_{\rm {21,gal}}(k,z)=P_{\rm {\rho,gal}}(k,z)+P_{\rm {x,gal}}(k,z)+P_{\rm {x\rho,gal}}(k,z). \end{equation} \tag{ 25 }$$

The individual terms on the right-hand side have the following physical interpretations. The first term, P_{ρ, gal}(k, z), is the cross-power spectrum between the matter and galaxy overdensity fields. The second term, P_{x, gal}(k, z), represents the cross-power spectrum between the neutral hydrogen fraction and the galaxy overdensity fields. The final term, P_{xρ, gal}, is a three-field term which contributes significantly only on small scales (Lidz et al. 2009), and hence is ignored in our calculation.

The halo model provides a simple calculation of the cross-power spectrum between matter and galaxies, and we have (Seljak 2000; Cooray & Sheth 2002)

$$\begin{eqnarray} P_{\rm {\rho,gal}}(k,z)&=&P^{1h}_{\rm {\rho,gal}}(k,z)+P^{2h}_{\rm {\rho,gal}}(k,z) \nonumber \\\ms P^{1h}_{\rm {\rho,gal}}(k,z)&=&\int dM \frac{dn(z,M)}{dM}\,\frac{M}{\bar{\rho }(z)}\,\frac{\langle N_{\rm {gal}}(M)\rangle }{\bar{n}_{\rm {gal}}(z)}\left [\!\frac{\rho _h(z,M,k)}{M}\!\right ]^p \nonumber \\ P^{2h}_{\rm {\rho,gal}}(k,z)&=&\int dM_1\,\frac{dn(z,M_1)}{dM_1}\,\frac{\rho _h(z,M_1,k)}{\bar{\rho }(z)}\,b_h(z,M_1) \nonumber \\\ms &&\times \int dM_2\,\frac{dn(z,M_2)}{dM_2}\, \frac{\langle N_{\rm {gal}}(M_2)\rangle }{\bar{n}_{\rm {gal}}(z)}\,\frac{\rho _h(z,M_2,k)}{M_2}\nonumber\\\ms &&\times b_h(z,M_2) P_{\rm {lin}}(k,z). \end{eqnarray} \tag{ 26 }$$

As described above, the matter–galaxy cross-correlation has the Poisson and halo-halo terms. $P^{1h}_{\rm {\rho,gal}}(k,z)$ denotes the correlation between galaxies and dark matter in the same halo and dominates on small scales. $P^{2h}_{\rm {\rho,gal}}(k,z)$ includes the correlations between galaxies and dark matter in neighboring halos and is dominant on large scales.

The next step is to compute the cross-correlation between the neutral hydrogen fraction and galaxy overdensity fields. In order to evaluate P_{x, gal} numerically under the analytic model, we focus on its behavior on scales much larger than bubble sizes. With this simplification, we can rewrite Equation (12) as

$$\begin{eqnarray} P_{\rm {i,i}}(k,z)&=&b^2_{\rm {eff}}\,P_{\rho,\rho }(k,z) \nonumber \\ &=&b^2_{\rm {eff}}\,\bar{b}^2_h\,P_{\rm {lin}}(k,z), \end{eqnarray} \tag{ 27 }$$

where $\bar{b}_h$ is the mean halo bias. Note that the filter function in Equation (12) is no longer included here. On large scales, the galaxy power spectrum can also be simplified to (Cooray & Sheth 2002)

$$\begin{equation} P_{\rm {gal,gal}}(k,z)=\bar{b}^2_{\rm {gal}}\,P_{\rm {lin}}(k,z), \end{equation} \tag{ 28 }$$

where $\bar{b}_{\rm {gal}}$ denotes the mean bias factor of the galaxy population. Then, the cross-correlation between the ionization fraction and the galaxy density fields, P_{i, gal}, can be given by

$$\begin{equation} P_{\rm {i,gal}}(k,z)=b_{\rm {eff}}\,\bar{b}_h\,\bar{b}_{\rm {gal}}\,P_{\rm {lin}}(k,z). \end{equation} \tag{ 29 }$$

However, as described in Equation (25), P_{21, gal}(k, z) depends on P_{x, gal}(k, z). This is obviously $\bar{x}_{{{\rm H}\,\scriptsize{I}}}\,P_{\rm {x,gal}}(k,z)=-\bar{x}_{\rm i}\,P_{\rm {i,gal}}(k,z)$. δ_x and δ_gal are therefore anti-correlated.

Figure 4 shows the 21 cm–galaxy cross-power spectrum as well as its components for z = (7, 10, 13), and the mean fractions of ionized hydrogen are $\bar{x}_{i}=(0.97,0.48,0.12)$, correspondingly. Here, the dimensionless cross-power spectrum is defined as $\Delta ^2_{\rm {a,b}}(k)=k^3P_{\rm {a,b}}/2\pi ^2$. The dashed lines represent the density–galaxy cross-power spectrum, $\Delta ^2_{\rm {\rho,gal}}$. We find that this power is always positive and increases by roughly four orders of magnitude over the scale range 0.1 < k < 100 Mpc⁻¹ at a redshift of z = 10. Compared to the matter power spectrum, $\Delta ^2_{\rm {\rho,\rho }}$, these high-redshift galaxies show very strong clustering even on large scales: the cross-power spectrum has an amplitude of $\Delta ^2_{\rm {\rho,gal}} \simeq 6.6\times 10^{-3}$ on a scale of k = 1 Mpc⁻¹ at z = 10. On the same scale, the amplitude of the matter power spectrum is $\Delta ^2_{\rm {\rho,\rho }} \simeq 8.1\times 10^{-4}$, which is ∼8 times smaller than that of $\Delta ^2_{\rm {\rho,gal}}$. On small scales, the cross-power spectrum $\Delta ^2_{\rm {\rho,gal}}$ grows exponentially, indicating that structure formation is highly biased to overdense regions during the EoR.

Zoom In Zoom Out Reset image size

We plot the absolute value of the cross-power spectrum between the neutral hydrogen fraction and galaxy density, $|\Delta ^2_{\rm {x,gal}}|$, as the dotted lines. It is worth noting that $\Delta ^2_{\rm {x,gal}}$ is negative at every value of k throughout the EoR. This occurs because galaxies form first and ionize their surroundings in overdense regions, whereas underdense regions are still mostly neutral and free of galaxies. Therefore, an anti-correlation between the neutral fraction and galaxy fields is naturally expected, especially on large scales. One can see that the strength of $|\Delta ^2_{\rm {x,gal}}|$ grows with the wavenumber. Moreover, as reionization proceeds, the cross-power is enhanced rapidly and finally exceeds the matter–galaxy cross-power spectrum.

From Equation (25), it is clear that the $\Delta ^2_{\rm {\rho,gal}}$ term and the $\Delta ^2_{\rm {x,gal}}$ term will together determine the quantity of $\Delta ^2_{\rm {21,gal}}$. We plot the absolute value of the 21 cm galaxy cross-power spectrum as the solid lines in Figure 4. Near the beginning of reionization, the 21 cm and galaxy density fields are positively correlated. However, these two fields quickly become anti-correlated on large scales. In our model, an interesting turnover appears when z ∼ 11 with $\bar{x}_{i}\sim 0.3$ and still persists at lower redshifts. Such a turnover does not exist in the 21 cm–galaxy cross-power spectrum at early times ($z=13,\,\bar{x}_{i}=0.12$). At z = 10, $\Delta ^2_{\rm {21,gal}}$ turns over on a scale of k ∼ 30 Mpc⁻¹. The turnover behavior occurs on progressively smaller scales as reionization proceeds. Since we use an approximate expression for $\Delta ^2_{\rm {x,gal}}$, the turnover presented here should be considered a qualitative feature. Actually, the specific shape of the turnover depends on the details of the reionization model (Lidz et al. 2009; Wiersma et al. 2013). Furthermore, one may note that the amplitudes of Δ²(k) shown here are slightly smaller than those shown in the simulations from Lidz et al. (2009), likely due to the smaller value of the low-mass cutoff M_min = 10⁶ M_☉ adopted in our model of galaxy formation (Scoccimarro et al. 2001).

In addition, we use the cross-correlation coefficient to quantify the relative strength of the cross-correlation between the 21 cm and NIR radiation backgrounds:

$$\begin{equation} r_{\rm {corr}}(\ell)=\frac{C^{\psi \phi }_{\ell }}{\sqrt{C^{\psi \psi }_{\ell }\,C^{\phi \phi }_{\ell }}}. \end{equation} \tag{ 30 }$$

In Figure 5, we plot predictions for the 21 cm–NIRB cross-power spectrum as well as the cross-correlation coefficient. We correlate the high-redshift component of the NIRB observed in the H band (left) and the L band (right) with 21 cm anisotropy at individual redshifts of z = (7, 8, 9, 10, 11, 12, 13, 14, 15), respectively. The absolute values of the cross-correlation signals are shown in the upper panels, where we have $\Delta ^2=\ell (\ell +1)C^{\psi \phi }_{\ell }/2\pi$. In our model, the 21 cm and NIR backgrounds are positively correlated in the early stages of reionization. At early times, galaxies are extremely rare objects and are only just starting to ionize their surroundings. The galaxies turn on first in large-scale overdense regions. Compared to the underdense regions, these areas contain more matter and initially more neutral hydrogen, and consequently glow more brightly in the 21 cm emission. This leads to a positive 21 cm–NIRB cross-correlation on large scales, as shown by the thin green lines in Figure 5. As reionization proceeds, the galaxies quickly ionize their overdense surroundings completely. Hence, these overdense regions dim in the 21 cm emission. On the other hand, the large-scale underdense regions are still mostly free of galaxies and roughly maintain their initial 21 cm brightness temperature. Therefore, an anti-correlation is demonstrated during the middle and late stages of reionization, as shown by the thick red lines in Figure 5. Overall, the amplitude of the cross-correlation grows even faster with redshift after the reionization is half completed, and the maximum is reached when z = 7 with $\bar{x}_{\rm {i}}\approx 0.97$. This marks a rise of the cosmic structure formation which is responsible for the reionization of neutral hydrogen at late times.

Zoom In Zoom Out Reset image size

In the lower panels, we show the cross-correlation coefficient for various 21 cm redshift slices. One can see that the absolute value of r_corr reaches a maximum of ∼0.05 when $\bar{x}_{\rm {i}}\sim 0.6$ for the two NIR bands considered here. Hence, we infer that only a small part of the NIRB cross-correlates with the 21 cm emission from the EoR. The result is easy to understand: the NIRB contains the remnant light of high-redshift galaxies. Specifically, the redshift intervals of the NIRB integral are 6 ⩽ z ⩽ 17 and 6 ⩽ z ⩽ 30 for the H band and the L band, respectively. However, the bandwidth of the 21 cm signal is significantly narrow Δν = 0.1 MHz, corresponding to Δz = 7 × 10⁻³ at z = 9. In the cross-correlation analysis, most of the NIR light comes from redshifts different from that of the measured 21 cm background. Fortunately, a cross-correlation does exist between these two radiation backgrounds. Any detection of this cross signal in future experiments could tell us about the properties of early stars and the physics of the reionization process.