© The Authors 2024

1. Introduction

Cold molecular gas is the material that fuels the galaxy machinery that eventually leads to star formation. Knowing the amount of gas available in galaxies, how efficiently it is converted into stars, and how it is replenished is crucial to our understanding of their evolutionary pathways. The cosmic history of the gas mass density resembles that of the star formation rate density (SFRD; Decarli et al. 2019; Riechers et al. 2019; Magnelli et al. 2020; Walter et al. 2020), peaking at z ∼ 2 and steadily decreasing until now. The gas mass (M_gas) content in galaxies at a fixed stellar mass (M_⋆) increases with redshift – at least at 0 < z < 3. At a fixed redshift, the gas fraction (f_gas = M_gas/(M_⋆ + M_gas)) decreases with M_⋆ (Genzel et al. 2010; Béthermin et al. 2015; Morokuma-Matsui & Baba 2015; Dessauges-Zavadsky et al. 2020; Tacconi et al. 2020; Magnelli et al. 2020; Wang et al. 2022).

The relation between M_gas and M_⋆ at different redshifts has been quantified in a variety of studies (e.g., Scoville et al. 2016; Tacconi et al. 2018, 2020; Liu et al. 2019a; Kokorev et al. 2021; Wang et al. 2022), covering 0 < z < 6. It is typically parameterized according to cosmic time or redshift, and the distance from the galaxies to the main sequence (MS) of star-forming galaxies (SFGs). The term MS refers to the tight correlation that exists between the SFR and M_⋆ (e.g., Noeske et al. 2007; Elbaz et al. 2011; Whitaker et al. 2012; Speagle et al. 2014; Tomczak et al. 2016; Santini et al. 2017; Pearson et al. 2018; Barro et al. 2019), which is seen to be present at least at 0 < z < 6.

Cold molecular gas can be studied directly using rotational lines of molecular hydrogen, H₂. However, the transition probabilities are very small, the line emission is weak, and transitions are sufficiently excited only in radiation or shock-warmed molecular gas-like photodissociation regions and outflows (Parmar et al. 1991; Richter et al. 1995). Common alternatives to studying the gas content in distant galaxies include the use of the low-transition ¹²CO millimeter rotational lines and dust continuum measurements.

For the first approach, it is typically assumed that the CO lower rotational lines are optically thick and the CO line luminosity is proportional to the total molecular gas mass (M_H2), using an empirical conversion factor (Dickman et al. 1986; Solomon et al. 1987; Bolatto et al. 2013).

For the second approach, the M_gas can be derived on the basis of the dust content by converting the dust mass (M_dust) obtained by fitting the infrared (IR) spectral energy distribution (SED; Draine & Li 2007) to M_gas, for which it is typically assumed a metallicity-dependent gas-to-dust ratio (δ_GDR; e.g., Magdis et al. 2012; Genzel et al. 2015). We can also use the photometry measured in the Rayleigh–Jeans (RJ) tail of the SED (e.g., Scoville et al. 2016; Hughes et al. 2017). The method of Scoville et al. (2016, S16, hereafter) works similarly to the previous one, assuming a constant δ_GDR with a mass-weighted dust temperature (T_dust) of 25 K. These approaches assume that zero-point calibrations based on z = 0 measurements are also valid at higher redshifts.

These methods have been previously used in other works to study the gas content in the local and distant universe (e.g., Saintonge et al. 2017; Decarli et al. 2019; Freundlich et al. 2019; Sanders et al. 2023 derived M_gas using the ¹²CO rotational lines; Schinnerer et al. 2016; Wiklind et al. 2019; Kokorev et al. 2021 used the dust emission; and Liu et al. 2019a; Aravena et al. 2020; Birkin et al. 2021 used both methods). However, despite the increasing number of studies in the field, most of the efforts so far focus on individual (sub-)millimeter detections of massive objects (> 10^{10 − 11} M_⊙). For instance, the Schinnerer et al. (2016) sample is made up of ALMA detections at 240 GHz, with M_⋆ > 10^10.7 M_⊙ at z ∼ 3.2. In Freundlich et al. (2019), they include CO emitters with 10^{10 − 11.8} M_⊙ within 0.5 < z < 3. The Liu et al. (2019a) sample contains galaxies at 0.3 < z < 6 that show high-confidence ALMA detections, with median M_⋆ = 10^10.7 M_⊙.

Alternatively, other studies have sought to extend this analysis to fainter galaxies and improve the completeness of the data sets by stacking the emission of similar sources, without imposing a flux criterion on the (sub-)millimeter emission of the sources. In Tacconi et al. (2020), part of their sample is based on this strategy, made up of stacks of Herschel far infrared (FIR) spectra. Their data set also includes individual CO emitters though. In Magnelli et al. (2020), the authors measured the cosmic density of dust and gas by stacking H-band selected galaxies above a certain M_⋆. In Garratt et al. (2021), they studied the evolution of the H₂ mass density back to z ≈ 2.5, measuring the average observed 850 μm flux density of near-infrared selected galaxies. In Wang et al. (2022), they employed stacking to derive the mean mass and extent of the molecular gas of a mass-complete sample down to 10¹⁰ M_⊙. In the latter study, they obtained M_H2 values that are generally lower than previous estimates, based on individual detections.

In this work, we use the emission at 1.1 mm measured with observations obtained by the Atacama Large Millimeter/submillimeter Array (ALMA) to infer the content of gas present in a mass-complete sample of galaxies at 1.0 < z < 3.0 by analyzing stacked ALMA images for subsamples in different redshift ranges and M_⋆ bins. Taking advantage of the galaxy catalog provided by the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) in Great Observatories Origins Deep Survey (GOODS; Dickinson et al. 2003; Giavalisco et al. 2004), specifically, in GOODS-S (Guo et al. 2013, G13 hereafter), we probed the 10^{10 − 11} M_⊙ stellar mass regime with a complete sample whose 80% completeness level reaches down to 10^8.6 M_⊙ at z = 1 (10^9.2 M_⊙ at z = 3.0) (Barro et al. 2019). Our analysis is aimed at removing potential biases at the high-mass end when using detections of individual galaxies. We aspire to check whether faint sources in ALMA follow the same scaling relations derived from brighter sources or, on the contrary, whether they present a distinct molecular gas content than that prescribed for their stellar masses. Moreover, this sample gives us the chance to explore the gas reservoirs of less massive galaxies, ∼10^9 − 10 M_⊙, for which previous scaling relations are still not well calibrated.

The structure of the paper is as follows. In Sect. 2 we present the data and sample selection. We then describe the physical properties of the sample and compare them with other catalogs in Sect. 3. In Sect. 4 we present our stacking and flux measurement methodology applied to the ALMA data. In Sects. 5 and 6, we present and discuss our results regarding the gas reservoirs of our sample, comparing them with previous scaling relations. The conclusions are summarized in Sect. 7. Throughout the paper we assume a flat cosmology with Ω_M = 0.3, Ω_λ = 0.7 and a Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹. We use AB magnitudes (Oke & Gunn 1983). All M_⋆ and SFR estimations refer to a Chabrier (2003) initial mass function (IMF).

2. Data and sample

2.1. Data

We base this work on the images provided by the GOODS-ALMA 1.1 mm galaxy survey (Franco et al. 2018; Gómez-Guijarro et al. 2022a) in GOODS-S, carried out in ALMA Band 6. This survey extends over a continuous area of 72.42 arcmin² (primary beam response level ≥20%) with a homogeneous average sensitivity. It is the result of two different array configurations. Cycle 3 observations (program 2015.1.00543.S; PI: D. Elbaz) were based on a more extended array configuration that provided a high-angular-resolution data set (Franco et al. 2018). Cycle 5 observations (program 2017.1.00755.S; PI: D. Elbaz) were based on a more compact array configuration, which resulted in a lower angular resolution data set (Gómez-Guijarro et al. 2022a). In this work, we use the low-resolution data set, which has a sensitivity of 95.2 μJy beam⁻¹ and an angular resolution of (synthesized beam full width at half maximum, FWHM, along the major × minor axis). This choice is motivated by our interest in detections and flux measurements, as opposed to resolving the extent of the sources.

2.2. Sample selection

We used the source catalog provided by CANDELS in GOODS-S (G13), which includes the redshifts, M_⋆, SFR, and other SED-derived parameters for the galaxies. We selected sources in the redshift range 1.0 ≤ z ≤ 3.0, where G13 cataloged 18 459 galaxies (out of the full sample of 34 930 galaxies). While the lower limit is chosen given our interest in high redshift galaxies, especially at cosmic noon, where the gas mass density reaches its peak, the upper limit is chosen considering that the completeness at 3 < z < 4 may be compromised. As pointed out in Mérida et al. (2023, M23 hereafter), H-band-based catalogs can be affected by significant selection effects when the H-band photometric point (i.e., the flux within the F160W HST filter) lies bluewards of the Balmer break. This feature is shifted to 320 nm at z = 3 (400 nm at z = 4).

We focused on galaxies with M_⋆ > 10⁸ M_⊙, taking into account the G13 mass completeness limits. The authors in G13 computed this limit by looking for the most massive galaxies whose flux is equal to the faint limit of the sample, which evolves with redshift (Fontana et al. 2004; Grazian et al. 2015). The limits we report were obtained considering a 90% completeness limit in the flux of H = 26 mag, which corresponds to the limit computed for the shallow part of GOODS-S. Additionally, we only kept those sources with M_⋆ < 10¹¹ M_⊙ for our sample. The number density of the G13 sample sharply decreases towards > 10¹¹ M_⊙, with only 14 objects remaining within our redshift range. These all have M_⋆ slightly above ∼10¹¹ M_⊙, right at the limit, so they cannot trace a higher mass bin than the one considered in this work (10^{10 − 11} M_⊙, see Sect. 4 for more information on the divisions in redshift and M_⋆ in further analysis). Additionally, considering the uncertainty in M_⋆ for these galaxies, 8 out of the 14 sources could fall into the 10^{10 − 11} M_⊙ range. If we include these 14 galaxies in the highest mass bin considered in this work, the f_gas values barely change, showing only a variation of the order of a hundredth.

We, moreover, restrict the sample to SFGs as indicated by the UVJ diagram (Whitaker et al. 2011), which allows us to classify the galaxies in question as quiescent or star-forming, according to their rest-frame colors. This UVJ selection guarantees that the M_gas we derive, based on stacking galaxies, is not biased to lower values because of the contribution of quiescent galaxies. In Fig. 1, we show the UVJ diagram depicting the galaxies that enter the selection and those that are discarded. All these criteria leave us with a sample of 15 236 sources.

Fig. 1. Rest-frame U − V versus V − J color for the G13 sample within 1 ≤ z ≤ 3 and 108 − 11 M⊙. In red, we show the sources classified as quiescent, according to this diagram. The star-forming sources are depicted in blue.

Finally, we discarded any source lying outside the GOODS-ALMA map coverage. This coverage is defined as the area where the noise is uniform across the map, excluding the edges of the outermost pointing, where there is no pointing overlap. Our final sample is thus composed of 5530 star-forming objects located at 1.0 ≤ z ≤ 3.0, with stellar masses ranging 10^8 − 11 M_⊙. This sample shows typical uncertainties in the SED-derived parameters of 0.11 for the redshifts, 0.07 dex for the M_⋆, and 0.05 dex for the SFRs.

Within this sample, we looked for the ALMA counterparts of our galaxies, as well as for galaxies in the vicinity of sources showing an ALMA counterpart, using the source catalog listed in Gómez-Guijarro et al. (2022a, GG22 hereafter). We select those galaxies closer than 5″ to any object included in GG22. This 5″ radius is chosen considering the growth curve of the low-resolution ALMA map point spread function (PSF) and the trade-off between the number of objects in the sample and the possible contamination of ALMA detected galaxies. This condition affects just ∼3% of the objects in the sample. The effect of the exclusion of these individually detected sources and their neighbors in the photometry and f_gas values is discussed in Sects. 4 and 5.1. We refer to the subsample obtained when excluding these galaxies as the “undetected data set” hereafter.

Additionally, we also looked for counterparts of these sources in other ALMA-based catalogs, namely, the ALMA twenty-six arcmin² survey of GOODS-S at one millimeter (ASAGAO; Hatsukade et al. 2018) and the ALMA Hubble Ultra Deep Field (ALMA-HUDF; Dunlop et al. 2017; Hill et al. 2024). In Yamaguchi et al. (2020), the authors included a list of those ASAGAO sources that have an optical counterpart in the Fourstar Galaxy Evolution Survey (ZFOURGE; Straatman et al. 2016); in Hill et al. (2024), they include the optical counterparts in G13 of the sources from the ALMA-HUDF catalog. Only four sources from the undetected data set coincide with objects from ASAGAO and another four with sources from ALMA-HUDF.

We also investigated whether any of these sources belonging to the undetected data set show significant emission at (sub-) millimeter wavelengths. We measured the photometry of each of these objects using the aperture photometry method provided by photutils from Python, selecting an aperture radius of . Only 18 out of the undetected sources show a signal-to-noise ratio of S/N > 3. Only 3 out of these 18 galaxies have S/N > 3.5. When considering the S/N at the peak, only 1% of the galaxies from the undetected data set exhibit values of S/N_peak > 3.5 (including the latter 18 sources). Based on the analysis carried out in GG22, none of the sources are massive enough for the ALMA emission excess to be regarded as real and, therefore, they are indistinguishable from random noise fluctuations (see GG22 for more details). Given the low S/N of our sources, we need to analyze stacked data to look into the gas reservoirs of these galaxies.

3. Properties of the sample and comparisons with other catalogs

In this section, we compare the properties of the galaxies from our sample with those from other catalogs that were used by previous studies that also aimed at inferring the gas content of galaxies. In particular, we refer to (i) the sources from the “super-deblended” catalogs, performed in GOODS-N and in the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007), which were used to derive the Kokorev et al. (2021) scaling relation; (ii) the galaxies from the Automated ALMA Archive mining in the COSMOS field (A³COSMOS), used to obtain the Liu et al. (2019a) scaling relation; (iii) objects from Tacconi et al. (2020), based on a compilation of individually detected galaxies from different surveys and also stacks, and used to derive their scaling relation; (iv) the galaxies from the COSMOS2020 catalog that lie in the A³COSMOS footprint (COSMOS2020^* hereafter), which is the sample the Wang et al. (2022) scaling relation is based on and, finally, (v) the sources from GG22. All these comparison samples are cut to only include galaxies within the same M_⋆ and redshift intervals that we are considering in this work (see Sect. 2). Later, in Sect. 5.2, we refer to these same data sets in the context of scaling relations. In Fig. 2, we show different diagrams that highlight the properties of the listed samples together with our data set, based on G13. In Table 1, we summarize the information displayed in Fig. 2.

Fig. 2. Stellar mass vs. redshift, a color vs. color diagram based on i − H and b − i, star formation rate versus stellar mass, histograms showing the distance to the main sequence in log scale, ΔMS, and ΔMS versus stellar mass (shown from left to right and up and down). We cut the comparison samples to only include galaxies within 1.0 < z < 3.0 and having 108 − 11 M⊙. In blue, we represent our sample. The GG22 galaxies are identified in yellow. The Tacconi et al. (2020) and A3COSMOS samples are shown in gray and green, respectively. The COSMOS2020* sample is represented in magenta and the galaxies from the super-deblended catalogs are displayed in maroon. In the stellar mass vs. redshift plot, the blue contours showing our sample enclose roughly 20%, 50%, 60%, 70%, 80%, and 90% of the data. In the color-color diagram, the contours roughly enclose the 10%, 20%, 40%, 60%, 80%, and 90% of the data. In these two panels, histograms of the quantities there represented are also included, following the same color code. Quartiles are represented as horizontal segments in all the histograms. In the z-histograms shown in the first panel, we artificially elevate the baselines for the sake of clarity. The Mérida et al. (2023) MS up to 1010 M⊙ and the Barro et al. (2019) MS above 1010 M⊙ are shown in the third panel as a solid black line. The dashed lines in the third, fourth, and fifth panels show the area enclosed within 3σ with respect to the main sequence, based on the typical scatter reported in Speagle et al. (2014, ∼0.30 dex). The last panel, showing ΔMS versus the stellar mass, is split for the sake of clarity, distinguishing between the super-deblended catalogs, Tacconi et al. (2020), A3COSMOS, and GG22 (top), and our sample and COSMOS2020* (bottom). The typical uncertainties for the redshifts, stellar mass, and SFRs of our galaxies are small, ∼0.11, 0.07 dex, and 0.05 dex, respectively. In the case of the i − H and b − i colors, these are 0.14 mag and 0.20 mag, respectively.

Given that some of these works do not report the values of the rest-frame colors of their sources, the i − H and b − i colors allow us to build a diagram that works similarly to an UVJ, but using apparent magnitudes. In Fig. 2, for the panel showing this color-versus-color diagram, we use the photometry measured within the F435W (b), the F775W (i), and the F160W (H) bands from HST for our sample and for GG22. For A³COSMOS and COSMOS2020^*, we use the photometry measured within the Subaru Prime Focus Camera (Suprime-Cam) b band, the Hyper Suprime-Cam (HSC) i band, and the UltraVISTA H band. For the super-deblended catalogs and the Tacconi et al. (2020) data set we also used the HST photometry, together with the Canada-France-Hawaii Telescope (CFHT) and Subaru observations in the absence of HST data.

In Fig. 2, we also make a comparison of the position of the samples in the SFR versus the M_⋆ plane. In that panel, we include the M23 fit, defined for 1.5 < z < 2.0, up to 10¹⁰ M_⊙, and the Barro et al. (2019) MS fit (B19, hereafter) above 10¹⁰ M_⊙. The distance of each point to the MS is re-scaled to its corresponding redshift. In the fourth panel, we show the difference with respect to the MS for each galaxy in these samples, ΔMS (ΔMS = log SFR−log SFR_MS, with ΔMS = 0 dex being equivalent to δMS = SFR/SFR_MS = 1). We use M23 to calculate the ΔMS of galaxies with M_⋆/M_⊙ < 10¹⁰ and B19 for sources with higher stellar masses, given that M23 focuses on the low-mass end of the MS.

It is important to mention that for our G13-based sample, Tacconi et al. (2020), and COSMOS2020^*, the SFR were computed following the ladder technique (Wuyts et al. 2011), which combines SFR indicators at UV, mid-infrared (MIR), and FIR. For the galaxies from the super-deblended catalogs, the SFR was computed from the integrated IR luminosity (LIR) using the Daddi et al. (2010) relation. The SFRs for A³COSMOS were computed from the IR luminosity using the Kennicutt (1998) calibration. For the GG22 galaxies, the SFRs were calculated as the sum of the SFR_IR (using the Kennicutt 1998 calibration) and the SFR_UV (using the Daddi et al. 2004 calibration). We checked that only 219 of our galaxies (4%, with only 10 galaxies with M_⋆ > 10¹⁰ M_⊙) are detected in the MIR and/or FIR using Spitzer MIPS and Herschel PACS and SPIRE. This means that the SFRs of our sample come mostly from the UV emission.

3.1. This work: G13-based data set

The sample evenly populates the redshift range considered in this work (first panel of the first row in Fig. 2), with median and quartiles . The low-mass coverage of G13 allows us to reach down to 10^8 − 9 M_⊙ (log (M_⋆/M_⊙) = ). In terms of the optical colors (second panel of the first row), our sample shows typical values of 0.81 mag for i − H and 0.12 mag for b − i.

The position of our galaxies in the SFR vs M_⋆ plane (first panel of the second row) is compatible with the MS, with only a minor population of galaxies above or below three times the typical scatter (∼0.04% and ∼1% of the galaxies above/below 3σ, respectively, with σ being ∼0.3 dex, according to Speagle et al. 2014). The median ΔMS (second panel of the second row) of our galaxies is ΔMS = dex. If we check the position of our galaxies with respect to the MS according to M_⋆ (third row), we see that this typical ΔMS is maintained over the whole M_⋆ range, including the high-mass end, where most of our comparison samples are located, thus making it difficult to see our galaxies in the SFR vs M_⋆ plane above 10^9.5 M_⊙.

3.2. Comparison data sets

3.2.1. “Super-deblended” catalogs

The “super-deblended” catalogs (Liu et al. 2018; Jin et al. 2018), performed in GOODS-N and COSMOS and constructed using FIR and sub-millimeter images, use the prior positions of sources from deep Spitzer/IRAC and Very Large Array (VLA) 20 cm observations to obtain the photometry of blended FIR/sub-millimeter sources. They also employ the SED information from shorter wavelength photometry as a prior to subtract lower redshift objects. In the case of the COSMOS super-deblended catalog, the authors additionally selected a highly complete sample of priors in the Ks-band from the UltraVista catalogs. Apart from selecting those galaxies satisfying our redshift and M_⋆ cuts, we only kept the galaxies exhibiting S/N > 3 in at least three FIR to sub-millimeter bands from 100 μm to 1.2 mm, following Kokorev et al. (2021). The optical photometry of these galaxies is obtained by looking for possible optical counterparts in the CANDELS catalog performed in GOODS-N (B19) and in the COSMOS2020 catalog (Weaver et al. 2022).

The median redshift and quartiles of the galaxies satisfying our selection criteria are: , with a higher concentration of lower redshift galaxies compared to our sample and in line with T20. This data set is biased toward more massive galaxies than our sample (log M_⋆/M_⊙ = 10.4), namely, ∼1.8 dex more massive than our galaxies. Their optical i − H and b − i colors are redder than the ones traced by our sources ( mag and mag, respectively), ∼1 mag redder in both colors.

In terms of the position of these galaxies with respect to the MS, these sources show ΔMS values compatible with being MS galaxies, showing ΔMS = dex. This value corresponds though to a more star-forming data set, more compatible with the upper envelope of the MS, given the typical scatter. 6% of the galaxies show values > 3σ. If we examine the evolution of ΔMS with M_⋆, we see that below 10¹⁰ M_⊙, the sample gets increasingly star-forming, with most of the sample below 10⁹ M_⊙ surpassing ΔMS = 1 dex.

3.2.2. A³COSMOS

The A³COSMOS data set (Liu et al. 2019b) contains ∼700 galaxies (0.3 < z < 6) with high-confidence ALMA detections in the (sub-)millimeter continuum. It consists of a blind extraction, imposing an S/N_peak > 5.40, and on a prior-based extraction, using the known positions of sources in the COSMOS field, cutting the final sample to S/N_peak > 4.35. We extracted the photometry of these sources from the COSMOS2020 catalog.

The A³COSMOS galaxies with redshifts and M_⋆ in common with this work are mainly located at higher redshifts () compared to our galaxies and are also biased toward more massive objects (log (M_⋆/M_⊙) = ), ∼2 dex more massive than our sources in this case. They also display redder optical colors, with values of 1.89 mag for i − H and 1.24 mag for b − i, ∼1 mag redder in both colors than our sample.

According to their position in the SFR versus M_⋆ plane, these objects are also compatible with the MS, but (as with the galaxies from the super-deblended catalogs, T20, and GG22), they are located in the upper envelope, showing values that are nearly twice the typical scatter (ΔMS = 0.58 dex, with 13% of the galaxies above 3σ).

3.2.3. Tacconi et al. (2020)

The Tacconi et al. (2020) sample (T20 hereafter) is based on the existing literature and ALMA archive detections for individual galaxies and stacks. It consists of 2052 SFGs, where 858 of the measurements are based on CO detections, 724 on FIR dust measurements, and 470 on ∼1 mm dust measurements. We extracted their photometry by looking for the counterparts of the individual objects in the CANDELS catalogs performed in the different cosmological fields, using the catalogs already specified together with the Stefanon et al. (2017) catalog for the Extended Groth Strip (EGS; Davis et al. 2007), and Galametz et al. (2013) for the Ultra Deep Survey (UDS; Lawrence et al. 2007; Cirasuolo et al. 2007). It is however true that, since part of their sample is based on stacking, our results regarding the colors will only reflect the nature of the individual detections that make up the sample.

We see that the Tacconi et al. (2020) galaxies meeting our redshift and M_⋆ criteria are centered at , which is in line with the super-deblended sample. In terms of M_⋆, this data set is made up mostly of massive objects (log (M_⋆/M_⊙) = ), which are 2 dex more massive than our sample. According to the optical colors, this sample traces redder values of i − H and b − i, typically mag and mag for each of these colors. This is more than 1 mag redder in i − H and ∼0.8 mag redder in b − i. These galaxies show greater star formation than our sources, exhibiting ΔMS = 0.33 dex, which is compatible with them being in the upper envelope of the MS (13% of the galaxies are located above 3σ).

3.2.4. COSMOS2020^*

The COSMOS2020 catalog comprises 1.7 million sources across the 2 deg² of the COSMOS field, with ∼966 000 of them measured with all available broad-band data. Compared to COSMOS2015 (Laigle et al. 2016), it reaches the same photometric redshift precision at almost one magnitude deeper. It goes down to 10^8.43 M_⊙ at z = 1 with 70% completeness (10^9.03 M_⊙ at z = 3). We kept those galaxies that lie within the A³COSMOS footprint (which we call COSMOS2020^*), consisting of 207 129 objects. This sample is not biased toward ALMA-detected galaxies, CO emitters, or high-mass systems, which makes it more similar to our data sample.

The median redshift of the galaxies within our redshift and M_⋆ intervals is , comparable to the values we retrieve for our sample. The COSMOS2020^* data set shows a typical M_⋆ of log (M_⋆/M_⊙) = , which is ∼0.5 dex more massive than our sample. According to the optical colors, these sources show similar i − H colors ( mag) to our galaxies and around ∼0.50 mag redder colors of b − i ( mag). These objects are located well within the MS typical scatter, with a ΔMS very similar to the one we obtain for our data set (ΔMS = dex, with 4% of the galaxies above 3σ and 7% of the galaxies below 3σ).

3.2.5. GG22

GG22 presented an ALMA blind survey at 1.1 mm and built a bona fide sample of 88 sources, comprising mostly massive dusty star-forming galaxies. Half of them are detected with a purity of 100% with a S/N_peak > 5 and half of them with 3.5 ≤ S/N_peak ≤ 5, aided by the Spitzer/IRAC and the VLA prior positions. We retrieved the optical fluxes of the GG22 ALMA-selected galaxies from ZFOURGE.

The GG22 sources compatible with our redshifts and M_⋆ cuts are also biased towards high redshifts compared to our sample, similarly to A³COSMOS (). As in the cases of the super-deblended data set of T20 and A³COSMOS, we can see that GG22 is mainly made up of massive galaxies, ∼2 dex more massive than our objects (log (M_⋆/M_⊙) = ). Their optical colors are also redder than the ones showed by our sample, with median and quartiles being mag for i − H (∼1 mag redder) and mag for b − i (0.7 mag redder). These galaxies are also MS objects but, as well as the sources from the super-deblended data set, T20, and A³COSMOS, they show greater star formation than our sources (ΔMS = 0.46 dex), located above the typical scatter of the MS (20% of the galaxies above 3σ).

3.3. Comparison remarks

The main differences between our data set and the comparison samples, with the exception of COSMOS2020^*, are the blue optical colors of our galaxies, their low M_⋆ coverage, and their closer proximity to the MS. However, the latter results concerning the blue optical colors of our galaxies can be a consequence of mixing different redshifts and M_⋆ when producing the color-color diagram. We thus decided to cut it in 0.5 redshift bins and select galaxies with > 10¹⁰ M_⊙, which allows for a direct comparison with the other catalogs. Below 10¹⁰ M_⊙, as highlighted in different panels in Fig. 2, we lack sources to make a full comparison. The red colors of the comparison samples are thus due to > 10¹⁰ M_⊙ systems.

In Fig. 3, we show the color-color diagram included in Fig. 2, divided into different redshift bins. We only show the super-deblended, A³COSMOS and COSMOS2020^* galaxies as comparison data sets since the number of objects in each redshift bin included in these catalogs still provides the means to obtain meaningful number statistics to compare with.

Fig. 3. Color-versus-color diagram based on i − H and b − i in different redshift bins. We only include galaxies with M⋆ ≥ 1010 M⊙ within our data set, the super-deblended catalogs, A3COSMOS, and COSMOS2020*. See Fig. 2 for the description of the color codes and markers here shown.

When restricting our sample to galaxies with > 10¹⁰ M_⊙, the difference in i − H diminishes and we retrieve similar values to those obtained for the comparison samples. We obtained mag at 1.0 ≤ z < 1.5, and mag at 2.5 ≤ z ≤ 3.0 for our data set.

For the b − i color, we traced bluer values than the super-deblended catalogs and A³COSMOS, while getting similar results to COSMOS2020^*. The difference between the set of our sample and COSMOS2020^*, and the super-deblended catalogs and A³COSMOS increases with redshift, and the color gets bluer as well. We obtained mag at 1.0 ≤ z < 1.5 ( mag at 2.5 ≤ z ≤ 3.0) according to our data set, compared to mag at 1.0 ≤ z < 1.5 ( mag at 2.5 ≤ z ≤ 3.0) for the super-deblended catalogs, and mag at 1.0 ≤ z < 1.5 ( mag at 2.5 ≤ z ≤ 3.0) for A³COSMOS. The similarity with COSMOS2020^* and the discrepancy with the super-deblended catalogs and A³COSMOS in this color are expected. The COSMOS2020^* includes all the galaxies at these M_⋆, regardless of their flux at (sub-)millimeter wavelengths, hence being mass-complete at 10¹⁰ M_⊙, similarly to our sample. On the contrary, the super-deblended catalogs use prior positions from deep Spitzer/IRAC and VLA observations, while the A³COSMOS only considers sources with high-confidence ALMA detections, which translates to redder colors of b − i and higher dust obscurations. Our galaxies show median optical extinctions, A(V), ranging from 1.03−1.71 mag, smaller as we increase in redshift, whereas these numbers are 2.08−2.28 mag for A³COSMOS.

4. Stacking analysis and flux measurements

To study the gas content of our galaxies, we stacked the emission of objects that are similar to each other. We group galaxies according to (1) redshift and (2) log M_⋆. We distinguished: (1) four redshift bins: 1.0 ≤ z < 1.5, 1.5 ≤ z < 2.0, 2.0 ≤ z < 2.5, and 2.5 ≤ z ≤ 3.0 and (2) three M_⋆ bins: 8≤log M_⋆/M_⊙ < 9, 9 ≤ log M_⋆/M_⊙ < 10, and 10 ≤ log M_⋆/M_⊙ ≤ 11.

These divisions in redshift and stellar mass are chosen as a result of an estimation used to evaluate and maximize the probability of obtaining detections according to different combinations of redshift and M_⋆ intervals. The estimation is based on the depth of the observations and the previous knowledge about the gas reservoirs in galaxies as given by the scaling relations derived in other works (see Sect. 5.2). If we consider the expected gas fractions provided by these relations and use the δ_GDR approach (see Sects. 1 and 5.1), we can calculate the typical flux density that corresponds to those gas fractions and we can roughly infer the number of objects necessary to obtain a measurement with S/N > 3. For this, we can quantify the relation (with σ being the resulting noise in the stacked map and N the number of objects) considering different combinations of redshift and M_⋆ bins and obtain that, for getting a measurement (S/N > 3) at 10 ≤ log M_⋆/M_⊙ ≤ 11, just a few objects (< 10) are required. For the 9 ≤ log M_⋆/M_⊙ < 10 bin, we require a number of objects of the order of hundreds. Finally, for the 8 ≤ log M_⋆/M_⊙ < 9 bin, we would need tens of thousands of objects to reach the necessary depth according to our current knowledge of the gas reservoirs in galaxies. We checked that the adopted redshift division guarantees these numbers for the 9 ≤ log M_⋆/M_⊙ < 10 and 10 ≤ log M_⋆/M_⊙ ≤ 11 mass bins, while for the 8 ≤ log M_⋆/M_⊙ < 9 bin, we lack objects, regardless of how we divide in redshift; this already offers a warning that the probability to obtain a measurement in this M_⋆ bin is very low. This estimation does not, however, ensure that we are obtaining measurements for the two remaining bins, given that scaling relations are not calibrated for the kind of objects we are considering in this work, but it can still be used as a starting point.

After defining the bins, we stacked 50 × 50 arcsec² (1000 × 1000 pix²) cutouts within the low-resolution ALMA mosaic, centered at each source and using the coordinates of the centroids provided by G13. Before the stacking, we corrected these centroids for a known offset between the HST and ALMA data, reported in different studies (e.g., Dunlop et al. 2017; Franco et al. 2018). We applied the correction from Whitaker et al. (2019), which corresponds to δRA (deg) = (0.011 ± 0.08)/3600 and δDec (deg) = (−0.26 ± 0.10)/3600. We opted for median stacking galaxies instead of mean stacking them. This choice is motivated by our aim to provide an estimate of the gas reservoirs of the bulk of the SFG population, not biased toward bright sources. Additionally, this method allows us to get closer to the detection threshold in the case of the 10^9 − 10 M_⊙ bin; the use of mean stacking gives a lower S/N in 3 of the 4 redshift bins.

Despite our choice, we computed the fluxes and further physical parameters from both mean stacking and median stacking measurements. In Appendix A, we offer a quantitative discussion of the effects introduced by mean or median stacking the galaxies in the gas content. We also include analogs of some of the figures and tables appearing in this paper showing the results obtained using mean stacking.

We checked that the centroids computed using the stacked emission in ALMA are compatible with those provided by G13, based on the HST imaging, within . Following GG22, the photometry is calculated within an aperture of . This radius provides the optimal trade-off between total flux retrieval and total S/N for the GG22 sample. Then, we apply the corresponding aperture correction by dividing this flux density by that enclosed within the synthesized dirty beam (normalized to its maximum value) using the same aperture radius (see GG22 for more details). This aperture correction is ∼1.67 for .

With an S/N < 3, we calculated an upper limit for the flux density based on the surrounding sky emission in the stacked image by placing 10 000 apertures at random positions across a 20 × 20 arcsec² cutout centered at the stacked source. We measured the photometry within each of these apertures and produce a histogram with all the values, fitting the resulting Gaussian distribution leftwards to the peak, to skip the possible emission of the source. We computed the upper limit as 3 times the standard deviation of the fit. We also checked that this approach is compatible with the standard deviation obtained by iteratively drawing N (equal to the number of sources in the stack) empty positions along the mosaic, stacking them, and measuring the flux within a aperture. The compatibility of the two methods is the result of the noise uniformity along the map.

With an S/N > 3 within the aperture, we repeated the measurement using an aperture radius r = 1″. We checked that this larger radius allows us to optimize the flux and/or S/N gain and/or loss, recovering ∼7% more flux. The aperture correction for r = 1″ is ∼1.28. The uncertainty associated with the measurements is the result of the combination of the error of the stacked data (which is used to compute the S/N) and the uncertainty linked to the underlying distribution of sources that contribute to the stack. The former component is calculated by placing 10 000 r = 1″ apertures at random positions across the 50 × 50 arcsec² stacked cutout. We measured the photometry within each aperture and fit the histogram leftwards to the peak, as done in the calculation of the upper limits in the previous case. The standard deviation provided by this fit is taken as the uncertainty. The uncertainty associated with the underlying distribution is computed via bootstrapping, considering the standard deviation of the bootstrap samples.

In Fig. 4, we show an image of all the stacks. In Table 2, we list the flux densities we measure, together with the derived uncertainties. In both Fig. 4 and Table 2, we also include the results for the undetected data set, defined in Sect. 2. In Sect. 5.1, we discuss the effects of the inclusion or exclusion of the GG22 sources and their neighbors in the f_gas. For both, the whole sample and the undetected data set, we obtained S/N > 3 flux density measurements for 10^{10 − 11} M_⊙ (high-mass bin) at all redshifts. The flux density enclosed within this M_⋆ bin increases towards higher redshifts. For the intermediate-mass (10^9 − 10 M_⊙) and the low-mass bins (10^8 − 9 M_⊙), we provide 3σ upper limits. In the case of the intermediate-mass bin, we obtain a signal close to our S/N threshold at 1.5 ≤ z < 2.0, with an S/N = 2.6 (S/N = 2.5 for the undetected data set), and again at 2.5 ≤ z ≤ 3.0, with an S/N = 2.2 (S/N = 2.1 for the undetected data set). Given the lack of detection in the lower mass bins, we tested whether regrouping the galaxies, stacking all sources within 10^8 − 10 M_⊙, would allow the threshold to be exceeded. We did this for each redshift bin, and also including all galaxies at 1 ≤ z ≤ 2 and 2 < z ≤ 3. However, these tests do not report any detections; the galaxies in the low-mass bin dominate the emission of the stack.

Fig. 4. Cutouts of 7 × 7 arcsec2 of the ALMA low-resolution map showing the median stacked galaxies in each redshift and mass bin. For each bin, we include cutouts that correspond to the whole sample and the undetected data set, as defined in Sect. 2. The apertures used to measure the photometry are displayed in green. The size of the beam is shown in the first panel of the figure in gray. The flux densities and the corresponding uncertainties for each stacked galaxy are included in Table 2.

The use of a certain aperture radius in our measurements, in this case, r = 1.0″, involves some flux loss. Departure from a point-like source may involve an additional flux correction based on the galaxy morphology (see Blánquez-Sesé et al. 2023). We consider two size estimations: the size of the dust component as prescribed by GG22 and the size of the stellar component as measured and reported in G13, based on H-band data.

As pointed out in several studies, the dust component is usually more concentrated than the stellar one (e.g., Kaasinen et al. 2020; Tadaki et al. 2020; Gómez-Guijarro et al. 2022a; Liu et al. 2023). However, it is currently uncertain if our stacks, based on a mass-complete sample including faint objects, follow the latter statement, given that previous size estimations of the dust component rely on individual detections of bright objects at (sub-)millimeter wavelengths. Due to this, we also include a size estimation based on the stellar component.

According to GG22, the effective radius R_eff (the radius that contains half of the total light) of the dust component of a source with z = 1.9 and M_⋆ = 10^10.5 M_⊙ is 010. At HST H-band resolution our galaxies are fitted by a Sérsic profile characterized by a median Sérsic index n = 1.36 and a median effective radius . In Fig. 5, we show the flux correction factor that should be taken for our measurements (i.e., at 10^{10 − 11} M_⊙) versus the R_eff. Focusing on the size estimation provided by GG22, the flux correction associated with our measurements is negligible. According to the size of the stellar component, for a Sérsic index n = 1.0 − 1.5, this correction ranges from 1.17−1.22. If the size of the dust component resembled that of the stellar component, this ∼20% correction would translate to 0.08 dex larger M_gas than those reported in Table B.1 and Fig. 6.

Fig. 5. Flux ratio between a point source emission and a modeled galaxy profile versus the effective radius. The blue line shows the profile for a Sérsic index, n = 0.5, the orange one for n = 1, the green one for n = 1.5, and the red one for n = 4. The vertical solid line shows the typical size of the dust component according to Gómez-Guijarro et al. (2022a) for a z = 1.9, log M⋆/M⊙ ∼ 10.5 galaxy. The vertical dashed line shows the typical size according to the median Sérsic parameters extracted from G13 for the M⋆ ≥ 1010 M⊙ galaxies in our sample.

5. Gas reservoirs

5.1. Observed evolution of the gas reservoir of our sample

We calculated the gas content of our sample based on two approaches. The first is based on the computation of a δ_GDR using a mass metallicity relation (MZR) and the second on the RJ dust continuum emission (see Sect. 1).

For the first one, we produce synthetic spectra of the dust emission of our galaxies, according to their median redshift and ΔMS, using the Schreiber et al. (2018) IR SED template library¹. This library contains 300 templates, divided into two classes: 150 dust continuum templates due to the effect of big dust grains, and 150 templates that include the MIR features due to polycyclic aromatic hydrocarbon molecules (PAHs). These templates (which can be co-added) correspond to the luminosity that is emitted by a dust cloud with a mass equal to 1 M_⊙. After scaling each template to the measured flux density of the stacked galaxy at 1.1 mm, we obtained the LIR by integrating the rest-frame template flux between 8 and 1000 μm. This luminosity is then translated to M_dust by multiplying the intrinsic M_dust/LIR of the template by the LIR that corresponds to the measured flux density. Schreiber et al. (2018) models use different dust grain composition and emissivity yielding lower M_dust by a factor of 2 (on average), when compared to the more widely used Draine & Li (2007) models. Therefore, to have comparable M_dust with the literature results and the prescriptions needed to convert them into M_gas, we re-scaled the results based on the Schreiber et al. (2018) by an appropriate factor at each source redshift (Leroy et al., in prep.). Then, M_gas was obtained through the dust emission using the δ_GDR − Z relation derived by Magdis et al. (2012), assuming the MZR from Genzel et al. (2015), using the median M_⋆ and z of the corresponding bin. The M_gas value that we get using this approach corresponds to the total gas budget of the galaxies, including the molecular and atomic phases. As explained in Gómez-Guijarro et al. (2022b) and references therein, the molecular gas dominates over the atomic one within the physical scales probed by the dust continuum observations at this wavelength. It is worth noting that this statement has been tested within the angular scales probed by dust continuum observations, but the HI may dominate at larger scales (Chowdhury et al. 2020).

We note that this approach assumes that the emissivity index (β) adopted in the Schreiber et al. (2018) templates (∼1.5, the average value for local dwarf galaxies; Lyu et al. 2016) is accurate for our galaxies since we did not have FIR data to better constrain this parameter. Leroy et al. (in prep.) perform stacking using ALMA and Herschel data to obtain the SED of typical MS galaxies. They obtained this β via SED fitting, obtaining values that are compatible with the β = 1.5, assumed in the Schreiber et al. (2018) models. In Shivaei et al. (2022), they used stacks of Spitzer, Herschel, and ALMA photometry to examine the IR SED of high-z subsolar metallicity (∼0.5 Z_⊙) luminous IR galaxies (LIRGs), adopting β = 1.5 for their sample. In this paper, they also discuss other possible values of this parameter, but still they performed the analysis using β = 1.5.

For the second approach, we followed S16, using the corrected version of Eq. (16) from that paper. In this paper, they affirm that the luminosity-to-mass ratio at 850 μm is relatively constant under a wide range of conditions in normal star-forming and starburst galaxies, at low and high redshifts. Thus, we are able to use the measurements of the RJ flux density, derive the luminosity, and estimate M_gas. They note that this approach is equivalent to a constant δ_GDR for high-stellar-mass galaxies. This approach is justified if the variation of the mass-weighted T_dust on galactic scales is small. This is true for galaxies in the vicinity of the MS, as reported in Magnelli et al. (2014), that uses stacked Herschel data up to z ∼ 2. However, higher and a wider range of temperatures are observed in systems further away from the MS (see Clements et al. 2010; Cochrane et al. 2022). The fact that the mass-weighted T_dust keeps constant at these redshifts for MS galaxies is also supported by simulations. Liang et al. (2019), using the high-resolution cosmological simulations from the Feedback in Realistic Environments (FIRE) project, report that the mass-weighted T_dust does not strongly evolve with redshift over z = 2 − 6 at fixed IR luminosity. At a fixed redshift, it is however tightly correlated with LIR, hence sources with very high LIR, normally starburst objects, show higher mass-weighted T_dust than 25 K. We do not expect such high LIR values for our galaxies.

It is also important to take into account that the RJ tail methods can be safely applied if λ_rest ≫ hc/(k_BT_dust), where k_B is the Boltzmann constant and T_dust refers to the mass-weighted T_dust. For a mass-weighted T_dust of 25 K, this requires λ_rest ≫ 580 μm. At higher redshifts and for lower mass-weighted T_dust, the rest-frame wavelength probed by ALMA band-6 observations is far from the RJ regime (Cochrane et al. 2022). Results at z ∼ 3 should therefore be interpreted with caution.

S16 insist that the calibration samples that they use are intentionally restricted to objects with high stellar mass (M_⋆ > 5 × 10¹⁰ M_⊙), hence, they do not probe lower metallicity systems. As a consequence, we only used this approach for the calculation of the M_gas in the high-mass bin. On the other hand, for the low- and intermediate-mass bins, the M_gas values are only computed following the δ_GDR method. We discuss the effect of S16 and other prescriptions in the calculation of the gas content in lower mass galaxies in Sect. 6. An offset between both approaches, S16 and δ_GDR, is expected though at high stellar masses, as reported in Gómez-Guijarro et al. (2022b), where they found a median relative difference ()/ between both measurements.

The errors are estimated through 10 000 Monte Carlo simulations perturbing the photometry randomly within the total uncertainties and assuming an error of 0.20 dex for the metallicity (Magdis et al. 2012). All these values are included in Table B.1. As done previously for the fluxes, we also computed the error due to the underlying distribution for the median log M_⋆, z, and ΔMS using bootstrapping. When propagating these uncertainties (of the order of the hundredth) to derive the error in the metallicity we obtained values lower than 0.20 dex and, hence, lower uncertainties for the metallicity-derived parameters. As a consequence, we decided to keep the 0.2 dex value for the metallicity uncertainty, which is more conservative.

To fully understand and check the consistency of our measurements, we include two additional calculations of the gas fractions for the high-mass bin: one for the undetected data set and another considering only the GG22 galaxies. This is intended both to compare the measurement of the gas fractions that we obtain when considering the full mass-complete sample with the results (excluding individually detected sources) and to check what the contribution is for these bright individually-detected galaxies in the stacked measurements. All these values are also included in Table B.1.

In Fig. 6, we show the evolution with M_⋆ of the gas fractions derived for each redshift bin. As mentioned in Sect. 4, we provide 3σ upper limits for the low- and intermediate-mass bins and measurements for the high-mass bin. For the low-mass bin, we get f_gas < 0.97 − 0.98 along 1 ≤ z ≤ 3. These numbers are f_gas < 0.69 − 0.77 for the intermediate-mass bin. For the high-mass bin, focusing first on the δ_GDR results, we obtained f_gas = 0.32 − 0.48. Looking at the two additional cases, we checked that removing the GG22 galaxies and their neighbors from our sample (i.e., considering only the undetected data set) causes a slight drop in the measurements (f_gas = 0.30 − 0.45). The number of detected objects is much lower compared to the contribution of the undetected galaxies (see Table B.1), which dominate the emission from the stack. If we only consider the detected galaxies from GG22 and stack them following the procedure described along Sect. 4, we get (as expected) a much higher gas fractions (f_gas = 0.46 − 0.66). Taking the individual gas fractions of the GG22 objects, provided in Gómez-Guijarro et al. (2022b), and calculating the mean for each redshift bin, we get values of f_gas = 0.53 − 0.69. If we turn to the gas fractions as derived using the S16 relation, we see that the latter provides higher values of the gas fractions, but both results are compatible within the uncertainties.

Fig. 6. Gas fractions versus stellar masses derived for each redshift range. Squares represent the gas fractions obtained for our sample, circles show the gas fractions derived for the undetected data set. Diamonds represent the gas fractions calculated for the GG22 galaxies only. The differences between samples highlight the effect of the exclusion and inclusion of individually detected galaxies in the gas fractions. Gas fractions from the colored filled markers are computed following the gas-to-dust ratio prescription whereas colored empty markers with a cross within represent the gas fractions as calculated using the S16 method (only for the high-mass bin, upper limits are only computed following the gas-to-dust ratio method). Both values are included in Table B.1. The 3σ upper limits at lower stellar masses are shown with smaller squares and vertical arrows. Uncertainties are also included for the measurements. We also show the Liu et al. (2019a, L19), Tacconi et al. (2020, T20), Kokorev et al. (2021, K21), and Wang et al. (2022, W22) scaling relations in green, gray, pink, and black, respectively. For each line, there is a dashed and a solid part. The solid part represents the mass range for which these relations are calibrated, whereas the dashed one shows an extrapolation of these relations for lower stellar masses. For the T20 relation, we also include the uncertainty as a shaded gray region. The scaling relations are computed using the median redshift of each bin and a ΔMS = 0 dex. The distance from each of the points to these scaling relations is re-scaled to their corresponding redshifts and ΔMS using T20.

5.2. Scaling relation framework and comparison with our sample

Based on galaxy samples such as the ones described in Sect. 3, several works have provided different scaling relations that allow us to obtain the gas content of galaxies given the redshift, the ΔMS, and the M_⋆. These include Liu et al. (2019a, L19), T20, Kokorev et al. (2021, K21), and Wang et al. (2022, W22). The L19 relation is based on the A³COSMOS project (introduced in Sect. 3), together with ∼1000 CO-observed galaxies at 0 < z < 4 (75% of them at z < 0.1). The galaxies from the A³COSMOS project probe the log M_⋆/M_⊙ ∼ 11 − 12 MS. Complementary sources (most of them belonging to Tacconi et al. 2018; Kaasinen et al. 2019) sample the log M_⋆/M_⊙ ∼ 10 − 11 MS at z > 1. For z < 0.03, the complementary sample also covers the log M_⋆/M_⊙ ∼ 9 − 10 MS, but these authors have insisted that the metallicity-dependent CO-to-H₂ conversion factor α_CO might be more uncertain and, consequently, the estimated M_gas.

Then, T20 is based on individually detected objects plus stacks of fainter galaxies, as pointed out in Sect. 3. This relation is an expansion of the results obtained in Tacconi et al. (2018). K21 uses ∼5000 SFGs at z < 4.5, drawn from the super-deblended catalogs introduced in Sect. 3. The median redshift of the sample K21 is based on is z ∼ 0.90, with a median value of M_⋆ ∼ 4.07 × 10¹⁰ M_⊙. The low-mass and low-redshift part of their sample is restricted to galaxies that lie above the MS. Nevertheless, 69% of the galaxies qualify as MS, while 26% are classified as starbursts and 5% as passive galaxies.

W22 is based on the COSMOS2015 galaxy catalog (in Sect. 3, we refer to COSMOS2020, which is an updated version of COSMOS2015). They selected star-forming MS galaxies with ALMA band 6 or 7 coverage in the A³COSMOS database and well within the ALMA primary beam, obtaining a final sample of 3037 sources. They stacked galaxies, binning in redshift and M_⋆, in the uv domain, covering the M_⋆ range 10^{10 − 12} M_⊙. They did not select galaxies according to a certain S/N threshold at (sub-) millimeter wavelengths, so the sample includes both detected and undetected ALMA sources. In Fig. 6, the values of the scaling relations there represented are derived using the median redshift of the bin and a ΔMS = 0 dex.

If we compare our results with the previous scaling relations, for the high-mass bin, we see that the measurements for our sample are more compatible with the W22 scaling relation than with any of the latter relations. These measurements are also within the uncertainties defined by T20 but below L19 and K21. In the case of the intermediate-mass bin, the upper limits lie on the level established by the W22 and T20 extrapolations, and well below L19 and K21. In the low-mass bin, the upper limits offer poor constraints and lie above most of the scaling relations.

6. Discussion

Pushing the limit to bluer, less dusty, more MS-like, and more mass-complete samples yields lower levels of gas than those prescribed by literature scaling relations based on redder and less complete data sets. The super-deblended catalogs and the A³COSMOS samples are highly dust-obscured, exhibiting red optical colors, and they are more star-forming than our galaxies. As a consequence, the L19 and K21 scaling relations yield f_gas − M_⋆ relations with a higher normalization. Going to mass-complete samples, like the one used in the W22 relation, leads to the inclusion of blue objects with low obscurations and SFRs compatible with being right on top of the MS. As a result, the W22 relation exhibits a lower normalization, better matching our results. T20 lies between the two regimes, presumably because it is based on a combination of individually detected red, dust-obscured objects, complemented by stacks of bluer and fainter galaxies.

Regarding our results for the high-mass bin, our low values of f_gas could still be contaminated to a certain extent by the presence of post-starburst galaxies in their way to quiescence. These galaxies could have passed the UVJ screening due to their blue U − V color and thereby pulling the f_gas to lower values. It is true, however, that this effect gains importance at z > 3 (D’Eugenio et al. 2020; Forrest et al. 2020; Valentino et al. 2020), out of the redshift range considered in this work. In Antwi-Danso et al. (2023), they tested the performance of the UVJ diagram selecting quiescent galaxies, including post-starbursts. According to their results based on the Prospector-α SED modeling, the UVJ selection reaches the ∼90% completeness at z ≤ 4. They defined this completeness as the number of quiescent galaxies that are selected divided by the total number of quiescent galaxies in the sample, with quiescence being defined as a specific SFR below the threshold of the green valley.

To quantify the effect of the possible contamination due to these post-starburst objects, we produced mock sources, based on sky positions where no galaxies have been cataloged, and introduce them in the stack, checking their imprint on the resulting f_gas. Considering the UVJ selection to be 90% complete at these redshifts, we see that after introducing these mock sources, we obtained between 5% and 7% less f_gas. This difference is smaller than the uncertainty we derive for this parameter.

Additionally, the fact that we are comparing our values of f_gas, obtained using a certain method, with the results provided by scaling relations whose measurements of f_gas come from different conversion prescriptions, might be another source of discrepancy.

The L19 and W22 scaling relations rely on the RJ-tail continuum method of Hughes et al. (2017). Using this prescription to compute the f_gas of our sample, adopting α_CO = 6.5 (K km s⁻¹ pc²)⁻¹, we get 7% larger values, similarly to what we get using S16. T20 use the Leroy et al. (2011)δ_GDR together with the Genzel et al. (2015) MZR. Using this prescription, we obtained 3% larger values of f_gas. K21 is based on the Magdis et al. (2012)δ_GDR prescription together with the Mannucci et al. (2010) fundamental metallicity relation (FMR) calibrated for Kewley & Dopita (2002) that they convert to the Pettini & Pagel (2004, PP04 N2) scale following Kewley & Ellison (2008). Using this method, we get 5% larger values of f_gas.

The discrepancy between some of the scaling relations and our data is therefore not a consequence of the methodology or other factors that might be artificially pulling down our values; rather, it is simply the end result of considering a mass-complete sample that includes bluer less-dusty objects in comparison with other samples progressively more and more biased to redder dustier galaxies.

Concerning our findings for the intermediate-mass bin, it is important to take into account that at low M_⋆, the link between metallicity and the α_CO or the δ_GDR is still not well constrained and can lead to a bad estimation of the gas content. Up to date, there is very little information about the gas content of galaxies with ∼10⁹ M_⊙ at high redshifts. Most efforts have been focused thus far on galaxies at z ∼ 0 (e.g., Jiang et al. 2015; Cao et al. 2017; Saintonge et al. 2017; Madden et al. 2020; Leroy et al. 2021). According to T20 and references therein, it is hard or impossible to detect low-mass galaxies with substantially subsolar metallicity and to determine their gas content quantitatively. They suggest that there might be an interstellar medium component that might be missed or overlooked with the current techniques, such as gas and/or dust at very low temperatures. Deeper observations would be required to provide a better constraint on the f_gas of these systems.

We also tested the effect of using different prescriptions to compute the f_gas in the intermediate-mass bin. Using RJ continuum methods such as S16 or Hughes et al. (2017) yields ∼10% lower values than the ones we report. Discrepancies are expected since these methods are calibrated for more massive galaxies. S16 relies on a sample of 0.2 − 4 × 10¹¹ M_⊙ galaxies, whereas the Hughes et al. (2017) sample comprises M_⋆ ranging from 6 − 11 × 10¹⁰ M_⊙. On the other hand, the Leroy et al. (2011) prescription provides values which are compatible with our results (they differ in less than 1%), whereas the use of the Magdis et al. (2012) prescription using the Mannucci et al. (2010) FMR yields similar f_gas at z < 2; however, it start to differ at higher redshifts, where this approach reports 8% less f_gas. This difference is compatible with the uncertainties but still could reflect the fact that the metallicity of low-mass galaxies at z > 2 deviates from that observed for local galaxies. This is contrary to what is seen in higher-mass systems, whose metallicity does not evolve with redshift until z > 2.5 (Mannucci et al. 2010). This highlights the need to re-calibrate these relations for less massive objects compatible with being MS galaxies. In most cases, the low-mass sample of this kind of studies are mainly made up of galaxies showing very high SFRs.

7. Summary and conclusions

Taking advantage of the CANDELS mass-complete catalog performed as part of the GOODS-S survey (Guo et al. 2013), we are able to explore the gas content of galaxies in ALMA, using Band-6 observations at 1.1 mm (Gómez-Guijarro et al. 2022a). Our sample is composed of 5530 star forming blue (⟨b − i⟩∼0.12 mag, ⟨i − H⟩∼0.81 mag) galaxies at 1.0 ≤ z ≤ 3.0, located on the MS. It has allowed us to explore the gas content of 10^{10 − 11} M_⊙ star-forming galaxies regardless of their emission at (sub-) millimeter wavelengths. Additionally, and thanks to the stellar mass coverage and completeness of the sample, we can provide an upper limit of the gas content of lower mass galaxies at ∼10^9 − 10 M_⊙. We report measurements at 10^{10 − 11} M_⊙ and 3σ upper limits for the gas fraction at 10^8 − 10 M_⊙.

At 10^{10 − 11} M_⊙, we tend to trace lower gas fractions, f_gas = 0.32 − 0.48, than those derived from other scaling relations that use samples of redder and dustier objects, on average; this reveals a bias toward individually-detected sources at (sub-) millimeter wavelengths which are more likely to be subject to higher attenuations and also undergoing more star formation than our galaxies. Relations based on more general mass-complete samples show values that are more compatible with to those we report here.

At 10^8 − 9 M_⊙, the values we retrieved lie well above the scaling relations extrapolation, whereas at 10^9 − 10 M_⊙, the upper limits, ranging from 0.69 to 0.77, are located well within the region defined by the Wang et al. (2022) and Tacconi et al. (2020) scaling relations. The position of the upper limits at these intermediate masses supports the idea that the extrapolation derived from these scaling relations is representative of the upper bound of the underlying f_gas − M_⋆ relation, as traced by the bulk of star-forming galaxies.

Acknowledgments

R.M.M. acknowledges support from Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant PGC2018-093499-B-I00, as well as from MDM-2017-0737 Unidad de Excelencia “Maria de Maeztu” – Centro de Astrobiología (CAB), CSIC-INTA, “ERDF A way of making Europe”, and INTA SHARDS^JWST project through the PRE-SHARDSJWST/2020 PhD fellowship. C.G.G. acknowledges support from CNES. P.G.P.-G. acknowledges support from grants PGC2018-093499-B-I00 and PID2022-139567NB-I00 funded by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033, FEDER, UE. P.S.B. acknowledges support from Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through the research projects with references PID2019-107427GB-C31 and PID2022-138855NB-C31. M.F. acknowledges NSF grant AST-2009577 and NASA JWST GO Program 1727. G.E.M. acknowledges the Villum Fonden research grant 13160 “Gas to stars, stars to dust: tracing star formation across cosmic time”, grant 37440, “The Hidden Cosmos”, and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under the grant No. 140. This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Centro de Astrobiología (CAB/INTA), partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC). We made use of the following ALMA data: ADS/JAO.ALMA#2015.1.00543.S and ADS/JAO.ALMA#2017.1.00755.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST 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.

References

Appendix A: Mean stacking

In table A.1 we include the fluxes and f_gas that we obtain from mean stacking the galaxies. In Fig. A.1 we include cutouts of the resulting stacked galaxies and Fig. A.2 is an analog of Fig. 6, showing the results obtained from mean stacking.

Fig. A.1. Cutouts of 7×7 arcsec2 of the ALMA low-resolution map showing the mean stacked galaxies in each redshift and mass bin (see Fig. 4) The flux densities and the corresponding uncertainties for each stacked galaxy are included in Table A.1.

Fig. A.2. Gas fractions versus stellar masses derived for each redshift range, derived from mean stacking (also for the GG22 sample). See Fig. 6 for a description of the markers and color codes shown here.

If we compare the values of f_gas in terms of the kind of stacking method chosen, we see that mean stacking provides slightly higher values of the gas content at 10^{10 − 11} M_⊙ (f_gas = 0.34 − 0.54), getting closer to the W22 and T20 relations. Overall, this variation is compatible within the uncertainties derived for f_gas.

Appendix B: Gas reservoirs of the different galaxy samples

In Table B.1 we present a summary of the physical properties obtained for the different data sets considered in this work.

All Tables

All Figures

	Fig. 1. Rest-frame U − V versus V − J color for the G13 sample within 1 ≤ z ≤ 3 and 108 − 11 M⊙. In red, we show the sources classified as quiescent, according to this diagram. The star-forming sources are depicted in blue.
In the text

Fig. 2. Stellar mass vs. redshift, a color vs. color diagram based on i − H and b − i, star formation rate versus stellar mass, histograms showing the distance to the main sequence in log scale, ΔMS, and ΔMS versus stellar mass (shown from left to right and up and down). We cut the comparison samples to only include galaxies within 1.0 < z < 3.0 and having 108 − 11 M⊙. In blue, we represent our sample. The GG22 galaxies are identified in yellow. The Tacconi et al. (2020) and A3COSMOS samples are shown in gray and green, respectively. The COSMOS2020* sample is represented in magenta and the galaxies from the super-deblended catalogs are displayed in maroon. In the stellar mass vs. redshift plot, the blue contours showing our sample enclose roughly 20%, 50%, 60%, 70%, 80%, and 90% of the data. In the color-color diagram, the contours roughly enclose the 10%, 20%, 40%, 60%, 80%, and 90% of the data. In these two panels, histograms of the quantities there represented are also included, following the same color code. Quartiles are represented as horizontal segments in all the histograms. In the z-histograms shown in the first panel, we artificially elevate the baselines for the sake of clarity. The Mérida et al. (2023) MS up to 1010 M⊙ and the Barro et al. (2019) MS above 1010 M⊙ are shown in the third panel as a solid black line. The dashed lines in the third, fourth, and fifth panels show the area enclosed within 3σ with respect to the main sequence, based on the typical scatter reported in Speagle et al. (2014, ∼0.30 dex). The last panel, showing ΔMS versus the stellar mass, is split for the sake of clarity, distinguishing between the super-deblended catalogs, Tacconi et al. (2020), A3COSMOS, and GG22 (top), and our sample and COSMOS2020* (bottom). The typical uncertainties for the redshifts, stellar mass, and SFRs of our galaxies are small, ∼0.11, 0.07 dex, and 0.05 dex, respectively. In the case of the i − H and b − i colors, these are 0.14 mag and 0.20 mag, respectively.

In the text

	Fig. 3. Color-versus-color diagram based on i − H and b − i in different redshift bins. We only include galaxies with M⋆ ≥ 1010 M⊙ within our data set, the super-deblended catalogs, A3COSMOS, and COSMOS2020*. See Fig. 2 for the description of the color codes and markers here shown.
In the text

Fig. 4. Cutouts of 7 × 7 arcsec2 of the ALMA low-resolution map showing the median stacked galaxies in each redshift and mass bin. For each bin, we include cutouts that correspond to the whole sample and the undetected data set, as defined in Sect. 2. The apertures used to measure the photometry are displayed in green. The size of the beam is shown in the first panel of the figure in gray. The flux densities and the corresponding uncertainties for each stacked galaxy are included in Table 2.

In the text

Fig. 5. Flux ratio between a point source emission and a modeled galaxy profile versus the effective radius. The blue line shows the profile for a Sérsic index, n = 0.5, the orange one for n = 1, the green one for n = 1.5, and the red one for n = 4. The vertical solid line shows the typical size of the dust component according to Gómez-Guijarro et al. (2022a) for a z = 1.9, log M⋆/M⊙ ∼ 10.5 galaxy. The vertical dashed line shows the typical size according to the median Sérsic parameters extracted from G13 for the M⋆ ≥ 1010 M⊙ galaxies in our sample.

In the text

Fig. 6. Gas fractions versus stellar masses derived for each redshift range. Squares represent the gas fractions obtained for our sample, circles show the gas fractions derived for the undetected data set. Diamonds represent the gas fractions calculated for the GG22 galaxies only. The differences between samples highlight the effect of the exclusion and inclusion of individually detected galaxies in the gas fractions. Gas fractions from the colored filled markers are computed following the gas-to-dust ratio prescription whereas colored empty markers with a cross within represent the gas fractions as calculated using the S16 method (only for the high-mass bin, upper limits are only computed following the gas-to-dust ratio method). Both values are included in Table B.1. The 3σ upper limits at lower stellar masses are shown with smaller squares and vertical arrows. Uncertainties are also included for the measurements. We also show the Liu et al. (2019a, L19), Tacconi et al. (2020, T20), Kokorev et al. (2021, K21), and Wang et al. (2022, W22) scaling relations in green, gray, pink, and black, respectively. For each line, there is a dashed and a solid part. The solid part represents the mass range for which these relations are calibrated, whereas the dashed one shows an extrapolation of these relations for lower stellar masses. For the T20 relation, we also include the uncertainty as a shaded gray region. The scaling relations are computed using the median redshift of each bin and a ΔMS = 0 dex. The distance from each of the points to these scaling relations is re-scaled to their corresponding redshifts and ΔMS using T20.

In the text

	Fig. A.1. Cutouts of 7×7 arcsec2 of the ALMA low-resolution map showing the mean stacked galaxies in each redshift and mass bin (see Fig. 4) The flux densities and the corresponding uncertainties for each stacked galaxy are included in Table A.1.
In the text

	Fig. A.2. Gas fractions versus stellar masses derived for each redshift range, derived from mean stacking (also for the GG22 sample). See Fig. 6 for a description of the markers and color codes shown here.
In the text