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

C. J. Anderson1, E. R. Switzer1, P. C. Breysse2
1NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
2Center for Cosmology and Particle Physics, Department of Physics, New York
University, 726 Broadway, New York, NY, 10003, U.S.A E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We perform a tomographic cross-correlation analysis of archival FIRAS data and
the BOSS galaxy redshift survey to constrain the amplitude of [C II]
${}^{2}P_{3/2}\rightarrow$ ${}^{2}P_{1/2}$ fine structure emission. Our
analysis employs spherical harmonic tomography (SHT), which is based on the
angular cross-power spectrum between FIRAS maps and BOSS galaxy over-densities
at each pair of redshift bins, over a redshift range of $0.24<z<0.69$. We
develop the SHT approach for intensity mapping, where it has several
advantages over existing power spectral estimators. Our analysis constrains
the product of the [C II] bias and [C II] specific intensity,
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$, to be $<0.31$ MJy/sr at
$z{\approx}0.35$ and $<0.28$ MJy/sr at $z{\approx}0.57$ at $95\%$ confidence.
These limits are consistent with most current models of the [C II] signal, as
well as with higher-redshift [C II] cross-power spectrum measurements from the
Planck satellite and BOSS quasars. We also show that our analysis, if applied
to data from a more sensitive instrument such as the proposed PIXIE satellite,
can detect pessimistic [C II] models at high significance.

###### keywords:

galaxies: evolution – methods: data analysis – infrared: general

††pubyear: 2022††pagerange: Constraining low redshift [C II] Emission by
Cross-Correlating FIRAS and BOSS Data–References

##  1 Introduction

Atomic and molecular line emissions are powerful diagnostic tracers of the
evolution of star formation over cosmic time. For instance, line emission
studies may help reveal the causes of the significant decline in star
formation rate after its peak at $z{\sim}2$, inferred from observations of
optical and infrared continuum radiation (Madau & Dickinson, 2014). One of the
most useful lines for understanding the context of star formation is the [C
II] ${}^{2}P_{3/2}\rightarrow$ ${}^{2}P_{1/2}$ fine structure line at $1901$
GHz, which is the brightest cooling line in the far-infrared (FIR) spectrum,
typically accounting for 0.1$\%$ to 1$\%$ of FIR energy (Stacey et al., 1991;
Malhotra et al., 1997; Díaz-Santos et al., 2017). The bulk of this [C II]
emission is expected to come from photodissociation regions (PDRs) on the
edges of molecular gas clouds (Stacey et al., 2010). This association suggests
that [C II] emission can trace the molecular gas available to fuel star
formation. A study of $z\sim 2$ galaxies detected by ALMA bore this out,
finding that [C II] emission is linearly related to molecular gas content
(Zanella et al., 2018). Measurements of [C II] emission and other lines (such
as [C I] and the rotational transitions of CO) can be used in concert with
galaxy formation models to constrain the physical properties of star-forming
regions within galaxies (Popping et al., 2019; Yang et al., 2021).

Studies of line emission from individual galaxies (e.g., Carilli & Walter
(2013); Hemmati et al. (2017); Zanella et al. (2018)) provide important
insights into the [C II] luminosity function. However, they are subject to
sample variance in small survey regions and are limited to galaxies above a
brightness threshold, so they may not capture the cosmic average of the
conditions of star formation. A complementary technique to individual galaxy
studies is the technique of intensity mapping, which aims to map large-scale
structure by detecting the aggregate redshifted line emission without
cataloging individual sources (Hogan & Rees, 1979; Scott & Rees, 1990; Madau
et al., 1997; Suginohara et al., 1999; Wyithe et al., 2008; Chang et al.,
2008; Visbal et al., 2011; Kovetz et al., 2017, 2019). Some of the advantages
of intensity mapping are that it captures all sources of emission rather than
a biased sample of only the brightest sources. It also puts significantly
lower requirements on telescope size, as the angular resolution need not be
sufficient for individual source detection. Similarly, because individual
objects do not need to be detected at a high signal-to-noise ratio, intensity
mapping surveys can quickly cover large cosmic volumes, providing a complete
census of emitting gas.

Intensity mapping has developed rapidly in the past fifteen years. Pathfinding
intensity mapping surveys have used pre-existing telescopes to detect
aggregate emission from the 21-cm line of neutral hydrogen (HI) (Pen et al.,
2009; Chang et al., 2010; Masui et al., 2013; Switzer et al., 2013; Wolz et
al., 2017; Anderson et al., 2018; Wolz et al., 2021) via cross-correlation
with optical galaxy surveys. Several dedicated 21-cm intensity mapping
experiments are now underway, targeting both the epoch of reionization (e.g.,
LOFAR, van Haarlem et al. (2013), and SKA precursors HERA and MWA, DeBoer et
al. (2017); McKinley et al. (2018)) as well as the era of dark energy
dominance (e.g., CHIME, Bandura et al. (2014), Tianlai, Chen (2012), HIRAX,
Crichton et al. (2021), SKA precursor MeerKAT, Wang et al. (2021), and BINGO,
Abdalla et al. (2021)). The initial focus on the 21-cm line has expanded to
include [C II], the rotational lines of CO, Lyman $\alpha$, H$\alpha$,
H$\beta$, and more (Kovetz et al., 2019).

The COPPS survey has utilized the Sunyaev-Zel’dovich Array (SZA) to make
tentative detections of CO emission at $z\sim 2.6$ in cross-correlation with
spectroscopic galaxy surveys (Keenan et al., 2021) and auto-correlation in the
shot-noise regime (Keating et al., 2016). A similar CO shot-noise detection
was made with Millimeter-wave Intensity Mapping Experiment (mmIME), using ALMA
and ACA facilities (Keating et al., 2020). Among a new generation of dedicated
ground-based intensity mapping instruments targeting CO and [C II] are COMAP
(Cleary et al., 2021), designed to measure CO(1-0) from $2.4<z<3.4$ and
CO(2-1) at $z=6-8$ , TIME (Crites et al., 2014), targeting [C II] emission
from the epoch of reionization and CO from $0<z<2$, CONCERTO (The CONCERTO
collaboration et al., 2020), focusing on [C II] emission from the epoch of
reionization, and the CCAT-prime receiver for the FYST (CCAT-Prime
collaboration et al., 2021), which targets [C II] at $3.5<z<8$ and [OIII] at
$z>7$. Two NASA balloon experiments designed to measure line emission in the
far-infrared (FIR) are TIM (Vieira et al., 2020), focusing on [C II] at
$0.52<z<1.67$, and EXCLAIM (Cataldo et al., 2021), focusing on CO and [C II]
in windows from $0<z<3.5$. The NASA MIDEX-class satellite mission, SPHEREx,
will produce intensity maps of multiple lines, including H$\alpha$, H$\beta$,
[OII], and [OIII] at $z<5$ (Gong et al., 2017).

The brightness of the [C II] line makes it an excellent candidate for the
technique of line intensity mapping. Pullen et al. (2018) recently
demonstrated the promise of [C II] intensity mapping through an analysis of
the cross-correlation between the angular distribution of high-redshift
($z{\sim}2.6$) quasars in the BOSS survey and the $353$, $545$, and $857$ GHz
Planck maps. They found a cross-correlation exceeding the expected thermal
continuum at $545$ GHz, consistent with [C II] emission correlated with BOSS
quasars. A refinement of the analysis (Yang et al., 2019) increased the
result’s significance. Still, the authors caution that greater spectral
resolution is required to verify that the excess cross-correlation is
explained by [C II] emission rather than the redshift evolution of the
correlated continuum emission.

Switzer (2017a) describes how future instruments for measuring spectral
distortions in the cosmic microwave background (CMB) can be employed for
intensity mapping. To anticipate the capabilities of future measurements, we
use data from the COBE-FIRAS instrument (Fixsen et al., 1994), in cross-
correlation with the BOSS CMASS and LOWZ galaxy catalogs, to make the first
tomographic intensity mapping constraint on [C II] emission. The FIRAS
instrument (Fixsen et al., 1994) was designed to precisely measure the
spectrum of the CMB, dust, and line emission from the Milky Way. It covers a
broad frequency range from $30$ GHz to $2910$ GHz with $13.6$ GHz spectral
channels. This, and the fact that FIRAS’ frequency range conveniently overlaps
the well-sampled LOWZ and CMASS galaxy catalogs, make FIRAS${\times}$BOSS a
natural candidate for a tomographic [C II] cross-correlation analysis.

Compared to the Planck data set used by Pullen et al. (2018), the FIRAS data
set has much higher thermal noise per pixel and much lower angular resolution
(the FIRAS beam has nearly a 7-degree full-width-half-maximum). At first
glance, it seems that these limitations give FIRAS${\times}$BOSS no hope of
approaching the sensitivity of Planck${\times}$QSO, but we achieve error bars
that are only about two times larger. One reason is that the effective number
of independent modes scales with the number of redshift bins. We use 14
redshift bins for FIRAS${\times}$CMASS and 16 for FIRAS${\times}$LOWZ, whereas
Planck${\times}$QSO only has one channel with correlated [C II] signal.
However, the number of modes also depends on the range of measurable angular
scales, and Planck${\times}$QSO more than makes up for its lack of redshift
resolution with its much larger range of observable angular modes.
Counteracting this, FIRAS${\times}$BOSS sees a larger cosmological signal due
to structure growth at lower redshifts, and, critically, the CMASS and LOWZ
galaxy catalogs are well-sampled enough that shot noise is subdominant at the
redshifts and angular scales we consider. The BOSS quasar sample used in
Planck${\times}$QSO, on the other hand, is dominated by shot noise due to the
sparser sampling.

###  1.1 Spherical harmonic tomography

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.

##  2 Parameter Estimation with the $C_{\ell}(z,z^{\prime})$ Statistic

###  2.1 Motivating the estimator

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). A recent analysis of BOSS
CMASS and LOWZ clustering using Spherical Harmonic Tomography (Loureiro et
al., 2019) 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. Another important feature of SHT 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)$
statistic insufficient. However, since $C_{\ell}(z,z^{\prime})$ does not
compress data along the line-of-sight direction, it describes redshift
evolution. 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)$
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.

###  2.2 Computing $C_{\ell}(z,z^{\prime})$ with Incomplete Sky Coverage

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

$D_{\ell}=\sum_{\ell^{\prime}}M^{AB}_{\ell\ell^{\prime}}\omega_{\ell^{\prime}}^{A}\omega_{\ell^{\prime}}^{B}C_{\ell^{\prime}},$
(1)

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

$M^{AB}_{\ell\ell^{\prime}}=\frac{2\ell^{\prime}+1}{4\pi}\sum_{\ell^{\prime\prime}}(2\ell^{\prime\prime}+1)\mathcal{W}^{AB}_{\ell^{\prime\prime}}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
0&0&0\end{pmatrix}^{2},$ (2)

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

$\begin{split}D_{b}&\equiv\sum_{\ell}B_{b\ell}D_{\ell}\\\
M^{AB}_{bb^{\prime}}&\equiv\sum_{\ell,\ell^{\prime}}B_{b\ell}M^{AB}_{\ell\ell^{\prime}}\omega^{A}_{\ell^{\prime}}\omega^{B}_{\ell^{\prime}}B^{u}_{\ell^{\prime}b^{\prime}}.\end{split}$
(3)

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

$\hat{C}_{b}=\sum_{b^{\prime}}\left(M^{AB}_{bb^{\prime}}\right)^{-1}\hat{D}_{b^{\prime}}.$
(4)

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

$\tilde{C}_{b}\equiv\sum_{b^{\prime}}\left(M^{AB}_{bb^{\prime}}\right)^{-1}\sum_{\ell}B_{b^{\prime}\ell}M^{AB}_{\ell\ell^{\prime}}\omega^{A}_{\ell^{\prime}}\omega^{B}_{\ell^{\prime}}C_{\ell^{\prime}}.$
(5)

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

$\begin{split}&\langle\Delta C_{b}^{A\times B}(z_{1},z_{2})\Delta
C_{b}^{A\times B}(z_{3},z_{4})\rangle\approx\\\ &\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}[C_{b}^{A\times A}(z_{1},z_{3})C_{b}^{B\times
B}(z_{2},z_{4})+\\\ &C_{b}^{A\times B}(z_{1},z_{4})C_{b}^{B\times
A}(z_{2},z_{3})],\end{split}$ (6)

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

$\begin{split}&\langle\Delta C_{b}^{\times}(z_{1},z_{2})\Delta
C_{b}^{\times}(z_{3},z_{4})\rangle\approx\\\ &\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm
IM}(z_{1},z_{3})C_{b}^{g}(z_{2},z_{4})].\end{split}$ (7)

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.

###  2.3 Implementing the estimator

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})$.

$\mathbf{x}\equiv vec[C_{b}(z,z^{\prime})].$ (8)

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

$\begin{split}-2\ln\mathcal{L}=&[\bf{x}(\Theta)-\hat{\bf{x}}]^{T}\bf{\Sigma}(\Theta)^{-1}[\bf{x}(\Theta)-\hat{\bf{x}}]\\\
&+\ln|\bf{\Sigma}(\Theta)|-2k\ln(2\pi),\end{split}$ (9)

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; 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.

##  3 The FIRAS and BOSS data sets

###  3.1 FIRAS Instrument and Data Set

Figure 1:  Each interferogram measured by FIRAS was multiplied by an
apodization function before Fourier transforming. The result is that the FIRAS
maps show the true sky spectrum convolved by the Fourier transform of that
apodization function. We call this Fourier transformed apodization function
the FIRAS frequency response function, $A(\Delta\nu)$, and plot it here.

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
111http://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$.

Figure 2:  Mollweide projection of the FIRAS inverse-noise weights for [C II]
emission at $z\sim 0.52$. We zero the weights at all regions that do not
overlap with the CMASS galaxy survey. Although this plot only shows the
weights for the single frequency bin centered at $z\sim 0.52$, the angular
distribution of the weights is identical for all frequencies. The frequency
dependence of the weights can be inferred from the orange curve of Figure 5.

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.

Figure 3:  The measured FIRAS beam (solid green) and our interpolation (dashed
blue). The beam was measured at a frequency of $1.5$ THz via observations of
the Moon. Figure 4:  The simulated FIRAS angular window function,
$\omega_{FIR}(\ell)$, which accounts for beam, scan, and pixelization effects.

###  3.2 BOSS Data Set

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).

Figure 5:  Galaxy counts as a function of redshift for the BOSS LOWZ (blue)
and CMASS (red) populations. Plotted in yellow is the expected FIRAS noise
angular power spectrum, $C_{\ell}$, at $\ell=0$ for the corresponding
frequencies of the [C II] fine-structure line. For most of the angular modes
accessible to FIRAS, the noise power spectrum will be higher, growing roughly
proportional to the inverse of the square of the FIRAS beam and scan window
function. For $z<0.2$, the FIRAS noise grows rapidly. The shaded blue and red
rectangles indicate the redshift ranges used for the cross-power analysis of
the FIRAS data with the LOWZ and CMASS samples. These regions were selected
because both the FIRAS thermal noise and galaxy shot noise are low.

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

$\delta^{g}(z,\theta)=\frac{n(z,\theta)}{\bar{n}(z,\theta)}-1,$ (10)

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$.

Figure 6:  Mollweide projection of the CMASS over-density map and selection
function at $z\sim 0.52$, near the redshift peak of CMASS, in Galactic
coordinates. There are several points in the over-density map well above 4,
but we choose to saturate the scale at 4 to show the broad clustering
features.

Previous analysis of the BOSS data has found evidence of systematic
contamination from stellar populations at low $\ell$ (Loureiro et al., 2019).
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).

##  4 Signal Model and Parameter Fit

The low spatial resolution of the FIRAS maps limits the range of scales to
$2<\ell<47$, or five bandpowers with $\Delta\ell=9$. Figure 7 shows the
expected variance on the cross-power signal for each $\ell$ and the
$\ell$-bins we use in this analysis. Since the first of the five available
$\ell$-bins shows signs of stellar contamination in the BOSS sample, we drop
it from our analysis. We also eliminate the second $\ell$-bin because, in that
bin, the mode mixing caused by partial sky coverage combines with the steep
angular index of Milky Way emission to mix the FIRAS auto-power negative. This
makes the empirical FIRAS model non-positive definite in that bin and
therefore unusable as verification for our cross-power fits. Consequently, we
use only the last three bins, as indicated in Figure 7, which contain most of
the sensitivity. In principle, extra information could be obtained at higher
$\ell$ by re-making the FIRAS maps from the raw data at higher $N_{\rm side}$,
but in practice, the 7-degree beam of the FIRAS instrument diminishes the
signal-to-noise ratio at higher $\ell$s (see the upward trend in noise-to-
signal ratio for $C_{\ell}^{\times}(z,z^{\prime})$ at $\ell>35$ in Figure 7).
To model the errors, we also restrict the FIRAS${\times}$FIRAS and
BOSS${\times}$BOSS analysis to these same 3 $\ell$-bins, even though much
higher $\ell$ information is available from the BOSS catalog. Similarly, we
must restrict the FIRAS${\times}$FIRAS analysis to the sky fraction covered by
the BOSS galaxy survey to avoid over-estimating the error bars by including
the bright Galactic plane, which has no overlap with the BOSS North or South
fields.

Figure 7:  The approximate expected variance in the cross-power diagonal at
$z{\sim}0.4$, according to the Gaussian error formula $\langle\Delta
C^{x}_{\ell}(z,z^{\prime})\Delta
C^{x}_{\ell}(z,z^{\prime})\rangle\approx(f_{\rm sky}(2\ell+1))^{-1}C^{\rm
IM}_{\ell}(z,z^{\prime})C_{\ell}^{g}(z,z^{\prime})$. Although the magnitude
varies somewhat with redshift, the shape is representative of all redshifts
studied. The small number of modes and large galactic foregrounds drive high
variance at low $\ell$. For $\ell>20$, the foregrounds have mostly subsided,
and thermal noise dominates. For $\ell>30$, the increase in the noise caused
by beam and scan convolution starts to overtake the advantage of extra modes
at higher $\ell$. Since the FIRAS data was originally mapped at $N_{\rm
side}=16$, spherical harmonics below $\ell=3N_{\rm side}=48$ form a complete
basis, and no further information can be extracted by considering higher
$\ell$. Dashed lines indicate the bounds of the 5 $\ell$-bins, and the shaded
regions show the three bins used in the cross-power analysis.

###  4.1 Dark Matter Model

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

$C_{\ell}^{A\times B}(z,z^{\prime})=\frac{2}{\pi}\int
k^{2}P^{\delta}(k,z{=}0)W_{A}^{\rm tot}(k,z)W_{B}^{\rm tot}(k,z^{\prime})dk,$
(11)

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),

$W_{A}(k,z)=b_{A}\int\phi_{z}(z^{\prime\prime})G(z^{\prime\prime},k)j_{\ell}[k\chi(z^{\prime\prime})]dz^{\prime\prime},$
(12)

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

$\begin{split}W_{A}^{\rm tot}(k,z)=&W_{A}(k,z)+W_{A}^{\rm RSD}(k,z)\\\
W_{A}^{\rm
RSD}(k,z)=&b_{A}\int\beta_{A}(z^{\prime\prime})\phi_{z}(z^{\prime\prime})\times\\\
\biggl{[}&\frac{2\ell^{2}+2\ell-1}{(2\ell+3)(2\ell-1)}j_{\ell}(k\chi(z^{\prime\prime}))\\\
&-\frac{\ell(\ell-1)}{(2\ell-1)(2\ell+1)}j_{\ell-2}(k\chi(z^{\prime\prime}))\\\
&-\frac{(\ell+1)(\ell+2)}{(2\ell+1)(2\ell+3)}j_{\ell+2}(k\chi(z^{\prime\prime}))\biggr{]}dz^{\prime\prime},\\\
\end{split}$ (13)

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

$C_{\ell}^{A\times
B}(z,z^{\prime})=b_{A}b_{B}C_{\ell}^{(2)}(z,z^{\prime})+\frac{b_{A}{+}b_{B}}{2}C_{\ell}^{(1)}(z,z^{\prime})+C_{\ell}^{(0)}(z,z^{\prime}).$
(14)

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.

###  4.2 FIRAS${\times}$BOSS

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

$\begin{split}&C_{\ell}^{\rm[C{\sc\,II}]\times g}(z,z^{\prime})=\\\
&I_{\rm[C{\sc\,II}]}(z)\cdot\Big{[}b_{g}b_{\rm[C{\sc\,II}]}C^{(2)}_{\ell}(z,z^{\prime})+\\\
&\frac{b_{g}+b_{\rm[C{\sc\,II}]}}{2}C^{(1)}_{\ell}(z,z^{\prime})+C^{(0)}_{\ell}(z,z^{\prime})\Big{]}.\end{split}$
(15)

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

$\begin{split}&C_{\ell}^{c\times g}(z,z^{\prime})=\\\
&\sum_{z^{\prime\prime}}\frac{dI_{\rm
CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta
z^{\prime\prime}\times\left[b_{g}b_{\rm[C{\sc\,II}]}C^{(2)}_{\ell}(z^{\prime\prime},z^{\prime})+\right.\\\
&\left.\frac{b_{g}+b_{\rm[C{\sc\,II}]}}{2}C^{(1)}_{\ell}(z^{\prime\prime},z^{\prime})+C^{(0)}_{\ell}(z^{\prime\prime},z^{\prime})\right],\end{split}$
(16)

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

$C_{\ell}^{\times}(z,z^{\prime})=A(z^{\prime\prime})\circledast\left[C_{\ell}^{\rm[C{\sc\,II}]\times
g}(z^{\prime\prime},z^{\prime})+C_{\ell}^{c\times
g}(z^{\prime\prime},z^{\prime})\right].$ (17)

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$.

Figure 8: Cross-power fits for the line and continuum amplitude from
FIRAS${\times}$BOSS, showing the CMASS analysis on the left and the LOWZ
analysis on the right. The top row shows the MCMC contours of our fit to the
data. Dark blue and light blue regions represent the 68 and 95 percent
contours, respectively. The MCMC analysis uses 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 a
minimal effect on our best-fit $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$
values, but it serves to prevent nonphysical negative values from counting
towards our upper limit constraints. The bottom row shows red histograms of
the distribution of best fit $\chi^{2}$ for 10,000 simulations in which we
draw $a_{\ell m}$ amplitudes from Gaussian distributions for the cosmological
signal, galaxy shot noise, and FIRAS auto-power (for more details of the
simulation, refer to the parametric simulations described in Appendix E). The
$\chi^{2}$ of the maximum likelihood fits to the data are plotted as vertical
dashed blue lines.

###  4.3 Modeling the FIRAS${\times}$BOSS Covariance

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.

####  4.3.1 [C II] and CIB signal

We simulate the [C II] and CIB variance models through map-space simulations
that include FIRAS instrumental effects. The [C II] and CIB signals are
painted onto maps drawn from the cosmological clustering signal with linear
bias. For the [C II] signal, this is accomplished by multiplying the drawn
maps by $I_{\rm[C{\sc\,II}]}(z)$. For the CIB signal, the maps are matrix-
multiplied by $\frac{dI_{\rm
CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta
z^{\prime\prime}$, in a map-space analogy to Equation 16. These maps are then
convolved by the FIRAS redshift response function, $A(z)$.

The magnitude of the portion of the covariance that comes from this [C II] and
CIB signal is a function of the [C II] and CIB amplitudes. In order to account
for this, our covariance is constructed from a linear combination of four
separate simulations (accounting for each cross-term). Appendix D describes
this process.

####  4.3.2 Clustering signal and shot noise from BOSS${\times}$BOSS

We model the variance in the BOSS survey as a tracer of the dark matter with
constant bias and linear RSD plus shot noise, or

$\begin{split}C_{\ell}^{g}(z,z^{\prime})&=b_{g}^{2}C_{\ell}^{(2)}(z,z^{\prime})+b_{g}C_{\ell}^{(1)}(z,z^{\prime})\\\
&+C_{\ell}^{(0)}(z,z^{\prime})+A_{SN}\frac{\Omega_{\rm
pixel}}{\bar{n}(z)}\delta(z,z^{\prime}),\end{split}$ (18)

where $\bar{n}(z)$ is the average number of BOSS galaxies per pixel in each
redshift slice, and $\Omega_{\rm pixel}$ is the angular size of a pixel in
steradians. While the complete BOSS${\times}$BOSS power spectrum requires
additional modeling to describe all scales (Loureiro et al., 2019), this
simple model with free parameters for a constant bias and shot-noise amplitude
is sufficient to describe the BOSS variance relevant to the angular scales
analyzed in FIRAS${\times}$BOSS. Due to the large number of galaxies in the
CMASS and LOWZ samples and the limited $\ell$-range of our analysis, the shot
noise is considerably smaller than the clustering signal, so $A_{SN}$ is
weakly constrained. The bias parameter is fit to a value of 1.81 and 1.82 for
the LOWZ and CMASS samples, consistent with previous work (Salazar-Albornoz et
al., 2017). We find reasonable $\chi^{2}$ values for both the CMASS
($\chi^{2}$ per degree of freedom of 1.02, PTE of 0.37) and LOWZ ($\chi^{2}$
per degree of freedom of 1.08, PTE of 0.12) fits. Figure 9 shows the measured
$C_{b}(z,z^{\prime})$ from CMASS versus the best-fit model for the three
$\ell$-bins we consider.

The covariance we use for the BOSS${\times}$BOSS parameter estimation is
computed using the approximate formula of Tristram et al. (2005) described in
Appendix C. Since the covariance is a function of the model, an MCMC analysis
must, in principle, recalculate the covariance for each different estimate of
the underlying model parameters. For BOSS${\times}$BOSS, we instead employ an
iterative approach. We make a best guess of the parameters, compute a
covariance for that guess, find a new maximum likelihood solution using that
covariance, and then repeat the procedure. After several iterations, this
procedure has converged, and the parameters assumed in the model for the
covariance equal the parameters at the maximum likelihood peak of our MCMC fit
to within a fractional tolerance of $3\times 10^{-3}$.

Figure 9:  The BOSS${\times}$BOSS binned angular power spectrum,
$C_{b}(z,z^{\prime})$, for the CMASS galaxies. The left and right columns show
the data and the best-fit model, respectively. The rows show the three angular
bins used in this analysis, each of size $\Delta\ell=9$ and centered at
$\ell=24$, 33, and 42. At the resolution of FIRAS, most of the cosmological
clustering signal occurs on the diagonal, where $z=z^{\prime}$, though there
is also a small correlation between neighboring redshift bins, visible just
off of the diagonal. The variation in amplitude along the diagonal is due to
the redshift evolution of the growth factor, the change in physical scales
being probed as a function of redshift, and changes in the CMASS galaxy
density.

####  4.3.3 Thermal noise and foregrounds from FIRAS${\times}$FIRAS

We model the FIRAS thermal noise from the inverse noise variance weights
provided by the FIRAS collaboration (FIRAS Collaboration, 1997). Let
$W(\theta,z)$ represent the FIRAS inverse noise variance maps at $N_{\rm
side}{=}128$. Thermal noise contributes constant variance in $\ell$-space,
given by $N(z){=}\Omega_{128}\frac{\langle
W(\theta,z)\rangle_{\theta}}{\langle W^{2}(\theta,z)\rangle_{\theta}}$, where
$\Omega_{128}$ is the pixel size in steradians, and the angular average is
taken only over pixels that overlap the BOSS galaxy survey. We then model the
FIRAS thermal noise as

$N(\ell,z,z^{\prime})\equiv
N(z)^{1/2}N(z^{\prime})^{1/2}A\left(|\nu(z)-\nu(z^{\prime})|\right),$ (19)

where $A\left(|\nu(z)-\nu(z^{\prime})|\right)$ accounts for the convolution of
the FIRAS spectrum by the Fourier Transform Spectrometer’s frequency response
function. The measured frequency correlations of the thermal noise agree with
the covariance model of the FIRAS collaboration (FIRAS Collaboration, 1997).
For this thermal noise component, we compute the expected binned angular power
spectrum (Equation 5) by applying the $N_{\rm side}=16$ pixel window function
to $N(\ell,z,z^{\prime})$, binning into bandpowers, and then unmixing with the
binned mixing matrix that uses the full simulated beam, scan, and pixel window
function. The result is that, after the binning and unmixing operator is
applied, the initially flat angular spectrum of the noise now rises roughly as
the inverse square of the FIRAS beam and scan window function.

We model the Milky Way foreground angular power spectrum as a simple power-law
with a free angular index $\gamma$ and amplitude $A_{\rm MW}^{2}$ that we fit
to the FIRAS auto-power spectrum through the form

$C^{\rm MW}_{\ell}(z,z^{\prime})=A_{\rm
MW}^{2}\ell^{-\gamma}D(z,T_{d})D(z^{\prime},T_{d}).$ (20)

The spectrum of the Galactic emission, $D(z,T_{d})$, is modeled as semi-
thermal dust emission, given by

$D(z,T_{d})\propto\nu^{\beta}B_{\nu}(T_{d}),$ (21)

where $\nu$ is converted to redshift assuming the [C II] line, $\beta=1.5$
(Planck Collaboration et al., 2014a), $B_{\nu}$ is the Planck function, and
the dust temperature $T_{d}$ is a free parameter. In principle, the Milky Way
emission is also convolved in the frequency direction by the frequency
response function of the FIRAS spectrometer, but the effect is negligible for
smooth spectral emission, so we do not include it. The full model we fit to
the data is

$C^{\rm IM}_{\ell}(z,z^{\prime})=A_{\rm
MW}^{2}\ell^{-\gamma}D(z,T_{d})D(z^{\prime},T_{d})+A_{N}N(\ell,z,z^{\prime}).$
(22)

There are four free parameters: the Milky Way amplitude $A_{\rm MW}$ at
$\ell\sim 1$ (units of MJy/sr), the dust temperature $T_{d}$ (units of K), the
unitless angular power-law index $\gamma$, and a unitless factor $A_{N}$
multiplying the expected noise signal ($A_{N}$ is expected to be near 1).
Figures 10 and 11 show color plots of the best-fit models and the data for
LOWZ and FIRAS respectively, over the full redshift and angular range of our
analysis. Figure 12 shows the redshift diagonal of the data and best-fit
models for both the LOWZ and CMASS redshift ranges, along with error bars
estimated from Monte Carlo simulations drawn from the best-fit model, fully
including the effects of FIRAS beam convolution, ecliptic scan convolution,
pixelization, and partial sky coverage.

The covariance we use to fit FIRAS${\times}$FIRAS to our parametric model is
simulated. We draw $a_{\ell m}$ amplitudes from a Gaussian distribution whose
variance is given by our model (Equation 22). We use this to produce 10,000
full-sky maps, to which we apply our FIRAS window function. We then use the
FIRAS inverse-noise weights and window function to compute a simulated
observed binned partial-sky $C_{b}(z,z^{\prime})$ for each of these 10,000
draws. The covariance computed from these simulated $C_{b}(z,z^{\prime})$
amplitudes accounts for the effects of beam convolution, ecliptic scan
smearing, and partial sky coverage. Since the covariance is a function of the
assumed parameters for the model, it should, in principle, be recalculated for
each different estimate of the underlying model parameters. For a simulated
covariance, this procedure is computationally expensive. As with
BOSS${\times}$BOSS, we instead employ an iterative approach for
FIRAS${\times}$FIRAS. We make a best guess of the parameters, simulate a
covariance for that guess, find a new maximum likelihood solution using that
covariance, and then repeat the procedure.

We repeat the above iteration until the $\chi^{2}$ per d.o.f. converges to
${<}1\%$. Figures 10, 11, and 12 show the parametric model compared to the
measured FIRAS${\times}$FIRAS power spectrum. The fit converges to a $1\%$
constraint on the thermal noise amplitude centered on 0.97 for LOWZ and 0.90
for CMASS. The three foreground parameters are correlated and weakly
constrained but yield dust temperatures consistent with the $15-25$ K measured
over the region of the sky we observe (Liu et al., 2017). Milky Way emission
contributes approximately half of the variance to the $\ell=24$ bandpower and
is negligible at the other two bandpowers. Because there are relatively few
spatial modes in the BOSS regions at $\ell=24$, and because the BOSS region is
relatively clear of Milky Way emission, the constraint on the Milky Way
contribution to $\ell=24$ is uncertain to $30\%$, but this uncertainty has a
small impact on the final bounds on line brightness from FIRAS${\times}$BOSS.
Since the foregrounds account for half of the redshift-diagonal variance at
$\ell=24$, a $30\%$ increase would cause a $15\%$ increase in total $\ell=24$
bin variance. The $\ell=24$ bin makes up roughly half of the total [C II]
signal-to-noise ratio, so this could at most increase the [C II] variance by
$7.5\%$. In amplitude, a $3.75\%$ increase in the final CII constraint is 0.01
MJy/sr. In reality, this is conservative because the foregrounds have long
frequency correlations that the [C II] signal does not have, so most of the
impact would instead be on the CIB constraint.

The best-fit parametric models have low $\chi^{2}$ per d.o.f. (0.74 for LOWZ
and 0.86 for CMASS), indicating that the error model may be overestimated. An
alternative to the parametric model approach is instead to use the measured
auto-power spectrum of the FIRAS data, $\hat{C}_{b}^{\rm IM}$, as the model of
variance in the FIRAS data cube. The attraction of this approach is that it
captures features that could be missing in any parametric model fit.
Unfortunately, because the FIRAS auto-power signal also contains the [C II]
and CIB signal, this approach introduces cosmic bias wherein high-scattering
modes with more [C II] and CIB signal are artificially down-weighted in the
cross-power analysis. This could bias both the measured [C II] signal and its
error bars low, an effect we have measured in simulations. Although the
measured FIRAS auto-power cannot be used in a covariance model acting on the
actual cross-power data, we use it on simulated cross-power signal to verify
that our parametric model produces similar error bars on the
FIRAS${\times}$BOSS parameters. Appendix E shows the results of two sets of
simulations, one with a covariance that uses the measured FIRAS auto-power and
one with a covariance that uses our best-fit parametric model. The two models
yield nearly identical error contours on the cross-power parameters.

We note that the complete FIRAS covariance has additional terms with complex
structure, described in FIRAS Collaboration (1997), which includes
correlations jointly between frequencies $\nu,\nu^{\prime}$ and position
vectors $\hat{n},\hat{n}^{\prime}$ on the sky. In the frequency range studied
here and the survey region outside the Galactic plane, this covariance is
dominated by thermal detector noise, which is diagonal in
$\hat{n},\hat{n}^{\prime}$, and by the correlated structure across frequencies
produced by Milky Way emission, at low multipoles. In measurements with future
instruments that achieve high significance, the absolute calibration error
must also be included. The FIRAS auto-power also contains the [C II] signal
and the continuum CIB, but at a low level that negligibly effects our
parametric thermal noise and foregrounds fits. In addition, there are several
prominent Galactic spectral lines (Fixsen et al., 1999). The only line that
falls in the frequencies of our analysis is the 205.3 $\mu{\rm m}$ NII line.
Although this line is visible in the full-sky FIRAS auto-power spectrum, it is
not detectable when we restrict the data to the BOSS survey regions, which are
out of the Galactic plane, so we do not include it in our model.

Figure 10:  The FIRAS${\times}$FIRAS binned angular power spectrum,
$C_{b}(z,z^{\prime})$, in the LOWZ redshift range for three angular bins of
width $\Delta\ell=9$, centered at $\ell=24$, 33, and 42. The left column shows
the data, and the right column shows the best-fit model, of the form of
Equation 22. The structure on the diagonal of the plots is due to thermal
noise, with small just-off-diagonal correlations due to the FIRAS frequency
window function, $A(\nu)$. The thermal noise increases with higher $\ell$,
roughly as the inverse of the FIRAS beam and scan window function squared. The
foregrounds are visible in the $\ell=24$ bin as a roughly constant offset to
all $z,z^{\prime}$ combinations. The foreground amplitude is comparable to the
thermal noise at $\ell=24$, but it drops at higher $\ell$ with a power-law
index of $\gamma\approx-2.3$. The foregrounds are negligible compared to the
thermal noise at $\ell=33$ and $\ell=42$.  Figure 11:  The
FIRAS${\times}$FIRAS binned angular power spectrum, $C_{b}(z,z^{\prime})$, in
the CMASS redshift range with layout and general properties similar to the
LOWZ region described in Figure 10. The vertical and horizontal lines in the
$\ell=24$ data at $z\sim 0.67$ are due to a spurious correlation between
thermal noise and Galactic foregrounds.  Figure 12:  The redshift diagonal of
the FIRAS${\times}$FIRAS binned angular power spectrum, $C_{b}(z=z^{\prime})$,
for both the CMASS and LOWZ redshift regions. The three angular bins, centered
at $\ell=24$, $\ell=33$, and $\ell=42$ are plotted in red, green, and blue,
respectively. The best-fit models are plotted as solid lines, and the data are
plotted as triangles or circles, with error bars computed from the best-fit
model. The $\ell=24$ data are artificially shifted backward by half a redshift
bin for visual clarity. Dashed lines show the thermal noise portion of the
best-fit model only. From the dashed lines, the effect of the FIRAS beam and
scan convolution can be seen, as the thermal noise increases from low to high
$\ell$. Only the $\ell=24$ bin has a significant foreground contribution.

##  5 Discussion

Figure 13 compares FIRAS${\times}$BOSS upper limits on
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ to several representative physical
models of [C II] brightness as a function of redshift. It also shows the
Planck${\times}$QSO intensity mapping constraint from Yang et al. (2019),
assuming that the excess power detected is [C II] emission. Models that do not
scale [C II] emission with star formation rate evolution are disfavored
because their low-redshift [C II] emission is too bright. Models that scale [C
II] emission with star formation rate can be calibrated with physical
parameters to be consistent with both FIRAS${\times}$BOSS and
Planck${\times}$QSO.

The yellow dotted region in Figure 13 shows the range of [C II] amplitudes
predicted by the collisional excitation model from Gong et al. (2012),
hereafter G12. In this model, the mean [C II] intensity is computed through a
simple radiative transfer model whose free parameters are the number density
$n_{e}$ and kinetic temperature $T_{e}^{K}$ of electrons within the emitting
galaxies. The yellow dotted region spans the range of $T_{e}^{K}$ and $n_{e}$
values considered by G12. Because this model predicts comparatively bright [C
II] emission, Yang et al. (2019) argue that it provides the best-fit to their
measurement and use their result to place constraints on the two free
parameters. FIRAS${\times}$BOSS upper limits at $z{\sim}0.5$ rule out the
brightest range of the G12 predictions. Since it is this bright end that is
consistent with the Yang et al. (2019) measurement at $z\sim 2.6$, the G12
model is disfavored to explain the FIRAS${\times}$BOSS and Planck${\times}$QSO
[C II] measurements.

The G12 model was originally created to forecast emission at much higher
redshifts, during the epoch of reionization. Since [C II] luminosity is
expected to be correlated with star formation, we qualitatively expect the [C
II] amplitude to follow the cosmic star formation history and come to a peak
around $z\sim 2-3$, declining to the present day (Madau & Dickinson, 2014).
This is not seen in the G12 model, which monotonically increases as we move to
lower redshift. Thus the bright G12 predictions, which are disfavored by the
combination of FIRAS${\times}$BOSS and Planck${\times}$QSO, may not be a
physically accurate estimate of the true [C II] evolution during more recent
epochs since $z{\sim}2$.

Next, we consider the scaling models from Silva et al. (2015) (hereafter S15),
shown as the hatched purple region in Figure 13. Silva et al. (2015) consider
four different empirically-calibrated power-law scaling relations between [C
II] luminosity and star formation rate (SFR). The hatched purple region shows
the full range of [C II] luminosity predictions from the four values of the
slope and amplitude of this power-law scaling from their Table 1. Although
these models were originally constructed to predict [C II] intensity at
reionization redshifts, they are correlated with the cosmic star formation
history and thus show the qualitative redshift evolution we expect, with a
peak at $z\sim 2$ and a decline at lower redshifts. These predictions fall a
factor of 100 or more below our FIRAS${\times}$BOSS upper limits at $z\sim
0.5$, putting them well below our ability to constrain. However, these models
also fall well below the Planck${\times}$QSO estimate, meaning they may also
be pessimistic predictions.

Finally, we examine the “semi-empirical" model from Sun et al. (2019)
(hereafter S19), plotted as a dashed green curve. This calculation falls
somewhat between the other two. It uses a physically motivated scaling between
the [C II] and infrared luminosity of a halo. This scaling is parameterized in
terms of the photoelectric heating efficiency from dust grains,
$\epsilon_{\rm{PE}}$. The distribution of galaxy infrared luminosities is
calibrated to empirical measurements of the cosmic infrared background (Planck
Collaboration et al., 2014b). This model also predicts a fiducial [C II]
amplitude well below our limits and the Yang et al. (2019) measurement.
However, as shown in Figure 10 of S19, the [C II] intensity can be brought
into agreement with the Planck${\times}$QSO measurement, if one increases
$\epsilon_{\rm{PE}}$ by a factor of ${\sim}6$ from their fiducial number. The
authors, however, note that this higher value may lead to tension with the
observed relation between [C II] luminosity and star formation rate in low-
redshift galaxies. The redshift evolution of this rescaled S19 model is
plotted as the dashed blue curve in Figure 13. In contrast with the G12 model,
we find that our FIRAS measurements are fully consistent with the redshift
evolution seen here, even scaling to the brighter emission at $z{\sim}2.6$
suggested by Yang et al. (2019).

Figure 13 also includes a point at $z{\sim}0$ from luminosity function
measurements (Hemmati et al., 2017), which integrates to $4.1\pm 2.7$ kJy/sr.
The measured galaxies cover the $10^{6.5}$ to $10^{9.5}$ solar luminosity
range, via a mix of direct detection of [C II] emission, at the bright end,
and inference of [C II] emission via FIR brightness and color measurements at
the dim end. An integral of this luminosity function yields the expected [C
II] brightness, $I_{\rm[C{\sc\,II}]}$, but in order to plot
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$, a model for the [C II] bias must be
assumed. We estimate a bias of 1, which corresponds to the value from the halo
model used in appendix B. This roughly matches the bias evolution of all the
models plotted, which predict decreasing bias at low redshift. This
convergence of bias toward unity at low redshift occurs because, as halo
masses increase at low redshift, the mass range for [C II]-hosting galaxies
becomes less biased towards the higher mass end of the halo distribution.

Even the models that scale with star formation rate struggle to decrease
quickly enough at low redshift to be consistent with both Planck${\times}$QSO
and Hemmati et al. (2017). There are several possible explanations. The
estimated $z{\sim}0$ [C II] bias of 1 could be too low. Or there could be an
undetected population of low-luminosity [C II] galaxies. In fact, Hemmati et
al. (2017) suggest the possibility of unmeasured faint low-metallicity dwarf-
galaxies, which would have a different [C II]-FIR calibration, contributing to
the faint end of the [C II] luminosity function. This possibility highlights a
trend where early intensity mapping results suggest unexpectedly bright
cumulative emission compared to direct detection of individual objects. For
instance, Breysse et al. (2021) suggests that potential bright CO emission
detected by mmIME (Keating et al., 2020) may be due to a large population of
dim CO galaxies. Another possibility is that the assumption of these models,
that [C II] emission is directly tied to the star formation rate, is too
simple.

Conclusions here are limited by thermal noise in FIRAS. To meaningfully
constrain the space of models, qualitatively new improvements in sensitivity
are required. To put FIRAS in the context of what could be achieved by a
modern mission, we compute the potential sensitivity achievable by cross-
correlating maps from the proposed PIXIE satellite with the CMASS galaxy maps
used in this paper. The Primordial Inflation Explorer (PIXIE) is a proposed
NASA Explorer-class mission to measure the polarized imprint of primordial
inflation on the CMB (Kogut et al., 2011), along with several other spectral
distortions of the CMB (Chluba et al., 2021). PIXIE will conduct an all-sky
survey with angular and spectral resolution and spectral coverage similar to
FIRAS. However, the beam will be slightly smaller, and the noise level will be
three orders of magnitude lower than FIRAS. These properties make PIXIE well-
placed to study star-formation and scale-dependent bias as a test of non-
Gaussianity via intensity mapping (Dizgah et al., 2019; Switzer, 2017b).

To assess PIXIE’s constraining power, we conduct a suite of 10,000 simulations
with no [C II] signal. These simulations project that a cross-correlation of
PIXIE maps with CMASS galaxies could place an upper limit of $100$ Jy/sr on
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ at $95\%$ confidence (see the red
transparent region in Figure 13). This result means that PIXIE${\times}$CMASS
is easily sensitive enough to detect even the most pessimistic [C II] models
at high significance.

Figure 13: Comparison of measured values of
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ to various models from the
literature. The gray boxes show the values allowed by the FIRAS${\times}$LOWZ
and FIRAS${\times}$CMASS 95% upper limits, with their widths representing the
redshift ranges used. The black circle shows the Planck${\times}$QSO
measurement from Yang et al. (2019). The horizontal error bars show the full
redshift range of that measurement, and the vertical error bars show the 95%
confidence intervals on their measurement of
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$. The black square shows the value
inferred from the $z\sim 0$ [C II] luminosity function (Hemmati et al., 2017),
assuming $b_{\rm[C{\sc\,II}]}=1$. The error bars come from their $1\sigma$
errors on luminosity function parameters. The purple hatched region shows the
range of brightnesses allowed by the scaling-relation models from Silva et al.
(2015), the yellow dotted region shows the range of brightnesses of the
collisional excitation models from Gong et al. (2012), and the green dashed
curve shows the semi-empirical model from Sun et al. (2019). The blue dashed
curve shows the Sun et al. (2019) model linearly rescaled to match the
Planck${\times}$QSO measurement at the appropriate redshift. The red region
below 100 Jy/sr in the FIRAS${\times}$CMASS redshift region shows the
projected sensitivity achievable with a cross-correlation between PIXIE and
CMASS galaxies in the same redshift range used for the FIRAS${\times}$CMASS
analysis. This would easily detect even the most pessimistic of the plotted [C
II] models at high sensitivity.

##  6 Conclusion

Our analysis constrains the product of the bias and specific intensity of [C
II], $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$, to be ${<}0.31$ MJy/sr at
$z{\sim}0.35$ and ${<}0.28$ MJy/sr at $z{\sim}0.57$ at $95\%$ confidence.
Through both FIRAS’ unique capability for tomographic measurement and the
depth of BOSS data at these redshifts, bounds on [C II] developed here are
competitive with those achieved in the recent analysis of Planck data (Yang et
al., 2019; Pullen et al., 2018). We rule out a swath of collisional excitation
models consistent with Planck${\times}$QSO, but note that those models were
not designed to predict low redshift [C II]. [C II] emission is expected to
correlate with the star formation rate, which appears to have peaked at $z\sim
2-3$ and then declined at later redshifts (Madau & Dickinson, 2014). Our
constraints are consistent with all models that include this expected
correlation of [C II] emission with star formation rate. If the high [C II]
intensity measured from Planck${\times}$QSO is accurate, then [C II] models
with realistic redshift evolution tied to the star formation rate predict
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}\approx 0.1$ MJy/sr at $z\sim 0.5$.
This expected value is nearly within reach of the FIRAS data. Measurements
described here are limited by FIRAS thermal noise and angular resolution
rather than the depth of the galaxy redshift survey, so substantial further
progress will require new intensity mapping measurements. We show that the
proposed PIXIE mission could detect even pessimistic models for [C II]
emission. The SHT approach developed here also applies well to wide-field
surveys such as HIRAX, CHIME, TianLai, HERA, and SPHEREx, both in auto-power
and crosspower. The SHT is broadly well-matched to intensity mapping analysis
through its ability to capture all Gaussian information directly in observed
map space and to easily account for survey effects such as the curved sky and
inhomogeneous noise.

## Acknowledgements

We would like to thank Dale Fixsen, Nathan Miller, and Nils Odegard for their
useful input regarding the FIRAS noise and beam model. Some of the results in
this paper have been derived using the healpy and HEALPix package. PCB was
supported by the James Arthur Postdoctoral Fellowship.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High Performance Computing at the University of Utah. The SDSS
website is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

## Data Availability

The FIRAS data underlying these results are available on NASA’s LAMBDA archive
at https://lambda.gsfc.nasa.gov/product/cobe/firas_prod_table.cfm. The BOSS
data underlying these results are on the BOSS data release 12 archive at
https://www.sdss.org/dr12/.

##  Appendix A Cross-power Approximate Sensitivity Forecast

The complete tomographic likelihood for $C_{b}^{\times}(z,z^{\prime})$ has a
complex structure. This appendix simplifies the sensitivity to its essential
elements to provide some intuition. For a survey where the intensity map auto-
power is noise and/or foreground dominated, the covariance on the cross-power
along the angular diagonal ($b=b^{\prime}$) is roughly

$\langle\Delta C_{b}^{\times}(z_{1},z_{2})\Delta
C_{b}^{\times}(z_{3},z_{4})\rangle\approx\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm IM}(z_{1},z_{3})C_{b}^{g}(z_{2},z_{4})].$
(23)

FIRAS lacks the redshift resolution for there to be a significant signal in
the [C II] or galaxy angular cross-power when $z\neq z^{\prime}$. Considering
only cross-correlations between the same redshift slice, the covariance
$\langle\Delta C_{b}^{\times}(z_{1},z_{2})\Delta
C_{b}^{\times}(z_{3},z_{4})\rangle$ is zero except along the redshift diagonal
$z_{1}{=}z_{2}{=}z_{3}{=}z_{4}$. We also assume that the bandpowers are broad
enough to have negligible correlations for $b\neq b^{\prime}$. Dropping
redshift subscripts, the binned [C II] cross-power signal can be modeled as
$C^{\times}_{b}=b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}C^{\delta}_{b}$, where
$C^{\delta}_{b}$ is the binned dark matter angular power spectrum,
$b_{\rm[C{\sc\,II}]}$ is the bias of [C II] galaxies, and
$I_{\rm[C{\sc\,II}]}$ is the specific intensity of [C II] emission. From this
signal model and Equation 23, the variance on a determination of
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ from the angular cross-power of any
single redshift and an angular bin is 222An interesting feature of this
formula is that, for a wide range of redshift bin sizes, this variance on the
line amplitude is roughly the same. For instance, if FIRAS had half of its
actual redshift resolution, the cosmological and galaxy clustering power
spectrum would wash out by a factor of 1/2, but the increased bandwidth would
also decrease the thermal noise signal in the FIRAS auto-power spectrum by
1/2. These effects would cancel out to keep the variance per angular bin and
redshift slice constant.

$\sigma^{2}_{bI}|_{\rm per-z,b}\approx\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}\frac{C_{b}^{\rm
IM}C_{b}^{g}}{C^{\delta}_{b}C^{\delta}_{b}}.$ (24)

The total constraining power is an inverse variance weighted sum over the
angular and redshift bins,

$\sigma^{2}_{bI}\approx\left[\sum_{b,z}\left(\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}\frac{C_{b}^{\rm
IM}C_{b}^{g}}{C^{\delta}_{b}C^{\delta}_{b}}\right)^{-1}\right]^{-1}.$ (25)

From this formula, we predict a sensitivity of 0.27 MJy/sr ($2\sigma$) on
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$. In computing this estimate, we have
neglected foregrounds and assumed that only thermal noise contributes to
$C_{b}^{\rm IM}$, since smooth foregrounds are de-weighted with only a small
penalty to sensitivity. This simple estimate closely matches the precision we
achieve with the complete calculation.

##  Appendix B Details of model for FIRAS${\times}$BOSS

To model the dust and [C II] emission, we employ a technique similar to Pullen
et al. (2018). Using the halo model, the redshift distribution of specific
intensity is

$\begin{split}\frac{dI_{\nu}}{dz}=\frac{d\chi}{dz}\frac{dI_{\nu}}{d\chi}&=\frac{d\chi}{dz}\int
dM\frac{L_{\nu}(M,z)}{4\pi}\frac{dn(z)}{dM}\\\
&=\frac{c}{(1+z)H(z)}\bar{j_{\nu}}(z),\end{split}$ (26)

where $L_{\nu}(M,z)$ represents the redshifted specific luminosity for a halo
of mass $M$, given by $L_{\nu}(M,z)=\frac{L_{\nu(1+z)}(M,z)}{1+z}$, where
$L_{\nu(1+z)}(M,z)$ represents the rest-frame specific luminosity evaluated at
$\nu(1+z)$. On the second line of Equation 26, we introduced the mean comoving
emission coefficient $\bar{j_{\nu}}(z)$, defined as

$\bar{j_{\nu}}(z)=\int dM\frac{L_{\nu(1+z)}(M,z)}{4\pi}\frac{dn(z)}{dM}.$ (27)

Following the analysis of extragalactic cosmic infrared background (CIB)
emission of Shang et al. (2012), the specific luminosity is parameterized as

$L_{\nu(1+z)}(M,z)=L_{0}\Phi(z)\Sigma(M)\Theta\left[\nu(1+z)\right],$ (28)

where $L_{0}$ is a normalization constant, $\Sigma(M)$ parameterizes the mass
dependence of the halo luminosity, and $\Phi(z)$ parameterizes the redshift
evolution of both the continuum and [C II] luminosity, which we assume evolve
together. $\Phi(z)$ is expected to have peaked around $z\sim 2-3$ (Madau &
Dickinson, 2014). A sufficiently precise intensity mapping experiment with a
large frequency bandwidth can, in principle, constrain the form of the
redshift dependence of $\Phi(z)$. FIRAS${\times}$BOSS, however, lacks the
precision to measure this evolution, so we choose to model it with a simple
power-law, $\Phi(z)=(1+z)^{2.3}$, following the parameterization and fit of
Pullen et al. (2018). $\Sigma(M)$ is assumed to follow a log-normal
distribution

$\Sigma(M)\propto\frac{M}{\sqrt{2\pi\sigma^{2}_{L/M}}}\exp{-\left[\frac{(\log_{10}M-\log_{10}M_{\rm
eff})^{2}}{2\sigma^{2}_{L/M}}\right]},$ (29)

where $\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). The integral over the
halo mass function $dn(z)/dM$ is computed with the python package hmf (Murray,
2014), assuming Planck 2015 (Ade et al., 2016) cosmological parameters and a
Tinker et al. (2008) functional form for $dn(z)/dM$. We include [C II]
emission in this specific luminosity by assuming that each CIB galaxy emits
continuum emission and a narrow [C II] line. Assuming that the ratio of total
power emitted by the line and the continuum is constant, and ignoring the
finite mass-dependent velocity width of the galaxy, we parameterize
$\Theta(\nu)$ as the sum of continuum emission ($\Theta^{c}(\nu)$, parametric
form from Shang et al. (2012)) and [C II] line emission

$\Theta(\nu)=\Theta^{c}(\nu)+\Theta^{\rm[C{\sc\,II}]}(\nu),$ (30)
$\Theta^{c}(\nu)=\left\\{\begin{array}[]{@{}ll@{}}(\nu/\nu_{0})^{\beta}[B_{\nu}(T_{d})/B_{\nu_{0}}(T_{d})]&\text{~{}~{}~{}~{}if}\
\nu<\nu_{0}\\\
(\nu/\nu_{0})^{-2}&\text{~{}~{}~{}~{}otherwise},\end{array}\right.$
$\Theta^{\rm[C{\sc\,II}]}(\nu)=\left[\int\Theta^{c}(\nu^{\prime})d\nu^{\prime}\right]\alpha\delta(\nu-\nu_{\rm[C{\sc\,II}]}),$

where we have defined $\alpha$ to be the ratio of the total power of the
continuum emission to the total power of [C II] line emission. We choose
$\beta=1.5$, in agreement with Planck Collaboration et al. (2014a). Serra et
al. (2014b) measured the dust temperature of the BOSS CIB to be 26 K, which
means all of the CIB emission in the FIRAS frequencies of this analysis is
below the turnover frequency $\nu_{0}$, so only the semi-thermal form of
$\Theta^{c}(\nu)$ is used.

The FIRAS instrument measures the average intensity in a given frequency bin.
Ignoring redshift space distortions and any contributions from thermal noise
and Milky Way foregrounds, the observed intensity is given by

$I_{\nu}(\theta,\phi)=\frac{1}{\Delta\nu}\int_{\nu-\Delta\nu/2}^{\nu+\Delta\nu/2}d\nu^{\prime}\int\left[\frac{dI^{c}_{\nu^{\prime}}}{dz}+\frac{dI^{\rm[C{\sc\,II}]}_{\nu^{\prime}}}{dz}\right]\cdot\left[1+b_{\rm
CIB}\delta\left(\chi\left(z\right),\theta,\phi\right)\right]dz,$ (31)

where we have explicitly separated the continuum and [C II] contributions. We
do this because the cross-correlation signal between FIRAS and BOSS has two
components: (1) a galaxy-[C II] cross-correlation component with narrow
redshift window functions for both the optical galaxies and the [C II]
sources, and (2) a galaxy-continuum component with a narrow redshift window
function for the optical galaxies but a broad redshift window function, given
by $\frac{dI^{c}_{\nu}}{dz}$, for the CIB continuum:

$C_{\ell}^{\times}(z,z^{\prime})=C_{\ell}^{\rm[C{\sc\,II}]\times
g}(z,z^{\prime})+C_{\ell}^{c\times g}(z,z^{\prime}).$ (32)

The angular power spectrum between signal types $i$ and $j$ is given by

$C_{\ell}^{i\times j}(z,z^{\prime})=\frac{2}{\pi}\int
W^{i}_{z}(k)W^{j}_{z^{\prime}}(k)k^{2}P(k,z=0)dk,$ (33)

where (ignoring RSD for now, and using equations 26 through 31 for the [C II]
and CIB signals) the Window functions for optical galaxies, [C II], and
continuum emission are as follows:

$W^{g}_{z}(k)=\frac{1}{\Delta z}\int_{z_{\rm min}}^{z_{\rm
max}}b_{g}G(z^{\prime},k)j_{\ell}(k\chi(z^{\prime}))dz^{\prime},$ (34)
$W^{\rm[C{\sc\,II}]}_{z}(k)=\left[\int\Theta^{c}(\nu)d\nu\right]\frac{\alpha}{\Delta\nu}\int_{z_{\rm
min}}^{z_{\rm
max}}b_{\rm[C{\sc\,II}]}G(z^{\prime},k)j_{\ell}(k\chi(z^{\prime}))\frac{cL_{0}\Phi(z^{\prime})}{4\pi(1+z^{\prime})^{2}H(z^{\prime})}\int\Sigma(M)\frac{dn(z^{\prime})}{dM}dMdz^{\prime},$
$W^{c}_{z}(k)=\int_{\nu-\Delta\nu/2}^{\nu+\Delta\nu/2}\frac{d{\nu^{\prime}}}{\Delta\nu}\int_{0}^{\infty}\frac{cL_{0}\Phi(z^{\prime})b_{\rm[C{\sc\,II}]}G(z^{\prime},k)j_{\ell}(k\chi(z^{\prime}))}{4\pi(1+z^{\prime})H(z^{\prime})}\Theta^{c}(\nu^{\prime}(1+z^{\prime}))\int\Sigma(M)\frac{dn(z^{\prime})}{dM}dMdz^{\prime},$

where $G(z,k)$ is the growth factor and $z_{\rm min}$ and $z_{\rm max}$ are
set by the spectral channels of FIRAS for $W^{\rm[C{\sc\,II}]}_{z}(k)$. As
stated in Section 3.2, the galaxy window function $W^{g}_{z}(k)$ was chosen to
have the same $z_{\rm min}$ and $z_{\rm max}$ as the FIRAS spectral channels.

If we make the reasonable approximation that all $z$-dependent quantities of
the dust and [C II] model vary slowly within the redshift bins determined by
the FIRAS frequency channels, then we can compute
$C_{\ell}^{\rm[C{\sc\,II}]\times g}(z,z^{\prime})$ and $C_{\ell}^{c\times
g}(z,z^{\prime})$ fully, including linear redshift space distortions, in terms
of $C^{(2)}_{\ell}(z,z^{\prime})$, $C^{(1)}_{\ell}(z,z^{\prime})$, and
$C^{(0)}_{\ell}(z,z^{\prime})$, which were defined in Section 4.1. For
$C_{\ell}^{\rm[C{\sc\,II}]\times g}(z,z^{\prime})$, we have

$C_{\ell}^{\rm[C{\sc\,II}]\times
g}(z,z^{\prime})=I_{\rm[C{\sc\,II}]}(z)\cdot\left[b_{g}b_{\rm[C{\sc\,II}]}C^{(2)}_{\ell}(z,z^{\prime})+\frac{b_{g}+b_{\rm[C{\sc\,II}]}}{2}C^{(1)}_{\ell}(z,z^{\prime})+C^{(0)}_{\ell}(z,z^{\prime})\right],$
(35)

where

$I_{\rm[C{\sc\,II}]}(z)=\left[\int\Theta^{c}(\nu^{\prime})d\nu^{\prime}\right]\frac{\alpha}{\nu_{\rm[C{\sc\,II}]}}\frac{cL_{0}\Phi(z)}{4\pi
H(z)}\int\Sigma(M)\frac{dn(z)}{dM}dM.$ (36)

In Equation 36, we have used $\Delta\nu=\frac{\nu_{\rm[C{\sc\,II}]}\Delta
z}{(1+z)^{2}}$. The [C II] intensity $I_{\rm[C{\sc\,II}]}(z)$ has a redshift
dependence, determined by both cosmological factors
$\left(\frac{dn(z)/dM}{H(z)}\right)$ and the assumed specific luminosity
parameterization $\left(L_{0}\Phi(z)\Sigma(M)\right)$ of Equation 28. They
lead to a [C II] intensity that increases by a factor of about 20 percent over
each of the CMASS and LOWZ redshift ranges. In our analysis, the shape of this
redshift evolution is fixed, and we fit only for the overall amplitude of
$b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}(z=z_{\rm center})$. We calculate the
CIB signal in a similar manner, as

$C_{\ell}^{c\times g}(z,z^{\prime})=\sum_{z^{\prime\prime}}\frac{dI_{\rm
CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta
z^{\prime\prime}\cdot\left[b_{g}b_{\rm[C{\sc\,II}]}C^{(2)}_{\ell}(z^{\prime\prime},z^{\prime})+\frac{b_{g}+b_{\rm[C{\sc\,II}]}}{2}C^{(1)}_{\ell}(z^{\prime\prime},z^{\prime})+C^{(0)}_{\ell}(z^{\prime\prime},z^{\prime})\right],$
(37)

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)$. Assuming the CIB spectrum does not change
significantly within the FIRAS frequency bin, then $\frac{dI_{\rm
CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta
z^{\prime\prime}$ is

$\frac{dI_{\rm
CIB}(\nu_{\rm[C{\sc\,II}]}^{z},z^{\prime\prime})}{dz^{\prime\prime}}\Delta
z^{\prime\prime}=\Delta
z^{\prime\prime}\cdot\left[\Theta^{c}\left(\nu_{\rm[C{\sc\,II}]}\frac{1+z^{\prime\prime}}{1+z}\right)\frac{cL_{0}\Phi(z^{\prime\prime})}{4\pi(1+z^{\prime\prime})H(z^{\prime\prime})}\int\Sigma(M)\frac{dn(z^{\prime\prime})}{dM}dM\right].$
(38)

The sum over $z^{\prime\prime}$ in Equation 37 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. Analogously to our [C II] analysis, 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$,
$T_{d}=26$ K, $\Phi(z)=(1+z)^{2.3}$, $\log_{10}(M_{\rm eff}/M_{\odot})=12.6$ ,
and $\sigma^{2}_{L/M}=0.5$. Because of this assumed spectral and redshift
evolution, the $\frac{dI_{\rm CIB}}{dz}$ parameter that we fit for is only the
exact continuum intensity per unit redshift at $z_{\rm center}$, the central
redshift of the survey, and $\nu_{\rm center}$, the central frequency of the
survey.

Equations 35 and 37 show how high spectral resolution allows a separation of
correlated line and continuum emission in intensity mapping. In the limit
where the spectral resolution is very low and $\Delta z^{\prime\prime}$ is
high, all the off-diagonal terms of $C_{\ell}(z^{\prime\prime},z^{\prime})$
are zero, collapsing the sum of Equation 37 to only the
$z^{\prime\prime}=z^{\prime}$ term. Dividing the [C II] signal by the CIB
signal, one finds that the ratio of [C II] to CIB cross-power scales as
$\frac{1}{\Delta z^{\prime}}$ (this point is subtle because one is free to
choose $\Delta z^{\prime\prime}$ to be an arbitrarily small value in Equation
37, but the redshift window for Equation 35 is fixed by the spectral
resolution of the intensity mapping survey. Therefore, one can choose the same
redshift window for the CIB calculation so that the $C_{\ell}(z,z^{\prime})$
terms cancel.). This simple $\frac{1}{\Delta z^{\prime}}$ scaling, where
higher redshift resolution continues to linearly improve the ability to
distinguish line emission from continuum emission, remains true until the
redshift resolution is high enough to resolve off-diagonal components of
$C_{\ell}(z^{\prime\prime},z^{\prime})$.

##  Appendix C Covariance angular coupling approximation

This section reviews a formula for approximating the $b-b^{\prime}$ coupling
in the $C_{b}$ covariance. The formula comes directly from the work of
Tristram et al. (2005), originally developed for pseudo-$C_{\ell}$ cross-
correlation studies of CMB maps from different instruments or different
detectors of the same instrument. We use the notation of Section 2, where
unbinned $C_{\ell}$ quantities retain the $\ell$ subscript, and binned
quantities are written with the subscript $b$. Let $C_{b}^{AB}$ represent the
binned angular power spectrum between maps of type $A$ and type $B$, where $A$
and $B$ can represent different redshifts or different types of maps
(intensity maps or binned galaxy surveys). Then, the covariance can be written

$\langle\Delta
C_{b}^{AB}C_{b^{\prime}}^{CD}\rangle=\sum_{b_{1},b_{2}}(M^{AB}_{bb_{1}})^{-1}\langle\Delta
D_{b_{1}}^{AB}\Delta D_{b_{2}}^{CD}\rangle(M^{CD}_{b^{\prime}b_{2}})^{-1},$
(39)

where $(M^{AB}_{bb_{1}})^{-1}$ represents the inverse of the binned mixing
matrix between maps $A$ and $B$, defined by equations 2 and 3, and where
$\langle\Delta D_{b_{1}}^{AB}\Delta D_{b_{2}}^{CD}\rangle$ is the covariance
of the pseudo-$C_{\ell}$, given by

$\langle\Delta D_{b_{1}}^{AB}\Delta
D_{b_{2}}^{CD}\rangle=\sum_{\ell_{1},\ell_{2}}B_{b_{1},\ell_{1}}B_{b_{2},\ell_{2}}\frac{\omega^{A}_{\ell_{1}}\omega^{B}_{\ell_{1}}\omega^{C}_{\ell_{2}}\omega^{D}_{\ell_{2}}}{2\ell_{2}+1}[M^{(2)}_{\ell_{1},\ell_{2}}(W^{AC,BD})C_{\ell_{1}}^{AC}C_{\ell_{2}}^{BD}+M^{(2)}_{\ell_{1},\ell_{2}}(W^{AD,BC})C_{\ell_{1}}^{AD}C_{\ell_{2}}^{BC}].$
(40)

In the above formula, $B_{b_{1},\ell_{1}}$ is a binning matrix,
$\omega^{A}_{\ell_{1}}$ is the combined beam and pixel window function for map
$A$ (and analogously for maps $B$, $C$ and $D$), and
$M^{(2)}_{\ell_{1},\ell_{2}}(W^{AC,BD})$ is the cross-correlation matrix,
defined by

$M^{(2)}_{\ell_{1},\ell_{2}}(W^{AC,BD})=\frac{2\ell_{2}+1}{4\pi}\sum_{\ell_{3}}(2\ell_{3}+1)W^{AC,BD}_{\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}^{2},$ (41)

where

$W^{AC,BD}_{\ell}=\frac{1}{2\ell+1}\sum_{m}w^{A}_{\ell,m}w^{C}_{\ell,m}(w^{B}_{\ell,m}w^{D}_{\ell,m})^{*}.$
(42)

In the above formula, $w^{A}_{\ell,m}$ represents a spherical harmonic
coefficient of the inverse noise weights for map $A$.

The covariance presented in this section uses an approximation that works best
when the sky coverage is large. We find it sufficient for our
BOSS${\times}$BOSS analysis, but for FIRAS${\times}$FIRAS and
FIRAS${\times}$BOSS, the formula appears to be inaccurate, possibly due to the
combination of the steep angular index on the foregrounds combined with the
small sky fraction. We instead use a simulated covariance for
FIRAS${\times}$FIRAS and FIRAS${\times}$BOSS.

##  Appendix D Updating the covariance for FIRAS${\times}$BOSS

We calculate a simulated covariance by drawing individual $a_{\ell m}$
amplitudes from our FIRAS foregrounds plus noise model, cosmological signal
model, and galaxy shot noise model. Each of these models is assumed to be
Gaussian and uncorrelated. We then compute simulated full-sky galaxy maps from
the sum of the cosmological and shot noise draws. We also compute simulated
full-sky FIRAS maps from the sum of the foreground and noise realization plus
a painting of [C II] and CIB signal onto the cosmological signal realization.
The simulated FIRAS signal is then convolved with the FIRAS beam, scan, and
pixel window function. Then, the simulated FIRAS and BOSS maps are multiplied
by their respective weights, and the cross-power spectrum between the FIRAS
and BOSS maps is computed. A simulated covariance is calculated from 10,000 of
these realizations. This procedure assumes a specific value for the [C II] and
CIB brightness, but the covariance should change as the [C II] and CIB
magnitude change. So, we use four sets of simulations to create a simulated
covariance that correctly updates its [C II] and CIB contribution to match the
current [C II] and CIB magnitude estimates.

From the assumption of Gaussian power spectra, the rough form of the
FIRAS${\times}$BOSS covariance is

$\langle\Delta C_{b}^{\times}(z_{1},z_{2})\Delta
C_{b}^{\times}(z_{3},z_{4})\rangle\approx\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm
IM}(z_{1},z_{3})C_{b}^{g}(z_{2},z_{4})+C_{b}^{\times}(z_{1},z_{4})C_{b}^{\times}(z_{3},z_{2})],$
(43)

where the partial sky also introduces some level of coupling between different
angular bins, which we compute via our simulations. Schematically, dropping
the redshift indices, the two terms in brackets can be broken down as follows.

$C_{b}^{\rm
IM}C_{b}^{g}\propto[C_{b}^{FG+N}+C_{b}^{\rm[C{\sc\,II}]\times[C{\sc\,II}]}+C_{b}^{\rm
CIB\times CIB}+C_{b}^{\rm[C{\sc\,II}]\times CIB}+C_{b}^{\rm
CIB\times[C{\sc\,II}]}]C_{b}^{g}.$ (44)
$C_{b}^{\times}C_{b}^{\times}\propto[C_{b}^{\rm[C{\sc\,II}]\times
g}+C_{b}^{\rm CIB\times g}][C_{b}^{\rm[C{\sc\,II}]\times g}+C_{b}^{\rm
CIB\times g}],$ (45)

where $FG+N$ refers to the model of Milky Way foregrounds plus noise. Due to
the large foregrounds and noise, the magnitude of the covariance is dominated
by the $C_{b}^{FG+N}C_{b}^{g}$ term. However, we find that the best-fit values
for $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ and $b_{\rm[C{\sc\,II}]}dI_{\rm
CIB}/dz$ are somewhat affected by the values one assumes for them in the
covariance model, so we need to be able to vary their magnitude in our
covariance model. This requires four simulations, which we describe in table
1.

Table 1: The four simulations required for computing the FIRAS${\times}$BOSS
covariance. All use the same BOSS${\times}$BOSS parameters, but the
foreground, thermal noise, [C II], and CIB amplitudes vary. Simulated
Covariances
---
Name |  FG, Noise Amplitude |  $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ |  $b_{\rm[C{\sc\,II}]}dI_{\rm CIB}/dz$ |  $b_{g},A_{SN}$
$\Sigma_{FG+N}$ |  FIRAS${\times}$FIRAS best-fit |  0 |  0 |  BOSS${\times}$BOSS best-fit
$\Sigma_{\rm[C{\sc\,II}]}$ |  0 |  1 |  0 |  BOSS${\times}$BOSS best-fit
$\Sigma_{\rm CIB}$ |  0 |  0 |  1 |  BOSS${\times}$BOSS best-fit
$\Sigma_{\rm[C{\sc\,II}]+CIB}$ |  0 |  1 |  1 |  BOSS${\times}$BOSS best-fit

These four simulations give the following approximate covariances:

$\Sigma_{FG+N}\approx\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}C_{b}^{FG+N}C_{b}^{g}$ (46)
$\Sigma_{\rm[C{\sc\,II}]}\approx\frac{1}{f_{\rm
sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm[C{\sc\,II}]\times[C{\sc\,II}]}C_{b}^{g}+C_{b}^{\rm[C{\sc\,II}]\times
g}C_{b}^{\rm[C{\sc\,II}]\times g}]$ (47) $\Sigma_{\rm
CIB}\approx\frac{1}{f_{\rm sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm CIB\times
CIB}C_{b}^{g}+C_{b}^{\rm CIB\times g}C_{b}^{\rm CIB\times g}]$ (48)
$\Sigma_{\rm[C{\sc\,II}]+CIB}=\Sigma_{\rm[C{\sc\,II}]}+\Sigma_{\rm
CIB}+\frac{1}{f_{\rm sky}\Delta\ell(2\ell+1)}[C_{b}^{\rm[C{\sc\,II}]\times
CIB}+C_{b}^{\rm CIB\times[C{\sc\,II}]}]C_{b}^{g}$ (49)

From the bottom three simulations, we can isolate the term that arises due to
cross-correlations between [C II] and CIB, as follows:

$\Sigma_{\rm[C{\sc\,II}]\times
CIB}=\Sigma_{\rm[C{\sc\,II}]+CIB}-\Sigma_{\rm[C{\sc\,II}]}-\Sigma_{\rm CIB}.$
(50)

The total covariance for arbitrary [C II] and CIB amplitudes can now be
constructed as follows:

$\Sigma=\Sigma_{FG+N}+b_{\rm[C{\sc\,II}]}^{2}I_{\rm[C{\sc\,II}]}^{2}\Sigma_{\rm[C{\sc\,II}]}+b_{\rm[C{\sc\,II}]}^{2}(dI_{\rm
CIB}/dz)^{2}\Sigma_{\rm
CIB}+b_{\rm[C{\sc\,II}]}^{2}I_{\rm[C{\sc\,II}]}(dI_{\rm
CIB}/dz)\Sigma_{\rm[C{\sc\,II}]\times CIB}.$ (51)

This is the covariance used for our cross-power estimator in Section 4.2 and
in our FIRAS validation simulations, described in Appendix E.

##  Appendix E Validating the parametric FIRAS model in the covariance for
FIRAS${\times}$BOSS

We compare two sets of simulations, each performed over both the CMASS and
LOWZ regions. In the first set of simulations, we draw the FIRAS $a_{\ell m}$
amplitudes from the parametric model fit described in Section 4.3.3. In the
second set of simulations, the FIRAS $a_{\ell m}$ amplitudes are drawn from an
empirical FIRAS model that is the measured FIRAS auto-power in bandpowers. In
all cases, the BOSS signal is drawn from the parametric fit of Section 4.3.2,
and both the [C II] and CIB amplitudes are set to zero. We draw 10,000
realizations for each set of simulations and compute the covariance directly
from the simulations. We then find the maximum likelihood solution for each
realization, using a Gaussian likelihood with this simulated covariance (and
also updating the covariance to account for the current [C II] and CIB
estimate according to Equation 51). Note that the empirical FIRAS-model
simulations are not subject to cosmic bias because the real [C II] and CIB
signal present in the empirical FIRAS model is not correlated with the
simulated cosmological clustering realization. Figure 14 shows that the
empirical (blue) and parametric (red) FIRAS simulations lead to near-identical
uncertainties on $b_{\rm[C{\sc\,II}]}I_{\rm[C{\sc\,II}]}$ and
$b_{\rm[C{\sc\,II}]}dI_{\rm CIB}/dz$.

Figure 14: Simulations of fits to the CIB and [C II] amplitude for two
different types of FIRAS models. The red contours show simulations that use
the parametric FIRAS model described in Section 4.3.3. The blue contours show
simulations that use the empirical FIRAS model, where the unbinned FIRAS data
are used directly as the FIRAS model. The left and right plots show the
FIRAS${\times}$CMASS and FIRAS${\times}$LOWZ simulations, respectively.

As a consistency check on our MCMC procedure, we also perform an MCMC analysis
on one of our realizations for the parametric covariance. We verify that the
95 and 68 percent MCMC contours are nearly identical in size and shape to the
95 and 68 percent contours of the histogram of maximum likelihood points from
our 10,000 simulations. This calculation confirms that the parametric variance
model for FIRAS${\times}$FIRAS yield errors consistent with those inferred
from the measured FIRAS${\times}$FIRAS power spectrum (that the parametric
FIRAS${\times}$FIRAS model is sufficiently accurate).

## References

  * Abdalla et al. (2021) Abdalla E., et al., 2021, arXiv e-prints, p. arXiv:2107.01633
  * Ade et al. (2016) Ade P. A. R., et al., 2016, Astronomy & Astrophysics, 594, A13
  * Alam et al. (2015) Alam S., et al., 2015, The Astrophysical Journal Supplement Series, 219, 12
  * Alonso (2018) Alonso D., 2018, NaMaster github, https://github.com/LSSTDESC/NaMaster
  * Alonso et al. (2019) Alonso D., Sanchez J., Slosar A., 2019, Monthly Notices of the Royal Astronomical Society, 484, 4127
  * Anderson et al. (2018) Anderson C. J., et al., 2018, Monthly Notices of the Royal Astronomical Society, 476, 3382
  * Asorey et al. (2012) Asorey J., Crocce M., Gaztañaga E., Lewis A., 2012, Monthly Notices of the Royal Astronomical Society, 427, 1891
  * Bandura et al. (2014) Bandura K., et al., 2014, in Stepp L. M., Gilmozzi R., Hall H. J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9145, Ground-based and Airborne Telescopes V. p. 914522 (arXiv:1406.2288), doi:10.1117/12.2054950
  * Breysse et al. (2021) Breysse P. C., Yang S., Somerville R. S., Pullen A. R., Popping G., Maniyar A. S., 2021, arXiv e-prints, p. arXiv:2106.14904
  * Bull et al. (2015) Bull P., Ferreira P. G., Patel P., Santos M. G., 2015, The Astrophysical Journal, 803, 21
  * CCAT-Prime collaboration et al. (2021) CCAT-Prime collaboration et al., 2021, arXiv e-prints, p. arXiv:2107.10364
  * Campagne et al. (2017) Campagne J. E., Neveu J., Plaszczynski S., 2017, Astronomy & Astrophysics, 602, A72
  * Carilli & Walter (2013) Carilli C. L., Walter F., 2013, Annual Review of Astronomy and Astrophysics, 51, 105
  * Cataldo et al. (2021) Cataldo G., et al., 2021, arXiv e-prints, p. arXiv:2101.11734
  * Chang et al. (2008) Chang T.-C., Pen U.-L., Peterson J. B., McDonald P., 2008, Physical Review Letters, 100, 091303
  * Chang et al. (2010) Chang T.-C., Pen U.-L., Bandura K., Peterson J. B., 2010, Nature, 466, 463
  * Chen (2012) Chen X., 2012, in International Journal of Modern Physics Conference Series. pp 256–263 (arXiv:1212.6278), doi:10.1142/S2010194512006459
  * Cheng et al. (2018) Cheng C., et al., 2018, The Astrophysical Journal, 868, 26
  * Chluba et al. (2021) Chluba J., et al., 2021, Experimental Astronomy, pp 1–40
  * Cleary et al. (2021) Cleary K. A., et al., 2021, arXiv e-prints, p. arXiv:2111.05927
  * Crichton et al. (2021) Crichton D., et al., 2021, arXiv e-prints, p. arXiv:2109.13755
  * Crites et al. (2014) Crites A. T., et al., 2014, in Holland W. S., Zmuidzinas J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII. p. 91531W, doi:10.1117/12.2057207
  * Dawson et al. (2012) Dawson K. S., et al., 2012, The Astronomical Journal, 145, 10
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Di Dio et al. (2013) Di Dio E., Montanari F., Lesgourgues J., Durrer R., 2013, Journal of Cosmology and Astroparticle Physics, 2013, 044
  * Di Dio et al. (2014) Di Dio E., Montanari F., Durrer R., Lesgourgues J., 2014, Annual Review of Astronomy and Astrophysics, 2014, 042
  * Díaz-Santos et al. (2017) Díaz-Santos T., et al., 2017, American Astronomical Society, 846, 32
  * Dizgah et al. (2019) Dizgah A. M., Keating G. K., Fialkov A., 2019, The Astrophysical Journal Letters, 870, L4
  * Doi et al. (2010) Doi M., et al., 2010, The Astronomical Journal, 139, 1628
  * FIRAS Collaboration (1997) FIRAS Collaboration 1997, FIRAS Explanatory Supplement, https://lambda.gsfc.nasa.gov/product/cobe/firas_exsupv4.cfm
  * Fisher et al. (1994) Fisher K. B., Scharf C. A., Lahav O., 1994, Monthly Notices of the Royal Astronomical Society, 266, 219
  * Fixsen et al. (1994) Fixsen D., et al., 1994, The Astrophysical Journal, 420, 457
  * Fixsen et al. (1999) Fixsen D. J., Bennett C., Mather J. C., 1999, The Astrophysical Journal, 526, 207
  * Fukugita et al. (1996) Fukugita M., Ichikawa T., Gunn J. E., Doi M., Shimasaku K., Schneider D. P., 1996, The Astronomical Journal, 111, 1748
  * Gong et al. (2012) Gong Y., Cooray A., Silva M., Santos M. G., Bock J., Bradford C. M., Zemcov M., 2012, The Astrophysical Journal, 745, 49
  * Gong et al. (2017) Gong Y., Cooray A., Silva M. B., Zemcov M., Feng C., Santos M. G., Dore O., Chen X., 2017, The Astrophysical Journal, 835, 273
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, The Astrophysical Journal, 622, 759
  * Gunn et al. (1998) Gunn J. E., et al., 1998, The Astronomical Journal, 116, 3040
  * Gunn et al. (2006) Gunn J. E., et al., 2006, The Astronomical Journal, 131, 2332
  * Guo et al. (2015) Guo H., et al., 2015, Monthly Notices of the Royal Astronomical Society, 446, 578
  * Hamimeche & Lewis (2008) Hamimeche S., Lewis A., 2008, Physical Review D, 77, 103013
  * Hemmati et al. (2017) Hemmati S., Yan L., Diaz-Santos T., Armus L., Capak P., Faisst A., Masters D., 2017, The Astrophysical Journal, 834, 36
  * Hivon et al. (2002) Hivon E., Górski K. M., Netterfield C. B., Crill B. P., Prunet S., Hansen F., 2002, The Astrophysical Journal, 567, 2
  * Hogan & Rees (1979) Hogan C. J., Rees M. J., 1979, Monthly Notices of the Royal Astronomical Society, 188, 791
  * Keating et al. (2016) Keating G. K., Marrone D. P., Bower G. C., Leitch E., Carlstrom J. E., DeBoer D. R., 2016, The Astrophysical Journal, 830, 34
  * Keating et al. (2020) Keating G. K., Marrone D. P., Bower G. C., Keenan R. P., 2020, The Astrophysical Journal, 901, 141
  * Keenan et al. (2021) Keenan R. P., Keating G. K., Marrone D. P., 2021, arXiv e-prints, p. arXiv:2110.02239
  * Kogut et al. (2011) Kogut A., et al., 2011, Journal of Cosmology and Astroparticle Physics, 2011, 025
  * Kovetz et al. (2017) Kovetz E. D., et al., 2017, arXiv preprint arXiv:1709.09066
  * Kovetz et al. (2019) Kovetz E., et al., 2019, Bulletin of the American Astronomical Society
  * Lanusse et al. (2015) Lanusse F., Rassat A., Starck J. L., 2015, Astronomy & Astrophysics, 578, A10
  * Lewis & Challinor (2011) Lewis A., Challinor A., 2011, CAMB: Code for Anisotropies in the Microwave Background (ascl:1102.026)
  * Liu et al. (2017) Liu H., von Hausegger S., Naselsky P., 2017, Physical Review D, 95, 103517
  * Loureiro et al. (2019) Loureiro A., et al., 2019, Monthly Notices of the Royal Astronomical Society, 485, 326
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Madau et al. (1997) Madau P., Meiksin A., Rees M. J., 1997, The Astrophysical Journal, 475, 429
  * Malhotra et al. (1997) Malhotra S., et al., 1997, The Astrophysical Journal Letters, 491, L27
  * Masui et al. (2013) Masui K. W., et al., 2013, The Astrophysical Journal Letters, 763, L20
  * Mather et al. (1993) Mather J. C., Fixsen D. J., Shafer R. A., 1993, in Infrared Spaceborne Remote Sensing. pp 168–180
  * McDonald & Seljak (2009) McDonald P., Seljak U., 2009, Journal of Cosmology and Astroparticle Physics, 10, 007
  * McKinley et al. (2018) McKinley B., et al., 2018, Monthly Notices of the Royal Astronomical Society, 481, 5034
  * Murray (2014) Murray S., 2014, HMF: Halo Mass Function calculator (ascl:1412.006)
  * Nicola et al. (2014) Nicola A., Refregier A., Amara A., Paranjape A., 2014, Physical Review D, 90, 063515
  * Oxholm & Switzer (2021) Oxholm T. M., Switzer E. R., 2021, Physical Review D, 104, 083501
  * Padmanabhan et al. (2007) Padmanabhan N., et al., 2007, Monthly Notices of the Royal Astronomical Society, 378, 852
  * Pen et al. (2009) Pen U.-L., Staveley-Smith L., Peterson J. B., Chang T.-C., 2009, Monthly Notices of the Royal Astronomical Society, 394, L6
  * Planck Collaboration et al. (2014a) Planck Collaboration et al., 2014a, Astronomy & Astrophysics, 566, A55
  * Planck Collaboration et al. (2014b) Planck Collaboration et al., 2014b, Astronomy & Astrophysics, 571, A30
  * Popping et al. (2019) Popping G., Narayanan D., Somerville R. S., Faisst A. L., Krumholz M. R., 2019, Monthly Notices of the Royal Astronomical Society, 482, 4906
  * Pullen et al. (2018) Pullen A. R., Serra P., Chang T.-C., Doré O., Ho S., 2018, Monthly Notices of the Royal Astronomical Society, 478, 1911
  * Salazar-Albornoz et al. (2017) Salazar-Albornoz S., et al., 2017, Monthly Notices of the Royal Astronomical Society, 468, 2938
  * Schöneberg et al. (2018) Schöneberg N., Simonović M., Lesgourgues J., Zaldarriaga M., 2018, Journal of Cosmology and Astroparticle Physics, 10, 047
  * Scott & Rees (1990) Scott D., Rees M. J., 1990, Monthly Notices of the Royal Astronomical Society, 247, 510
  * Serra et al. (2014a) Serra P., Lagache G., Doré O., Pullen A., White M., 2014a, Astronomy & Astrophysics, 570, A98
  * Serra et al. (2014b) Serra P., Lagache G., Doré O., Pullen A., White M., 2014b, Astronomy & Astrophysics, 570, A98
  * Serra et al. (2016) Serra P., Doré O., Lagache G., 2016, The Astrophysical Journal, 833, 153
  * Shang et al. (2012) Shang C., Haiman Z., Knox L., Oh S. P., 2012, Monthly Notices of the Royal Astronomical Society, 421, 2832
  * Silva et al. (2015) Silva M., Santos M. G., Cooray A., Gong Y., 2015, The Astrophysical Journal, 806, 209
  * Smee et al. (2013) Smee S. A., et al., 2013, The Astronomical Journal, 146, 32
  * Smith et al. (2003) Smith R., et al., 2003, Monthly Notices of the Royal Astronomical Society, 341, 1311
  * Stacey et al. (1991) Stacey G. J., Geis N., Genzel R., Lugten J. B., Poglitsch A., Sternberg A., Townes C. H., 1991, The Astrophysical Journal, 373, 423
  * Stacey et al. (2010) Stacey G. J., Hailey-Dunsheath S., Ferkinhoff C., Nikola T., Parshley S. C., Benford D. J., Staguhn J. G., Fiolet N., 2010, The Astrophysical Journal, 724, 957
  * Suginohara et al. (1999) Suginohara M., Suginohara T., Spergel D. N., 1999, The Astrophysical Journal, 512, 547
  * Sun et al. (2019) Sun G., Hensley B. S., Chang T.-C., Doré O., Serra P., 2019, The Astrophysical Journal, 887, 142
  * Switzer (2017a) Switzer E. R., 2017a, The Astrophysical Journal, 838, 82
  * Switzer (2017b) Switzer E. R., 2017b, The Astrophysical Journal, 838, 82
  * Switzer et al. (2013) Switzer E. R., et al., 2013, Monthly Notices of the Royal Astronomical Society, 434, L46
  * Switzer et al. (2015) Switzer E. R., Chang T.-C., Masui K. W., Pen U.-L., Voytek T. C., 2015, The Astrophysical Journal, 815, 51
  * Switzer et al. (2019) Switzer E. R., Anderson C. J., Pullen A. R., Yang S., 2019, The Astrophysical Journal, 872, 82
  * The CONCERTO collaboration et al. (2020) The CONCERTO collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.14246
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, The Astrophysical Journal, 688, 709
  * Tristram et al. (2005) Tristram M., Macías-Pérez J., Renault C., Santos D., 2005, Monthly Notices of the Royal Astronomical Society, 358, 833
  * Vieira et al. (2020) Vieira J., et al., 2020, arXiv e-prints, p. arXiv:2009.14340
  * Visbal et al. (2011) Visbal E., Trac H., Loeb A., 2011, Annual Review of Astronomy and Astrophysics, 8, 010
  * Wang et al. (2021) Wang J., et al., 2021, Monthly Notices of the Royal Astronomical Society, 505, 3698
  * Wolz et al. (2017) Wolz L., et al., 2017, Monthly Notices of the Royal Astronomical Society, 464, 4938
  * Wolz et al. (2021) Wolz L., et al., 2021, arXiv e-prints, p. arXiv:2102.04946
  * Wyithe et al. (2008) Wyithe J. S. B., Loeb A., Geil P. M., 2008, Monthly Notices of the Royal Astronomical Society, 383, 1195
  * Yang et al. (2019) Yang S., Pullen A. R., Switzer E. R., 2019, arXiv preprint arXiv:1904.01180
  * Yang et al. (2021) Yang S., Somerville R. S., Pullen A. R., Popping G., Breysse P. C., Maniyar A. S., 2021, The Astrophysical Journal, 911, 132
  * Zanella et al. (2018) Zanella A., et al., 2018, Monthly Notices of the Royal Astronomical Society, 481, 1976
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, Astronomy & Astrophysics, 556, A2