A Foreground Masking Strategy for [C ii] Intensity Mapping Experiments Using Galaxies Selected by Stellar Mass and Redshift

simstack is an algorithm that takes galaxy positions from an external catalog, splits them into subsets (typically, but not necessarily, by stellar mass and redshift), and generates mock map layers that are simultaneously regressed with the real-sky map to estimate the mean flux density of each subset. Formally, it is an extension of simple thumbnail stacking (Marsden et al. 2009), the difference being that the off-diagonal entries in the subsets' covariance matrix are not assumed to be zero, so as to account for galaxy clustering. The simultaneous fitting provides a solution to the limitations of stacking in highly confused maps (i.e., biased flux density estimates due to the clustering of sources at angular scales comparable to that of the FIR/sub-mm beam), such that in the theoretical limit where the catalog is complete it naturally leads to a completely unbiased estimator (see Appendix A for some justification). Viero et al. (2013) show that simstack yields unbiased results at any beam size, while conventional thumbnail stacking (e.g., "median" or "mode" stacking, etc.), without additional corrections, inevitably leads to wavelength-dependent biases in the presence of galaxy clustering.

The first step in measuring ${L}_{\mathrm{IR}}({M}_{* },z)$ is to split the catalog into subsets of star-forming and quiescent galaxies based on their U − V versus V − J colors (UVJ, e.g., Williams et al. 2009), and then again into bins of stellar mass (five and three layers for star-forming and quiescent galaxies, respectively) and redshift (eight layers), determined by their optical and near-infrared photometry. We developed an algorithm to calculate the optimized locations of the 5 × 8 + 3 × 8 = 64 stellar-mass/redshift bins so that each bin contains at least 100 (10) star-forming (quiescent) galaxies, as illustrated in Figure 1.

Zoom In Zoom Out Reset image size

Next, the average FIR/sub-mm flux density for each bin and at each wavelength is estimated with simstack. Uncertainties on the mean flux densities are estimated with an extended bootstrap technique which takes into account the uncertainties in the photometric redshift and stellar-mass estimates of individual sources. L_IR for each bin is estimated by first fitting a modified blackbody (or graybody ) with emissivity index β = 2, and the Wien side approximated as a power law with slope α = −2 (Blain et al. 2002), to the full spectrum of intensities $\nu {I}_{\nu }(\lambda )$, and then integrating under the best-fit graybody from ${\lambda }_{\mathrm{rf}}=8$ to 1000 μm. The final step is to fit the full set of mean L_IR with multiple linear regression as described by Viero et al. (2013):

$$\begin{eqnarray}&&\mathrm{log}{L}_{\mathrm{IR}}({M}_{* },z)=\displaystyle \sum _{p=0}^{n}\left\{\left[\displaystyle \sum _{q=0}^{n}{A}_{p,q}{(\mathrm{log}{M}_{* })}^{q}\right]{z}^{p}\right\},\end{eqnarray} \tag{ 1 }$$

where n = 2 and 1 for star-forming and quiescent galaxies, respectively. The coefficient matrices ${A}_{p,q}$ are found to be

$$\begin{eqnarray}{A}_{p,q}^{\mathrm{sf}}=\left(\begin{array}{ccc}2.417 & 0.733 & 0.004\\ -38.84 & 8.080 & -0.406\\ 4.947 & -1.223 & -0.069\end{array}\right)\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}{A}_{p,q}^{\mathrm{qt}}=\left(\begin{array}{cc}0.845 & 0.820\\ 4.556 & -0.354\end{array}\right).\end{eqnarray} \tag{ 3 }$$

Figure 2 shows two sample SED fittings to the stacked fluxes, together with the best-fit polynomials to the mean ${L}_{\mathrm{IR}}({M}_{* },z)$ relations of star-forming and quiescent galaxies separately. As demonstrated in Viero et al. (2012), the modified blackbody approximation produces mean SEDs consistent with best-fit templates such as Chary & Elbaz (2001) and the derived mean L_IR is largely insensitive to the exact choice of the Wien side slope α.

Zoom In Zoom Out Reset image size

At this point, we have modeled the mean infrared luminosity as a function of stellar mass and redshift, but naturally we expect L_IR of individual galaxies to depart from this model, with some characteristic scatter. The question we aim to answer now is: what is the degree of scatter of the full ensemble of sources?

The answer lies in the standard deviation of the residual map (see Table 2 for the definition), which is the difference between the real-sky map at each wavelength and a synthetic map made by applying the ${L}_{\mathrm{IR}}({M}_{* },z)$ model to the original catalog (i.e., the actual stellar masses, redshifts, and sky positions). In a universe where (i) objects are perfectly described by the mean model with no scatter, (ii) catalogs are 100% complete, and (iii) maps have no noise, the residual map would be completely blank. In practice, the actual residual map will have structure due to the intrinsic stochasticity of the galaxy populations, catalog incompleteness, as well as instrumental noise.

Table 2.  Summary of the Terms Used in Our Discussion of Methodology in Sections 2 and 3

Term	Description	Reference
Bins	Stellar mass or redshift intervals used to divide galaxies into sub-populations for stacking analysis.	S2.1
Layer	A subset of a real/mock sky image (or map) attributed to only the sources in the corresponding stellar mass or redshift bin.	S2.1, S2.2
Scatter	In this paper, we exclusively define "scatter" as the standard deviation of flux density or luminosity in the source population, which is characterized and represented by σS in Equation (4).	S2, Appendix A, Equation (4)
(Un)perturbed	Fluxes being assigned to the sources in a specific layer are drawn from a distribution with the mean equal to the best-fit value given by SIMSTACK and some (zero) nonzero width defined by the scatter.	S2.2, Appendix A, Equation (4)
Real/Mock	"Real" refers to the actual sky image, whereas "mock" refers to the image reconstructed using source locations and perturbed mean fluxes from SIMSTACK. More specifically, in our analysis we construct the mock sky image by merging (1) a layer of interest perturbed according to a distribution with a tunable scatter and (2) background layers perturbed by a distribution with a fiducial scatter of 0.3 dex.	S2.2, Equations (5), (6)
Base	The "base" map, different from the mock image, is obtained by merging (1) an unperturbed layer of interest and (2) background layers perturbed by a distribution with a fiducial scatter of 0.3 dex.	S2.2, Equations (5), (6)
Residual	The difference between the real or noise-added mock sky image and a "base" one.	S2.2, Equations (5), (6)
${D}_{\mathrm{real}}$, ${D}_{\mathrm{mock}}$	A small cutout image a few pixels by side, where each pixel measures the standard deviation of a data cube obtained by thumbnail-stacking the residual map at the positions of the sources in each i, j layer.	S2.2, Equations (5), (6)

Download table as:  ASCIITypeset image

We now introduce a method to formally characterize the scatter about the mean ${L}_{\mathrm{IR}}({M}_{* },z)$ relation by leveraging the structure in the residual map. Although our method has similarities with the "scatter stacking" method described in Schreiber et al. (2015), our use of residual maps—estimated by taking the difference with "base" maps generated with simstack-derived luminosities—makes us less susceptible to clustering contamination, and provides a more robust estimate of the scatter in each $({M}_{* },z)$ bin, validated through an extensive set of end-to-end simulations (see Appendix A). Due to the layered structure of our maps, the interplay between individual layers (often with a root-mean-square, or rms, amplitude below the confusion limit of the real map) must be investigated through simulations to estimate the scatter in each layer. For simplicity, we assume that the scatter is dominated by the stochasticity of the star formation activity and therefore is independent of wavelength. We perform the scatter calibration with the 250 μm SPIRE/HerMES (Griffin et al. 2010; Oliver et al. 2012) map which covers the entire COSMOS/UltraVISTA field.

We assume that for a given redshift z_i and stellar mass ${M}_{* ,j}$, the actual flux density S (and therefore the total infrared luminosity) is log-normally distributed about the mean value with a scatter σ_S, an assumption that is motivated by the observed scatter in the star formation main sequence (SFMS, e.g., Sargent et al. 2012). Namely,

$$\begin{eqnarray}&&P(x| \mu ,{\sigma }_{S})=\displaystyle \frac{1}{{\sigma }_{S}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{{(x-\mu )}^{2}}{2{\sigma }_{S}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where $\mu =\mathrm{log}\langle S({M}_{* },z)\rangle$ is the log of the mean flux density measured by simstack. Intuitively, σ_S can be estimated by examining the statistics of source fluxes (e.g., standard deviation) in each stellar-mass/redshift bin of interest. However, for the highly confused far-infrared maps we use, clustering could render measured statistics biased by the contribution from sources in other bins, whose scatter must also be properly accounted for.

Therefore we assign a fiducial scatter to the "background" sources. As will be shown in Section 3.2, the actual scatter can be measured without bias using our method, as long as it is not drastically different from the fiducial value. In particular, the scatter being investigated here is analogous to that of the SFMS, which can be explained as an application of the central limit theorem (Kelson 2014) and is measured to be around ∼0.3 dex (Behroozi et al. 2013; Sparre et al. 2015). We therefore adopt 0.3 dex as the fiducial population scatter for the "background" sources in our mock maps and demonstrate in Section 3 that it is indeed a reasonable choice.

We hereafter refer to the actual sky image as the "real map" and the synthetic map based on the ${L}_{\mathrm{IR}}$ model as the "mock map." In addition, we call a layer unperturbed when the flux density of its sources is constant and equal to the average value μ found using simstack, while a layer is perturbed when each source has been assigned a flux density according to a log-normal distribution of mean μ and scatter σ_S. Finally, as anticipated before, the residual map is either a real or mock map from which the "base" map (the layer of interest, unperturbed, plus the background layers perturbed with the fiducial scatter) is subtracted. This nomenclature is summarized in Table 2.

The crucial step of the algorithm is that we take the standard deviation ${D}_{\mathrm{real}}^{k}$, computed over the positions of all cataloged sources in the $({M}_{* },z)$ bin of interest, of the residual real map, and compare it to its counterpart ${D}_{\mathrm{mock}}^{k}$, which is the standard deviation of a residual mock map obtained by adding up all perturbed layers, plus noise floor (e.g., Nguyen et al. 2010).

Mathematically, at a given pixel of interest k, we have

$$\begin{eqnarray}&&{D}_{\mathrm{real}}^{k}=\mathrm{SD}\left[{S}_{\mathrm{real}}^{k}-\displaystyle \sum _{i,j}{S}_{\mathrm{base}}^{k}({z}_{i},{M}_{* ,j})\right]\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}{D}_{\mathrm{mock}}^{k} & = & \mathrm{SD}\left[\displaystyle \sum _{i,j}{S}_{\mathrm{mock}}^{k}({z}_{i},{M}_{* ,j})+{S}_{\mathrm{noise}}^{k}+\right.\\ & & -\left.\displaystyle \sum _{i,j}{S}_{\mathrm{base}}^{k}({z}_{i},{M}_{* ,j})\right],\end{eqnarray} \tag{ 6 }$$

where SD stands for taking the standard deviation of the thumbnail-stacked cube at each pixel k, and i, j are the indices of redshift and stellar-mass bins (see Table 2 for a reminder of the definitions).

Note that the layer of interest in the mock map is perturbed with different, adjustable levels of σ_S (while all other layers are perturbed with the constant, fiducial value 0.3 dex) to provide a "calibration curve" to compare with the real map. The idea is that the mock map is our best representation of the real map, including the positional source clustering and the scatter in luminosity that may be present in the actual galaxy populations, as well as instrument and confusion noise. From each of these, we want to subtract our best estimate of the average flux density in each layer, i.e., the unperturbed simstack values. At this point, the residual real map will contain, at the positions of the sources in the layer of interest, information on the layer's intrinsic scatter in flux density. The magnitude of this scatter is then simply measured by gauging which level of σ_S in the residual mock map matches the scatter in the residual real map. This is illustrated in Figure 3, where we show how the thumbnail-stacked mock data cube, ${D}_{\mathrm{mock}}^{k}$, compares to the real one, ${D}_{\mathrm{real}}^{k}$, as we tune up the level of the scatter σ_S. For the purpose of measuring the scatter, we focus only on the central pixel of ${D}_{\mathrm{mock}}^{k}$ and ${D}_{\mathrm{real}}^{k}$.

Zoom In Zoom Out Reset image size

As shown in Figure 4, the standard deviation in the mock thumbnail cubes gradually increases with increasing input scatter. The horizontal dashed line represents the standard deviation measured in the real map. The scatter in the real maps can be consequently inferred from the intersection points, marked as squares in the figure. Since the maps are confusion-noise limited, the calibration curves do not start at zero, but rather at some noise floor equivalent to the standard deviation obtained by thumbnail stacking on random, non-source positions. A more detailed justification of this method based on end-to-end simulations is provided in Appendix A.

Zoom In Zoom Out Reset image size

Table 3 lists the results of our extended simstack procedure, i.e., the number of galaxies in each bin, their mean total infrared luminosity and their scatter about the mean. In particular, we find an average logarithmic scatter, of $\langle {\sigma }_{{\rm{L}}}\rangle \sim 0.35$ dex, with no evidence for systematic dependence on redshift or stellar mass, which is consistent with both observations (e.g., Whitaker et al. 2012) and theoretical expectations (e.g., Kelson 2014; Sparre et al. 2015) of the dispersion about the SFMS. Dutton et al. (2010) investigate the origin of such small, roughly constant scatter in the SFMS using a semi-analytic model for disk galaxies based on smooth mass accretion onto dark matter halos and show that the scatter is mainly dominated by the variations in the gas accretion history and therefore does not evolve strongly with time or mass. Note that the method fails to give a reliable estimate of the scatter when the source population's flux density is too close to the noise floor.

Table 3.  Number of Galaxies, Mean Total Infrared Luminosity, and Scatter about the Mean

	Number of Galaxies (${N}_{\mathrm{gal}}$), Luminosity ($\mathrm{log}[{L}_{\mathrm{IR}}/{L}_{\odot }]$), and Scatter (${\sigma }_{{\rm{L}}}\,[\mathrm{dex}]$)
	0 < z < 0.3	0.3 < z < 0.5	0.5 < z < 0.7	0.7 < z < 1a
Star-forming Galaxies
${10}^{10.5}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{13}$	117, ${10.94}_{-0.01}^{+0.01}$, 0.33	296, ${11.04}_{-0.01}^{+0.01}$, 0.34	360, ${11.25}_{-0.01}^{+0.01}$, 0.35	849, ${11.42}_{-0.01}^{+0.01}$, 0.33
${10}^{10.2}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{10.5}$	154, ${10.78}_{-0.01}^{+0.01}$, 0.35	298, ${10.84}_{-0.01}^{+0.01}$, 0.29	338, ${11.11}_{-0.02}^{+0.01}$, 0.34	926, ${11.29}_{-0.02}^{+0.01}$, 0.35
${10}^{10}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{10.2}$	188, ${10.64}_{-0.01}^{+0.01}$, 0.37	367, ${10.68}_{-0.02}^{+0.02}$, 0.33	494, ${10.90}_{-0.02}^{+0.01}$, 0.33	1018, ${11.08}_{-0.01}^{+0.02}$, 0.42
${10}^{9.5}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{10}$ b	691, ${10.28}_{-0.01}^{+0.01}$, 0.44	1095, ${10.43}_{-0.02}^{+0.02}$, 0.43	1561, ${10.61}_{-0.03}^{+0.02}$, ≲0.6	3461, ${10.69}_{-0.02}^{+0.02}$, ≲0.7
Quiescent Galaxies
${10}^{11}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{13}$	89, ${10.08}_{-0.04}^{+0.05}$, ⋯	138, ${10.26}_{-0.06}^{+0.07}$, ⋯	114, ${10.33}_{-0.05}^{+0.05}$, ⋯	255, ${10.36}_{-0.06}^{+0.09}$, ⋯
${10}^{10}\lt \tfrac{{M}_{* }}{{M}_{\odot }}\lt {10}^{11}$	450, ${10.00}_{-0.03}^{+0.03}$, ⋯	774, ${9.85}_{-0.07}^{+0.09}$, ⋯	684, ${9.96}_{-0.09}^{+0.10}$, ⋯	1591, ${10.00}_{-0.08}^{+0.05}$, ⋯

Notes. The scatters of faint, star-forming galaxies in the two low-mass, high-redshift bins are shown as upper limits since in these cases the noise floor (both instrument and confusion) dominates the variance of the residual data cube.

^aRedshift bins are only shown up to z ∼ 1 where the majority of CO foreground comes from. ^bThe lowest-mass layers of faint star-forming and quiescent galaxies are not shown here as their mean and variance are less well-constrained.

Download table as:  ASCIITypeset image

TIME is a high-throughput, millimeter-wave imaging spectrometer array, designed to measure the 3D [C ii] power spectrum. The clustering amplitude constrains the aggregate luminosity of [C ii] emission from EoR galaxies. The instrument parameters of the proposed experiment are summarized in Table 4.

Table 4.  TIME Specifications

TIME Instrument Parameters
Dish size	12 m
Instantaneous FOV	$14^{\prime} \times 0\buildrel{\,\prime}\over{.} 43$
Survey area	$1\buildrel{\circ}\over{.} 3\times 0\buildrel{\,\prime}\over{.} 43$ (1 × 180 beams)
Number of spectrometers	32 (total), 16 per polarization
Spectral range	183–326 GHz
Spectral resolution	90–120
Survey volume	194 Mpc × 1.1 Mpc × 1240 Mpc

Download table as:  ASCIITypeset image

CO emission is derived for each object in the catalog by converting infrared luminosity to CO line strength with the well-established L_IR–${L}_{\mathrm{CO}}^{{\prime} }$ correlation (e.g., Carilli & Walter 2013, hereafter CW13; Greve et al. 2014, hereafter G14; Dessauges-Zavadsky et al. 2015). We can then use the measured stellar mass functions of galaxies at 0 < z < 2 to calculate the power of CO line foregrounds, with the ability of monitoring different subsets (e.g., different stellar mass bins, quiescent versus star-forming galaxies, etc.) Now the total mean intensity of CO contamination can be expressed as

$$\begin{eqnarray}{\bar{I}}_{\mathrm{CO}} & = & \displaystyle \sum _{J}{\displaystyle \int }_{{M}_{* }^{\min }}^{{M}_{* }^{\max }(z)}{{dM}}_{* }{\rm{\Phi }}({M}_{* },z)\\ & & \times \displaystyle \frac{{L}_{\mathrm{CO}}^{J}({L}_{\mathrm{IR}}({M}_{* },z))}{4\pi {D}_{L}^{2}}{y}^{J}(z){D}_{A}^{2},\end{eqnarray} \tag{ 7 }$$

where ${M}_{* }^{\min }=1.0\times {10}^{8}\,{M}_{\odot }$, 3 ≤ J_upp ≤ 6, representing all CO transitions acting as foregrounds and ${\rm{\Phi }}({M}_{* },z)$ being the stellar mass function measured by the COSMOS/UltraVISTA Survey (Muzzin et al. 2013b). ${M}_{* }^{\max }(z)$ represents an evolving mass cut that measures the depth of foreground masking (see Section 4.3 for a detailed discussion) and is set to ${M}_{* ,0}^{\max }=1.0\times {10}^{13}\,{M}_{\odot }$ when no masking is applied. The factor ${y}^{J}{(z)=d\chi /d{\nu }_{\mathrm{obs}}^{J}={\lambda }_{\mathrm{rf}}^{J}(1+z)}^{2}/H(z)$ accounts for the mapping of frequency into distance along the line of sight (Visbal & Loeb 2010). The comoving radial distance χ, the comoving angular diameter distance D_A, and the luminosity distance D_L are related by $\chi ={D}_{A}={D}_{L}/(1+z)$. In the presence of scatter (as is always the case), the expectation value of a function ${ \mathcal F }$ of CO luminosity at a best-fit ${L}_{\mathrm{IR}}({M}_{* },z)$ given by simstack can be written as

$$\begin{eqnarray}&&E[{ \mathcal F }({L}_{\mathrm{CO}}^{J})]={\int }_{-\infty }^{\infty }{dx}{ \mathcal F }({10}^{x})P(x| \mathrm{log}{L}_{\mathrm{IR}}),\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&P(x| \mathrm{log}{L}_{\mathrm{IR}})=\displaystyle \frac{1}{{\sigma }_{\mathrm{tot}}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{{(x-\mu )}^{2}}{2{\sigma }_{\mathrm{tot}}^{2}}\right),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}\mu & = & \mathrm{log}{L}_{\mathrm{CO}}^{J}(\mathrm{log}{L}_{\mathrm{IR}})\\ & = & {\alpha }^{-1}(\mathrm{log}{L}_{\mathrm{IR}}-\beta )+\mathrm{log}{r}_{J}\end{eqnarray} \tag{ 10 }$$

is derived from the best-fit ${L}_{\mathrm{CO}}^{{\prime} }\mbox{--}{L}_{\mathrm{IR}}$ correlation, with r_J being some scaling factor for different J. Consequently, in the presence of scatter, Equation (7) becomes $\langle {\bar{I}}_{\mathrm{CO}}\rangle \equiv E[{\bar{I}}_{\mathrm{CO}}({L}_{\mathrm{CO}}^{J})]$, which describes the expectation value of the total CO mean intensity, averaged over the probability distribution of ${L}_{\mathrm{CO}}^{J}$ as specified by μ and ${\sigma }_{\mathrm{tot}}$. In our calculation, we consider two prescriptions: (i) CW13, who give α = 1.37 ± 0.04, β = −1.74 ± 0.40, and scaling relations appropriate for sub-millimeter galaxies which are used to convert to transitions higher than $J=1\to 0$, and (ii) G14, who provide α and β coefficients for each individual J transition (i.e., r_J = 1) based on samples of low-z ultra-luminous infrared galaxies and high-z dusty star-forming galaxies comparable to CW13. Since the total CO foreground consists of multiple J transitions, we deem the G14 prescription more appropriate for our purposes, because it treats both the slope and intercept as free parameters when fitting to galaxies observed in different J. Henceforth, we present our results based on the G14 model unless otherwise stated.

It is worth noting that the ${L}_{\mathrm{CO}}^{{\prime} }\mbox{--}{L}_{\mathrm{IR}}$ relation is usually determined using a compilation of galaxy samples which collectively spans the stellar mass range ${\mathrm{log}}_{10}({M}_{* })\sim 9.5\mbox{--}11.5$. Simply extrapolating this relation to lower stellar masses without considering possible changes in the ISM in this regime likely overestimates the predicted CO emission, as observations suggest that local galaxies with ${M}_{* }\lt {10}^{9}\,{M}_{\odot }$ are deficient in CO, due to lower molecular gas content and low metallicities (e.g., Bothwell et al. 2014). Dessauges-Zavadsky et al. (2015) find a 0.38 dex scatter in L_IR for a given ${L}_{\mathrm{CO}(1-0)}^{{\prime} }$, which corresponds to 0.32 dex scatter when converting infrared luminosity into CO luminosity. Comparable levels of scatter have also been identified by CW13 and G14 using galaxy samples of similar types. We note that, different from our assumption in Equation (4), Li et al. (2016) re-normalize the log-normal distribution so that the linear mean remains constant and that the level of scatter only affects the shot-noise component of the power spectrum. Instead, we choose to fix the logarithmic mean in this work to best represent the distribution about the best-fit line for the observed $\mathrm{log}{L}_{\mathrm{IR}}\mbox{--}\mathrm{log}{L}_{\mathrm{CO}}^{{\prime} }$ correlation. Also, for simplicity, we ignore any potential correlation between the ∼0.3 dex scatter intrinsic to the total infrared luminosity and the comparable ∼0.3 dex scatter in the infrared-to-CO conversion¹⁰ and combine them orthogonally (i.e., adding in quadrature; see also Li et al. 2016), yielding a total scatter of σ_tot ∼ 0.5 dex, which is what our reference model assumes hereafter.

Figure 5 shows the CO(1 − 0) power spectra predicted by our CO model at z = 1, compared with the best-fit model from Padmanabhan (2018) which is derived from abundance matching the halo mass function to the CO luminosity function observed at 0 < z < 3. The overall power spectrum of the CO foreground can be written as the sum of the clustering and shot-noise terms

$$\begin{eqnarray}&&{P}_{\mathrm{CO}}^{\mathrm{tot}}({z}_{f},{k}_{f})={P}_{\mathrm{CO}}^{\mathrm{clust}}({z}_{f},{k}_{f})+{P}_{\mathrm{CO}}^{\mathrm{shot}}({z}_{f},{k}_{f}),\end{eqnarray} \tag{ 11 }$$

where the clustering component can be derived from the mean intensity I_CO, the average bias ${\bar{b}}_{\mathrm{CO}}(z)$ (Visbal & Loeb 2010) and the nonlinear matter power spectrum ${P}_{\delta \delta }^{\mathrm{nl}}$ (computed with the CAMB-based HMFcalc code; Murray et al. 2013) as

$$\begin{eqnarray}&&{P}_{\mathrm{CO}}^{\mathrm{clust}}({z}_{f},{k}_{f})=\displaystyle \sum _{J}{\bar{b}}_{\mathrm{CO}}^{2}{({\bar{I}}_{\mathrm{CO}}^{J})}^{2}{P}_{\delta \delta }^{\mathrm{nl}}({z}_{f},{k}_{f}),\end{eqnarray} \tag{ 12 }$$

and the shot-noise or Poisson component is given by

$$\begin{eqnarray}{P}_{\mathrm{CO}}^{\mathrm{shot}}({z}_{f}) & = & \displaystyle \sum _{J}{\displaystyle \int }_{{M}_{* }^{\min }}^{{M}_{* }^{\max }(z)}{{dM}}_{* }{\rm{\Phi }}({M}_{* },{z}_{f})\\ & & \times {\left\{\displaystyle \frac{{L}_{\mathrm{CO}}^{J}({L}_{\mathrm{IR}}({M}_{* },{z}_{f}))}{4\pi {D}_{L}^{2}}{y}^{J}({z}_{f}){D}_{A}^{2}\right\}}^{2}.\end{eqnarray} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

Estimating the CO contamination for any given observed [C ii] power spectrum also requires rescaling (i.e., projecting) the corresponding CO comoving power spectrum at low redshift to the redshift of [C ii]. Following Visbal & Loeb (2010) and Gong et al. (2014), the projected CO power spectrum can be written as

$$\begin{eqnarray}&&{P}_{\mathrm{obs},\mathrm{CO}}^{J}({z}_{s},{k}_{s})={P}_{\mathrm{CO}}^{J}({z}_{f},{k}_{f})\times {\left(\displaystyle \frac{\chi ({z}_{s})}{\chi ({z}_{f})}\right)}^{2}\displaystyle \frac{{y}^{J}({z}_{s})}{{y}^{J}({z}_{f})},\end{eqnarray} \tag{ 14 }$$

where $| {{\boldsymbol{k}}}_{f}| =\sqrt{{(\chi {({z}_{s})/\chi ({z}_{f}))}^{2}{k}_{\perp }^{2}+({y}^{J}({z}_{s})/{y}^{J}({z}_{f}))}^{2}{k}_{\parallel }^{2}}$ is the 3D k-vector at the redshift of CO foreground. Here we assume ${k}_{\perp }=\sqrt{{k}_{1}^{2}+{k}_{2}^{2}}$ and ${k}_{1}={k}_{2}={k}_{\parallel }$ for the 3D k-vector $| {{\boldsymbol{k}}}_{s}| =\sqrt{{k}_{\perp }^{2}+{k}_{\parallel }^{2}}$ at the redshift of [C ii] signal.

The left panel of Figure 6 shows our predicted CO power spectra projected into the frame of [C ii] at redshift z = 6.5 or an equivalent observing frequency of ${\nu }_{\mathrm{obs}}=250\,\mathrm{GHz}$. Contributions from different CO transitions (green curves) to the total CO power (gray and black solid curves) are shown by different line styles. For comparison, we also show two alternative CO models from Silva et al. (2015). The simulation-based model (purple line) is derived from fitting to the simulated ${L}_{\mathrm{CO}}^{J}-{M}_{\mathrm{halo}}$ relations (Obreschkow et al. 2009a, 2009b), whereas the observational CO model (red line) is based on rescaling the observed infrared luminosity function (Sargent et al. 2012) with the ratios given by CW13. Finally, we note that Breysse et al. (2015) assume "Model A" of Pullen et al. (2013), which models the CO luminosity at a given dark matter halo mass with a simple scaling relation and predict a much higher level of CO foreground for the [C ii] signal given by the "m2" model of Silva et al. (2015). However, we note that the Pullen et al. (2013) "Model A" is only optimized for observations at z ∼ 2 and fails to capture the transition to sub-linear scaling of the L_CO–M_h relation at halo masses ${M}_{h}\gt {10}^{11}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

The variation in modeling the [C ii] intensity is illustrated by the predicted signals from Gong et al. (2012) and Silva et al. (2015, "m2" model), shown as blue and yellow lines, respectively. Recent ALMA observations of several typical star-forming galaxies at 5 ≲ z ≲ 8 have tentatively suggested a high [C ii]-to-infrared luminosity ratio (Capak et al. 2015; Aravena et al. 2016). The [C ii] luminosity function derived from these observations is similar to that of Gong et al. (2012), implying a high clustering amplitude. In terms of the cumulative number density of z ∼ 6 galaxies, the Gong et al. (2012) model is also supported by recent observations (Aravena et al. 2016; Hayatsu et al. 2017), which suggest a cumulative number density more than 10 times higher for galaxies with ${L}_{[{\rm{C}}{\rm{II}}]}\gt 2\times {10}^{8}\,{L}_{\odot }$ compared to Silva et al. (2015).

3D positional information (R.A., decl., z) from the galaxy catalog allows us to remove spectral–spatial elements (voxels) in the survey. Namely, after a 3D intensity map consisting of ∼8000 voxels has been measured, we discard the voxels contaminated by at least one foreground CO line falling into TIME's spectral range. We specifically use stellar mass as a measure of the masking depth because it is directly provided by the galaxy catalog, and CO power spectra are conveniently parameterized in terms of it. This approach is different from the blind, bright-voxel masking approach (e.g., Breysse et al. 2015), which does not exploit spectral information to identify and mask the voxels contaminated by faint CO sources, and thus fails to reduce the CO foreground sufficiently.

Provided that the catalog is complete between the integration limits (i.e., ${M}_{* }^{\min }\lt {M}_{* }\lt {M}_{* }^{\max }$), it is possible to estimate the loss of survey volume at a given masking depth by simply counting the number of voxels contaminated by the CO lines emitted from galaxies to be masked. Laigle et al. (2016) lists the 90% completeness levels for the COSMOS/UltraVISTA (UltraDeep, or "UD") catalog under consideration here (also shown in Figure 1); the stellar mass limits are ${M}_{* }^{90 \% }\leqslant {10}^{8.9}\,{M}_{\odot }$ for all CO transitions of interest. Since galaxies with ${M}_{* }\leqslant {10}^{8.9}\,{M}_{\odot }$ contribute a negligible fraction (≲0.5%) of the total CO power, the loss fraction is essentially dominated by the choice of masking strategy.

We optimize the masking sequence using an "evolving mass" cut, as shown in Figure 7. Instead of masking galaxies with a simple, universal stellar-mass cut, which results in removing more voxels containing higher-redshift, relatively faint CO-emitters, in order to mask equally massive, lower-redshift, CO-bright counterparts, we define a function ${M}_{* }^{\max }(z)\leqslant {M}_{* ,0}^{\max }=1.0\times {10}^{13}\,{M}_{\odot }$ that is designed to follow a threshold of constant CO(4–3) flux (in W m⁻², assuming G14). Motivated by the range of uncertainty in [C ii] models, we show two examples here corresponding to an extensive masking scheme (Case A) for the Silva et al. (2015) model as well as a moderate masking scheme (Case B) for the Gong et al. (2012) model. Masking essentially reduces the amplitude of CO power spectrum by varying the integration limit of the first and second CO-luminosity moments of the stellar mass function ${\rm{\Phi }}({M}_{* })$, which correspond to the mean intensity and shot-noise power, respectively. Namely,

$$\begin{eqnarray}&&\displaystyle \frac{\langle {\bar{I}}_{\mathrm{CO}}{\rangle }_{{\rm{m}}}}{\langle {\bar{I}}_{\mathrm{CO}}{\rangle }_{\mathrm{um}}}=\displaystyle \frac{E[{\int }_{{M}_{* }^{\min }}^{{M}_{* }^{\max }(z)}{L}_{\mathrm{CO}}({M}_{* }){\rm{\Phi }}({M}_{* }){{dM}}_{* }]}{E[{\int }_{{M}_{* }^{\min }}^{{M}_{* ,0}^{\max }}{L}_{\mathrm{CO}}({M}_{* }){\rm{\Phi }}({M}_{* }){{dM}}_{* }]}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\langle {P}_{\mathrm{CO}}^{\mathrm{shot}}{\rangle }_{{\rm{m}}}}{\langle {P}_{\mathrm{CO}}^{\mathrm{shot}}{\rangle }_{\mathrm{um}}}=\displaystyle \frac{E[{\int }_{{M}_{* }^{\min }}^{{M}_{* }^{\max }(z)}{L}_{\mathrm{CO}}^{2}({M}_{* }){\rm{\Phi }}({M}_{* }){{dM}}_{* }]}{E[{\int }_{{M}_{* }^{\min }}^{{M}_{* ,0}^{\max }}{L}_{\mathrm{CO}}^{2}({M}_{* }){\rm{\Phi }}({M}_{* }){{dM}}_{* }]},\end{eqnarray} \tag{ 16 }$$

where the angle brackets indicate that values are averaged over the log-normal distribution described by Equations (4) and (8).

Zoom In Zoom Out Reset image size

The constant versus evolving stellar-mass cut approaches are explicitly compared in Figure 8, where we show the predicted CO(4–3) power at k = 0.1 h/Mpc (after scale-projecting into the frame of [C ii] at z ∼ 6.5) for the two different masking strategies, and for the two different L_IR–${L}_{\mathrm{CO}}^{{\prime} }$ prescriptions considered in this work, namely CW13 and G14. One can see that there is a clear advantage in using the evolving mass cut strategy, yielding a CO contamination almost an order of magnitude lower than the constant mass cut (at equal masking fractions). We show this masking scheme for TIME voxels individually in Figure 9, where they are positioned according to their spatial (x-axis) and spectral channel (y-axis) indices. For the more extensive Case A masking, of which the depth decreases from $\sim {10}^{9}\,{M}_{\odot }$ at z ∼ 0.3 to $\sim {10}^{10}\,{M}_{\odot }$ at z ∼ 2, about 8% of voxels need to be masked, as indicated by the yellow strips.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the right panel of Figure 6, we show how this masking strategy can effectively bring down the CO contamination to levels that are sub-dominant to the clustering [C ii] power. The power of total CO emission is calculated only from unmasked galaxies with an evolving stellar-mass cut, which results in a ∼8% loss of the total survey volume (Case A). We note that, although our analysis disfavors a total scatter larger than ∼0.5 dex, the uncertainty in the scaling relations of different rotational J transitions among galaxy populations may also affect the predictions of CO power. Such uncertainty can be readily absorbed into the total scatter in our model by examining a broader range of scatter.

The effect of voxel masking on the [C ii] power spectrum is to essentially remove a small fraction of voxels from the survey volume in a nearly random (i.e., uncorrelated) pattern. In Appendix B, we demonstrate using a simulated light cone that simply discarding the CO-contaminated voxels would only cause a change in the raw, measured [C ii] power spectrum of the order of the masked fraction (≲10%), which is already small compared to the expected measurement uncertainty and thus will not affect our predictions for the [CII]/CO power ratio.

In order to obtain the true power spectrum, though, one must correct for the artifact arising from the coupling between Fourier modes due to windowing (i.e., masking) in real space (Hivon et al. 2002; Zemcov et al. 2014). Specifically, individual k modes are propagated through the mask to characterize how their powers are mixed into other modes $k^{\prime}$. A mode-coupling matrix ${M}_{{kk}^{\prime} }$ can be constructed consequently, whose inverse provides the appropriate transformation from a masked power spectrum to an unmasked one. Provided that mode mixing and other systematics such as instrument beam and experimental noise are properly corrected, the [C ii] power spectrum should be measured in an unbiased way in the presence of voxel masking. Alternatively, the correlation information can also be extracted from the two-point correlation function, which is formally the Fourier transform of the power spectrum. It has the advantage of being less affected by the complicated survey geometry and incomplete sky coverage due to masking, albeit making the theoretical interpretation less straightforward. A detailed discussion of such corrections and alternatives is beyond the scope of this paper and thus left for future work.

We show in Figure 10 the evolution of CO power at scale $k=0.1\,h/\mathrm{Mpc}$, where the clustering term dominates, with the masking depth expressed in K-band magnitude ${m}_{{\rm{K}}}^{\mathrm{AB}}$, infrared luminosity L_IR and stellar mass M_*. Two dominant CO transitions (3–2 and 4–3; see Figure 6) are displayed separately here because the conversion between different masking depth expressions is redshift dependent. For our reference model, masking out voxels containing galaxies with ${m}_{{\rm{K}}}^{\mathrm{AB}}\lesssim 22$ at z < 1 renders a total CO power small enough compared with the [C ii] clustering power with a moderate ∼8% loss of total survey volume.

Zoom In Zoom Out Reset image size

The accuracy of masking depends on the error in photometric redshift estimates with respect to instrument spectral resolution; for COSMOS DR2, ${\sigma }_{z}^{\mathrm{phot}}/(1+{z}^{\mathrm{phot}})$ is less than 1% (Laigle et al. 2016), comparable to TIME's typical voxel size in redshift space. While for simplicity the presented masked fractions are calculated assuming the maximum-likelihood photometric redshift, one may perform an even more conservative masking by accounting for the 68% confidence interval of the photometric redshift distribution, which would approximately double the masking fraction. Compared with the uncertainty in masking fraction due to fitting errors in the CO–infrared relation shown in Figure 9, photometric redshift errors would likely dominate the uncertainty in the predicted masking fraction.

As illustrated in Figure 8, we expect to be masking at most ∼700 voxels at 0 < z < 2 to reduce the level of CO contamination to a level required for a solid [C ii] detection; hence, a follow-up campaign to measure spectroscopic redshifts is straightforward, if deemed necessary. For moderate masking (Case B; ∼350 voxels), a typical z ∼ 1 star-forming galaxy close to our masking threshold ${m}_{{\rm{K}}}^{\mathrm{AB}}\sim 21$ requires about three hours of integration to obtain a robust spectroscopic redshift measurement with a multi-object spectrometer like MOSFIRE, which amounts to a total exposure time of about 30 hours for all ∼200 galaxies¹¹ that need to be masked within TIME's survey volume. For the more extensive masking (Case A; down to ${m}_{{\rm{K}}}^{\mathrm{AB}}\sim 22$), spectroscopic confirmation becomes more costly ($\gt 60\,\mathrm{hr}$), so the masking of these fainter sources will be guided solely by photometric redshifts.

Finally, we note that this masking formalism is flexible enough that it can be further optimized in multiple ways. First of all, stacking using more information on the sources (e.g., by including dust extinction, see M. P. Viero et al. 2018, in preparation) than the mass–redshift plane could improve the total infrared luminosity model by reducing the scatter. Moreover, although here the masking depth is chosen quite arbitrarily to roughly trace a constant level of observed CO flux, it can be more formally optimized based on the properties of the foreground emitters, including the level of scatter.

Given the uncertainties in the strength of the [C ii] signal and the CO contamination (see Figure 6), it is desirable to probe the level of remaining CO foreground after the voxel masking technique is applied in order to determine whether the foreground has been removed sufficiently. Silva et al. (2015) discuss the usefulness of cross-correlation as a way to constrain the degree of post-masking foreground. Specifically, cross-correlation can be done either between a foreground CO line and another dark matter tracer (e.g., a known population of galaxies) at the same redshift, or between two foreground CO lines (e.g., $J=4\to 3$ and $J=3\to 2$) emitted from the same redshift but contaminating the intensity maps observed at two different frequencies. The CO-galaxy cross-correlation requires an external data set like COSMOS. The correlation can be checked as the masking depth increases. The CO–CO cross-correlation can be done within the experiment's own data set, albeit at the expense of a potentially lower sensitivity after masking. The cross power in this case serves as a tracer of the degree of contamination as a function of masking depth. Since [C ii] signals from different redshifts are uncorrelated, they do not contribute to the overall cross-correlation power. It is worth noting that these methods can test whether the CO foreground has been removed satisfactorily, although without indicating which sources must be further removed. In Appendix C, we present a more detailed discussion of the usefulness of cross-correlating CO lines from the same redshift, including how it can be used to measure CO lines themselves and thus constrain the cosmic molecular gas content.