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

Intensity mapping provides a unique means to probe the epoch of reionization (EoR), when the neutral intergalactic medium was ionized by the energetic photons emitted from the first galaxies. The [CII] 158$\mu$m fine-structure line is typically one of the brightest emission lines of star-forming galaxies and thus a promising tracer of the global EoR star-formation activity. However, [CII] intensity maps at $6 \lesssim z \lesssim 8$ are contaminated by interloping CO rotational line emission ($3 \leq J_{\rm upp} \leq 6$) from lower-redshift galaxies. Here we present a strategy to remove the foreground contamination in upcoming [CII] intensity mapping experiments, guided by a model of CO emission from foreground galaxies. The model is based on empirical measurements of the mean and scatter of the total infrared luminosities of galaxies at $z<3$ and with stellar masses $M_{*}>10^{8}\,\rm M_{\rm \odot}$ selected in $K$-band from the COSMOS/UltraVISTA survey, which can be converted to CO line strengths. For a mock field of the Tomographic Ionized-carbon Mapping Experiment (TIME), we find that masking out the"voxels"(spectral-spatial elements) containing foreground galaxies identified using an optimized CO flux threshold results in a $z$-dependent criterion $m^{\rm AB}_{\rm K} \lesssim 22$ (or $M_{*} \gtrsim 10^{9} \,\rm M_{\rm \odot}$) at $z<1$ and makes a [CII]/CO$_{\rm tot}$ power ratio of $\gtrsim 10$ at $k=0.1$ $h$/Mpc achievable, at the cost of a moderate $\lesssim 8\%$ loss of total survey volume.


Introduction
C[II] intensity mapping is a powerful tool to study the global impact of star-forming galaxies on the intergalactic medium (IGM) during cosmic reionization (e.g., Gong et al. 2012;Silva et al. 2015;Serra et al. 2016). However, an accurate measurement of the C[II] power spectrum from fluctuations in background intensities relies on the fact that the foreground contamination, which is dominated by CO rotational lines from intermediate redshift galaxies (3 ≤ J upp ≤ 6, at 0 < z < 2), can be appropriately identified and subtracted, or masked.
A variety of foreground removal techniques for general line intensity mapping experiments have been proposed. The masking of pixels (or 3-dimensional "voxels") approach relies on an external galaxy catalog from which positional information and physical properties (notably stellar mass) guide the identification of bright, candidate sources to be masked. Foreground contamination is reduced by masking the survey geometry containing those sources (Gong et al. 2014;Breysse et al. 2015;Silva et al. 2015). Cross-correlating with an alternative tracer of the same cosmic volume has also been widely suggested as a useful way to distinguish the signal and its foreground (Visbal & Loeb 2010;Gong et al. 2012Gong et al. , 2014Silva et al. 2015). Another promising tool of line deconfusion, which makes use of the anisotropy of power spectra projected onto a common frame, is introduced by Visbal & Loeb (2010) and studied in detail recently by Lidz & Taylor (2016) and Cheng et al. (2016). This paper will focus on the optimization of what is inherently the simplest of these approaches: voxel masking. The optimal level at which CO-emitting sources should be masked is an open question, mostly because direct measurements of the CO line emission power spectra themselves is still an emerging field (e.g., Keating et al. 2016). An alternative approach is to use proxies for CO, such as recent star formation (traced by bolometric infrared luminosity) to model the foreground empirically. However, while direct measurements of the mean infrared luminosities of galaxies has been successfully measured by multiple groups (Viero et al. 2013;Schreiber et al. 2015;Viero et al. 2015), the scatter about the mean, which is just as critical for optimal masking, is still uncertain. In this paper we present a novel method to estimate the CO foreground level -complete with mean and scatter -and explore the effect that different levels of masking has on the resulting power spectrum. The method uses bolometric infrared luminosity (L IR ) as the star-formation tracer; models the mean and variance (or scatter) in L IR as a function of stellar mass and redshift; and uses the prescriptions of Carilli & Walter (2013) to convert L IR to CO. The scatter in the model background is estimated by using simulations to calibrate the variance in thumbnail stacks on far-infrared/sub-millimeter residual maps, where a residual map is defined as the difference between the modeled background and the data. We then use the estimated levels of CO foreground to determine the degree of masking necessary to significantly detect the C[II] power spectrum with the Tomographic Ionized-carbon Mapping Experiment (TIME, Crites et al. 2014).
The paper is arranged as follows. In Section 2, we present how we model the mean bolometric infrared luminosity of galaxies as a function of stellar mass and redshift with the simultaneous stacking formalism and algorithm developed by Viero et al. (2013, SIMSTACK 1 ). We also describe in detail the innovative technique of thumbnail stacking on residual maps, used to characterize the scatter. Section 3 is dedicated to the validation of our stacking analysis through simulations. We discuss the observational implications for C[II] intensity mapping experiments in Section 4 and briefly conclude in Section 5. Throughout this paper, we assume a flat, ΛCDM cosmology consistent with the most recent measurement by the Planck Collaboration et al. (2016). Base-10 logarithms are assumed, and all distances are quoted in comoving units unless stated otherwise.

Modeling L IR (M * , z) and its Scatter
The work discussed in this section is an application and an extension of the SIMSTACK method introduced by Viero et al. (2013), which was based on the UDS field. The measurements presented in this work are performed in the COSMOS field (Scoville et al. 2007), where deeper imaging is now publicly available; combining a catalog derived using imaging from Laigle et al. (2016) as processed by the Muzzin et al. (2013b) pipeline, with maps spanning the full far-infrared/sub-millimeter spectral range of the thermal spectral energy distribution (SED) from warmed dust. Note that, in addition to the maps used in Viero et al. (2013), we use maps at 450µm and 850µm from deep SCUBA-2 observations made available by Casey et al. (2013). While the SCUBA-2 maps cover only a small fraction of the COSMOS field (394 arcmin 2 ), they provide critical constraints on the lowenergy end of the SED (for details on the fitting routine, see Moncelsi et al. 2011;Viero et al. 2012). The full data-set is described in detail in Viero et al. (in prep.).
2.1. Estimating the Mean L IR (M * , z) with SIM-STACK SIMSTACK is a simple algorithm which takes positions from an external catalog, split into subsets (typically, but not necessarily, by stellar mass and redshift), and uses them to generate mock map layers that are simultaneously regressed with the real sky map to estimate the mean flux density of each subset. It is formally equivalent to conventional "thumbnail stacking", with the exception that the off-diagonals in the covariance are not excluded. The reason for stressing the term simultaneous is that it provides a solution to the limitation of stacking in highly confused maps, i.e., biased estimates from clustered sources, such that in the theoretical limit where the catalog is complete it naturally leads to a completely unbiased estimator. We refer interested readers to Viero et al. (2013) for further details.
The first step in measuring L IR (M * , z) is to split the catalog into subsets of star-forming and quiescent galaxies based on their U − V vs. V − J colors (UVJ, e.g., Williams et al. 2009), and then again into bins of stellar mass and redshift. Next, average flux densities are estimated for each bin of stellar mass and redshift simultaneously with SIMSTACK, and that is repeated at each wavelength. 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 the individual sources. L IR for each bin is estimated by first fitting a modified blackbody (or graybody) with emissivity index β = 1.8, and the Wien side approximated as a power-law with slope α = −2 (Blain et al. 2002), to the full spectrum of intensities νI ν (λ) (converted from flux density estimates); and then integrating under the best-fit graybody from (rest-frame) λ = 8 to 1000 µm. The final step is to fit the full set of L IR 's with a simple model using linear regression as described in Viero et al. (2013).

Characterization of the Scatter in
To this point, we have described a model for the mean infrared luminosity of galaxies as a function of stellar mass and redshift, but naturally we expect the L IR of individual galaxies to depart from this model, with some characteristic scatter. The question we aim to answer now is by how much; or what is the degree of scatter of the full ensemble of sources?
The key insight is that the answer lies in the variance of the residual map, which is the difference of the actual sky image at each wavelength and a synthetic map made by applying the L IR (M * , z) model to the original catalog (i.e., the actual stellar masses, redshifts, and sky positions). In a universe where the objects are perfectly described by the mean model with no scatter, the residual map would be zero everywhere, but of course the actual residual map will have some structure.
We now introduce a method to formally characterize the scatter about the mean L IR (M * , z) relation by leveraging the structure in the residual map. We note that although our method has similarities with the "scatter stacking" method described in Schreiber et al. (2015), our use of residual maps -estimated from mocks generated with SIMSTACK-derived luminosities -makes us measurably 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-toend simulations (see Section 3). Due to the layered structure of our maps, the potential interactions among layers and their scatter must be investigated through simulations, if a credible measurement of the scatter in each layer is to be performed. For simplicity, we assume that this scatter is independent of wavelength and 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. Viero et al. (2013) show explicitly that SIMSTACK yields unbiased results at any beam size, while alternative techniques (e.g., "median" or "mode" stacking, etc.) with no rigorous mathematical justification may lead to wavelengthdependent biases, in the presence of galaxy clustering.
We assume that for a given redshift z i and stellar mass M * ,j , the actual flux density is log-normally distributed about the mean value with a scatter σ S . Namely, where µ = log S (M , z) is the mean flux density measured by SIMSTACK. Intuitively, we might expect that such a scatter should be measurable for each specific bin by looking at how the standard deviation of source fluxes changes with the value of σ S assigned in a mock map constructed using the positional information from a galaxy catalog, and then compare that with the same measurement for the real map. However, this simple comparison is subject to biases due to the presence of clustered sources in other bins, and must be adequately dealt with. This can be achieved by assigning a fiducial scatter value from the literature to those "background" sources. As we will show explicitly in Section 3.2, the actual scatter can be unbiasedly measured with 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 star formation main sequence (SFMS), which can be understood as a consequence of the central limit theorem (Kelson 2014) and is measured to be around ∼ 0.3 dex from both observations (see Table 5 and 9 in Behroozi et al. 2013 for a compilation of observed values) and numerical simulations (Sparre et al. 2015). We therefore adopt 0.3 dex as the fiducial population scatter for the "background" sources in our mock maps. We hereafter refer to the actual sky image as the "real map" and the synthetic map based on the L 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 . The crucial step of the algorithm is as follows: we take the standard deviation D real,k , computed over the positions of all cataloged sources in the (M * , z) bin of interest, of the real map minus the sum of the unperturbed mock layers, and compare it to its counterpart D mock,k , which is the standard deviation of a mock map obtained by adding up all perturbed layers, plus an instrumental noise floor (e.g., Nguyen et al. 2010), again minus the sum of the unperturbed mock layers.
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), so as 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 scatter that may be present in the actual galaxy populations, as well as instrumental noise. From each of these, we want to subtract our best effort at measuring the average flux density in each layer, i.e. the unperturbed SIMSTACK values. At this point, the two residual maps will contain, at the positions of the sources in the layer of interest, information on the layer's intrinsic dispersion in flux density. The magnitude of this dispersion is then simply measured by gauging which level of σ S in the mock map matches the dispersion in the real map (see Figure 1).
Mathematically, at a given pixel of interest k, we have and As shown in Figure 2, the standard deviation of 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. On average, we find a scatter of ∼ 0.33 ± 0.04 dex, with no obvious systematic dependence on redshift or stellar mass. Note that the maps are confusion-noise limited and thus the calibration curves start at the level of the noise floor rather than zero.

Consistency Check of the Thumbnail Stacking Analysis with Simulated Maps
In order to demonstrate the validity of our scatter characterization through thumbnail stacking analyses, and to identify any potential bias due to clustering, we apply the procedures from the previous section to simulated maps where the input mean and scatter are known. Since C[II] intensity mapping observations are most likely to occur in well-studied fields to facilitate the treatment Fig. 2.-Calibration curves for the scatter in the derived L IR (M * , z) relation at 250 µm. The top and bottom panels correspond to galaxies with stellar mass 10 10 − 10 10.5 M and 10 10.5 − 10 11 M , respectively. The x-axis is the level of input scatter injected into the simulations (labeled "sigma" in the panels of Figure 1), and the y-axis is the measured standard deviation level in the residual maps (illustrated by the color bar in Figure 1). The horizontal, dashed lines are the measured standard deviation levels in the real residual map. The solid curves are the measured output standard deviation levels, with increasing input scatter, of the mock residual maps. The intersecting squares indicate the level of scatter in the real sky images.
of foreground issues, we base our simulations on the COS-MOS catalog, and defer a more thorough analysis involving mock clustered catalogs to future work.
Simulated maps are generated in the same way as the total mock maps for scatter characterization were created, as described in Section 2.2. We note that different realizations of random flux assignments ensure that we obtain distinct maps with similar statistics. Flux densities are assigned to sources in each individual bin {i, j} according to a log-normal distribution, with scatter σ in and mean S in . Viero et al. (2013) demonstrated, using simulations of varying source densities and beam sizes, that flux densities estimated with SIMSTACK are unbiased. Here we Fig. 3.-Robustness of measuring the mean flux densities with SIMSTACK at 250 µm. Open circles represent the mean flux densities before perturbation, and closed circles are after a constant 0.3 dex perturbation. SIMSTACK measurements before and after perturbation are shown as thin and thick lines, respectively, which pass through their respective circles, indicating good agreement with all simulated input values.
further verify that SIMSTACK estimates are robust for subsets with different, nontrivial, levels of scatter. As shown in Figure 3, estimated flux densities S out before and after a 0.3 dex scatter (thin and thick lines, respectively) are consistent with the mean inputs (thin and thick circles, where filled circles are the scattered bins) for a simulated map. We also test that SIMSTACK estimates are unbiased for perturbations of up to 0.5 dex on different combinations of bins, which provides confidence that the mean flux density distribution in the actual ∼ 0.3 dex measurement is correctly estimated.
As a final validation test, we estimate the scatter in simulated maps in which all layers have been collapsed into one image. Figure 4 shows the ratio of calibrated scatter to the input scatter of the simulated map for different subsets (i.e. bins), where the error bars indicate the 68% confidence intervals. Our thumbnail stacking analysis recovers the input scatter to within ∼ 20%. For reference, we also show cases where bins at z ∼ 0.23 and z ∼ 0.75 are perturbed by σ real = 0.2 dex and 0.4 dex, shown respectively as down and up-pointing triangles, to demonstrate that our method is able to correctly identify Fig. 4.-The effectiveness of scatter characterization by thumbnail stacking on the residual maps. Filled markers represent the ratio of the scatter measured by our stacking formalism to the real scatter input (0.3 dex) of the simulated maps. This ratio is evaluated for different stellar mass and redshift bins, as indicated by the x-axis and marker colors. Also shown at z ∼ 0.23 and z ∼ 0.75 are cases where the bin under examination is perturbed by σ real = 0.2 dex and 0.4 dex (down and up-pointing triangles, respectively). Note that the data points are slightly shifted along the x-axis for clarity. the scatter of varying inputs.

Implications for the Masking Strategy of C[II] Intensity Mapping Experiment
We now apply our estimates of the mean and scatter in the L IR (M * , z) relation to model CO emission in the z < 2 sky, for the purpose of guiding foreground masking strategies of future C[II] intensity mapping experiments such as the Tomographic Ionized-carbon Mapping Experiment (TIME, Crites et al. 2014). 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 CO correlation (e.g., Carilli & Walter 2013, Greve et al. 2014, Dessauges-Zavadsky et al. 2015, and 3D positional information (i.e., RA, Dec, z) from the catalog allows us to perform masking at the intra-voxel level. Note that accuracy of the masking will depend on the error in photometric redshift estimates with respect to instrument spectral resolution; for COSMOS DR2, σ z /(1 + z) is less than 1% (Laigle et al. 2016). As illustrated in Figure 7, we expect to be masking less than 1,000 sources at 0 < z < 2 to reduce the level of CO contamination to a level required for a solid C[II] detection; detection; hence, a follow-up campaign to measure spectroscopic redshifts would be trivial, if deemed necessary.
This allows us to look at the masking problem from a different angle; we can use the measured stellar mass functions of galaxies at 0 < z < 2 to calculate the power of CO foregrounds, with the ability of monitoring different subsets (e.g., different stellar mass bins, quiescent vs star-forming galaxies, etc.). Now the total mean intensity of CO contamination can be expressed as where M min * = 1.0 × 10 8 M , y(z) = dχ/dν = λ rf (1 + z) 2 /H(z), 3 ≤ J upp ≤ 6, representing all CO transitions acting as foreground and Φ(M * , z) being the stellar mass function measured by the COSMOS/UltraVISTA Survey (Muzzin et al. 2013a). Note that in the presence of scatter (as is always the case), the expectation of CO luminosity at a best-fit L IR (M * , z) given by SIMSTACK can be written as and where µ = log L J CO (log L IR ) = (1/α)(log L IR − β) + log r J is derived from the best-fit L CO -L IR correlation with α = 1.37 ± 0.04 and β = −1.74 ± 0.40, together with the scaling relations appropriate for sub-millimeter galaxies (Carilli & Walter 2013). Dessauges-Zavadsky et al. (2015) find a 0.38 dex scatter in L IR for a given L CO(1−0) , which corresponds to 0.3 dex scatter when converting infrared luminosity into CO luminosity. Hence, for simplicity, we ignore the correlation between the two sources of scatter 2 and combine them orthogonally (i.e. adding in quadrature), yielding a total scatter of σ tot ∼ 0.5 dex. Consequently, the total foreground CO power spectrum can be expressed written as the sum of the clustering and shot noise terms 2 This is a somewhat arbitrary choice given the potentially similar physics (star formation, dust attenuation, etc.) that leads to the observed scatters in both cases. As it is difficult to accurately determine this correlation, we simply assume here that 0.5 dex is a relatively conservative estimate of the total scatter. where and We note that estimating the CO contamination for any given observed C[II] power spectrum 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 where is the 3D k-vector at the redshift of CO foreground. Here we assume k ⊥ = k 2 1 + k 2 2 and k 1 = k 2 = k for the 3D k-vector |k s | = k 2 ⊥ + k 2 at the redshift of C[II] signal. Figure 5 shows our predicted CO power spectra projected into the frame of C[II] at redshift z = 6.5 or an equivalent observing frequency of ν obs = 250 GHz. A total conservative scatter of 0.5 dex is assumed and contributions from different CO transitions (thin gray curves) to the total CO power (thick black curve) are shown by different line styles. For comparison, we also include two CO alternative models from Silva et al. (2015). The simulation-based model is derived from fitting to the simulated L J CO −M halo relations (Obreschkow et al. 2009a,b), whereas the observational CO model (teal line) is based on rescaling the observed infrared luminosity function (Sargent et al. 2012) with the ratios given by Carilli & Walter (2013). The large variation in modeling the C [II] intensity is on the other hand illustrated by the predicted signals from Gong et al. (2012) and Silva et al. (2015, "m2" model), shown as blue and yellow crosses, respectively. We however note that recent ALMA observations of several typical star-forming galaxies at 5 z 8 have tentatively suggested a high C[II]-to-IR luminosity ratio (Capak et al. 2015;Aravena et al. 2016), favoring a more optimistic scheme similar to Gong et al. (2012).
In the right panel of Figure 5, we demonstrate how masking can effectively bring down the level of CO contamination by showing the power of total CO emission only from unmasked galaxies with stellar mass below a given threshold (M * < 10 10 M and 10 9 M ). For completeness, the power spectra corresponding to lower levels of scatter are also shown by the dashed (0.3 dex) and at ν obs ∼250 GHz (z C[II] ∼ 6.5). Our model illustrates the relative contribution to the total CO power spectrum (thick black line) from different CO transitions (thin gray lines), with a simulated input scatter of σ tot = 0.5 dex. Also shown are two CO models from Silva et al. (2015): the first based on simulated L J CO − M halo relations (purple line labeled 'sim'); and the second from rescaling the observed infrared luminosity function (teal line labeled 'obs'). We note that although our model neatly splits the two Silva et al. (2015) lines, it is by coincidence, not by construction. C[II] signals predicted by Gong et al. (2012) and Silva et al. (2015, "m2" model), shown as blue and yellow crosses, respectively, highlight the large range of existing C[II] predictions. Right: The same comparison as in the left panel, but after masking. The black and red sets of curves denote masking of all galaxies with stellar masses above 10 9 M and 10 10 M , respectively, and z < 2. For illustrative purposes we show the effect that varying σ tot has on the model power spectrum.
dash-dotted (0.1 dex) curves. We note that although our model based on scatter characterization disfavors a total scatter larger than ∼ 0.5 dex, the uncertainty in the scaling relations of different J's among galaxy populations may also affect the predictions of CO power. As such uncertainty can be readily absorbed into the total scatter in our model, we ignore its effect on our discussion of foreground masking by examining a broader range of scatter than what typical characterizations suggest.
Therefore, provided that the catalog is complete between the integration limits (i.e. M min * < M * < M thres * ), 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. Table 6 in Laigle et al. (2016) shows the 90% completeness levels for the COSMOS/UltraVISTA (UltraDeep, or "UD") catalog under consideration here; the mass limits are M 90% * ≤ 10 8.9 M for all CO transitions of interest.
We show in Figure 6 the evolution of CO power at scale k = 0.1 h/Mpc, where the clustering term dominates, with the masking depth expressed in K-band magnitude m AB K , infrared luminosity L IR , stellar mass M * and the masked fraction f mask . Two dominant CO transitions (3-2 and 4-3, see the left panel of Figure 5) are displayed separately here as the conversion between different masking depth expressions is redshift dependent. For our reference model with a total scatter of 0.5 dex, masking out galaxies with stellar mass greater than 10 9 M (or equivalently m AB K 22) will suffice to make the CO power subordinate to the physically interesting C[II] clustering power with a moderate 15% loss of total survey volume.
This masking scheme is shown for TIME voxels indivually in Figure 7, where they are positioned according their spatial (x-axis) and spectral channel (y-axis) indices. For a masking depth of 10 9 M (10 10 M ), about 13% (4%) of voxels need to be masked, whose locations are indicated by the yellow strips. We note that in reality TIME makes a 180 × 1 beam = 1 • × 35 line scan with Fig. 7.-Voxel masking as a method of attenuating the CO foreground in C[II] intensity mapping experiments. The two panels show the masked voxels at masking depths of 10 10 M (left) and 10 9 M (right). By masking all the voxels that are contaminated by CO emission lines (3 ≤ J upp ≤ 6) from low-z galaxies with stellar mass down to 10 9 M , we lose only a moderate fraction ( 15%) of our survey volume. The exact voxels being masked are illustrated in terms of their channel indices (44 spectral and 180 spatial channels) and are calculated from a mock TIME field chosen from COSMOS/UltraVISTA field. Note that the spectral-tospatial aspect ratio of the voxels here is set to 10 for visual clarity, while TIME's will be roughly 30. 44 scientific, spectral channels between ∼ 200 − 300 GHz, which corresponds to an aspect ratio of 30 rather than 10 shown here for visual clarity.
We note that this masking formalism is flexible enough so that the exact masking depths can be optimized for individual subsets in any given actual survey to obtain a desired signal-to-noise. As this paper focuses on the demonstration of foreground removal with galaxies selected by stellar mass and redshfit in general, rather than strategies for specific experiment setups, we will defer a detailed discussion of such optimization frameworks to future work.

Summary
We have presented a novel method to estimate the mean and scatter of CO line emission from measurements of the bolometric infrared luminosity, L IR , and showed how it can be applied as a foreground removal strategy for C[II] intensity mapping experiments. Combining empirical estimates with information about their abundance and spatial distribution provided by COS-MOS survey, we optimized the trade-off between the relative strength of CO/C[II] power being detected and the loss of survey volume. We found that even in the most conservative scenario, by masking galaxies at z < 1 with log(M/M ) ) > 9 -which amounts to K-band magnitudes of m AB 22, or 15% of all voxels, and is already easy to do in current deep and wide-area survey fields such as COSMOS -a C[II]/CO power ratio of greater than 20 achievable. , and CO(4-3) (bottom panel), at redshift z = 6.5. Multiple x-axes are shown to illustrate how the masking depth in AB magnitude projects to L IR , stellar mass M , and mask fraction f mask . Note that the x-axis m AB K and L IR scales differ from top to bottom because the interloping lines in the top and bottom panels originate from different redshifts: z = 0.36 and 0.82 for CO(3-2) and CO(4-3), respectively. Solid horizontal lines represent model predictions from Gong et al. (2012, blue) and Silva et al. (2015, red). Orange solid, dashed, and dash-dotted curves represent the CO power levels vs. masking fraction assuming a scatter of σ tot = 0.1, 0.5, and 1.0 dex, respectively. For a total conservative scatter of 0.5 dex and masking level of log(M/M ) > 9, the power ratio of C[II]/CO(3-2) is between 25 and 300, and of C[II]/CO(4-3) is between 200 and 2000, depending on the comparison model. If we assume the more pessimistic Silva et al. (2015, "m2") C[II] model, and if we include all CO transitions, the total C[II]/CO 3≤Jupp≤6 power ratio is roughly 20.