<figure><img src="image_1.png" /></figure>

# The ALMA Spectroscopic Survey in the HUDF: Constraining Cumulative CO Emission at 1 ≲ z ≲ 4 with Power Spectrum Analysis of ASPECS LP Data from 84 to 115 GHz

Bade D. Uzgil¹,² ID, Chris Carilli¹,³ ID, Adam Lidz⁴ ID, Fabian Walter¹,² ID, Nithyanandan Thyagarajan¹ ID, Roberto Decarli⁵ ID, Manuel Aravena⁶ ID, Frank Bertoldi⁷ ID, Paulo C. Cortes⁸,⁹ ID, Jorge González-López⁶,¹⁰, Hanae Inami¹¹,¹², Gergö Popping² ID, Dominik A. Riechers¹³,¹⁴,<fnref n="18" /> ID, Paul Van der Werf¹⁵, Jeff Wagg¹⁶, and Axel Weiss¹⁷ ID

¹National Radio Astronomy Observatory, Pete V. Domenici Array Science Center, P.O. Box 0, Socorro, NM 87801, USA; buzgil@nrao.edu

²Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

³Cavendish Laboratory, University of Cambridge, 19 J J Thomson Avenue, Cambridge CB3 0HE, UK

⁴Department of Physics & Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA

⁵INAF-Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy

⁶Núcleo de Astronomia de la Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile

⁷Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany

⁸Joint ALMA Observatory—ESO, Av. Alonso de Córdova, 3104, Santiago, Chile

⁹National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

¹⁰Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile

¹¹Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon (CRAL) UMR5574, F-69230 Saint-Genis-Laval, France

¹²Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

¹³Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA

¹⁴Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

¹⁵Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile

¹⁶SKA Organization, Lower Withington Macclesfield, Cheshire SK11 9DL, UK

¹⁷Max-Planck-Institut für Radioastronomie, Auf dem Hügel 71, D-53121 Bonn, Germany

Received 2019 July 29; revised 2019 October 18; accepted 2019 October 25; published 2019 December 9

## Abstract

We present a power spectrum analysis of the ALMA Spectroscopic Survey Large Program (ASPECS LP) data from 84 to 115 GHz. These data predominantly probe small-scale fluctuations ($k = 10–100 h \text{ Mpc}^{-1}$) in the aggregate CO emission in galaxies at $1 \lesssim z \lesssim 4$. We place an integral constraint on CO luminosity functions (LFs) in this redshift range via a direct measurement of their second moments in the three-dimensional (3D) autopower spectrum, finding a total CO shot-noise power $P_{CO,CO}(k_{CO(2-1)}) \le 1.9 \times 10^2 \mu\text{K}^2 (\text{Mpc h}^{-1})^3$. This upper limit ($3\sigma$) is consistent with the observed ASPECS CO LFs in Decarli et al. but rules out a large space in the range of $P_{CO,CO}(k_{CO(2-1)})$ inferred from these LFs, which we attribute primarily to large uncertainties in the normalization $\Phi_*$ and knee $L_*$ of the Schechter-form CO LFs at $z > 2$. Also, through power spectrum analyses of ASPECS LP data with 415 positions from galaxies with available optical spectroscopic redshifts, we find that contributions to the observed mean CO intensity and shot-noise power of MUSE galaxies are largely accounted for by ASPECS blind detections. Finally, we sum the fluxes from individual blind CO detections to yield a lower limit on the mean CO surface brightness at 99 GHz of $\langle T_{CO} \rangle = 0.55 \pm 0.02 \mu\text{K}$, which we estimate represents 68%–80% of the total CO surface brightness at this frequency.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Molecular gas (1073)*

## 1. Introduction

The formation of molecular clouds from atomic hydrogen gas and their subsequent consumption as fuel for star formation are important transitions linking the early stages in the life cycle of the interstellar medium (ISM) to the evolution of galaxies. Obtaining an unbiased and complete measure of the cold gas content and star formation activity of galaxies as functions of cosmic time provides insight into the underlying physical processes that regulate this evolution. With the cosmic star formation rate density (SFRD) well characterized out to $z \sim 3–4$ and rest-frame UV observations setting constraints on the SFRD into the first billion years after the Big Bang (see, e.g., Madau & Dickinson 2014 for a review), there are ongoing efforts to complement this understanding with concurrent trends in the atomic (H I; e.g., Neeleman et al. 2016) and molecular (H$_2$; e.g., Decarli et al. 2019) gas history, particularly during the epoch of galaxy mass assembly at $z \sim 2$, when cosmic star formation activity was approximately 10 times

higher than in the present epoch and more than half of the stellar mass in the universe was accumulated (Madau & Dickinson 2014).

While the cosmic H I gas density has been inferred from observations of damped Ly$\alpha$ systems in quasar spectra at $z \lesssim 5$ (Wolfe et al. 2005), a more direct method is to observe the H I gas in emission via the 21 cm hyperfine transition. The 21 cm experiments have constrained the atomic gas density in cosmic volumes out to $z \sim 0.8$ (Chang et al. 2010; Switzer et al. 2013), and extended surveys are underway (e.g., CHIME, Tianlai, HIRAX, BINGO, and the Ooty Wide Field Array) to push this redshift limit in a continuous range out to $z \sim 3.5$. As the primary science goal of many of these experiments is to use H I as a tracer of large-scale structure in order to measure the imprint of baryon acoustic oscillations, these experiments utilize an observational technique known as line intensity mapping to survey large areas of sky ($\mathcal{O}(10^3–10^4 \text{ deg}^2)$) with coarse angular resolution ($\mathcal{O}(10 \text{ arcmin})$) in a spectral line across a wide fractional bandwidth (30%–60%), resulting in 3D maps of spatially confused line emission throughout

<footnote n="18">¹⁸ Humboldt Research Fellow.</footnote>
cosmological volumes. Line intensity mapping experiments measure the surface brightness fluctuations in the targeted spectral line, as well as any additional line or continuum emission contributing to the aggregate surface brightness at the observed frequencies, via the power spectrum.

Owing to the large collecting areas and wide bandwidths available in existing facilities such as the Atacama Large Millimeter Array (ALMA), Karl G. Jansky Very Large Array (JVLA), and IRAM NOrthern Extended Millimeter Array (NOEMA), the cosmic evolution of molecular gas density has already been measured out to $z \sim 4$—well into the epochs of galaxy mass assembly and peak cosmic star formation history—with various surveys targeting different rotational $J$ transitions of the CO molecule as an $H_2$ gas tracer (Walter et al. 2014, 2016; Pavesi et al. 2018). Unlike the $H\alpha$ intensity mapping experiments, the CO surveys performed blind spectral scans, or so-called molecular deep fields, to build a census of galaxies’ gas content by detecting emission from individual CO-emitting sources that are brighter than the survey’s flux limit. Given the relatively small fields of view and longer baselines of the telescopes employed in these efforts, the molecular deep fields are characterized by survey areas $O(10^0-10^1 \text{ arcmin}^2)$) and angular resolutions ($O(10^0)$ arcsec) well-suited for observing individual galaxies.

The ALMA Spectroscopic Survey Large Program (ASPECS LP) in the *Hubble* Ultra Deep Field (HUDF) is the latest example of a blind spectral scan that has resulted—with 68 hr of total telescope time to scan the full ALMA Band 3 from 84 to 115 GHz—in the tightest blind constraints to date on the evolution of CO luminosity functions (LFs), which directly translate to measurements on the cosmic molecular gas density (Decarli et al. 2019, hereafter D19) over ~12 Gyr of the universe’s history, revealing the levels of accumulation and consumption of molecular gas in galaxies from $z \sim 4$ to the present day. The ASPECS LP targeted an ~4.6 arcmin$^2$ field in a region of the HUDF containing the deepest near-infrared (near-IR) photometric data on the sky (Illingworth et al. 2013; Koekemoer et al. 2013) and ~1500 spectroscopic redshifts for rest-frame optically/UV-selected galaxies (Inami et al. 2017), which facilitated the confirmation and redshift identification of blindly detected line candidates and further enabled the characterization of physical properties such as molecular gas mass, stellar mass, active galactic nucleus (AGN) fraction, metallicity, IR luminosity, and star formation rate (SFR) for all secure detections in the field, as well as for hundreds of fainter sources. Papers from the ASPECS team discuss key results from the ASPECS LP scan in ALMA Band 3, including the observed CO LFs (D19), which also contain a detailed description of the ASPECS survey and ancillary data sets, blind searches for spectral line and continuum detections (González-López et al. 2019, hereafter GL19), MUSE-based CO identifications and demographics of the ASPECS CO sample from spectral energy distribution (SED) modeling (Boogaard et al. 2019, hereafter B19), theoretical perspectives on the cosmic molecular gas density evolution (Popping et al. 2019), ISM properties (Aravena et al. 2019), and stacking analysis with MUSE galaxies in the UDF.

In this paper, we consider the ASPECS LP Band 3 data in the context of a power spectrum analysis. Although we adopt the power spectrum approach used in line intensity mapping experiments, our data set is inherently distinct from those produced by the aforementioned line intensity mapping

experiments, given the marked differences in sky coverage and angular resolution. Our overarching goals, however—to (1) probe the flux from galaxies below the survey’s sensitivity threshold for individual line detections and (2) improve the constraints on the observed cumulative emission (from CO, in our case) by measuring surface brightness fluctuations within the survey volume—are akin to the objectives common throughout the line intensity mapping experimental landscape (see, e.g., Kovetz et al. 2017, for a review), including experiments with goals of mapping the CO intensity field at $2 < z < 3$ to determine the cosmic molecular gas density (Keating et al. 2016; Li et al. 2016). Furthermore, the parallel analysis by the ASPECS team to extract individual CO detections, along with the rich multiwavelength data sets available in the HUDF, provide valuable information to aid in the interpretation of the power spectrum results and enable an exploration of the complementarities between the two approaches.

The organization of this work is as follows. In Section 2, we place ASPECS in the context of a power spectrum analysis, identifying, e.g., the relevant scales that the survey covers in Fourier space. In Section 3, we describe the Band 3 data, as well as details regarding our approach to measuring the power spectrum. Our results on lower limits from blindly detected sources, the three-dimensional (3D) CO autopower spectrum, and statistical analysis of the CO fluctuation data including information from galaxy catalogs are presented in Section 4. In Section 5, we briefly discuss our findings in the framework of a comparison between the power spectrum analysis and the blind line search in recovering the true CO power and comment on the capability of current facilities to measure the CO power at high redshift. Finally, we summarize our findings in Section 6.

Throughout this work, we adopt a cosmological model with $\Omega_M = 0.7$, $\Omega_\Lambda = 0.3$, $\Omega_k = 0$, and $h = H_0/100 = 0.70$.

## 2. ASPECS in the Context of a Power Spectrum Analysis

### 2.1. Mapping Survey Dimensions from Real Space to Fourier Space

The original goal of the ASPECS LP—to reach a sensitivity such that the predicted “knee” of the CO LF could be reached at $z \sim 2''$ (Walter et al. 2016)—was a key driver of the chosen survey parameters,<fnref n="19" /> including total observing time (and, hence, rms sensitivity per beam per channel), spectral resolution $\Delta\nu_{\text{chn}}$, survey bandwidth $\Delta\nu_{\text{BW}}$, array configuration (or synthesized beam size $\Delta\theta_b$), and survey width $\Delta\theta_s$. We do not—and, in some cases, cannot—alter these experimental parameters for the purposes of the power spectrum analysis, except when redefining $\Delta\nu_{\text{chn}}$ and $\Delta\theta_s$, to be explained in more detail below. Upon adopting a target redshift for the observations, $\Delta\nu_{\text{chn}}$, $\Delta\nu_{\text{BW}}$, $\Delta\theta_b$, and $\Delta\theta_s$ can be translated to physical comoving length scales via the standard cosmological relations, and, thus, set the range of distances where statistical correlations between galaxies can be probed. Given that 11 of the 16 secure, blindly detected sources in GL19 with known redshifts in ASPECS correspond to CO(2–1) emitters, we have adopted a target redshift $z_{\text{cen,CO(2-1)}} = 1.315$ to represent the redshift of CO(2–1) emission observed at the band center, $\nu_{\text{cen}} = 99.572$ GHz. We discuss the issue of redshift ambiguities in the CO line emission in Section 2.2.1.

<footnote n="19">¹⁹ See Walter et al. (2016) for more on the rationale behind the opted survey design.</footnote>
### 2.1.1. Real-space Dimensions

The largest physical scale, then, accessible in real space in the line-of-sight dimension, $r_{\parallel,max}$, is determined by the survey's frequency coverage, $\Delta\nu_{BW}$,

$$r_{\parallel,max} = \frac{c}{H_0} \int_{z_{\min}}^{z_{\max}} \frac{dz}{\sqrt{\Omega_{\mathrm{M}}(1+z)^3 + \Omega_{\Lambda}}} \\ = \chi(z_{\max}) - \chi(z_{\min}), \qquad (1)$$

where $\chi(z)$ is the comoving line-of-sight distance to redshift $z$. In the above expression, $z_{\min}$ and $z_{\max}$ correspond to the minimum and maximum redshifts observed at, respectively, the highest and lowest frequencies, $\nu_{\max} = 114.750$ GHz and $\nu_{\min} = 84.278$ GHz, of the survey bandwidth, so that $z_{\min} = \nu_{\mathrm{rest,CO(2-1)}/\nu_{\max} - 1} = 1.009$, $z_{\max} = \nu_{\mathrm{rest,CO(2-1)}/\nu_{\min} - 1} = 1.735$, and $r_{\parallel,max} = 1054.8$ Mpc h$^{-1}$.

Similarly, the channel resolution establishes the smallest physical scale probed in the line of sight, $r_{\parallel,min}$:

$$r_{\parallel,min} = \chi(z_{\mathrm{chn},i+1}) - \chi(z_{\mathrm{chn},i}). \qquad (2)$$

Here $z_{\mathrm{chn},i}$ and $z_{\mathrm{chn},i+1}$ correspond to redshifts of the $i$th and $i+1$th channels at observed frequencies $\nu_i$ and $\nu_{i+1} = \nu_i + \Delta\nu_{\mathrm{chn}}$, so that $z_{\mathrm{chn},i} = \nu_{\mathrm{rest,CO(2-1)}/\nu_i - 1$ and $z_{\mathrm{chn},i+1} = \nu_{\mathrm{rest,CO(2-1)}/(\nu_i + \Delta\nu_{\mathrm{chn}}) - 1$. Because $\Delta\nu_{\mathrm{chn}}$ is a constant across the band, the physical separation between channels increases gradually with redshift. In practice, then, to facilitate computing the Fourier transform, we do not use Equation (2) when converting channel widths from frequency to physical distance. Rather, we define channel separations of equal width in space, dividing the total line-of-sight distance, $r_{\parallel,max}$, by the number of channels, $N_{\mathrm{chn}} = 196$, in the band, yielding $r_{\parallel,min} = 5.38$ Mpc h$^{-1}$. This value is equal to the channel width at $\nu_{\mathrm{cen}}$ and is a reasonable substitute for the true $r_{\parallel,min}$ per channel, due to the modest 12.6% relative change in $r_{\parallel,min}$ from either band edge to the band center. Note that the choice of $N_{\mathrm{chn}} = 196$ reflects the fact that we have imaged the ASPECS data cube while rebinning the native ALMA channel resolution by a factor of 40 (i.e., $\Delta\nu_{\mathrm{chn}} = \Delta\nu_{40\mathrm{chn}}$), compared to the factor of 2 rebinning ($\Delta\nu_{\mathrm{chn}} = \Delta\nu_{2\mathrm{chn}}$) used when imaging the data cube for purposes of CO line searches, etc. The coarser spectral resolution $\Delta\nu_{40\mathrm{chn}} = 0.156$ GHz, with a velocity width $\Delta\nu_{40\mathrm{chn}} \sim 470$ km s$^{-1}$ at $\nu_{\mathrm{cen}}$, ensures that most CO emission is spectrally unresolved throughout the data cube; the median FWHM of the Gaussian line flux profiles of the blindly detected lines in ASPECS is 355 km s$^{-1}$, and the full range of observed FWHMs spans 40.0–617 km s$^{-1}$ (GL19). We favor the larger channel width to avoid significant contributions to the power spectrum from emission lines with FWHM $\gg \Delta\nu_{\mathrm{chn}}$, since we do not attempt to characterize the effect of falsely elongating the observed flux density of these lines from the $\sim$kpc scales of localized emission within the CO-bright galaxy to $\sim$Mpc scales when converting channel widths to cosmological line-of-sight distances.

In the transverse, or on-sky, dimensions, the largest and smallest physical scales accessible to probe CO fluctuations in real space, $r_{\perp,max}$ and $r_{\perp,min}$, are determined by the survey width and synthesized beam size,

$$r_{\perp,max} = D_{\mathrm{A,co}}(z_{\mathrm{cen,CO(2-1)}}) \Delta\theta_{\mathrm{S}}, \qquad (3)$$

$$r_{\perp,min} = D_{\mathrm{A,co}}(z_{\mathrm{cen,CO(2-1)}}) \Delta\theta_{\mathrm{b}}, \qquad (4)$$

where the units of $\Delta\theta_S$ and $\Delta\theta_b$ are in radians, and $D_{\mathrm{A,co}}(z)$, the comoving angular diameter distance at redshift $z$, is equal to $\chi(z)$ for $\Omega_k = 0$.

As the antenna primary beam size grows with observed wavelength, $\Delta\theta_b$ and $\Delta\theta_S$ gradually increase with redshift across the survey bandwidth. As with $r_{\parallel,min}$, we adopt fixed values for each quantity calculated at $\nu_{\mathrm{cen}}$, where change is modest to either band edge. At this frequency, the synthesized beam—an ellipse described by the FWHMs of its major and minor axes, $\Delta\theta_{b,maj}$ and $\Delta\theta_{b,min}$, respectively—for the full data set is $\Delta\theta_b = \Delta\theta_{b,maj} \times \Delta\theta_{b,min} = 1.80 \times 1.48$, corresponding to comoving transverse distances 0.0250 Mpc h$^{-1} \times 0.0205$ Mpc h$^{-1}$. (For reference, $\Delta\theta_b = 2.11 \times 1.67$ at $\nu_{\min}$, and $\Delta\theta_b = 1.51 \times 1.32$ at $\nu_{\max}$.) Expressing the beam area as $A_b = \Delta\theta_{b,maj}\Delta\theta_{b,min}\pi/(4\ln 2)$, we let $r_{\perp,min} = \sqrt{A_b} = 0.0240$ Mpc h$^{-1}$. Note that the data cube for an interferometric image is gridded with rectangular cell sizes $\Delta\theta_{\mathrm{cell}}$ a factor of a few times smaller than $\Delta\theta_b$, chosen such that $\Delta\theta_{\mathrm{cell}}$ represents Nyquist sampling of the longest-baseline visibility data (Taylor et al. 1999). Thus, the smallest transverse dimension present in the data set is actually $\Delta\theta_{\mathrm{cell}} = 0.36 (=0.005)$ Mpc h$^{-1}$ at $\nu_{\mathrm{cen}}$, though there is no information on CO fluctuations contained within physical scales smaller than $\Delta\theta_b$.

At $\nu_{\mathrm{cen}}$, the full width of the survey spans roughly $\Delta\theta_{S,tot} = 2.83$ at a primary beam response cutoff of 20%. The sensitivity profile of the mosaic primary beam implies, however, that the antenna response drops to 50% at $\Delta\theta_{S,HPBW} \approx 2.15$. Beyond this threshold, the rms noise statistics deteriorate rapidly, as indicated by the noise map$^{20}$ in Figure 1 for $\nu_{\mathrm{cen}}$. Not only is the overall rms noise higher for the survey field past the half-power point of the mosaic primary beam, the spatial variation of the rms is also significantly greater in this region, compared to the central $\sim 4$ arcmin$^2$, e.g., where the rms remains mostly between 0.13 and 0.18 Mjy beam$^{-1}$ channel$^{-1}$, gradually reaching 0.3 Mjy beam$^{-1}$ channel$^{-1}$ at the half-power point;$^{21}$ beyond the half-power point, the rms increases from 0.3 to 0.8 Mjy beam$^{-1}$ at the outermost edge of the mosaic, defined by the primary beam cutoff at 20% antenna response. We note as well that results from the blind search for individual CO emitters in ASPECS data (GL19) suggest a lower fidelity (i.e., higher probability of false identification) of line candidates in the survey volume corresponding to $<50\%$ antenna response, so other studies (e.g., CO LF measurements presented in D19) within the ASPECS collaboration have excluded data that lie outside $\Delta\theta_{S,HPBW}$. Thus, we limit our analysis to a square region (shown as the black dotted square in Figure 1) with area $\Delta\theta_S^2 = (1.84 \mathrm{arcmin})^2 = (1.53 \mathrm{Mpc h}^{-1})^2$, chosen to lie within the 50% power threshold at all observed frequencies. For reference, this region encompasses roughly 85% of the volume contained within $\Delta\theta_{S,HPBW}$ and 55% of the volume within $\Delta\theta_{S,tot}$.

### 2.1.2. Fourier-space Dimensions

Since we will be characterizing the CO fluctuation field by its power spectrum, we must relate the relevant physical scales probed in real space (Equations (1)–(4)) to the wavevectors

<footnote>$^{20}$ Noise maps were generated by calculating the rms, in units of mJy beam$^{-1}$, for each pixel, using all surrounding data within a 70 pixel $\times$ 70 pixel box.</footnote>

<footnote>$^{21}$ Here the quoted rms values refer to channel widths with $\Delta\nu_{\mathrm{chn}} = \Delta\nu_{2\mathrm{chn}}$, or Mjy beam$^{-1}$ per 7.81 MHz channel.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Noise maps at $v_{min} = 84.3$ and 91.9 GHz, $v_{cen} = 99.6$ and 107.3 GHz, and $v_{max} = 114.8$ GHz. Solid contours indicate curves of constant rms noise (in units of mJy beam⁻¹). Notable changes in the respective rms at different frequencies are due to an overlap of the frequency bands in the observations (top right panel) and decreasing atmospheric transmission at higher frequencies (bottom row; GL19). The black and red dashed contours show, respectively, the mosaic primary beam</figcaption></figure>

with magnitude $k = \sqrt{k_{\parallel}^2 + k_{\perp}^2}$ that are accessible to ASPECS in Fourier space:

$$k_{\parallel, \min} = \frac{2\pi}{r_{\parallel, \max}} \quad \text{and} \quad k_{\perp, \min} = \frac{2\pi}{r_{\perp, \max}}, \qquad (5)$$

$$k_{\parallel, \max} = \frac{2\pi}{2r_{\parallel, \min}} \quad \text{and} \quad k_{\perp, \max} = \frac{2\pi}{2r_{\perp, \min}}, \qquad (6)$$

where $k_{\parallel, \min}$ and $k_{\perp, \min}$ represent the lowest $k$ modes available in the line-of-sight and transverse dimensions, respectively; note that these fundamental modes map to the largest scales accessible in real space, with frequencies spanning a single oscillation across $\Delta v_{BW}$ and $\Delta \theta_s$. The highest $k$ modes probed by the survey, $k_{\parallel, \max}$ and $k_{\perp, \max}$, correspond to Nyquist frequencies, $k_{\parallel, \text{Nyq}} \propto 1/(2r_{\parallel, \min})$ and $k_{\perp, \text{Nyq}} \propto 1/(2r_{\perp, \min})$, mapping these modes to the smallest physical scales in the survey.

Table 1 summarizes the ASPECS survey parameters adopted for the power spectrum analysis and their mappings to physical dimensions in real and $k$ space. In Figure 2, we indicate the location of $k_{\parallel}$ and $k_{\perp}$ values for ASPECS relative to the predicted total CO(2–1) power—including both clustering and shot-noise contributions—at $z = 1$ from Sun et al. (2018). The ranges of transverse and line-of-sight $k$ modes have important

<table><caption>**Table 1**
Mapping ASPECS Survey Parameters to Real- and Fourier-space Dimensions</caption><tr><td>Survey bandwidth, Δv<sub>BW</sub></td><td>84.278–114.750 GHz</td></tr><tr><td>Channel resolution, Δv<sub>chn</sub></td><td>0.156 GHz</td></tr><tr><td>Survey width, Δθ<sub>s</sub></td><td>1.84</td></tr><tr><td>Beam size, Δθ<sub>b</sub></td><td>1.80 × 1.48</td></tr><tr><td>Central redshift, z<sub>cen,CO(2-1)</sub></td><td>1.315</td></tr><tr><td>r<sub>||,min</sub> ≤ r<sub>||</sub> ≤ r<sub>||,max</sub></td><td>5.38 Mpc h<sup>-1</sup> &lt; r<sub>||</sub> &lt; 1054.8 Mpc h<sup>-1</sup></td></tr><tr><td>r<sub>⊥,min</sub> ≤ r<sub>⊥</sub> ≤ r<sub>⊥,max</sub></td><td>0.0240 Mpc h<sup>-1</sup> &lt; r<sub>⊥</sub> &lt; 1.53 Mpc h<sup>-1</sup></td></tr><tr><td>k<sub>||,min</sub> ≤ k<sub>||</sub> ≤ k<sub>||,max</sub></td><td>0.00596 h Mpc<sup>-1</sup> &lt; k<sub>||</sub> &lt; 0.584 h Mpc<sup>-1</sup></td></tr><tr><td>k<sub>⊥,min</sub> ≤ k<sub>⊥</sub> ≤ k<sub>⊥,max</sub></td><td>4.107 h Mpc<sup>-1</sup> &lt; k<sub>⊥</sub> &lt; 130.900 h Mpc<sup>-1</sup></td></tr></table>

implications to be considered when computing the power spectrum.

For example, the large bandwidth and relatively narrow survey area of ASPECS dictate that the CO brightness fluctuations on physical scales larger than the survey width, $r_{\perp, \max} = 1.53$ Mpc h⁻¹ (i.e., for $k < k_{\perp, \min} = 4.107$ h Mpc⁻¹), will be probed exclusively by $k_{\parallel}$ modes. Thus, the power spectrum measured at $k < 4.107$ h Mpc⁻¹ will be an inherently one-dimensional (1D) measurement, dominated by power from the shorter wavelength, high-$k_{\perp}$ modes projected into the line of sight. These physical scales are important, however, for extracting information about large-scale clustering. Based on models (e.g., Pullen et al. 2013; Sun et al. 2018) for the total
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Total predicted CO(2–1) power from Sun et al. (2018). Here CW13 and G14 refer to different prescriptions used for relating CO luminosities to estimated IR luminosities in the model, based on the compilations in Carilli & Walter (2013) and Greve et al. (2014). Power from clustering and shot noise dominates at $k < 1$ and $k > 1$ h Mpc$^{-1}$, respectively. Vertical dashed lines indicate $k$ values corresponding to ASPECS survey parameters, namely, survey bandwidth $\Delta\nu_{BW}$, channel width, survey width $\Delta\theta_s$ (1′.84), antenna primary beam (0′.98), mosaic spacing (25′′.4), and synthesized beam size $\Delta\theta_b$ (1′′.8 × 1′′.5).</figcaption></figure>

CO power spectrum at the redshift range relevant to this study, we expect any power from galaxy clustering between dark matter halos, $P_{\text{CO,CO}}^{\text{clust}}(k, z)$, to dominate the total CO power spectrum, $P_{\text{CO,CO}}^{\text{tot}}(k, z)$, up to $k \lesssim 1$ h Mpc$^{-1}$ compared to contributions from small-scale clustering of galaxies that share a common host dark matter halo or shot-noise power, $P_{\text{CO,CO}}^{\text{shot}}$ (see Figure 2). If the CO surface brightness fluctuations, $\langle T_{\text{CO}} \rangle$, trace the large-scale clustering of galaxies with some mean bias factor, $\langle b_{\text{CO}} \rangle$, that offsets CO-emitting galaxies from the underlying dark matter distribution, i.e., if

$$P_{\text{CO,CO}}^{\text{clust}}(k) = \langle T_{\text{CO}}(z) \rangle^2 \langle b_{\text{CO}}(z) \rangle^2 P_{\text{m,m}}(k, z), \quad (7)$$

where $P_{\text{m,m}}(k, z)$ is the linear matter power spectrum (appropriate for $k < 0.1$ h Mpc$^{-1}$), then the low-$k$ component of the power spectrum is, in principle, useful for constraining the aggregate CO emission within a given cosmological volume. Note that the units of $P_{\text{CO,CO}}^{\text{clust}}(k, z)$ in Equation (7) are in $\mu\text{K}^2$ (Mpc h$^{-1}$)$^3$ for $\langle T_{\text{CO}}(z) \rangle$ in $\mu\text{K}$, $P_{\text{m,m}}(k, z)$ in (Mpc h$^{-1}$)$^3$, and a dimensionless $b_{\text{CO}}(z)$.

For physical scales smaller than $r_{\perp,\text{max}} = 1.53$ Mpc h$^{-1}$ (i.e., for $k \ge k_{\perp,\text{min}} = 4.107$ h Mpc$^{-1}$), it is clear from Table 1 that $k$ can have contributions from both $k_{\perp}$ and $k_{\parallel}$ as long as $k = \sqrt{k_{\perp}^2 + k_{\parallel}^2}$ is within the range $k = 4.107-130.900$ h Mpc$^{-1}$—we are measuring a full 3D power spectrum in this regime—though the $k_{\perp}$ modes, with wavelengths on the order of $\Delta\theta_b$ up to $\Delta\theta_s$, provide most of the information on the power at these scales. At $z \sim 1$, we expect any power from galaxy–galaxy clustering at $k \gtrsim 4$ h Mpc$^{-1}$ to be buried under a Poissonian shot-noise component, $P_{\text{CO,CO}}^{\text{shot}}(k)$, dominated by bright CO emitters in the survey volume. By restricting our analysis to $k \gtrsim 10$ h Mpc$^{-1}$, where the true power spectrum is expected to be flat, the power spectrum measurement is unaffected by the highly anisotropic ASPECS survey window function.

## 2.2. Inherent Challenges to the Autopower Spectrum Measurement

One of the intrinsic benefits of the autopower spectrum measurement is its sensitivity to intensity fluctuations from all sources, faint and bright, contained within the survey volume. The inclusion in the power spectrum analysis of all flux densities present in the data cube also presents specific challenges to the interpretation of the measured power. For the purpose of the ASPECS power spectrum analysis, a primary concern is the redshift ambiguity of CO emission within the observed survey bandwidth. We also briefly discuss the effects of a possible contribution from continuum emission.

### 2.2.1. Redshift Ambiguity of CO Emission

As with blindly detected individual line candidates, where the redshift of a line candidate without spectroscopically or photometrically confirmed counterparts can be ambiguous, the true redshifts of sources contributing to the intensity fluctuations contained within the ASPECS survey volume are unknown. However, in defining a real- and Fourier-space grid to perform our power spectrum calculations, we have assumed a specific target redshift $z_{\text{cen,CO(2-1)}} = 1.315$—corresponding to the central redshift of CO(2–1) in the ASPECS bandwidth—for the emission. While CO(2–1) is expected to dominate the mean surface brightness at 99 GHz, based on the number of blind detections in ASPECS relative to other line transitions, for example, it is not the only source of spectral line emission present in the survey volume. Figure 2 of D19 illustrates the number of different spectral lines that are, in principle, observable within the ASPECS frequency coverage, emitted from galaxies within the local universe (such as CO(1–0) at $z < 0.37$) and at high redshift (such as any CO transition from $J > 2$ at $z > 2$). Thus, the measured autopower spectrum of the ASPECS data should be interpreted as a sum of the power from surface brightness fluctuations from all relevant CO transitions,

$$P_{\text{CO,CO}}(k) = P_{\text{CO(1-0),CO(1-0)}}(k_{\text{CO(1-0)}}) \\ + P_{\text{CO(2-1),CO(2-1)}}(k_{\text{CO(2-1)}}) \\ + P_{\text{CO(3-2),CO(3-2)}}(k_{\text{CO(3-2)}}) \\ + P_{\text{CO(4-3),CO(4-3)}}(k_{\text{CO(4-3)}}) + \dots, \quad (8)$$

where we truncate the sum, in practice, to include contributions from CO(1–0), CO(2–1), CO(3–2), and CO(4–3), because the ASPECS survey has only resulted in CO blind detections up to $J = 4$ (GL19). Each term in the right-hand side of the above equation is expressed as a function of wavenumber $k_{\text{CO}(J-(J-1))}$, corresponding to the Fourier space defined at the emitted redshift of the respective $J$ transition at the band center. For reference, at $\nu_{\text{cen}} = 99.572$ GHz, CO(1–0), CO(3–2), and CO(4–3) can be emitted from $z = 0.157, 2.470$, and $3.629$, respectively. The $k$ appearing on the left-hand side of Equation (8) is intentionally ambiguous; ultimately, we would like to write $P_{\text{CO,CO}}(k)$ as $P_{\text{CO,CO}}(k_{\text{CO(2-1)}})$. We can convert $k_{\text{CO}(J-(J-1))}$ to $k_{\text{CO(2-1)}}$, defined at $z_{\text{cen,CO(2-1)}}$, using so-called
“distortion” factors (as given in, e.g., Lidz & Taylor 2016),

$$ 
\begin{aligned}
\alpha_{\perp}(z_{\text{cen,CO}(J-(J-1))}) &= \frac{k_{\perp,\text{CO}(2-1)}}{k_{\perp,\text{CO}(J-(J-1))}} \\
&= \frac{D_{\text{A,co}}(z_{\text{cen,CO}(J-(J-1))})}{D_{\text{A,co}}(z_{\text{cen,CO}(2-1)})}
\end{aligned}
\qquad (9) $$

and

$$ 
\begin{aligned}
\alpha_{\parallel}(z_{\text{cen,CO}(J-(J-1))}) &= \frac{k_{\parallel,\text{CO}(2-1)}}{k_{\parallel,\text{CO}(J-(J-1))}} \\
&= \frac{H(z_{\text{cen,CO}(2-1)})}{H(z_{\text{cen,CO}(J-(J-1))})} \\
&\quad \times \frac{(1+z_{\text{cen,CO}(J-(J-1))})}{(1+z_{\text{cen,CO}(2-1)})},
\end{aligned}
\qquad (10) $$

that relate transverse and line-of-sight comoving distances, respectively, between the true emitted redshift of intensity fluctuations, $z_{\text{cen,CO}(J-(J-1))}$, and the adopted redshift $z_{\text{cen,CO}(2-1)}$. In Equation (10), the expression $H(z)$ refers to the Hubble parameter at redshift $z$. Finally, combining Equations (8)–(10), we obtain the total CO power in terms of $k_{\text{CO}(2-1)}$:

$$ P_{\text{CO,CO}}(k_{\text{CO}(2-1)}) = P_{\text{CO}(2-1),\text{CO}(2-1)}(k_{\text{CO}(2-1)} + \sum_{J=1,3,4} \left[ \frac{1}{\alpha_{\perp}(z_{\text{cen,CO}(J-(J-1)})^2 \alpha_{\parallel}(z_{\text{cen,CO}(J-(J-1)}) \times P_{\text{CO}(J-(J-1)),\text{CO}(J-(J-1))} \left( \frac{k_{\perp,\text{CO}(2-1)}}{\alpha_{\perp}(z_{\text{cen,CO}(J-(J-1)})}, \times \frac{k_{\parallel,\text{CO}(2-1)}}{\alpha_{\parallel}(z_{\text{cen,CO}(J-(J-1)})} \right) \right] \qquad (11) $$

The multiplicative pre-factor $1/(\alpha_{\perp}^2 \alpha_{\parallel})$ in the second term on the right-hand side of the above equation represents the ratio of volume probed by the survey in CO(2–1) relative to the other J transitions. The ratio $1/(\alpha_{\perp}^2 \alpha_{\parallel}) > 1$ for $J = 1$, reflecting the fact that the volume probed by CO(1–0) is less than the volume probed by CO(2–1), while the opposite is true for the $J = 3$ and 4 transitions, where $\frac{1}{\alpha_{\perp}^2 \alpha_{\parallel}} < 1$. Contributions to the total measured power from different CO transitions are correspondingly magnified or demagnified when projected into the CO(2–1) frame (Equation (11)).

Note that Equation (11) is only valid for large separations in redshift between the sources of CO emission; if there is overlap between the redshift ranges of the different transitions in the ASPECS survey volume, then Equation (11) will contain cross-terms that represent the cross-power spectrum between the CO transitions that overlap in redshift. For the ASPECS spectral coverage, we point out that there is a small overlap in redshift for the CO(3–2) and CO(4–3) in the survey at $z = 3.011–3.107$ (see Table 1 in D19), but the cross-terms here will be negligible given that the mean redshifts $z_{\text{cen,CO}(3–2)} = 2.470$ and $z_{\text{cen,CO}(4–3)} = 3.629$ are widely separated.

### 2.2.2. Continuum Emission

A search for continuum emission in the ASPECS LP Band 3 data was presented in GL19. This study identified six

continuum sources, with the brightest emission on the order of ~10 µJy, indicating that the continuum level in each channel of the 3 mm cube is negligible for our purposes. To ensure, however, that our power spectrum measurements reflect power from spectral line (CO) fluctuations only and do not contain contributions from the continuum, we perform continuum subtraction on the cube with a linear baseline fit, described in Section 3. This continuum-subtracted cube is used for all power spectrum and cross-power spectrum analyses.

## 3. Data and Methods

The ASPECS LP survey consisted of two blind frequency scans at 1.2 and 3 mm in the HUDF (Beckwith et al. 2006). The methods and subsequent power spectrum analysis presented here utilize the 3 mm observations, which consist of a 17-pointing mosaic over ~4.7 arcmin² at five frequency tunings spanning the full extent of ALMA Band 3. As described in D19 (see their Section 2.2 for more details) of this series, the resulting visibility data were imaged using the CASA task tclean—with natural weighting applied in the uv-plane and frequency rebinning over two of the native 3.91 MHz spectral resolution elements—to produce an image cube with mean rms = $\langle\sigma_{N,2chn}\rangle = 1.96 \times 10^{-4}$ Jy beam⁻¹ channel⁻¹ across all 3935 channels in the cube. Here the rms per channel of the data cube has been inferred by computing the rms in the central 70 × 70 pixels for each channel map; given the lack of known sources in this 70 × 70 pixel wide skewer through the data cube, the rms in this region is expected to be a valid representation of the noise in the cube. (Recall that Figure 1 shows the spatial variation of the rms noise at a number of representative frequencies.) At $\nu_{\text{cen}}$, for example, this mean rms translates to a mean surface brightness sensitivity in units of Jy sr⁻¹ and, via the Rayleigh–Jeans law, a mean brightness temperature $\langle\sigma_{N,2chn}\rangle = 2.76 \times 10^6$ Jy sr⁻¹ = $9.07 \times 10^3$ µK. In the context of a power spectrum analysis, this mean rms can give rise to a spectrally featureless (i.e., “white”) noise power, $P_N$,

$$ P_N = \langle\sigma_{N,2chn}\rangle^2 V_{\text{vox,2chn}}, \qquad (12) $$

where $V_{\text{vox,2chn}}$ refers to the voxel volume defined by the beam area and channel width. At $z_{\text{cen}}$, $V_{\text{vox,2chn}} = \Delta\theta_b^2 \Delta\nu_{2chn} = (0.024 \text{ Mpc h}^{-1})^2 \times (0.27 \text{ Mpc h}^{-1}) = 1.57 \times 10^{-4}$ (Mpc h⁻¹)³ and implies $P_N = 1.29 \times 10^4$ µK² (Mpc h⁻¹)³. This image cube, referred to in this work as $T_{0,4chn}$, has served as the principal data product exploited in a variety of analysis efforts by the ASPECS team, including the identification of blindly detected individual line emitters. The catalog of reliable blind detections—specifically, where the probability that the line is due to noise has been determined in GL19 to be less than 10%—is reproduced in Table 2.

As already discussed in Section 2.1.1, the 7.81 MHz channel width is too fine a spectral resolution for the purposes of the power spectrum analysis. Thus, unless otherwise noted, we have imaged with a frequency rebinning over 40 native spectral resolution elements to obtain an image cube $T_{0,4chn}$ (or $T_0$, hereafter, for brevity) characterized by a lower mean rms, $\langle\sigma_{N,4chn}\rangle = 4.55 \times 10^{-5}$ Jy beam⁻¹ channel⁻¹ $\approx \langle\sigma_{N,2chn}\rangle/\sqrt{20}$, as depicted in the bottom right panel of Figure 3. Note that the noise power remains unchanged as the larger channel width counteracts the change in $\langle\sigma_{N,4chn}\rangle^2$ (see Equation (12)).
<table><caption>**Table 2**
Blind CO Detections in ASPECS LP 3 mm Survey</caption><thead><tr><th>ID<br>(1)</th><th>Line<br>(2)</th><th>ν<sub>obs</sub><br>(GHz)<br>(3)</th><th>Flux<br>(Jy km s<sup>-1</sup>)<br>(4)</th><th>Redshift<br>(5)</th></tr></thead><tbody><tr><td>ASPECS LP-3mm.01<sup>a</sup></td><td>CO(3–2)</td><td>97.584</td><td>1.02 ± 0.04</td><td>2.543</td></tr><tr><td>ASPECS LP-3mm.02<sup>a</sup></td><td>CO(2–1)</td><td>99.513</td><td>0.47 ± 0.04</td><td>1.317</td></tr><tr><td>ASPECS LP-3mm.03<sup>a</sup></td><td>CO(3–2)</td><td>100.131</td><td>0.41 ± 0.04</td><td>2.454</td></tr><tr><td>ASPECS LP-3mm.04<sup>a</sup></td><td>CO(2–1)</td><td>95.501</td><td>0.89 ± 0.07</td><td>1.414</td></tr><tr><td>ASPECS LP-3mm.05<sup>a</sup></td><td>CO(2–1)</td><td>90.393</td><td>0.66 ± 0.06</td><td>1.550</td></tr><tr><td>ASPECS LP-3mm.06<sup>a</sup></td><td>CO(2–1)</td><td>110.038</td><td>0.48 ± 0.06</td><td>1.095</td></tr><tr><td>ASPECS LP-3mm.07<sup>a</sup></td><td>CO(3–2)</td><td>93.558</td><td>0.76 ± 0.09</td><td>2.696</td></tr><tr><td>ASPECS LP-3mm.08<sup>a</sup></td><td>CO(2–1)</td><td>96.778</td><td>0.16 ± 0.03</td><td>1.382</td></tr><tr><td>ASPECS LP-3mm.09</td><td>CO(3–2)</td><td>93.517</td><td>0.40 ± 0.04</td><td>2.698</td></tr><tr><td>ASPECS LP-3mm.10</td><td>CO(2–1)</td><td>113.192</td><td>0.59 ± 0.07</td><td>1.037</td></tr><tr><td>ASPECS LP-3mm.11<sup>a</sup></td><td>CO(2–1)</td><td>109.966</td><td>0.16 ± 0.03</td><td>1.096</td></tr><tr><td>ASPECS LP-3mm.12</td><td>CO(3–2)</td><td>96.757</td><td>0.14 ± 0.02</td><td>2.574</td></tr><tr><td>ASPECS LP-3mm.13</td><td>CO(4–3)</td><td>100.209</td><td>0.13 ± 0.02</td><td>3.601</td></tr><tr><td>ASPECS LP-3mm.14</td><td>CO(2–1)</td><td>109.877</td><td>0.35 ± 0.05</td><td>1.098</td></tr><tr><td>ASPECS LP-3mm.15</td><td>CO(2–1)</td><td>109.971</td><td>0.21 ± 0.03</td><td>1.096</td></tr><tr><td>ASPECS LP-3mm.16</td><td>CO(2–1)</td><td>100.503</td><td>0.08 ± 0.01</td><td>1.294</td></tr></tbody></table>

Note. (1) Catalog ID. (2) Identified line transition. (3) Observed frequency at line center. (4) Integrated line flux from Table 6 of GL19. (5) Redshift of observed CO transition.

After imaging and applying a primary beam correction to correct flux densities for the effect of the mosaic sensitivity pattern, we estimate and subtract any possible continuum emission by running CASA task imcontsub in the full cube to ensure that all surface brightness fluctuations in $T_0$ are due to spectral emission. The continuum was approximated using a linear baseline fit across all channels to prevent introducing artificial spectral structure. We inspected rms levels and spectra at random positions in the cube before and after continuum subtraction, finding negligible (less than 0.1%–1%) change in both quantities, confirming our expectations based on GL19 (see Section 2.2.2).

For the purposes of the power spectrum analysis, however, we do not work directly with $T_0$ or $T_{0,2chn}$. That is, we do not assess the level of astrophysical signal in the 3 mm data set by taking the (auto)power spectrum of $T_0$, $P_{T_0, T_0}(k)$, defined as

$$ \langle T_0^*(\mathbf{k}) T_0(\mathbf{k}') \rangle \equiv (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_{T_0, T_0}(k), \quad (13) $$

where $(2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') = \int d^3\mathbf{x} e^{-i(\mathbf{k} - \mathbf{k}') \cdot \mathbf{x}}$ is a Dirac delta function. Explicitly, $P_{T_0, T_0}(k)$ is the 3D average of the Fourier transform of the two-point correlation function, $\xi(|\mathbf{x} - \mathbf{x}'|) = \xi(r)$: $\int d^3\mathbf{r} e^{-i\mathbf{k}\cdot\mathbf{r}} \xi(r)$. We note, however, that the power spectrum in Equation (13) can include contributions from astrophysical signal at the target and/or other redshifts, as well as instrument noise, characterized by $P_N$. In Sections 2.2.1 and 2.2.2, we explained why the main astrophysical source of surface brightness fluctuations in the data cube is expected to be the CO line transitions (from $J = 1, 2, 3, 4$). As our estimated noise power is 2–3 orders of magnitude greater than the predicted CO power spectrum signal, the instrument noise introduces a significant bias in the measured power spectrum throughout all $k$ probed by the survey. Thus, in this low signal-to-noise (S/N) regime, we seek a way to remove the noise bias from the data in order to accurately measure the CO signal; we

avoid subtracting this noise-bias term from the data based on independent estimates of $P_N$ that may not reflect the true noise amplitude in the data or exhibit deviations from Gaussianity that would, for example, invalidate Equation (12), which assumes that the measured rms describes a white-noise random field.

Dillon et al. (2014) demonstrated that it is possible—and, indeed, preferable when the expected S/N is subunity, as in the current generation of reionization-era 21 cm intensity mapping experiments discussed in their paper—to remove the noise bias by computing the cross-power spectrum of two data cubes, e.g., $T_I$ and $T_{II}$, that are derived as subsets of the original data cube, $T_0$, in a manner that preserves the real and Fourier spaces sampled by $T_0$. For ASPECS data, we can split the original CASA measurement set$^{22}$ $\tilde{T}_0$ into two subsets, such that the sum of visibilities in $\tilde{T}_I$ and $\tilde{T}_{II}$ gives $\tilde{T}_0$, i.e., $\tilde{T}_I + \tilde{T}_{II} = \tilde{T}_0$, in order to produce the corresponding image cubes $T_I$ and $T_{II}$. The working assumption here is that the cross-power spectrum between $T_I$ and $T_{II}$, $P_{T_I, T_{II}}(k)$, given by

$$ \langle T_I^*(\mathbf{k}) T_{II}(\mathbf{k}') \rangle \equiv (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_{T_I, T_{II}}(k), \quad (14) $$

contains only astrophysical signal and, in principle, any residual correlated noise; random noise present in each cube will be uncorrelated and produce zero mean signal in the cross. Hereafter, we refer to $P_{T_I, T_{II}}(k)$ as the noise-bias-free power spectrum.

Errors on $P_{T_I, T_{II}}(k)$, $\delta P_{T_I, T_{II}}(k)$, can be similarly evaluated by first creating additional subsets $\tau$ of two visibility data sets from each parent visibility data set, $\tilde{T}_I$ or $\tilde{T}_{II}$, such that $\tilde{\tau}_I + \tilde{\tau}_2 = \tilde{T}_I$ and $\tilde{\tau}_3 + \tilde{\tau}_4 = \tilde{T}_{II}$. Then, the cross-power spectrum (performed, again, in the image domain) between the mathematical differences of each pair is computed to yield the error on $P_{T_I, T_{II}}(k)$,

$$ \delta P_{T_I, T_{II}}(k) = \gamma_N \langle (\tau_1(\mathbf{k}) - \tau_2(\mathbf{k}))^* (\tau_3(\mathbf{k}) - \tau_4(\mathbf{k})) \rangle. \quad (15) $$

In this scheme, the purpose of differencing the two data cubes derived from either $T_I$ or $T_{II}$ is to remove signal from $T_I$ or $T_{II}$, respectively, such that the mathematical difference represents noise-only data. Then, one can compute the cross-power spectrum between the pair of differences to remove the noise bias in the noise-only data, yielding a so-called noise-bias-free error on the noise-bias-free power spectrum, $P_{T_I, T_{II}}(k)$. Furthermore, we can create two additional realizations of $\delta P_{T_I, T_{II}}(k)$ by reordering the differences and obtain a final average error, $\langle \delta P_{T_I, T_{II}}(k) \rangle$, as follows:

$$ \langle \delta P_{T_I, T_{II}}(k) \rangle = \frac{\gamma_N}{3} \left( \langle (\tau_1(\mathbf{k}) - \tau_2(\mathbf{k}))^* (\tau_3(\mathbf{k}) - \tau_4(\mathbf{k})) \rangle \right. \\ \left. + \langle (\tau_1(\mathbf{k}) - \tau_3(\mathbf{k}))^* (\tau_2(\mathbf{k}) - \tau_4(\mathbf{k})) \rangle \right. \\ \left. + \langle (\tau_1(\mathbf{k}) - \tau_4(\mathbf{k}))^* (\tau_2(\mathbf{k}) - \tau_3(\mathbf{k})) \rangle \right). \quad (16) $$

The pre-factor, $\gamma_N$, is determined by relating the expected (and actual) noise properties of $\tau_1$, $\tau_2$, $\tau_3$, and $\tau_4$ to the parent data cubes $T_I$ and $T_{II}$ and grandparent data cube $T_0$. Specifically, in the case that only thermal noise is present in the data, then—as long as the number of visibilities in $\tilde{T}_0$ is divided equally among $\tilde{T}_I$ and $\tilde{T}_{II}$, and the weights (determined by antenna system temperatures) on the visibility data are not dramatically

<footnote>a Source is classified as extended in GL19.</footnote>

<footnote>$^{22}$Throughout this paper, we use the tilde ($\sim$) symbol to denote a visibility data set $\tilde{T}$ used to generate an image cube $T$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top: tree diagram illustrating the derivation of visibility data subsets $\tilde{T}_I$, $\tilde{T}_{II}$, $\tilde{\eta}_1$, $\tilde{\eta}_2$, $\tilde{\eta}_3$, and $\tilde{\eta}_4$ from the original visibility data set $\tilde{T}_0$, used to determine the noise-bias-free power spectrum, $P_{\tilde{\eta}_1, \tilde{\eta}_2}(k)$, and corresponding error $\delta P_{\tilde{\eta}_1, \tilde{\eta}_2}(k)$. The expected mean rms of the resulting image generated from each visibility data set is also labeled. Bottom left: beam major (solid curves) and minor (dashed curves) axes as a function of observed frequency for image subsets $\tau_1$, $\tau_2$, $\tau_3$, and $\tau_4$. Numbers in the upper x-axis refer to the channel number of the data cube, which has been imaged with a factor of 40 rebinning in frequency. Bottom right: rms as a function of observed frequency. The upper x-axis is the same as in the bottom left panel.</figcaption></figure>

different in one subset compared to the other—the resulting mean rms in $T_I$ and $T_{II}$ will be equal to $\sqrt{2} \langle \sigma_{N,40chn} \rangle$, where $\langle \sigma_{N,40chn} \rangle$ refers to the mean rms in $T_0$, derived earlier in this section; and, the cross-power spectrum (Equation (14)) will have a noise covariance equal to $2 \langle \sigma_{N,40chn} \rangle^2$. Similarly, if the same conditions hold for the splitting of CASA measurement sets corresponding to $\tilde{T}_I$ into $\tilde{\eta}_1$ and $\tilde{\eta}_2$ and $\tilde{T}_{II}$ into $\tilde{\eta}_3$ and $\tilde{\eta}_4$, then each image $\tau_1$ through $\tau_4$ will have a mean rms = $2 \langle \sigma_{N,40chn} \rangle$. Figure 3 (bottom right panel), which shows that the rms at all $\nu_{\text{obs}}$ in $\tau_1$ through $\tau_4$ (colored curves) matches well the rms level representing $2 \langle \sigma_{N,40chn} \rangle$ (black dotted curve), confirms the assumptions that our noise is predominantly thermal, the visibility data sets have been divided equally among the subsets, and visibility weights do not differ significantly among the subsets. (We point out that the rms in $\tau_4$ is slightly higher than the black dotted curve at some frequency intervals and have verified that this discrepancy is due to a smaller number of visibilities entering into this subset at those frequencies.) Then, the images representing mathematical differences ($\tau_1 - \tau_2$) and ($\tau_3 - \tau_4$) will each have a mean rms = $2\sqrt{2} \langle \sigma_{N,40chn} \rangle$, and the

resulting cross-power spectrum (Equation (15)) will have a noise covariance equal to $8 \langle \sigma_{N,40chn} \rangle^2$. Therefore, we find $\gamma_N = 0.25$. The process is visualized as a tree diagram in the top panel of Figure 3.

In practice, we forgo the step of dividing $\tilde{T}_0$ into subsets $\tilde{T}_I$ and $\tilde{T}_{II}$ and begin by dividing $\tilde{T}_0$ into quarters $\tilde{\eta}_1$ through $\tilde{\eta}_4$. The CASA measurement set represented by $\tilde{T}_0$ contains visibility data corresponding to 11 frequency ranges, or spectral windows,<fnref n="23" /> that have been stitched together from all available data so that the spectral windows in $T_0$ are ordered from lowest to highest observed frequency, with no overlapping regions. (Please refer to Section 2.4 of Walter et al. 2016 for more details on the construction of $T_0$.) For each spectral window in $T_0$, there are multiple blocks of visibility data corresponding to the various execution blocks scheduled by ALMA for observing. Each block, in turn, is typically comprised of eight to nine scans that repeat 17 times to cover the entire spatial area of the mosaic. Using the CASA task *split*, we select and

<footnote n="23">²³ We do not refer to the four 1.875 GHz spectral windows of the ALMA sidebands.</footnote>
distribute these scans evenly among four subsets. We repeat these steps for every block of scans and every spectral window and merge visibilities in the four subsets using the CASA task `concat`—followed by `statwt`, for a homogenous weighting system in the concatenated data—to produce the subsets $\tilde{\tau}_1, \tilde{\tau}_2, \tilde{\tau}_3$, and $\tilde{\tau}_4$. Our choice of dividing $T_0$ this way guards against possible frequency and/or temporal biases, as each subset contains the full range of frequency (84–115 GHz)—required, moreover, so that all image cubes probe the same line-of-sight distance—and time (2016 December 2–21) covered by the observations. The splittings have also resulted in small (i.e., less than the 0.′/36 cell size in the gridded image cubes) variations in beam sizes between subsets (see bottom left panel in Figure 3), which is necessary when taking the cross-power spectrum or performing other mathematical operations, like subtraction, between any two images.<fnref n="24" />

# 4. Results

## 4.1. Limits from Detected Sources

### 4.1.1. Mean Surface Brightness at 99 GHz

A direct measurement of the mean CO surface brightness, $\langle T_{CO} \rangle$, across observed frequencies $\nu_{obs} = 84.3–114.8$ GHz provides an empirical point of comparison to model predictions in the context of CO intensity mapping experiments at moderate redshifts, and also places a constraint on foreground emission for cosmic microwave background (CMB) experiments aiming to map spectral distortions at high redshift. We repeat the analysis performed for the ASPECS-Pilot program by Carilli et al. (2016) to place a lower limit on $\langle T_{CO} \rangle$ at $\nu_{obs} = 99$ GHz, based on blindly detected sources in the ASPECS LP 3 mm survey.

Following Carilli et al. (2016), we consider the aggregate emission from all observed CO transitions that contribute to the mean sky brightness at 99 GHz. We begin by summing the line fluxes of the 16 blindly detected CO emission line candidates reported in GL19 to obtain a total CO flux of $6.91 \pm 0.19$ Jy km s$^{-1}$ or, equivalently, $(2.28 \pm 0.06) \times 10^6$ Jy Hz. Dividing this total flux by $\Delta\nu_{BW}$ yields a total mean CO flux density $\langle S_\nu \rangle = (7.53 \pm 0.21) \times 10^{-5}$ Jy. Finally, to derive a mean surface brightness in units of $\mu$K, we apply the Rayleigh–Jeans approximation, $\langle T_{CO} \rangle \sim 1360 \langle S_\nu \rangle \lambda_{obs}^2 / A_S^{\text{blind}} = 0.55 \pm 0.02 \ \mu$K, where $\lambda_{obs}$ is the observed wavelength (in units of cm) and $A_S^{\text{blind}}$ is the survey area (in units of arcsec$^2$) utilized in the line search, corresponding to the region of the mosaic where primary beam attenuation is less than 20%: $1.69 \times 10^4$ arcsec$^2$ at 99 GHz. Following the prescription in Moster et al. (2011), we estimate a 19.5% relative uncertainty on $\langle T_{CO} \rangle$ due to cosmic variance in the pencil beam survey by combining the fractional uncertainties<fnref n="25" /> calculated for each identified line entering into the above flux sum, given the survey depth in stellar mass, mean redshift, and survey volume probed by the respective J transition. Because Moster et al. (2011) estimated cosmic variance as the product of galaxy bias

and the dark matter cosmic variance, their prescription is strictly applicable here in the case where the galaxy bias is identical to the bias $b_{CO}$ of CO emission with respect to the matter density field.

For ASPECS-Pilot, which consisted of a single $\sim 1$ arcmin$^2$ pointing with the same spectral coverage as ASPECS LP, $\langle T_{CO} \rangle$ was found to be $0.94 \pm 0.09$ $\mu$K at 99 GHz (Carilli et al. 2016), which is a factor of 1.72 times greater than reported here for ASPECS LP. Since the time of publication of that analysis, however, four of the 10 line candidates reported by ASPECS-Pilot (namely, 3 mm.4, 3 mm.7, 3 mm.8, and 3 mm.9 in Table 2 of Walter et al. 2016) have been reclassified as “unconfirmed”—i.e., likely spurious, given their narrow line widths—sources based on the improved line search algorithms developed in GL19 and are excluded from the present analysis. Additionally, two of the ASPECS-Pilot line candidates (3 mm.6 and 3 mm.10) are outside the ASPECS LP survey coverage and similarly excluded. Thus, when including emission from only the four remaining confirmed sources from the original 10 sources listed in Walter et al. (2016), one finds that the total observed CO flux scales linearly with the decrease in observed survey area, resulting in a revised $\langle T_{CO} \rangle = 0.55 \pm 0.05$ $\mu$K for ASPECS-Pilot, consistent with our new measurement.

It is important to note that the measurement of $\langle T_{CO} \rangle$ presented here is considered a lower limit because the blind detections represent only a fraction of the total CO emission in the ASPECS LP survey volume; the fraction recovered by blind detections is determined by the sensitivity limit of the survey and the shape of the relevant CO LFs. We compute the mean CO surface brightness based on the observed CO(2–1), CO(3–2), and CO(4–3) LFs for ASPECS LP presented in D19 as follows:

$$ \langle T_{CO(J-(J-1))} \rangle = \int d \log_{10} L_{CO(J-(J-1))} \times \Phi(L_{CO(J-(J-1))}) \frac{L_{CO(J-(J-1))}}{4\pi D_L^2} y D_{A,\text{co}}^2, \quad (17) $$

where $D_L$, y, and $D_{A,co}$ refer, respectively, to the luminosity distance, the derivative of the comoving radial distance with respect to the observed frequency (i.e., $y = d\chi/d\nu = \lambda_{\text{rest}}(1+z)^2/H(z)$), and the comoving angular diameter distance and are all evaluated at $z_{\text{cen,CO(J-(J-1))}}$. Here $\Phi(L_{CO(J-(J-1))})$ is originally expressed as a function of the integrated source brightness temperature, $L'_{CO(J-(J-1))}$ (in units of K km s$^{-1}$ pc$^2$), $\Phi(L_{CO(J-(J-1))}')$, and is given in the logarithmic Schechter form

$$ \log_{10} \Phi(L'_{CO(J-(J-1))}) = \log_{10} \Phi_* + \alpha \log_{10} \left( \frac{L_{CO(J-(J-1))}}{L'_{CO(J-(J-1))*}} \right) \\ \times \frac{1}{\ln 10} \frac{L_{CO(J-(J-1))}'}{L'_{CO(J-(J-1))*}} + \log_{10}(\ln 10). \quad (18) $$

We convert from $L'_{CO(J-(J-1))}$ to $L_{CO(J-(J-1))}$ (in units of solar luminosity) via

$$ L_{CO(J-(J-1))} = 3 \times 10^{-11} \nu_{\text{rest,CO(J-(J-1))}}^3 L'_{CO(J-(J-1))} \quad (19) $$

from Carilli & Walter (2013). Fits to the LF data have yielded Schechter parameters $\alpha$, $\Phi_*$, and $L'_{CO(J-(J-1))*}$, with

<footnote n="24">²⁴ Splitting $\tilde{T}_0$ must be done in a way that preserves the real and Fourier spaces probed by $T_0$. For example, if $\tilde{T}_0$ were split into two sets $\tilde{T}_1$ and $\tilde{T}_{11}$ that contained visibilities from the first and second half of the channels, respectively, in $T_0$, then $\tilde{T}_1 + \tilde{T}_{11} = \tilde{T}_0$ would still hold, but the images $T_1$ and $T_{11}$ would each probe only half of the volume in $T_0$.</footnote>

<footnote n="25">²⁵ The relative uncertainty on the mean CO(2–1), CO(3–2), and CO(4–3) surface brightnesses due to cosmic variance is 17%, 23%, and 54% for minimum stellar masses probed of $6 \times 10^9$, $2 \times 10^{10}$, and $3 \times 10^{10}$ M$_\odot$, respectively.</footnote>
<table><caption>**Table 3**
CO LF Schechter Parameters from D19</caption><thead><tr><th>Line</th><th>Redshift</th><th>α</th><th>log<sub>10</sub>Φ<sub>*</sub><br/>[log<sub>10</sub><br/>(Mpc<sup>-3</sup> dex<sup>-1</sup>)]</th><th>log<sub>10</sub>L'<sub>*</sub><br/>[log<sub>10</sub><br/>(K km s<sup>-1</sup> pc<sup>2</sup>)]</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th></tr></thead><tbody><tr><td>CO(2–1)</td><td>1.43</td><td>-0.2 (fixed)</td><td>-2.79<sup>+0.09</sup><sub>-0.09</sub></td><td>10.09<sup>+0.10</sup><sub>-0.09</sub></td></tr><tr><td>CO(3–2)</td><td>2.61</td><td>-0.2 (fixed)</td><td>-3.83<sup>+0.13</sup><sub>-0.12</sub></td><td>10.60<sup>+0.20</sup><sub>-0.15</sub></td></tr><tr><td>CO(4–3)</td><td>3.80</td><td>-0.2 (fixed)</td><td>-3.43<sup>+0.19</sup><sub>-0.22</sub></td><td>9.98<sup>+0.22</sup><sub>-0.14</sub></td></tr></tbody></table>

**Note.** (1) Line transition. (2) Mean redshift of LF redshift bin. (3) Faint-end slope parameter in Equation (18). (4) Normalization parameter in Equation (18). (5) Characteristic luminosity parameter in Equation (18).

uncertainties summarized in Table 3. Note that the faint-end slope, α, has been fixed at α = -0.2 for all LFs.

Integrating the LFs (Equation (17)) from an upper luminosity limit $L'_{\text{upp}} = 10^{12}$ K km s$^{-1}$ pc$^2$ down to the mean $7\sigma$ line sensitivity$^{26}$ $L'_{\text{min},7\sigma}$ in the respective redshift interval covered by each CO transition, which reflects the ASPECS LP detection threshold,$^{27}$ yields a mean total surface brightness $\langle T_{\text{CO}} \rangle_{\text{LF},7\sigma} = 0.49-1.78$ µK, where the quoted range reflects the uncertainty in the LF parameters; please see Table 4 for a breakdown of the inferred $\langle T_{\text{CO}} \rangle$ by J transition. Extending the lower limit of integration down to $L'_{\text{min}} = 10^8$ K km s$^{-1}$ pc$^2$ at all redshifts implies a total mean surface brightness of $\langle T_{\text{CO}} \rangle_{\text{LF}} = 0.72-2.24$ µK. Therefore, we estimate that our blind detections represent $\langle T_{\text{CO}} \rangle_{\text{LF},7\sigma} / \langle T_{\text{CO}} \rangle_{\text{LF}} = 68.1\% - 79.5\%$ of the total CO surface brightness at this observed frequency.

### 4.1.2. CO Shot-noise Power

As with the limit on mean CO surface brightness, the blindly detected sources in Table 2 can also be used to place a lower limit on the expected CO shot-noise power.

The total CO shot-noise power from only the detected sources, $[P_{\text{CO,CO}}^{\text{shot}}(k_{\text{CO}(2-1)})]_{\text{det}}$, will contain contributions from galaxies emitting in the observed transitions $J = 2, 3$, and $4$,

$$[P_{\text{CO,CO}}^{\text{shot}}(k_{\text{CO}(2-1)})]_{\text{det}} = [P_{\text{CO(2-1),CO(2-1)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}} \\ + [P_{\text{CO(3-2),CO(3-2)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}} \\ + [P_{\text{CO(4-3),CO(4-3)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}}, \quad (20)$$

where $[P_{\text{CO(3-2),CO(3-2)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}}$ and $[P_{\text{CO(4-3),CO(4-3)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}}$ have been converted to the CO(2–1) frame using Equation (11). Each term on the right-hand side of Equation (20) can be determined analytically by summing the $N$ individual line fluxes per the expression

$$\sum_{i=1}^{N} \frac{1}{V_S} \left( \frac{L_{\text{CO}(J-(J-1))}}{4\pi D_L^2} y D_{\text{A,co}}^2 \right)^2, \quad (21)$$

<table><caption>**Table 4**
Mean CO Surface Brightness Inferred from Schechter-form LFs</caption><thead><tr><th>Line</th><th>⟨T<sub>CO(J-(J-1))</sub>⟩<br/>⟨L'<sub>min</sub> = L'<sub>min,7σ</sub>⟩<br/>(μK)</th><th>⟨T<sub>CO(J-(J-1))</sub>⟩<br/>⟨L'<sub>min</sub> = 10<sup>8</sup> K km s<sup>-1</sup> pc<sup>2</sup>⟩<br/>(μK)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th></tr></thead><tbody><tr><td>CO(2–1)</td><td>0.53<sup>+0.32</sup><sub>-0.21</sub></td><td>0.72<sup>+0.38</sup><sub>-0.25</sub></td></tr><tr><td>CO(3–2)</td><td>0.25<sup>+0.31</sup><sub>-0.12</sub></td><td>0.29<sup>+0.33</sup><sub>-0.13</sub></td></tr><tr><td>CO(4–3)</td><td>0.12<sup>+0.25</sup><sub>-0.08</sub></td><td>0.20<sup>+0.32</sup><sub>-0.11</sub></td></tr></tbody></table>

**Notes.** (1) Line transition. (2) Mean CO surface brightness calculated by integrating Equation (17) with lower and upper limits of integration $L'_{\text{min}} = L'_{\text{min},7\sigma}$ and $L'_{\text{upp}} = 10^{12}$ K km s$^{-1}$. (3) Same as column (2) but for $L'_{\text{min}} = 10^8$ K km s$^{-1}$ pc$^2$.

where $V_S$ refers to the survey volume at $z_{\text{cen,CO}(J-(J-1))}$. Note that the above expression for shot-noise power has units of surface brightness squared times volume ($\mu K^2$ (Mpc h$^{-1}$)$^3$) and is equal to the same value at all $k$, appropriate for a Poisson sampling of galaxies.

Starting with the CO(2–1), CO(3–2), and CO(4–3) source fluxes from GL19, reported in Table 2, we find, for the entire ASPECS LP 3 mm survey volume used in the blind search, lower limits on the expected shot-noise power of $[P_{\text{CO(2-1),CO(2-1)}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}} = 63.64$, $[P_{\text{CO(3-2),CO(3-2)}}^{\text{shot}}(k_{\text{CO(3-2)}})]_{\text{det}} = 98.49$, and $[P_{\text{CO(4-3),CO(4-3)}}^{\text{shot}}(k_{\text{CO(4-3)}})]_{\text{det}} = 1.05 \mu K^2$ (Mpc h$^{-1}$)$^3$, respectively. For the cropped region (black dotted square in Figure 1) corresponding to the volume used in the power spectrum analysis, we find that the CO(2–1), CO(3–2), and CO(4–3) line emitters each give rise to a respective shot-noise power of 73.99, 71.06, and $1.21 \mu K^2$ (Mpc h$^{-1}$)$^3$. The slightly higher shot-noise power predicted in the cropped region for CO(2–1) and CO(4–3) is due to the decrease in volume after the crop; the CO(3–2) shot-noise power decreases due to the fact that two of the five detected sources are located outside of the boundary of the cropped region. After converting the CO(3–2) and CO(4–3) shot-noise power into the CO(2–1) frame, we sum each contribution to arrive at a total shot-noise power at $z_{\text{cen}} = 1.315$ arising from the blind detections: $[P_{\text{CO,CO}}^{\text{shot}}(k_{\text{CO(2-1)}})]_{\text{det}} = 118.45$ and $113.24 \mu K^2$ (Mpc h$^{-1}$)$^3$ for the full survey and cropped region, respectively.

Finally, we estimate the expected shot-noise power based on the D19 CO LFs,

$$P_{\text{CO(J-(J-1)),CO(J-(J-1))}}^{\text{shot}} = \int d \log_{10} L_{\text{CO(J-(J-1))}} \\ \times \Phi(L_{\text{CO(J-(J-1))}}) \left[ \frac{L_{\text{CO(J-(J-1))}}}{4\pi D_L^2} y D_{\text{A,co}}^2 \right]^2, \quad (22)$$

and find that the detected sources (i.e., integrating Equation (22) down to the relevant $7\sigma$ line sensitivity limit) recover 95.2%–97.7% of $P_{\text{CO(2-1),CO(2-1)}}^{\text{shot}}(k_{\text{CO(2-1)}})$, 98.6%–99.7% of $P_{\text{CO(3-2),CO(3-2)}}^{\text{shot}}(k_{\text{CO(3-2)}})$, and 84.9%–96.0% of $P_{\text{CO(4-3),CO(4-3)}}^{\text{shot}}(k_{\text{CO(4-3)}})$. In total, the recovered fraction is 96.0%–98.6% of $P_{\text{CO,CO}}^{\text{shot}}(k_{\text{CO(2-1)}})$.

<footnote>$^{26}$ For reference, the mean $7\sigma$ line sensitivity for ASPECS in CO(2–1), CO(3–2), and CO(4–3), is $2.68 \times 10^9$ K km s$^{-1}$ pc$^2$ (9.85 $\times 10^5$ L$_{\odot}$) $3.70 \times 10^9$ K km s$^{-1}$ pc$^2$ (4.58 $\times 10^6$ L$_{\odot}$), and $3.93 \times 10^9$ K km s$^{-1}$ pc$^2$ (1.15 $\times 10^7$ L$_{\odot}$), respectively.</footnote>

<footnote>$^{27}$ The S/N threshold S/N $\ge 6.8$ applied to the catalog of all possible line candidates (including candidates down to low S/N) yields the 16 high-fidelity detections presented in Table 2.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Measurement of the noise-bias-free CO autopower spectrum (solid black curve), $P_{CO,CO}(k_{CO(2-1)})$, in the ASPECS LP Band 3 survey. Error bars on $P_{CO,CO}(k_{CO(2-1)})$ represent values from a polynomial fit to the raw errors, $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle$, calculated from Equation (16). The inverse-variance-weighted mean CO power is plotted for two bins averaging modes in the upper and lower halves (open red diamonds) of probed $k_{CO(2-1)}$, as well as for a bin (filled red diamond) containing the inverse-variance-weighted mean CO power for all $k_{CO(2-1)} \sim 10-100 \text{ h Mpc}^{-1}$. The $3\sigma$ upper limit, calculated using the uncertainty on the latter binned power spectrum, is plotted as the red dotted-dashed line. For comparison, the gray swath bounds the $1\sigma$ confidence region for the noise-bias-free power spectrum of a single realization of a noise-only simulated data cube from the CASA task simobserve, with the corresponding error bars calculated in the same way as for the real data. Theoretical predictions for $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})$ from Popping et al. (2016; cyan dashed line) and a model based on SED fitting of known sources in the ASPECS survey field ("MAGPHYS model"; magenta dashed line) are also plotted. A dotted black line that illustrates where the measured power is zero is drawn for reference.</figcaption></figure>

## 4.2. Measurement of CO Autopower Spectrum at 0.001 ≲ z ≲ 4.5

The noise-bias-free autopower spectrum, $P_{CO,CO}(k_{CO(2-1)})$, is presented in Figure 4. We have averaged the power spectrum in linear bins of width $dk_{CO(2-1)} = 2\pi/r_{\perp,max} = 4.1 \text{ h Mpc}^{-1}$, measuring CO fluctuations on scales from $k_{CO(2-1)} \sim 10$ to $100 \text{ h Mpc}^{-1}$. Formally, the ASPECS survey volume provides access to 3D modes (i.e., modes containing both $k_{\perp}$ and $k_{\parallel}$ components) down to the fundamental mode $k_{CO(2-1)} = 4.1 \text{ h Mpc}^{-1}$ (see Table 1), though the number of independent modes $N_m (=196, or one mode per every channel in the cube)$ in this lowest wavenumber bin is small, and the resulting S/N on the power spectrum is low; it has been discarded in this analysis.

Errors on the power spectrum at each $k_{CO(2-1)}$ bin, $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle$, have been calculated using a 6° polynomial fit to the raw values calculated per Equation (16), as we would expect $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle$ to approach a smooth function as the number of realizations of the noise-only cubes—($\tau_1(k) - \tau_2(k)$), ($\tau_3(k) - \tau_4(k)$), etc.—approaches infinity.

As an independent check on our error estimation, we also compute the noise-bias-free power spectrum of noise-only simulated data cubes, $P_{N,N}(k_{CO(2-1)})$, created with the CASA task simobserve. The output of simobserve is the CASA measurement sets, $\tilde{\tau}_{1,N}$, $\tilde{\tau}_{2,N}$, $\tilde{\tau}_{3,N}$, and $\tilde{\tau}_{4,N}$, that have been generated to mock the ASPECS observational setup, including an identical mosaic pointing pattern and antenna configuration,

which determine the mosaic power pattern and synthesized beam sizes, respectively. We then produce dirty image data cubes with the same parameters (e.g., 40-channel rebinning in frequency) adopted for the real data and normalize the flux densities in each cube so that the rms of each frequency slice (or channel map) at a given $\nu_{\text{obs}}$ for a given simulated cube (e.g., $\tau_{1,N}$) is identical to the rms noise of the corresponding data cube (e.g., $\tau_1$) at the same $\nu_{\text{obs}}$ (Figure 3). In this way, we have constructed noise-only simulated image cubes, $\tau_{1,N}-\tau_{4,N}$, with noise properties similar to the real data cubes $\tau_1-\tau_4$. The resulting noise-bias-free power spectrum of this simulated noise-only data set is shown alongside $P_{CO,CO}(k_{CO(2-1)})$ in Figure 4.

To improve the S/N on $P_{CO,CO}(k_{CO(2-1)})$, we have averaged the power within individual wavenumber bins into two wider bins containing the first and second halves of the full $k_{CO(2-1)}$ range and a third set containing all $k_{CO(2-1)}$ bins in the available range. We then report the inverse-variance-weighted mean and corresponding inverse-variance-weighted error for the bin representing the power spectrum averaged across all $N_b = 23$ bins from $k_{CO(2-1)} = 9.55$ to $100.05 \text{ h Mpc}^{-1}$,

$$ \langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}} = -45 \pm 77 \text{ } \mu\text{K}^2 \text{ (Mpc h}^{-1})^3. $$

We compute similar quantities for the bin containing the lower half ($9.55 \text{ h Mpc}^{-1} \le k_{CO(2-1)} \le 54.98 \text{ h Mpc}^{-1}$) of the modes only, $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{low}}$, and the upper half ($59.20 \text{ h Mpc}^{-1} \le k_{CO(2-1)} \le 100.05 \text{ h Mpc}^{-1}$) of the modes
<table><caption>**Table 5**

Predictions for CO Shot-noise Power, $P_{CO(J-(J-1),CO(J-(J-1))}^{\text{shot}}(k_{CO(J-(J-1),z_{cen,CO(J-(J-1))}})[\mu\text{K}^2 (\text{Mpc h}^{-1})^3]$</caption><thead><tr><th>Line<br>(1)</th><th>Popping et al. (2016)<br>(2)</th><th>MAGPHYS Model<br>(3)</th><th>ASPECS LF<br>(4)</th><th>COLDz<sup>a</sup> LF<br>(5)</th><th>COPSS II<sup>b</sup><br>(6)</th></tr></thead><tbody><tr><td>CO(1–0)</td><td>0.55</td><td>0.04</td><td>0.040<sup>+0.31</sup><sub>-0.03</sub></td><td>...</td><td>...</td></tr><tr><td>CO(2–1)</td><td>8.15</td><td>89.23</td><td>80<sup>+71</sup><sub>-38</sub></td><td>170<sup>+509</sup><sub>-120</sub><sup>a</sup></td><td>1600 ± 700<sup>a</sup></td></tr><tr><td>CO(3–2)</td><td>17.82</td><td>8.61</td><td>130<sup>+320</sup><sub>-82</sub></td><td>...</td><td>...</td></tr><tr><td>CO(4–3)</td><td>8.77</td><td>6.33</td><td>27<sup>+90</sup><sub>-19</sub></td><td>...</td><td>...</td></tr><tr><td>Total<sup>b</sup></td><td>43.31</td><td>98.13</td><td>170<sup>+290</sup><sub>-92</sub></td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** (1) Line transition. (2) Popping et al. (2016). (3) Model based on SED fitting and the $L'_{CO}-L_{IR}$ relation (see text for details). (4) CO LFs from D19. The CO LFs refer to Schechter fits integrated from $L'_{min} = 10^8$ to $L'_{max} = 10^{12}$ K km s$^{-1}$ pc$^2$ except for $J = 1$, where no Schechter fit was performed. For CO(1–0), we have used tabulated LF data for each luminosity bin from $L'_{min} = 10^8$ to $L'_{max} = 10^9$ K km s$^{-1}$ pc$^2$. (5) Riechers et al. (2019). (6) Keating et al. (2016).

<sup>a</sup> Total area surveyed: ~60 arcmin$^2$.

<sup>b</sup> Total area surveyed: ~0.7 deg$^2$.

<sup>c</sup> CO(2–1) power inferred from CO(1–0) LF measurement. See text in Section 4.2.2 for details regarding the applied conversion.

<sup>d</sup> Sum of CO(1–0), CO(2–1), CO(3–2), and CO(4–3) shot-noise power, converted to the CO(2–1) frame.

only, $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{high}}$, finding, overall,

$$ \langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{low}} = -260 \pm 170 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 $$

$$ \langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{high}} = +10 \pm 86 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 $$

The measurements above are generally consistent with nondetections (i.e., $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle = 0$) at the quoted $1\sigma$ or $1.5\sigma$ (in the case of $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{low}}$) level and comparable to the noise-bias-free power spectrum measured for the simulated noise-only cubes:

$$ \langle P_{N,N}(k_{CO(2-1)}) \rangle_{\text{tot}} = +41 \pm 87 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 $$

$$ \langle P_{N,N}(k_{CO(2-1)}) \rangle_{\text{low}} = -47 \pm 180 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 $$

$$ \langle P_{N,N}(k_{CO(2-1)}) \rangle_{\text{high}} = +70 \pm 100 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 $$

Furthermore, the reported power spectra in each total, low, and high $k_{CO(2-1)}$ bin agree within ~$1\sigma$ of each other, suggesting that our measurement does not discern any spectral structure. Thus, we adopt $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}}$ as representative of the measured flat power spectrum and use it to place a $3\sigma$ upper limit on the noise-bias free CO power spectrum, $P_{CO,CO}(k_{CO(2-1)}) \le 190 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3 = \langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}} + 3 \times \langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}}$.

The analytically estimated shot-noise power based only on ASPECS LP blind detections (see Section 4.1.2), $[P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})]_{\text{det}} = 113.24 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3$, lies roughly a factor of 2 below our upper limit. (We refer here to $[P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})]_{\text{det}}$ calculated for the cropped region used in the power spectrum analysis, with CO flux values from GL19.)

We note that the effect of the CMB radiation is likely negligible on the cumulative measurement presented here. Assuming gas kinetic temperatures close to ~40 K, which are appropriate for the CO-emitting sources in ASPECS based on the output dust temperatures from MAGPHYS SED fits (Boogaard et al. 2019), corrections due to CMB on the CO line fluxes contributing to the emission are expected to be less than 25% at $z \le 4$, with the possible exception of CO(4–3) at $z \sim 4$ in low-density gas ($n_H = 10^{3.2} \text{ cm}^{-3}$) and non-LTE conditions (Sections 3.3 and 3.4 in da Cunha et al. 2013). In that case, the intrinsic CO(4–3) flux can be up to a factor of 2 higher than observed. However, since the product of distortion factors converting CO(4–3) power into the CO(2–1) frame is small, we do not expect this to have a significant effect on the

aggregate emission measured by the power spectrum, even if such conditions were representative of the ISM.

The two theoretical models we include in our comparison in Figure 4 (MAGPHYS and Popping) are also factors of ~2 to ~4.5 below our upper limit. The Popping model refers to the semi-analytic model described in Popping et al. (2016). The MAGPHYS model refers to a model CO catalog for a subset of ~1000 sources in the HUDF with (1) a (photometric, grism, or spectroscopic) redshift that is in the allowable ranged covered by CO transitions in ASPECS LP up to $J = 4$ and a reliable fit to the SED using MAGPHYS from the rest-frame UV/optical wavelengths out to the IR based on (2) a 1.6 $\mu$m flux density greater than 0.1 $\mu$Jy and (3) detections in at least five photometric bands. (For more details on the SED fitting of HUDF sources within the ASPECS field, see D19 and B19.) The CO(1–0) luminosities for the sources in this catalog have been estimated by scaling to the output IR luminosity from MAGPHYS according to Carilli & Walter (2013) and then using excitation corrections from Daddi et al. (2015) to predict higher $J$ CO line luminosities. We estimated the total CO shot-noise power for the model CO catalog per Equation (21), finding that—if we include all sources with stellar masses $M_* \ge 10^7 \ M_\odot$—$P_{CO,CO}^{\text{shot}}(k_{CO(2-1)}) = 98.13 \ \mu\text{K}^2 \ (\text{Mpc h}^{-1})^3$, which is consistent with the $3\sigma$ upper limit measured from the power spectrum and comparable to $[P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})]_{\text{det}}$.

Expectations for $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})$—separated into the constituent power from specific $J$ transitions, when possible—based on the above models, empirically derived CO LFs (Section 4.2.1), and higher-redshift CO(1–0) measurements (Section 4.2.2) have been gathered, for reference, in Table 5. Note that predictions based on semi-analytic simulations from Popping et al. (2016) are known to underestimate the number density of bright sources compared to observed CO LFs (D19) and thus may represent an unrealistically low prediction for $P_{CO,CO}(k_{CO(2-1)})$. (Please see the following section, particularly Equation (22), for details connecting the shot-noise power spectrum measurement to the LF.) The modeling framework in that paper has since been updated in Popping et al. (2019), now including, e.g., predictions using hydrodynamical simulations, but the authors found that the models again underpredict the bright end ($L'_{CO} > 10^{10} \text{ K km s}^{-1} \text{ pc}^2$) of the LFs by factors of 1–3 dex. (We direct the reader to Popping et al. 2019 for a discussion of the potential origins of the discrepancy between their models and observations.) Also, we point out that the
discrepancy between the Popping et al. (2016) predictions for CO(1–0) shot-noise power and the expected CO(1–0) shot noise from the MAGPHYS model and the ASPECS observations can be ascribed to the small volume coverage of the ASPECS survey volume at $z \sim 0.28$, which limits the number of bright sources in the MAGPHYS model catalog and precludes a Schechter function fit for CO(1–0) LF data in the D19 analysis. (See also the discussion on cosmic variance in Popping et al. 2019.) Specifically, in Table 5, the expected CO(1–0) shot-noise power from the ASPECS LF is calculated by integrating the LF only in luminosity bins where LF data are available, i.e., in the range from $L'_{CO(1-0)} = 10^8–10^9$ K km s$^{-1}$ pc$^2$, instead of the full range of $L'_{CO} (=10^8–10^{12}$ K km s$^{-1}$ pc$^2$) encapsulated by the Schechter function fits provided at higher $z$.

### 4.2.1. Comparison to CO LFs Derived from ASPECS LP

Integrating the LFs measured by ASPECS LP (with lower and upper limits of integration equal to $10^8$ and $10^{12}$ K km s$^{-1}$ pc$^2$) per Equation (22) yields the following estimates of shot-noise power for CO(2–1), CO(3–2), and CO(4–3), respectively: $80^{+71}_{-38}$, $130^{+320}_{-82}$, and $27^{+90}_{-19} \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$. In order to compare to our measurement of the total CO shot-noise power, we convert the individual shot-noise powers estimated for each transition into the CO(2–1) frame and sum, finding $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)}) = 170^{+290}_{-92} \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$. The $3\sigma$ upper limit on $P_{CO,CO}(k_{CO(2-1)}) \le 190 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$ determined via the power spectrum places a more stringent constraint on the total CO shot-noise power; it rules out a significant fraction of the allowable range, $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)}) = 78$ to $460 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$, obtained from the LF fits. The low sensitivity, however, on the power spectrum measurement precludes us from determining the amplitude of individual contributions from the different CO transitions to the aggregate value, so translating our limit on $P_{CO,CO}(k_{CO(2-1)})$ to constraints on individual LFs is highly speculative—as would be attempting to constrain the Schechter parameters of the individual CO LFs. We note, however, that the fixed faint-end slope $\alpha = -0.2$ in D19 implies that the shot-noise power spectrum is relatively insensitive to the low-luminosity systems in the CO LF. In this case, the shot-noise power is more sensitive to the Schechter parameters $\Phi_*$ and $L_*$, and our measurement suggests either lower normalizations $\Phi_*$ or a knee in the Schechter function that occurs at lower luminosities $L'_*$. This may be particularly applicable to the observed CO(3–2) and CO(4–3) LFs at $z > 2$, where the empirical constraints span a limited range in luminosity compared to the CO(2–1) LF at $z \sim 1$, and the resulting uncertainties on the Schechter parameters are large.

### 4.2.2. Comparison to Higher-redshift CO(1–0) Observations

In this section, we compare our measured $P_{CO,CO}(k_{CO(2-1)})$ to independent observational constraints on CO(1–0) shot-noise power at higher redshift.

First, we consider the CO(1–0) LFs determined at $z = 2–3$ by the VLA COLDz program (Riechers et al. 2019), which targeted survey volumes in COSMOS (~9 arcmin$^2$) and GOODS-N (~51 arcmin$^2$) across 8 GHz of bandwidth in the Ka band ($\nu_{\text{obs}} \approx 30–38$ GHz) with typical synthesized beam sizes ~3". We compute the probability distribution of $P_{CO(1-0),CO(1-0)}(k_{CO(1-0)}, z_{\text{cen}} = 2.4)$, using Equation (22) and the Schechter function parameter samples from the posterior distributions obtained with the approximate Bayesian computation

method for the merged COSMOS and GOODS-N data set (see Figure 6 in Riechers et al. 2019). We find that the distribution has a median $P_{CO(1-0),CO(1-0)}(k_{CO(1-0)}, z_{\text{cen}} = 2.4) = 276.58 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$, with a probable range of $P_{CO(1-0),CO(1-0)}(k_{CO(1-0)}, z_{\text{cen}} = 2.4) = 74.75 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$ (5th percentile) to $1119.35 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$ (95th percentile).

The COLDz constraints on the CO(1–0) shot-noise power spectrum at $z_{\text{cen,CO(1-0)}} = 2.4$ can be converted to a constraint on the CO(2–1) luminosity density at $z \sim 1$ (and, thus, the CO(2–1) shot-noise power spectrum at the same redshift) assuming that (1) the CO(2–1) line is thermalized and has the same brightness temperature as CO(1–0) at $z_{\text{cen}} = 2.4$ and (2) there is no evolution in the CO(2–1) luminosity density from $z_{\text{cen,CO(1-0)}} = 2.4$ to $z_{\text{cen,CO(2-1)}} = 1.3$. The first assumption is reasonable for the low $J$ transitions relevant here and is supported by observations of a variety of high-$z$ systems (see, e.g., Table 2 in Carilli & Walter 2013, which shows $L'_{CO(2-1)}/L'_{CO(1-0)} \sim 0.9$ for submillimeter galaxies, color-selected star-forming galaxies, etc.). The latter assumption is likely false in detail, but, given that the cosmic SFR density is relatively flat between $z = 1$ and 3, it is a reasonable first approximation; note also that the cosmic molecular gas densities at $z \sim 1$ and 3 are indistinguishable by current empirical standards. In any case, assumption (2) results in an underestimation of $P_{CO(2-1),CO(2-1)}$ if the CO(2–1) luminosity density at $z = 2.4$ is, in reality, lower than the luminosity density at $z = 1$, and vice versa. So, assuming that items (1) and (2) are valid, we have $P_{CO(1-0),CO(1-0)}(k_{CO(1-0)}, z_{\text{cen}} = 2.4) = P_{CO(2-1),CO(2-1)}(k_{CO(2-1)}, z_{\text{cen}} = 2.4)$, which implies a CO(2–1) luminosity density $\rho_{CO(2-1)}(z_{\text{cen}} = 2.4) = \rho_{CO(2-1)}(z_{\text{cen}} = 1.3) = 5.63 \times 10^{10} L_{\odot} \text{ Mpc}^{-3}$, yielding $P_{CO(2-1),CO(2-1)}^{\text{shot}}(k_{CO(2-1)}, z_{\text{cen}} = 1.3) = 167.03 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$; the 5th and 95th percentiles similarly give lower and upper bounds on $P_{CO(2-1),CO(2-1)}^{\text{shot}}(k_{CO(2-1)}, z_{\text{cen}} = 1.3) = 45.14$ and $675.98 \mu\text{K}^2$ (Mpc $h^{-1}$)$^3$, respectively. This range is consistent with both our upper limit on the measured $P_{CO,CO}(k_{CO(2-1)}, z_{\text{cen}} = 1.3)$ and the $P_{CO(2-1),CO(2-1)}^{\text{shot}}(k_{CO(2-1)}, z_{\text{cen}} = 1.3)$ estimated from the CO(2–1) LF fit derived from ASPECS LP data (see Section 4.2.1).$^{28}$ A direct comparison of the CO(1–0) LF inferred$^{29}$ by ASPECS at $z_{\text{cen}} = 2.6$ and that measured by COLDz at nearby redshifts also reveals an excellent agreement across the overlapping CO(1–0) luminosity range of $L'_{CO(1-0)} = (0.1–4) \times 10^{11}$ K km s$^{-1}$ pc$^2$, which suggests that our assumptions (1) and (2) are reasonable. Importantly, the apparent agreement at the bright end of the CO(1–0) LFs also suggests that the impact of cosmic variance on the shot-noise power spectrum presented in this study is modest, since the shot-noise measurement is inherently more sensitive to high-luminosity systems (because of the $\propto L^2$ dependence in Equation (22)).

We also compare the measured CO power spectrum of ASPECS LP data to the power spectrum measurement of CO(1–0) at wavenumbers $k_{CO(1-0)} \sim 1–10 \text{ h Mpc}^{-1}$ at $z \sim 2–3$.

<footnote>$^{28}$ Note that larger uncertainties on shot-noise power derived from COLDz are largely due to differences in fitting the measured LFs: the authors in Riechers et al. (2019) treated the faint-end slope of the LF as a free parameter, while Decarli et al. (2019) kept it fixed during their fitting.</footnote>

<footnote>$^{29}$ Decarli et al. (2019) used the CO(3–2) LF measured at $z_{\text{cen}} = 2.6$ to infer a CO(1–0) LF after accounting for CO excitation. The fiducial prescription to convert from CO(3–2) to CO(1–0) luminosities was based on Daddi et al. (2015), but the authors there explored alternative lower- and higher-excitation scenarios as well, finding that their results are qualitatively robust to the adopted CO line ratio.</footnote>
from COPSS II (Keating et al. 2016). COPSS II was a dedicated intensity mapping experiment carried out with the Sunyaev-Zel'dovich Array to observe noncontiguous fields (totaling a survey area of ~0.7 deg²) on the sky at ~2' spatial resolution. Keating et al. (2016) presented a marginal 2σ detection of $P_{CO(1-0),CO(1-0)}^{\text{shot}}(k_{CO(1-0)}) = 3000 \pm 1300 \ \mu\text{K}^2$ (Mpc h⁻¹)³ at $z_{\text{cen}} = 2.8$, which is at least a factor of 6 larger than the COLDz measurement at similar redshift. Following the same procedure as outlined above, we estimate the CO(2–1) shot noise at $z_{\text{cen}} = 1.3$ based on the COPSS II measurement, finding $1600 \pm 700 \ \mu\text{K}^2$ (Mpc h⁻¹)³, which is 4.8–12 times higher than our current 3σ upper limit, though the ASPECS measurement probes CO power at higher wavenumbers (>10 h Mpc⁻¹) compared to COPSS II.

### 4.3. CO-galaxy Statistics

While the noise-bias-free autopower spectrum (Section 4.2) provides an unbiased view of the aggregate CO-emitting galaxy population, the detection of the autopower spectrum is inherently challenging given that the measurement weights noisy voxels<fnref n="30" /> equally with voxels containing CO emission. Adjusting these weights by including information from another field that is correlated in spatial distribution with the target field of CO fluctuations could yield a higher-fidelity measurement of the CO properties (e.g., shot-noise power, mean surface brightness) of the secondary population represented in the additional field, as long as this additional data set does not contribute significantly to the noise.

Given that we expect the CO fluctuations to originate within the ISM of galaxies, we consider the available galaxy catalogs with spectroscopic redshifts (spec-z) in the HUDF as potential data sets to perform power spectrum analyses with ASPECS data. Spectroscopic redshifts are required in order to match the spectral precision of the ASPECS LP data cube, where astrophysical emission is expected to be contained within the channel width $\Delta\nu_{40chn}$, corresponding to a redshift resolution $\Delta z_{40chn}/(1+z_{\text{cen,CO(2-1)}}) = 0.0013$. In contrast, photometric redshifts measured via SED template fitting procedures, e.g., are typically characterized by uncertainties $\Delta z_{\text{phot}}/(1+z_{\text{phot}}) \sim 0.05$ (e.g., Coe et al. 2006; Brammer et al. 2008) and are unsuitable for use in the power spectrum analysis.

The MUSE Ultra Deep Field (UDF) survey (Bacon et al. 2017) has yielded ~1500 spectroscopic redshifts within a $3' \times 3'$ field in the HUDF (Inami et al. 2017), which includes the ASPECS LP areal footprint. The redshift range covered by the MUSE UDF extends beyond the ASPECS LP Band 3 spectral coverage, but there is overlap at $z < 1.5$ and $z > 3$ due primarily to the identification in MUSE spectra of rest-frame optical and UV emission features, respectively, such as [O II] and Lyα emission; from $z \sim 2$ to 3, MUSE spectroscopic redshifts are determined with the presence of various UV emission and/or absorption features.<fnref n="31" /> There are, in total, 680 sources identified by MUSE that are available in the ASPECS LP areal and redshift coverage, and 415 of these lie within the cropped region used here for the power spectrum analysis; unless otherwise noted, we include all spec-z reported

in the MUSE catalog.<fnref n="32" /> Of these 415, 24 MUSE sources fall within the redshift interval and volume probed by CO(1–0) in ASPECS LP and can be used to search for CO(1–0) emission in the ASPECS LP data; for CO(2–1), there are 128 MUSE sources in the relevant redshift range and area; for CO(3–2), 64 sources; and, for CO(4–3), 199 sources. Note that the relative number of MUSE sources available to correlate with the different CO J transitions is affected by the availability and intrinsic strength of the various aforementioned spectral features used to determine MUSE spec-z within each redshift interval.

Since galaxies in the MUSE spec-z catalog have been selected at shorter rest wavelengths than the millimeter-wave emission observed in ASPECS, it is not a priori known how closely they correlate with the CO fluctuations. However, based on the counterpart analysis of the 16 secure ASPECS LP CO detections, which all have optical/near-IR counterparts, we expect a nonzero positive correlation. In fact, as described in B19, all 10 CO(2–1) blind detections in ASPECS LP have a counterpart MUSE spec-z, with velocity offsets typically less than 100 km s⁻¹ between the CO and MUSE redshift (see Table 1 in B19); two of the five CO(3–2) blind detections in ASPECS LP have MUSE spec-z; and the single CO(4–3) blind detection does not have a MUSE spec-z.

In the following sections, we explore two statistics, namely, a masked autopower spectrum $P_{CO,CO}^{\text{gal}}(k_{CO(2-1)})$ (Section 4.3.1) and cross-shot-noise power spectrum $P_{CO,gal}^{\text{gal}}(k_{CO(2-1)})$ (Section 4.3.2), to assess the contributions to the observed shot-noise power and mean CO surface brightness from MUSE-selected galaxies, respectively.

#### 4.3.1. Masked Noise-bias-free Autopower Spectrum

We obtain $P_{CO,CO}^{\text{gal}}(k_{CO(2-1)})$ by adopting the same methods to measure the noise-bias-free autopower spectrum, $P_{CO,CO}^{\text{gal}}(k_{CO(2-1)})$, with one key difference: masking was performed on $\tau_1$, $\tau_2$, $\tau_3$, and $\tau_4$ to remove flux densities from voxels beyond a spatial radius of ~1" and spectral width of 0.165 GHz (or 1 channel) from the center of known positions of MUSE galaxies listed in Inami et al. (2017), prior to estimating the power spectrum. The ~1" radius was chosen so that any enclosed source flux would be encompassed by a full beamwidth; the true source flux is not recovered in the case of extended emission (see Table 2 for sources identified as extended and corresponding fluxes extracted per GL19). Explicitly, we are evaluating a statistic defined as

$$ \langle T_{CO}^{\text{gal}*}(\mathbf{k}) T_{CO}^{\text{gal}}(\mathbf{k}') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_{CO,CO}^{\text{gal}}(\mathbf{k}). \quad (23) $$

In the above equation, $T_{CO}^{\text{gal}}$ represents the masked ASPECS data cube, which has been weighted at every ith voxel according to

$$ (T_{CO}^{\text{gal}})_i = w_i (T_{CO})_i, \quad (24) $$

where $w_i = 1$ for voxels containing a MUSE galaxy, and $w_i = 0$ otherwise. Thus, in the same way that the noise-bias-free autopower spectrum measurement at $k = 10-100 \ h \ \text{Mpc}^{-1}$ yielded constraints on the second moments of the CO LFs (per

<footnote n="30">³⁰ Voxels refer to 3D resolution elements defined by the beam area and channel width.</footnote>

<footnote n="31">³¹ Please see Figure 1 of B19 for the redshift distribution—color-coded by the respective spectral feature(s) used for redshift determination—of MUSE galaxies within the ASPECS LP survey volume.</footnote>

<footnote n="32">³² As discussed in Inami et al. (2017), spec-z has been assigned confidences ranging from 1 to 3; redshifts with CONFID = 2 or 3 are considered “secure,” and those with CONFID = 1 have been determined as a “possible” redshift from the presence of a spectral line with uncertain identification.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Masked noise-bias-free autopower spectrum, $P_{CO,CO}^{MUSE}(k_{CO(2-1)})$, of ASPECS LP data cubes with MUSE 3D source positions from Inami et al. (2017). All 415 MUSE sources with spec-z that fall within redshift ranges observable by ASPECS in CO(1-0), CO(2-1), CO(3-2), and CO(4-3) are included. The solid black curve represents the total power measured using ASPECS LP data and all MUSE sources with potential CO emission up to $J=4$. The solid and dashed red curves show the power using all MUSE sources with and without, respectively, a previously detected CO counterpart from the ASPECS line search. The y-axis has been scaled by a factor $1/(P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{all})$, so that each curve represents the contribution of the respective subset of galaxies to the total power measured when using all MUSE galaxies.</figcaption></figure>

Equation (22)), the masked noise-bias-free autopower spectrum measurement at the same $k$ range constrains the second moments of the LFs of CO-emitting MUSE galaxies.

Errors on the masked noise-bias-free autopower spectrum, $\delta P_{CO,CO}^{gal}(k_{CO(2-1)})$, were obtained from Equation (16) and include an additional contribution estimated from 100 simulations of random MUSE source positions. This Poisson term was deemed necessary to prevent underestimating the error based only on Equation (16); after removing significant source flux from the cubes in the masking step, the pre-factor $\gamma = 0.25$ that appears in this equation might no longer accurately describe the relation of noise properties in $\tau_1$, $\tau_2$, $\tau_3$, and $\tau_4$ to those in $T_1$ and $T_2$.

Figures 5 and 6 show the measured masked noise-bias-free CO autopower spectrum, $P_{CO,CO}^{MUSE}(k_{CO(2-1)})$, between CO fluctuations in the ASPECS LP data cube and 3D positions from the MUSE spec-z catalog. In these figures, we compare the power measured when including (1) the MUSE positions of the CO blind detections (solid red curves), (2) all available MUSE positions (solid black curves), and (3) all available MUSE positions, excluding those corresponding to the positions of the CO blind detections (dashed red curves). To facilitate comparison of the relative contributions of each sample described by items (1)–(3), we have scaled the y-axis in each plot by a factor $1/(P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{all})$, where the denominator represents the total noise-bias-free cross-power using all MUSE positions and including potential CO emission from $J=1-4$ transitions. In the case where item (3) yields $P_{CO,CO}^{MUSE} = 0$, we can be confident that there is no measured “excess” power from MUSE galaxies with previously undetected CO emission. Note that, although our power spectrum measurements cannot probe fluctuations on scales of $k_{CO(2-1)} \gtrsim 100 \text{ h Mpc}^{-1}$, which correspond to the sub-beam size pixel gridding, we include these wavenumbers in the plot to show the effect of the beam size on the measured power,

which follows an exponential drop-off as expected; the sensitivity on the autopower spectrum measurement was insufficient to “detect” the beam roll-off, so we truncated the spectrum in Figure 4 at $k_{CO(2-1)} \approx 100 \text{ h Mpc}^{-1}$.

We find that the noise-bias-free cross-power spectrum is detected at high significance with an amplitude suggesting that the majority of CO emission contained within the survey volume is accounted for by the ASPECS blind detections and that the rest-frame optical/UV galaxies in the same field closely trace the observed CO emission (Figure 5). The latter result is unsurprising, given that 12/16 of CO blind detections had known MUSE counterparts. Only a small level of excess power—$\langle P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{\text{blind removed}} \rangle / \langle P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{\text{all}} \rangle = 7.05 \pm 3.2\%$, where we have averaged over wavenumber bins from $k_{CO(2-1)} = 9.55$ to $54.98 \text{ h Mpc}^{-1}$ to avoid the effects of the exponential beam roll-off$^{33}$ seen at higher $k_{CO(2-1)}$—is observed to originate from the positions of MUSE galaxies in the same field that do not have previously detected ASPECS counterparts, amounting to $3.8\%–10.7\%$ of the total masked CO autopower measured when including all MUSE positions with potential CO $J=1, 2, 3$, or 4 emission. In other words, the percentage of power recovered by the MUSE sources with previously detected ASPECS counterparts is 89\%–96\%.

The high S/N on $P_{CO,CO}^{MUSE}(k_{CO(2-1)})$ enables further decomposition of the total masked CO autopower spectrum into the individual $J$ transitions contributing to the aggregate signal (Figure 6). Here we identify emission from the CO(3-2) transition as the principal source of residual power after removing the ASPECS blind detections in the total masked CO autopower spectrum in Figure 6, while the power derived from all other $J$ transitions is zero. The level of residual power detected in the masked CO(3-2) power spectrum,$^{34}$ $P_{CO(3-2),CO(3-2)}^{MUSE}(k_{CO(2-1)}) = 1.8 \pm 0.56 \text{ µK}^2 \text{ (Mpc h}^{-1})^3$, represents $5.2\% \pm 1.6\%$ of the total masked CO autopower $\langle P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{\text{all}} \rangle$ (including all $J$ transitions), or $14\%–22\%$ of the total masked CO(3-2) power spectrum amplitude.

Next we consider the nature of the sources contributing to the excess in CO(3-2) emission relative to the expected power from the CO(3-2) blind detections only. The level of observed power could be attributed to a single bright source, e.g., with a CO(3-2) flux of the order of $\sim 0.10 \text{ Jy km s}^{-1}$, or multiple fainter sources. Since a CO(3-2) source flux of $\sim 0.10 \text{ Jy km s}^{-1}$ at $2 \le z \le 3$ implies a line luminosity of $4.2 \times 10^6 \text{ L}_\odot$ that is just below the mean sensitivity limit for the survey for that redshift range, it is possible that a single bright source with this flux would have been previously undetected in the ASPECS line search. The scenario where the excess power originates from CO(3-2) emitters below the individual detection threshold is also plausible, however, given that the MUSE catalog contains galaxies with lower stellar mass $M_*$ and SFRs than probed by the CO blind detections. Specifically, at $2 \le z \le 3$, MUSE sources probe down to $M_* \sim 10^9 \text{ M}_\odot$ and SFR $\sim 0.3 \text{ M}_\odot \text{ yr}^{-1}$, while the ASPECS-detected galaxies have $M_* \gtrsim 10^{10} \text{ M}_\odot$ and SFR $\gtrsim 10 \text{ M}_\odot \text{ yr}^{-1}$ (B19). We can test these scenarios by masking the ASPECS data cube down to progressively lower flux thresholds and

<footnote>$^{33}$ If we average over all $k_{CO(2-1)}$ from 9.55 to 100.05 h Mpc$^{-1}$, we find that the excess power is $\langle P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{\text{blind removed}} \rangle / \langle P_{CO,CO}^{MUSE}(k_{CO(2-1)})^{\text{all}} \rangle = 7.1 \pm 1.2\%$.</footnote>

<footnote>$^{34}$ Here $P_{CO(3-2),CO(3-2)}^{MUSE}(k_{CO(2-1)}) = 1.3 \pm 0.23 \text{ µK}^2 \text{ (Mpc h}^{-1})^3$, when averaging over all $k_{CO(2-1)}$ up to 100 h Mpc$^{-1}$, instead of the lower half of the modes, as adopted in the main text.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cross-power spectrum of ASPECS LP data cubes with MUSE 3D positions, separated into contributions by CO J transition. Color-coding is the same as in Figure 5. The magenta dashed curve in the bottom left panel denotes masked CO(3–2) autopower measured upon removing MUSE ID 24, in addition to MUSE IDs corresponding to ASPECS blind detections.</figcaption></figure>

determine the flux level where the excess power vanishes. Figure 7 shows the results of this analysis, where voxels with flux densities with $|F_ν| > 1.0, 0.75, 0.50, 0.25$, and $0.1$ mJy have been blanked so that the remaining emission in the data cubes is due to sources (or noise) fainter than the masking threshold; note that we consider the absolute value of the flux densities because we are working with dirty image cubes that contain $F_ν < 0$. From Figure 7, it is clear that roughly 75% of the excess power originates from voxels with flux densities greater than 0.25–0.50 mJy, which implies fluxes $F_ν \Delta\nu_{40\text{chn}} > 0.12–0.24$ Jy km s⁻¹, suggesting that a single relatively bright emitter is responsible for the majority of the observed power from galaxies with previously undetected CO sources. We note that this flux level is fainter than any of the previous blind CO(3–2) detections in Table 2 and comparable to a 0.17 Jy km s⁻¹ ASPECS CO(3–2) detection identified with a MUSE spec-z prior at z = 2.028; this source cannot be responsible for the excess power discussed here, however, as it lies outside the region of sky used in the power spectrum analysis. Nonzero power due to voxels with $F_ν < 0.25$ mJy implies that very faint sources with fluxes less than 0.10 Jy km s⁻¹ may exist in the data cube and contribute to the observed excess in the CO-galaxy power spectrum, but their overall contribution is small (i.e., <25%).

Since the masking analysis suggests that the observed excess CO(3–2) power could be due to a single source, we can try to identify this source by masking (i.e., setting to zero) one by one each MUSE source position; note that only one MUSE position

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Contribution of voxels with different flux density thresholds ($|F_ν| < 1.00, 0.75, 0.50, 0.25$, and $0.1$ mJy) to the masked noise-bias-free CO(3–2) autopower spectrum using MUSE source positions that lack a previous ASPECS blind CO(3–2) detection. The uppermost red dashed curve is identical to the red dashed curve in the bottom left panel in Figure 6; please refer to that figure for error bars on the measured $P_{\text{CO(3-2),gal}}(k_{\text{CO(2-1)}})$.</figcaption></figure>

is masked at a time. If the excess power is due to a single source, then the power will remain unchanged (within the measured uncertainties) until the MUSE position corresponding to the CO(3–2) emission is masked and the power goes to zero. Following this procedure, we observe the power drop
to zero (magenta curve in Figure 6) upon masking source MUSE ID 24. Masking of all other MUSE sources resulted in negligible changes to the power spectrum. Examination of the data cube $T_0$ reveals no significant flux at the source position corresponding to MUSE ID 24 (R.A. = 53°160088, decl. = -27.776356, $v_{\text{obs}}$ = 97.57 GHz). However, this location is within a beam's width of a known blind detected source, ASPECS LP-3mm.01 or MUSE ID 35 (at R.A. = 53°160587, decl. = -27.776120, $v_{\text{obs}}$ = 97.58), so the ~1" radius used to extract emission from MUSE positions when computing the masked autopower spectrum encompasses flux from the blindly detected source. Thus, we conclude that the observed excess in CO(3-2) power is not due to any MUSE sources with previously undetected CO(3-2) emission.

We thus recompute the masked autopower spectrum using all MUSE positions with potential CO $J = 1, 2, 3$, or 4 emission, now excluding MUSE ID 24, in order to revise the estimate of the fraction of total masked autopower recovered by MUSE positions corresponding to ASPECS blind detections. We find that the MUSE positions corresponding to ASPECS blind detections recover 106% ± 6.5% of the total CO shot-noise power.

### 4.3.2. Cross-shot-noise Power Spectrum

Additionally, we measure the cross-power spectrum between the ASPECS data and MUSE position field and refer to this quantity as the cross-shot-noise power spectrum, $P_{\text{CO,gal}}(k_{\text{CO}(2-1)})$. In this case, we work directly with the ASPECS data cube $T_0$, which has the lowest rms compared to the subsets $T_1, T_{II}$, etc. Since the noise in $T_0$ has not been quantified, the error on $P_{\text{CO,gal}}(k_{\text{CO}(2-1)})$ is derived via simulation of random MUSE positions only and is a lower estimate of the true error.

We normalize the grid of MUSE galaxy positions to have units corresponding to the dimensionless density fluctuation field $\delta_{\text{gal}}(\mathbf{x}_i) = (n_{\text{gal}}(\mathbf{x}_i) - \langle n_{\text{gal}} \rangle) / \langle n_{\text{gal}} \rangle$, where $n_{\text{gal}}(\mathbf{x}_i)$ refers to the number density of galaxies at position $\mathbf{x}_i$ in the cube, and $\langle n_{\text{gal}} \rangle$ is the mean number density of galaxies in the full volume. The cross-shot-noise power spectrum between $T_0$ and the dimensionless density fluctuation cube, $G$, is then

$$ \langle T_0^*(\mathbf{k}) G(\mathbf{k}') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_{\text{CO,gal}}(\mathbf{k}), \quad (25) $$

with units of $\mu\text{K} (\text{Mpc h}^{-1})^3$.

As derived in, e.g., Wolz et al. (2017) and Breysse & Alexandroff (2019), $P_{\text{CO,gal}}(k)$ is proportional to the mean CO surface brightness of CO-emitting MUSE galaxies, $\langle T_{\text{CO,gal}} \rangle$,

$$ P_{\text{CO,gal}}(k) = \frac{\langle T_{\text{CO,gal}} \rangle}{\langle n_{\text{gal}} \rangle}. \quad (26) $$

Since the factor $1/\langle n_{\text{gal}} \rangle$ is equal to the amplitude of the shot-noise term of the galaxy autopower spectrum, $P_{\text{gal,gal}}^{\text{shot}}$ (units of $(\text{Mpc h}^{-1})^3$), we rearrange Equation (26) to write

$$ \langle T_{\text{CO,gal}} \rangle = \frac{P_{\text{CO,gal}}(k)}{P_{\text{gal,gal}}^{\text{shot}}}, \quad (27) $$

which has units of $\mu\text{K}$.

Constraints on $\langle T_{\text{CO,gal}} \rangle$ are obtained from the cross-shot-noise power spectrum, $P_{\text{CO,gal}}(k_{\text{CO}(2-1)})$, according to Equation (27). As illustrated in Figure 8, we strongly detect the cross-shot-noise power spectrum in both considered $J = 3$ and

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mean CO(2–1) (top panel) and CO(3–2) (bottom panel) surface brightness of MUSE galaxies. Color-coding is the same as in Figure 5. Values on the y-axis in each panel have been scaled by the total mean CO(2–1) and CO(3–2) surface brightnesses, $\langle T_{\text{CO(2-1),MUSE}} \rangle_{\text{all}}$ and $\langle T_{\text{CO(3-2),MUSE}} \rangle_{\text{all}}$, so that curves plotted represent the relative contribution of MUSE galaxies with (red solid curve) and without (red and orange dashed curves) ASPECS blind detections. In the top panel for CO(2–1), the orange curve represents MUSE galaxies with potential CO(2–1) emission, excluding spectra with low confidence (CONFID = 1).</figcaption></figure>

$J = 2$ transitions, finding no significant contribution to the observed mean CO surface brightness from MUSE sources with previously undetected CO(2–1) or CO(3–2) emission. It is interesting to note, however, that the average value of $\langle T_{\text{CO(2-1),gal}} \rangle$ is slightly nonzero after removing the blindly detected sources from the sample (red dashed curve, top panel). After we discard the 18 MUSE positions with potential CO(2–1) emission that have poorly characterized MUSE spectra (or CONFID = 1, i.e., a nonsecure redshift based on a singly unidentified line), then there is marginal (2.5$\sigma$) detection of excess CO(2–1) surface brightness from MUSE emitters without previous ASPECS CO detections (orange dashed curve) that can contribute $0.07 \pm 0.02 \ \mu\text{K}$, or $19\% \pm 7.8\%$ of the total observed CO(2–1) emission.<fnref n="35" /> Not shown in Figure 8, we measure nondetections of CO(1–0) and CO(4–3) surface brightnesses ($-0.0020 \pm 0.0082$ and $-0.022 \pm 0.028 \ \mu\text{K}$, respectively) in MUSE galaxies.

**Implications for CO LFs.** If there is a one-to-one correlation between CO-emitting galaxies and MUSE-selected galaxies, then $\langle T_{\text{CO,gal}} \rangle = \langle T_{\text{CO}} \rangle$ (and $P_{\text{CO,gal}}^{\text{NBF}}(k_{\text{CO}(2-1)}) = P_{\text{CO,CO}}(k_{\text{CO}(2-1)})$), and one can use the above measurements on $\langle T_{\text{CO,gal}} \rangle$ (and $P_{\text{CO,gal}}^{\text{NBF}}(k_{\text{CO}(2-1)})$) to constrain the CO LF via Equations (17) and (22). For example, if the observed ~20% excess surface brightness in CO(2–1) from galaxies without ASPECS detections is real, and the MUSE galaxies represent the complete population of total CO(2–1) emitters, then one can deduce that the ASPECS

<footnote n="35">³⁵ When averaging all bins from $k = 10$ to $100 \ h \ \text{Mpc}^{-1}$, we find that the significance of this “excess” improves, and the previously undetected sources contribute 23% ± 5.8% of the total observed CO(2–1) emission.</footnote>
survey recovers ≥80% of the CO(2–1) surface brightness at its sensitivity threshold. Keeping other Schechter parameters in Table 3 fixed, this suggests a relatively flat faint-end slope $\alpha \lesssim -0.1$ for the CO(2–1) LF at $z \sim 1$.36 Of course, if there are CO(2–1) emitters that do not have MUSE spec-z, then $\langle T_{CO,gal} \rangle \neq \langle T_{CO} \rangle$, and we cannot reliably infer constraints on the CO(2–1) LF. Given the high percentage (100%) of ASPECS CO(2–1) detections with MUSE spec-z counterparts, it is possible that $\langle T_{CO(2-1),gal} \rangle$ is not dramatically different from $\langle T_{CO(2-1)} \rangle$. For CO(3–2), however, the percentage is much lower (40%), so we do not attempt to draw conclusions about the CO(3–2) LF based on $\langle T_{CO(3-2),gal} \rangle$.

## 5. Beyond ASPECS: Detecting the CO Power Spectrum

Based on the analytic estimate $[P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})]_{\text{det}} = 113.24 \ \mu\text{K}^2 (\text{Mpc h}^{-1})^3$ from blindly detected sources (Section 4.1.2) and our $3\sigma$ upper limit $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)}) \lesssim 187.29 \ \mu\text{K}^2 (\text{Mpc h}^{-1})^3$, we consider the following question: what would be needed in an ALMA line survey to obtain a significant ($\ge 5\sigma$) detection on the CO shot-noise power spectrum?

The empirically determined uncertainty (Equation (16)) on the measured shot-noise power $\langle P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}} = -45 \ \mu\text{K}^2 (\text{Mpc h}^{-1})^3$ was found to be $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle_{\text{tot}} = 77 \ \mu\text{K}^2 (\text{Mpc h}^{-1})^3$. If the ASPECS blind detections can account for the bulk of the observed power, i.e., if $[P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})]_{\text{det}}$ closely approximates $P_{CO,CO}(k_{CO(2-1)})$ in our measured $k$ range, which is supported by the cross-power spectrum analysis in Section 4.3.1, then $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle$ must be reduced by a factor of 3.4 in order to obtain a $5\sigma$ detection on the total CO shot-noise power, or $\sigma_{N,40chn}$ must be reduced by a factor of $\sqrt{3.4} = 1.8$, since $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle \propto \sigma_{N,40chn}^2$. Since $\langle \delta P_{CO,CO}(k_{CO(2-1)}) \rangle \propto \sigma_{N,40chn}^2$, this can be achieved by increasing the integration time $t_{\text{int}}$ by a factor of 3.4, as $\sigma_{N,40chn} \propto t_{\text{int}}^{-1/2}$.

Alternatively, the S/N on $P_{CO,CO}^{\text{shot}}(k_{CO(2-1)})$ can be improved by increasing the number of independent $k$ modes in the survey. Because $S/N = P_N/(N_m)^{1/2}$, the total number of modes would need to be increased by a factor of (3.4)$^2 = 12$. In order to cover the same redshift range in CO, this would require enlarging the survey area to $12\Delta\theta_S^2 = 12(1.8)^2 = 39 \ \text{arcmin}^2$ while scanning in frequency across the same 30 GHz bandwidth, which is substantially more expensive than the ~threefold increase in $t_{\text{int}}$ estimated above. However, increased areal coverage would enable observations of more massive galaxies that are not captured in ASPECS LP.

The community is considering the Next Generation VLA (ngVLA) project as an order-of-magnitude improvement in observational capabilities in the 1–115 GHz regime over existing facilities, such as ALMA and the JVLA (Murphy et al. 2018; Selina et al. 2018). The ngVLA core array will have 100 × 18 m diameter antennas within an ~1 km diameter and a minimum bandwidth of 20 GHz. While the primary beam, and thus the instantaneous field of view, is a factor of 2.25 smaller than that of ALMA at a given frequency, the collecting area is about 4.5 times larger, and the bandwidth is at least 2.5 times larger. The implied time to cover the same cosmic

volume to the same sensitivity is then a factor of about 20 shorter than for the current ASPECS program.37

### 5.1. Power Spectrum versus Individual Galaxy Detection

Given that the ASPECS blind detections recover ~100% of the observed noise-bias-free cross-power between CO and MUSE galaxies, it is difficult to motivate longer integration times to obtain a detection on the power spectrum. To explain why the ASPECS LP 3mm survey is more efficient at recovering the CO shot-noise power by detecting individual galaxies, we consider the relationship between the significance of the individual detections and the expected uncertainty on the shot-noise power spectrum.

Suppose the survey can detect galaxies down to some minimum luminosity $L_{\text{min,det}}$. Explicitly, the relation between the mean surface brightness sensitivity $\sigma_N$ (in units of Jy sr$^{-1}$, for example) and $L_{\text{min,det}}$ is

$$ \frac{L_{\text{min,det}}}{4\pi D_L^2 \Delta\nu_{\text{chn}} \Delta\theta_b^2} = \epsilon\sigma_N, \qquad (28) $$

where $\Delta\nu_{\text{chn}}$ and $\Delta\theta_b^2$ are in units of Hertz and steradian, respectively, and $\epsilon$ is the required significance for a detection. For example, in Section 4.1.1, we set the ASPECS blind detection threshold as $\epsilon = 7$. Then, we can rewrite $P_N = \sigma_N^2 V_{\text{vox}}$ (Equation (12)) in terms of $L_{\text{min,det}}$ and $\epsilon$,

$$ \begin{aligned} P_N &= \frac{1}{\epsilon^2} \left( \frac{L_{\text{min,det}}}{4\pi D_L^2 \Delta\nu_{\text{chn}} \Delta\theta_b^2} \right)^2 y D_A^2 \Delta\nu_{\text{chn}} \Delta\theta_b^2 \\ &= \frac{A^2 L_{\text{min,det}}^2}{\epsilon^2 V_{\text{vox}}}, \end{aligned} \qquad (29) $$

where $A = yD_A^2 = cD_L^2/(H(z)\nu_{\text{rest}})$.

For the simple toy case, where all galaxies have identical luminosity $L = L_{\text{min,det}}$ and number density $n$, the shot-noise power (see Equation (22)) can be written as $P_{\text{shot}} = A^2 L_{\text{min,det}}^2 n$ and the S/N on $P_{\text{shot}}$ as

$$ \begin{aligned} S/N &= \frac{P_{\text{shot}}}{P_N} \sqrt{N_m} \\ &= \epsilon^2 \sqrt{N_m} n V_{\text{vox}} \\ &= \frac{\epsilon^2 N_{\text{gal}}}{\sqrt{N_m}}. \end{aligned} \qquad (30) $$

Note that we have used $nV_{\text{vox}} = nV_{\text{survey}} \times (V_{\text{vox}}/V_{\text{survey}}) = N_{\text{gal}}/N_m$, where $N_{\text{gal}}$ is the total number of galaxies in the survey, to arrive at Equation (30).

For the more realistic case, where all galaxies follow a Schechter-form LF, $P_{\text{shot}} = A^2 L_*^2 \Phi_* \Gamma(3 + \alpha)$ (e.g., Lidz & Taylor 2016), and we obtain

$$ \begin{aligned} S/N &= \frac{\epsilon^2 L_*^2 \Gamma(3 + \alpha)}{\sqrt{N_m} L_{\text{min,det}}^2} \Phi_* V_{\text{survey}} \\ &= \frac{\epsilon^2 N_{\text{gal}}}{\sqrt{N_m}}, \end{aligned} \qquad (31) $$

<footnote n="37">37 Note that we have assumed an observed wavelength range centered at 3 mm, identical to ASPECS LP, for this comparison, and ngVLA will also offer powerful centimeter-wave capabilities.</footnote>

<footnote n="36">36 This constraint on $\alpha$ is robust to changes in $\Phi_*$ and $L'_*$ within the quoted uncertainties in Table 3.</footnote>
where $N_{\text{gal}} = L_*^2 \Gamma(3 + \alpha) / L_{\text{min,det}}^2 \Phi_* V_{\text{survey}}$ is the effective number of detectable galaxies within the survey volume.

For the ASPECS survey and relevant Schechter function LFs, $\epsilon = 7$ and $N_{\text{gal}} = 15$ CO(2–1) + 7 CO(3–2) + 5 CO(4–3) = 27$ galaxies. Here $N_m$ is $1.18 \times 10^7$ modes, or the fraction of independent modes that fall within $k = 10–100$ h Mpc$^{-1}$. Note that, for mode counting, we refer to the $T_{0,2chn}$ cube defined by a grid of roughly 64 beams $\times$ 64 beams $\times$ 3935 channels, because this is the cube where the threshold for individual detection was defined. Then, S/N = 0.39, and the CO power spectrum is not expected to be detected at high significance;<fnref n="38" /> the shot-noise signal from individual galaxies is “diluted” by the factor $1/\sqrt{N_m}$. In this regime of small $\sigma_N$ and survey volume, the “traditional galaxy survey analysis,” where one identifies individual sources above the flux limit, provides a better S/N on the shot-noise term than the power spectrum analysis; the number and formulae here demonstrate this.

## 6. Conclusions

We presented a power spectrum analysis of the ASPECS LP 3 mm data set. Key results from this analysis are as follows.

(i) We derive a lower limit on the mean total CO surface brightness $\langle T_{\text{CO}} \rangle$ at 99 GHz by summing the observed fluxes from the 16 ASPECS blind detections: $\langle T_{\text{CO}} \rangle \ge 0.55 \pm 0.02$ µK. As the CO clustering power is proportional to $\langle T_{\text{CO}} \rangle^2$, this information from individually detected sources provides information on the CO clustering amplitude at large physical scales ($k < 1.0$ h Mpc$^{-1}$), as well as an indication of foreground contamination in CMB spectral distortion mapping experiments.

(ii) We derive an upper limit ($3\sigma$) on the 3D CO autopower spectrum at $10$ h Mpc$^{-1} \lesssim k \lesssim 100$ h Mpc$^{-1}$, $P_{\text{CO,CO}}(k_{\text{CO}(2-1)}) \lesssim 187.29$ µK$^2$ (Mpc h$^{-1}$)$^3$, which is broadly consistent with the ASPECS-observed CO LFs presented in D19. The upper limit measured in this study places a significantly tighter constraint on total CO power than predictions based on the Schechter-form LFs, due to the large uncertainties in parameters $\Phi_*$ and $L_*$ for the individual LFs at $z > 2$. Extrapolating the observed CO(1–0) LF at $z = 2–3$ in the COLDz survey to a CO(2–1) power in the ASPECS redshift range is consistent with our result, while the same procedure applied to the COPSS II measured shot-noise power yields a CO(2–1) power that is >5 times greater than our upper limit.

(iii) We report that detections of the masked noise-bias-free autopower spectrum, $P_{\text{CO,CO}}^{\text{gal}}(k_{\text{CO}(2-1)})$, and cross-shot-noise power spectrum, $P_{\text{CO,gal}}(k_{\text{CO}(2-1)})$, between CO and rest-frame optical/UV tracers from the MUSE spectroscopic redshift catalog yield constraints on the second and first moments, respectively, of the CO-emitting MUSE galaxy LFs. We found that $106\% \pm 6.5\%$ of the measured power in total CO shot noise is comprised of surface brightness fluctuations from MUSE galaxies

with previously detected CO from the ASPECS blind line search. We also constrained the contribution of ASPECS blind detections to the observed CO mean intensity of MUSE emitters, finding that up to ~20% of $\langle T_{\text{CO(2-1),MUSE}} \rangle$ is attributed to emission from galaxies below the threshold for individual detection. With the assumption that all observed CO(2–1) emission originates from MUSE emitters, $P_{\text{CO,gal}}(k_{\text{CO(2-1)}})$ can be used to place a direct constraint on the mean CO(2–1) surface brightness and, thus, faint-end slope of the LF, suggesting $\alpha \le -0.1$ at $z \sim 1$.

B.D.U. would like to thank Riccardo Pavesi for providing data to aid in the comparison of the ASPECS shot-noise power spectrum measurement with COLDz luminosity functions and Guochao Sun for providing model data to use in Figure 2. We thank Patrick Breysse for useful discussions pertaining to the interpretation of the CO-galaxy cross-power spectrum. F.W. acknowledges support from ERC grant “Cosmic Gas” (740246). D.R. acknowledges support from the National Science Foundation under grant No. AST-1614213 and the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. The authors jointly thank the anonymous referee for a constructive report that helped to improve this paper.

**Facilities:** ALMA. **data:** 2016.1.00324.L. ALMA is a partnership of the ESO (representing its member states), NSF (USA), and NINS (Japan), together with the NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, AUI/NRAO, and NAOJ.

## ORCID iDs

Bade D. Uzgil https://orcid.org/0000-0001-8526-3464

Chris Carilli https://orcid.org/0000-0001-6647-3861

Adam Lidz https://orcid.org/0000-0002-3950-9598

Fabian Walter https://orcid.org/0000-0003-4793-7880

Nithyanandan Thyagarajan https://orcid.org/0000-0003-1602-7868

Roberto Decarli https://orcid.org/0000-0002-2662-8803

Manuel Aravena https://orcid.org/0000-0002-6290-3198

Frank Bertoldi https://orcid.org/0000-0002-1707-1775

Paulo C. Cortes https://orcid.org/0000-0002-3583-780X

Gergő Popping https://orcid.org/0000-0003-1151-4659

Dominik A. Riechers https://orcid.org/0000-0001-9585-1462

Axel Weiss https://orcid.org/0000-0003-4678-3939

## References

Aravena, M., Decarli, R., González-López, J., et al. 2019, ApJ, 882, 136

Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1

Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729

Boogaard, L. A., Decarli, R., González-López, J., et al. 2019, ApJ, 882, 140

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Breysse, P. C., & Alexandroff, R. M. 2019, MNRAS, 490, 260

Carilli, C. L., Chluba, J., Decarli, R., et al. 2016, ApJ, 833, 73

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Natur, 466, 463

Coe, D., Benítez, N., Sánchez, S. F., et al. 2006, AJ, 132, 926

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46

Decarli, R., Walter, F., González-López, J., et al. 2019, ApJ, 882, 138

<footnote n="38">³⁸ Note that $N_{\text{gal}}$ and, thus, the derived S/N contain contributions from different J transitions. If we project the S/N into the CO(2–1) frame, as we have done for the power spectrum analysis, we must scale the relative contributions to the total S/N by the appropriate distortion factors, resulting in S/N = 0.29.</footnote>
Dillon, J. S., Liu, A., Williams, C. L., et al. 2014, PhRvD, 89, 023002

González-López, J., Decarli, R., Pavesi, R., et al. 2019, ApJ, 882, 139

Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142

Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6

Inami, H., Bacon, R., Brinchmann, J., et al. 2017, A&A, 608, A2

Keating, G. K., Marrone, D. P., Bower, G. C., et al. 2016, ApJ, 830, 34

Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2013, ApJS, 209, 3

Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, arXiv:1709.09066

Li, T. Y., Wechsler, R. H., Devaraj, K., & Church, S. E. 2016, ApJ, 817, 169

Lidz, A., & Taylor, J. 2016, ApJ, 825, 143

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415.

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Murphy, E. J., Bolatto, A., Chatterjee, S., et al. 2018, in ASP Conf. Ser. 517, Science with a Next Generation Very Large Array, ed. E. Murphy (San Francisco, CA: ASP), 3

Neeleman, M., Prochaska, J. X., Ribaudo, J., et al. 2016, ApJ, 818, 113

Pavesi, R., Sharon, C. E., Riechers, D. A., et al. 2018, ApJ, 864, 49

Popping, G., Pillepich, A., Somerville, R. S., et al. 2019, ApJ, 882, 137

Popping, G., van Kampen, E., Decarli, R., et al. 2016, MNRAS, 461, 93

Pullen, A. R., Chang, T.-C., Doré, O., & Lidz, A. 2013, ApJ, 768, 15

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7

Selina, R. J., Murphy, E. J., McKinnon, M., et al. 2018, in ASP Conf. Ser. 517, Science with a Next Generation Very Large Array, ed. E. Murphy (San Francisco, CA: ASP), 15

Sun, G., Moncelsi, L., Viero, M. P., et al. 2018, ApJ, 856, 107

Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, MNRAS, 434, L46

Taylor, G. B., Carilli, C. L., & Perley, R. A. (ed.) 1999, ASP Conf. Ser. 180, Synthesis Imaging in Radio Astronomy II, A Collection of Lectures from the Sixth NRAO/NMIMT Synthesis Imaging Summer School (San Francisco, CA: ASP)

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Walter, F., Decarli, R., Sargent, M., et al. 2014, ApJ, 782, 79

Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, ARA&A, 43, 861

Wolz, L., Blake, C., & Wyithe, J. S. B. 2017, MNRAS, 470, 3220