Hanae Inami

Roberto Decarli

Fabian Walter

Axel Weiss

Chris Carilli

Manuel Aravena

Leindert Boogaard

Jorge Gonzaĺez-López

Gergö Popping

Elisabete Da Cunha

Jorge González-López

Roland Bacon

Franz Bauer

Thierry Contini

Paulo C Cortes

Pierre Cox

Emanuele Daddi

Tanio Díaz-Santos

Melanie Kaasinen

Dominik A Riechers

Jeff Wagg

Paul Van Der Werf

Lutz Wisotzki

The ALMA Spectroscopic Survey in the Hubble Ultra Deep Field: Constraining the Molecular Content at log( M * / M ) ∼ 9.5 with CO Stacking of MUSE-detected z ∼ 1.5 Galaxies

Keywords: galaxies, high-redshift -galaxies, evolution -galaxies, ISM -galaxies, star formation

INTRODUCTION

Stars form inside dense molecular gas clouds. It is thus critical to reveal how much molecular gas exists in galaxies to characterize galaxy formation and evolution (e.g., Kennicutt & Evans 2012;Carilli & Walter 2013;Tacconi et al. 2020;Hodge & da Cunha 2020). It is well established that the cosmic star formation rate density increased from the early stages of the Universe towards its peak around z ∼ 1 -3, after which it progressively decreased until the current epoch [START_REF] Madau | [END_REF]. A broad consensus is emerging on the cause of this growth, peak, and decline of the star formation history over cosmic time via measurements of the gas that fuels star formation (e.g., Walter et al. 2014;Decarli et al. 2016a;Scoville et al. 2017;Decarli et al. 2019;Riechers et al. 2019;Liu et al. 2019;Magnelli et al. 2020;Lenkić et al. 2020;Tacconi et al. 2020). This evolution could be due to either the available supply of molecular gas for forming stars, a mechanism that causes high efficiency in star formation such as galaxy mergers, or a combination of these processes.

Most of the star formation in the Universe occurs in galaxies residing on the so-called "main sequence", a tight correlation between the star formation rate (SFR) and stellar mass (M * ) of star-forming galaxies (e.g., Noeske et al. 2007;Salmi et al. 2012;Whitaker et al. 2012Whitaker et al. , 2014;;Schreiber et al. 2015;Popesso et al. 2019a,b). This correlation is observed at redshifts up to at least z ∼ 6.5 (e.g., Speagle et al. 2014;Salmon et al. 2015). In other words, most star formation in the Universe is longlasting and evolves steadily, supporting the idea that it is predominantly regulated by the gas accretion of the available fuel supply and feedback processes, rather than stochastic events like galaxy mergers (Lilly et al. 2013;Dekel & Krumholz 2013;Tacchella et al. 2016). Such stochastic events cause enhanced star formation activity, which elevates SFRs significantly above the mainsequence relation. However, these galaxies, referred to as starbursts are in the minority (e.g., Rodighiero et al. 2011;Sargent et al. 2012;Lutz 2014). To understand these two star formation modes and the efficiency of the star formation process, it is essential to determine the H 2 gas supply and deficiency. The most common tracer of the molecular gas, the fuel of star formation, is line emission from carbon monoxide ( 12 CO) rotational transitions at low excitation (Bolatto et al. 2013;Carilli & Walter 2013).

The connection between molecular gas content, stellar mass, and SFR for z > 1 galaxies on and above the main-sequence relation has been investigated in various targeted studies (e.g., Magdis et al. 2012;Santini et al. 2014;Scoville et al. 2016;Tacconi et al. 2013Tacconi et al. , 2018)). However, the targets have so far been limited to massive (> 10 10 M ) and highly star-forming ( 50M yr -1 ) galaxies. Little is known about the molecular gas reservoirs in galaxies either below the main sequence or at modest stellar masses, mostly because of sensitivity limits.

The ALMA Spectroscopic Survey in the Hubble Ultra-Deep Field (ASPECS, Walter et al. 2016) has been conducted as a spectroscopic survey over the entire frequency range of ALMA Bands 3 and 6 in the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006) to perform an unbiased search for multiple rotational transitions of CO emission (Walter et al. 2016;González-López et al. 2019). Spectral line scans have an advantage in assessing the molecular gas content based on a complete line flux-limited sample without any target preselections. González-López et al. (2019) conducted a blind search of line and continuum sources directly in the ASPECS data cube (Band 3) and evaluated its completeness. Using these reliable CO detections, Decarli et al. (2019) constructed CO luminosity functions and presented the evolution of the cosmic gas mass density. A census of the molecular gas content of galaxies that have direct CO detections is shown and discussed in Aravena et al. (2019). These gas measurements were compared with model predictions from cosmological sim-ulations and semi-analytical models in Popping et al. (2019). Uzgil et al. (2019) performed a power spectrum analysis and probed CO emission at 1 z 4 below the sensitivity limit of individual detections and gave a constraint on missing CO emission from individually undetected galaxies.

Here, we maximize the sensitivity of the ASPECS data to detect CO emission by stacking ALMA spectra using the Band 3 data. Our stacking analysis is based on optical spectroscopic redshifts from another large, unbiased, blind spectroscopic survey in the HUDF carried out by the integral field unit (IFU) instrument MUSE (Multi Unit Spectroscopic Explorer) on the Very Large Telescope (Bacon et al. 2017). The combination of the three-dimensional (3D) data obtained by both the AS-PECS and MUSE surveys not only made the stacking analysis possible, but also enabled a direct comparison between the molecular gas properties and rest-frame optical/ultraviolet properties (Boogaard et al. 2019).

This paper is structured as follows: we first introduce the observations and data taken with the ASPECS and MUSE HUDF surveys and the ancillary data in §2. We describe the method of the stacking analysis along with the sample used for stacking in §3. In §4, the stacked CO spectra are presented. We then convert the measured CO emission to molecular gas mass and discuss the gas mass content in §5. We summarize and conclude our findings in §6. A flat ΛCDM cosmology with H 0 = 70 km s -1 Mpc -1 , Ω Λ = 0.7, and Ω m = 0.3 is adopted throughout this paper.

THE ASPECS AND MUSE DATASETS

The MUSE survey covered the entire area of the HUDF (Bacon et al. 2017), whereas the ASPECS LP observed almost the entire region of the Hubble eXtremely Deep Field (XDF; Illingworth et al. 2013). MUSE is an optical IFU (Bacon et al. 2015) on the Very Large Telescope (VLT) Yepun (UT4) of the European Southern Observatory (ESO) with a wide fieldof-view (FoV, 1 × 1 ), high sensitivity, wide wavelength coverage (4650 -9300 Å), and high spectral resolution (R ∼ 3000). In particular, the large size of its FoV facilitates spectroscopic redshift surveys without requiring any target preselection, achieving a spatially homogeneous spectroscopic completeness. It offers redshift determination at z = 0-6.5 based on rest-frame ultraviolet and optical emission and absorption lines.

Complementing the MUSE spectroscopic survey, the ASPECS line scan survey in ALMA Band 3 (3mm) can detect multiple CO rotational transition lines from J = 1 -0 to J = 6 -5 in the redshift range of z = 0 -7 (with gaps 1 at 0.37 < z < 1.00 and 1.74 < z < 2.00, Walter et al. 2016). The MUSE redshift coverage overlaps well with the redshift coverage of the ASPECS Band 3 line scan survey at z > 1 (Figure 1), except for a small gap at z = 1.74 -2.00.

In this paper, we concentrate on the ASPECS Band 3 data where the full MUSE redshift coverage can be exploited (see §3.2). The significant redshift overlap of these two unbiased spectroscopic surveys in the same field is beneficial for performing stacking analyses on ASPECS CO spectral lines based on the optically determined MUSE redshifts. In this section, we will briefly present the survey designs of the HUDF conducted by ASPECS and MUSE.

ASPECS observations and data

The detailed survey strategy and data reduction of the ASPECS Pilot and Large Program surveys are presented in Walter et al. (2016) and González-López et al. (2019), respectively. The observational setups of the Large Program (LP) survey were the same as the Pilot survey, except its coverage is extended to ∼ 5 arcmin 2 . In this work, we used the combined data of the ASPECS Pilot and LP surveys taken with ALMA Band 3. AS-PECS LP carried out a full frequency scan in Band 3 (84 -115 GHz) over the XDF, which resides in the HUDF. With a total of 17 pointings centered at (RA, Dec) = (03:32:38.5, -27:47:00), which completely cover the Pilot region, the total area with a primary beam response > 50% in the LP survey is 4.6 arcmin 2 at ≈ 99.5 GHz (the central frequency of Band 3).

The CASA (Common Astronomy Software Applications) software was used to calibrate and image the data. With the C40-3 array configuration, we obtained a synthesized beam size of 1.75 × 1.49 with a position angle 91.5 deg at ≈ 99.5 GHz by using natural weighting in CASA when imaging the data. The frequency channel was rebinned to 7.813 MHz (23.5 km s -1 at 99.5 GHz). The sensitivity for this channel bin size was ∼ 0.2 mJy beam -1 throughout the entire scanned frequency range. The 5σ CO(2-1) line sensitivity is > 1.4 × 10 9 K km s -1 pc 2 assuming a line width of 200 km s -1 . Based on assumptions made in this work (see §5.1), the corresponding molecular gas (H 2 ) limit is 6.8 × 10 9 M . For our main target emission line, CO(2-1), the ALMA Band 3 scan offers redshift coverage of z = 1.0059 -1.7387.

Using the same data cube, González-López et al. ( 2019) reported CO emission line detections from an unbiased blind search without prior knowledge of source positions and observed CO frequencies. They found 16 high significance CO emitting sources, among which 11 were identified as CO(2-1) emission. We refer the readers to González-López et al. (2019) for comprehensive discussions on the detection methods. Aravena et al. (2019) assessed the molecular gas properties of these sources, which will be used for comparison in this paper.

MUSE observations and data

The MUSE-HUDF deep survey was conducted as a two layer survey of different depths (Bacon et al. 2017). The 3 ×3 deep survey observed the entire HUDF region, whereas the 1 × 1 ultra deep survey was carried out near the center of the deep survey area. The ≈ 10 h and ≈ 31 h exposure times, respectively, reached 3σ emission line flux limits of 3.1 and 1.5 × 10 -19 erg s -1 cm -2 at ∼ 7000 Å.

The MUSE spectra were extracted with two methods: prior extractions based on the UVUDF catalog (Rafelski et al. 2015) and a blind search for emission lines in the data cubes. In the former case, the coordinates for the source extraction were from UVUDF 2 , whereas for the latter, the coordinates were determined from the MUSE data. For a more detailed description of the survey and data treatments, see Bacon et al. (2017). Following Dunlop et al. (2017), we corrected the known systematic offset of the Hubble Space Telescope (HST) positions to match the radio astrometric reference frame by applying +0.279 in Declination (Dec) and -0.076 in Right Ascension (R.A.).

The MUSE spectroscopic redshifts were measured from the spectral features in the extracted spectra. Each redshift has an associated confidence level of 3 (secure redshift, determined by multiple spectral lines), 2 (secure redshift, determined by a single spectral line), or 1 (possible redshift, determined by a single spectral line whose spectral identification remains uncertain). The typical MUSE redshift uncertainty is σ z = 0.00012(1+z) or σ v ≈ 40 km s -1 , which is smaller than the typical line width of CO emission. We refer to Inami et al. (2017) for details about the redshift determination and redshift catalogs of the MUSE-HUDF field. In this work, we only use the reliable MUSE redshifts of confidence levels 2 and 3.

The MUSE-HUDF survey obtained 1338 reliable redshifts in total, a factor of eight increase over the previously available spectroscopic redshifts in this field. The simultaneous wavelength coverage of 4650 -9300 Å and the spectral resolution of R ∼ 3000 of MUSE allowed detections and unambiguous identifications of major restframe ultraviolet and optical emission lines, including Hα, [O II], [O III], and Lyα. These lines were used to determine redshifts over the range 0 < z < 6.5. Although more difficult to probe, the redshift range z ∼ 1.5 -3 can be covered by C III] and absorption features (e.g., C IV, Fe II, Mg II; see Figure 13 of Inami et al. 2017).

Among the spectroscopic redshifts assessed in the MUSE HUDF field, 503 sources with a confidence level of 2 or higher lie within the ASPECS survey region (LP primary beam response > 50%). Out of these sources, 107 sources are [O II] emitters covering z = 0.25 -1.5 and 363 sources are Lyα emitters covering z = 2.8 -6.6. For the analysis presented in the main part of this paper, we took advantage of the prevalent [O II] line detections and their redshift overlap with the CO(2-1) selection function in ALMA Band 3 (Figure 1) to perform a stacking analysis. In the redshift range where CO(2-1) can be detected, the MUSE spectroscopic redshift sources with absorption features and C III] emission also contributed to the stacking, although the number was small (Figure 1). We also attempted to stack spectra to detect higher-J CO lines. The MUSE redshifts used in higher-J CO stacked spectra were mostly measured using the Lyα line, which is known to be offset from the systemic redshift for a few hundred km s -1 (e.g., Shapley et al. 2003). We have applied a correction to this offset ( §4.2).

Ancillary data and physical parameters derived from SED fitting

Owing to the same HUDF coverage of the MUSE and ASPECS observations, there are abundant ancillary data available. We assembled optical and near-infrared photometric data from Skelton et al. (2014). These photometry catalogs and data include ultraviolet to infrared from the HST, various ground-based telescopes, and all of the Spitzer IRAC channels.

Based on this photometric dataset, we inferred physical parameters such as star formation rate (SFR) and stellar mass (M * ) via modeling with the high redshift extension of the Spectral Energy Distribution (SED) fitting code MAGPHYS (da Cunha et al. 2008, 2015). In addition to Spitzer/MIPS and Herschel/PACS, the ALMA 1.2 and 3 mm photometry from the ASPECS data (González-López et al. 2019, 2020) were also used for the SED fitting for a subset of the sources whose CO or continuum emission was detected. The same procedure was carried out in the other ASPECS work (Aravena et al. 2020;Boogaard et al. 2020 submitted). The SED fits were computed based on the MUSE spectroscopic redshift. See Boogaard et al. (2019) for the detailed process of the SED fitting.

METHODS

ALMA spectral extraction and stacking

From the ALMA data, we first extracted a sub-cube centered at the position of each MUSE source with a secure spectroscopic redshift. This sub-cube has a size of 11 ×11 ×3000 km s -1 . The primary beam correction using the combined Pilot and LP primary beam response was applied after the sub-cube extraction.

The uncertainty of the extracted spectra was calculated based on the region in the data cube where the combined beam response is > 99% of the peak sensitivity at the phase center. The uncertainty for each spectrum is then scaled by the combined primary beam response of the location of the objects under consideration.

To perform stacking, the channel frequencies of the extracted spectra were converted to the rest-frame, then the spectra were resampled onto a common frequency grid. In the rest-frame, a CO line detected at the lower frequency end of the spectrum (i.e., galaxies at higher redshift) has a coarser spectral sampling due to the (1 + z) correction of the redshift. Thus, the common frequency grid was given at the coarsest sampling, corresponding to the CO line detected at the lowest frequency end of Band 3. In the case of CO(2-1), we set the final spectral range and the channel width of the extracted spectra to be 3000 km s -1 , centered on the CO(2-1) rest-frame frequency (230.538 GHz), and 27.8 km s -1 (21.4 MHz), respectively. This resulted in 108 channels in the extracted spectra. The final extracted rest-frame frequency range is 229.395 -231.682 GHz (±1500 km s -1 ). The same procedure is adopted for stacking high-J CO lines.

This conversion was implemented by first shifting the channel frequencies of the spectra to the rest-frame, according to the MUSE redshift. We then applied a Gaussian decimation filter to the rest-frame spectra via the Fourier plane. This filter removes fine scale (high sampling frequency) noise in the spectrum which would otherwise be aliased into spurious noise when sub-sampled to the lower target resolution [START_REF] Lyons | Understanding Digital Signal Processing[END_REF]. A better signal-to-noise ratio (SNR) was thus obtained when we binned all spectra onto the coarser common frequency grid in the later step. In the Fourier plane, the width of the Gaussian filter was set to give an attenuation of 40 dB at the Nyquist folding frequency of half a cycle per channel of the new channel width (0.5 cycles per 27.8 km s -1 ). This value is a reasonable compromise between the desirable suppression of aliased noise and the corresponding (minor) reduction in the resolution of the sub-sampled spectra. The Gaussian kernel was scaled to be unity at the origin of the Fourier plane to guarantee flux conservation. In addition, the kernel was given an odd number of elements, placed symmetrically around the origin, to prevent spectral shifts during the convo- lution. The effect of this filter was equivalent to performing a linear convolution of the spectrum with an area-normalized Gaussian of FWHM 63.2 km s -1 .

Next, we interpolated the filtered spectrum onto the frequency grid of the stacked spectra. The mean of the resulting spectra was then calculated to obtain the final stacked spectrum3 . We did not apply any weighting when performing the stack.

We also carried out two different random extractions following the same procedure described above, to produce random stacked spectra for comparisons. One random extraction involved assigning random redshifts to each spectrum before combining them at the location of the known MUSE position. The other involved using the correct MUSE redshifts, but extracting the spectra at random positions within the region where the LP primary beam response is > 50%, instead of at the known source position. The number of randomized extractions is the same as the number of CO emission samples. Neither of these methods should produce stacked spectra with real features. This random spectral extraction highlights which features in the real stacked spectra should be discounted as noise.

In Figure 2, we show the standard deviation of randomly stacked spectra against the number of the stacks.

For this plot, we randomize both redshift and position for the spectral extractions. The standard deviation of stacked spectra roughly decreases with 1/ √ N .

Sample selection

We took advantage of the significant redshift overlap between the ASPECS and MUSE surveys to carry out a stacking analysis on CO line emission: we focused on the ASPECS Band 3 data where the full MUSE redshift coverage can be exploited. In particular, over the range 1.0059 ≤ z ≤ 1.7387, the ASPECS survey has the highest sensitivity for detecting CO(2-1) among all of the observable CO emission features (Figure 9 of Boogaard et al. 2019). In addition, this redshift range is where MUSE is efficient at detecting spectral features (up to z ∼ 1.5 for [O II]) as shown in Figure 1.

We first selected a subset of the MUSE sources for stacking CO(2-1) emission in the ASPECS Band 3 data cube. The criteria were the following: (1) 1.0059 ≤ MUSE spec-z ≤ 1.7387, (2) MUSE spec-z confidence levels ≥ 2, and (3) location within ALMA Band 3 LP primary beam response ≥ 50%. For high-J CO emission, we used the same selection criteria, except for the redshift range. The ranges were z = 2.0088 -3.1080 for CO(3-2), z = 3.0115 -4.4771 for CO(4-3) z = 4.0142 -5.8460 for CO(5-4) z = 5.0166 -7.2146 for CO .

The resulting number of MUSE sources for stacking CO(2-1) was 111 in total4 . For 104 sources, the MUSE spectroscopic redshift identifications were based on [O II], five sources used absorption features, one source used C III], and one quasar had strong Mg II emission. Among these 111 sources, 10 sources5 were identified with CO(2-1) emission by the ASPECS blind search: MUSE IDs 8 (ASPECS-LP-3mm.06), 16 (3mm.11), 924 (3mm.14), 925 (3mm.16), 996 (3mm.02), 1001 (3mm.05), 1011 (3mm.10), 1117 (3mm.04), 6415 (3mm.08), and 6870 (3mm.15) (Boogaard et al. 2019;Aravena et al. 2019).

We performed the stacking in each group of the sample galaxies classified by the stellar mass and specific SFR (SSF R = SF R/M * ) derived from the MAGPHYS SED fits (see §2.3). The SF R -M * relation of our sam-ple is shown in Figure 3. The SSFR classification was based on the main-sequence (MS) relation from Eq.11 of Boogaard et al. (2018), using the mean redshift of the CO(2-1) line detectable in Band 3 (z = 1.43; Walter et al. 2016). Galaxies lying within, above, and below the intrinsic scatter (0.44 dex6 ) of this relation are referred to as the "MS", "above" the MS, and "below" the MS, respectively, throughout this paper. We here use the main-sequence relation from Boogaard et al. (2018) because their SF R -M * correlation is assessed with the objects whose spectral features and redshifts were measured based on the MUSE data. Compared with earlier studies such as Whitaker et al. (2014) and Schreiber et al. (2015) whose samples are mostly massive galaxies, Boogaard et al. (2018) better constrained the low-mass end of the relation. The numbers of galaxies in each group for stacking are presented in Table 1.

The main-sequence relation of Boogaard et al. (2018) was also adopted for classifying galaxies which expected to have CO J > 2 lines. Similar to CO(2-1), we used the mean redshifts of z = 2.61, 3.80, 4.99, and 6.18 for CO(3-2), CO(4-3), CO(5-4), and CO(6-5) lines, respectively (Walter et al. 2016).

CO line flux and upper limit measurements

For each stacked spectrum, we performed a best fit on the 2D image (moment-0) with a 2D Gaussian function to estimate the line flux or upper limit. The x-and y-positions are allowed to vary in the vicinity of the MUSE positions within a radius of 1/3 of the ALMA beam size. The widths were fixed to the beam size. We considered a CO emission line to have been detected when the amplitude of the fitted Gaussian function was higher than 3σ of the local fluctuations in the 2D image. The line flux was estimated by integrating the fitted Gaussian function. When the CO line was not detected, the standard deviation in the central area of the image was used to evaluate the 3σ upper limit, assuming a point source.

Uncertainties in MUSE redshifts can cause some flux losses when we stack spectra. We assumed an emission line with a width of 300 km s -1 (full width at half maximum, FWHM) and performed a bootstrap simulation to assess the flux loss resulting from σ z = 0.00012(1 + z) ( §2.2). The total flux was measured within ∆v = 528 km s -1 (corresponding to 19 slices in frequency; see §4. the rest-frequency of CO(2-1). This velocity range was the same as the one we used to create coadded 2D images for flux measurements (see §4). With 10000 realizations, we found that 92% of them resulted in less than 3% flux loss. We do not scale up our measured fluxes in this work to take account of this flux loss, because its impact on our line flux measurements is not significant.

CO EMISSION FROM STACKED SPECTRA

We performed CO emission line stacking from CO(2-1) to CO . In this section, we will present the resulting stacked spectra.

CO(2-1)

Stacking the entire sample

We first performed stacking without binning in stellar mass or SSFR to obtain a constraint on the average properties of all galaxies. The entire sample (111 spectra) and a subset of the sample (101 spectra) that excluded 10 objects whose CO emission were detected by the blind search ( §3.2) were used to search for signals in the stacks (below the noise threshold of individual galaxies).

For these two sets of samples, we inspected 2D images (moment-0) that were coadded between 230.335 ≤ ν rest /GHz ≤ 230.741 (∆v = 528 km s -1 ) centered at the CO(2-1) rest frequency (230.538 GHz) in the corresponding 1D spectrum. This width is consistent with 2σ of the mean line width of the ASPECS blind CO detections (308 km s -1 in FWHM; Aravena et al. 2019) to recover 95% of the emission.

The total line fluxes were obtained by fitting a 2D Gaussian function in the 2D images. This measure helps to avoid problems such as a possible slight positional offset between the optical and CO emission. All of the Gaussian parameters, the peak position, and amplitude were set to be free parameters, whereas the widths were fixed to the mean beam size. We obtained total line fluxes (upper limit) of 0.031 ± 0.007 Jy km s -1 (< 0.012 Jy km s -1 ) for the stacked spectra including (excluding) the sources with direct CO(2-1) detections.

Stacking in stellar mass and SSFR bins

We carried out the stacking on both the entire sample and the sample that excluded the direct CO detections. To investigate the dependence of the molecular gas content on fundamental properties of galaxies, we divided the sample into ranges of stellar mass and SSFR to perform the stacking. As shown in Figure 3, we set the bins to have steps of log (M * /M ) = 1.0 over the range 8.0 < log (M * /M ) ≤ 12.0 (4 bins) with stellar mass, and galaxies "above" the MS, on the "MS", and "below" the MS based on the main-sequence relation at z = 1.43 discussed in §3.2. In the redshift range where CO(2-1) can be detected with the ASPECS Band 3 survey (z = 1.0 -1.7), the main-sequence relation does not evolve significantly.

The stacked 2D and 1D spectra are shown in Figure 4 and the CO line flux measurements are summarized in Table 2. When the directly detected CO(2-1) emission is included, there are solid detections (> 3σ) at log (M * /M ) > 10.0, regardless of their SSFRs (except in the 10.0 < log (M * /M ) ≤ 11.0 bin where there is only a single source). On the other hand, at log (M * /M ) ≤ 10.0, we only obtain a signal in the bin above the MS with 9.0 < log (M * /M ) ≤ 10.0. When we exclude the known CO(2-1) emission from the stacks, a significant detection is seen in the MS bin with 10.0 < log (M * /M ) ≤ 11.0 (seven sources) 7 .

Although the bins of the MS galaxies in the ranges of 8.0 < log (M * /M ) ≤ 9.0 and 9.0 < log (M * /M ) ≤ 10.0 contain the largest numbers of galaxies, we do not detect a CO(2-1) line. To attempt to detect stacked CO emission for galaxies in these low mass bins, we adopt a ∆ log (M * /M ) = 0.5 stellar mass bin size in case some hidden emission from a small number of sources has been diluted by averaging too many sources without any emission. This smaller bin size results in 14, 24, 23, and 9 MS sources from 8.0 < log (M * /M ) ≤ 8.5 to 9.5 < log (M * /M ) ≤ 10.0 in steps of log (M * /M ) = 0.5. No detection is identified even with these finer bins. The estimated 3σ upper limits with the 2D data are 0.028, 0.024, and 0.017, 0.030 Jy km s -1 , from lower to higher stellar mass bins, respectively.

Among the bins with stacked detections, although the censoring fraction is unity for the bin of below the MS with stellar mass 10.0 < log (M * /M ) ≤ 11, all of the remaining bins are ≤ 0.75. In particular, the fraction is ≤ 0.6 for the MS galaxies.

The CO(2-1) lines detected in individual galaxies

We inspect individual spectra of galaxies in the MS bin of 10 < log (M * /M ) ≤ 11, where the stacked spectrum has the highest SNR (the spectra are depicted in Appendix A). Excluding the emission already found by the blind search, it is possible to identify the CO(2-1) line with lower SNR in MUSE IDs 879, 985, and 1308. We note that MUSE IDs 879 and 985 have already been reported by Boogaard et al. (2019) as the MUSE priorbased sample (see their Table 2 andFigure 4). A further potential stacked CO(2-1) detection is seen in the range of 10.0 < log (M * /M ) ≤ 11.0 for the mean of two galaxies lying below the MS. A detected CO(2-1) line is dominated by MUSE ID 928.

Higher-J CO stacked spectra

Along with the CO(2-1) stacking analysis, we attempt to detect higher excitation CO emission, up to J = 6 -5, with the same stacking method. The census of the stacked sources is shown in Table 3.

For CO(3-2), among our sample selection ( §3.2), MUSE IDs 35 and 11248 are the only sources that were already known from the blind search (ASPECS-LP-3mm.01 and 3mm.12). The former is a galaxy above the MS and the latter is on the MS with log (M * /M ) = 10.39 and 10.64, respectively. In addition, we consider MUSE ID 35's pair galaxy MUSE ID 24 (log (M * /M ) = 9.45) to have a CO(3-2) detection in our analysis. This is because the CO spatial extent, which peaks at the location of ID 35, covers ID 24 and it is difficult to distinguish the flux contributions from these two sources. When the stacking was performed in the same way as for CO(2-1), these two sources dominated the detections in their bins (Appendix B). If we discard them, no significant emission remains. This result agrees with the finding of Uzgil et al. (2019) who also did not find additional CO(3-2) emission with a masked auto-power spectrum using all MUSE positions with LP primary beam response ≥ 20%. There is also no stack detection for J up ≥ 4 CO emission in our sample with the blind search9 , nor the stacking.

We note that the MUSE spectroscopic redshifts of galaxies at z 3, which allow detections of J up > 4 CO emission, were mostly determined using the Lyα line. Its redshift was measured using the peak emission, which can have a few hundred km s -1 offset from the systemic redshift (e.g., Shapley et al. 2003). When the redshifts of the stacked sources were measured with Lyα, we tried to recover their systemic redshifts based on an empirical correlation between the Lyα emission peak and the Lyα line width (Verhamme et al. 2018). The intrinsic scatter of the relation is ±73 km s -1 measured with 13 sources that have both Lyα and C III] detections. We do not find any detection for J up > 4 CO emission either. This is likely either due to dilution of the stacked signal, owing to uncertainties in the velocity offset, or to Lyα emitters on average having smaller SFRs and molecular gas content.

We also visually investigate individual CO spectra with the MUSE redshifts as priors to search for possible CO emission which is washed out by stacking. Two potential CO(5-4) emission lines were identified as shown in Figure 10 in Appendix C. If confirmed, they may be one of the first cases of high-J CO detection at z > 4.5 for Lyα emitters. Deeper observations are needed to confirm these detections.

DISCUSSION

In this section, we will concentrate on discussing the molecular gas content at z ∼ 1.5, as we have obtained the most reliable measurements at this redshift with our analysis.

CO Luminosities and Molecular (H 2 ) Gas Masses

We employed the following equation to obtain the CO luminosities from the measured CO(2-1) emission (Solomon et al. 1997;Solomon & Vanden Bout 2005):

L COJ-[J-1] K km s -1 pc 2 = 3.25 × 10 7 F line Jy km s -1 D 2 L (1 + z) 3 ν 2 obs (1)
where F line is the integrated line flux density, D L is the luminosity distance in Mpc, and ν obs is the observed frequency in GHz. For the redshift (z), we used the mean redshift of the galaxies in each stacking bin.

We then adopted the following equation and the same conditions as Decarli et al. (2016b) to infer molecular gas masses (M ) from the line luminosities (K km s -1 pc 2 ):

M gas = α CO r J1 L COJ-[J-1] (2)
where J is the upper level of the CO excitation and α CO is the CO luminosity to gas mass conversion factor. We assume r 21 = 0.76 ± 0.09 following Daddi et al. (2015). Based on galaxies detected with the AS-PECS survey, Boogaard et al. (2020 submitted) estimated r 21 = 0.75 ± 0.11, which is comparable to the value from Daddi et al. (2015). Here we adopt the value from the former to be consistent with the molecular mass measurements published in Aravena et al. (2019), which are used for comparison in this paper. For α CO , we used 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the stellar mass bins of log (M * /M ) > 10 which is one of the best estimates for high-redshift main-sequence star-forming galaxies (Daddi et al. 2010). Among the ASPECS directly detected CO sources with metallicity estimates available from the MUSE spectra, all have roughly solar metallicity (Boogaard et al. 2019), which further justifies our choice of α CO (Aravena et al. 2019). For the galaxies with log (M * /M ) ≤ 10.0, because they are likely to have lower metallicity, α CO could be higher (e.g., Bolatto et al. 2013). In this work, we assume that the metallicities of galaxies with log (M * /M ) ∼ 9.5 and ∼ 8.5 are ∼ 2/3Z and ∼ 1/2Z (e.g., Savaglio et al. 2005;Erb et al. 2006;Mannucci et al. 2009;Zahid et al. 2011;Sargent et al. 2014;Sanders et al. 2020). Thus, α CO would increase by a factor of ∼ 2 and ∼ 3, respectively, compared to galaxies with log (M * /M ) > 10.0 (e.g., Genzel et al. 2012;Schruba et al. 2012). We adopt α CO = 7.2 and 10.8 M (K km s -1 pc 2 ) -1 for galaxies with log (M * /M ) ∼ 9.5 and log (M * /M ) ∼ 8.5, respectively. In these lower stellar mass bins, the M gas estimates with a constant α CO = 3.6 M (K km s -1 pc 2 ) -1 are also shown as dotted symbols in figures for reference. The calculated line luminosities and molecular gas masses from CO(2-1) are summarized in Table 4.

In the following subsections, we discuss the gas content of galaxies at z ∼ 1.5, together with their stellar masses and SFRs derived from the SED fitting ( §2.3). The CO(2-1) emission used here stems from the measurements based on the 2D Gaussian fits in the images presented in §4.1.2 and Figure 4. Because we are interested in global properties of the molecular gas, we use the measurements that include the direct CO(2-1) detections from the blind search.

Molecular gas and stellar mass scaling relation of MS galaxies

Previous observational and theoretical studies of molecular gas in high redshift galaxies have established a set of scaling relations that relate the molecular gas content to galaxy properties such as stellar masses, SFRs, source sizes (e.g., Tacconi et al. 2013Tacconi et al. , 2018;;Freundlich et al. 2019;Aravena et al. 2019;Genzel et al. 2015;Magdis et al. 2012;Santini et al. 2014;Scoville et al. 2016;Somerville et al. 2012;Davé et al. 2017;Popping et al. 2019). These relations have been key to understanding how the galaxy growth process has taken place through cosmic time. We first compare average molecular gas properties to the stellar mass of distant galaxies.

In conjunction with the MUSE deep survey in the same field, the stacking analysis facilitates the exploration of the gas scaling relations below stellar masses of log (M * /M ) = 10, a regime that is rather uncharted in CO emission for individual detections at z > 1 (see also Appendix D). We focus on discussing the mainsequence galaxies, which have the largest numbers of stacked spectra in our CO(2-1) sample.

In Figure 5 (left), we plot the molecular gas mass from the stacking analysis against stellar mass of the mainsequence galaxies, along with individual ASPECS CO detections (González-López et al. 2019;Aravena et al. 2019). The majority of directly detected galaxies have log (M * /M ) > 10. At log (M * /M ) ≤ 10, there is only one direct detection. Although stacking also leads to no detection, the stacks of 38 and 32 spectra provide tight constraints on the average gas mass in the range of 8 < log (M * /M ) ≤ 10 at z ∼ 1.5. Taking upper limits into account, an increase of gas mass is found with increasing stellar mass from at least log (M * /M ) ∼ 9.0 to 12.0.

For comparison, we also show the contours of the molecular gas measurements from the PHIBSS2 survey (included both the CO-and dust-based measurements; Tacconi et al. 2013Tacconi et al. , 2018;;Freundlich et al. 2019). For further comparison with the literature, we show the distributions of the PHIBSS2 galaxies at lower redshifts in Appendix D. The gas masses of the PHIBSS2 sources were derived via a metallicity-based prescription for this parameter (Genzel et al. 2012), resulting in α CO ∼ 4 -7 M (K km s -1 pc 2 ) -1 at z ∼ 1.5.

We also compare our measurements with cosmological galaxy formation model calculations of molecular gas mass presented in Popping et al. (2019). In the same diagram we show three different model predictions: the IllustrisTNG hydrodynamical simulations with the "3.5arcsec" and "Grav" apertures and the Santa Cruz semi-analytic model (SC SAM). The former corresponds to all the H 2 within a radius of 3.5 of the source center (similar to ASPECS) and the latter corresponds to all the H 2 gravitationally bound to the galaxy. The H 2 properties of galaxies were derived based on the molecular hydrogen fraction recipe of Gnedin & Kravtsov (2011) 10 .

We use the gas mass model predictions at z = 1.43 (Popping et al. 2019), which is the mean redshift of the CO(2-1) line detectable in Band 3. At the higher stellar mass end, all of the three models predicted lower gas mass than the observed gas mass. At log (M * /M ) ∼ 10.5, about half of the sources with ASPECS direct detections lie within the 2σ scatter of the predictions. At log (M * /M ) 10.0, where the predicted gas mass is below the ASPECS detection limit, the only CO constraint is from stacking. The upper limit at log (M * /M ) ∼ 9.5 derived from the stacking is consistent with the models. This result also implies that at least a factor of 10 increase in the sample size is needed to confirm or rule out the model predictions.

Molecular gas-to-stellar mass ratio of MS galaxies as a function of stellar mass

We depict this plot with a different presentation in Figure 5 (right) to show a more common presentation of molecular gas mass normalized by the stellar mass (molecular gas-to-stellar mass ratio, µ gas ≡ M gas /M * ) as a function of stellar mass. In the literature, a decrease of the gas-to-stellar mass ratio with increasing stellar mass at log (M * /M ) > 10.5 has been reported 10 These models are the same as the ones shown in Figure 2 in Popping et al. (2019), but without the ASPECS observational selection effects (i.e. we show here the entire population of galaxies predicted by the models). [Left] Molecular gas mass as a function of stellar mass of the main-sequence galaxies in this study. We use αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively, to obtain molecular gas mass (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The large green circles are from the stacked CO detections and the downward-pointing arrows are the 3σ upper limits. The squares indicate the ASPECS CO detections from the blind search (Aravena et al. 2019). The PHIBSS2 galaxies (Tacconi et al. 2013(Tacconi et al. , 2018;;Freundlich et al. 2019), classified as MS galaxies with the MS relation of Boogaard et al. (2018) over 1.00 ≤ z ≤ 1.74, are shown as the gray contours for comparison to our results. We also display the model predictions of molecular gas mass from SC SAM (pink line and shading, which indicate the median and 2σ scatter, respectively), IllustrisTNG with the 3.5arcsec aperture (blue), and IllustrisTNG with the Grav aperture (beige) from Popping et al. (2019). [Right] Molecular gas-to-stellar mass ratio is presented instead as a function of stellar mass. The symbols are the same as in the left panel. The black line with dots shows the running median of the gas-to-stellar mass ratios of the PHIBSS2 sources lying on the MS of Boogaard et al. (2018) (gray contours; see also Tacconi et al. 2018).

at z ∼ 1.5 (e.g., Tacconi et al. 2013). However, no observational CO constraints exist at lower stellar masses below 10 10 M . Our stacking analysis facilitates exploration of this lower stellar mass regime.

For galaxies with log (M * /M ) ∼ 10.5 and 11.5, we estimate µ gas = 0.49 ± 0.08 and 0.24 ± 0.02, respectively. Similarly to Figure 5 (left), these values are in agreement with the measurements of the PHIBSS2 galaxies lying on the MS relation of Boogaard et al. (2018) (grey contours; see Appendix D for a version using the MS relation of Speagle et al. 2014). In the lower stellar mass bin, log (M * /M ) ∼ 9.5, we constrain µ gas to have a 3σ upper limit of < 1.25 with α CO = 7.2 M (K km s -1 pc 2 ) -1 (or < 0.63 if α CO = 3.6 M (K km s -1 pc 2 ) -1 ).

Beyond log (M * /M ) ∼ 10.5, a decrease of the gasto-stellar mass ratio is discerned with increasing stellar mass in our stack results, as well as in previous work reporting both CO-based and dust-based gas estimates (Magdis et al. 2012;Tacconi et al. 2013Tacconi et al. , 2018;;Genzel et al. 2015;Scoville et al. 2017;Aravena et al. 2019;Liu et al. 2019). As an example, we show in Figure 5 the result from PHIBSS2 with a steady decrease of µ gas with an increase of stellar mass (black line).

This decline of µ gas does not seem to be as steep at log (M * /M ) < 10.0. If a constant α CO is also adopted for these lower mass bins, the gas-to-stellar mass ratio is consistent with being constant from log (M * /M ) ∼ 9.5 to 10.5, then turning down around log (M * /M ) ∼ 10.5. We note that the models from Popping et al. (2019), despite a discrepancy with the observed results at the high mass side, show a similar constant gas-to-stellar mass ratio up to log (M * /M ) ∼ 10.5. This is consistent with our finding if we assume a constant conversion factor.

Given the limitations of our data, and the unknown dependence of the α CO conversion factor, we cannot determine at which stellar mass a possible downturn of the gas-to-stellar mass ratio occurs. We, however, note that such a "plateau" in the low stellar mass regime has also been identified for local galaxies with direct detections of CO emission: here µ gas stays constant from at least log (M * /M ) ∼ 9.0 and starts decreasing around log (M * /M ) ∼ 10.5 (e.g., Saintonge et al. 2017;Bothwell et al. 2014). A similar result is found in Tacconi et al. (2018) who compared local star-forming galaxies, where CO measurements for low mass galaxies are available (including Saintonge et al. 2017), to distant galaxies after removing assumed redshift effects on their gas mass content (see also Appendix D).

The declining µ gas at the high stellar mass end can be attributed to stellar feedback (e.g., Tacconi et al. 2013;Genzel et al. 2015;Davé et al. 2011). In contrast, the flatter µ gas at log (M * /M ) 10.5 may imply that the effects of feedback are weaker. The drop in gas-to-stellar mass ratio seems to appear around log (M * /M ) ∼ 10.5, which is relevant to the characteristic stellar mass of star-forming galaxies (e.g., Duncan et al. 2014) where mass-quenching is becoming dominant (Peng et al. 2010).

Furthermore, the gas-to-stellar mass ratio below ∼ 1 at 3σ of low stellar mass galaxies suggests that we may have retrieved most of the CO emission at z ∼ 1.5. As shown in Decarli et al. (2019), CO luminosity functions derived from the ASPECS data are assumed to have a fixed faint-end slope, which is consistent with the faintend slope of the stellar mass function of star-forming galaxies. Thus, the value of µ gas 1 indicates that the assumed faint-end slope of the CO luminosity function is at least consistent or could be flatter than that of the stellar mass function. If this is the case, most of the CO emission at z ∼ 1.5 has been recovered (see also Uzgil et al. 2019).

The estimated molecular gas density ρ(H 2 ) at z ∼ 1.5, based on our CO(2-1) stacking measurements, is (0.49 ± 0.09) × 10 8 M Mpc -3 . This value is lower than, but formally consistent with, the value derived in Decarli et al. (2019) with a different CO selection. The CO emission used in Decarli et al. (2019) included sources that did not enter the present analysis because the lack of a counterpart with high-quality MUSE redshift. Our total CO flux estimate from stacking is also in line with the result from the CO auto-power spectrum analysis using the MUSE positions (Uzgil et al. 2019).

Dependence of the molecular gas content on SSFR

In the previous section we only considered mainsequence galaxies; in the remaining discussion, we will include galaxies above and below the main-sequence relation of Boogaard et al. (2018) to discuss the dependence of molecular gas content on SSFR. Below, we keep the same assumed α CO conversion factors as above: 3.6, 7.2, and 10.8 M (K km s -1 pc 2 ) -1 for galaxies with log (M * /M ) > 10, ∼ 9.5, and 8.5, respectively ( §5.1).

We include the constraints using a constant α CO = 3.6 throughout the stellar mass bins in the figures for reference.

Gas-to-stellar mass ratios

As shown earlier, the gas-to-stellar mass ratio (µ gas ) provides information on the supply and depletion of gas reservoirs in galaxies. We compare µ gas against SSFRs in the left panel of Figure 6. We find that µ gas increases with increasing SSFRs, which is also seen in earlier studies of high-redshift galaxies, including the dust-based measurements of the molecular gas mass (e.g., Tacconi et al. 2013Tacconi et al. , 2018;;Genzel et al. 2015;Magdis et al. 2012;Scoville et al. 2016;Liu et al. 2019). With the stacking analysis, we show that this trend holds even for galaxies below the main-sequence at z ∼ 1.5.

In the right panel of Figure 6, we show µ gas of stacked sources color-coded by different SSFRs to compare with their evolution (Geach et al. 2011;Magdis et al. 2012). Because our data points are derived from the CO(2-1) line, all of them lie between z = 1.0-1.5. A wide spread is related to the variations of µ gas across the mainsequence relation seen in the left panel (see also e.g., Tacconi et al. 2018). The lower µ gas of low SSFR sources is also found by Spilker et al. (2018) who investigated galaxies lying below the MS with log (M gas /M ) ∼ 11 at z ∼ 0.7. The depletion times of our sample below the MS are on average comparable with Spilker et al.'s sample. Following the same assumption as Spilker et al. (2018), if these galaxies continuously consume the existing gas with the observed SFR, then the µ gas of the galaxies with log (M gas /M ) ∼ 10.5 (11.5) would reduce to 1/10 at z ≈ 0.5 (0.2) and 1/100 at z ≈ 0.1 (0.0). Hence, by z = 0 their µ gas would be comparable to those of passive galaxies at the current epoch.

Gas depletion time

The gas depletion time, t depl = M gas /SF R, is another way of examining the gas content in galaxies. It estimates the time taken for the gas to be fully consumed at the current SFR without accounting for additional fueling. We show the gas depletion time as a function of SSFR in Figure 7. Our results from the stacking analysis are in good agreement with previous studies. Both the stacked and directly detected sources show decreasing t depl with increasing SSFR, except that the 9.0 < log (M * /M ) ≤ 10.0 bin above the MS shows an elevated t depl value. The data points from these two sets of measurements occupy the same range in the diagram, following the constant gas-to-stellar mass ratio (M gas /M * ) from around 0.1 to 1.0. The decreasing gas depletion timescale with SSFR demonstrates that galaxies with more extreme star formation consume their gas [Left] Gas-to-stellar mass ratios as a function of SSFR. The symbols are the same as in Figure 5, but now we show the entire sample, below, on, and above the main-sequence relation, color-coded with their stellar masses as indicated in the legend. The same as the other figures involving molecular gas mass estimates, we use αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The data points are centered on the mean values of the parameters on the x-axis. The three diagonal gray dotted lines represent constant depletion time scales of 0.1, 1.0, and 10 Gyr from the bottom to the top. The background pink, green, and blue filled colors indicate the below, on, and above the main-sequence relation of Boogaard et al. (2018) for a galaxy at log (M * /M ) = 10.0. A version using the main-sequence relation of Speagle et al. (2014) is shown in Appendix D. [Right] The evolution of gas-to-stellar mass ratios. The symbols are the same as Figure 5 but now we include galaxies above and below the main-sequence relation: above (blue), below (red), and on the MS relation (green). The values inside each symbol indicate the corresponding stellar mass bins in log scale. The PHIBSS2 galaxies classified as above, below, and on the MS by the main-sequence relation of Boogaard et al. (2018) are also depicted as red, green, and blue contours, respectively. The gray curve in the background and its shaded region show the evolutionary track of main-sequence galaxies from Tacconi et al. (2018).

at a higher rate. This may be the major cause of the scatter as shown in Figure 6.

SUMMARY AND CONCLUSIONS

Based on the accurate redshifts from the MUSE IFU survey in the Hubble Ultra Deep Field, we perform a CO emission stacking analysis with the ASPECS Band 3 data. When we split the sample into stellar mass and SSFR bins (on, below, and above the main-sequence relation), we detect CO(2-1) emission (z ∼ 1.43) down to log (M * /M ) = 10.0, even after removing previously reported CO detections. We do not recover any higher-J CO emission at higher redshift (z 3) with stacking when excluding the sources with direct CO detections.

The 3σ upper limits on CO(2-1) emission at log (M * /M ) < 10.0 in the redshift range 1.0 < z < 1.7 provide meaningful upper limits on molecular gas mass estimates of main-sequence star-forming galaxies in this stellar mass range, which has been poorly explored at z ∼ 1.5. Under the assumption of a metallicity-based α CO conversion factor, we observe an increase in gas mass with increasing stellar mass from log (M * /M ) ∼ 9.0 to 11.0. The upper limits at the low mass end are consistent with the model predictions, but to confirm or Depletion time (t depl ) plotted against SSFR. The symbols are the same as Figure 6 (left). The same as the earlier figures, the molecular gas masses are estimated based on αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The colors of the filled symbols, shown in the legend, indicate each mass bin used for the stacking. The three diagonal gray dotted lines represent constant gas-to-stellar mass ratios of 0.1, 1.0, and 10 from the bottom to the top.

rule out these predictions, at least a factor of 10 increase in the sample size is needed. The gas-to-stellar mass ratio (µ gas ≡ M gas /M * ) from log (M * /M ) ∼ 9.0 to 10.5 declines at a slower rate compared with the known decrease of µ gas at higher stellar masses. If a fixed α CO = 3.6 M (K km s -1 pc 2 ) -1 conversion factor is assumed across the stellar mass range explored here, µ gas is consistent with being constant from at least log (M * /M ) ∼ 9.0 to 10.5, in agreement with predictions by models and observations at lower redshifts.

Furthermore, the gas-to-stellar mass ratios of ∼ 0.5 and 1 at the stellar masses of log (M * /M ) ∼ 10.5 and log (M * /M ) ∼ 9.5, respectively, imply that the faint end slope of the CO luminosity function is at least consistent or could be flatter than the stellar mass function. We have successfully recovered the majority of the CO emission at z ∼ 1.5 with the stacking analysis. The molecular gas density ρ(H 2 ) = (0.49 ± 0.09) × 10 8 M Mpc -3 from this stacking analysis is comparable with the one inferred from a CO-driven selection (Decarli et al. 2019).

When we compare the gas-to-stellar mass ratio (µ gas ) against SSFR, we confirmed the known correlations to also hold for galaxies with low SSFRs at z ∼ 1.2. The scatter in the µ gas -SSFR correlation seems to be related to the decrease of µ gas at log (M * /M ) 10.5 for starforming galaxies. We also show that the gas-to-stellar mass ratios of massive galaxies (log (M * /M ) ∼ 11) below the MS at z ∼ 1.2 are comparable to those at z ∼ 0.7 (Spilker et al. 2018;Tacconi et al. 2018).

Our stacking analysis of the combined volumetric surveys of ALMA and MUSE has let us explore regimes that were uncharted before and that will remain challenging for investigations that rely on direct CO detections.

The authors would like to thank the referee whose constructive comments helped improving the manuscript. Note-The units are in 10 9 K km s -1 pc 2 and 10 9 M for (1) CO(2-1) line luminosities (the first row in each MS bin) and (2) molecular masses (Mgas; the second row), respectively. The CO-to-H2 conversion factors are assumed to be αCO = 3.6, 7.2, and 10.8 M (K km s -1 pc 2 ) -1 for galaxies in the bins of log (M * /M ) > 10, 9 < log (M * /M ) ≤ 10, and 8 < log (M * /M ) ≤ 9, respectively (see §5.1). The values in parentheses are after excluding the galaxies which have the CO(2-1) line identified by the blind search (see §3.2).

Table 5. On the MS (1) < 14.12 (< 14.12) < 1.25 (< 1.25) 0.49 ± 0.08 (0.46 ± 0.12) 0.24 ± 0.02 (-) (2) < 11.27 (< 11.27 Note-The first and second rows in each MS bin are (1) molecular-to-stellar mass ratios (Mgas/M * ) and (2) depletion times (Mgas/SF R in Gyr), respectively. The CO-to-H2 conversion factors are assumed to be αCO = 3.6, 7.2, and 10.8 M (K km s -1 pc 2 ) -1 for galaxies in the bins of log (M * /M ) > 10, 9 < log (M * /M ) ≤ 10, and 8 < log (M * /M ) ≤ 9, respectively (see §5.1). The values in parentheses are after excluding the galaxies which have the CO(2-1) line identified by the blind search (see §3.2). 5, to show our results in the context of previous studies of the molecular gas content. Note that the MS galaxies here are classified based on the prescription of Speagle et al. (2014) following Tacconi et al. (2018). Thus, the data points and the MS boundaries are slightly shifted compared to Figure 5 in the main text. The distributions of PHIBSS2 galaxies at lower redshifts are displayed in the left panels (color-coded by redshift as given in the legend). Most of the PHIBSS2 galaxies at z ∼ 0 are from xCOLD GASS (Saintonge et al. 2017). 2014) is used here instead. Thus, the data points and the MS boundaries are slightly shifted compared to Figures 6 and7. The light pink, green, and blue background colors indicate the regions of below, on, and above the main-sequence relation at log (M * /M ) = 10, respectively. The orange dashed lines are the best-fit lines for the PHIBSS2 sample (Tacconi et al. 2018).

Figure 1 .

 1 Figure1. The redshift coverages of the MUSE (blue) and ASPECS (green) HUDF surveys. The numbers of ASPECS directly detected sources and MUSE sources with a secure redshift measurement are listed on the right hand side (There are also nine and two sources classified as "stars" and "others", respectively,Inami et al. 2017). The numbers in parentheses are the MUSE sources whose locations are in the region of LP primary beam response > 50% for each CO emission. The absorption features detected with MUSE include C IV, Fe II, Mg II. The analysis presented in this work focuses on the stacking frequencies that correspond to the CO(2-1), CO(3-2), CO(4-3), CO(5-4), and CO(6-5) lines.

Figure 2 .

 2 Figure2. Reduction of standard deviation of stacked spectra (100 realizations) with increasing number of stacked spectra (red region). The stacked spectra used here were extracted randomly from the ASPECS data cube. The red dotted line represents 1/ √ N , where N is the number of the stacked spectra.

  1.1 for this choice of ∆v) centered at

Figure 3 .

 3 Figure 3. SFR-M * relation of our CO(2-1) stacking sample. The gray filled band shows the main-sequence relation and its intrinsic scatter for MUSE-detected galaxies (at z = 1.43, Boogaard et al. 2018). The gray solid, dashed and dashed dotted lines are the main-sequence relations from Speagle et al. (2014) (which is used for PHIBSS2 Tacconi et al. 2018, see also Appendix D), Whitaker et al. (2014) andSchreiber et al. (2015), respectively. The blue vertical and diagonal dotted lines show the divisions of the bins that we used for stacking. The numbers at the top left of each grid are the counts of galaxies in the bin. The filled circles indicate the galaxies which have CO(2-1) line detections in the ASPECS blind search (see §3.2).

Figure 4 .

 4 Figure 4. Resulting stacked 2D images (left) and 1D spectra (right) in the bins of SSFRs (above, on, and below the galaxy main-sequence, as discussed in §3.2; y-axis) and stellar masses (x-axis). The upper and lower two panels show the results of including and excluding the individually reported direct CO(2-1) detections, respectively (Aravena et al. 2019; Boogaard et al. 2019). The 2D images are coadded between 230.335 ≤ νrest/GHz ≤ 230.741 (∆v = 528 km s -1 ) and the white ellipses indicate the mean of the beam size. The 1D spectrum is extracted from the central spaxel (the red solid line). Spectra extracted with random redshifts in place of the MUSE position and with random positions at the correct MUSE redshift are shown by the light blue and green lines, respectively ( §3.1).

  Figure 5.[Left] Molecular gas mass as a function of stellar mass of the main-sequence galaxies in this study. We use αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively, to obtain molecular gas mass (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The large green circles are from the stacked CO detections and the downward-pointing arrows are the 3σ upper limits. The squares indicate the ASPECS CO detections from the blind search(Aravena et al. 2019). The PHIBSS2 galaxies(Tacconi et al. 2013(Tacconi et al. , 2018;;Freundlich et al. 2019), classified as MS galaxies with the MS relation ofBoogaard et al. (2018) over 1.00 ≤ z ≤ 1.74, are shown as the gray contours for comparison to our results. We also display the model predictions of molecular gas mass from SC SAM (pink line and shading, which indicate the median and 2σ scatter, respectively), IllustrisTNG with the 3.5arcsec aperture (blue), and IllustrisTNG with the Grav aperture (beige) fromPopping et al. (2019). [Right] Molecular gas-to-stellar mass ratio is presented instead as a function of stellar mass. The symbols are the same as in the left panel. The black line with dots shows the running median of the gas-to-stellar mass ratios of the PHIBSS2 sources lying on the MS ofBoogaard et al. (2018) (gray contours; see alsoTacconi et al. 2018).

Figure 6 .

 6 Figure 6. [Left] Gas-to-stellar mass ratios as a function of SSFR. The symbols are the same as in Figure5, but now we show the entire sample, below, on, and above the main-sequence relation, color-coded with their stellar masses as indicated in the legend. The same as the other figures involving molecular gas mass estimates, we use αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The data points are centered on the mean values of the parameters on the x-axis. The three diagonal gray dotted lines represent constant depletion time scales of 0.1, 1.0, and 10 Gyr from the bottom to the top. The background pink, green, and blue filled colors indicate the below, on, and above the main-sequence relation ofBoogaard et al. (2018) for a galaxy at log (M * /M ) = 10.0. A version using the main-sequence relation ofSpeagle et al. (2014) is shown in Appendix D.[Right] The evolution of gas-to-stellar mass ratios. The symbols are the same as Figure5but now we include galaxies above and below the main-sequence relation: above (blue), below (red), and on the MS relation (green). The values inside each symbol indicate the corresponding stellar mass bins in log scale. The PHIBSS2 galaxies classified as above, below, and on the MS by the main-sequence relation ofBoogaard et al. (2018) are also depicted as red, green, and blue contours, respectively. The gray curve in the background and its shaded region show the evolutionary track of main-sequence galaxies fromTacconi et al. (2018).

  Figure 7.Depletion time (t depl ) plotted against SSFR. The symbols are the same as Figure6(left). The same as the earlier figures, the molecular gas masses are estimated based on αCO = 3.6 M (K km s -1 pc 2 ) -1 for galaxies in the log (M * /M ) > 10 bins and 7.2 and 10.8 for the log (M * /M ) ∼ 9.5 and 8.5 bins, respectively (see §5.1). We also show the estimates based on a constant αCO = 3.6 M (K km s -1 pc 2 ) -1 in the lower mass bins with the dashed symbols for reference. The colors of the filled symbols, shown in the legend, indicate each mass bin used for the stacking. The three diagonal gray dotted lines represent constant gas-to-stellar mass ratios of 0.1, 1.0, and 10 from the bottom to the top.

  Summary of molecular-to-stellar mass ratios (µgas = Mgas/M * ) and depletion times (t depl = Mgas/SF R)

Figure 8 .

 8 Figure 8. Individual spectra used for stacking in the bins of 10.0 < log (M * /M ) ≤ 11.0 for the galaxies on the MS galaxies (left) and below the MS (right). The 2D image (moment-0) and 1D spectrum are displayed for each galaxy with the MUSE ID number at the top left in each panel of the 1D spectrum. If the emission has been detected by the blind (González-López et al. 2019) or prior searches (Boogaard et al. 2019), then its ASPECS ID is also shown in parentheses. The lines and symbols are the same as in Figure 4.

Figure 9 .

 9 Figure 9. The stacked results of CO(3-2), CO(4-3), CO(5-4), and CO(6-5), shown in the same manner as the left panel in Figure 4, at top left, top right, bottom left, and bottom right, respectively. No significant detections in the stacks are reported.

Figure 10 .

 10 Figure 10. Tentative detections of CO(5-4) of two Lyα emitters. The left panels show the 2D images of CO(5-4) created over the velocity range highlighted in the middle panels. The middle panels depict the ALMA 1D spectra. In the right panels, the MUSE 1D spectra are presented along with the MUSE redshift in the top right corner. The velocity is defined relative to the MUSE redshift with an empirical correction for Lyα systemic redshift (see §4.2; Verhamme et al. 2018).

Figure 11 .

 11 Figure11. The same diagrams as Figure5, to show our results in the context of previous studies of the molecular gas content. Note that the MS galaxies here are classified based on the prescription ofSpeagle et al. (2014) followingTacconi et al. (2018). Thus, the data points and the MS boundaries are slightly shifted compared to Figure5in the main text. The distributions of PHIBSS2 galaxies at lower redshifts are displayed in the left panels (color-coded by redshift as given in the legend). Most of the PHIBSS2 galaxies at z ∼ 0 are from xCOLD GASS(Saintonge et al. 2017).

Figure 12 .

 12 Figure 12. The same diagrams as Figures 6 and 7 in the main text, but lower redshift bins are also shown, following Figure 11. Again, the main-sequence relation of Speagle et al. (2014) is used here instead. Thus, the data points and the MS boundaries are slightly shifted compared to Figures6 and 7. The light pink, green, and blue background colors indicate the regions of below, on, and above the main-sequence relation at log (M * /M ) = 10, respectively. The orange dashed lines are the best-fit lines for the PHIBSS2 sample(Tacconi et al. 2018).

Table 1 .

 1 The authors thank Ian Smail for useful suggestions to improve this study. This work was supported by JSPS KAKENHI Grant Number JP19K23462 (HI). Numbers of galaxies in the stellar mass and SSFR bins for

	Este trabajo cont ó con el apoyo de CONICYT + PCI
	+ INSTITUTO MAX PLANCK DE ASTRONOMIA
	MPG190030. FW acknowledges support by ERC Ad-
	vanced Grant 740246 (Cosmic Gas). T.D-S. acknowl-
	edges support from the CASSACA and CONICYT
	fund CAS-CONICYT Call 2018. D.R. acknowledges
	support from the National Science Foundation under
	grant numbers AST-1614213 and AST-1910107 and
	from the Alexander von Humboldt Foundation through
	a Humboldt Research Fellowship for Experienced Re-
	searchers. MK acknowledges support from the Interna-
	tional Max Planck Research School for Astronomy and
	Cosmic Physics at Heidelberg University (IMPRS-HD).
	ALMA is a partnership of ESO (representing its member
	states), NSF (USA) and NINS (Japan), together with
	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
	ESO, AUI/NRAO and NAOJ.
	Facilities: ALMA, VLT:Yepun
	Software: astropy (Astropy Collaboration et al. 2013;
	The Astropy Collaboration et al. 2018), CASA (McMullin
	et al. 2007)

Table 4 .

 4 Summary of CO(2-1) line luminosities (L CO2-1 ) and molecular gas masses (Mgas)

	log (M * /M )

The ASPECS Band 6 (1mm) survey covers most of these gaps by observing higher-J CO lines.

In cases where MUSE could not spatially resolve the HSTdetected sources, these sources were "merged" into a single MUSE object. Its new coordinates are the HST F775W flux-weighted center of all the merged objects.

The median stacking also produced similar results, but here we adopt the mean stacking for simpler noise estimates and for a better treatment when there are only two samples to stack.

There are three galaxies(MUSE IDs 6314, 6450, and 6530) which meet the selection criteria, but their stellar masses and SFRs cannot be constrained because they have no optical counterpart on which to perform an SED fit. These galaxies are found by MUSE (no HST counterpart in the UVUDF catalog due to blending). These objects are not included in our sample of 111.

5 Four of them detected by the CO blind search, ASPECS-LP-3mm.02, 04, 05, and 08, have a MUSE redshift confidence level of 1, but they are included in this work.

This intrinsic scatter is found by assuming a Gaussian function in a model of the star formation sequence. As noted inBoogaard et al. (2018), this value is higher than the average value reported in previous work.

Note that a detection in the bin below the MS with 10.0 < log (M * /M ) ≤ 11.0 has already been observed when the known CO(2-1) emission were included, because none of the galaxies in this bin include the directly detected CO(2-1) emission (none of the galaxies have been excluded to make any changes).

This source has a new additional MUSE redshift z = 2.5739 whose foreground galaxy is identified at z = 1.098(Boogaard et al. 2019).

There is one CO(4-3) detection in the blind search, but this source does not have a MUSE redshift(Aravena et al. 2019;Boogaard et al. 2019).

Table 2. Summary of measured stacked CO(2-1) line fluxes (F line ) log (M * /M ) 8.0 -9.0 9.0 -10.0 10.0 -11.0 11.0 -12.0 Above the MS < 0.032 (< 0.032) 0.110 ± 0.027 (< 0.041) < 0.081 (< 0.081) -(-) On the MS < 0.022 (< 0.022) < 0.019 (< 0.019) 0.130 ± 0.021 (0.095 ± 0.024) 0.453 ± 0.042 (-) Below the MS -(-) < 0.051 (< 0.051) 0.119 ± 0.039 (0.119 ± 0.039) 0.152 ± 0.027 (< 0.054)

Note-The units are in Jy km s -1 . Values in parentheses are after excluding the galaxies which have the CO(2-1) line identified by the blind search (see §3.2).

Table 3. Total number of galaxies used for stacking high-J CO emission log (M * /M ) CO 7.0 -8.0 8.0 -9.0 9.0 -10.0 10.0 -11.0 11.0 -12.0

Above the MS 3-2 1 (1) 5 ( 5) 5 ( 4) 1 (-) -4-3 4 (4) 4 (4) 5 ( 5) 16) 14 ( 14) 7 ( 7)

Note-The numbers in parentheses are after excluding the galaxies which have the CO line identified by the blind search (see §3.2).

APPENDIX

A. THE INDIVIDUAL SPECTRA IN THE STACKING

The individual spectra for the stacking of the galaxies on and below the MS with 10.0 < log (M * /M ) ≤ 11.0 are displayed in Figure 8 ( §4.1.3).

B. THE STACKED IMAGES, SPECTRA, AND MEASURED UPPER LIMITS FOR HIGH-J CO LINES

The stacked images and spectra for CO(3-2), CO(4-3), CO(5-4), and CO are shown in Figure 9. The measured line flux upper limits are summarized in Table 6.

C. TENTATIVE DETECTIONS OF CO(5-4)

In Figure 10, we present the two cases of potential CO(5-4) detections. As discussed in § 4.2, further observations are required to verify these tentative detections.

Note-The units are in Jy km s -1 .

D. COMPARISONS WITH EARLIER STUDIES

Here we demonstrate our results in the context of earlier studies including molecular gas measurements for galaxies at lower redshifts. We adopt the main-sequence relation of Speagle et al. (2014), which was also used in the analysis of PHIBSS2 (Tacconi et al. 2018), to be consistent with the literature. Note that in the main text of this paper, the main-sequence relation of Boogaard et al. (2018) is utilized throughout, because this is calibrated based on the MUSE-detected sources whose average stellar masses and SFRs are lower than the sample used in Speagle et al. (2014). Therefore, the data points and the boundaries of the main-sequence relation shown in this Appendix are slightly shifted compared to the ones shown in the main text. As presented in the rightmost panels of Figures 11 and12, the results are consistent even when the MS relation of Speagle et al. (2014) is adopted.

In this study, we put constraints on the molecular gas content through CO emission at log (M * /M ) ∼ 9.5, which has only been explored at lower redshift, as illustrated in the two left panels of Figure 11. The majority of the PHIBSS2 sample at z ∼ 0 are from xCOLD GASS (Saintonge et al. 2011a(Saintonge et al. ,b, 2016(Saintonge et al. , 2017)). We find that at z ∼ 1.5, the molecular gas mass of the MS galaxies is increasing with increasing stellar mass, at least from log (M * /M ) ∼ 9.5 (top right). At lower redshift, this trend is seen starting from log (M * /M ) ∼ 9.0, but at smaller gas masses (top left panel; see also Saintonge et al. 2017). In the bottom panels, a constant gas-to-stellar mass ratio between log (M * /M ) ∼ 9.5 and ∼ 10.5 is shown at both z ∼ 0 and z ∼ 1.5 for galaxies lying on the MS. Then the gas-to-stellar mass ratio declines towards higher stellar masses.

In Figure 12, we show the same plots as in Figures 6 and7 in the main text, but now the main-sequence relation of Speagle et al. ( 2014) is adopted. We also include the best-fit lines from Tacconi et al. (2018).