Open Access
Issue
A&A
Volume 674, June 2023
Article Number A166
Number of page(s) 7
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202345977
Published online 16 June 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Understanding why galaxies die, in other words, why they become quiescent, is one of the most sought milestones in the current study of galaxy evolution. The existence of a population of dead galaxies, in contrast to the main sequence (MS) of star-forming galaxies (SFGs), has been confirmed and identified both photometrically (e.g. Daddi et al. 2005, 2007; Toft et al. 2005; Kriek et al. 2006) and spectroscopically (Toft et al. 2012; Whitaker et al. 2013; D’Eugenio et al. 2020) up to z ∼ 4 (Valentino et al. 2020) and even beyond (Carnall et al. 2023). Quenched galaxies are characterized by very low levels or even total absence of star formation activity and red colours due to the old age of their evolved stellar populations. Their total comoving mass density has remarkably increased since z ∼ 2.5 (Tomczak et al. 2014), an epoch where the most massive galaxies (log(M*/M) > 10.7) have seen their star formation rate (SFR) drastically reduced (Muzzin et al. 2013; Davidzon et al. 2017), signifying the onset of the decline of the star formation rate density of the Universe (Madau & Dickinson 2014).

The physical mechanisms responsible for the truncation of star formation in galaxies still remain uncertain, though many competing scenarios have been proposed, including: morphological quenching (Martig et al. 2009; Cornuault et al. 2018), active galactic nuclei (AGN; Di Matteo et al. 2005) and stellar feedback (Ciotti et al. 1991), virial shock heating due to massive halos (Rees & Ostriker 1977), or cosmological starvation (Feldmann & Mayer 2015) (see Man & Belli 2018 for further information).

The critical role of the interstellar medium (ISM) and more importantly of the cold molecular gas in the regulation of star formation and therefore the quenching of galaxies has also been established, playing a crucial part in all the proposed mechanisms. Measuring the amount of gas left in the ISM of quiescent galaxies (QGs) is pivotal to progress in this field. Furthermore, it is essential to observe these systems at high-z where they offer an opportunity to study QGs closer to the quenching episode, making them prime laboratories to test the various quenching processes.

The molecular gas mass (Mgas) and subsequently the gas fraction (fgas = Mgas/M*) of local QGs has been measured and found to range between 0.3−1%, typically ∼100 times lower than that of local SFGs (Young et al. 2011; Cappellari et al. 2013; Boselli et al. 2014; Davis et al. 2014; Lianou et al. 2016). Quite naturally, the detection of the observable signatures of such low levels of gas reservoirs is becoming progressively more difficult as we move to higher redshifts, posing an observational challenge for the study of the ISM of distant QGs. Nevertheless, there have been several attempts to measure the gas reservoir in high-z QGs and trace the evolution of their gas fraction with redshift, utilising CO (Sargent et al. 2015; Suess et al. 2017; Hayashi et al. 2018; Spilker et al. 2018; Bezanson et al. 2019; Belli et al. 2021; Williams et al. 2021), [CII] (Schreiber et al. 2018) and [CI] (Suzuki et al. 2022) line observations of individual galaxies as well as dust continuum emission of stacked ensembles (Gobat et al. 2018; Magdis et al. 2021, hereafter, G18 and M21) or gravitationally lensed systems (Caliendo et al. 2021; Whitaker et al. 2021a).

Intriguingly, these nascent studies have reached inconclusive and somewhat contradicting results. For example, CO(2−1) observations for a handful of very massive (log(M*/M)> 11.55) QGs at 1.2 <  z <  1.5 point towards negligible fgas <  1% remaining gas reservoirs in their ISM (Williams et al. 2021). A similar result was found with dust continuum observations of 1 <  z <  3 lensed QGs (Whitaker et al. 2021a) obtained with the Atacama Large Millimiter/submillimeter Array (ALMA). Despite the aid of lensing amplification, the majority of the targets remained undetected at 1.1 mm, placing stringent upper limits on their fgas comparable to that of local QGs (although see G18). On the other hand, far-infrared/millimeter (FIR/mm) stacking analysis of massive (log(M*/M)> 10.7), 0.5 <  z <  2.5 QGs (G18; M21) recovered dust continuum emission indicative of the presence of substantial amounts of gas and relatively large gas fractions, fgas ∼ 5 − 10%, that is, about fifty times larger than those of local ellipticals.

Although these contradictory results reveal a fully depleted versus partially depleted ISM in distant QGs and possibly point towards a range of quenching mechanisms, any attempt to draw robust conclusions is hampered by selection effects and systematic uncertainties in the conversion of the observables (i.e. line luminosities, dust continuum emission, etc.) to Mgas. Spectroscopic observations of the gas tracers (CO, [CII], and [CI]) of high-z QGs are still scarce and focused exclusively on the very massive end of the population, while the analysis of gravitationally lensed systems might also suffer from surface brightness limitations, small sample sizes, and other caveats that may introduce systematic uncertainties (Gobat et al. 2022). On the other hand, the coarse resolution of the FIR/mm observations employed in the stacking analysis of G18 and M21 (Herschel, JCMT/SCUBA-2, and ASTE/AzTEC) can result in a blending and contamination of the recovered signal and, thus, of the inferred Mgas (and fgas). The tension between the two approaches (stacking vs individual galaxies) could also stem from selection bias, with individual observations targeting the most exotic, massive high-z QGs and stacking representing the average, less massive QGs.

In this work, we aim to push the field a step further and overcome the main caveat of the stacking results presented in G18 and M21, namely, the poor resolution of the FIR/mm observations. In particular, we present a stacking analysis of 1 <  z <  3 QGs by exploiting high angular resolution ALMA observations at 1.1 mm and we present new measurements of the Mgas and fgas of the population. In Sects. 2 and 3, we present the data and the sample selection used in our study. In Sect. 4, we describe the stacking analysis and the recovered signal. In Sect. 5, we present the measured dust mass (Mdust) and resulting Mgas and we trace the evolution of fgas with redshift, while in Sect. 6, we discuss the implication of our results and explore a potential fgas − M* anti-correlation for high-z QGs. In Sect. 7, we summarise our conclusions. Throughout this work, we assume a standard ΛCDM cosmology with ΩM = 0.3, ΩΛ = 0.7 and H0 = 70 km s−1 Mpc−1. We also adopt the Salpeter initial mass function (IMF; Salpeter 1955) and the AB magnitude system.

2. ALMA data

In this work, we focus on QGs lying in GOODS-S (Dickinson et al. 2003), taking advantage of the rich multi-wavelength ancillary data in the field that we combine with 1.1 mm data from the GOODS-ALMA survey (Franco et al. 2018; Gómez-Guijarro et al. 2022). As a second step, in order to maximise the sensitivity of our stacking analysis, we complement our ALMA dataset with the full suite of available ALMA Band 6 archival data in GOODS-S.

2.1. GOODS-ALMA

GOODS-ALMA is a blind 1.1 mm galaxy survey in the GOODS-S field, centered at α  =  3h32m30s, δ  =   − 27° 48′00. The survey covers a contiguous area of 72.42 arcmin2 with ALMA Band 6 and it is a combination of two observing campaigns carried out at different angular resolutions at a homogenous sensitivity. The high resolution dataset has an average sensitivity of 89.0 μJy beam−1 and angular resolution of 0.251″ × 0.232″ (synthesised beam FWHM). The low-resolution dataset has an average sensitivity of 95.2 μJy beam−1 and angular resolution of 1.33″ × 0.935″. The combined mosaic reaches a factor 1.5 deeper sensitivity of 68.4 μJy beam−1 at angular resolution of 0.447″ × 0.418″. The high resolution mosaic was presented in Franco et al. (2018, GOODS-ALMA 1.0). The low-resolution mosaic and the final combination of the two mosaics was presented in Gómez-Guijarro et al. (2022, GOODS-ALMA 2.0). We refer to this study for further details on the observations and data processing.

2.2. ALMA archive

We downloaded all ALMA Band-6 projects that cover (i.e. with a primary beam correction better than 0.2) 1 of our 121 QGs (see Sect. 3) from the archive. This was done using the Python package astroquery. We then calibrated all these projects using the calibration scripts provided in the archive by the ALMA observatory (i.e. the so-called scriptForPI.py). Each QG is covered by at least two ALMA pointings (GOODS-ALMA 1.0 and GOODS-ALMA 2.0) with four ALMA pointings, on average. This implies that in the final stacking analysis, a total of 477 ALMA pointings were combined, each with a different depth and spatial resolution. Stacking in the uv-plane with appropriate weightings (see Sect. 4) makes it possible to combine all these different pointings.

3. Sample selection

In order to select QGs in the GOODS-S field, we employed the ZFOURGE catalogue (Straatman et al. 2016), which contains multi-wavelength ultraviolet (UV) to near-infrared (NIR) photometry of 30 911 K-band detected sources. The detection image is produced from a combination of Fourstar/Ks-band observations with previous K-band images in the field (Retzlaff et al. 2010; Hsieh et al. 2012; Fontana et al. 2014), reaching a 5σ detection limit that varies varies between 26.2 and 26.5 mag across the field. The catalogue was complemented with photometric redshifts (photo-z) derived with the EAZY fitting code (Brammer et al. 2008) as well as with the fundamental physical properties of the sources (M*, SFR, AV, and age) inferred by FAST (Kriek et al. 2009), after adopting the photo-z estimates from EAZY with the stellar population models of Bruzual & Charlot (2003), an exponentially declining star formation history, fixed solar metallicity, and a Calzetti et al. (2000) dust attenuation law (Straatman et al. 2016). As a sanity check, we also re-ran FAST with the same configuration and replaced photo-z estimates with spectroscopic redshifts for sources that are included in the most recent spectroscopic catalogue in GOODS-S (Garilli et al. 2021).

We first applied a selection on the ZFOURGE catalogue based on the (use = 1) flag, removing 17 612 photometrically uncertain sources or those that might suffer flux contamination from a nearby star, among other criteria (see Straatman et al. 2016, for further information). Since here we are interested in typical QGs around cosmic noon, we selected 644 sources with 1 <  z <  3 and in the stellar mass range of 10.20 < log(M*/M) < 11.50, as depicted in Fig. 1. We further narrowed down our sample for sources that fall within the GOODS-ALMA footprint, excluding galaxies lying at the edges of the map, due to their poorer sensitivity. The selected parent sample consists of 435 galaxies that meet our selection criteria.

thumbnail Fig. 1.

Galaxy selection. Density plot of the used ZFOURGE catalogue in the redshift vs stellar mass plane, consisting of a total 13 299 galaxies. The red shaded area represents the region covered by our redshift and stellar mass selection criteria (1 <  z <  3 and 10.20 < log(M*/M) < 11.50) embedding 852 sources. Orange circles correspond to the QGs that constitute our final selection.

We then proceeded to select QGs from the parent sample using the UVJ criterion (Williams et al. 2009) after inferring the rest-frame colour from the best fit FAST spectral energy distributions (SEDs). For our purposes, we adopted the slightly modified colour selection introduced by Schreiber et al. (2015):

{ U V > 1.3 , V J < 1.6 , U V > 0.88 × ( V J ) + 0.49 , $$ \begin{aligned} {\left\{ \begin{array}{ll} U-V > 1.3,\\ V-J < 1.6,\\ U-V > 0.88 \times (V-J) + 0.49, \end{array}\right.} \end{aligned} $$(1)

to take into account the different photometric coverage and the uncertainties in the zero-point corrections. The UVJ colour diagram along with the 140 QGs that meet the adopted colour criteria are shown in Fig. 2.

thumbnail Fig. 2.

UVJ colour diagram. Distribution of the parent sample of 435 galaxies that meet our selection criteria in the U − V, V − J colour–colour space, colour-coded by their log(SFR). The red box represents the quiescent region limits defined in Schreiber et al. (2015), which enclose the 140 QGs from which we drew our final sample.

In order to validate the quiescent nature of the selected galaxies, we also performed some additional tests. First, we compared the inferred SFR of the galaxies in our sample to that of the MS galaxies at the corresponding redshift and stellar mass, quantified as δMS = SFR/SFRMS, using the MS prescription of Schreiber et al. (2015). Adopting a δMS <  1/5 cut, we found that our entire sample is also quiescent according to their SFR levels. In the process, we performed a quality check by visually inspecting the best-fit SEDs. We also confirmed that our selected galaxies are undetected (down to 2σ) in the Herschel/PACS and SPIRE bands, adopting the GOODS-Herschel catalogues presented in Elbaz et al. (2011). They are also undetected in the 1.1 mm GOODS-ALMA 2.0 map (Gómez-Guijarro et al. 2022). However, we chose to include galaxies that are detected at MIPS 24 μm (but not in PACS/SPIRE/ALMA) to avoid biasing our sample against QGs with strong mid-IR emission originating from evolved old stellar populations or intermediate AGN activity (Fumagalli et al. 2014). This decision does not affect the recovered signal in our stacking analysis. Finally, to avoid contamination from neighbouring sources in our stack, we removed QGs that lie within a distance of 3″ from a GOODS-ALMA detected source. The final sample consist of 121 QGs, with a median z of 1.53 ± 0.03 and a median M* of (5.5 ± 0.2) × 1010M.

4. Stacking analysis

We performed a stacking analysis of the final sample of 121 QGs at hand in the uv-plane (see methodology in Gómez-Guijarro et al. 2022; Wang et al. 2022), which has been shown to improve the signal-to-noise ratio (S/N) of the stacked image compared to a stacking analysis performed in the image plane (Lindroos et al. 2015); most importantly, it allows us to stack datasets of different angular resolution. In particular, we applied a mean uv-stacking to the QGs sample, concatenating the visibilities from the GOODS-ALMA observations. The main advantage of the GOODS-ALMA map is the homogeneous coverage of the field that results in a uniform contribution of the stacked sources in the final image.

Subsequent aperture photometry at the phase center yielded a 3σ upper limit flux of 47.3 μJy. While the measurement is below a formal detection limit (> 3σ), the curve of growth of the emission at phase center closely follows that of the PSF from the stack, hinting that the recovered signal could originate from a real source, which might be detected if deeper observations were available. This curve of growth could also be caused by a noise peak located right at the phase center. However the probability of this scenario is extremely low (< 0.05%).

In order to improve the sensitivity of the stack, we decide to concatenate visibilities from the entire ALMA archival observations in Band 6 in the GOODS-S field, as described in Sect. 2.2 (see Fig. 3). Once more, we perform aperture photometry at the phase center. We applied a minor re-centering of the aperture position according to the uncertainty in the astrometric accuracy of the ALMA pointings, which, according to the ALMA technical handbook, is calculated as:

Δ S = Beam FWHP 0.9 × S / N · $$ \begin{aligned} {\Delta }S = \frac{{\mathrm{Beam} }_{\mathrm{FWHP} }}{0.9 \times S/N}\cdot \end{aligned} $$(2)

thumbnail Fig. 3.

Stacked image. 7″ × 7″ image cutout of the ALMA Band 6 archival uv-stack centered at the phase center. Blue and red contours represent the 2σ and 3σ flux emission levels, respectively.

Approximating the beam size to BeamFWHP ∼ 0.5″, we estimated an astrometric uncertainty of ΔS ∼ 0.18″. We then relocated the aperture center position by Δα = −0.05″ and Δδ = −0.10″, these being the changes in right ascension and declination, respectively, in order to optimize the S/N.

To estimate the uncertainty of the measured flux density, we randomly placed 3000 empty apertures around the phase center and adopted the standard deviation of their measured flux density distribution as the error. The diameter of the placed apertures is the same as the one we selected to make the source flux density measurement. We selected an aperture size that maximizes the S/N of the recovered flux and is also consistent with the expected dust component size of QGs at z ∼ 1.5. Recent results indicate that the dust continuum emission is expected to be smaller than the stellar emission (Whitaker et al. 2021a). Therefore, using a prior conservative estimate of the stellar component size for QGs at z ∼ 1.5 and log(M*/M) = 10.7, namely, Reff ∼ 0.2″ (van der Wel et al. 2014), we defined an aperture with Dap = 0.9″, allowing for the capture of the entire dust emission of the stack. The flux is then corrected by the appropriate aperture correction to account for the flux losses outside the aperture, calculated by dividing the flux within the aperture of Dap = 0.9″ by the flux enclosed in the synthesized dirty beam of the stack (BeamFWHM ≈ 0.5″) within the same aperture (normalised to its maximum value, see Fig. 4), yielding a correction factor of 3.26.

thumbnail Fig. 4.

Image analysis. Curve of growth at the phase center of the ALMA archival uv-stacked image (black), with its associated uncertainty shown as grey shaded area, and of the PSF (green). The blue dotted line depicts the selected aperture diameter of 0.9″.

The applied methodology ultimately yields a flux measurement of 11.70 ± 3.59 μJy with S/N = 3.25 and we apply this value for the remainder of this study. Concerning the physical properties of the stacked QGs sample, the heterogeneity in the integration time of each galaxy from the ALMA archive weighs differently each individual source. The analysis of the weight distribution of the stack shows three distinct objects with weights an order of magnitude higher than the mean value. These three galaxies bias the entire stack to a misleading higher M* and a lower z when compared to the remaining sample. Therefore, we remove these outliers and perform the same stacking analysis. The resulting S/N is hardly affected (only from the third significant digit on). The weighted average physical properties of the quiescent sample vary slightly from ⟨z⟩ = 1.62 ± 0.04 and ⟨M*⟩ = (6.55 ± 0.39) × 1010M, associated to the homogeneously weighted GOODS-ALMA stack, to ⟨z⟩ = 1.47 ± 0.03 and ⟨M*⟩ = (5.00 ± 0.26) × 1010M, associated to this latter ALMA archival stack.

5. Gas fraction

The observed frame 1.1 mm emission from galaxies at z ∼ 1 − 3 samples the Rayleigh–Jeans (R-J) tail of dust emission. In order to obtain a dust mass estimate from the measured stacked flux density, we utilised the QG SED model shown in M21, which is calibrated to make valid predictions for the average QG population. For a more detailed explanation of the utilised data and techniques to construct this QG SED model, we refer to M21. This empirical model is characterised by a dust temperature (Tdust) of 20 K. The observed weak evolution in the Tdust of massive QG makes this average temperature a robust estimate for their ISM. This has lately been found consistent by simulations in Cochrane et al. (2022), who has used the FIRE (Feedback In Realistic Environments) project (Hopkins et al. 2014) to study the dependence of dust mass on FIR flux and Tdust.

Fitting a unique model for a stack ensemble covering a large redshift range intrinsically assumes that the FIR flux density at the observed 1.1 mm does not evolve significantly with cosmic time. As shown in M21, redshift has a negligible effect in the flux density in ALMA Band 6 at 1 <  z <  3. This is attributed to the negative K-correction, namely, the counter-action between the flux dimming due to increasing cosmic distance and the flux boosting due to the redshifting of the galaxy light together with the negative slope of the R-J tail.

For the calculation of the Mdust of the stacked ensemble we compiled, we brought our flux density measurement to the rest frame and compared it to that of the model, resulting in an scaling factor of N:

N = f ν f ν 0 , $$ \begin{aligned} {N} = \frac{f_{\nu }}{f^{0}_{\nu }}, \end{aligned} $$(3)

where fν and f ν 0 $ f^{0}_{\nu} $ are the measured flux density and the normalised SED model flux density at observed frame 1.1 mm, respectively. Since the cold gas mass scales linearly with the Mdust, and the latter can be linearly scaled from the original SED, we can calculate the corresponding gas fraction as:

f gas = N × GDR ( Z ) × M dust 0 M , $$ \begin{aligned} f_{\mathrm{gas} } = \frac{N \times \mathrm{GDR}(Z) \times M^{0}_{\rm dust}}{M_{*}}, \end{aligned} $$(4)

where M dust 0 $ M^{0}_{\mathrm{dust}} $ is the dust mass of the normalised SED model and GDR(Z) is the gas-to-dust ratio, which corresponds to GDR(Z) = 92 when adopting a universal solar metallicity (Leroy et al. 2011; Magdis et al. 2012). Following this prescription, we obtained an estimate for Mdust = (2.82 ± 0.86) × 107 M , which we convert to an Mgas = (2.63 ± 0.79) × 109 M. We then calculated the corresponding fdust = 0.05 ± 0.02% and fgas = 5.3 ± 1.8% (see Fig. 5). For completeness, we also calculated Mgas using the monochromatic gas mass estimate method presented in Scoville et al. (2017), which yields a Mgas = (1.8 ± 0.56) × 109 M and a corresponding value of fgas = 3.5 ± 1.19%. This decrease is fully consistent with the higher dust temperature adopted in this method to convert the FIR flux density to gas mass. Even though estimating Mgas through Mdust (i.e. using a metallicity-dependent GDR) introduces an additional uncertainty in our analysis, the adopted value yields the most conservative estimate. We note that a recent work by Morishita et al. (2022) on a gravitationally lensed QG has shown that the normally assumed GDR(Z) could underestimate the actual gas mass by a factor of ×1.6, while the SIMBA cosmological simulation indicates a range across four orders of magnitude in the GDR of QGs (Whitaker et al. 2021b).

thumbnail Fig. 5.

Gas and dust fractions of QGs. Top: selection of fdust and fgas measurements as a function of redshift for QGs. Circles correspond to dust derived gas fractions: this work and previous stacks studies (G18; M21) are shown in red and blue, respectively. The white circles show two different estimates (connected by a grey dotted line) for a sample of individually observed lensed galaxies. The lower values correspond to those presented in Whitaker et al. (2021a) and the upper values show the new estimates provided in Gobat et al. (2022). The grey diamonds represent CO derived fgas estimates (Sargent et al. 2015; Bezanson et al. 2019; Williams et al. 2021). The red dashed area embeds fgas measurements of local QGs obtained for the ATLAS3D sample (Young et al. 2011; Cappellari et al. 2013; Davis et al. 2014). The blue shaded area and the purple dashed line represent the best fit to the M21 data and the Gobat et al. (2020) model respectively. For reference, we add the fgas evolution of main sequence galaxies according to Liu et al. (2019). Bottom: dust and gas fraction as a function of stellar mass for measurements at z ∼ 1.5. Symbols are the same as in the top panel. The dotted line shows the fgas prediction according to the Davé et al. (2012) galaxy evolution models, color coded as a function of Mhalo. For reference, we add the fgas − M* trend measured by Magdis et al. (2012), Liu et al. (2019). Light yellow scattered diamonds and arrows mark the fgas detections and upper limits for local QGs, with the corresponding best fit plotted as a black dashed line.

In Fig. 5, we show the derived fdust and fgas as a function of redshift along with a collection of estimates from the literature including local QGs (Young et al. 2011; Cappellari et al. 2013; Boselli et al. 2014; Davis et al. 2014), high-z stacked ensembles (G18; M21), and dust continuum observations of high-z QGs (Whitaker et al. 2021a). For consistency, we derived Mdust, Mgas, fgas using the same methodology as that described in Sect. 4 to the dust continuum measurements of Whitaker et al. (2021a). This results in a mean 30% increase in fgas compared to their presented values. Nevertheless, this only reinforces the conclusions presented below and in Sect. 6. The depleted gas reservoirs of local QG with fgas < 1% is in direct contrast with the increasing fgas values of the stack ensembles estimate from z = 0 to z ∼ 1. The steep increase is followed by a flat evolution towards higher redshift. This behaviour is partially mirrored by the evolution of fgas in MS SFGs, approximately 1.5 dex higher. On the other hand, the individual observations of QGs show, on average, a much weaker evolution in the gas budget of QGs with redshift, predominantly yielding upper limits in fgas that are in tension with the estimates of stacks.

Our measurement appears to be consistent with previous stacking analyses (G18; M21), indicating a non-negligible amount of gas reservoir in the average population of QGs at z ∼ 1.5 and an increase by a factor of ×10 in fgas with respect to their local counterparts. We recall that previous studies were based on the stacking of approximately 1500 sources in Herschel/SPIRE 250, 350, 500 μm, JCMT/SCUBA-2 850 μm, ASTE/AzTEC 1.1 mm maps with a resolution that ranges from 15″ to 36″, yielding a 0.12 mJy upper limit (3σ) at 1.1 mm. The sub-arcsecond resolution of the ALMA observations considered here, yielding a clear detection at 1.1 mm, seems to suggest that the reported tension between stacked results and the studies of individual QGs reported in the literature, does not (at least not fully) originate from the poorer resolution of the former.

6. Discussion

The analysis described in the previous section yields a gas fraction of fgas = 5.3 ± 1.8% for the average population of quenched galaxies at z ∼ 1.5. This estimate is consistent with the values derived in the stacking analysis of G18 and M21, indicating that high-z galaxies (on average) are quenched with non-negligible amounts of gas still being present.

6.1. Gas fraction as a function of stellar mass

We also wanted to explore the fgas − M* relation for high-z QGs. For the star forming galaxy population, an anticorrelation between these two parameters has been found at both the local universe (Cicone et al. 2017; Saintonge et al. 2017) and high-z (Sargent et al. 2014; Liu et al. 2019). The same has also been confirmed for local QGs (Boselli et al. 2014; Saintonge & Catinella 2022). This trend yet remains to be investigated for high-z QGs due to the narrow dynamical range in M* of the data obtained up to date. So far, the scarce fgas measurements of distant QGs have essentially focused on the very massive population (log(M*/M) > 11.2) (Williams et al. 2021) or gravitationally lensed systems (Whitaker et al. 2021a). Our stacking analysis, capturing the intermediate mass population of QGs, allows us to probe this trend.

In Fig. 5, we show for the first time the fgas as a function of M* at z ∼ 1.5. As a reference, we also plot the well established trend for MS SFGs at the same redshift as well as for local QGs. First, we observe that all QGs measurements lie well below the MS (fgas ≲ 15%). More importantly, all the upper limits, consistent with the gas depletion picture, are clustered at the massive end of the plot. Indeed, spectroscopic studies and imaging of individual high-z QGs have so far been heavily biased towards the most extreme and massive systems (log M* > 11.3). Conversely, the stacks, supportive of a gas retention scenario, reflect an average gas fraction that is representative of the more typical, moderately massive QGs at z ∼ 1.5. This suggests that the apparent tension between the studies of stack ensembles and of individual sources can be attributed to selection effects, with more massive QGs having experienced a larger depletion of their gas mass reservoir at the same redshift.

Nevertheless, the large uncertainties of our estimates together with the lack of data at the low-mass end of high-z QG population prevent us from fully uncovering the dependence of fgas on M*. Indeed, whether the fgas remains flat or increases at lower stellar masses, as is the case for SFGs or local QGs, cannot be constrained with the current data.

Finally, it is worth considering some theoretical predictions regarding the dependence of fgas on M*. For that purpose, we considered the analytic model presented by Davé et al. (2012), which is based on the assumption that galaxy growth is regulated by an equilibrium between star formation, as well as gas inflows and outflows. We note that the model can only reproduce QGs at z ∼ 1.5 that reside in the most massive dark matter halos (log(Mhalo/M) > 13.5) and, thus, the comparison to our data is only valid in the high-mass end. For log(M*/M) > 11.2, the model predicts a sharp decline in the fgas of QGs, consistent with the observed upper limits (Fig. 5), further supporting the scenario where more massive QGs that reside in more massive halos exhibit lower fgas than those of lower M*.

6.2. Caveats

It is acknowledged that the stacking analysis and the subsequent derivation of Mdust and Mgas from dust-continuum observations come with a suite of inherent uncertainties, caveats, and implicit assumptions. As we discussed and addressed in the previous sections, these include the selection of QGs based on the UVJ diagram, assumption of homogeneity in the properties of the population (inherent to any stacking technique), and adopted methodology to convert the observed flux densities to physical parameters (e.g. adopted template, dust temperature, and GDR). On top of these, the spatial extent of stacked sources should also be considered.

The adopted aperture correction applied to measured flux density assumes that the object at the phase center is either a point source or marginally resolved. We test the robustness of this assumption by estimating how much flux would not be measured should the object not be a point source. To do so we convolved a number of Sérsic profiles (Sérsic 1963) with the ALMA-archive synthesized dirty beam and studied how the flux changes as a function of the source size. In Fig. 6, we show the ratio of missed flux between a point source and galaxy profiles with varying Sérsic indexes, n, as a function of the galaxy effective radius, Reff. Measurements of the stellar sizes of QGs at z ∼ 1.5 with log(M*/M) = 10.7, as traced by the optical wavelengths, yields a Reff ∼ 0.18″ (van der Wel et al. 2014), which we used as a reference. We notice that, for a typical QGs n = 4, using a conservative Reff ∼ 0.18″, the flux would be underestimated by ∼5%, yielding a fgas ∼ 6%, thus reinforcing the gas retention scenario and subsequent conclusions in Sect. 5 and in the discussion.

thumbnail Fig. 6.

Size assumption effect. Flux correction factor ratio between the assumed point source emission and a modelled galaxy profile as a function of effective radius. Each colored line shows the ratios for profiles with different Sérsic indexes. The black vertical dashed line represents a size reference of Reff based on the stellar component size of QGs at z ∼ 1.5 with log(M*/M) = 10.7 (van der Wel et al. 2014).

7. Conclusions

We conducted a 1.1 mm stacking analysis of a carefully selected sample of 1 <  z <  3 QGs, using both the GOODS-ALMA survey and the full ensemble of ALMA archival data in the GOODS-S field. We recovered, at the phase center of the image, a signal of S1.1 mm = 11.70 ± 3.59 μJy, arising from the stacked ensemble of QGs with weighted mean values of ⟨z⟩ = 1.47 and log⟨M*⟩ = 10.70 M. Using the empirical templates of M21, we estimated an average log(Mdust/M) = 7.47 ± 0.13 and log(Mgas/M) = 9.42 ± 0.14 that correspond to fdust = 0.05 ± 0.02% and fgas = 5.3 ± 1.8%. With these estimates, along with with previous stacking results and studies of individual high-z QGs, we reached the following conclusions:

  • The average population of QGs at z = 1 − 3 appears to retain a non-negligible amount of gas, with an average of fgas,= 5.3 ± 1.8%. This is in agreement with previous stacking results, but considerably larger compared to the values – or upper limits inferred by studies of individual high-z QGs.

  • While still poorly constrained due to the scarcity of data in the low mass regime, the fgas of QGs at ⟨z⟩ = 1.5 appears to remain flat up to log(M*/M) ∼ 11 and then rapidly drop towards the largest stellar masses. This could alleviate the tension between the stacking results and the studies of individual sources that so far have predominantly drawn targets from the very massive end of the high-z population of QGs.

The sub-arcsecond resolution of the ALMA observations exploited in this study has been instrumental into overcoming the main caveat of the previous stacking attempts of high-z QGs, namely, the uncertainties introduced by the coarse resolution of Herschel and other ground-based sub-millimetre and millimetre facilities. However, the debate between gas depletion versus gas retention in high-z QGs and the physical mechanisms involved in the truncation of the star formation of galaxies is still on and requires deeper continuum observations or resorting to alternative gas mass tracers (e.g. [CII]; Zanella et al. 2018).

Acknowledgments

G.E.M. and D.B.S. acknowledge financial support from the Villum Young Investigator grant 37440 and 13160 and the Cosmic Dawn Center (DAWN), funded by the Danish National Research Foundation under grant No. 140 P.D. and M.T. acknowledge support from the NWO grant 016.VIDI.189.162 (“ODIN”).

References

  1. Belli, S., Contursi, A., Genzel, R., et al. 2021, ApJ, 909, L11 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bezanson, R., Spilker, J., Williams, C. C., et al. 2019, ApJ, 873, L19 [NASA ADS] [CrossRef] [Google Scholar]
  3. Boselli, A., Cortese, L., & Boquien, M. 2014, A&A, 564, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [Google Scholar]
  5. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  6. Caliendo, J. N., Whitaker, K. E., Akhshik, M., et al. 2021, ApJ, 910, L7 [NASA ADS] [CrossRef] [Google Scholar]
  7. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862 [NASA ADS] [CrossRef] [Google Scholar]
  9. Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023, MNRAS, 520, 3974 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cicone, C., Bothwell, M., Wagg, J., et al. 2017, A&A, 604, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Ciotti, L., D’Ercole, A., Pellegrini, S., & Renzini, A. 1991, ApJ, 376, 380 [CrossRef] [Google Scholar]
  12. Cochrane, R. K., Hayward, C. C., & Anglés-Alcázar, D. 2022, ApJ, 939, L27 [Google Scholar]
  13. Cornuault, N., Lehnert, M. D., Boulanger, F., & Guillard, P. 2018, A&A, 610, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680 [NASA ADS] [CrossRef] [Google Scholar]
  15. Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156 [NASA ADS] [CrossRef] [Google Scholar]
  16. Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98 [Google Scholar]
  17. Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Davis, T. A., Young, L. M., Crocker, A. F., et al. 2014, MNRAS, 444, 3427 [NASA ADS] [CrossRef] [Google Scholar]
  19. D’Eugenio, C., Daddi, E., Gobat, R., et al. 2020, ApJ, 892, L2 [Google Scholar]
  20. Dickinson, M., Giavalisco, M., & GOODS Team 2003, in The Mass of Galaxies at Low and High Redshift, eds. R. Bender, & A. Renzini, 324 [Google Scholar]
  21. Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 [NASA ADS] [CrossRef] [Google Scholar]
  22. Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Feldmann, R., & Mayer, L. 2015, MNRAS, 446, 1939 [Google Scholar]
  24. Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Fumagalli, M., Labbé, I., Patel, S. G., et al. 2014, ApJ, 796, 35 [Google Scholar]
  27. Garilli, B., McLure, R., Pentericci, L., et al. 2021, A&A, 647, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Gobat, R., Daddi, E., Magdis, G., et al. 2018, Nat. Astron., 2, 239 [NASA ADS] [CrossRef] [Google Scholar]
  29. Gobat, R., Magdis, G., D’Eugenio, C., & Valentino, F. 2020, A&A, 644, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Gobat, R., D’Eugenio, C., Liu, D., et al. 2022, A&A, 668, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022, A&A, 658, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Hayashi, M., Tadaki, K.-I., Kodama, T., et al. 2018, ApJ, 856, 118 [Google Scholar]
  33. Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hsieh, B.-C., Wang, W.-H., Hsieh, C.-C., et al. 2012, ApJS, 203, 23 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kriek, M., van Dokkum, P. G., Franx, M., et al. 2006, ApJ, 649, L71 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kriek, M., van Dokkum, P. G., Franx, M., Illingworth, G. D., & Magee, D. K. 2009, ApJ, 705, L71 [NASA ADS] [CrossRef] [Google Scholar]
  37. Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [Google Scholar]
  38. Lianou, S., Xilouris, E., Madden, S. C., & Barmby, P. 2016, MNRAS, 461, 2856 [NASA ADS] [CrossRef] [Google Scholar]
  39. Lindroos, L., Knudsen, K. K., Vlemmings, W., Conway, J., & Martí-Vidal, I. 2015, MNRAS, 446, 3502 [NASA ADS] [CrossRef] [Google Scholar]
  40. Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40 [Google Scholar]
  41. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  42. Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]
  43. Magdis, G. E., Gobat, R., Valentino, F., et al. 2021, A&A, 647, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Man, A., & Belli, S. 2018, Nat. Astron., 2, 695 [NASA ADS] [CrossRef] [Google Scholar]
  45. Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250 [Google Scholar]
  46. Morishita, T., Abdurro’uf, Hirashita, H., et al. 2022, ApJ, 938, 144 [Google Scholar]
  47. Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18 [NASA ADS] [CrossRef] [Google Scholar]
  48. Rees, M. J., & Ostriker, J. P. 1977, MNRAS, 179, 541 [NASA ADS] [CrossRef] [Google Scholar]
  49. Retzlaff, J., Rosati, P., Dickinson, M., et al. 2010, A&A, 511, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Saintonge, A., & Catinella, B. 2022, ARA&A, 60, 319 [NASA ADS] [CrossRef] [Google Scholar]
  51. Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22 [Google Scholar]
  52. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  53. Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19 [NASA ADS] [CrossRef] [Google Scholar]
  54. Sargent, M. T., Daddi, E., Bournaud, F., et al. 2015, ApJ, 806, L20 [NASA ADS] [CrossRef] [Google Scholar]
  55. Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018, A&A, 611, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 [Google Scholar]
  58. Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41 [Google Scholar]
  59. Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103 [NASA ADS] [CrossRef] [Google Scholar]
  60. Straatman, C. M. S., Spitler, L. R., Quadri, R. F., et al. 2016, ApJ, 830, 51 [NASA ADS] [CrossRef] [Google Scholar]
  61. Suess, K. A., Bezanson, R., Spilker, J. S., et al. 2017, ApJ, 846, L14 [NASA ADS] [CrossRef] [Google Scholar]
  62. Suzuki, T. L., Glazebrook, K., Schreiber, C., et al. 2022, ApJ, 936, 61 [Google Scholar]
  63. Toft, S., van Dokkum, P., Franx, M., et al. 2005, ApJ, 624, L9 [NASA ADS] [CrossRef] [Google Scholar]
  64. Toft, S., Gallazzi, A., Zirm, A., et al. 2012, ApJ, 754, 3 [NASA ADS] [CrossRef] [Google Scholar]
  65. Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85 [Google Scholar]
  66. Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93 [Google Scholar]
  67. van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [Google Scholar]
  68. Wang, T.-M., Magnelli, B., Schinnerer, E., et al. 2022, A&A, 660, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Whitaker, K. E., van Dokkum, P. G., Brammer, G., et al. 2013, ApJ, 770, L39 [NASA ADS] [CrossRef] [Google Scholar]
  70. Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021a, Nature, 597, 485 [NASA ADS] [CrossRef] [Google Scholar]
  71. Whitaker, K. E., Narayanan, D., Williams, C. C., et al. 2021b, ApJ, 922, L30 [NASA ADS] [CrossRef] [Google Scholar]
  72. Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879 [NASA ADS] [CrossRef] [Google Scholar]
  73. Williams, C. C., Spilker, J. S., Whitaker, K. E., et al. 2021, ApJ, 908, 54 [NASA ADS] [CrossRef] [Google Scholar]
  74. Young, L. M., Bureau, M., Davis, T. A., et al. 2011, MNRAS, 414, 940 [NASA ADS] [CrossRef] [Google Scholar]
  75. Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976 [Google Scholar]

All Figures

thumbnail Fig. 1.

Galaxy selection. Density plot of the used ZFOURGE catalogue in the redshift vs stellar mass plane, consisting of a total 13 299 galaxies. The red shaded area represents the region covered by our redshift and stellar mass selection criteria (1 <  z <  3 and 10.20 < log(M*/M) < 11.50) embedding 852 sources. Orange circles correspond to the QGs that constitute our final selection.

In the text
thumbnail Fig. 2.

UVJ colour diagram. Distribution of the parent sample of 435 galaxies that meet our selection criteria in the U − V, V − J colour–colour space, colour-coded by their log(SFR). The red box represents the quiescent region limits defined in Schreiber et al. (2015), which enclose the 140 QGs from which we drew our final sample.

In the text
thumbnail Fig. 3.

Stacked image. 7″ × 7″ image cutout of the ALMA Band 6 archival uv-stack centered at the phase center. Blue and red contours represent the 2σ and 3σ flux emission levels, respectively.

In the text
thumbnail Fig. 4.

Image analysis. Curve of growth at the phase center of the ALMA archival uv-stacked image (black), with its associated uncertainty shown as grey shaded area, and of the PSF (green). The blue dotted line depicts the selected aperture diameter of 0.9″.

In the text
thumbnail Fig. 5.

Gas and dust fractions of QGs. Top: selection of fdust and fgas measurements as a function of redshift for QGs. Circles correspond to dust derived gas fractions: this work and previous stacks studies (G18; M21) are shown in red and blue, respectively. The white circles show two different estimates (connected by a grey dotted line) for a sample of individually observed lensed galaxies. The lower values correspond to those presented in Whitaker et al. (2021a) and the upper values show the new estimates provided in Gobat et al. (2022). The grey diamonds represent CO derived fgas estimates (Sargent et al. 2015; Bezanson et al. 2019; Williams et al. 2021). The red dashed area embeds fgas measurements of local QGs obtained for the ATLAS3D sample (Young et al. 2011; Cappellari et al. 2013; Davis et al. 2014). The blue shaded area and the purple dashed line represent the best fit to the M21 data and the Gobat et al. (2020) model respectively. For reference, we add the fgas evolution of main sequence galaxies according to Liu et al. (2019). Bottom: dust and gas fraction as a function of stellar mass for measurements at z ∼ 1.5. Symbols are the same as in the top panel. The dotted line shows the fgas prediction according to the Davé et al. (2012) galaxy evolution models, color coded as a function of Mhalo. For reference, we add the fgas − M* trend measured by Magdis et al. (2012), Liu et al. (2019). Light yellow scattered diamonds and arrows mark the fgas detections and upper limits for local QGs, with the corresponding best fit plotted as a black dashed line.

In the text
thumbnail Fig. 6.

Size assumption effect. Flux correction factor ratio between the assumed point source emission and a modelled galaxy profile as a function of effective radius. Each colored line shows the ratios for profiles with different Sérsic indexes. The black vertical dashed line represents a size reference of Reff based on the stellar component size of QGs at z ∼ 1.5 with log(M*/M) = 10.7 (van der Wel et al. 2014).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.