Intensity Mapping in the Presence of Foregrounds and Correlated Continuum Emission

A major goal of modern astrophysics is to understand the evolution of galaxies in their cosmological context. Two epochal shifts of great interest are the formation of the galaxies that reionize the universe and the subsequent dramatic decline in the star formation rate (SFR) from $z\sim 2$ to the present (Madau & Dickinson 2014) despite the continued growth of dark matter halos. Studies of CO, [C ii], and 21 cm line emission are particularly valuable as tracers of star formation and its precursors in the interstellar medium during these eras (Carilli & Walter 2013).

Most studies of line emission to date measure the properties of individual galaxies to draw broader conclusions about galaxy evolution. These galaxies are often selected from surveys targeting very luminous objects, such as quasars and luminous IR galaxies (Carilli & Walter 2013). Such catalogs provide a biased sample of gas tracers rather than measuring the average population. Even when line emission is detected blindly (e.g., Decarli et al. 2016), these detections may be limited to the brightest objects, so they are not representative of most galaxies. Blind surveys are often in cosmologically small volumes (especially at low redshift), making the source counts subject to cosmic variance (Robertson 2010). Furthermore, surveys for individual objects are expensive because they must detect galaxies at high significance, which requires large apertures or interferometers to gain high flux sensitivity and avoid confusion.

Intensity mapping (IM; Hogan & Rees 1979; Scott & Rees 1990; Madau et al. 1997; Suginohara et al. 1999; Chang et al. 2008; Wyithe et al. 2008; Visbal & Loeb 2010; Visbal et al. 2011) overcomes many of these challenges, so it is attracting increasing interest and investment, as detailed in the recent white paper by Kovetz et al. (2017). Rather than identifying individual objects as in a galaxy redshift survey, IM measures the integrated emission at a given (observed) frequency from galaxies or the intergalactic medium (IGM). It requires only modest aperture sizes to reach the smallest scales of cosmological interest, regardless of confusion. The method was originally developed for 21 cm radiation from reionization (Hogan & Rees 1979; Scott & Rees 1990; Madau et al. 1997) but has expanded to numerous lines and science goals. Observations of CO, [C ii], and 21 cm can constrain the luminosity function of galaxies too faint to observe individually at high redshift. These observations will be critical for constraining the most uncertain aspects of current galaxy formation models, namely, the role of star formation and "feedback" from stars, supernovae, and active black holes.

Several measurements that bear on galaxy evolution have already been provided by IM. Cross-correlation of 21 cm emission and WiggleZ (Switzer et al. 2013) has determined the H i abundance at $z\approx 0.8$, which is consistent with feedback from active galactic nuclei (Villaescusa-Navarro et al. 2016; Padmanabhan et al. 2017). A preliminary indication of [C ii] emission in cross-correlation between Planck and the Baryon Oscillation Spectroscopic Survey (BOSS) quasar sample for $2\lt z\lt 3.2$ (Pullen et al. 2018) favors collisional excitation models at these redshifts that are different from local scaling relations. COPSS-II (Keating et al. 2016) observed CO to provide the first inferences of the molecular gas abundance $2.3\lt z\lt 3.3$, ruling out several models. Croft et al. (2018) found diffuse $\mathrm{Ly}\alpha$ emission in cross-correlation with the $\mathrm{Ly}\alpha$ forest and quasars at a level that is consistent with hydrodynamic simulations. Kovetz et al. (2017) described several dedicated IM experiments, both ongoing and future.

An interpretation of an autonomous IM survey must argue that there is no bias from residual bright continuum radiation (Milky Way and extragalactic), interloper lines from other redshifts, or instrumental response. Cross-correlation with overdensities inferred from a galaxy redshift survey provides an approach to circumvent this foreground contamination. Here, uncorrelated foreground variance increases errors but does not bias the determination. Cross-correlation also allows breakdown by galaxy type or environment to probe the impact of those factors on the interstellar medium (ISM) (Wolz et al. 2017).

We describe several challenges and approaches related to cross-correlation in the specific context of the measurement in Pullen et al. (2018). Here the cross-power of Planck 545 GHz and BOSS quasars is used to constrain [C ii] across $2\lt z\lt 3.2$, and [C ii] is a tracer of great interest for IM. Star formation excites the $157.7\,\mu {\rm{m}}$ (1900 GHz) ${}^{2}{P}_{3/2}\to {}^{2}{P}_{1/2}$ fine-structure transition of [C ii] in several phases of the ISM. Emitting as much as $0.65 \%$ of the typical ${L}_{\mathrm{FIR}}\sim {10}^{12}{L}_{\odot }$ far-IR luminosity of the host (Stacey et al. 2010), [C ii] is the brightest far-IR cooling line of star-forming galaxies. There is a well-established log-linear relation between [C ii] and the SFR (De Looze et al. 2014; Pineda et al. 2014; Herrera-Camus et al. 2015) that becomes more complex at high redshift (Capak et al. 2015; Vallini et al. 2015; Lagache et al. 2018) and as a function of metallicity (Croxall et al. 2017; Smith et al. 2017).

The analysis in Pullen et al. (2018) estimates the [C ii] surface brightness as ${I}_{\nu }={6.6}_{-4.8}^{+5.0}\times {10}^{4}$ Jy sr^–1 but finds inconclusive support for a [C ii] emission component from a Bayesian model test. This is a tantalizing suggestion of [C ii] at $2\lt z\lt 3.2$, and one of our goals here is to attempt to upgrade it to a secure detection or more stringent limit. The BOSS-North region is well out of the galactic plane, but the 545 GHz map does have significant Milky Way dust emission. This contamination suggests that additional effort toward suppressing the foregrounds may be fruitful. With this as a backdrop, we will describe the general problem of rejecting bright Milky Way continuum foregrounds while constraining correlated line and continuum emission.

Typically, IM is described in terms of one isolated line, such as [C ii] above. However, the target galaxies also contribute a characteristic spectral energy distribution (SED), which includes bright dust or a synchrotron continuum contribution at most frequencies of interest for IM. The overdensity field traced by a companion galaxy redshift survey will, therefore, have nonzero cross-correlation with the intensity survey across many frequencies. The cross-correlation of the intensity survey at ${\nu }_{i}$ with galaxies at z_j therefore measures the SED of galaxies at z_j emitting into ${\nu }_{i}$, $\mathrm{SED}({\nu }_{i},{z}_{j})$ (Serra et al. 2014; Switzer 2017; Chiang et al. 2018; Pullen et al. 2018).

A determination of the line brightness must marginalize over the continuum part of the SED. Traditional foreground cleaning that removes uncorrelated continuum emission (such as from the Milky Way) will also suppress the correlated continuum radiation. The cross-correlation in a cleaned map no longer traces $\mathrm{SED}({\nu }_{i},{z}_{j})$ but has residual correlated continuum emission due to the slight differences in SED between uncorrelated (largely Milky Way) and correlated continuum radiation. While foreground cleaning may be successful in revealing the redshifted line emission, it may muddle the correlated continuum and produce bias. We argue that IM analyses should apply a self-consistent approach to deweighting uncorrelated continuum foregrounds and estimating correlated continuum in the target galaxies.

Estimates of the emission of the Milky Way that are independent of extragalactic sources have utility in removing foregrounds without impacting the extragalactic correlated continuum signal. We find that existing spatial templates of reddening to stars (Green et al. 2015) and inferred from Hi (Lenz et al. 2017) do not significantly aid the line amplitude estimation in the case of Planck 545 GHz × BOSS quasars. Similar templates may have great utility in the future for separating continuum terms.

A fundamental challenge in the Planck analysis is the number of available frequency channels relative to the number of estimated parameters. We consider a future tomographic spectral survey with many frequency bands and argue that low ${k}_{\parallel }$ bins isolate the correlated continuum emission (barring instrumental effects). Cutting these bins has a minor impact except in reducing sensitivity at low $| k|$. These conclusions bode well for future spectral IM surveys (Kovetz et al. 2017), which anticipate numerous frequency channels. Throughout, we will use the example of [C ii], but considerations here apply to 21 cm and other targets of IM.

Section 2 reviews Milky Way foregrounds in the Planck data and describes a simple map-space cleaning approach. This approach neglects the joint interpretation of uncorrelated continuum emission from the Milky Way and correlated continuum emission from the target galaxies. Section 3 develops some general results for IM with both correlated and uncorrelated continuum emission in a simplified setting. It argues for the utility of Milky Way–only templates, which we apply to the Planck × BOSS analysis in Section 5. Marginalizing over a correlated continuum model amounts to throwing out those degrees of freedom (dof) along the line of sight. An intensity survey should therefore have many more bands than correlated continuum model parameters. Section 6 considers this true tomographic case and shows that correlated continuum bias lies at low ${k}_{\parallel }$.

Pullen et al. (2018) presented a preliminary indication of [C ii] emission from $2\lt z\lt 3.2$ that is correlated with large-scale structure (LSS) traced by quasars in the BOSS survey. Planck High Frequency Instrument (HFI) data (Lamarre et al. 2010; Planck HFI Core Team et al. 2011) at 545 GHz serve as the intensity map, and measurements at 353 and 857 GHz help to isolate dust continuum and thermal Sunyaev-Zel'dovich (tSZ) emission from the [C ii] of interest. Figure 1 shows significant contamination from the Milky Way in the BOSS-North region of the Planck 545 GHz map. This Milky Way emission is uncorrelated with the [C ii] signal from $2\lt z\lt 3.2$, so it contributes variance but not bias to the cross-correlation with quasars. Figure 2 shows that the [C ii] signal from Pullen et al. (2018) remains $100\times$ lower than the total 545 GHz power spectrum, which is dominated by foregrounds at all scales.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Pullen et al. (2018) performed no foreground cleaning in the map domain, so the bright Milky Way emission contributes significant variance to each of the cross-powers of 353, 545, and 857 GHz with BOSS. Our primary goal here is to assess whether some form of foreground cleaning could upgrade the modest indication of [C ii] in Pullen et al. (2018) to a secure detection of redshifted line emission. We will start with a standard approach of cleaning Milky Way (uncorrelated) continuum emission in map space and argue that this is complicated by correlated dust continuum emission from the target galaxies.

A simple approach to cleaning galactic emission is to form a linear combination of adjacent bands, as

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{\mathrm{clean}}={{\boldsymbol{x}}}_{545}-{\kappa }_{353}{{\boldsymbol{x}}}_{353}-{\kappa }_{857}{{\boldsymbol{x}}}_{857}.\end{eqnarray} \tag{ 1 }$$

The goal of finding linear combination parameters $\{{\kappa }_{353},{\kappa }_{857}\}$ is to minimize the auto-power or rms variation of the cleaned map (and hence the error bars of the cross-correlation). Stack the 353 and 857 GHz maps into a matrix ${\boldsymbol{A}}$. Then the linear combination coefficients that minimize variance are given by

$$\begin{eqnarray}&&{\boldsymbol{\kappa }}={\left({{\boldsymbol{A}}}^{T}{{\boldsymbol{N}}}^{-1}{\boldsymbol{A}}\right)}^{-1}{{\boldsymbol{A}}}^{T}{{\boldsymbol{N}}}^{-1}{{\boldsymbol{x}}}_{545},\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{N}}$ is the covariance of the 545 GHz map. Figure 3 shows the map after linear combination cleaning, with some residual Milky Way dust emission still visible but greatly reduced.

Zoom In Zoom Out Reset image size

Figure 4 shows the cross-power of quasars with the cleaned map compared to the input 545 GHz map.³ The dramatic change in the cross-power from linear combination cleaning motivates a reconsideration of foregrounds in the Pullen et al. (2018) analysis. Not only do the errors diminish, but the level of the cross-power also drops considerably. This is because much of the cross-correlation between Planck 545 GHz and BOSS quasars is from correlated continuum emission.

Zoom In Zoom Out Reset image size

The linear combination that cleans the Milky Way from the 545 GHz map also removes much of the correlated continuum emission. If the dust continuum emission shared the same SED, the linear combination cleaning would be highly effective in removing Milky Way foregrounds and allow a unique determination of the line brightness independent of correlated continuum emission. The fact that the Milky Way dust SED (which drives the linear combination coefficients) is not identical to the correlated continuum SED means that the cleaned cross-power in Figure 4 has some poorly determined residual correlated continuum emission.

The inference of clustered [C ii] has two parts: (1) a determination of the amplitude of [C ii] plus dust continuum (CIB) clustering and (2) a marginalization over the correlated dust continuum to separate the [C ii] contribution. The clustering of [C ii] emission with quasars in $2\lt z\lt 3.2$ is localized to the 545 GHz band, so a linear combination of bands that cleans 545 GHz is an effective way to isolate [C ii]. In contrast, the marginalization over the correlated dust continuum depends on information from 353 and 857 GHz, so it becomes entangled with attempts to use those bands to clean Milky Way emission at 545 GHz.

This section builds up a set of analytic results for multitracer analysis of IM data with correlated and uncorrelated continuum emission contamination. It has become standard to clean continuum contamination in IM data in the map domain and simulate the impact of the resulting signal loss on the power spectrum (e.g., Switzer et al. 2015 for single-dish instruments and Cheng et al. 2018 for interferometers). In contrast, the multitracer approach developed here jointly (1) deweights Milky Way foregrounds, (2) marginalizes over correlated continuum emission, and (3) estimates the line amplitude. In this context, signal loss does not need to be estimated separately because it is included self-consistently in the parameter estimate. Rather than calculating signal loss through simulations, expressions here will focus on the error of the line estimate and how continuum emission causes constraints to degrade from thermal noise limits.

For the expressions to remain tractable and illustrative, we will make several simplifying assumptions: (1) shot noise in the galaxy redshift survey is negligible, and the linear bias is known; (2) inter-bandpower correlations are negligible, so we can consider expressions for a single bandpower/multipole; and (3) there is no stochasticity between the density traced in the intensity and galaxy redshift survey. The line intensity and galaxy fields should trace the same underlying density on linear scales, making this approximation reasonable on large-area surveys such as used in Section 5.

This regime is applicable to single-dish surveys that can access linear modes, such as BINGO (Battye et al. 2016), GBT (Switzer et al. 2013), Parkes (Anderson et al. 2018), and MeerKAT/SKA (Bull et al. 2015) for 21 cm, and most of the proposed experiments for high-frequency lines: CCAT-prime (Vavagiakis et al. 2018), COMAP (Tveit Ihle et al. 2018), CONCERTO (Dumitru et al. 2018), EXCLAIM (Padmanabhan 2018), SPHEREx (Doré et al. 2014), STARFIRE (Uzgil et al. 2014), and HETDEX (Fonseca et al. 2017). Some of these simplifying assumptions can be easily relaxed in this framework. Appendix A describes the impact of stochasticity and shot noise, and Section 4.2 describes low ${k}_{\parallel }$ modes. Additionally, the spherical harmonic approach here is well-matched to surveys on large (curved) regions of the sky.

3.1. Uncorrelated Foregrounds in a Multitracer Setting

Foregrounds that are uncorrelated with cosmological density will increase the errors of a determination of the line brightness from cross-correlation but not bias the result. This section quantifies the impact of uncorrelated foregrounds in a multitracer (McDonald & Seljak 2009) setting, reviewed in Appendix A.

Let ${\boldsymbol{\delta }}$ (a vector) be a spatial map of the cosmological overdensity in a redshift slice that corresponds to a frequency channel in the intensity survey. Model the observed galaxy redshift survey map ${{\boldsymbol{x}}}_{g}$ as a tracer with linear bias b_g and shot noise ${{\boldsymbol{n}}}_{g}$. Model the intensity map as a line surface brightness amplitude S_L times overdensity ${\boldsymbol{\delta }}$ plus noise ${{\boldsymbol{n}}}_{I}$. (Here ${S}_{L}={I}_{L}{b}_{L}$ refers to the line brightness I_L multiplied by the linear cosmological bias b_L of the line emission.) Emission from the Milky Way and galaxies not in the selected redshift slice will contribute continuum emission that is uncorrelated with ${\boldsymbol{\delta }}$. Let the uncorrelated continuum foregrounds have a spatial template ${{\boldsymbol{x}}}_{u}$ and spectral dependence $f(\nu )$. These appear in the intensity map at frequency ${\nu }_{I}$ as $f({\nu }_{I}){{\boldsymbol{x}}}_{u}$. All vectors refer to 2D map slices at constant frequency or redshift. The observed galaxy redshift survey (${{\boldsymbol{x}}}_{g}$) and intensity (${{\boldsymbol{x}}}_{I}$) maps are

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{g} & = & {b}_{g}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{g},\\ {{\boldsymbol{x}}}_{I} & = & {S}_{L}{\boldsymbol{\delta }}+f({\nu }_{I}){{\boldsymbol{x}}}_{u}+{{\boldsymbol{n}}}_{I}.\end{array}\end{eqnarray} \tag{ 3 }$$

Form a data vector ${\boldsymbol{d}}=\{{C}_{{\ell }}^{{gI}},{C}_{{\ell }}^{{gg}}\}$ from the power spectrum of the galaxy redshift survey ${C}_{{\ell }}^{{gg}}$ and the cross-power between the intensity and galaxy redshift survey maps ${C}_{{\ell }}^{{gI}}$. Model the observed two-point functions as

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{{gI}} & = & {b}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta },\\ {C}_{{\ell }}^{{gg}} & = & {b}_{g}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{g},\end{array}\end{eqnarray} \tag{ 4 }$$

where the galaxy bias and shot noise are assumed to be known, for simplicity. The model then has parameters ${\boldsymbol{\theta }}=\{{S}_{L},{C}_{{\ell }}^{\delta \delta }\}$, where ${C}_{{\ell }}^{\delta \delta }$ is the power spectrum of the underlying density ${\boldsymbol{\delta }}$ fluctuations. For illustrative purposes, we consider the constraint based on a single multipole or bandpower. The log-likelihood is

$$\begin{eqnarray}&&2{ \mathcal L }=\mathrm{ln}\,\det {\boldsymbol{\Sigma }}({\boldsymbol{\theta }})+{\left[{\boldsymbol{d}}-{\boldsymbol{\mu }}({\boldsymbol{\theta }})\right]}^{T}{{\boldsymbol{\Sigma }}}^{-1}({\boldsymbol{\theta }})[{\boldsymbol{d}}-{\boldsymbol{\mu }}({\boldsymbol{\theta }})].\end{eqnarray} \tag{ 5 }$$

The terms in the covariance ${\boldsymbol{\Sigma }}$ under the assumption of Gaussian fluctuations in the data are given by

$$\begin{eqnarray}&&\mathrm{Cov}({C}_{{\ell }}^{{AB}},{C}_{{\ell }^{\prime} }^{{CD}})=\displaystyle \frac{1}{{v}_{{\ell }}}{\delta }_{{\ell },{\ell }^{\prime} }({C}_{{\ell }}^{{AD}}{C}_{{\ell }}^{{BC}}+{C}_{{\ell }}^{{AC}}{C}_{{\ell }}^{{BD}}),\end{eqnarray} \tag{ 6 }$$

where ${C}_{{\ell }}^{{AB}}$ is the cross-power of fields A and B, and ${v}_{{\ell }}$ is the number of modes in the measurement, which is roughly $(2{\ell }+1){f}_{\mathrm{sky}}{\rm{\Delta }}{\ell }$ for a 2D map on the fraction ${f}_{\mathrm{sky}}$ of the sky and a bandpower of width ${\rm{\Delta }}{\ell }$. The delta function ${\delta }_{{\ell },{\ell }^{\prime} }$ is exact for all-sky, and for simplicity, we assume that it is approximately the case for broad bandpowers on a partial sky.

Following the approach of Appendix A, the variance of a determination of the line amplitude S_L from a single multipole ℓ is (marginalized over ${C}_{{\ell }}^{\delta \delta }$)

$$\begin{eqnarray}&&{\sigma }_{{S}_{L}}^{2}=\displaystyle \frac{{N}_{I}+f{\left({\nu }_{I}\right)}^{2}{C}_{{\ell }}^{{uu}}}{{C}_{{\ell }}^{\delta \delta }{v}_{{\ell }}},\end{eqnarray} \tag{ 7 }$$

where ${C}_{{\ell }}^{{uu}}$ is the power spectrum of the uncorrelated emission (the map ${{\boldsymbol{x}}}_{u}$), N_I is the noise variance of the intensity map, and ${C}_{{\ell }}^{\delta \delta }$ is the power spectrum of the density fluctuations ${\boldsymbol{\delta }}$. This expression is exact (neglecting shot noise), rather than an expansion about small-intensity map noise (Bernstein & Cai 2011). From Equation (7), foregrounds that are uncorrelated with the cosmological signal play the same role as noise in the intensity map.

Cosmic variance would appear in ${\sigma }_{{S}_{L}}^{2}$ as a term $\propto {S}_{L}^{2}$. As described in Appendix A, this is avoided by adding covariance to the galaxy auto-power ${C}_{{\ell }}^{{gg}}$ in the likelihood (assuming negligible stochasticity between the intensity and galaxy redshift survey).

3.2. Deweighting Uncorrelated Continuum Emission

Intensity measurements at two frequencies can separate uncorrelated continuum and line emission. Include a "veto" map ${{\boldsymbol{x}}}_{V}$ at a different frequency ${\nu }_{V}$. The uncorrelated continuum foreground $f({\nu }_{I}){{\boldsymbol{x}}}_{u}$ in the intensity map will appear at this frequency as $f({\nu }_{V}){{\boldsymbol{x}}}_{u^{\prime} }$, where $f({\nu }_{V})$ accounts for the SED and $u^{\prime}$ in ${{\boldsymbol{x}}}_{u^{\prime} }$ denotes that the spatial pattern of uncorrelated foregrounds in the veto map may not be fully coherent with the spatial pattern ${{\boldsymbol{x}}}_{u}$ in the intensity map. Model this set of maps as

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{g} & = & {b}_{g}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{g},\\ {{\boldsymbol{x}}}_{I} & = & f({\nu }_{I}){{\boldsymbol{x}}}_{u}+{S}_{L}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{I},\\ {{\boldsymbol{x}}}_{V} & = & f({\nu }_{V}){{\boldsymbol{x}}}_{u^{\prime} }+{{\boldsymbol{n}}}_{V}.\end{array}\end{eqnarray} \tag{ 8 }$$

The veto map will also contain line signal from a different redshift, given as ${S}_{L}({z}_{V}){\boldsymbol{\delta }}({z}_{V})$, but we will assume that this has a negligible correlation with ${\boldsymbol{\delta }}({z}_{I})$ in the primary intensity survey. That is an excellent approximation for, e.g., [C ii] emission in Planck 353 GHz versus 545 GHz, which originates from considerably different redshifts. In this case, the redshifted line emission in the veto band amounts to uncorrelated noise, so it is accommodated in the assumed noise power N_V. Two bands are sufficient to fit for a line amplitude and one continuum parameter.

Intensity surveys with greater adjacency in the bands must account for signal correlations in the likelihood. Standard codes (Challinor & Lewis 2011; Di Dio et al. 2013) can calculate interband signal correlations, and several analyses of LSS (Gaztañaga et al. 2012; Di Dio et al. 2014; Nicola et al. 2014; Salazar-Albornoz et al. 2017) have used similar approaches that could be adopted for IM. Appendices B and C describe the multiband case, which can accommodate more general SEDs of the continuum correlations and interband signal correlations. Sections 4.2 and 6 describe low-${k}_{\parallel }$ modes in the ${C}_{{\ell }}$ and $P({k}_{\parallel },{k}_{\perp })$ context, respectively.

To characterize the lack of coherence between the spatial distribution of uncorrelated foreground contamination in the intensity map (${{\boldsymbol{x}}}_{u}$) and veto map (${{\boldsymbol{x}}}_{u^{\prime} }$), define "foreground" stochasticity ${r}_{F}={C}_{{\ell }}^{{uu}^{\prime} }/\sqrt{{C}_{{\ell }}^{{uu}}{C}_{{\ell }}^{u^{\prime} u^{\prime} }}$, where ${C}_{{\ell }}^{{uu}^{\prime} }$ is the cross-power between ${{\boldsymbol{x}}}_{u}$ and ${{\boldsymbol{x}}}_{u^{\prime} }$.

Expand the observation vector in Equation (5) to include the cross-power of the veto map and galaxy redshift survey ${\boldsymbol{d}}=\{{C}_{{\ell }}^{{gI}},{C}_{{\ell }}^{{gV}},{C}_{{\ell }}^{{gg}}\}$ and estimate the parameters ${\boldsymbol{\theta }}=\{{S}_{L},{C}_{{\ell }}^{\delta \delta }\}$. In the limit that the veto band measures the foregrounds well, e.g., ${N}_{V}\ll f{({\nu }_{V})}^{2}{C}_{{\ell }}^{{uu}}$. Then,

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{S}_{L}}^{2} & = & \displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left[{N}_{I}+\displaystyle \frac{f{\left({\nu }_{I}\right)}^{2}}{f{\left({\nu }_{V}\right)}^{2}}{N}_{V}\right]\\ & & +\ \displaystyle \frac{f{\left({\nu }_{I}\right)}^{2}{C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}(1-{r}_{F}^{2}).\end{array}\end{eqnarray} \tag{ 9 }$$

To interpret this result, make a cleaned map that removes uncorrelated continuum emission using the veto map, as ${{\boldsymbol{x}}}_{\mathrm{clean}}={{\boldsymbol{x}}}_{I}-[f({\nu }_{I})/f({\nu }_{V})]{{\boldsymbol{x}}}_{V}$. In the limit that the foreground in the veto band perfectly traces the intensity map (r_F = 1), ${{\boldsymbol{x}}}_{\mathrm{clean}}={S}_{L}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{I}-[f({\nu }_{I})/f({\nu }_{V})]{{\boldsymbol{n}}}_{V}$. The noise variance in this cleaned map is ${N}_{I}+{N}_{V}[f{({\nu }_{I})}^{2}/f{({\nu }_{V})}^{2}]$, which is the numerator of the first term in Equation (9). The first term therefore describes the impact of thermal noise in a cleaned map. To interpret the second term, note that the numerator is $f{({\nu }_{I})}^{2}{C}_{{\ell }}^{{uu}}$, the foreground covariance in the intensity map. If r_F = 0, the veto band provides no cleaning in the intensity band, and the variance is increased by the full foreground brightness analogous to Equation (7). Alternately, if r_F = 1, the veto band is a perfect tracer of uncorrelated foregrounds in the intensity map, and the second term goes to zero. In this limit, the line brightness can be determined without impact from the foreground variance.

We have not explicitly cleaned the map anywhere in the likelihood. The appearance of an underlying cleaned map ${{\boldsymbol{x}}}_{\mathrm{clean}}$ arises self-consistently with the parameter estimate through the action ${{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{d}}$ in the likelihood, Equation (5). Foregrounds common to both the intensity and veto maps show up as correlated noise between the cross-powers ${C}_{{\ell }}^{{gI}}$ and ${C}_{{\ell }}^{{gV}}$. Specifically, these off-diagonal terms in ${\boldsymbol{\Sigma }}$ appear from the cross-variance ${C}_{{\ell }}^{{IV}}=f({\nu }_{I})f({\nu }_{V}){C}_{{\ell }}^{{uu}^{\prime} }$. For a simple demonstration in the limit that r_F = 1 and there is no noise in the veto map (N_V = 0), ${{\boldsymbol{\Sigma }}}^{-1}$ takes a linear combination of the cross-powers ${C}_{{\ell }}^{{gI}}$ and ${C}_{{\ell }}^{{gV}}$ to form a cleaned cross-power analogous to the cleaned map, as

$$\begin{eqnarray}&&{C}_{{\ell }}^{{gI}}{| }_{\mathrm{clean}}={C}_{{\ell }}^{{gI}}-[f({\nu }_{I})/f({\nu }_{V})]{C}_{{\ell }}^{{gV}}.\end{eqnarray} \tag{ 10 }$$

An additional map cleaned through a linear combination in map space will not add information to the likelihood. Furthermore, taking a linear combination such as Equation (1) outside of the context of the likelihood can induce bias from, e.g., spurious correlation of foregrounds and signal. This necessitates signal loss simulations to account for the impact of any map operations done before the parameter likelihood in power-spectrum space (Switzer et al. 2015).

In summary, bright Milky Way continuum foregrounds can be deweighted through observations at multiple frequencies, limited by the coherence r_F between frequencies.

3.3. Uncorrelated and Correlated Continuum Emission

In addition to line emission, the host galaxies also emit a dust continuum that correlates with cosmological overdensity. To determine the line amplitude independently of the continuum in these galaxies, the continuum contribution must be modeled and marginalized over. Model this correlated continuum component as $g(\nu ){S}_{C}{\boldsymbol{\delta }}$, where $g(\nu )$ describes the SED of the correlated continuum and S_C is the amplitude of the continuum. Extending the map model in Equation (8),

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{g} & = & {b}_{g}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{g},\\ {{\boldsymbol{x}}}_{I} & = & f({\nu }_{I}){{\boldsymbol{x}}}_{u}+g({\nu }_{I}){S}_{C}{\boldsymbol{\delta }}+{S}_{L}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{I},\\ {{\boldsymbol{x}}}_{V} & = & f({\nu }_{V}){{\boldsymbol{x}}}_{u^{\prime} }+g({\nu }_{V}){S}_{C}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{V}.\end{array}\end{eqnarray} \tag{ 11 }$$

Continue to use the observation vector ${\boldsymbol{d}}=\{{C}_{{\ell }}^{{gI}},{C}_{{\ell }}^{{gV}},{C}_{{\ell }}^{{gg}}\}$ but expand the parameters to include an estimate of the correlated continuum amplitude S_C, as ${\boldsymbol{\theta }}=\{{S}_{L},{S}_{C},{C}_{{\ell }}^{\delta \delta }\}$. The variance of the estimate of S_L not marginalized over the other parameters (${({F}_{{S}_{L},{S}_{L}})}^{-1}$ for Fisher matrix F) is the same as Equation (9). That is, we can determine the amplitude of ${\boldsymbol{\delta }}$ in the intensity map as before. However, now the correlated amplitude is ${S}_{L}+g({\nu }_{I}){S}_{C}$. To determine the line brightness S_L independently of the correlated continuum $g({\nu }_{I}){S}_{C}$, we need to marginalize over S_C, giving

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{S}_{L}}^{2} & = & \displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left[{N}_{I}+\displaystyle \frac{g{\left({\nu }_{I}\right)}^{2}}{g{\left({\nu }_{V}\right)}^{2}}{N}_{V}\right]\\ & & +\ \displaystyle \frac{{C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}{\left[f({\nu }_{I})-f({\nu }_{V})\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}\\ & & +\ \displaystyle \frac{2{C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}(1-{r}_{F})f({\nu }_{I})f({\nu }_{V})\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}.\end{array}\end{eqnarray} \tag{ 12 }$$

To interpret this result, make a cleaned map that removes the correlated continuum emission using the veto map, as ${{\boldsymbol{x}}}_{\mathrm{clean}}={{\boldsymbol{x}}}_{I}-[g({\nu }_{I})/g({\nu }_{V})]{{\boldsymbol{x}}}_{V}$. Unlike the previous section, this map cleans the correlated continuum emission using a linear combination of bands. The first term in Equation (12) is just the thermal noise in this cleaned map. The term in brackets in the second line is the residual of the uncorrelated continuum foreground emission after the correlated emission has been cleaned. If both the correlated and uncorrelated continua have the same SED, $g(\nu )\propto f(\nu )$, both are deweighted equally. The third term is a cross-term and arises from the stochasticity between the spatial shape of uncorrelated continuum foregrounds in maps ${{\boldsymbol{x}}}_{I}$ and ${{\boldsymbol{x}}}_{V}$, which contributes variance to ${\sigma }_{{S}_{L}}^{2}$ even if the correlated and uncorrelated continua have the same SED.

3.4. Bias from Incorrect Correlated Continuum Models

If the correlated continuum spectrum is not accurately modeled, residuals may be spuriously interpreted as line emission. In the linear model for two intensity bands (Equation (11)), let the true correlated continuum SED be $g^{\prime} (\nu )$ and model it as $g(\nu )$. In this case, the estimated line amplitude is

$$\begin{eqnarray}&&{\hat{S}}_{L}={S}_{L}+{S}_{C}\left[g^{\prime} ({\nu }_{I})-\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}g^{\prime} ({\nu }_{V})\right],\end{eqnarray} \tag{ 13 }$$

which is biased by the difference between the true SED and the model.

In the case where the number of instrument bands ${N}_{\mathrm{band}}$ equals the number of parameters ${N}_{\mathrm{param}}$, the spectral model uses all dof, and there is no additional handle to identify bias in the line amplitude. The model may have an excellent goodness of fit because it is using S_L to fit a residual from an insufficient continuum model. If ${N}_{\mathrm{band}}\gg {N}_{\mathrm{param}}$, channels that have no expected line emission correlation should be consistent with zero, providing a test.

In summary, in the case of relatively few bands, the line amplitude may be undetectably compromised by an incomplete continuum SED model and subject to variance from bright Milky Way emission that cannot be independently downweighted. Section 4 develops a model to forecast an experiment with many bands in a ${C}_{{\ell }}$ approach. Section 6 describes correlated continuum emission in the context of $P({k}_{\perp },{k}_{\parallel })$ in a survey with many bands.

The preceding toy models give some intuition for the impact of uncorrelated and correlated continuum emission. The analytic Fisher matrix approach presented there becomes cumbersome to treat more than two intensity channels, requiring analytic ${N}_{\mathrm{band}}\times {N}_{\mathrm{band}}$ inverses for ${N}_{\mathrm{band}}$ maps for each parameter, followed by an analytic ${N}_{\mathrm{param}}\times {N}_{\mathrm{param}}$ inverse over parameters to find the marginalized error on S_L.

Appendix C builds a linear model for the joint deweighting of uncorrelated continuum and fit to both correlated line and continuum emission. Here we model the cross-power at each ℓ as a set of correlated continuum and line parameters ${\boldsymbol{\theta }}$ that linearly describe the measured cross-powers as a function of frequency ${{\boldsymbol{d}}}_{{\ell }}$ as ${{\boldsymbol{d}}}_{{\ell }}={\boldsymbol{M}}{\boldsymbol{\theta }}+{{\boldsymbol{n}}}_{\times }$, where the bandpowers have noise ${{\boldsymbol{n}}}_{\times }$ described by covariance ${{\boldsymbol{N}}}_{\times }$. The parameters ${\boldsymbol{\theta }}=\{{S}_{L},{{\boldsymbol{S}}}_{C}\}$, where S_L is the line amplitude and ${{\boldsymbol{S}}}_{C}{| }_{j}={S}_{C,j}$ is the amplitude of linear spectral templates ${g}_{j}({\nu }_{i})$ that describe the correlated continuum emission as a function of frequency ${\nu }_{i}$. The matrix ${\boldsymbol{M}}\equiv [{\boldsymbol{\xi }},{{\boldsymbol{g}}}_{1},\; ...{{\boldsymbol{g}}}_{{N}_{\mathrm{comp}}}]$ holds the spectral templates for the line correlation ${\boldsymbol{\xi }}$ and ${N}_{\mathrm{comp}}$ component spectral modes of the correlated continuum, and ${{\boldsymbol{\xi }}}_{i}=\xi ({\nu }_{i})$ describes the correlation of the line intensity at each frequency ${\nu }_{i}$ with the galaxies at redshift z. In the case of Planck, the bands are widely spaced in redshift so that $\xi ({\nu }_{i})$ is a delta function at ${\nu }_{i}=545\,\mathrm{GHz}$ to an excellent approximation. Section 4.2 describes the impact of correlations at low ${k}_{\parallel }$ on an intensity survey with many narrow bands. Appendix B describes how the Milky Way produces strong, low-rank correlations in ${{\boldsymbol{N}}}_{\times }^{-1}$, and Appendix D describes the case of foregrounds that have the same spatial template at all frequencies.

The variance ${{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}$ on the correlated line and continuum amplitude parameters ${\boldsymbol{\theta }}$ in this simple linear model is

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}={\left({{\boldsymbol{M}}}^{T}{{\boldsymbol{N}}}_{\times }^{-1}{\boldsymbol{M}}\right)}^{-1}.\end{eqnarray} \tag{ 14 }$$

Appendix C shows that this simple form reproduces the results from previous sections under a slightly more restrictive set of assumptions.

4.1. How Many Bands Are Needed?

The linear formulation in Equation (14) emphasizes an accounting for overall dof. Let the line signal have one dof (the line amplitude S_L), the correlated continuum model have ${D}_{\mathrm{corr}}$ dof, and the uncorrelated continuum emission model have ${D}_{\mathrm{uncorr}}$ dof. Assuming these modes are mutually independent, the total number of dof needed to account for all terms is ${D}_{\mathrm{corr}}+{D}_{\mathrm{uncorr}}+1$. If a survey has ${D}_{\mathrm{corr}}+1$ bands, it can separate correlated line and continuum emission, but it has no dof remaining to suppress uncorrelated (Milky Way) emission. In the model of Equation (11), there are two intensity map channels and one dof for the line and correlated continuum emission, respectively. This leaves no additional freedom to remove uncorrelated foregrounds, and these show up in the second term of Equation (12).

Equation (14) permits a rapid simulation of an instrument with ${N}_{\mathrm{band}}$ uniform bands for a hypothetical [C ii] survey from $2\lt z\lt 2.5$ (543 and 633 GHz) to illustrate dof. For correlated emission, we fit an unknown amplitude for a ${T}_{\mathrm{dust}}=27.2$ K, β = 1.5 dust source at z = 2.25 (${I}_{\nu }^{\mathrm{dust}}={\nu }^{\beta }B[{T}_{\mathrm{dust}},(1+z)\nu ]$, where $B({T}_{\mathrm{dust}},\nu )$ is the Planck law) and for the amplitude of the derivative of intensity with respect to temperature ${{dI}}_{\nu }^{\mathrm{dust}}/{dT}$, corresponding to a fit to the dust temperature to first order. There are two correlated dof to fit the continuum and one signal dof, which gives three dof in total. For Milky Way dust emission, we take emission at 19 K and β = 1.5 with an amplitude that is $10\times$ the thermal noise (after marginalizing over the correlated continuum). We take the signal correlation to exist only in one spectral band, but Section 4.2 describes the impact of signal correlations along the line of sight. Figure 5 shows the inflation of errors (relative to thermal) as a function of the number of bands. In a survey with three bands, there is no remaining latitude to deweight uncorrelated continuum, and the error on the line amplitude is penalized by $10\times$ thermal, consistent with the variance of the contamination. Adding one additional channel provides a dof to suppress the Milky Way, independent of the correlated continuum determination.

Zoom In Zoom Out Reset image size

The simulation in this section is set up to emphasize the impact of dof in a general survey. Section 5.6 uses Equation (14) to interpret the Planck × BOSS and shows that there is not a strong penalty for the three bands in this particular case, given the measured amplitude of the Milky Way foreground and similarity in the correlated and uncorrelated spectral indices.

4.2. Line-of-sight Correlations

Cosmological structure is also correlated along the line of sight, so a slice of density from the galaxy redshift survey at redshift z will correlate not only with line emission from redshift z but also nearby redshifts $z^{\prime}$. Correlations at low ${k}_{\parallel }$ (and consequently at low ℓ) can extend over many spectral bins in the intensity map. In the case of Planck, this effect is negligible between the bands 353, 545, and 857 GHz (at the ℓ reported here), but the effect will be significant for future tomographic surveys with many narrow bands.

Section 6 develops a complete view of this effect by considering the 3D power spectrum of a tomographic survey. It shows that signal at low ${k}_{\parallel }$ is corrupted by correlated continuum emission. However, the linear framework in Equation (14) provides another approach to understand the impact of this effect. Appendix C describes line-of-sight signal correlations in terms of a correlation kernel $\xi ({\nu }_{i})$ between the line intensity at ${\nu }_{i}$ and a galaxy survey at fixed redshift z. This kernel can be calculated as described in Challinor & Lewis (2011) and Di Dio et al. (2013).

Figure 6 demonstrates the inflation of errors at low ${k}_{\parallel }$ in the same simulation setup as Section 4.1 and taking 200 channels. At low ${k}_{\parallel }$, the clustered line emission is spectrally smooth. It is, therefore, less distinguishable from correlated and uncorrelated continuum emission, so its errors rapidly grow. Here modes with spatial wavenumbers ${k}_{\parallel }$ $10\times$ higher than the largest modes in the volume gain immunity by varying rapidly spectrally compared to the continuum.

Zoom In Zoom Out Reset image size

It has been recognized since Hogan & Rees (1979) that bright and spectrally smooth Milky Way emission contaminates long ${k}_{\parallel }$ modes. A new point emphasized here is that correlated continuum emission also contaminates these modes.

4.3. Application to Auto-power Analysis

4.4. Adding a Template Map for the Milky Way

The measurement of Planck × BOSS quasars is in the challenging regime where ${N}_{\mathrm{band}}\approx {N}_{\mathrm{param}}$. Given that this is the best intensity map for [C ii] currently available, we can try to improve the foreground cleaning with a Milky Way template.

Extend the map model in Equation (11) to include a Milky Way template with amplitude A_M and noise variance N_M, as

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{g} & = & {b}_{g}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{g},\\ {{\boldsymbol{x}}}_{I} & = & f({\nu }_{I}){{\boldsymbol{x}}}_{u}+g({\nu }_{I}){S}_{C}{\boldsymbol{\delta }}+{S}_{L}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{I},\\ {{\boldsymbol{x}}}_{V} & = & f({\nu }_{V}){{\boldsymbol{x}}}_{u}+g({\nu }_{V}){S}_{C}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{V},\\ {{\boldsymbol{x}}}_{M} & = & {A}_{M}{{\boldsymbol{x}}}_{u}+{{\boldsymbol{n}}}_{M}.\end{array}\end{eqnarray} \tag{ 15 }$$

Equation (14) gives

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{S}_{L}}^{2} & = & \displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left[{N}_{I}+{\left[\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}{N}_{V}\right]\\ & & +\ \displaystyle \frac{{N}_{M}{A}_{M}^{-2}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}{\left[f({\nu }_{I})-f({\nu }_{V})\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}.\end{array}\end{eqnarray} \tag{ 16 }$$

This is just Equation (12) where the bright Milky Way variance ${C}_{{\ell }}^{{uu}}$ has been traded with signal-to-noise ratio in the tracer map, ${N}_{M}/{A}_{M}^{2}$. Appendix E includes the impact of stochasticity r_M between the tracer of the Milky Way emission and the continuum emission in the observed bands. In this case, an additional term of the variance is

$$\begin{eqnarray}&&{\sigma }_{{S}_{L}}^{2}{| }_{\mathrm{stoc}}=\displaystyle \frac{{C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}{\left[f({\nu }_{I})-f({\nu }_{V})\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}(1-{r}_{M}^{2}).\end{eqnarray} \tag{ 17 }$$

If the Milky Way template traces the dust emission poorly, ${r}_{M}\to 0$, this reverts to the full foreground contamination in Equation (12). Appendix D describes the extension to many bands.

The previous section argued that Milky Way templates can provide leverage to suppress foregrounds, especially in the case where the number of spectral channels is similar to the number of correlated continuum parameters. The analysis of Planck × BOSS for [C ii] in Pullen et al. (2018) is in this regime. This section considers the application of Milky Way templates in this case to boost the signal-to-noise ratio of the correlated line emission.

5.1. Existing Milky Way Templates

There are approaches to using broadband measurements to model thermal dust emission in the Milky Way (Planck Collaboration et al. 2014a, 2016b; Meisner & Finkbeiner 2015). However, these galactic templates are prone to also containing extragalactic thermal dust emission, including correlated continuum emission. We instead seek models of thermal dust emission that do not rely on sensitivity to dust emission directly.

Thermal dust emission depends on both the column depth and the temperature. Existing inferences of temperature measure dust emission around the peak of the SED using broadband observations. These maps may, therefore, correlate with extragalactic radiation and cannot be used. Both reddening and H i abundance trace the dust column depth and can be inferred independently of extragalactic emission.

We use the direct inference of reddening to Pan-STARRS stars from Green et al. (2015), shown in Figure 7. This map is noisy, capturing only the brightest features in the BOSS-North field. To reach higher sensitivity, we also consider the reddening map from Lenz et al. (2017) inferred from HI4PI (HI4PI Collaboration et al. 2016) measurements of the H i column through 21 cm radiation (Figure 8). Both Green et al. (2015) and Lenz et al. (2017) are imperfect tracers of Milky Way dust emission. The H i tracers for dust density become biased in high-density regions where ${{\rm{H}}}_{2}$ has formed (Lenz et al. 2017), so the H i tracer is best in areas out of the galactic plane considered here. Reddening and H i are imperfect tracers of the dust column depth (Cardelli et al. 1989; Fitzpatrick & Massa 2007; Schlafly et al. 2016), and the dust column is an imperfect tracer of the thermal dust emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the reduction in bandpower errors when the Green et al. (2015) and Lenz et al. (2017) maps are projected out of the Planck 545 GHz channel. This linear combination in map space is only for visual demonstration. The Milky Way emission template from Lenz et al. (2017) decreases variance on the largest scales in the measured cross-power, while the template from Green et al. (2015) produces higher noise in the cross-power for ${\ell }\gt 300$.

Zoom In Zoom Out Reset image size

5.2. Planck and BOSS Data

5.3. Measured Cross-powers

Figure 10 shows the six angular cross-powers used to estimate [C ii] emission and continuum nuisance parameters. These cross-spectra were performed over three Planck bands: 353, 545, and 857 GHz. These maps were cross-correlated with overdensity maps of quasars and LRGs cataloged by the BOSS survey. We perform the cross-correlations over angular scales ${\ell }=100\mbox{--}1000$. These spectra are identical to those presented in Pullen et al. (2018).

Zoom In Zoom Out Reset image size

We also consider two Milky Way templates to marginalize over contamination in the Planck maps. Specifically, we use the Hi–derived dust map from Lenz et al. (2017) and the reddening map from Green et al. (2015). As observables, we cross-correlate these dust maps with both the quasar and LRG overdensity maps (Figure 11). Section 5.4 investigates evidence for correlation with LSS.

Zoom In Zoom Out Reset image size

5.4. Models for CIB and C ii Signals

We briefly review the CIB and [C ii] emission models used in the [C ii] constraints from Pullen et al. (2018). We begin with the angular cross-power spectrum between CIB emission in a Planck map at frequency ν and an LSS tracer overdensity field, an expression given by

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{\mathrm{LSS}-\mathrm{CIB}} & = & \int \displaystyle \frac{{dz}}{\chi {\left(z\right)}^{2}}{\left(\displaystyle \frac{d\chi }{{dz}}\right)}^{-1}{b}_{\mathrm{LSS}}{b}_{\mathrm{CIB}}(k,z)\\ & & \times \displaystyle \frac{{dN}}{{dz}}(z)\displaystyle \frac{{dS}}{{dz}}(z,\nu ){P}_{\mathrm{DM}}\left(\displaystyle \frac{{\ell }}{\chi (z)},z\right),\end{array}\end{eqnarray} \tag{ 18 }$$

where χ is the comoving distance, dN/dz is the redshift distribution of the LSS tracer, ${P}_{\mathrm{DM}}(k,z)$ is the matter power spectrum computed using CAMB, and ${b}_{\mathrm{LSS}}$ and ${b}_{\mathrm{CIB}}(k,z)$ are the clustering biases for the LSS tracer and CIB emitters, respectively. For the LSS tracers we consider, we fix ${b}_{\mathrm{LRG}}=2.1$, and we allow ${b}_{\mathrm{QSO}}$ to float in the range 3.2–3.8.

The redshift source distribution of CIB sources can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{dS}}_{\nu }}{{dz}}=\displaystyle \frac{c}{H(z)(1+z)}\int \,{dL}\displaystyle \frac{{dn}}{{dL}}(M,z)\displaystyle \frac{{L}_{\nu (1+z)}}{4\pi },\end{eqnarray} \tag{ 19 }$$

where dn/dL is the infrared galaxy luminosity function and ${L}_{\nu (1+z)}$ is a model for the CIB luminosity emitted at rest-frame frequency $\nu (1+z)$. Here ${b}_{\mathrm{CIB}}(k,z)$ and dn/dL are predicted using the same halo model as in Shang et al. (2012). The model for the luminosity depends on several parameters; however, we only allow three parameters in our fit to vary: a luminosity amplitude L₀, a redshift evolution parameter δ, and the dust temperature T_d.

Following Pullen et al. (2018), we extend the CIB spectrum model in Shang et al. (2012) to include [C ii] emission in the host galaxies with amplitude ${A}_{{\rm{C}}{\rm{II}}}$. We also add correlated thermal SZ emission, which we model with an amplitude A_tSZ multiplying a spectral template. See Pullen et al. (2018) for more details.

The bandpower covariance that appears in the likelihood is described in Pullen et al. (2018) and includes ℓ–${\ell }^{\prime}$ correlations from masked sky regions.

5.5. CIB and C ii Parameter Estimates

We perform our Markov chain Monte Carlo (MCMC) analysis using CosmoMC (Lewis & Bridle 2002) for the six parameters

$$\begin{eqnarray}&&{\rm{\Xi }}\equiv \{{T}_{d},\delta ,{L}_{0},{A}_{{\rm{C}}{\rm{II}}},{A}_{\mathrm{tSZ}},{b}_{\mathrm{QSO}}\}.\end{eqnarray} \tag{ 20 }$$

As in Pullen et al. (2018), we fit the parameters using the six cross-power spectra between the three high-frequency Planck bands and the quasars and LRGs presented in Figure 10; these are measured in nine bins over $100\lt {\ell }\lt 1000$. We label this fit the "Pullen + 2018" fit. We include the full covariance matrix between all cross-power spectra. Also as before, we fit the mean levels of the CIB emission in the three Planck bands (Béthermin et al. 2012), as well as 10 SFR density measurements (Madau & Dickinson 2014).

The Milky Way templates show some evidence of correlation with the quasar overdensity. This could bias the determination of the correlated [C ii] amplitude, so we additionally marginalize over a nuisance clustering anisotropy based on the CIB × quasar clustering model with amplitude α.

Figure 12 shows the parameter fits before and after including Milky Way emission templates. The posteriors for the CIB parameters and thermal SZ amplitude do not change significantly from those of Pullen et al. (2018). The posterior distribution of ${A}_{{\rm{C}}{\rm{II}}}$ is only modestly impacted by marginalizing over the Lenz et al. (2017) template. We also find that the A_Cii measurement is not degenerate with α for either Milky Way template. The MCMC constrains $-0.092\lt {\alpha }_{\mathrm{Lenz}}\lt -0.025$ and $-0.121\,\lt {\alpha }_{\mathrm{Green}}\lt -0.042$ at 95% confidence. This indication of a correlation between the Milky Way–only tracers and the quasar overdensity warrants future investigation. Synchrotron emission from the quasars could correlate with the anisotropy of the synchrotron background in the 21 cm maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.6. Interpreting the MCMC Results

We can interpret the lack of impact of Milky Way templates in Planck × BOSS quasar measurements using the multiband formalism in Section 4. To do this, recast the relevant parts of the full nonlinear model as linear modes that fit the correlated continuum emission as in Equation (14). Let the data vector be the cross-powers of the Planck 353, 545, and 857 GHz maps and BOSS quasars. Take a simple model for the SED of the correlated continuum with ${I}_{\nu }^{\mathrm{dust}}={\nu }^{\beta }B[{T}_{\mathrm{dust}},(1+z)\nu ]$, where $B({T}_{\mathrm{dust}},\nu )$ is the Planck law, and we fix β = 1.5. To build a linear model, let ${\boldsymbol{M}}$ have a mode that is the best-fit dust with ${T}_{\mathrm{dust}}=27.2$ K, evaluated at z = 2.25 (near the peak of the quasar number density). To model changes in temperature, also add a mode ${{dI}}_{\nu }^{\mathrm{dust}}/{dT}(\nu )$.

For the covariance in Equation (14), use the bandpower covariance of the Planck data. Figure 13 shows this principal mode of the bandpower covariance, which is consistent with emission from the Milky Way (β = 1.5, ${T}_{\mathrm{dust}}=19$ K from z = 0). The principal mode is spectrally similar to emission with ${T}_{\mathrm{dust}}=27.2$ K from z = 2.25 in the correlated continuum model. The similarity in spectral indices means that the operation that marginalizes over correlated continuum emission is also very effective at suppressing the Milky Way emission. In the setting of Equation (16), the correlated and uncorrelated continuum terms have similar SEDs, so $f(\nu )\approx g(\nu )$ and the second term of the variance are highly suppressed. The value of a Milky Way template map is diminished in this case.

Zoom In Zoom Out Reset image size

In the linear model, the Lenz et al. (2017) template increases the signal-to-noise ratio in the line amplitude estimate by only 8%. Adding a derivative mode ${{dI}}_{\nu }^{\mathrm{dust}}/{dT}$ (as a linearized temperature parameter) reduces the improvement of the Lenz et al. (2017) template to 3%. Most of this improvement comes from ${\ell }\lt 500$, where the Milky Way template is best correlated to Planck 545 GHz (Figure 14). In contrast, if the extragalactic dust temperature were the same as the Milky Way's $\approx 19$ K (but redshifted to z = 2.25), the spectral match between the correlated SED modes and the galactic dust would not be nearly as good. For a 19 K emitter at z = 2.25, the dust spectrum crests in the 545 GHz band, while at z = 0, the dust emission in 857 GHz is brighter. In this scenario, the Lenz et al. (2017) template would lead to a 27% improvement in the line constraint. Hence, the Planck × BOSS quasar measurements benefit from a lucky coincidence between the galactic and redshifted extragalactic emission SED. This similarity allows the marginalization over the extragalactic SED to also effectively deweight the Milky Way emission.

Section 4 argues that an intensity survey needs at least as many frequencies as dof in the correlated continuum mode and argued that suppression of continuum contamination is equivalent to throwing out ${k}_{\parallel }$ (line-of-sight wavevector) information. This section describes how correlated continuum emission appears in the full 3D power spectrum ${P}_{\times }({k}_{\perp },{k}_{\parallel })$ of tomographic surveys with numerous spectral bins.

Here we will mock up a future tomographic survey to show the correlated continuum and line emission in the full 3D power spectrum ${P}_{\times }({k}_{\perp },{k}_{\parallel })$ for wavenumbers perpendicular, ${k}_{\perp }$, and parallel, ${k}_{\parallel }$, to the line of sight. Correlated continuum emission enters at low ${k}_{\parallel }$. Cutting ${k}_{\parallel }\lt {k}_{\parallel }^{\mathrm{cut}}$ recovers the line correlation independently of the continuum but results in some loss of sensitivity at low $| k|$.

6.1. Impact of Correlated Continuum Emission in the 3D Cross-power

Before doing a numerical simulation, we can get intuition for the impact of the SED on the IM cross-power by working in the Hubble approximation. Here the SED is a direct convolution of the overdensity field in the line-of-sight direction, which becomes a multiplication in ${k}_{\parallel }$ space.

The surface brightness at frequency ν is the integral of comoving specific emission intensity $j(\nu ,z)$ across all redshifts in a tomographic survey slice,

$$\begin{eqnarray}&&{I}_{\nu }=\int \displaystyle \frac{{dz}}{1+z}\displaystyle \frac{d\chi }{{dz}}j(\nu ,z),\end{eqnarray} \tag{ 21 }$$

where $\chi (z)$ is the comoving distance. For a thin tomographic redshift slice $\delta z$, the brightness of the line emission will not depend on $\delta z$. In contrast, the correlated component of the continuum emission will be the integral in Equation (21) over the redshift slice, so it will depend linearly on the redshift width $\delta z$.

Let the continuum SED at frequency ν from galaxies at redshift z be ${\rm{\Theta }}[\nu (1+z)]$ and take a graybody dust (Blain et al. 2002). Evaluate the SED on the comoving grid points and find its discrete Fourier transform $\tilde{{\rm{\Theta }}}({k}_{\parallel })$, normalized such that $\tilde{{\rm{\Theta }}}({k}_{\parallel }=0)=1$. If the overdensity in 3D k-space is ${{\boldsymbol{\delta }}}_{k}$, then the convolution with the full line plus continuum SED is

$$\begin{eqnarray}&&{{\boldsymbol{\delta }}}_{k}^{{IM}}={S}_{L}\left[1+{\rm{\Delta }}{{zf}}_{\mathrm{cont}}\tilde{{\rm{\Theta }}}({k}_{\parallel })\right]{{\boldsymbol{\delta }}}_{k},\end{eqnarray} \tag{ 22 }$$

where ${\rm{\Delta }}z$ is the total redshift range, or ${N}_{\mathrm{chan}}\delta z$ for ${N}_{\mathrm{chan}}$ channels, and ${f}_{\mathrm{cont}}$ is the fraction of continuum to line intensity from galaxies per unit redshift. Relating this to the Planck–BOSS cross-correlation model in Section 5.4, ${f}_{\mathrm{cont}}=\tfrac{2.5}{{A}_{{\rm{C}}{\rm{II}}}}$. Here S_L is the line brightness times the cosmological linear bias ${I}_{{\rm{C}}{\rm{II}}}{b}_{{\rm{C}}{\rm{II}}}$, and we assume that the dust emission shares the same bias for simplicity.

The cross-power ${{\boldsymbol{\delta }}}_{k}^{{IM}}{{\boldsymbol{\delta }}}_{k}^{g* }$ binned onto ${k}_{\parallel }$ and ${k}_{\perp }$ is then

$$\begin{eqnarray}&&{P}_{\times }({k}_{\perp },{k}_{\parallel })={S}_{L}{b}_{g}\left[1+{\rm{\Delta }}{{zf}}_{\mathrm{cont}}\tilde{{\rm{\Theta }}}({k}_{\parallel })\right]{P}_{\delta \delta }({k}_{\perp },{k}_{\parallel }),\end{eqnarray} \tag{ 23 }$$

where ${P}_{\delta \delta }({k}_{\perp },{k}_{\parallel })$ is the power spectrum of the underlying dark matter field ${\boldsymbol{\delta }}$. This is just the expected IM cross-power ${P}_{\times }({k}_{\perp },{k}_{\parallel })={S}_{L}{b}_{g}{P}_{\delta \delta }({k}_{\perp },{k}_{\parallel })$ plus a term that multiplies the ${k}_{\parallel }$ direction. In the case where galaxies have a constant SED over this spectral range, ${\rm{\Theta }}({k}_{\parallel })$ is a delta function, and the continuum only appears at ${k}_{\parallel }=0$. In practice, continuum emission will contaminate higher ${k}_{\parallel }$ because the SED has additional spectral structure, the Hubble approximation will break down, and there may be mixing between ${k}_{\parallel }$ modes because of the survey geometry.

We simulate the effect of correlated continuum emission on the cross-power spectrum by drawing Gaussian sample maps from a z = 2.25 power spectrum (Di Dio et al. 2013). The maps are $500\,\mathrm{Mpc}/h$ cubes with a pixel size of $1\,\mathrm{Mpc}/h$, which corresponds to $2\lt z\lt 2.5$ in the line-of-sight direction. The continuum SED is then simulated by adding a graybody dust SED whose magnitude is proportional to the overdensity at each voxel and redshifted according to that voxel's redshift. To illustrate the sensitivity of a future survey, we take ${A}_{{\rm{C}}{\rm{II}}}$ a factor of 20 smaller than the posterior in Section 5.5. Since our focus is only on the ratio of correlated continuum to line emission, we do not apply the multiplicative prefactor ${S}_{L}{b}_{g}$ (Equation (23)) to our simulated intensity maps. Therefore, the units of the power spectra in Figures 15 and 16 are ${\mathrm{Mpc}}^{3}{h}^{-3}$. Figure 15 shows the two-dimensional power spectrum of one realization of this simulated model. (The impact of the continuum convolution is partly stochastic, so averaging over many simulations also averages down part of the correlated continuum.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Only the two lowest nonzero k_∥ modes are noticeably biased by the continuum SED. This is shown more clearly in Figure 16, in which we plot the power spectrum averaged over k_⊥ as a function of k_∥. Figure 17 shows the decrease in the number of modes for $| k|$ by cutting the two lowest nonzero k_∥ modes from the binned one-dimensional power spectrum.

Zoom In Zoom Out Reset image size

We can also revisit the question from Section 4.1 about the number of bands that a survey needs to treat correlated continuum emission without significant signal loss. In the tomographic case, the number of ${k}_{\parallel }$ modes corresponds roughly to the number of frequency channels in the spectrometer. Removing the five lowest ${k}_{\parallel }$ bins rejects the majority of the correlated continuum. Schematically, if there are a total of 500 channels, this results in a $5/500=1 \%$ loss in signal. Conversely, if there are 10 channels, and five ${k}_{\parallel }$ modes must be removed, 50% of the signal is lost. The calculations here have made simple assumptions for the number of dof in the correlated continuum SED. In practice, both the correlated and uncorrelated continua may be more complex through their astrophysics or interaction with the instrument.

Cross-correlation between IM data sets and galaxy redshift surveys has become a gold standard for robustness against foreground contamination. This cross-correlation is sensitive to not only the line emission in target galaxies but also the complete, correlated SED. The inference of line brightness needs to marginalize over nonline contributions in the SED. The influence of correlated continuum emission has only recently been appreciated (Serra et al. 2014; Switzer 2017; Chiang et al. 2018; Pullen et al. 2018). Here we have expanded the consideration to a more general survey to understand the circumstances when correlated continuum emission can be a significant pitfall for IM.

Correlated continuum emission should be estimated jointly (self-consistently) with bright Milky Way foreground cleaning. For most IM surveys, a synchrotron or dust continuum dominates the correlated SED. This SED is similar but not identical to the bright continuum emission from the Milky Way. We show that a properly formed likelihood will deweight foregrounds in ${C}_{{\ell }}$ space analogously to map-space cleaning through a linear combination. Doing the deweighting as part of the parameter estimation from the power spectra self-consistently treats the uncorrelated and correlated continuum emission. In contrast, an analysis that cleans in map space before the likelihood analysis requires simulations to quantify residual correlated continuum emission and signal loss.

The IM observations with a limited number of bands have unique challenges that we consider in the case of Planck × BOSS (Pullen et al. 2018). When the number of channels is approximately the number of correlated continuum parameters, there is limited latitude to (1) deweight the bright uncorrelated foregrounds while also marginalizing over the correlated continuum and (2) verify that the line amplitude estimate is not biased because of an incomplete continuum model.

We reconsider the analysis of Planck × BOSS with the addition of Milky Way–only templates. These templates have the potential to provide leverage to separate and clean Milky Way emission independently of the correlated dust continuum. However, in practice, we find that existing templates (Green et al. 2015; Lenz et al. 2017) do not yet provide significant benefit, due both to the stochasticity with which they trace the Milky Way and to a similarity in the correlated and Milky Way dust continuum spectra.

Our analysis of Planck × BOSS will have shortcomings for future IM data, which move from simple detection to astrophysical characterization. Rather than taking a parametric model for the clustered emission, future work could directly solve for the $\mathrm{SED}({z}_{i},{\nu }_{j})$ and anisotropy structure ${C}_{{\ell }}$ of the correlation, similar to a "CMB-only" approach (Dunkley et al. 2013). Direct inference of the clustering anisotropy could measure strong scale-dependent bias and stochasticity of the tracers that may not be understood well in a simple parametric model. Proposed IM experiments (Kovetz et al. 2017) have numerous frequency channels and should permit a unique measurement of the clustering of line and continuum radiation as a function of redshift.

We acknowledge David Spergel for recommendations for cleaning correlated foregrounds. This work was supported under NASA 16-ADAP16-0033, and is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Recent work (Alonso & Ferreira 2015; Bull et al. 2015; Fonseca et al. 2015; Kovetz et al. 2017; Switzer 2017) argues that the multitracer approach (McDonald & Seljak 2009; Seljak 2009; Bernstein & Cai 2011) can apply fruitfully to IM and allow a determination of line amplitudes without cosmic variance. This appendix describes this approach in map and power-spectrum space (Abramo & Leonard 2013; Blake et al. 2013) without foregrounds.

Following the definitions of Section 3, a model of the observed galaxy redshift survey and intensity map, respectively, is ${{\boldsymbol{x}}}_{g}={b}_{g}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{g}$ and ${{\boldsymbol{x}}}_{I}={S}_{L}{\boldsymbol{\delta }}+{{\boldsymbol{n}}}_{I}$. We will start with the more standard derivation of the Fisher matrix in map space (e.g., Bernstein & Cai 2011). Per ℓ mode, the covariance of the galaxy redshift survey ${x}_{g,{\ell }}$ and intensity map ${x}_{I,{\ell }}$ is

$$\begin{eqnarray}{{\boldsymbol{\Sigma }}}_{{gI}}=\left(\begin{array}{ll}{b}_{g}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{g,{\ell }} & {b}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta }\\ {b}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta } & {S}_{L}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{I,{\ell }}\end{array}\right).\end{eqnarray} \tag{ 24 }$$

The Fisher matrix element for the determination of the covariance parameters ${\theta }_{i}$ and ${\theta }_{j}$ given the observables $\{{{\boldsymbol{x}}}_{g},{{\boldsymbol{x}}}_{I}\}$ is then (Tegmark et al. 1997)

$$\begin{eqnarray}&&{F}_{{\theta }_{i},{\theta }_{j}}=\displaystyle \frac{1}{2}\mathrm{Tr}\left({{\boldsymbol{\Sigma }}}_{{gI}}^{-1}\displaystyle \frac{d{{\boldsymbol{\Sigma }}}_{{gI}}}{d{\theta }_{i}}{{\boldsymbol{\Sigma }}}_{{gI}}^{-1}\displaystyle \frac{d{{\boldsymbol{\Sigma }}}_{{gI}}}{d{\theta }_{j}}\right).\end{eqnarray} \tag{ 25 }$$

In this setup, the parameters to constrain are ${\boldsymbol{\theta }}=\{{S}_{L},{C}_{{\ell }}^{\delta \delta },{N}_{A}\}$. From the inverse of the Fisher matrix, the minimum variance of the determination of S_L for each multipole (suppressing the ℓ label) marginalized over ${C}_{{\ell }}^{\delta \delta }$ and N_I yields

$$\begin{eqnarray}&&{\sigma }_{{S}_{L}}^{2}=\displaystyle \frac{1}{{v}_{{\ell }}}\left(\displaystyle \frac{{N}_{I}}{{C}_{{\ell }}^{\delta \delta }}+\displaystyle \frac{{N}_{g}}{{C}_{{\ell }}^{\delta \delta }}\displaystyle \frac{{S}_{L}^{2}}{{b}_{g}^{2}}+\displaystyle \frac{{N}_{I}{N}_{g}}{{\left({C}_{{\ell }}^{\delta \delta }\right)}^{2}{b}_{g}^{2}}+2\displaystyle \frac{{N}_{g}^{2}}{{\left({C}_{{\ell }}^{\delta \delta }\right)}^{2}}\displaystyle \frac{{S}_{L}^{2}}{{b}_{g}^{4}}\right),\end{eqnarray} \tag{ 26 }$$

where ${v}_{{\ell }}$ is the number of modes in the measurement, which is roughly $(2{\ell }+1){f}_{\mathrm{sky}}$ for a 2D map on the fraction ${f}_{\mathrm{sky}}$ of the sky. This uses no expansion in the approximation of small noise (e.g., Bernstein & Cai 2011).

Sample variance produces a term linear in S_L², but here all factors of S_L appear with orders of ${N}_{g}/{C}_{{\ell }}^{\delta \delta }$. To simplify expressions throughout to provide some analytic insight, we will assume that galaxy survey shot noise is negligible, or ${N}_{g}/{C}_{{\ell }}^{\delta \delta }\ll 1$ (taking no terms beyond order zero). In this case, ${\sigma }_{{S}_{L}}^{2}={N}_{I}/({v}_{{\ell }}{C}_{{\ell }}^{\delta \delta })$. This has a simple interpretation as the signal-to-noise ratio per mode, divided by the number of modes.

The standard discussion above aligns well with a likelihood analysis of pixel space. Pixel-space likelihoods are expensive in practical data analysis, requiring an ${N}_{\mathrm{pix}}\times {N}_{\mathrm{pix}}$ matrix inverse for ${N}_{\mathrm{pix}}$ pixels. We instead prefer to formulate the analytic Fisher matrix in power-spectrum space, where the likelihood for parameters in Section 5 is performed. Working with a likelihood of two-point information also allows more flexibility in the parameter constraint. For example, a likelihood using Equation (24) has access to information from both the cross-power ${b}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta }$ and the auto-power ${S}_{L}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{I}$ of the intensity survey. In practice, the auto-power may not be usable because of additive bias from residual foregrounds, but it can simply be excluded from the two-point likelihood.

Here the two-point functions are

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{{II}} & = & {S}_{L}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{I},\\ {C}_{{\ell }}^{{gI}} & = & {b}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta },\\ {C}_{{\ell }}^{{gg}} & = & {b}_{g}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{g}.\end{array}\end{eqnarray} \tag{ 27 }$$

Let the observable be the vector ${\boldsymbol{\mu }}=\{{C}_{{\ell }}^{{II}},{C}_{{\ell }}^{{gI}},{C}_{{\ell }}^{{gg}}\}$ and estimate the parameters $\{{S}_{L},{C}_{{\ell }}^{\delta \delta },{N}_{I}\}$. In this case, the Fisher matrix for determination of the mean parameters (Tegmark et al. 1997) is

$$\begin{eqnarray}&&{F}_{{\theta }_{i},{\theta }_{j}}=\displaystyle \frac{1}{2}\mathrm{Tr}\left[\mathrm{Cov}{\left({\boldsymbol{\mu }}\right)}^{-1}\left(\displaystyle \frac{d{\boldsymbol{\mu }}}{d{\theta }_{i}}\displaystyle \frac{d{{\boldsymbol{\mu }}}^{T}}{d{\theta }_{j}}+\displaystyle \frac{d{\boldsymbol{\mu }}}{d{\theta }_{j}}\displaystyle \frac{d{{\boldsymbol{\mu }}}^{T}}{d{\theta }_{i}}\right)\right],\end{eqnarray} \tag{ 28 }$$

where $\mathrm{Cov}({\boldsymbol{\mu }})$ is the covariance of the two-point functions in the observation vector ${\boldsymbol{\mu }}$. The matrix elements of this covariance matrix are given by Equation (6).

Marginalizing over $\{{C}_{{\ell }}^{\delta \delta },{N}_{I}\}$ and ignoring galaxy redshift survey shot noise, one recovers the estimate of S_L with variance ${\sigma }_{{S}_{L}}^{2}={N}_{I}/({v}_{{\ell }}{C}_{{\ell }}^{\delta \delta })$, as above. Sample variance is absorbed in ${C}_{{\ell }}^{\delta \delta }$ as the explanatory variable, allowing S_L to be determined independently. Indeed, when the observable is only the cross-power ${\boldsymbol{\mu }}=\{{C}_{{\ell }}^{{gI}}\}$, and the line amplitude S_L is the only parameter, the variance ${\sigma }_{{S}_{L}}^{2}=(2{S}_{L}^{2}+{N}_{I}/{C}_{{\ell }}^{\delta \delta })/{v}_{{\ell }}$ has a term $\propto {S}_{L}^{2}$.

To simplify the analytics throughout and avoid use of the intensity auto-power, we reduce the observables to the cross-powers and the galaxy redshift survey auto-power, or ${\boldsymbol{\mu }}=\{{C}_{{\ell }}^{{gI}},{C}_{{\ell }}^{{gg}}\}$ here. This provides two equations allowing a fit of two parameters $\{{S}_{L},{C}_{{\ell }}^{\delta \delta }\}$, so it assumes that N_I is fixed and known. In this case, the variance on S_L marginalized over ${C}_{{\ell }}^{\delta \delta }$ is again ${\sigma }_{{S}_{L}}^{2}={N}_{I}/({v}_{{\ell }}{C}_{{\ell }}^{\delta \delta })$. This setup is a starting point for Section 3.

Both approaches above assume that the galaxy redshift survey and line emission are perfect tracers of the same overdensity. In practice, there is some stochasticity r between the two tracers, yielding ${C}_{{\ell }}^{{gI}}={{rb}}_{g}{S}_{L}{C}_{{\ell }}^{\delta \delta }$. In either the map or power-spectrum space approach, and assuming negligible shot noise, the line amplitude constraint with stochasticity between the tracers becomes

$$\begin{eqnarray}&&{\sigma }_{{S}_{L}}^{2}=\displaystyle \frac{{N}_{I}}{{r}^{2}{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}+\displaystyle \frac{{S}_{L}^{2}}{{v}_{{\ell }}}\left(\displaystyle \frac{1}{{r}^{2}}-1\right).\end{eqnarray} \tag{ 29 }$$

Note that this does contain cosmic variance terms $\propto {S}_{L}^{2}$ when $r\ne 1$. When there is stochasticity between the populations, the component of variance that is not correlated between the tracers adds cosmic variance. Furthermore, the thermal noise term ${N}_{I}/({r}^{2}{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta })$ scales as $1/{r}^{2}$. On the large angular scales considered in Section 5, the galaxy redshift survey and intensity trace the same underlying linear perturbations, where $r\approx 1$ is a reasonable approximation.

Evasion of cosmic variance in the line amplitude estimate can be understood in both map and power-spectrum space. In map space, the galaxy redshift survey provides a spatial template of overdensity whose amplitude can be fit in the intensity survey without cosmic variance (it is a determination of a mean amplitude rather than of variance). In power-spectrum space, cosmic variance is correlated between the intensity and galaxy redshift survey. Fitting jointly for ${C}_{{\ell }}^{\delta \delta }$ and S_L makes ${C}_{{\ell }}^{\delta \delta }$ the explanatory variable of the cosmic covariance. Also, S_L only appears as an amplitude in the intensity survey, so it is determined independently of the cosmic covariance.

Extend the model above to two intensity bands, A and B, at frequencies ${\nu }_{A}$ and ${\nu }_{B}$. Let the amplitude of the correlated signal in maps A and B be ϕ and ψ. Add foreground covariance with spatial anisotropy ${C}_{{\ell }}^{{uu}}$ and amplitude $f(\nu )$ to both intensity bands. The two-point functions are

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{{AA}} & = & {\phi }^{2}{C}_{{\ell }}^{\delta \delta }+f{\left({\nu }_{A}\right)}^{2}{C}_{{\ell }}^{{uu}}+{N}_{A}\\ {C}_{{\ell }}^{{BB}} & = & {\psi }^{2}{C}_{{\ell }}^{\delta \delta }+f{\left({\nu }_{B}\right)}^{2}{C}_{{\ell }}^{{uu}}+{N}_{B}\\ {C}_{{\ell }}^{{AB}} & = & \phi \psi {C}_{{\ell }}^{\delta \delta }+f({\nu }_{A})f({\nu }_{B}){C}_{{\ell }}^{{uu}}\\ {C}_{{\ell }}^{{Ag}} & = & \phi {b}_{g}{C}_{{\ell }}^{\delta \delta }\\ {C}_{{\ell }}^{{Bg}} & = & \psi {b}_{g}{C}_{{\ell }}^{\delta \delta }\\ {C}_{{\ell }}^{{gg}} & = & {b}_{g}^{2}{C}_{{\ell }}^{\delta \delta }+{N}_{g},\end{array}\end{eqnarray} \tag{ 30 }$$

In the limit of negligible shot noise, the covariance of determinations $\hat{\phi }$ and $\hat{\psi }$ is

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Cov}(\hat{\phi },\hat{\phi }) & = & \displaystyle \frac{f{\left({\nu }_{A}\right)}^{2}{C}_{{\ell }}^{{uu}}+{N}_{A}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }},\\ \mathrm{Cov}(\hat{\phi },\hat{\psi }) & = & \displaystyle \frac{f({\nu }_{A})f({\nu }_{B}){C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }},\\ \mathrm{Cov}(\hat{\psi },\hat{\psi }) & = & \displaystyle \frac{f{\left({\nu }_{B}\right)}^{2}{C}_{{\ell }}^{{uu}}+{N}_{B}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}.\end{array}\end{eqnarray} \tag{ 31 }$$

Extending this to many bands gives covariance the form

$$\begin{eqnarray}&&{\boldsymbol{N}}=\displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}({{\boldsymbol{N}}}_{\mathrm{th}}+{C}_{{\ell }}^{{uu}}{{\boldsymbol{ff}}}^{T}),\end{eqnarray} \tag{ 32 }$$

where ${{\boldsymbol{N}}}_{\mathrm{th}}$ is the diagonal matrix of thermal noise in each band (N_A and N_B above) and ${{\boldsymbol{f}}}_{i}=f({\nu }_{i})$ is the foreground spectrum evaluated in each of the bands. In summary, if the Milky Way adds common spatial structure to the intensity maps, it will appear as a rank-1 term in the bandpower covariance as a function of frequency.

Here we develop a simple linear model for an intensity survey with many channels. Let the multiband IM survey have frequencies labeled ${\nu }_{i}$ and find the cross-power with the galaxy redshift survey at fixed redshift z. Following the notation in Section 3, model this power spectrum as

$$\begin{eqnarray}&&{C}_{{\ell }}[{{\boldsymbol{x}}}_{I}({\nu }_{i})\times {{\boldsymbol{x}}}_{g}(z)]=[{S}_{L}\xi ({\nu }_{i})+{S}_{C}({\nu }_{i})]\cdot {b}_{g}{C}_{{\ell }}^{\delta \delta }(z).\end{eqnarray} \tag{ 33 }$$

In the case of Planck, quasars in $2\lt z\lt 3.2$ only correlate with one frequency channel, so $\xi ({\nu }_{i})=1$ in ${\nu }_{i}=545\,\mathrm{GHz}$ and $\xi ({\nu }_{i})=0$ at the other Planck frequencies. In a more typical intensity survey with narrow spectral channels, the clustering of LSS at low ${k}_{\parallel }$ produces correlations between channels in the kernel $\xi ({\nu }_{i})$ (Challinor & Lewis 2011; Di Dio et al. 2013).

In the illustrative limit of negligible shot noise, the galaxy redshift survey auto-power fixes ${C}_{{\ell }}^{\delta \delta }(z)$, including cosmic variance, making the model for the cross-powers in Equation (33) approximately linear in the amplitude $[{S}_{L}\xi ({\nu }_{i})+{S}_{C}({\nu }_{i})]$. We can therefore use a simple analytic model of the cross-power as a function of frequency to understand the joint constraint of correlated line and continuum emission. In practice with BOSS quasars, shot noise is not negligible, but the linear approach gives intuition for the geometry of the joint analysis of correlated and uncorrelated continua. Future IM experiments would be complemented well by galaxy redshift surveys with small shot noise on the scales of interest.

Form a new data vector ${{\boldsymbol{d}}}_{{\ell }}$ from the measured cross-powers at frequency ${\nu }_{i}$ divided by ${b}_{g}{\hat{C}}_{{\ell }}^{\delta \delta }(z)$ inferred from the galaxy redshift survey

$$\begin{eqnarray}&&{d}_{{\ell }}^{\times }({\nu }_{i})={\hat{C}}_{{\ell }}[{{\boldsymbol{x}}}_{I}({\nu }_{i})\times {{\boldsymbol{x}}}_{g}(z)]/({b}_{g}{\hat{C}}_{{\ell }}^{\delta \delta }(z)).\end{eqnarray} \tag{ 34 }$$

This is an estimator for ${S}_{L}\xi ({\nu }_{i})+{S}_{C}({\nu }_{i})$. Model the SED of the correlated continuum as the sum of linear components ${g}_{j}({\nu }_{i})$ and amplitudes ${S}_{C,j}$,

$$\begin{eqnarray}&&{S}_{C}({\nu }_{i})=\displaystyle \sum _{j}^{{N}_{\mathrm{comp}}}{S}_{C,j}{g}_{j}({\nu }_{i}).\end{eqnarray} \tag{ 35 }$$

Pack each spectral component into a vector ${{\boldsymbol{g}}}_{j}{| }_{i}={g}_{j}({\nu }_{i})$ and the signal correlation kernel into a vector ${\boldsymbol{\xi }}| i=\xi ({\nu }_{i})$. In the case of widely separated bins in the Planck analysis, this is a delta function at 545 GHz. Pack the line signal and ${N}_{\mathrm{comp}}$ correlated continuum modes into ${\boldsymbol{M}}\equiv [{\boldsymbol{\xi }},{{\boldsymbol{g}}}_{1},\; ...{{\boldsymbol{g}}}_{{N}_{\mathrm{comp}}}]$. The linear amplitude parameter vector is then ${\boldsymbol{\theta }}\equiv [{S}_{L},{{\boldsymbol{S}}}_{C}]$, where ${{\boldsymbol{S}}}_{C}$ are the ${S}_{C,j}$ linear parameters for the SED. With these substitutions, the model for multiband cross-powers in Equation (33) is

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{{\ell }}^{\times }={\boldsymbol{M}}{\boldsymbol{\theta }}+{{\boldsymbol{n}}}_{\times },\end{eqnarray} \tag{ 36 }$$

where ${{\boldsymbol{n}}}_{\times }$ is the noise of the cross-power measurements drawn from covariance ${{\boldsymbol{N}}}_{\times }$ that includes both thermal noise (${{\boldsymbol{N}}}_{\mathrm{th}}$) and foreground (${{\boldsymbol{N}}}_{\mathrm{fg}}$) covariance as

$$\begin{eqnarray}&&{{\boldsymbol{N}}}_{\times }=\displaystyle \frac{1}{{\nu }_{{\ell }}{C}_{{\ell }}^{\delta \delta }}({{\boldsymbol{N}}}_{\mathrm{th}}+{{\boldsymbol{N}}}_{\mathrm{fg}}).\end{eqnarray} \tag{ 37 }$$

(Note that the data vector defined in Equation (34) is an estimator for the amplitudes for ϕ and ψ in Appendix B but applied to many bands. This results in the same form for the covariance as Equation (32)).

This model solves jointly for the line amplitude S_L and spectral baseline amplitudes ${{\boldsymbol{S}}}_{C}$ that describe the correlated continuum emission (Switzer 2017). The linear estimate $\hat{{\boldsymbol{\theta }}}$ and its covariance ${{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}$ are

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{\theta }}} & = & {{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}{{\boldsymbol{M}}}^{T}{{\boldsymbol{N}}}_{\times }^{-1}{{\boldsymbol{d}}}_{{\ell }},\\ {{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}} & = & {\left({{\boldsymbol{M}}}^{T}{{\boldsymbol{N}}}_{\times }^{-1}{\boldsymbol{M}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 38 }$$

We have written this expression for a single multipole ℓ for simplicity, but the model could be extended as a likelihood on the cross-powers at all ℓ in practice.

This linear form recovers our primary results so far. Equation (38) recovers Equation (9) in the case where ${\boldsymbol{M}}={[[1,0]}^{T}]$ (only S_L is fit from the data) and r_F = 1. Here ${{\boldsymbol{N}}}_{\mathrm{fg}}={C}_{{\ell }}^{{uu}}{[f({\nu }_{I}),f({\nu }_{V})]}^{T}[f({\nu }_{I}),f({\nu }_{V})]$ and ${{\boldsymbol{N}}}_{\mathrm{th}}\,=\mathrm{diag}[{N}_{I},{N}_{V}]$. Adding one correlated SED component as ${\boldsymbol{M}}={[{[1,0]}^{T},[g({\nu }_{I}),g({\nu }_{V})]}^{T}]$ recovers Equation (12), again with r_F = 1.

We can now consider the effect of Milky Way foregrounds by letting the covariance ${{\boldsymbol{N}}}_{\times }^{-1}$ be the rank-1 foreground model Equation (32), equivalent to a spatial pattern common across all frequencies. In practice, Milky Way foregrounds will have higher rank, which can be included in ${{\boldsymbol{N}}}_{\times }^{-1}$. In this simplified setting, we can continue analytically by applying the Sherman–Morrison formula twice to give

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}} & = & \displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left({{\boldsymbol{\Sigma }}}_{\mathrm{th}}+\displaystyle \frac{{C}_{{\ell }}^{{uu}}{{\boldsymbol{\alpha }}}_{f}{{\boldsymbol{\alpha }}}_{f}^{T}}{1+{C}_{{\ell }}^{{uu}}{\rm{\Delta }}}\right)\end{array}\,,\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }} & \equiv & {{\boldsymbol{f}}}^{T}{{\boldsymbol{N}}}_{\mathrm{th}}^{-1}({\boldsymbol{f}}-{\boldsymbol{M}}{{\boldsymbol{\alpha }}}_{f}),\\ {{\boldsymbol{\Sigma }}}_{\mathrm{th}} & \equiv & {\left({{\boldsymbol{M}}}^{T}{{\boldsymbol{N}}}_{\mathrm{th}}^{-1}{\boldsymbol{M}}\right)}^{-1},\\ {{\boldsymbol{\alpha }}}_{f} & \equiv & {{\boldsymbol{\Sigma }}}_{\mathrm{th}}{{\boldsymbol{M}}}^{T}{{\boldsymbol{N}}}_{\mathrm{th}}^{-1}{\boldsymbol{f}}.\end{array}\end{eqnarray} \tag{ 40 }$$

Here ${{\boldsymbol{\Sigma }}}_{\mathrm{th}}$ is the variance for the estimate of ${\boldsymbol{\theta }}$ with only thermal noise and no foregrounds, and ${{\boldsymbol{\alpha }}}_{f}$ is the optimal linear estimator, with thermal noise only, applied to the foreground SED. Hence, ${C}_{{\ell }}^{{uu}}{{\boldsymbol{\alpha }}}_{f}{{\boldsymbol{\alpha }}}_{f}^{T}$ is the amplitude of the foregrounds in the ${\boldsymbol{\theta }}$ parameter basis. Finally, Δ is a discriminant that describes how much information about the uncorrelated foregrounds is left after fitting for the correlated spectral components in ${\boldsymbol{M}}$.

We can now consider the geometry of the joint correlated and uncorrelated continuum separation problem. The key quantity here is the discriminant Δ. This is a measure of the remaining rms of the uncorrelated foregrounds ${\boldsymbol{f}}$ after fitting out the correlated continuum modes ${{\boldsymbol{g}}}_{i}$. (To see this, note that ${{\boldsymbol{\alpha }}}_{f}$ projects ${\boldsymbol{f}}$ onto the parameter basis ${\boldsymbol{\theta }}$, and then ${\boldsymbol{M}}$ projects onto a spectrum differenced with ${\boldsymbol{f}}$ in Δ.)

If ${\boldsymbol{M}}$ is a complete spectral basis, then ${\rm{\Delta }}=0$ regardless of the spectral shape of ${\boldsymbol{f}}$. This will be the case any time the number of independent parameters in the correlated SED model is equal to the number of bands, such as in the case of Planck × BOSS. In this regime,

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}=\displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left({{\boldsymbol{\Sigma }}}_{\mathrm{th}}+{C}_{{\ell }}^{{uu}}{{\boldsymbol{\alpha }}}_{f}{{\boldsymbol{\alpha }}}_{f}^{T}\right).\end{eqnarray} \tag{ 41 }$$

These two terms correspond to the thermal noise and foreground components in Equation (12) with r_F = 1. The second term represents foregrounds that contaminate the line amplitude estimate after marginalizing the correlated continuum spectral modes ${{\boldsymbol{g}}}_{j}$. An observation with few bands may not provide sufficient dof to deweight uncorrelated continuum foregrounds and model the correlated continuum foreground.

Adding a template leads to a qualitative improvement in the foreground cleaning. When ${N}_{\mathrm{band}}={N}_{\mathrm{param}}$, Equation (41) shows that the line amplitude variance always has a term proportional to ${C}_{{\ell }}^{{uu}}$. By adding a template, ${C}_{{\ell }}^{{uu}}{\rm{\Delta }}$ in the denominator of Equation (39) can remain large. In the limit of large foregrounds,

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}=\displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left({{\boldsymbol{\Sigma }}}_{\mathrm{th}}+\displaystyle \frac{{{\boldsymbol{\alpha }}}_{f}{{\boldsymbol{\alpha }}}_{f}^{T}}{{\rm{\Delta }}}\right).\end{eqnarray} \tag{ 42 }$$

Now the foreground brightness ${C}_{{\ell }}^{{uu}}$ drops out, and the numerator ${{\boldsymbol{\alpha }}}_{f}{{\boldsymbol{\alpha }}}_{f}^{T}$ is suppressed by Δ. Again in this limit, the variance of the S_L estimate is not compromised by bright foregrounds, making this a promising approach. This is the extension of Equation (16) to multiple bands.

Extend the template model of Equation (15) to allow for imperfect correlation ${r}_{M}={C}_{{\ell }}^{{Mu}}/\sqrt{{C}_{{\ell }}^{{MM}}{C}_{{\ell }}^{{uu}}}$ with the actual foreground ${{\boldsymbol{x}}}_{u}$ in the intensity map. The correlations between this map and itself (MM), the intensity map (MI), and veto map (MV) are

$$\begin{eqnarray}&&{C}_{{\ell }}^{{MM}}={A}_{M}^{2}{C}_{{\ell }}^{{uu}}+{N}_{M},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{C}_{{\ell }}^{{MI}}={r}_{M}{A}_{M}f({\nu }_{I}){C}_{{\ell }}^{{uu}},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{C}_{{\ell }}^{{MV}}={r}_{M}{A}_{M}f({\nu }_{V}){C}_{{\ell }}^{{uu}}.\end{eqnarray} \tag{ 45 }$$

For simplicity, let ${N}_{M}\to 0$, so the template has no noise but is still an imperfect tracer of the true foreground. Then,

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{S}_{L}}^{2} & = & \displaystyle \frac{{C}_{{\ell }}^{{uu}}}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}{\left[f({\nu }_{I})-f({\nu }_{V})\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}(1-{r}_{M}^{2})\\ & & +\displaystyle \frac{1}{{v}_{{\ell }}{C}_{{\ell }}^{\delta \delta }}\left[{N}_{I}+{\left[\displaystyle \frac{g({\nu }_{I})}{g({\nu }_{V})}\right]}^{2}{N}_{V}\right].\end{array}\end{eqnarray} \tag{ 46 }$$

In this paper, we argue that redshift-space distortions (RSD) can be neglected for the bulk of our analysis. To prove this, we first derive the angular power spectrum, which is directly related to the CIB-LSS angular cross-power spectrum we use in our analysis, under the Limber approximation (Loverde & Afshordi 2008). We begin with the full expression for the angular power spectrum in redshift space given by ${C}_{{\ell }}={C}_{{\ell }}^{00}+2\beta {C}_{{\ell }}^{0r}+{\beta }^{2}{C}_{{\ell }}^{{rr}}$, where β is the RSD parameter given by the growth rate-clustering bias ratio,

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{00} & = & \displaystyle \frac{2}{\pi }\int {dk}\,{k}^{2}P(k){[{W}_{{\ell }}^{0}(k)]}^{2},\\ {C}_{{\ell }}^{0r} & = & \displaystyle \frac{2}{\pi }\int {dk}\,{k}^{2}P(k){W}_{{\ell }}^{0}(k){W}_{{\ell }}^{r}(k),\\ {C}_{{\ell }}^{{rr}} & = & \displaystyle \frac{2}{\pi }\int {dk}\,{k}^{2}P(k){[{W}_{{\ell }}^{r}(k)]}^{2},\end{array}\end{eqnarray} \tag{ 47 }$$

where P(k) is the matter power spectrum at redshift z = 0 and the window function ${W}_{{\ell }}(k)={W}_{{\ell }}^{0}(k)+{W}_{{\ell }}^{r}(k)$. Given the selection function $\phi (r)$, we can write (Padmanabhan et al. 2007)

$$\begin{eqnarray}&&\begin{array}{l}{W}_{{\ell }}^{0}(k)=\int {dr}\,\phi (r){j}_{{\ell }}({kr})\\ {W}_{{\ell }}^{r}(k)=\int {dr}\,\phi (r)\left[\displaystyle \frac{2{{\ell }}^{2}+2{\ell }-1}{(2{\ell }-1)(2{\ell }+3)}{j}_{{\ell }}({kr})\right.\\ \left.-\displaystyle \frac{{\ell }({\ell }-1)}{(2{\ell }-1)(2{\ell }+1)}{j}_{{\ell }-2}({kr})-\displaystyle \frac{({\ell }+1)({\ell }+2)}{(2{\ell }+1)(2{\ell }+3)}{j}_{{\ell }+2}({kr})\right].\end{array}\end{eqnarray} \tag{ 48 }$$

We wish to perform operations similar to those in Loverde & Afshordi (2008) to get the Limber approximation for this expression. Using Equations (7)–(11) in Loverde & Afshordi (2008), we can show that the Limber approximation at first order is

$$\begin{eqnarray}\begin{array}{rcl}\sqrt{\displaystyle \frac{2}{\pi }}\int {dr}\,\phi (r){j}_{{\ell }}({kr}) & \simeq & \displaystyle \frac{1}{k\sqrt{{\ell }+1/2}}\phi \left(\displaystyle \frac{{\ell }+1/2}{k}\right)\\ & \simeq & \displaystyle \frac{1}{k\sqrt{{\ell }}}\phi \left(\displaystyle \frac{{\ell }}{k}\right),\end{array}\end{eqnarray} \tag{ 49 }$$

where we approximate ${\ell }+1/2\simeq {\ell }$ to make the following equations less cluttered. Inserting this into Equation (47), we can write

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{00} & = & \displaystyle \frac{1}{{\ell }}\int {dk}\,{\phi }^{2}\left(\displaystyle \frac{{\ell }}{k}\right)P(k),\\ {C}_{{\ell }}^{0r} & = & \displaystyle \frac{1}{\sqrt{{\ell }}}\int {dk}\,\phi \left(\displaystyle \frac{{\ell }}{k}\right)F({\ell },k)P(k),\\ {C}_{{\ell }}^{{rr}} & = & \int {dk}\,{F}^{2}({\ell },k)P(k),\end{array}\end{eqnarray} \tag{ 50 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}F({\ell },k) & = & \displaystyle \frac{2{{\ell }}^{2}+2{\ell }-1}{(2{\ell }-1)(2{\ell }+3)}\displaystyle \frac{\phi ({\ell }/k)}{\sqrt{{\ell }}}\\ & & -\ \displaystyle \frac{{\ell }({\ell }-1)}{(2{\ell }-1)(2{\ell }+1)}\displaystyle \frac{\phi [({\ell }-2)/k]}{\sqrt{{\ell }-2}}\\ & & -\ \displaystyle \frac{({\ell }+1)({\ell }+2)}{(2{\ell }+1)(2{\ell }+3)}\displaystyle \frac{\phi [({\ell }+2)/k]}{\sqrt{{\ell }+2}}.\end{array}\end{eqnarray} \tag{ 51 }$$

Now we wish to use this expression to show when RSD can be neglected. In this argument, we will assume $\phi (r)$ is a top-hat distribution centered at r = R with full width ${\rm{\Delta }}R$. In this case, it is evident that if ${\ell }/k-R\ll {\rm{\Delta }}R$, then $\phi [({\ell }\pm 2)/k]\,\simeq \phi ({\ell }/k)$. Additionally, assuming ℓ is large enough such that $\sqrt{{\ell }\pm 2}\simeq \sqrt{{\ell }}$, we can show that in this case, $F({\ell },k)\,\simeq 1/(4{{\ell }}^{5/2})$, which decreases rapidly to zero, eliminating the RSD effect. Taking the approximation that ${\ell }/k\sim R$, we argue that all the ϕs are nonzero when $2/k\lt {\rm{\Delta }}R/2$, which implies that the RSD effect vanishes when ${\ell }\gtrsim 4R/{\rm{\Delta }}R$. For the quasars, $4R/{\rm{\Delta }}R\sim 17$, while for the LRGs, $4R/{\rm{\Delta }}R\sim 10$. Thus, we are justified in neglecting RSD for our analysis.