Constraining low redshift [C II] Emission by Cross-Correlating FIRAS and BOSS Data

We cross-correlate FIRAS data with cosmological overdensity inferred from the BOSS spectroscopic galaxy redshift survey to search for extragalactic [C II] emission. Two competing effects make the BOSS CMASS ($0.41<z<0.75$) and LOWZ ($0.06<z<0.49$) especially well-tuned for this goal. The growth of large-scale structure and geometric factors yield increasing density contrast on large angular scales where FIRAS is most sensitive. However, by $z<0.25$ ($\nu>1500$ GHz), the FIRAS noise rises dramatically. For our analysis, we bin the CMASS and LOWZ galaxy maps into redshift slices that correspond to the FIRAS frequency channels for the [C II] line. We then use the technique of spherical harmonic tomography (SHT) (see, for example, Asorey et al. (2012); Nicola et al. (2014); Lanusse et al. (2015)) to compute the angular cross-power spectrum, $C_{\ell}^{\times}(z,z^{\prime})$, between the FIRAS maps and BOSS galaxy over-densities at each pair of redshifts.

With a fine enough binning in redshift, SHT captures the complete information available in the power spectrum (Asorey et al., 2012), and it is well-matched to large area surveys such as BOSS, where a flat-sky approximation would introduce significant distortions. In order to fit our data to a model, the expected angular clustering signal, $C_{\ell}^{\delta}(z,z^{\prime})$, must be computed by integrating a 3D power spectrum model over highly oscillatory Bessel functions (details included in Section 4.1). Standard anisotropy codes, such as CAMB (Lewis & Challinor, 2011) and CLASS (Di Dio et al., 2013, 2014), include methods for performing this integration (with or without the Limber approximation).

We describe the analysis over several sections. Section 2 motivates the $C_{\ell}(z,z^{\prime})$ statistic, explains how to use the pseudo-$C_{\ell}$ technique to compute $C_{\ell}(z,z^{\prime})$ for surveys with partial sky coverage, and describes the likelihood function we use for parameter estimation with $C_{\ell}(z,z^{\prime})$. Section 3 describes the FIRAS dataset, viewed as a [C II] intensity map, and the BOSS CMASS and LOWZ data sets, viewed as binned galaxy over-density maps. Section 4 describes our dark matter model and its use in fitting the [C II] and CIB amplitude to our FIRAS${\times}$BOSS data. It also describes how we model the FIRAS${\times}$BOSS covariance by fitting parametric models to the BOSS${\times}$BOSS and FIRAS${\times}$FIRAS data. Section 5 discusses our [C II] constraints and how they relate to other measurements and astrophysical models. We conclude in Section 6.

We perform our analysis with Spherical Harmonic Tomography (SHT), a two-point statistic wherein each redshift slice of data is decomposed into spherical harmonics, and the angular power spectrum, $C_{\ell}(z,z^{\prime})$, is calculated between each pair of redshift slices. SHT captures the complete information available in the more typically used three-dimensional power spectrum statistic, $P(k)$ (Asorey et al., 2012). The SHT technique has already been used in several clustering analyses and forecasts for galaxy surveys (e.g. Thomas et al. (2011); Leistedt et al. (2013); Balaguera-Antolínez et al. (2018); Loureiro et al. (2019a); Xavier et al. (2019); Fonseca et al. (2019); Nicola et al. (2020); Viljoen et al. (2020); Tanidis & Camera (2021)). Among these, a recent analysis of BOSS CMASS and LOWZ clustering using Spherical Harmonic Tomography (Loureiro et al., 2019a) found equivalent or better constraints on cosmological parameters compared to standard power spectrum analysis techniques. As the set of two-point cross-correlations, the SHT contains all information in the data cubes that is statistically isotropic and Gaussian.

SHT has some inherent geometrical advantages, especially for large-angle and deep surveys. A significant advantage is that the spherical coordinates apply to wide-angle surveys like BOSS without any flat-sky approximation. A traditional $P(k)$ analysis relies on the flat-sky approximation to distinguish between transverse and line-of-sight modes, which is critical for accurately representing redshift space distortions (RSDs) and distinguishing continuum foregrounds from line signal. By contrast, in the SHT formalism, both foregrounds and linear redshift space distortions take exact, simple forms. This geometric advantage is shared with the related analysis technique of spherical Fourier-Bessel (SFB) decomposition (Fisher et al., 1995; Heavens & Taylor, 1995; Rassat & Refregier, 2012; Leistedt et al., 2012; Yoo & Desjacques, 2013; Lanusse et al., 2015; Liu et al., 2016; Grasshorn Gebhardt & Doré, 2021), which, in addtion to spherical harmonic decomposition, involves a further Fourier-Bessel transform and data-compression along the line-of-sight. However, SFB does not share another important feature of SHT, which is its ability to capture redshift-dependent change over cosmological time in deep surveys. For deep surveys, structure growth and changes in star formation rate break the assumption of translational invariance in the line-of-sight direction, rendering the $P(k)$ or SFB statistic insufficient. However, since $C_{\ell}(z,z^{\prime})$ does not compress data along the line-of-sight direction, it describes redshift evolution. A study (Mondal et al., 2022) of simulated 21-cm data from the Epoch of Reionization (EoR) found that an implementation of the SHT technique, MAPS (Mondal et al., 2018; Mondal et al., 2019; Mondal et al., 2020), obtains $\sim$2 times more stringent error bars on model parameters than techniques that fail to capture redshift evolution due to data compression along the line-of-sight, such as $P(k)$ or SFB. A final geometric advantage is that $C_{\ell}(z,z^{\prime})$ describes the data in observing coordinates of angle and frequency (or, equivalently, redshift) rather than re-gridding the data onto cosmological distances in an assumed cosmological model. An MCMC likelihood analysis that constrains cosmological parameters would therefore not need to recompute the data statistic at each step, which in principle would be needed for a $P(k)$ or SFB analysis.

There are several practical advantages to parameter estimation with SHT. Due to the scan strategy, intensity mapping data generally have inhomogeneous noise and partial sky coverage in the angular direction. Multiplication by noise weights in real space couples transverse modes in Fourier space, which can be easily accounted for in SHT analysis, thanks to pre-existing work (Hivon et al., 2002; Tristram et al., 2005) on the pseudo-$C_{\ell}$ technique from CMB analysis. Intensity maps will also often have significant variation in the noise level in the frequency (line-of-sight) direction because of variations in spectrometer noise or chromatic contamination such as terrestrial radio frequency interference. Inhomogeneous noise in real space produces coupling of line-of-sight Fourier modes in analyses of $P(k)$. In contrast, because $C_{\ell}(z,z^{\prime})$ does not perform any Fourier or Bessel transform in the redshift direction, no coupling occurs, and the noise can be expressed as a simple function of the redshift slice. Additionally, chromatic beam effects can be modeled per redshift slice without introducing any flat-sky approximation.

In addition, SHT provides an avenue for self-consistently handling foregrounds and other continuum signals. In contrast, other approaches that remove foregrounds in map space before the two-point analysis must contend with signal loss (Switzer et al., 2015; Cheng et al., 2018). With SHT, the covariance of $C_{\ell}(z,z^{\prime})$ can be constructed to include information about foregrounds. The inverse covariance weights the data in the likelihood analysis for the line brightness, suppressing modes contaminated by foregrounds self-consistently with the parameter estimation. Additionally, the $C_{\ell}(z,z^{\prime})$ model can include signal terms such as the cross-correlation of a galaxy redshift sample with the correlated continuum emission of the galaxies in the IM volume (Serra et al., 2014a; Pullen et al., 2018; Switzer et al., 2019).

SHT also introduces several new challenges to the analysis. First, it places significant requirements on memory for the computation of the likelihood. This is because the size of the covariance scales as $N_{b}^{2}N_{z}^{4}$, where $N_{z}$ is the number of redshift bins and $N_{b}$ is the number of angular bins. For FIRAS, this is fairly manageable because we analyze only three angular bins and about 15 redshift slices for each of FIRAS${\times}$CMASS and FIRAS${\times}$LOWZ. The size of the covariance will be more challenging for future instruments with higher spectral and angular resolution.

The visualization of $C_{\ell}(z,z^{\prime})$ and its errors presents an additional challenge. The extra dimensionality of $C_{\ell}(z,z^{\prime})$, which allows it to measure redshift evolution, also means that angular and redshift information cannot be shown in a single plot. We develop several approaches for displaying the data and describing the goodness of fit to high-dimensional data with complex covariance. Lastly, the evaluation of the $C_{\ell}(z,z^{\prime})$ model is computationally expensive. Since the current generation of IM observations focuses on detection and measurements of line brightness (Kovetz et al., 2019), this is not an issue. In this regime, the $C_{\ell}(z,z^{\prime})$ model can be constructed from linear combinations of cosmological clustering and shot noise templates that need only be calculated once. In contrast, studies of large-scale structure acquire information from the nonlinear dependence of $C_{\ell}(z,z^{\prime})$ on cosmological parameters, so can require expensive recalculation of the anisotropy. Recent work has accelerated the integrals that convert $P(k,z)$ to $C_{\ell}(z,z^{\prime})$ without using the Limber approximation (Campagne et al., 2017; Schöneberg et al., 2018), making cosmological parameter estimation with $C_{\ell}(z,z^{\prime})$ more practical.

Extensive work in CMB data analysis has developed an approach (Hivon et al., 2002) for dealing with incomplete sky coverage and an approximate formula for the covariance it induces between angular scales (Tristram et al., 2005). The code package NaMaster (Alonso et al., 2019; Alonso, 2018) includes this full functionality, along with contaminant removal for polarized and unpolarized maps, on both large curved-sky maps and small maps, using a flat-sky approximation.

For observed maps A and B with known, possibly different, inverse noise weights and beams, the angular power spectrum of the inverse noise weighted maps is the pseudo-$C_{\ell}$ spectrum (Hivon et al., 2002), which we label $D_{\ell}$, following the notation of Tristram et al. (2005). The pseudo-$C_{\ell}$ spectrum is related to the true full-sky angular power spectrum $C_{\ell}$ by

where $\omega^{A}_{\ell^{\prime}}$ is the product of the beam and pixel window function of map A, $\omega^{B}_{\ell^{\prime}}$ is the product of the beam and pixel window function of map B, and $M^{AB}_{\ell\ell^{\prime}}$ is the mixing matrix, computed via

where $\mathcal{W}^{AB}_{\ell^{\prime\prime}}$ is the angular power spectrum of the inverse noise spatial weights of maps $A$ and $B$, and $\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\ 0&0&0\end{pmatrix}$ represents a Wigner 3-j symbol. These formulas work for both auto and cross-correlations. The mixing matrix is not invertible with partial sky coverage, and it is convenient to define the binned pseudo-$C_{\ell}$ spectrum and binned mixing matrix. Here, a range of multipoles $\ell$ are combined into a bandpower $b$, for which the mixing matrix of bandpowers $b$ and $b^{\prime}$ is invertible using suitable binning. Define $B_{b\ell}$ as a binning matrix that averages blocks of consecutive $\ell$s into bins indexed by $b$ and define $B^{u}_{\ell b}$ as an equivalent unbinning matrix that assigns the average value of a bin $b$ to each $\ell$ within that bin. Then, define the binned quantities

If the bin size is wide enough, $M^{AB}_{bb^{\prime}}$ is invertible even for partial sky coverage. One can then form an estimate, $\hat{C}_{b}$, of the binned angular power spectrum from the observed binned pseudo-$C_{b}$ spectrum, $\hat{D}_{b}$, via

This procedure should only recover the true angular power spectrum if $C_{\ell}$ is constant within each bin $b$. Since this is not generally the case, to compare $\hat{C}_{b}$ to a model $C_{\ell}$, as we do in Section 4, we always compute or simulate the expected binned angular power spectrum, $\tilde{C}_{b}$, which is related to the model via

If one assumes that the true $C_{\ell}(z,z^{\prime})$ distribution is Gaussian, then the $b=b^{\prime}$ terms of the covariance can be roughly approximated as

where $f_{\rm sky}$ is the fraction of the sky seen by the survey and $\Delta\ell$ is the width of the bin centered at $\ell$. For intensity mapping tomography, the first term is the product of the autopower of the IM data $C_{b}^{\rm IM}$ and the galaxy survey $C_{b}^{g}$, and the second term is the product of cross-powers $C_{b}^{\times}$ between the IM data and the galaxy survey. For the FIRAS${\times}$BOSS analysis, the first term dominates the covariance due to the high thermal noise and large (at low $\ell$) Milky Way foreground signal present in $C_{b}^{\rm IM}(z,z^{\prime})$. So, as a rough rule of thumb, the $b=b^{\prime}$ terms of the cross-power covariance are

Appendix A applies this equation to derive an approximate formula for the expected sensitivity that a cross-power survey could achieve on the line intensity. The covariance can be more accurately approximated under the assumption of large sky coverage (Tristram et al. (2005), formula included in Appendix C) or computed via simulated draws from an assumed $C_{\ell}(z,z^{\prime})$ model and repeated pseudo-$C_{\ell}$ computation of the resulting $C_{b}(z,z^{\prime})$. For our BOSS${\times}$BOSS analysis, we use the approximate covariance of Tristram et al. (2005), as it matches our simulations well. For the FIRAS${\times}$FIRAS and FIRAS${\times}$BOSS analysis, the approximation fails to match simulations, so we instead use a fully simulated covariance. We suspect that the inconsistency between the simulations and the approximate formula is due to the combination of small sky coverage and the steep angular index of the Milky Way foregrounds. Further details on the covariance used for each analysis are included in Section 4.3.

For a range of binned multipoles indexed by bandpower $b$, the computed $C_{b}(z,z^{\prime})$ is a rank-3 tensor, indexed by $b$, $z$, and $z^{\prime}$. To perform a likelihood analysis, we define $\mathbf{x}$ as a flattened vector of all the unique elements of $C_{b}(z,z^{\prime})$.

Note that for auto-powers, there are only $N_{b}\times N_{z}\times(N_{z}+1)/2$ unique elements, since $C_{b}^{A\times A}(z,z^{\prime})=C_{b}^{A\times A}(z^{\prime},z)$. For the cross-power, all elements are unique, and $\mathbf{x}$ has a length of $N_{b}\times N_{z}^{2}$.

The likelihood formed from $C_{b}(z,z^{\prime})$ is well-studied in CMB literature (see e.g., Hamimeche & Lewis (2008) for a survey of likelihood forms depending on assumptions). If the amplitudes of the spherical harmonics $a_{\ell m}(z)$ of the maps are drawn from an underlying Gaussian $C_{\ell}(z,z^{\prime})$ distribution, then the resulting measured vector of anisotropies, $\mathbf{\hat{x}}$, is Wishart-distributed. With our bin size of $\Delta\ell=9$, and for the $\ell$-range with signal sensitivity (see Figure 7), the number of modes contributing to each bin is high enough that the Wishart distribution is reasonably well-approximated by a Gaussian distribution

where $\bf{\Sigma}(\Theta)$ is the bandpower covariance matrix of the flattened data vector and includes binned angular $b{-}b^{\prime}$ coupling induced by incomplete sky coverage, $\hat{\bf{x}}$ is a flattened vector of the $\hat{C}_{b}(z,z^{\prime})$ estimate computed from the data, and $\bf{x}(\Theta)$ is a flattened vector of the $C_{b}(z,z^{\prime})$ model, which depends on the science parameters of interest, $\bf{\Theta}$.

When analyzing multiple datasets, such as galaxy surveys and intensity maps, the likelihood can be expanded to apply to a data vector that concatenates the auto- and cross-powers of all the datasets. So, for the FIRAS and BOSS data, the data vector $\hat{\bf{x}}$ could contain flattened FIRAS${\times}$FIRAS, BOSS${\times}$BOSS, and FIRAS${\times}$BOSS $C_{b}(z,z^{\prime})$ estimates. If the thermal noise and galaxy shot noise were small, then the high correlation between the galaxy data and [C II] data, which trace the same matter fluctuations, could be exploited to remove cosmic variance from the measurement of the [C II] line intensity (McDonald & Seljak, 2009; Bull et al., 2015; Abramo et al., 2016; Switzer, 2017b; Switzer et al., 2019; Oxholm & Switzer, 2021). For this FIRAS${\times}$BOSS analysis, thermal noise and foregrounds, rather than cosmic variance, are dominant sources of error, so the benefits of such an approach are negligible. However, future intensity mapping experiments may benefit from this approach. For simplicity, we focus on [C II] constraints from the cross-power, FIRAS${\times}$BOSS, and we only use FIRAS${\times}$FIRAS and BOSS${\times}$BOSS to validate the cross-power covariance model.

FIRAS is a rapid-scan polarizing Michelson interferometer (Mather et al., 1993) on the COBE satellite that mapped the frequency spectrum of the full infrared sky at a coarse angular resolution. The frequency spectrum for each pointing was obtained via an inverse Fourier transform of the interferogram of measured powers over a discrete range of instrument path length differences. The resulting measurements of the sky spectrum are equal to the true spectrum convolved by the inverse Fourier transform of an apodization function (Fixsen et al., 1994). Figure 1 plots this frequency response function, which we shall denote $A(\Delta\nu)$.

The published FIRAS maps are binned in the HEALPix ¹¹1http://healpix.sourceforge.net (Górski et al., 2005; Zonca et al., 2019) format with resolution parameter $N_{\rm side}{=}16$, corresponding to 3072 angular pixels, sufficient to sample the 7-degree beam. In addition to the sky maps, inverse-noise weight maps were produced based on fluctuations of the different interferograms contributing to each pixel. We upgrade the map binning to $N_{\rm side}{=}128$ for analysis. This regridding does not gain any angular information from the FIRAS maps, but it allows finer, more accurate noise weights for the galaxy over-density maps. Figure 2 shows the inverse-noise weights at $N_{\rm side}=128$ for [C II] emission from $z\sim 0.52$.

The FIRAS beam is formed by a non-imaging parabolic concentrator, which creates a near-tophat beam response with 7-degree FWHM. This beam was measured by in-flight scans of the Moon (Figure 3). In addition to this intrinsic beam convolution, the finite time required to complete an interferogram combined with the FIRAS scan motion causes the maps to be further smoothed in the ecliptic scan direction by a 2.4-degree tophat function. We account for all beam, scan, and pixelization smoothing effects on the power spectrum via simulation: realizations of a model power spectrum are drawn at $N_{\rm side}=128$, convolved by the FIRAS beam model, convolved by a 2.4-degree tophat function in the ecliptic direction, degraded to $N_{\rm side}=16$, and re-grid onto $N_{\rm side}=128$. We then compute an angular transfer function in $\ell$ by calculating the ratio of the angular power spectrum calculated from these modified maps to the original angular power spectrum. The square root of this transfer function, denoted $\omega_{FIR}(\ell)$ (Figure 4), represents the combined beam, scan, and pixel window function for the FIRAS maps.

The BOSS survey (Dawson et al., 2012) extends the Sloan Digital Sky Survey (SDSS). It was designed to measure the Baryon Acoustic Oscillation feature in the galaxy matter power spectrum. Precise spectroscopic redshifts were obtained for approximately 1.5 million galaxies in the redshift range $0{<}z{<}0.8$, selected to have approximately constant stellar mass. Details about the telescope and instruments of SDSS can be found in Fukugita et al. (1996), Gunn et al. (1998), Gunn et al. (2006), Doi et al. (2010), and Smee et al. (2013).

BOSS data release 12 (Alam et al., 2015) includes 100 mock unclustered realizations of CMASS and LOWZ galaxies. We bin both the real catalogs and the unclustered mocks onto HealPix maps with $N_{\rm side}=128$, in redshift bins corresponding precisely to the FIRAS frequency bins under the assumption that the FIRAS signal is redshifted [C II] emission. We construct the CMASS and LOWZ galaxy selection functions, denoted $\bar{n}(z,\theta)$, by averaging their mock catalogs, assuming a selection function that is separable in angle and redshift. The separability assumption reduces shot noise in the selection function and should be sufficiently accurate for the cross-correlation, limited by FIRAS noise. We define the boundary of the survey by zeroing pixel weights where the selection function is more than $1.3$ standard deviations below the angular average of the sample, which additionally de-weights some regions with lower coverage. The galaxy over-density field is then formed via

where $n(z,\theta)$ denotes the binned galaxy maps and $\bar{n}(z,\theta)$ denotes the selection function. Figure 6 displays the galaxy over-density maps and selection function for a single redshift slice at $z\sim 0.52$.

Previous analysis of the BOSS data has found evidence of systematic contamination from stellar populations at low $\ell$ (Loureiro et al., 2019a). We also find evidence of systematic low $\ell$ contamination, with spurious $z-z^{\prime}$ correlations visible at the lowest $\ell$-bin of our analysis ($2\leq\ell\leq 10$), indicative of a common systematic component across redshifts. $\chi^{2}$ tests of BOSS data compared to our fitted model are high when including this bin. However, they drop to expected values when excluding this bin. We cut this lowest angular bin and the next bin ($11\leq\ell\leq 19$) from our analysis (see Figure 7).

Both our BOSS${\times}$BOSS and FIRAS${\times}$BOSS models require the dark matter angular power spectrum of the overdensity field, $C_{\ell}^{\delta}(z,z^{\prime})$. We calculate this angular power spectrum with the Boltzmann code CLASS (Di Dio et al., 2013, 2014), using cosmological parameters inferred from the Planck 2015 (Ade et al., 2016) temperature and low-$\ell$ polarization maps (TT+LowP). The Halofit routine (Smith et al., 2003) provides nonlinear corrections to the power spectrum. However, on the several-degree scales of this analysis, the fluctuations are well-described by linear perturbation theory ($k_{\rm max}\sim\frac{\ell_{\rm max}}{\chi(z)_{\rm min}}\sim\frac{50}{880}$ h/Mpc $\sim 0.06$ h/Mpc), and nonlinear corrections are small. CLASS computes the angular power spectrum from the 3D power spectrum, $P(k)$, according to the equation

where $P^{\delta}(k,z{=}0)$ is the dark matter power spectrum at the current epoch, and, if there are no redshift space distortions (RSDs),

where $b_{A}$ is the bias for dark matter tracer $A$, $\phi_{z}(z^{\prime\prime})$ is a tophat redshift selection function that is non-zero only over the range of the redshift slice centered at redshift $z$ and normalized to integrate to 1, $j_{\ell}$ is a spherical Bessel function of the first kind with parameter $\ell$, $G(z^{\prime\prime},k)$ is the growth factor, and $\chi(z^{\prime\prime})$ is the radial comoving distance to the shell at redshift $z^{\prime\prime}$. Linear redshift space distortions (RSDs) can be included (Fisher et al., 1994; Padmanabhan et al., 2007) by replacing $W_{A}(k,z)$ with $W_{A}^{\rm tot}(k,z)$, where

where $\beta_{A}(z)=f(z)/b_{A}$, with $f(z)$ being the logarithmic growth rate of linear perturbations in the matter power spectrum. The $C_{\ell}^{A\times B}(z,z^{\prime})$ that results from using $W_{A}^{\rm tot}(k,z)$ and $W_{B}^{\rm tot}(k,z)$ contains cosmological terms proportional to $b_{A}b_{B}$, proportional to $(b_{A}{+}b_{B})/2$, and independent of both $b_{A}$ and $b_{B}$. We label these terms $C_{\ell}^{(2)}(z,z^{\prime})$, $C_{\ell}^{(1)}(z,z^{\prime})$, and $C_{\ell}^{(0)}(z,z^{\prime})$ respectively. They are calculated from CLASS via linear combinations of $C_{\ell}(z,z^{\prime})$ computations without RSD and bias 1, with RSD and bias 1, and with RSD and bias 0. With this formalism, the cross-power spectrum of two biased matter tracers is given by

The bias dependence of these terms is reminiscent of the Kaiser correction in power spectrum space. Indeed, these equations can be derived by including the Kaiser enhancement term in a plane-wave expansion of the power spectrum and integrating along the line-of-sight (Padmanabhan et al., 2007).

We do not model finger-of-God effects or redshift smearing due to spectroscopic survey errors since our redshift bin size of $\Delta z\sim 0.02$ is more than 10 times larger than the 400 km/s satellite galaxy velocity dispersion fit and 30 km/s spectroscopic error fit found by Guo et al. (2015)’s analysis of CMASS galaxies. Analyses with sufficient redshift resolution to resolve the finger-of-God effect can include it in the signal model by smoothing the radial window function in the $C_{\ell}(z,z^{\prime})$ calculation (Loureiro et al., 2019b; Grasshorn Gebhardt & Jeong, 2020).

The cross-correlation signal model consists solely of correlated continuum emission (dust) and [C II] emission from the BOSS galaxies. Because they are uncorrelated with the cosmological overdensity field, the thermal noise and foregrounds contribute zero average cross-power and thus will not factor into the mean signal model (although the variance caused by their spurious correlation with the galaxy survey will be included through the FIRAS auto-power in the covariance).

In Appendix B, we derive the functional form of our cross-power model. The [C II] part of the model is

The redshift dependence of $I_{\rm[C{\sc\,II}]}(z)$ is shown in Equation 36. An intensity mapping experiment with sufficient sensitivity can fit the parameters that control the redshift evolution of $I_{\rm[C{\sc\,II}]}(z)$. However, due to the high noise and large beam of FIRAS, we fix all of those evolution parameters with reasonable values from Pullen et al. (2018) (see details in Appendix B). Those values lead to a modest evolution in which the brightness increases ${\approx}20\%$ toward higher redshift over each of the CMASS and LOWZ redshift ranges. Our MCMC analysis assumes this redshift shape and fits for the overall amplitude of $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}(z{=}z_{\rm center})$ at the central redshift of each region (CMASS and LOWZ, respectively).

The CIB portion of the cross-power model is

where $\frac{dI_{\rm CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta z^{\prime\prime}$ is the intensity of the CIB that is emitted from sources in a redshift bin of size $\Delta z^{\prime\prime}$, centered at redshift $z^{\prime\prime}$, and measured at a frequency of $\nu_{\rm[C{\sc\,II}]}/(1+z)$. The sum over $z^{\prime\prime}$ in Equation 16 should, in principle, be carried out over all redshifts, even those outside of the galaxy survey. In practice, for the redshift bins and $\ell$ bins considered in this analysis, the bracketed $C_{\ell}(z^{\prime\prime},z^{\prime})$ kernel is dominated by the $z^{\prime\prime}=z^{\prime}$ term, with the neighboring off-diagonal terms around 20% of the diagonal term, and the rest of the terms are negligible. As with our [C II] analysis, because of the limited sensitivity of the FIRAS data, we do not attempt to constrain the redshift evolution of the CIB brightness or the spectral shape of the CIB emission. Instead, these are fixed by the assumed values of $\beta=1.5$ (Planck Collaboration et al., 2014b), $T_{d}=26$ K (Serra et al., 2014a), $\Phi(z)=(1+z)^{2.3}$ (Pullen et al., 2018), $\log_{10}(M_{\rm eff}/M_{\odot})=12.6$ (Planck Collaboration et al., 2014b; Serra et al., 2016), and $\sigma^{2}_{L/M}=0.5$ (Shang et al., 2012; Planck Collaboration et al., 2014b; Serra et al., 2016). Because of this assumed spectral and redshift evolution, the parameter we constrain is $\frac{dI_{\rm CIB}(\nu_{\rm[C{\sc\,II}]}^{z}=\nu_{\rm center},z^{\prime\prime}=z_{\rm center})}{dz^{\prime\prime}}$, where $z_{\rm center}$ is the central redshift of the analysis region, and $\nu_{\rm center}$ is the corresponding central frequency of the analysis region.

Finally, our signal model must account for the FIRAS data being convolved by the FIRAS frequency response function, $A(\nu)$. To do this, we convert the frequency response function $A(\nu)$ to a function of redshift, $A(z)$, and convolve $C_{\ell}^{\times}(z,z^{\prime})$ by $A(z)$, resulting in the final signal model

In our MCMC analysis, we fix the value of the galaxy bias, $b_{g}$, to the best-fit value from our BOSS${\times}$BOSS analysis (Section 4.3.2) and also fix the [C II] and CIB bias to be identical to the galaxy bias ($b_{\rm[C{\sc\,II}]}=b_{g}$). This fixing of the [C II]/CIB bias has almost no effect on our fits since only the small RSD terms ($C^{1}_{\ell}(z,z^{\prime})$ and $C^{0}_{\ell}(z,z^{\prime})$) can break the degeneracy between $b_{\rm[C{\sc\,II}]}$ and $I_{\rm[C{\sc\,II}]}$, and FIRAS${\times}$BOSS lacks the precision to break that degeneracy. The next section, 4.3 develops the covariance model employed by the MCMC.

Figure 8 shows the MCMC contours for our parameter fits to the FIRAS${\times}$BOSS data with the CMASS and LOWZ galaxies. Also shown, on the bottom row, are the $\chi^{2}$ of the CMASS and LOWZ maximum likelihood fits, compared to a simulated distribution of best-fit $\chi^{2}$ values. For the MCMC analysis, we apply a simple flat prior that restricts both $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ and $b_{\rm[C{\sc\,II}]}dI_{\rm CIB}/dz$ to positive values. This prior has almost no effect on our best-fit $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ values, but since we do not have enough sensitivity for detection, it prevents unphysical negative values from counting towards our quoted upper limit constraints. We find, at 95 percent confidence, that $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}<0.31$  MJy/sr at $z\sim 0.35$ and $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}<0.28$  MJy/sr at $z\sim 0.57$.

The required pieces for the covariance used in the MCMC parameter estimation are models for: 1) the [C II] and CIB signal associated with cosmological clustering, 2) the BOSS signal associated with cosmological clustering, plus shot noise, and 3) the FIRAS thermal noise and foregrounds. In the following three subsections, we describe each of these models in turn. Appendix D describes our method of simulating the covariance from these three models.