<figure><img src="image_28.png" /><figcaption>Figure 17. Estimated purity as a function of candidate S/N. The blue points indicate the higher-significance candidates reported in this work. The purity estimation utilizes an extended candidate list down to an S/N of 5. The red shaded area indicates the S/N range that we consider for calculating purities. We also show the more conservative prescription adopted in applying purity corrections, i.e., assuming a uniform distribution for the purity, treating the calculated value as an upper limit. Our alternative prescription instead uses the calculated purities as having a Gaussian uncertainty of 100% (the upper $1\sigma$ limit is indicated by the dashed line). We do not show the highest-S/N line detections corresponding to AzTEC-3 and GN10 because they correspond to CO(2–1) line emission, and they were not included in the purity estimation because of their high S/N (see the text for details).</figcaption></figure>

to be real (Figure 17). In the GOODS-N field, there is no gap in the S/N distribution of the line candidates (the highest-S/N source is GN10 at $z > 5$). Therefore, we include all candidates in the range $5 < S/N < 6.4$.

The procedure we have described would attribute a purity of 70% for the candidate COLDz.GN.1, 50% to GN19 (which we manually correct to be 100% because we know it to be a real line), and in the 30%–50% range for the other S/N ~ 6 candidates, subsequently decreasing to about 7% at S/N = 5.5 (Figure 17).

When utilizing these purities to assemble the CO luminosity function, we consider two possible alternative strategies that allow us to estimate the effects of the systematic uncertainties introduced by our purity computation. In the first approach, we treat these purities as having 100% uncertainty; i.e., we will draw purities (for the Monte Carlo sampling used to estimate the allowed range of the luminosity function) as independent random numbers normally distributed around the estimated

values, with standard deviation equal to the purity estimate themselves and truncating at zero. The alternative approach is to implement these purities as upper limits and to draw purities from a uniform distribution between zero and the calculated values. The latter provides a more conservative purity estimation. Therefore, the luminosity function constraints are somewhat lower in this method, although they are compatible between the two methods. This conservative approach attempts to implement the additional information coming from the lack of clear multiwavelength counterparts to our moderate-S/N candidates. We will present the detailed results of both approaches in Paper II.

The S/N thresholds adopted in Section 4 and Table 3 correspond to approximate purities of ~4% and ~7% for COSMOS and GOODS-N, respectively. We emphasize that previously employed definitions of purity have differed significantly. In particular, we attempt to assess the fidelity defined by Walter et al. (2016) for our candidate selection. The comparison is not straightforward because the definition of fidelity used in that work relies on the details of their line search algorithm, but an approximate implementation of their method indicates an equivalent fidelity of approximately 80%–90% for COSMOS and 50%–60% for GOODS-N in their method.

## F.2. Estimating Noise Tail Extent from Data Cube Sizes

Due to the short-scale noise correlation intrinsic to interferometric noise (over the synthesized beam length scale), the calculation of the highest expected S/N due to noise (both positive and negative) is not straightforward, as the counting of “independent elements” is nontrivial. Vio & Andreani (2016) and Vio et al. (2017) independently discussed a similar analysis of this case. We have reached the same conclusions, although we take a slightly different approach, as we describe below. A detailed analysis of extreme value statistics in the case of smooth Gaussian random fields (which is a good approximation for interferometric noise) was developed by Bardeen et al. (1986) and Bond & Efstathiou (1987), among others, and expanded upon by Colombi et al. (2011). Here we only summarize the main results as relevant to our data and discuss the implications. The objective is a description of the probability distribution function for the highest S/N in a data cube that is uniquely due to noise and how this varies as a function of cube “size.” If we consider the original data cube, then noise is not correlated across different channels, and a noise realization is equivalent to a 2D case with spatial correlation only and an area equivalent to the total area across the full cube (i.e., the sum of the areas over the independent channels). In this case, the approximate cumulative distribution function for the highest S/N ($\nu$) to be expected from such a noise realization is given by

$$P(\nu_{\max} < \nu) \approx \exp\left(-\frac{1}{4\sqrt{2\pi}} N_n \nu e^{-\nu^2/2}\right), \quad (15)$$

where $N_n$ is the “naive” counting given by the total area divided by the “beam area” (defined by a radius equal to the beam standard deviation). The second case regards the case where correlation of the noise across channels has been introduced (for example, by convolution with a spectral template in order to match-filter) and is also relevant to line searches in the form of the noise properties of matched-filtered cubes. In this case, the correlation takes place in 3D, and a slightly different
approximate formula describes the cumulative distribution function for the highest S/N ($\nu$) to be expected,

$$P(\nu_{\text{max}} < \nu) \approx \exp\left(-\frac{1}{6\pi\sqrt{2}}N_n\nu^2e^{-\nu^2/2}\right), \quad (16)$$

where $N_n$ is the “naive” counting given by the total cube volume divided by the “effective beam volume” (an ellipsoid with a radius equal to the beam standard deviation in the spatial dimension and the standard deviation of the template used in the spectral dimension). We have verified that the highest-significance noise peaks measured as negative features in our data are compatible with these probabilistic predictions. We note that these distribution functions are quite broad and only predict the highest S/N expected due to noise to be approximately in the S/N = 5.5–6 range for our deeper mosaic and S/N = 5.7–6.4 for our wider mosaic. This is a manifestation of the strong intrinsic stochasticity of the noise tails. We also note that the “effective number of independent elements” implied by these estimates is ~10 and ~20 times higher than the naive counting in the 2D and 3D cases, respectively, and that these ratios are themselves increasing functions of data-cube size. The naive counting of independent elements would therefore lead to a significant underestimation of the extent of the noise tails. This conclusion is compatible with the results of Vio & Andreani (2016) and Vio et al. (2017). However, we here report equations that explicitly describe the distribution of the maximum S/N to be expected from noise rather than implicitly through the probability distribution function of local maxima. The analysis above is only approximately equivalent to the analysis presented in Vio et al. (2017), because they express the distribution function of interest as a function of $N_p$, i.e., the number of local maxima in the noise realization, which is itself a random variable with its own probability distribution.

### F.3. Artificial Source Analysis

In order to estimate the completeness and biases introduced by our line search and flux extraction methods, we perform an extensive probabilistic analysis of artificially injected sources into our maps. The main goal of this analysis is to establish a probabilistic connection between recovered candidate properties and intrinsic properties such as spatial size, velocity width, and line flux. This will provide some control over the uncertainties that affect the analysis of the CO luminosity function (Paper II). We also develop a method to correct the luminosity function by the completeness of our line search, which avoids a purely “per-source” completeness estimation as far as possible (due to the bias of “per-source” corrections) while avoiding assumptions that would significantly affect the result.

Since the large majority of the data cube contains very little signal, we use the data themselves as our model for the noise and inject artificial sources of varying size, velocity width, and flux at random positions in the data cube (Table 8). We inject sources in each cube (500 in COSMOS and 2500 in GOODS-N), estimating that this will not cause crowding of the field—therefore not causing overlaps between different sources—during the line search and effectively simulating the recovery of each injected source individually. We then analyze each injected cube following the same steps of MF3D that we applied to the real data, and in the end, we search for the injected sources to determine the recovered S/N and the line parameters that would have been

<table><caption>**Table 8**
Artificial Source Injected Sizes</caption><thead><tr><th></th><th>COSMOS</th><th>GOODS-N</th></tr></thead><tbody><tr><td>Intrinsic spatial size (arcsec)</td><td>~0.5, ~3.0, ~4.7</td><td>1.0, 2.5, 4.5</td></tr><tr><td>Convolved spatial size (arcsec)</td><td>2.6, 4.0, 5.3</td><td>2–3, 3–4, 4.8–5.4</td></tr><tr><td>Frequency width (MHz)</td><td>23.2, 46.8, 70</td><td>23.2, 46.8, 70</td></tr><tr><td>Velocity width (km s<sup>-1</sup>)</td><td>~200, ~400, ~600</td><td>~200, ~400, ~600</td></tr></tbody></table>

**Note.** Gaussian sizes are utilized for the injected artificial sources. All sizes refer to the Gaussian FWHM. The convolved sizes are the injected sizes in the Natural-mosaic. These are fixed in COSMOS, while in GOODS-N, because of the larger beam differences across the mosaic, we injected sources of fixed intrinsic sizes and convolved them to the local beam size appropriate for each mosaic position.

measured. We define the “flux-factor” as the ratio of the measured line flux to the injected flux. Therefore, the distribution of flux-factors captures both flux corrections and uncertainties on our flux estimations. The purpose of the flux-factor analysis is not just to correct for potential biases in our flux extraction procedure but also to estimate the uncertainty of the flux recovery. We subsequently utilize these flux probability distributions to inform our luminosity function estimates (Paper II). We ignore dependencies on frequency or position of the injected source flux-factors (determined as a function of local S/N) and completeness, thereby obtaining average values that correctly sample the data for an approximately uniform distribution of real sources in our cube.

Both the completeness and the flux-factors are dependent on the source size and line FWHM. Since we only inject sources of three spatial sizes and three frequency widths, we develop a probabilistic framework to relate each line candidate to the different injected sizes (Table 8). For each detected candidate, based on the template size and velocity width where their S/N peaks, we determine a probability distribution of belonging to each category of “injected” spatial size and line width, therefore matching in a continuous and probabilistic way the measured sizes to a discrete grid of intrinsic properties, as explained in detail in the following.

#### F.3.1. Flux-factors Based on Artificial Sources

The artificial source analysis allows us to estimate how well our measured fluxes correspond to the injected flux for candidates of different S/N, spatial size, and velocity width. The objective of the flux-factor analysis is to characterize the uncertainty and bias of our flux estimates in order to correctly estimate the uncertainty of our luminosity function measurement.

We confirm that aperture fluxes have a slight bias toward higher fluxes because positive noise adjacent to a candidate source tends to enlarge the fitted sizes and therefore contribute spurious flux to the candidate (Condon 1997). We thus need to estimate the magnitude of this bias at the S/N range of interest (~5–6) to correct the measured fluxes accordingly. A correction factor relies on an estimate of how likely a measured extended source is to be due to noise rather than real extended structure. In order to determine this bias, we need to assume an expected approximate size distribution for our sources to be combined with information from the artificial sources regarding how the flux is affected by the interplay of real and measured sizes.
<figure><img src="image_29.png" /><figcaption>Figure 18. Probability of injected spatial size (proxy for “real” size) as a function of S/N for different measured sizes (injected sizes are color-coded: smallest in blue, intermediate in green, significantly extended in red; sizes of the artificial sources listed in Table 8). Measured sizes are indicated by the spatial size of the peak template (see Table 5). Left: results for the COSMOS field. Right: results for the GOODS-N field. A measured (peak template) size of −1" corresponds to a source that achieves its peak S/N in the Natural-mosaic, i.e., before any smoothing, at the native resolution, while the 0" template is to a point-source template applied to the Smoothed-mosaic.</figcaption></figure>

We use the artificial sources to determine the probabilities of spatial extension (probability of being like spatial bin 0, 1, or 2 of the injected sources; Table 8) given the measured size as traced by the size of the spatial template that gives the highest S/N. This probability is used in the following to relate the measured properties of each line candidate to the flux-factor and completeness, which are computed for the bins of injected properties. We use the Bayes theorem to relate the probability of a given real size conditional to a measured size: p(real-injected size|measured size) to the probability distribution that we can measure from the artificial sources, which is the probability of measuring a given size conditional to a certain injected size p(measured size|injected size), by employing a prior on the expected real size distribution.33 We must employ a prior for the probabilities of the real sizes that captures our expectation that most sources would be unresolved, while allowing for a fraction of resolved sources coming from extended gas reservoirs and merging and/or blended objects. We adopt 88%, 10%, and 2% for the size bins in Table 8, respectively. We stress that although these relative fractions are uncertain, their precise choice does not significantly affect any of the results, because their effect is simply to modulate our assignment between line candidates and injected sources. The main effect of this assignment is in the estimation of completeness corrections, where the uncertainty introduced by these priors is small compared to the systematic uncertainty introduced by choosing a weighting based on the detected size distribution.

We estimate p(measured size|injected size) by measuring the fraction of injected sources of a given size recovered at different match-filter spatial template sizes. In this way, we compute the final posterior probability for a given measured size to originate from an unknown “injected” size by combining this with the prior (Figure 18):

$$ \begin{gathered} p(\text{real/injected size}|\text{measured size, S/N}) \\ \propto p(\text{real size}) \times p(\text{measured size}|\text{injected size, S/N}). \end{gathered} \tag{17} $$

We follow the same procedure for relating the measured velocity width (from the peak template) to the injected line widths by assuming a flat prior for the line width over the three

injected bins (Table 8). These are required to compute the completeness of line candidates by relating their measured properties to the injected sources. The derived probability distribution functions are shown in Figure 19.

In order to calculate flux-factors from the artificial sources, we employ an analogous technique. We calculate probability distributions of the flux-factors, given source spatial template size and S/N, by weighing the distribution of flux-factors found for given measured sizes by the probability that the given measured size originates from the different possible injected sizes. Specifically, the correction ratio depends on both the injected and the measured size, as a larger ratio is needed to correct for a compact source that appears extended:34

$$ \begin{aligned} p(\text{flux ratio}|\text{measured size, S/N}) &= \sum_{\text{injected}} p(\text{flux ratio}|\text{measured, injected, S/N}) \\ & \quad \times p(\text{injected}|\text{measured, S/N}). \end{aligned} \tag{18} $$

These flux-factor distributions (Figure 20) can be approximated by lognormal distributions peaking near a factor of one but with a tail to larger ratios to correct for the bias toward larger spatial size, which is introduced by including positive noise as part of the candidate source. We stress that this is just a convenient parameterization of the measured distributions. We calculate the probability distribution of flux ratios for each measured spatial size in bins of S/N. Since the distributions are noisy due to the limited number of artificial sources and do not strongly depend on S/N in the narrow 5–6 range of interest for our candidates, we consider the mean distribution over the 5 < S/N < 6 range (Figure 20). The fitted lognormal curves to these mean distributions, which provide a good interpolation to the noisy distribution estimates, will be utilized in the construction of the luminosity function. We can understand the general trend seen in these shapes as follows. The smallest template selects point sources. Therefore, the distribution peaks near 1. Slightly extended sources have a larger mean flux correction, reflecting the finding that they are most likely noise-smearied point sources, and so their flux needs to be reduced.

<footnote n="33">33 Note that the artificial sources can only be used to estimate distributions conditional to a given injected size, because the relative frequency of injected sizes that was utilized (uniform) is not representative of the expected distribution of real sizes.</footnote>

<footnote n="34">34 The only distribution that can be directly estimated through counting artificial sources is conditional to injected size. We cannot marginalize over those sizes without specifying a distribution of expected source sizes first.</footnote>
<figure><img src="image_30.png" /><figcaption>**Figure 19.** Probability of injected velocity width (proxy for “real” width) as a function of S/N for different measured widths (injected sizes are color-coded: narrowest in blue, intermediate in green, wide in red; FWHMs of the artificial sources are listed in Table 8) in the COSMOS (left) and GOODS-N (right) fields. The measured FWHMs are indicated by the frequency size of the peak template (see Table 5).</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 20.** Probability distribution functions of the “flux-factor” (i.e., the ratio of measured-to-injected line flux) for the range $5 < \text{S/N} < 6$ (in steps of $0.1\sigma$) conditional to the measured spatial size, as indicated by the spatial size of the peak template. We also show (black lines) the mean distribution taken over the full S/N range to obtain a more representative estimate. We fit this probability distribution by a lognormal distribution and show the best-fit model in red. On the left, we show the results for the COSMOS field, and on the right is the GOODS-N field. A measured (peak template) size of $-1''$ corresponds to a source that achieves its peak S/N in the Natural-mosaic, i.e., before any smoothing, at the native resolution, while the $0''$ template is to a point-source template applied to the Smoothed-mosaic.</figcaption></figure>

Slightly larger (intermediate-size) sources then require less correction because it becomes more likely that they are somewhat extended in reality. Even larger sources show a long tail of larger flux-factors because it is extremely unlikely that the real source is very extended. Therefore, their fluxes need to be significantly corrected (or rather, there is significant uncertainty as to their real flux, and we need to account for this in constructing the luminosity function).

### F.3.2. Completeness

In order to estimate the completeness of our detection process, we utilize the artificial sources to measure the fraction of the injected sources that are detected. The

objective of the completeness correction is to account for the fraction of the mosaic volume where a given line candidate would be detectable and for the fraction of objects of a given intrinsic line luminosity that would be missed by a fixed S/N threshold.

We assume that the fraction of detected lines (a proxy for the probability of detection) only depends on the integrated line flux, spatial size, and velocity width. By injecting artificial sources that uniformly sample random positions within the edges of the mosaic, we derive completeness corrections that account for the effects of the spatial and frequency variation of the sensitivity, as previously adopted by Walter et al. (2016).
<figure><img src="image_32.png" /><figcaption>**Figure 21.** Completeness corrections calculated as a function of the injected flux of the artificial sources. Top: completeness in the COSMOS field. Bottom: completeness in the shallower GOODS-N field. Colors distinguish the three different velocity widths of the artificial sources, and the line (marker) style distinguishes different spatial sizes. Markers represent the measured completeness in bins of line flux, and lines represent the best-fitting interpolating functions.</figcaption></figure>

While completeness is a property of the overall number counts in a luminosity bin, we partially adopt the spirit of the $1/V_{max}$ method of calculating a “per-source” completeness only insofar as this depends on the line velocity width and especially spatial size of the detection. While this may potentially introduce a bias<fnref n="35" /> (see, e.g., Hogg et al. 2010), it saves us from additional assumptions about the size distribution. We therefore point out the important caveat that our luminosity function estimate does not correct for missed objects due to poor sampling of the size distribution. The effect can be seen by noticing, for example, that the completeness for extended objects at low flux values drops very quickly. This is reflected in the fact that all of our low-flux objects are point sources (therefore, the completeness correction in the lowest line luminosity bins misses the potential contribution from undetected extended sources).

The completeness is measured from the artificial sources as a function of “injected” properties, i.e., injected integrated flux,

and spatial size and frequency width (in three bins each; see Table 8) as the ratio of injected sources recovered with S/N above a threshold value of $5\sigma$ to the total number of injected sources (the precise choice of a threshold does not change the result appreciably). These measured completeness values are shown in Figure 21, together with the interpolating functions that we use in deriving the luminosity function: the two-parameter ($I_0$ and $I_1$) family of functions $1 - e^{-I/I_0}/(I + I_1)$, where $I$ is the integrated line flux. While this chosen family of interpolating functions has no specific significance, we found it to provide an appropriate description of the measured completeness.

The optimal way to correct for the completeness of a luminosity function bin would require calculating the mean completeness within the bin over the full “internal” parameter space (in our case, these are spatial and velocity line sizes, frequency, and precise luminosity within the bin) weighted by the model expectation for the distribution of sources within this space. We adopt an intermediate approach between this “mean completeness” approach and a purely “per-candidate” approach. We assume that the distribution of line luminosity within a bin and the frequency distribution are uniform and average over this subset of the internal space by randomly sampling it. On the other hand, we do not assume an intrinsic spatial size and velocity width distribution in order to avoid biasing our result. We therefore adopt the sizes and line widths of the candidates to calculate the appropriate completeness, thereby letting the data determine the size and velocity width distributions.

In summary, for each detected line candidate, the size properties (i.e., the spatial size and frequency width) of the “real” underlying source are estimated probabilistically based on the measured size and width (from the peak template). Then, probability-weighted completenesses are the factors that enter the evaluation of the luminosity function. The dependences of these completenesses on the line flux are mean values evaluated for the full luminosity function bin, rather than depending on the candidate line flux measurement. This hybrid approach helps us mitigate the bias that would derive from a purely “per-source” correction (e.g., Hogg et al. 2010).

### F.4. Implementation of the Statistical Corrections

In order to assemble the luminosity function, we use a variation on the method used by Decarli et al. (2014, 2016). We weigh the contribution of each line candidate to the luminosity bin by its purity and inversely by its completeness and use the total cosmic comoving volume covered within the edges of the mosaic. The completeness correction converts this volume to an effective $V_{max}$ for each galaxy, also accounting for the spatial variation in sensitivity. In order to estimate the systematic uncertainty introduced by our assumptions, we evaluate the luminosity function with many random realizations of flux-factor, purity assignment, and luminosity bin widths and boundaries. One of the differences between our approach and the approach employed by Decarli et al. (2014, 2016) consists of including a larger number of candidates. We have also calculated the luminosity function using the same method employed by Decarli et al. (2016), and the result is consistent with our more extensive method. The advantage of our approach consists of relying less heavily on the properties of the few moderate-S/N individual candidates that happen to be located near

<footnote n="35">³⁵ Regions of parameter space that have very low completeness tend to be poorly sampled and hence cannot properly be accounted for in a completeness calculation that is weighted by the detected candidates.</footnote>
the top of the S/N list but that still have a limited probability of being real. Although a fraction of the moderate-S/N candidates are expected to be real, it is not clear that those near the top (of the set of uncertain candidates) of the S/N list have a significantly greater likelihood of being real given the limited range in S/N considered (~5–6). By utilizing a larger sample of candidates in deriving constraints on the luminosity function down-weighted by appropriate purities, we do not introduce additional bias but rather better explore the implications of the systematic uncertainties. In particular, the statistical justification for a “per-source” purity and completeness correction (which we cannot fully avoid) only holds for large enough samples. By better sampling the “internal” space of possible candidate sizes and line widths and adopting average per-bin completenesses (i.e., not using the uncertain measured fluxes), we aim to achieve a more accurate completeness correction and evaluation of the systematic uncertainties introduced by these factors.

In detail, for each $L'_{CO}$ bin, we calculate a completeness factor appropriate for each of the nine bins in the spatial size-frequency width grid by averaging over 1000 random realizations of values of $L'_{CO}$ in each bin, using random redshifts to calculate the corresponding integrated flux and hence the appropriate completeness correction for each. We average over this frequency distribution and precise $L'_{CO}$ within each bin. We therefore enforce a uniform prior and maintain the dependence on the spatial-frequency size information separate. We subsequently apply them for each line candidate as weighted by the candidate probability distribution for its spatial-frequency size assignment. In this way, we use the measured relative occurrence of different spatial and velocity sizes as weights for the appropriate completeness for each luminosity function bin.

In order to explore the range of luminosity function values allowed by our systematic uncertainty, we use 10,000 Monte Carlo realizations of the luminosity function calculation for each bin width and shift, where we vary the purity assignment independently for each candidate and the flux-factor to be applied. We therefore “move candidates around” among adjacent luminosity bins, simulating the effect of the uncertainty in their intrinsic fluxes. We separately implement the purity in one of the two ways we described in Appendix F.1, either as normally or uniformly distributed. We also implement the flux correction as taking a random value drawn from the appropriate lognormal distribution, which was derived for each spatial size in the previous section. We also add a normal uncertainty of 20% to the measured flux to reproduce the uncertainty in our flux calibration.

Finally, in order to describe the range of values for the luminosity function (in log comoving volume density space) spanned by our 10,000 Monte Carlo realizations, we calculate the median value for each bin and a measure of the scatter around the median. We evaluate the scatter conservatively by quoting luminosity function ranges that include 90% of the probability. We also evaluate the statistical Poisson uncertainty as appropriate for each bin as a relative uncertainty of $1/\sqrt{N}$, where N corresponds to the number of candidates in an $L'_{CO}$ bin.

# Appendix G
## Additional Stacking Results of Galaxies with Optical Redshifts

In addition to the stacking of individual sets of galaxies described in the main text, we have attempted to stack three more sets of galaxies based on optical redshift information. We do not detect significant signal in these stacked spectra.

### G.1. COSMOS Protocluster at $z = 5.3$

We search for CO(2–1) from nine member galaxies of the AzTEC-3 protocluster, as identified through the Lyman Break (LBG) technique and color selection by Capak et al. (2011). While two of them show a hint of a CO emission line signal (~2.5σ), all others (including LBG-1, the one with the strongest [C II] emission; Riechers et al. 2014 and 2018 Capak et al. 2015) are consistent with noise. The positions of the LBGs with tentative CO detections are J2000 10:00:21.96+02:36:08.5 and 10:00:20.13+02:35:53.9. Assuming a line FWHM of 250 km s⁻¹ (i.e., the width of the [C II] line in LBG-1), we derive CO line fluxes of approximately 0.03 ± 0.012 Jy km s⁻¹, corresponding to $L'_{CO21} = 7 \pm 3 \times 10^9$ K km s⁻¹ pc². We do not claim any detections due to the low significance and only quote these as approximate limits for reference. We also stack all spectra at the positions of LBGs (limited to the seven galaxies that are not contaminated by emission from the bright CO(2–1) line in AzTEC-3 at the resolution of our survey) and do not detect a significant signal (Figure 22). We therefore place a 3σ limit of <0.012 Jy km s⁻¹ on their average CO(2–1) emission, corresponding to $L'_{CO21} < 3 \times 10^9$ K km s⁻¹ pc². The average stellar mass for the stacked LBGs, as reported by the COSMOS2015 catalog, is $4 \times 10^9 M_{\odot}$. If we assume an $\alpha_{CO} = 3.6 M_{\odot}$ (K km s⁻¹ pc²)⁻¹, we obtain a gas mass upper limit of <1.1 × 10¹⁰ $M_{\odot}$ and therefore a gas fraction of $M_{gas}/M_{*} < 3$, which is not strongly constraining. None of the LBGs are detected in 3 GHz radio continuum emission (Smolčić et al. 2017). This would place a limit of <165 $M_{\odot}$ yr⁻¹ on their star formation rate if we adopted the redshift evolution of the radio–FIR correlation measured by Delhaize et al. (2017). This may be converted to an expected limit of $L'_{CO} < 2.2 \times 10^{10}$ K km s⁻¹ pc² by assuming the star formation law (Daddi et al. 2010b; Genzel et al. 2010). This limit is higher than what we derive from our CO nondetection, implying that our deep observations provide strong constraints on the CO luminosity of $z > 5$ LBGs.

### G.2. Grism Redshifts in the GOODS-N Field

The 3D-HST catalog (Momcheva et al. 2016) contains 694 galaxies in GOODS-N with grism redshifts for which the CO(1–0) line is covered by our data. Nevertheless, the majority of these grism spectra do not significantly improve the redshift determination over the photometric redshift and are therefore not usable for stacking. We search for matches to our line search candidates (down to $4\sigma$) within a radius of 2" and ~500 km s⁻¹. We find 16 potential matches and assess the contamination by chance association by also matching our blind detection catalog to 694 random positions within the signal data cube. We find that the distribution due to random associations is well described by a Gaussian with a mean 8 associations and a standard deviation of 3.5 associations. Therefore, the majority of our line associations are likely to be random, but some may be expected to be real. The measured
<figure><img src="image_33.png" /><figcaption>**Figure 22.** Additional spectral stacks of sets of galaxies with spectroscopic or grism redshifts. Top: CO(1–0) stacks of subsets of galaxies with the highest-quality grism spectra with redshift uncertainty <500 (left) and <700 km s⁻¹ (right). Bottom: CO(2–1) stacks for samples of potential z > 5 galaxies belonging to previously identified overdensities in our fields, i.e., the AzTEC-3 protocluster at z = 5.3 (left) and the overdensity around HDF850.1 (right; Capak et al. 2011; Walter et al. 2012). The spectral resolution is the same as in Figure 5.</figcaption></figure>

line fluxes of these candidate counterparts would imply gas masses that are sometimes larger than the stellar masses. This is possible, but we consider it more likely that those cases may be random noise associations.

We also stack the spectra extracted at the positions of the galaxies with high-quality grism redshifts (Figure 22). Thirty-seven of them have redshift uncertainties less than 500 km s⁻¹ (based on the 95th percentile of the redshift probability distribution reported by Momcheva et al. 2016), and 35 more have redshift uncertainties less than 700 km s⁻¹. At this level of uncertainty, it would be likely that a fraction of the line signal contributes to the stacked spectrum. Both of these stacks show no detection. Assuming a line FWHM of 300 km s⁻¹, this implies 3σ upper limits of <0.011 and <0.008 Jy km s⁻¹, corresponding to ≲3 and 2 × 10⁹ K km s⁻¹ pc², respectively.

$$G.3. \quad \text{HDF850.1 } z = 5.2 \quad \text{Galaxy Overdensity in GOODS-N}$$

We also search for potential CO(2–1) emission from galaxies in the z ~ 5.2 overdensity around the SMG HDF850.1 (Walter et al. 2012), taking advantage of the abundance of available spectroscopic redshifts in this region. Twenty-four of the 105 galaxies with spectroscopic redshifts presented in that work fall within our data, but none of them are individually detected. We stack the spectra to obtain an average spectrum (Figure 22). No significant emission is found after stacking. Assuming a line FWHM of

300 km s⁻¹ implies a 3σ upper limit of <0.015 Jy km s⁻¹, and L'CO21 < 3.5 × 10⁹ K km s⁻¹ pc² at z ~ 5.2. We match these 105 galaxies to galaxies within 1" in the photometric catalog by Skelton et al. (2014), finding 83 matches, but we do not adopt their stellar mass estimates because the redshifts of these galaxies were often greatly underestimated.

## ORCID iDs

Riccardo Pavesi https://orcid.org/0000-0002-2263-646X
Chelsea E. Sharon https://orcid.org/0000-0002-6250-5608
Jacqueline A. Hodge https://orcid.org/0000-0001-6586-8845
Roberto Decarli https://orcid.org/0000-0002-2662-8803
Fabian Walter https://orcid.org/0000-0003-4793-7880
Chris L. Carilli https://orcid.org/0000-0001-6647-3861
Emanuele Daddi https://orcid.org/0000-0002-3331-9590
Ian Smail https://orcid.org/0000-0003-3037-257X
Rob J. Ivison https://orcid.org/0000-0001-5118-1313
Mark Sargent https://orcid.org/0000-0003-1033-9684
Elisabete da Cunha https://orcid.org/0000-0001-9759-4797
Manuel Aravena https://orcid.org/0000-0002-6290-3198
Jeremy Darling https://orcid.org/0000-0003-2511-2060
Vernesa Smolčić https://orcid.org/0000-0002-3893-8614
Nicholas Z. Scoville https://orcid.org/0000-0002-0438-3323
References

Aravena, M., Carilli, C., Daddi, E., et al. 2010, ApJ, 718, 177

Aravena, M., Carilli, C. L., Salvato, M., et al. 2012, MNRAS, 426, 258

Aravena, M., Hodge, J. A., Wagg, J., et al. 2014, MNRAS, 442, 558

Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2013, ApJ, 769, 80

Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, PhR, 369, 111

Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2015, ApJ, 809, 175

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bond, J. R., & Efstathiou, G. 1987, MNRAS, 226, 655

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bourne, N., Dunne, L., Bendo, G. J., et al. 2013, MNRAS, 436, 479

Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Natur, 522, 455

Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Natur, 470, 233

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

Carilli, C. L., McKinnon, M., Ott, J., et al. 2015, arXiv:1510.06438

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

Casey, C. M., Hodge, J. A., Lacy, M., et al. 2015, arXiv:1510.06411

Casey, C. M., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Chen, C.-C., Hodge, J. A., Smail, I., et al. 2017, ApJ, 846, 108

Colombi, S., Davis, O., Devriendt, J., Prunet, S., & Silk, J. 2011, MNRAS, 414, 2436

Condon, J. J. 1997, PASP, 109, 166

Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

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

da Cunha, E., Walter, F., Decarli, R., et al. 2013b, ApJ, 765, 9

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Bournaud, F., Walter, F., et al. 2010a, ApJ, 713, 686

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

Daddi, E., Dannerbauer, H., Elbaz, D., et al. 2008, ApJL, 673, L21

Daddi, E., Dannerbauer, H., Krips, M., et al. 2009, ApJL, 695, L176

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

Daddi, E., Elbaz, D., Walter, F., et al. 2010b, ApJL, 714, L118

Dannerbauer, H., Walter, F., & Morrison, G. 2008, ApJL, 673, L127

Darling, J., & Zeiger, B. 2012, ApJL, 749, L33

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

Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJ, 782, 78

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

Eales, S., Smith, M. W. L., Auld, R., et al. 2012, ApJ, 761, 168

Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJL, 600, L93

González-López, J., Bauer, F. E., Aravena, M., et al. 2017, A&A, 608, A138
Gregory, P. 2010, Bayesian Logical Data Analysis for the Physical Sciences
(Cambridge: Cambridge Univ. Press)

Gregory, P. C., & Loredo, T. J. 1992, ApJ, 398, 146

Greve, T. R., Bertoldi, F., Smail, I., et al. 2005, MNRAS, 359, 1165

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Groves, B. A., Schinnerer, E., Leroy, A., et al. 2015, ApJ, 799, 96

Guidetti, D., Bondi, M., Prandoni, I., et al. 2017, MNRAS, 471, 210

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Hogg, D. W., Myers, A. D., & Bovy, J. 2010, ApJ, 725, 2166

Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913

Keating, G. K., Bower, G. C., Marrone, D. P., et al. 2015, ApJ, 814, 140

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

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kogut, A., Chuss, D. T., Dotson, J., et al. 2014, Proc. SPIE, 9143, 91431E

Kogut, A., Fixsen, D. J., Chuss, D. T., et al. 2011, JCAP, 7, 025

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Lentati, L., Wagg, J., Carilli, C. L., et al. 2015, ApJ, 800, 67

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

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

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6

Mangum, J. G., Darling, J., Henkel, C., & Menten, K. M. 2013, ApJ, 766, 108

Mangum, J. G., Darling, J., Menten, K. M., & Henkel, C. 2008, ApJ, 673, 832

McKinnon, M., Carilli, C., & Beasley, T. 2016, Proc. SPIE, 9906, 990627

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Miettinen, O., Novak, M., Smolčić, V., et al. 2017, A&A, 602, A54

Molnár, D. C., Sargent, M. T., Delhaize, J., et al. 2018, MNRAS, 475, 827

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Morrison, G. E., Owen, F. N., Dickinson, M., Ivison, R. J., & Ibar, E. 2010, ApJS, 188, 178

Narayanan, D., & Krumholz, M. R. 2014, MNRAS, 442, 1411

Pope, A., Scott, D., Dickinson, M., et al. 2006, MNRAS, 370, 1185

Popping, A., Jurek, R., Westmeier, T., et al. 2012, PASA, 29, 318

Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398

Popping, G., van Kampen, R., 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., Capak, P. L., Carilli, C. L., et al. 2010a, ApJL, 720, L131

Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84

Riechers, D. A., Carilli, C. L., Maddalena, R. J., et al. 2011a, ApJL, 739, L32

Riechers, D. A., Carilli, L. C., Walter, F., et al. 2011b, ApJL, 733, L11

Riechers, D. A., Carilli, C. L., Walter, F., & Momjian, E. 2010b, ApJL, 724, L153

Riechers, D. A., Hodge, J., Walter, F., Carilli, C. L., & Bertoldi, F. 2011c, ApJL, 739, L31

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2018, ApJ, submitted (arXiv:1808.04371)

Riechers, D. A., Walter, F., Carilli, C. L., et al. 2006, ApJ, 650, 604

Righi, M., Hernández-Monteagudo, C., & Sunyaev, R. A. 2008, A&A, 489, 489

Rudnick, G., Hodge, J., Walter, F., et al. 2017, ApJ, 849, 27

Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5

Schinnerer, E., Groves, B., Sargent, M. T., et al. 2016, ApJ, 833, 112

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Selina, R., & Murphy, E. 2017, ngVLA Reference Design Development and Performance Estimates, http://library.nrao.edu/public/memos/ngvla/NGVLA_17.pdf

Sharon, C. E., Riechers, D. A., Hodge, J., et al. 2016, ApJ, 827, 18

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 807, 128

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Spilker, J. S., Bezanson, R., Marrone, D. P., et al. 2016, ApJ, 832, 19

Steidel, C. C., Shapley, A. E., Pettini, M., et al. 2004, ApJ, 604, 534

Steinhardt, C. L., Speagle, J. S., Capak, P., et al. 2014, ApJL, 791, L25

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Natur, 463, 781

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246

Tacconi, L. J., Neri, R., Chapman, S. C., et al. 2006, ApJ, 640, 228

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Vio, R., & Andreani, P. 2016, A&A, 589, A20

Vio, R., Vergès, C., & Andreani, P. 2017, A&A, 604, A115

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

Walter, F., Decarli, R., Carilli, C., et al. 2012, Natur, 486, 233

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

Westmeier, T., Popping, A., & Serra, P. 2012, PASA, 29, 276

Whiting, M. T. 2012, MNRAS, 421, 3242

Zeiger, B., & Darling, J. 2010, ApJ, 709, 386