Canceling out intensity mapping foregrounds

TL;DR

This paper introduces a novel analysis method combining intensity maps with galaxy surveys to effectively cancel foreground contamination in 21 cm intensity mapping, enabling more precise measurements of the HI mass function.

Contribution

The paper proposes a new one-point statistic that is unbiased by foregrounds and leverages combined data to improve HI mass function measurements.

Findings

Method can measure HI mass function with high precision

Foreground bias is effectively canceled

Applicable to current 21 cm observations

Abstract

21 cm intensity mapping has arisen as a powerful probe of the high-redshift universe, but its potential is limited by extremely bright foregrounds and high source confusion. In this Letter, we propose a new analysis which can help solve both problems. From the combination of an intensity map with an overlapping galaxy survey we construct a new one-point statistic which is unbiased by foregrounds and contains information left out of conventional analyses. We show that our method can measure the HI mass function with unprecedented precision using observations similar to recent 21 cm detections.

Equations

$$\varphi (M_{\mathrm{HI}})=\ln (10){\varphi }_{\ast }{\left(\frac{M_{\mathrm{HI}}}{M_{\ast }}\right)}^{1+\alpha }{e}^{-M_{\mathrm{HI}}/M_{\ast }-M_{\min }/M_{\mathrm{HI}}}.$$

$${\mathcal{P}}_{1+2}(T)=\int \int {\mathcal{P}}_1(T^{\prime }){\mathcal{P}}_2(T^{\prime \prime }){\delta }_D(T-T^{\prime }-T^{\prime \prime })d{T}^{\prime }d{T}^{\prime \prime },$$

$${\mathcal{P}}_{1+2}(T)=\int {\mathcal{P}}_1(T^{\prime }){\mathcal{P}}_2(T-T^{\prime })d{T}^{\prime }=({\mathcal{P}}_1\circ {\mathcal{P}}_2)(T),$$

$${\tilde{\mathcal{P}}}_{1+2}(\mathcal{T})={\tilde{\mathcal{P}}}_1(\mathcal{T}){\tilde{\mathcal{P}}}_2(\mathcal{T}),$$

$$\frac{\tilde{\mathcal{P}}(\mathcal{T}|N_{\det }^1)}{\tilde{\mathcal{P}}(\mathcal{T}|N_{\det }^2)}=\frac{{\tilde{\mathcal{P}}}_S(\mathcal{T}|N_{\det }^1){\tilde{\mathcal{P}}}_{\mathrm{FG}}(\mathcal{T})}{{\tilde{\mathcal{P}}}_S(\mathcal{T}|N_{\det }^2){\tilde{\mathcal{P}}}_{\mathrm{FG}}(\mathcal{T})}=\frac{{\tilde{\mathcal{P}}}_S(\mathcal{T}|N_{\det }^1)}{{\tilde{\mathcal{P}}}_S(\mathcal{T}|N_{\det }^2)}.$$

$${\tilde{\mathcal{R}}}_i^{N_{\det }^1{N}_{\det }^2}=\frac{{\tilde{B}}_i^{N_{\det }^1}}{{\tilde{B}}_i^{N_{\det }^2}}.$$

$$\text{var}(B_i)=⟨B_i⟩(1-⟨B_i⟩/N_{\mathrm{vox}})\approx ⟨B_i⟩.$$

$$\tilde{\mathcal{P}}(\mathcal{T}|N_{\det })={\tilde{\mathcal{P}}}_{\det }(\mathcal{T}|N_{\det }){\tilde{\mathcal{P}}}_{\mathrm{un}}(\mathcal{T}|N_{\det }).$$

$${\varphi }_{\mathrm{un}}(M_{\mathrm{HI}})=\varphi (M_{\mathrm{HI}})e^{-M_{\mathrm{HI}}/M_{\mathrm{cut}}},$$

$${\varphi }_{\det }(M_{\mathrm{HI}})=\varphi (M_{\mathrm{HI}})\left(1-e^{-M_{\mathrm{HI}}/M_{\mathrm{cut}}}\right),$$

$${\tilde{\mathcal{P}}}_{\det }(\mathcal{T}|N_{\det })={\left[{\tilde{\mathcal{P}}}_{\det }(\mathcal{T}|N_{\det }=1)\right]}^{N_{\det }}.$$

$${\tilde{\mathcal{P}}}_{\mathrm{un}}(\mathcal{T}|N_{\det })=\sum _{N_{\mathrm{un}}}{\left[{\tilde{\mathcal{P}}}_{\mathrm{un}}(\mathcal{T}|N_{\mathrm{un}}=1)\right]}^{N_{\mathrm{un}}}\mathcal{P}(N_{\mathrm{un}}|N_{\det }).$$

$${\mathcal{P}}_{\mathrm{LN}}(\mu )=\frac{1}{\mu \sqrt{2\pi {\sigma }_G^2}}\exp \left\{-\frac{1}{2{\sigma }_G^2}{\left[\ln \left(\frac{\mu }{\overline{N}}\right)+\frac{{\sigma }_G^2}{2}\right]}^2\right\},$$

$$\mathcal{P}(N_{\mathrm{un}}|N_{\det })=\int \mathcal{P}(N_{\mathrm{un}}|f_{\mathrm{un}}{\mu }_{\det })\mathcal{P}({\mu }_{\det }|N_{\det })d{\mu }_{\det }.$$

$$\mathcal{P}({\mu }_{\det }|N_{\det })\propto \mathcal{P}(N_{\det }|{\mu }_{\det }){\mathcal{P}}_{\mathrm{LN}}({\mu }_{\det }),$$

$$\begin{gathered} \mathcal{P}(N_{\mathrm{un}}|N_{\det })\propto \int {\mathcal{P}}_{\mathrm{Poiss}}(N_{\mathrm{un}}|f_{\mathrm{un}}{\mu }_{\det }){\mathcal{P}}_{\mathrm{LN}}({\mu }_{\det }) \\ \times {\mathcal{P}}_{\mathrm{Poiss}}(N_{\det }|{\mu }_{\det })d{\mu }_{\det }, \end{gathered}$$

Full text

Many experiments are studying the evolution of the universe with the redshifted 21 cm line Tingay et al. (2013); van Haarlem et al. (2013); Bandura et al. (2014); Ali et al. (2015); Xu et al. (2015); Newburgh et al. (2016); DeBoer et al. (2017). These experiments seek to map vast swaths of the sky from the local universe to cosmic dawn using line intensity mapping Suginohara et al. (1999); Morales and Wyithe (2010); Pritchard and Loeb (2012); Kovetz et al. (2017). Intensity maps do not resolve individual emitters, but instead map fluctuations in the density of neutral hydrogen. This gives them sensitivity to the aggregate emission from all galaxies, as well as the neutral intergalactic medium. By targeting a narrow emission line, maps can be made in three dimensions by observing at many closely-spaced frequencies. Intensity mapping surveys can quickly access large volumes of space, allowing unprecedented constraints on cosmology and fundamental physics Shaw et al. (2015); Liu et al. (2016); Ewall-Wice et al. (2016); Liu and Parsons (2016); Foreman et al. (2018). In addition, they can provide statistical information about faint objects below detection thresholds of conventional surveys. Theory predicts that the number of detected objects is dwarfed by these systems, which should then accrete and merge to form galaxies such as our Milky Way Somerville and Davé (2015).

Currently, the largest limiting factor in 21 cm cosmology is the presence of foreground emission which is typically orders of magnitude brighter than the signal Oh and Mack (2003); Wang et al. (2006); Liu and Tegmark (2011). These foregrounds have limited any attempts to observe HI in autocorrelation. There have, however, been detections of cross-spectra between 21 maps and optical galaxy surveys, at redshift $z\sim 0.08$ Anderson et al. (2018) with Parkes telescope data and the 2dF galaxy survey Colless et al. (2001), and another at $z\sim 0.8$ with Green Bank Telescope (GBT) data Masui et al. (2013) and the WiggleZ survey Parkinson et al. (2012). Cross-correlation is robust against foreground contamination because galaxies and HI trace the same underlying large-scale structure, while the foregrounds do not.

Even when the signal can be detected it is challenging to interpret. Intensity maps are typically analyzed using power spectra. A density field is fully described by its power spectrum if that field is Gaussian. However, 21 cm maps weight galaxies by their gas content, which is determined by complex, nonlinear baryon dynamics within dark matter halos. As a result, the resulting intensity field is non-Gaussian.

For example, consider an experiment like Parkes/2dF or GBT/WiggleZ. At $z<1$, virtually all of the neutral hydrogen is found within halos, and we can describe how it is distributed with the HI mass function (HIMF) $\phi_{\rm{HI}}(M_{\rm{HI}})$, which gives the number density of halos with HI masses between $M_{\rm{HI}}$ and $M_{\rm{HI}}+dM_{\rm{HI}}$. Let us consider an HIMF with a slightly-modified Schechter form:

[TABLE]

This basic form has been used many times in the literature Schechter (1976); Zwaan et al. (2005); Jones et al. (2018). We have added an additional low-mass cutoff at $M_{\rm{min}}$, necessary as intensity maps lack hard detection thresholds. A power spectrum, sensitive to only the Gaussian part of a field, can only access the first two moments of the HIMF Lidz et al. (2011), while it would take at least four numbers, $\left(\phi_{*},M_{*},\alpha,M_{\rm{min}}\right)$, to fully determine Eq. (1). Furthermore, there are substantial degeneracies between $\phi$ and fundamental cosmological parameters that intensity maps might otherwise be ideally situated to measure, such as halo bias, the growth rate of density fluctuations, or the amplitude of primordial non-Gaussianity Obuljen et al. (2018); Chen et al. (2018); Castorina and White (2019).

It was suggested in Refs. Breysse et al. (2016, 2017) that this non-Gaussianity could be accessed using one-point statistics, as opposed to two-point statistics like power spectra. A technique called P(D) analysis, which has seen use for decades across the electromagnetic spectrum Scheuer (1957); Barcons et al. (1994); Windridge and Phillipps (2000); Lee et al. (2009); Patanchon et al. (2009); Glenn et al. (2010), allows mapping between the HIMF and the probability distribution function (PDF) $\mathcal{P}(T)$ of a voxel (or three-dimensional pixel) having intensity between $T$ and $T+dT$. This statistic, termed the Voxel Intensity Distribution (VID), has been shown to add significantly to the astrophysical information which can be gained from an intensity map Ihle et al. (2019).

Unfortunately, as with the auto-spectrum, the VID of a 21 cm map would be dangerously contaminated with foregrounds. Even after foreground cleaning, there can be substantial residual contamination. Rather than expose ourselves to this large potential source of bias, we introduce here a one-point analogue to the cross-spectrum, extending the VID formalism to remove foreground bias. We will then forecast how this technique can be applied to the Parkes/2dF and GBT/WiggleZ data.

The approach we propose here relies on a simple fact: for independent random variables $T_{1}$ and $T_{2}$, the PDF $\mathcal{P}_{1+2}(T)$ of their sum $T=T_{1}+T_{2}$ is the convolution of their individual PDFs $\mathcal{P}_{1}(T)$ and $\mathcal{P}_{2}(T)$. This is straightforward to prove. We can write

[TABLE]

where $\delta_{D}$ is a Dirac delta function. Evaluating one integral leaves a convolution,

[TABLE]

We can further make use of the Fourier convolution theorem to turn this into a product in Fourier space:

[TABLE]

where $\mathcal{T}\equiv 2\pi/T$ is the Fourier conjugate of the intensity $T$ and $\tilde{\mathcal{P}}(\mathcal{T})$ is the Fourier transform of $\mathcal{P}(T)$ (also known as the characteristic function). In an intensity map, the observed $T$ in a given voxel is the sum of a signal component $T_{S}$ with PDF $\mathcal{P}_{S}$ and a contribution from noise and foregrounds which we will abbreviate as $T_{\rm{FG}}$ with PDF $\mathcal{P}_{\rm{FG}}$. The full VID of a map is then the convolution of $\mathcal{P}_{S}$ and $\mathcal{P}_{\rm{FG}}$.

As stated above, $\mathcal{P}_{\rm{FG}}$ is hard to model to sufficient precision, so it will be difficult to apply the VID statistic directly to 21 cm data. We will therefore use a separate data set with different systematics to isolate our signal. Assume that our volume contains both a 21 cm map and an optical galaxy survey. In each voxel we know the total radio intensity $T_{S}+T_{FG}$ and the number $N_{\rm{det}}$ of detected optical galaxies. $T_{S}$ and $N_{\rm{det}}$ will be correlated, due both to large-scale structure and the HI content of the optical galaxies. We can therefore construct conditional PDFs $\mathcal{P}(T|N_{\rm{det}})$. We refer to these PDFs as conditional Voxel Intensity Distributions, or CVIDs.

Crucially, each CVID will be a convolution of a signal part, which depends on $N_{\rm{det}}$, and a noise/foreground part, which does not. If we compare voxels with different $N_{\rm{det}}$, we can write,

[TABLE]

From the second equality, it is clear that the above ratio is unbiased by foregrounds, as the deconvolution cancels out the component which is common to both CVIDs.

In practice, we do not compute continuous PDFs directly from maps. Instead, we estimate PDFs using histograms $B_{i}\approx\mathcal{P}(T_{i})\Delta TN_{\rm{vox}}$, where $B_{i}$ is the number of voxels in a bin of width $\Delta T$ centered at $T_{i}$, and $N_{\rm{vox}}$ is the total number of voxels. If we separate our map by $N_{\rm{det}}$ values, compute histograms $B_{i}^{N_{\rm{det}}}$ from each part, then compute their Fourier transforms $\tilde{B}_{i}^{N_{\rm{det}}}=\Delta T\sum_{j}B_{j}^{N_{\rm{det}}}\exp(i\kappa_{i}T_{j})$, we can write down the CVID Ratio (CVR),

[TABLE]

It is easy to show that the expectation value of $\tilde{\mathcal{R}}^{N_{\rm{det}}^{1}N_{\rm{det}}^{2}}_{i}$ is proportional to the ratio from Eq. (5), and thus is unbiased by foregrounds.

Figure 1 shows an example of how this works. The left panel shows predicted histograms for a Parkes/2dF data set with a Gaussian $\mathcal{P}_{\rm{FG}}$, with the signal contribution computed as described below. The right panel shows the CVR computed only from the signal. Both panels include, for illustration, a toy data realization. If we assume that $T_{S}$ and $T_{\rm{FG}}$ in each voxel are independent draws from $\mathcal{P}_{S}$ and $\mathcal{P}_{\rm{FG}}$, then each $B_{i}$ is an independent draw from a binomial distribution with mean $B_{i}$ and variance

[TABLE]

Based on the nearly-diagonal correlation matrices shown in Figure 2 of Ref. Ihle et al. (2019), this approximation seems reasonable, at least for low signal-to-noise data. We can therefore create a sample histogram by drawing randomly from binomial distributions, with errors following Eq. (7). Detailed testing of this binomial approximation will require actual mock maps, which we leave for future work.

From our toy data, it is clear that the simulation, which includes foregrounds, gives the same CVR as the foreground-free theory. Though we used a Gaussian $\mathcal{P}_{\rm{FG}}$, the same would hold for any general PDF. If we can estimate the error on $B_{i}^{N_{\rm{det}}}$, either with Eq. (7) or with simulations, then we can propagate this error through to get the error on the CVR directly from the data. This has two important implications. First, it means that we do not need a model of the foregrounds to estimate errors on $\tilde{\mathcal{R}}^{10}_{i}$. Second, it accounts for instabilities which appear when the denominator of Eq. (6) approaches zero, as seen at high-$\mathcal{T}$ in Figure 1. Though the measured CVR deviates significantly from the expectation, the large error bar means that these points get correspondingly little weight.

Now we need to connect the CVR to the HIMF. We can relate the two using a modified $P(D)$ analysis Lee et al. (2009); Breysse et al. (2017); Chen et al. (2017). In each voxel, there are $N_{\rm{det}}$ detected optical galaxies and an unknown number $N_{\rm{un}}$ of unresolved HI emitters, so:

[TABLE]

To compute $\tilde{\mathcal{P}}_{\rm{det}}$ and $\tilde{\mathcal{P}}_{\rm{un}}$, we need separate HIMFs for the detected and undetected galaxy populations. For now, we assume an exponential separation,

[TABLE]

[TABLE]

with free parameter $M_{\rm{cut}}$.

If $N_{\rm{det}}=0$, there is no contribution from detected galaxies and $\mathcal{P}_{\rm{det}}(T|0)=\delta_{D}(T)$. If $N_{\rm{det}}=1$, then $\mathcal{P}_{\rm{det}}(T|1)$ is proportional to $\phi_{\rm{det}}$, with appropriate normalization (see Breysse et al. (2017) for details). For higher values of $N_{\rm{det}}$, we can recursively apply Eq. (4) to get

[TABLE]

Note that in Eq. (11) we have implicitly assumed that the HIMF in a voxel is independent of how many galaxies it contains. This is known to be inaccurate, as more massive objects will be more strongly biased, shifting the HIMF to larger masses in dense voxels. Accurately modeling this effect will likely require simulations, so we neglect it for now.

For undetected galaxies, we do not know the value of $N_{\rm{un}}$, so we have

[TABLE]

where $\mathcal{P}_{\rm{un}}(\mathcal{T}|N_{\rm{un}}=1)$ is proportional to $\phi_{\rm{un}}$, $\mathcal{P}(N_{\rm{un}}|N_{\rm{det}})$ describes the correlation between optical galaxies and unresolved HI emitters. If there is no clustering, $\mathcal{P}_{\rm{un}}$ has no dependence on $N_{\rm{det}}$, and will cancel out of the CVR. However, we know that HI and optical galaxies should be at least somewhat correlated, though that correlation may be color- and scale-dependent Hess and Wilcots (2013); Papastergis et al. (2013); Anderson et al. (2018). Here we will employ a simple model to capture the leading order effects.

Following Breysse et al. (2017); Agrawal et al. (2017), assume each voxel has expected counts $\mu_{\rm{un}}$ and $\mu_{\rm{det}}$ of detected and undetected galaxies, and that the actual counts are Poisson draws with these means. We take $\mu_{\rm{un}}=f_{\rm{un}}\mu_{\rm{det}}$, where $f_{\rm{un}}$ is the ratio $\overline{N}_{\rm{un}}/\overline{N}_{\rm{det}}$ of the mean number counts. Galaxy count distributions are known to be reasonably-approximated by lognormal distributions Coles and Jones (1991); Kayo et al. (2001),

[TABLE]

where the $\sigma_{G}$ parameter acts as the “bias” in this model. The uncorrelated case corresponds here to $\sigma_{G}=0$. We then can write

[TABLE]

We can use Bayes’ Theorem to state that,

[TABLE]

where $\mathcal{P}_{\rm{LN}}$ acts as our “prior”. With our assumption of Poisson statistics, we then have

[TABLE]

where $\mathcal{P}_{\rm{Poiss}}(N|\mu)$ is the Poisson distribution with mean $\mu$.

We can now finally predict a CVR from our HIMF model. We will now examine what information could be gained from such an analysis. Consider a model of the Parkes/2dF maps as described in Ref. Anderson et al. (2018), with free parameters $(\phi_{*},M_{*},\alpha,M_{\rm{min}},M_{\rm{cut}},\sigma_{G})$. Assume that HI evolves negligibly from $z=0-0.05$, and use the best-fit HIMF from the ALFALFA survey Jones et al. (2018) as a model, with $\phi_{*}=(4.5\pm 0.8)\times 10^{-3}$ Mpc*-3* dex*-1*, $\log(M_{*}/M_{\odot})=9.94\pm 0.05$, and $\alpha=-1.25\pm 0.1$. We arbitrarily choose $M_{\rm{min}}=10^{5}\ M_{\odot}$ so that it falls below the ALFALFA detection threshold, and $M_{\rm{cut}}=3\times 10^{8}\ M_{\odot}$ to get the correct number of optically detected galaxies Colless et al. (2001). We compute $\sigma_{G}$ following Ref. Breysse et al. (2017). We assume a Gaussian $\mathcal{P}_{\rm{FG}}$, but note again that this procedure would work regardless of the assumed form. As the strength of the contamination depends on how much prior foreground cleaning is assumed, we consider a pessimistic case with $\sigma_{\rm{FG}}=35$ mK and an optimistic case with $\sigma_{N}=10$ mK, arbitrarily chosen to very roughly match the Parkes maps from Figure 1 of Ref. Anderson et al. (2018). Note that we do not allow $\phi_{*}$ and $\alpha$ to vary between the detected and undetected populations, but assume a single value for each.

We can forecast constraints on our model using the Fisher matrix formalism Fisher (1935); Tegmark et al. (1997). We only consider $N_{\rm{det}}=0$ and 1 here, as few data voxels have $N_{\rm{det}}>1$. As we neglect evolution from the $z=0$ ALFALFA galaxies, we can use their quoted systematic errors as priors on $\phi_{*}$, $M_{*}$ and $\alpha$. Though we do not directly forecast the higher-redshift GBT/WiggleZ survey, we can get a rough idea of its performance using a weaker set of priors, which we take as ten times worse than the ALFALFA priors. In both cases, we assume $10\%$ prior knowledge of $\sigma_{G}$, which would have to come from simulations, and uninformative fractional priors of 10 on both $M_{\rm{cut}}$ and $M_{\rm{min}}$. We assume a cosmology consistent with the Planck 2015 results Planck Collaboration et al. (2016).

Figure 2 shows the results of our Fisher forecasts. We plot $1\sigma$ confidence intervals around our fiducial models for the total, detected, and undetected populations. Even in the most pessimistic case, with strong foregrounds and weak (10x ALFALFA) priors, we get a good measurement of the bright end of the HIMF and the detected galaxy HIMF. If we know enough to trust the strong ALFALFA priors, then the CVID adds important constraints on the HIMF of the optical galaxies and a modest measurement of that of the unresolved galaxies, both of which cannot be obtained from ALFALFA alone. With 10 mK foregrounds, the CVID dramatically improves on the priors. The very brightest end of the HIMF is still dominated by ALFALFA, but the CVID has added a wealth of information about faint galaxies which cannot be obtained conventionally. Unfortunately, even the intensity mapping data loses sensitivity at the very brightest end of the HIMF. Even stronger foreground cleaning would be needed to measure $M_{\rm{min}}$.

These forecasts clearly demonstrate the utility of this method. However, we have made a number of assumptions here that deserve further study. The lognormal form of the galaxy number count PDF is likely overly simplistic, and could be replaced by a more sophisticated prescription Leicht et al. (2019). In Eqs. (9) and (10), we assumed that optical galaxies host the brightest 21 cm emitters. This should hold for blue, gas-rich galaxies, but a number of bright optical galaxies are red and gas-poor Cole et al. (2005), and therefore have weaker HI emission. As mentioned above, we have entirely neglected luminosity-dependent bias in our forecasting. Though the Parkes and GBT surveys are single-dish, many 21 cm experiments are interferometric, which may introduce extra systematics. Finally, we assumed that we can apply some degree of foreground cleaning to our data without affecting our signal. In the Parkes foreground cleaning, the signal was suppressed along with the foregrounds Anderson et al. (2018), which if uncorrected would bias our CVR measurements. These caveats motivate additional study of this method using mock data sets, as in Ref. Ihle et al. (2019).

Though we have focused here on the HIMF, the potential utility of this method extends much further. For example, one could study the HI-color relation mentioned above Wolz et al. (2017). With intensity maps of other lines Kovetz et al. (2019); Crites et al. (2014); Keating et al. (2016); Li et al. (2016); Gong et al. (2017); Croft et al. (2018); The OST mission concept study team (2018); Padmanabhan (2018); Stacey et al. (2018); Cooray et al. (2019); Dumitru et al. (2019), one could, for example, measure the molecular gas properties of Lyman-alpha emitters Chung et al. (2019), or study AGN feedback Breysse and Alexandroff (2019). The CVR can also be modified to combine intensity maps of different lines, another common experimental target Lidz et al. (2011); Breysse and Rahman (2016); Fonseca et al. (2018); Beane et al. (2019).

Cross-correlations have long been a powerful tool for cosmologists, and will only become more critical as more intensity mapping surveys come online. With this work we have demonstrated a one-point cross-correlation method that can be used to clean foregrounds and probe astrophysics inaccessible to conventional surveys. With some refinement, this will be a valuable tool for many future experiments.

The authors would like to thank Ue-Li Pen and Hamsa Padmanabhan for useful conversations. Philippe Berger was supported by Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Copyright 2019. All rights reserved.