Intensity Mapping in the Presence of Foregrounds and Correlated Continuum Emission

TL;DR

This paper discusses methods for separating line emission from foreground and continuum contamination in intensity mapping surveys, with a focus on cross-correlation techniques and future multi-band survey capabilities.

Contribution

It introduces a joint estimation approach for correlated continuum emission and foreground cleaning, applied to Planck and BOSS data for high-redshift CII detection.

Findings

Effective separation of line and continuum emission in cross-correlation.

Challenges of foreground contamination in low-band surveys.

Potential of multi-band surveys to improve emission separation.

Abstract

Intensity mapping has attracted significant interest as an approach to measure the properties of the interstellar medium in typical galaxies at high redshift. Intensity mapping measures the statistics of surface brightness as a function of frequency, making it sensitive not only to all line emission of interest but also radiation from all other sources. Significant effort has gone into developing approaches that reject foreground contamination. Additionally, the target galaxies have multiple sources of emission that can complicate the interpretation of the line brightness. We describe the problem of jointly estimating correlated continuum emission and cleaning uncorrelated continuum emission, such as from the Milky Way. We apply these considerations to a cross-correlation of Planck data with BOSS quasars for a determination of CII for 2 < z < 3.2. Intensity mapping surveys with few bands have unique challenges for treating foregrounds and avoiding bias from correlated continuum emission. We show how a future intensity mapping survey with many bands can separate line from continuum emission in cross-correlation.

Figures

Equations

$${\mathbf{x}}_{\mathrm{clean}}={\mathbf{x}}_{545}-{\kappa }_{353}{\mathbf{x}}_{353}-{\kappa }_{857}{\mathbf{x}}_{857}.$$

$$\mathbf{\kappa }=({\mathbf{A}}^T{\mathbf{N}}^{-1}\mathbf{A}{)}^{-1}{\mathbf{A}}^T{\mathbf{N}}^{-1}{\mathbf{x}}_{545}$$

$${\mathbf{x}}_g$$

$${\mathbf{x}}_I$$

$$C_{\ell }^{gI}$$

$$C_{\ell }^{gg}$$

$$2\mathcal{L}=\ln \det \mathbf{\Sigma }(\mathbf{\theta })+[\mathbf{d}-\mathbf{\mu }(\mathbf{\theta }){]}^T{\mathbf{\Sigma }}^{-1}(\mathbf{\theta })[\mathbf{d}-\mathbf{\mu }(\mathbf{\theta })].$$

$$\mathrm{Cov}(C_{\ell }^{AB},C_{{\ell }^{\prime }}^{CD})=\frac{1}{v_{\ell }}{\delta }_{\ell ,{\ell }^{\prime }}(C_{\ell }^{AD}{C}_{\ell }^{BC}+C_{\ell }^{AC}{C}_{\ell }^{BD}),$$

$${\sigma }_{S_L}^2=\frac{N_I+f({\nu }_I{)}^2{C}_{\ell }^{uu}}{C_{\ell }^{\delta \delta }{v}_{\ell }},$$

$${\mathbf{x}}_g$$

$${\mathbf{x}}_I$$

$${\mathbf{x}}_V$$

$${\sigma }_{S_L}^2$$

$$C_{\ell }^{gI}{|}_{\mathrm{clean}}=C_{\ell }^{gI}-[f({\nu }_I)/f({\nu }_V)]C_{\ell }^{gV}.$$

$${\mathbf{x}}_g$$

$${\mathbf{x}}_I$$

$${\mathbf{x}}_V$$

$${\sigma }_{S_L}^2$$

$${\hat{S}}_L=S_L+S_C\left[g^{\prime }({\nu }_I)-\frac{g({\nu }_I)}{g({\nu }_V)}{g}^{\prime }({\nu }_V)\right],$$

$${\mathbf{\Sigma }}_{\mathbf{\theta }}=({\mathbf{M}}^T{\mathbf{N}}_{\times }^{-1}\mathbf{M}{)}^{-1}.$$

$${\mathbf{x}}_g$$

$${\mathbf{x}}_I$$

$${\mathbf{x}}_V$$

$${\mathbf{x}}_M$$

$${\sigma }_{S_L}^2$$

$${\sigma }_{S_L}^2{|}_{\mathrm{stoc}}=\frac{C_{\ell }^{uu}}{v_{\ell }{C}_{\ell }^{\delta \delta }}{\left[f({\nu }_I)-f({\nu }_V)\frac{g({\nu }_I)}{g({\nu }_V)}\right]}^2(1-r_M^2).$$

$$C_{\ell }^{\mathrm{LSS}-\mathrm{CIB}}$$

$$\frac{d{S}_{\nu }}{dz}=\frac{c}{H(z)(1+z)}\int dL\frac{dn}{dL}(M,z)\frac{L_{\nu (1+z)}}{4\pi },$$

$$\Xi \equiv \{T_d,\delta ,L_0,A_{\mathrm{CII}},A_{\mathrm{tSZ}},b_{\mathrm{QSO}}\}.$$

$$I_{\nu }=\int \frac{dz}{1+z}\frac{d\chi }{dz}j(\nu ,z),$$

$${\mathbf{\delta }}_k^{IM}=S_L\left[1+\Delta z{f}_{\mathrm{cont}}\tilde{\Theta }(k_{\parallel })\right]{\mathbf{\delta }}_k,$$

$$P_{\times }(k_{\perp },k_{\parallel })=S_L{b}_g\left[1+\Delta z{f}_{\mathrm{cont}}\tilde{\Theta }(k_{\parallel })\right]P_{\delta \delta }(k_{\perp },k_{\parallel }),$$

$${\mathbf{\Sigma }}_{gI}=\left(\begin{array}{cc}{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)$$

$$F_{{\theta }_i,{\theta }_j}=\frac{1}{2}Tr\left({\mathbf{\Sigma }}_{gI}^{-1}\frac{d{\mathbf{\Sigma }}_{gI}}{d{\theta }_i}{\mathbf{\Sigma }}_{gI}^{-1}\frac{d{\mathbf{\Sigma }}_{gI}}{d{\theta }_j}\right).$$

$${\sigma }_{S_L}^2=\frac{1}{v_{\ell }}\left(\frac{N_I}{C_{\ell }^{\delta \delta }}+\frac{N_g}{C_{\ell }^{\delta \delta }}\frac{S_L^2}{b_g^2}+\frac{N_I{N}_g}{(C_{\ell }^{\delta \delta }{)}^2{b}_g^2}+2\frac{N_g^2}{(C_{\ell }^{\delta \delta }{)}^2}\frac{S_L^2}{b_g^4}\right),$$

Full text

1. Introduction

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 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 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; Wyithe et al., 2008; Chang et al., 2008; Visbal & Loeb, 2010; Visbal et al., 2011) overcomes many of these challenges, so is attracting increasing interest and investment, as detailed in the recent whitepaper 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 IGM. IM 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.

Intensity mapping has already provided several measurements that bear on galaxy evolution. Cross-correlation of $21$ cm emission and WiggleZ (Switzer et al., 2013) has determined the HI abundance at $z{\approx}0.8$, which is consistent with feedback from active galactic nuclei (Villaescusa-Navarro et al., 2016; Padmanabhan et al., 2017). Preliminary indication of [C ii] emission in cross-correlation between Planck and the BOSS quasar sample for $2<z<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<z<3.3$, ruling out several models. Croft et al. (2018) find diffuse ${\rm Ly}\alpha$ emission in cross-correlation with the ${\rm Ly}\alpha$ forest and quasars at a level that is consistent with hydrodynamic simulations.Kovetz et al. (2017) describes several dedicated IM experiments, both ongoing and future.

An interpretation of an autonomous intensity mapping 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 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<z<3.2$. [C ii] is a tracer of great interest for intensity mapping. Star formation excites the $157.7\,\mu{\rm m}$ ($1900$ GHz) ${}^{2}P_{3/2}\rightarrow\,^{2}P_{1/2}$ fine structure transition of [C ii] in several phases of the ISM. [C ii] is the brightest far-IR cooling line of star-forming galaxies, emitting as much as $0.65\%$ of the typical $L_{\rm FIR}{\sim}10^{12}{\rm L_{\odot}}$ far-IR luminosity of the host (Stacey et al., 2010). There is a well-established log-linear relation between [C ii] and the star formation rate (SFR) (Pineda et al., 2014; De Looze et al., 2014; Herrera-Camus et al., 2015) which 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^{+5.0}_{-4.8}\times 10^{4}$ Jy/sr, 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<z<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.

Intensity mapping is typically described in terms of one isolated line, such as [C ii] above. However, the target galaxies also contribute a characteristic spectral energy distribution, which includes bright dust or synchrotron continuum contribution at most frequencies of interest for intensity mapping. The overdensity field traced by a companion galaxy redshift survey will, therefore, have non-zero 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 spectral energy distribution (SED) of galaxies at $z_{j}$ emitting into $\nu_{i}$, ${\rm SED}(\nu_{i},z_{j})$ (Serra et al., 2014; Pullen et al., 2018; Switzer, 2017; Chiang 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 ${\rm 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 intensity mapping 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 $\times$ 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 intensity mapping 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 intensity mapping.

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 intensity mapping with both correlated and uncorrelated continuum emission in a simplified setting. It argues for the utility to Milky Way-only templates, which we apply to the Planck cross BOSS analysis in Section 5. Marginalizing over a correlated continuum model amounts to throwing out those degrees of freedom 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}$.

2. Milky Way foregrounds in the Planck data

Pullen et al. (2018) present a preliminary indication of [C ii] emission from $2<z<3.2$ that is correlated with large-scale structure 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 GHz and 857 GHz help to isolate dust continuum and tSZ emission from the [C ii] of interest. Fig. 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<z<3.2$, so contributes variance but not bias to the cross-correlation with quasars. Fig. 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.

Pullen et al. (2018) perform 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

[TABLE]

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$ GHz and $857$ GHz maps into a matrix ${\bf A}$. Then the linear combination coefficients which minimize variance are given by

[TABLE]

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

Fig. 4 shows the cross-power of quasars with the cleaned map compared to the input 545 GHz map.111An internal linear combination can be biased by spurious correlations between signal and foreground, especially on large angular scales. Simulations show $<1\%$ bias for $\ell>200$ (range reported here), but as we will argue, this map-domain cleaning is primarily illustrative and is not used in parameter estimates. 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.

The linear combination that cleans the Milky Way from the 545 GHz map also removes much of the correlated continuum emission. If the dust continua 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 independently 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 Fig. 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, 2) a marginalization over the correlated dust continuum to separate the [C ii] contribution. The clustering of [C ii] emission with quasars in $2<z<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 correlated dust continuum depends on information from 353 and 857 GHz, so becomes entangled with attempts to use those bands to clean Milky Way emission at 545 GHz.

3. Sensitivity with continuum emission

This section builds up a set of analytic results for multi-tracer analysis of intensity mapping data with correlated and uncorrelated continuum emission contamination. It has become standard to clean continuum contamination in intensity mapping data in the map domain, and to simulate the impact of 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 multi-tracer 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, 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), EXCLAIM222Eric Switzer, private communication, 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 multi-tracer 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 multi-tracer (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 ${\bf x}_{g}$ as a tracer with linear bias $b_{g}$ and shot noise ${\bf n}_{g}$. Model the intensity map as a line surface brightness amplitude $S_{L}$ times overdensity ${\boldsymbol{\delta}}$ plus noise ${\bf 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 ${\bf x}_{u}$ and spectral dependence $f(\nu)$. These appear in the intensity map at frequency $\nu_{I}$ as $f(\nu_{I}){\bf x}_{u}$. All vectors refer to 2D map slices at constant frequency or redshift. The observed galaxy redshift survey (${\bf x}_{g}$) and intensity (${\bf x}_{I}$) maps are

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $C^{AB}_{\ell}$ is the cross-power of fields $A$ and $B$. $v_{\ell}$ is the number of modes in the measurement, which is roughly $(2\ell+1)f_{\rm sky}\Delta\ell$ for a 2D map on the fraction $f_{\rm sky}$ of the sky and a bandpower of width $\Delta\ell$. The delta function $\delta_{\ell,\ell^{\prime}}$ is exact for all-sky, and for simplicity we assume 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 $\ell$ is (marginalized over $C^{\delta\delta}_{\ell}$)

[TABLE]

where $C^{uu}_{\ell}$ is the power spectrum of the uncorrelated emission (the map ${\bf x}_{u}$), $N_{I}$ is the noise variance of the intensity map, and $C^{\delta\delta}_{\ell}$ 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 Eq. 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^{2}_{S_{L}}$ as a term $\propto S_{L}^{2}$. As described in Appendix A, this is avoided by adding covariance to the galaxy auto-power $C^{gg}_{\ell}$ 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 ${\bf x}_{V}$ at a different frequency $\nu_{V}$. The uncorrelated continuum foreground $f(\nu_{I}){\bf x}_{u}$ in the intensity map will appear at this frequency as $f(\nu_{V}){\bf x}_{u^{\prime}}$, where $f(\nu_{V})$ accounts for the SED and $u^{\prime}$ in ${\bf x}_{u^{\prime}}$ denotes that the spatial pattern of uncorrelated foregrounds in the veto map may not be fully coherent with the spatial pattern ${\bf x}_{u}$ in the intensity map. Model this set of maps as

[TABLE]

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 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 originate from considerably different redshifts. In this case, the redshifted line emission in the veto band amounts to uncorrelated noise, so 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 (Di Dio et al., 2013; Challinor & Lewis, 2011) can calculate inter-band signal correlations, and several analyses of large-scale structure (Gaztañaga et al., 2012; Nicola et al., 2014; Di Dio et al., 2014; Salazar-Albornoz et al., 2017) have used similar approaches that could be adopted for intensity mapping. Appendices B and C describe the multi-band case, which can accommodate more general spectral energy distributions of the continuum correlations and inter-band signal correlations. Section 4.2 and Section 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 (${\bf x}_{u}$) and in the veto map (${\bf x}_{u^{\prime}}$), define “foreground" stochasticity $r_{F}=C^{uu^{\prime}}_{\ell}/\sqrt{C^{uu}_{\ell}C^{u^{\prime}u^{\prime}}_{\ell}}$, where $C^{uu^{\prime}}_{\ell}$ is the cross-power between ${\bf x}_{u}$ and ${\bf x}_{u^{\prime}}$.

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

[TABLE]

To interpret this result, make a cleaned map that removes uncorrelated continuum emission using the veto map, as ${\bf x}_{\rm clean}={\bf x}_{I}-[f(\nu_{I})/f(\nu_{V})]{\bf x}_{V}$. In the limit that the foreground in the veto band perfectly traces the intensity map ($r_{F}=1$), then ${\bf x}_{\rm clean}=S_{L}{\boldsymbol{\delta}}+{\bf n}_{I}-[f(\nu_{I})/f(\nu_{V})]{\bf 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 Eq. 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^{uu}_{\ell}$, 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 analogously to Eq. 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 ${\bf x}_{\rm clean}$ arises self-consistently with the parameter estimate through the action ${\boldsymbol{\Sigma}}^{-1}{\bf d}$ in the likelihood, Eq. 5. Foregrounds common to both the intensity and veto map show up as correlated noise between the crosspowers $C^{gI}_{\ell}$ and $C^{gV}_{\ell}$. Specifically, these off-diagonal terms in ${\boldsymbol{\Sigma}}$ appear from the cross-variance $C^{IV}_{\ell}=f(\nu_{I})f(\nu_{V})C^{uu^{\prime}}_{\ell}$. 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^{gI}_{\ell}$ and $C^{gV}_{\ell}$ to form a cleaned cross-power analogous to the cleaned map, as

[TABLE]

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 Eq. 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 Eq. 8,

[TABLE]

Continue to use the observation vector ${\bf d}=\{C^{gI}_{\ell},C^{gV}_{\ell},C^{gg}_{\ell}\}$, but expand the parameters to include an estimate of the correlated continuum amplitude $S_{C}$, as ${\boldsymbol{\theta}}=\{S_{L},S_{C},C^{\delta\delta}_{\ell}\}$. 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 Eq. 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

[TABLE]

To interpret this result, make a cleaned map that removes correlated continuum emission using the veto map, as ${\bf x}_{\rm clean}={\bf x}_{I}-[g(\nu_{I})/g(\nu_{V})]{\bf x}_{V}$. Unlike the previous section, this map cleans the correlated continuum emission using a linear combination of bands. The first term in Eq. 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 continuum 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 ${\bf x}_{I}$ and ${\bf x}_{V}$, which contributes variance to $\sigma^{2}_{S_{L}}$ even if the correlated and uncorrelated continuum 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 (Eq. 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

[TABLE]

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

In the case where the number of instrument bands $N_{\rm band}$ equals the number of parameters $N_{\rm param}$, the spectral model uses all degrees of freedom and there is no additional handle to identify bias in the line amplitude. The model may have a excellent goodness of fit because it is using $S_{L}$ to fit a residual from an insufficient continuum model. If $N_{\rm band}\gg N_{\rm 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 which 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.

4. Extension to 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_{\rm band}\times N_{\rm band}$ inverses for $N_{\rm band}$ maps for each parameter, followed by an analytic $N_{\rm param}\times N_{\rm 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 $\ell$ as a set of correlated continuum and line parameters ${\boldsymbol{\theta}}$ which linearly describe the measured cross-powers as a function of frequency ${\bf d}_{\ell}$ as ${\bf d}_{\ell}={\bf M}{\boldsymbol{\theta}}+{\bf n}_{\times}$, where the bandpowers have noise ${\bf n}_{\times}$ described by covariance ${\bf N}_{\times}$. The parameters ${\boldsymbol{\theta}}=\{S_{L},{\bf S}_{C}\}$, where $S_{L}$ is the line amplitude and ${\bf 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 ${\bf M}\equiv[{\boldsymbol{\xi}},{\bf g}_{1},...{\bf g}_{N_{\rm comp}}]$ holds the spectral templates for the line correlation ${\boldsymbol{\xi}}$ and $N_{\rm comp}$ component spectral modes of the correlated continuum. ${\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$ 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 ${\bf 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

[TABLE]

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 Eq. 14 emphasizes an accounting for overall degrees of freedom (DOF). Let the line signal have one degree of freedom (the line amplitude $S_{L}$), the correlated continuum model have $D_{\rm corr}$ DOF, and the uncorrelated continuum emission model have $D_{\rm uncorr}$ DOF. Assuming these modes are mutually independent, the total number of degrees of freedom needed to account for all terms is $D_{\rm corr}+D_{\rm uncorr}+1$. If a survey has $D_{\rm corr}+1$ bands, it can separate correlated line and continuum emission, but it has no degrees of freedom remaining to suppress uncorrelated (Milky Way) emission. In the model of Eq. 11, there are two intensity map channels, and one degree of freedom for line and correlated continuum emission respectively. This leaves no additional freedom to remove uncorrelated foregrounds, and these show up in the second term of Eq. 12.

Equation 14 permits a rapid simulation of an instrument with $N_{\rm band}$ uniform bands for a hypothetical [C ii] survey from $2<z<2.5$ ($543$ GHz and $633$ GHz) to illustrate degrees of freedom. For correlated emission, we fit an unknown amplitude for a $T_{\rm dust}=27.2$ K, $\beta=1.5$ dust source at $z=2.25$ ($I^{\rm dust}_{\nu}=\nu^{\beta}B[T_{\rm dust},(1+z)\nu]$ where $B(T_{\rm dust},\nu)$ is the Planck law) and for the amplitude of the derivative of intensity with respect to temperature $dI_{\nu}^{\rm dust}/dT$, corresponding to a fit to the dust temperature to first order. There are two correlated degrees of freedom 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 $\beta=1.5$ with an amplitude that is $10\times$ thermal noise (after marginalizing over the correlated continuum). We take the signal correlation to exist only in one spectral band, but Sec. 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 degree of freedom to suppress the Milky Way, independently of the correlated continuum determination.

The simulation in this section is set up to emphasize the impact of DOF in a general survey. Sec. 5.6 uses Eq. 14 to interpret the Planck cross 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}$ (consequently at low-$\ell$) can extend over many spectral bins in the intensity map. In the case of Planck, this effect is negligible between the bands $353$ GHz, $545$ GHz, and $857$ GHz (at the $\ell$ reported here), but the effect will be significant for future tomographic surveys with many narrow bands.

Sec. 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 Eq. 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 Di Dio et al. (2013); Challinor & Lewis (2011).

Fig. 6 demonstrates the inflation of errors at low $k_{\parallel}$, in the same simulation setup as Sec. 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 continuum.

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

The formalism and discussion here have applied to the cross-power between an intensity survey and galaxy redshift survey tracer. It is useful to contrast this to considerations for the intensity survey auto-power. In the cross-power, the uncorrelated continuum (Milky Way) adds variance, so the aim is simply to reduce the Milky Way variance. In the auto-power, residual Milky Way variance translates directly into a bias that is difficult to model because the Milky Way is 1) much brighter than the correlated continuum (and so instrumental response is even more critical), 2) is not statistically isotropic, and 3) does not have a uniform spectrum. Our discussion is therefore limited to the cross-power, with possible application to the auto-power deferred to future work.

4.4. Adding a template map for the Milky Way

The measurement of Planck $\times$ BOSS quasars is in the challenging regime where $N_{\rm band}\approx N_{\rm 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 Eq. 11 to include a Milky Way template with amplitude $A_{M}$ and noise variance $N_{M}$, as

[TABLE]

Eq. 14 gives

[TABLE]

This is just Eq. 12 where the bright Milky Way variance $C^{uu}_{\ell}$ has been traded with noise-to-signal 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

[TABLE]

If the Milky Way template traces the dust emission poorly $r_{M}\rightarrow 0$, this reverts to the full foreground contamination in Eq. 12. Appendix D describes the extension to many bands.

5. Milky Way template cleaning applied to Planck x BOSS

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 cross 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 signal to noise of the correlated line emission.

5.1. Existing Milky Way templates

There are approaches to use broadband measurements to model thermal dust emission in the Milky Way (Planck Collaboration et al., 2016b, 2014a; Meisner & Finkbeiner, 2015). However, these galactic templates are prone to also contain 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 HI 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 Fig. 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 HI column through $21$ cm radiation (Fig. 8). Both Green et al. (2015) and Lenz et al. (2017) are imperfect tracers of Milky Way dust emission. HI tracers for dust density become biased in high-density regions where ${\rm H}_{2}$ has formed (Lenz et al., 2017), so the HI tracer is best in areas out of the galactic plane considered here. Reddening and HI are imperfect tracers of the dust column depth (Cardelli et al., 1989; Fitzpatrick & Massa, 2007; Schlafly et al., 2016), and dust column is an imperfect tracer of the thermal dust emission.

Fig. 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. Lenz et al. (2017) decreases variance on the largest scales, and Green et al. (2015) adds noise for $\ell>300$.

5.2. Planck and BOSS data

The Planck 353, 545 and 857 GHz maps have full-width at half-max (FWHM) of 4.41, 4.47, and 4.23 arcmin, respectively, and are binned in HEALPix (Górski et al., 2005) pixelization with $N_{\rm side}=2048$. We use a point source mask that is the union of published masks (Planck Collaboration et al., 2016a) across the three frequencies and apodized across $0.5^{\circ}$, and a galactic emission mask with $2^{\circ}$ apodization. In total, this leaves $f_{\rm sky}=0.332$ in the BOSS survey region. Throughout, we will work in MJy/sr units, and convert between measurements in $K_{\rm CMB}$ units to MJy/sr using the mean coefficients 287.45 (353 GHz), 58.04 (545 GHz), and 2.27 (857 GHz) $({\rm MJy/sr})/K_{\rm CMB}$ (Planck Collaboration et al., 2014b).

As a tracer of the underlying overdensity, we use the CORE uniform sample (Ross et al., 2012) of the Baryon Oscillation Spectroscopic Survey (BOSS) (Dawson et al., 2013) spectroscopic quasar sample (Pâris et al., 2017) from Data Release 12 (DR12) (Alam et al., 2015). We develop masks and overdensity as described in Pullen et al. (2018), identical to approaches for quasar clustering (White et al., 2012; Pâris et al., 2014; Eftekharzadeh et al., 2015). For DR12, this procedure yields a catalog of 178,622 quasars. Restricting to mask weight greater than $90\%$ gives 82,522 quasars over $8294\,{\rm deg}^{2}$ and an overlap with the Planck map of $6483\,{\rm deg}^{2}$ with 75,244 quasars.

It is also beneficial to have a galaxy redshift sample that has no cross-correlation with the [C ii] emission, but does trace the CIB and SZ clustered emission, as part of a parametric model for the components. Here we use the BOSS DR12 CMASS luminous red galaxy (LRG) sample (Alam et al., 2017; Reid et al., 2016; Alam et al., 2015), which consists of 862,735 galaxies over $9376\,{\rm deg}^{2}$ and mean $z=0.57$, and is designed to be stellar-mass limited at $z>0.45$. Sample selection here requires spectroscopic sectors (Aihara et al., 2011) with completeness $>70\%$ and redshift completeness $>80\%$. Additionally requiring $0.43<z<0.7$ and pixel coverage $>90\%$ gives 777,202 galaxies over $10,229\,{\rm deg}^{2}$.

In both the quasar and LRG sample, we form the overdensity ${\boldsymbol{\delta}}=({\bf n}-\bar{\bf n})/\bar{\bf n}$ using systematic weights (Anderson et al., 2014). The LRG and CMASS mask regions are additionally apodized by a $0.5^{\circ}$ FWHM kernel.

5.3. Measured cross-powers

Fig. 10 shows the six angular cross-powers used to estimate [C ii] emission and continuum nuisance parameters. These cross-spectra were performed over 3 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-1000$. These spectra are identical to those presented in Pullen et al. (2018).

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 (Fig. 11). Section 5.4 investigates evidence for correlation with large-scale structure.

5.4. Models for CIB and CII 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 $\nu$ and a large-scale structure (LSS) tracer overdensity field, an expression given by

[TABLE]

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

The redshift source distribution of CIB sources can be written as

[TABLE]

where $dn/dL$ is the infrared galaxy luminosity function and $\mathrm{L_{\nu(1+z)}}$ is a model for the CIB luminosity emitted at rest frame frequency $\nu(1+z)$. $b_{\rm 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 3 parameters in our fit to vary: a luminosity amplitude $L_{0}$, a redshift evolution parameter $\delta$, 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 CII}$. 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$-$\ell^{\prime}$ correlations from masked sky regions.

5.5. CIB and CII parameter estimates

We perform our MCMC analysis using CosmoMC (Lewis & Bridle, 2002) for the 6 parameters

[TABLE]

As in Pullen et al. (2018), we fit the parameters using the 6 cross-power spectra between the 3 high-frequency Planck bands and the quasars and LRGs presented in Fig. 10; these are measured in 9 bins over $100<\ell<1000$. We label this fit the “Pullen+2018” fit. We include the full covariance matrix between all the cross-power spectra. Also as before, we fit the mean levels of the CIB emission in the 3 Planck bands (Béthermin et al., 2012) as well as 10 star formation rate 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$\times$quasar clustering model with amplitude $\alpha$.

Fig. 12 shows the parameter fits before and after including Milky Way emission templates. The posteriors for CIB parameters and thermal SZ amplitude do not change significantly from Pullen et al. (2018). The posterior distribution of $A_{\rm CII}$ 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 $\alpha$ for either Milky Way template. The MCMC constrains $-0.092<\alpha_{\rm Lenz}<-0.025$ and $-0.121<\alpha_{\rm Green}<-0.042$ at 95% confidence. This indication of 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.

5.6. Interpreting the MCMC results

We can interpret the lack of impact of Milky Way templates in Planck$\times$BOSS-quasar measurements using the multi-band formalism in Sec. 4. To do this, recast relevant parts of the full nonlinear model as linear modes that fit the correlated continuum emission as in Eq. 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^{\rm dust}_{\nu}=\nu^{\beta}B[T_{\rm dust},(1+z)\nu]$ where $B(T_{\rm dust},\nu)$ is the Planck law, and we fix $\beta=1.5$. To build a linear model, let ${\bf M}$ have a mode which is the best-fit dust with $T_{\rm 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^{\rm dust}_{\nu}/dT(\nu)$.

For the covariance in Eq. 14, use the bandpower covariance of the Planck data. Fig. 14 shows this principal mode of the bandpower covariance, which is consistent with emission from the Milky Way ($\beta=1.5$, $T_{\rm dust}=19$ K from $z=0$). The principal mode is spectrally similar to emission with $T_{\rm 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 Eq. 16, the correlated and uncorrelated continuum terms have similar SED so $f(\nu)\approx g(\nu)$ and the second term of the variance is highly suppressed. The value of a Milky Way template map is diminished in this case.

In the linear model, the Lenz et al. (2017) template increases the signal-to-noise in the line amplitude estimate by only $8\%$. Adding a derivative mode $dI^{\rm dust}_{\nu}/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<500$, where the Milky Way template is best correlated to Planck 545 GHz (Fig. 13). 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 $\times$ 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.

6. A true tomographic survey

Section 4 argued that an intensity survey needs at least an many frequencies as degrees of freedom 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}<k_{\parallel}^{\rm 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 intensity mapping 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 $\nu$ is the integral of comoving specific emission intensity $j(\nu,z)$ across all redshifts in a tomographic survey slice

[TABLE]

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 Eq. 21 over the redshift slice, and so it will depend linearly on the redshift width $\delta z$.

Let the continuum SED at frequency $\nu$ from galaxies at redshift $z$ be $\Theta[\nu(1+z)]$, and take a greybody dust (Blain et al., 2002). Evaluate the SED on the comoving grid points and find its discrete Fourier transform $\tilde{\Theta}(k_{\parallel})$, normalized such that $\tilde{\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

[TABLE]

where $\Delta z$ is the total redshift range, or $N_{\rm chan}\delta z$ for $N_{\rm chan}$ channels, and $f_{\rm 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_{\rm cont}=\frac{2.5}{A_{\rm CII}}$. Here $S_{L}$ is the line brightness times cosmological linear bias $I_{\rm CII}b_{\rm CII}$, 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

[TABLE]

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 intensity mapping 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, $\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$\,{\rm Mpc}/h$ cubes with a pixel size of 1$\,{\rm Mpc}/h$, which corresponds to $2<z<2.5$ in the line-of-sight direction. The continuum SED is then simulated by adding a greybody 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 CII}$ 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}$ (Eq. 23) to our simulated intensity maps. Therefore, the units of the power spectra in Figs. 15 and 16 are ${\rm Mpc}^{3}h^{-3}$. Fig. 15 shows the 2-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.)

Only the two lowest non-zero $k_{\parallel}$ 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_{\perp}$, as a function of $k_{\parallel}$. Figure 17 shows the decrease in the number of modes for $|k|$ by cutting the two lowest non-zero $k_{\parallel}$ modes from the binned one-dimensional power spectrum.

We can also revisit the question from Sec. 4.1 of 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 5 lowest $k_{\parallel}$ bins rejects the majority of 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 5 $k_{\parallel}$ modes must be removed, $50\%$ of the signal is lost. Calculations here have made simple assumptions for the number of degrees of freedom in the correlated continuum SED. In practice, both the correlated and uncorrelated continuum may be more complex through their astrophysics, or through interaction with the instrument.

7. Discussion

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

Correlated continuum emission should be estimated jointly (self-consistently) with bright Milky Way foreground cleaning. For most intensity mapping surveys, a synchrotron or dust continuum dominate 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.

Intensity mapping observations with a limited number of bands have unique challenges that we consider in the case of Planck cross 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, 2) verify that the line amplitude estimate is not biased because of an incomplete continuum model.

We reconsider the analysis of Planck cross 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 cross BOSS will have shortcomings for future intensity mapping 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 $SED(z_{i},\nu_{j})$ and anisotropy structure $C_{\ell}$ of the correlation, similarly 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 intensity mapping 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.

8. Acknowledgements

We acknowledge David Spergel for recommendations for cleaning correlated foregrounds.

Appendix A The multi-tracer approach

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $v_{\ell}$ is the number of modes in the measurement, which is roughly $(2\ell+1)f_{\rm sky}$ for a 2D map on the fraction $f_{\rm 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}^{2}$, but here all factors of $S_{L}$ appear with orders of $N_{g}/C^{\delta\delta}_{\ell}$. To simplify expressions throughout to provide some analytic insight, we will assume that galaxy survey shot noise is negligible, or $N_{g}/C^{\delta\delta}_{\ell}\ll 1$ (taking no terms beyond order 0). In this case, $\sigma^{2}_{S_{L}}=N_{I}/(v_{\ell}C^{\delta\delta}_{\ell})$. This has the simple interpretation as the noise-to-signal 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_{\rm pix}\times N_{\rm pix}$ matrix inverse for $N_{\rm 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 Eq. A1 has access to information from both the cross-power $b_{g}S_{L}C^{\delta\delta}_{\ell}$ and the auto-power $S_{L}^{2}C^{\delta\delta}_{\ell}+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

[TABLE]

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

[TABLE]

where ${\rm 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 Eq. 6.

Marginalizing over $\{C^{\delta\delta}_{\ell},N_{I}\}$, and ignoring galaxy redshift survey shot noise, one recovers the estimate of $S_{L}$ with variance $\sigma^{2}_{S_{L}}=N_{I}/(v_{\ell}C^{\delta\delta}_{\ell})$, as above. Sample variance is absorbed in $C^{\delta\delta}_{\ell}$ as the explanatory variable, allowing $S_{L}$ to be determined independently. Indeed, when the observable is only the cross-power ${\boldsymbol{\mu}}=\{C^{gI}_{\ell}\}$ and the line amplitude $S_{L}$ is the only parameter, the variance $\sigma^{2}_{S_{L}}=(2S_{L}^{2}+N_{I}/C^{\delta\delta}_{\ell})/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^{gI}_{\ell},C^{gg}_{\ell}\}$ here. This provides two equations allowing a fit of two parameters $\{S_{L},C^{\delta\delta}_{\ell}\}$, so assumes that $N_{I}$ is fixed and known. In this case, the variance on $S_{L}$, marginalized over $C^{\delta\delta}_{\ell}$ is again $\sigma^{2}_{S_{L}}=N_{I}/(v_{\ell}C^{\delta\delta}_{\ell})$. 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^{gI}_{\ell}=rb_{g}S_{L}C^{\delta\delta}_{\ell}$. In either the map or power-spectrum space approach, and assuming negligible shot noise, the line amplitude constraint with stochasticity between the tracers becomes

[TABLE]

Note that this does contain cosmic variance terms $\propto S^{2}_{L}$ when $r\neq 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^{\delta\delta}_{\ell})$ 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^{\delta\delta}_{\ell}$ and $S_{L}$ makes $C^{\delta\delta}_{\ell}$ the explanatory variable of the cosmic covariance. $S_{L}$ only appears as an amplitude in the intensity survey so is determined independently of the cosmic covariance.

Appendix B A multiband model with separable foregrounds

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

[TABLE]

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

[TABLE]

Extending this to many bands gives covariance the form

[TABLE]

where ${\bf N}_{\rm th}$ is the diagonal matrix of thermal noise in each band ($N_{A}$ and $N_{B}$ above) and ${\bf 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.

Appendix C Multiband SED reconstruction

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

[TABLE]

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

In the illustrative limit of negligible shot noise, the galaxy redshift survey auto-power fixes $C^{\delta\delta}_{\ell}(z)$, including cosmic variance, making the model for the cross-powers in Eq. C1 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 continuum. Future intensity mapping experiments would be complemented well by galaxy redshift surveys with small shot noise on the scales of interest.

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

[TABLE]

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}$,

[TABLE]

Pack each spectral component into a vector ${\bf g}_{j}|_{i}=g_{j}(\nu_{i})$, and pack 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_{\rm comp}$ correlated continuum modes into ${\bf M}\equiv[{\boldsymbol{\xi}},{\bf g}_{1},...{\bf g}_{N_{\rm 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 Eq. C1 is

[TABLE]

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

[TABLE]

(Note that the data vector defined in Eq. C2 in an estimator for the amplitudes for $\phi$ and $\psi$ in App. B, but applied to many bands. This results in the same form for the covariance as Eq. B3).

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

[TABLE]

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

This linear form recovers our primary results so far. Eq. C6 recovers Eq. 9 in the case where ${\bf M}=[[1,0]^{T}]$ (only $S_{L}$ is fit from the data) and $r_{F}=1$. Here, ${\bf N}_{\rm fg}=C^{uu}_{\ell}[f(\nu_{I}),f(\nu_{V})]^{T}[f(\nu_{I}),f(\nu_{V})]$ and ${\bf N}_{\rm th}={\rm diag}[N_{I},N_{V}]$. Adding one correlated SED component as ${\bf M}=[[1,0]^{T},[g(\nu_{I}),g(\nu_{V})]^{T}]$ recovers Eq. 12, again with $r_{F}=1$.

Appendix D Geometry of the joint correlated and uncorrelated continuum estimation

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

[TABLE]

where

[TABLE]

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

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

If ${\bf M}$ is a complete spectral basis, then $\Delta=0$ regardless of the spectral shape of ${\bf 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 cross BOSS. In this regime,

[TABLE]

These two terms correspond to the thermal noise and foreground components in Eq. 12 with $r_{F}=1$. The second term represents foregrounds that contaminate the line amplitude estimate after marginalizing the correlated continuum spectral modes ${\bf g}_{j}$. An observation with few bands may not provide sufficient degrees of freedom 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_{\rm band}=N_{\rm param}$, Eq. D3 shows that the line amplitude variance always has a term proportional to $C^{uu}_{\ell}$. By adding a template, $C^{uu}_{\ell}\Delta$ in the denominator of Eq. D1 can remain large. In the limit of large foregrounds,

[TABLE]

Now the foreground brightness $C^{uu}_{\ell}$ drops out, and the numerator ${\boldsymbol{\alpha}}_{f}{\boldsymbol{\alpha}}_{f}^{T}$ is suppressed by $\Delta$. 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 Eq. 16 to multiple bands.

Appendix E The impact of stochasticity between Milky Way tracers and dust emission

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

[TABLE]

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

[TABLE]

Appendix F Redshift-space distortions on the Limber-approximated angular power spectrum

In the paper we argue that redshift-space distortions 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 $\beta$ is the RSD parameter given by the growth rate-clustering bias ratio,

[TABLE]

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)

[TABLE]

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

[TABLE]

where we approximate $\ell+1/2\simeq\ell$ to make the following equations less cluttered. Inserting this into Eq. F, we can write

[TABLE]

where

[TABLE]

Now we wish to use this expression to show when RSD can be neglected. In this argument we will assume $\phi(r)$ is a tophat distribution centered at $r=R$ with full width $\Delta R$. In this case, it is evident that if $\ell/k-R\ll\Delta R$, then $\phi[(\ell\pm 2)/k]\simeq\phi(\ell/k)$. Additionally assuming $\ell$ 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 $\phi$s are nonzero when $2/k<\Delta R/2$, which implies the RSD effect vanishes when $\ell\gtrsim 4R/\Delta R$. For the quasars $4R/\Delta R\sim 17$, while for the LRGs $4R/\Delta R\sim 10$. Thus, we are justified in neglecting RSD for our analysis.