ALMA Measures Rapidly Depleted Molecular Gas Reservoirs in Massive Quiescent Galaxies at z ∼ 1.5

The five galaxies confirmed with Keck were originally selected for spectroscopy from the NEWFIRM Medium Band Survey (NMBS; Whitaker et al. 2011). NMBS includes multi-wavelength photometry from the UV to 24 μm, and in particular, medium-band near-IR filters that sample the Balmer/4000 Å break at our target redshifts. The Subaru target from the sample of Onodera et al. (2012) was selected from the BzK color-selected catalog published in McCracken et al. (2010).

To identify our target galaxies, we used the stellar masses measured from the photometric spectral energy distribution (SED) as published in the original studies. These works fit the photometry using FAST (Kriek et al. 2009) to estimate stellar masses assuming Bruzual & Charlot (2003) stellar population models with exponentially declining star formation history (SFH), a Chabrier (2003) IMF and Calzetti et al. (2000) dust attenuation. Where relevant, we convert literature measurements to Chabrier IMF for comparison to our measurements. Not all targets had stellar ages measured in the literature (based on the UV to IR photometry), but where available they indicate old stellar ages (1–1.5 Gyr), with the exception of our pilot galaxy 21434, whose published stellar age is 800 Myr (Bezanson et al. 2019). The pilot galaxy's rest-frame colors are also the closest to the bluer post-starburst region of the UVJ-quiescent diagram (Whitaker et al. 2012; Belli et al. 2019).

For results presented herein, we re-fit the UV to near-IR photometry uniformly using the SED-fitting code prospector (Johnson et al. 2019, 2020), which uses the Flexible Stellar Populations Synthesis code (Conroy et al. 2009; Conroy & Gunn 2010). We fit using the prospector-α model framework (Leja et al. 2017), which includes a non-parametric SFH that has been shown to be more realistic and physically representative of massive galaxies (Leja et al. 2019a). For the purpose of fitting the stellar population properties probed by the UV to near-IR photometry, we augment the prospector-α model by removing emission due to active galactic nuclei (AGNs; which contribute primarily at mid-IR wavelengths) and the dust emission. Re-fitting the galaxies uniformly in this way also enables us to measure the mass-weighted age, which is more directly comparable to the cosmological simulations we present in Section 4.

We use the NMBS photometric catalog that includes medium-band near-infrared photometry for all galaxies, with the exception of one (ID 307881), which lies outside the NMBS footprint. For this galaxy we use the UltraVISTA catalog with broad-band photometry in the near-infrared (Muzzin et al. 2013b). We present the stellar population properties measured with prospector using the modified prospector-α model and default priors in Table 1. Using these fits instead of the literature values results in an average difference of ∼0.1 dex higher stellar mass. We find a similar difference in stellar masses between using a non-parametric SFH and an exponentially declining SFH within prospector. This difference in mass with assumed SFH is consistent with that characterized for massive log₁₀ ${M}_{* }$/${M}_{\odot }$ > 11 galaxies (Leja et al. 2019b). Our results do not significantly depend on the choice of SFH or its impact on measured stellar mass, which affect our measurements of ${f}_{{{\rm{H}}}_{2}}$ by less than a factor of 1.5.

Table 1. Properties of ALMA Targets

ID a	R.A.	Decl.	zspec	Mass	SFR${}_{\mathrm{UV}+\mathrm{IR}}$ b	SFR${}_{30\mathrm{Myr}}$ c	Re(kpc) d	Age e	Age f	References
22260	149.818229	2.561610	1.240	11.51 ${}_{-0.03}^{+0.04}$	3.6	5.3 ${}_{-1.91}^{+3.41}$	7.6	3.4	4.6	Bezanson+2013
20866	149.800931	2.537990	1.522	11.46 ${}_{-0.03}^{+0.03}$	12.8	0.7 ${}_{-0.68}^{+2.69}$	2.4	2.4	1.7	Bezanson+2013
34879	150.131380	2.523800	1.322	11.32 ${}_{-0.04}^{+0.04}$	22.9	1.4 ${}_{-1.20}^{+2.30}$	5.5	2.5	2.1	Belli+2015
34265	150.170160	2.481100	1.582	11.51 ${}_{-0.03}^{+0.03}$	7.4	0.3 ${}_{-0.34}^{+1.61}$	0.9	2.1	1.3	Belli+2015
21434	149.816230	2.549250	1.522	11.39 ${}_{-0.03}^{+0.03}$	19.1	0.5 ${}_{-0.49}^{+1.79}$	1.9	2.1	1.2	Bezanson+2013,2019
307881	150.648487	2.153990	1.429	11.63 ${}_{-0.03}^{+0.03}$	5.0	0.7 ${}_{-0.66}^{+1.73}$	2.7	2.7	3.2	Onodera+2012

Notes.

^a We adopt IDs as published in the source reference. ID for galaxy 34265 is from Belli et al. (2015) but is referred to as NMBS-COSMOS18265 in van de Sande et al. (2013). ^b SFR${}_{\mathrm{UV}+\mathrm{IR}}$ values correspond to those published by the UltraVISTA survey (Muzzin et al. 2013b). ^c Corresponds to the average SFR over the past 30 Myr as derived from our SED fitting with prospector. ^d Circularized half-light radius, defined as ${R}_{e}={r}_{e}\sqrt{b/a}$ where r_e is the semimajor axis and b/a is the axis ratio. Measured with GALFIT (Peng et al. 2002) ^e Mass-weighted stellar age as derived from fitting with non-parametric SFH prospector (in Gyr). ^f Mass-weighted stellar age assuming an exponentially declining SFH (in Gyr).

Download table as:  ASCIITypeset image

A second impact of the non-parametric SFH is that the mass-weighted age of the galaxies are typically older than that derived using parametric SFH (Leja et al. 2019b). While the ages measured assuming an exponentially declining model are typically of order 1–3 Gyr, the non-parametric model returns ages of order 2–3 Gyr. These imply that the major star formation episodes in our sample happened above z > 3. We list both values in Table 1. In the rest of this work, stellar age will refer to mass-weighted age, and we adopt the older ages from the non-parametric model because it is the more conservative constraint, as we will discuss in Section 4.2.

2.1.1. Estimation of the Star Formation Rates

2.1.2. Hubble Space Telescope Imaging

The ALMA observations of our target galaxies were carried out in project 2018.1.01739.S (PI: Williams) in separate observing sessions from 2018 December 18 to 2019 January 17 using the Band 3 (3 mm) receivers. We combine the results from this program with ALMA data for one similar galaxy from a previous pilot program (2015.1.00853.S; see Bezanson et al. 2019, for details).

The correlator was configured to center the CO(2–1) line within a spectral window of 1.875 GHz width, which provides ∼5500 km s⁻¹ of bandwidth centered on the expected frequency of the CO line, ≈89.3–102.9 GHz. Three additional spectral windows were used for continuum observations. Targets 20866, 22260, and 307881 were observed for a total of ∼90–100 minutes on-source, while 34265 and 34879 were each observed for about twice as long. The array was in a compact configuration yielding synthesized beam sizes ∼15–25 so as not to spatially resolve the target galaxies. Bandpass and flux calibrations were performed using J1058+0133 and gain calibration using J0948+0022. The data were reduced using the standard ALMA pipeline and the reductions checked manually. Our cleaning procedure involved first masking regions with clearly detected emission (S/N > 5) and then we used a stopping criterion of 3× the image rms.

Images of both the continuum and line emission were created using natural weighting of the visibilities to maximize sensitivity, with pixel sizes chosen to yield 5–10 pixels across the synthesized beam. The spectral cubes have a typical noise of 50–65 μJy beam⁻¹ in a 400 km s⁻¹ channel measured near the rest-frequency of the CO(2–1) line. The continuum data combined all available spectral windows, and reach a typical sensitivity of ∼5–9 μJy beam⁻¹, calculated as the rms of the non-primary-beam corrected image. All target galaxies are undetected in the continuum. However, the continuum imaging yielded several serendipitous 3 mm sources in these deep data. Two of these continuum sources were previously unknown galaxies and are presented in Williams et al. (2019).

To extract CO(2–1) spectra for each source, we used the uvmultifit package (Martí-Vidal et al. 2014) to fit pointlike sources to the visibilities, averaging together channels in order to produce a number of resulting spectra with channel widths ranging from 50 to 800 km s⁻¹. Given the low spatial resolution of the data and the compact galaxy sizes as measured in the available HST imaging, the pointlike source approximation is likely valid. For most sources, we fixed the position of the point-source component to the phase center of the ALMA data, with two exceptions detailed below, leaving only the flux density at each channel as a free parameter. The spectra of each target are shown in Figure 2.

Zoom In Zoom Out Reset image size

In source 22260, we detect a weak emission line at the correct frequency for the galaxy's redshift, but offset ∼12 ± 03 from the expected position of the target galaxy, a marginally significant offset given the signal-to-noise of these 2'' resolution data. It is not clear if this offset is spurious, due to an astrometric offset with respect to the HST imaging (although unlikely given our careful registration to GAIA), or reflective of a more complex physical scenario with a gas-rich region within this galaxy or a very nearby secondary galaxy as has been seen in high-redshift quiescent galaxies (Schreiber et al. 2018). For this source, we fit two point sources to the visibilities, fixing the position of one to the phase center (where we find no detection) and the other to the position of the slightly offset source, which is shown in Figure 2. We subsequently treat this as a real detection of CO(2–1) from our target. As can be seen in the HST imaging shown in Figure 2 this galaxy has a secondary optical/near-IR component (seen most prominently at HST/ACS 814W shown as inset), possibly indicating a recent minor merger. Deeper high-resolution ALMA data would be necessary to conclusively determine if the origin of the CO(2–1) emission is the secondary optical-IR component.

34879 is a similar scenario, although in that case the line emitter is brighter, offset in velocity from the redshift of our target (${\rm{\Delta }}v\sim 600$ km s⁻¹), and clearly identifiable with a nearby galaxy in the HST imaging (Figure 2). We again extract spectra by fitting multiple point sources to the visibility data for this field, fixing the positions of the sources to the phase center and the observed position of the line emitter, respectively. After this procedure, the spectra extracted at the phase centers of each field show no evidence of CO emission, although we note that the channel fluxes in these spectra are now slightly correlated with the spectra of the offset sources due to the small sky separations compared to the synthesized beam sizes.

For the detected source, 22260, we measure the line flux by fitting a simple Gaussian to the CO(2–1) spectrum. For galaxies that are not detected in CO(2–1), we set upper limits to the line flux. The undetected galaxies have velocity dispersions, σ, measured from the rest-frame optical stellar absorption features in the range $\sigma \sim 200\mbox{--}370$ km s⁻¹ (Bezanson et al. 2013; Belli et al. 2014) with the exception of 307881 for which it was not measured (Onodera et al. 2012). To measure upper limits to the integrated CO(2–1) line flux of undetected galaxies, we assume similar line widths for the stars as any molecular gas, and adopt typical values for the FWHM of the CO(2–1) of $2.355\times \sigma \,\sim$ 500–600 km s⁻¹. We use these line widths and the channel noise to set upper limits on the integrated line fluxes of each target. We note that large line widths are conservative, and that these upper limits scale with velocity interval ${\rm{\Delta }}v$ as $\sqrt{{\rm{\Delta }}v}$. Assuming a smaller line width would decrease our limiting integrated flux. The CO(2–1) line luminosity for our detected galaxy 22260, and the 1σ upper limits in the case of non-detections, are reported in Table 2.

Table 2. Molecular Gas Properties

ID	${S}_{\nu }$ a , b	${S}_{\nu }$ dvb	${L}_{{\rm{CO}}(2-1)}^{{\rm{{\prime} }}}$ b	${M}_{{{\rm{H}}}_{2}}$ c	${f}_{{{\rm{H}}}_{2}}$ c
	(μJy)	(mJy km s−1)	(108 K km s−1 pc2)	(109 M⊙)	(%)
22260 d	180 ± 38	90 ± 19	19 ± 4	10.5 ± 2.2	3.2 ± 0.7
20866	47.4	23.7	7.5	12.3	<4.3
34879 d	27.5	13.8	3.3	5.5	<2.6
34265 d	35.1	17.6	5.9	9.8	<3.0
21434	69.6	34.8	8.0	13.7	<5.5
307881	37.8	18.9	5.3	8.8	<2.1
Stack	⋯	10.3 e	2.89	4.7	<1.6

Notes.

^a Line flux is measured in a 500 km s⁻¹ channel. ^b 1σ upper limits. ^c 3σ upper limits, and assuming r₂₁ = 0.8 in temperature units, α_CO = 4.4. Molecular gas masses can be rescaled under different assumptions as ${M}_{{{\rm{H}}}_{2}}$×(0.8/r₂₁)(${\alpha }_{\mathrm{CO}}$/4.4). ^d O ii λ3727 detected in emission with Keck. ^e Assumes 400 km s⁻¹ bin. Can be scaled to width 500 km s⁻¹ by multiplying by $\sqrt{500/400}$.

Download table as:  ASCIITypeset image

To convert our measurements of CO(2–1) line luminosities into molecular gas mass (${M}_{{{\rm{H}}}_{2}}$), we make the following standard assumptions about the molecular gas conditions. We first assume a CO excitation (namely the luminosity ratio between the CO(2–1) and CO(1–0) transitions, r₂₁) following observations from the local universe. Although local galaxies with low sSFR are observed to exhibit a range of values ${r}_{21}=0.7-1$ (e.g., Saintonge et al. 2017), bulges and the central nuclei of galaxies that are thought to be similar to high-redshift compact quiescent galaxies exhibit near-thermalized excitation, ${r}_{21}=1$ (e.g., Leroy et al. 2009). For our analysis we assume ${r}_{21}=0.8$ following Spilker et al. (2018), which results in more conservative (higher) molecular gas mass measurements and limits than the assumption of thermalized emission. For comparisons to other measurements in the literature of z > 1 passive galaxies we therefore rescale other values to this excitation as ${M}_{{{\rm{H}}}_{2}}$ × (0.8/r₂₁) (including the object from Bezanson et al. 2019, which we convert to our value of r₂₁ = 0.8). Assuming a larger value for r₂₁ (e.g., 1, as is done for other studies of passive systems across redshifts) does not significantly change our results, and instead would imply even lower molecular gas fractions that further strengthen our conclusions. This assumption has a minimal impact on our ${M}_{{{\rm{H}}}_{2}}$ uncertainty budget (10%–20%) compared to e.g., a factor of ≳2 uncertainty due to the assumed value of the CO–H₂ conversion factor to translate the measured CO luminosity to ${M}_{{{\rm{H}}}_{2}}$.

In this work we assume a Milky Way–like value of ${\alpha }_{\mathrm{CO}}$ = 4.4 ${M}_{\odot }$ (K km s⁻¹ pc²)⁻¹, which is a reasonable assumption for massive galaxies with presumably high metallicities (e.g., Narayanan et al. 2012b; see also the review by Bolatto et al. 2013).

With five out of six galaxies undetected in CO(2–1) (including 21434; Bezanson et al. 2019), we perform a stacking analysis of the five non-detected galaxies. We calculate the weighted average (mean) to account for the slight differences in map rms, and use the non-primary beam corrected maps, which have Gaussian noise properties. Since the nearby companion of galaxy 34879 has significantly detected CO(2–1) emission offset by 600 km s⁻¹, but with roughly width of 200 km s⁻¹, we restrict our exploration of stacked CO(2–1) emission using image cubes with velocity resolution ≲400 km s⁻¹ to prevent the flux from the companion entering the stack, given the companion's location within 15 of the target galaxies in the stack. We construct image cubes at 400 km s⁻¹ velocity resolution centered at the rest-frequency of CO(2–1) of each galaxy, and stack the velocity bin that contains the CO(2–1) line.

We do not detect any CO(2–1) from the stack of individually undetected sources, with an rms noise limit of 25.6μJy beam⁻¹ for the 400 km s⁻¹ channel width of the stack. We use the mean redshift of the non-detected galaxies ($\langle z\rangle$ = 1.476) to put a 1σ upper limit to the average CO luminosity of ${L}_{\mathrm{CO}}^{\prime}$ ${}_{(2-1)}\lt 2.9\times {10}^{8}$ K km s⁻¹ pc². We make the same assumptions listed in Section 2.3 to convert this measurement to a molecular gas mass and find ${M}_{{{\rm{H}}}_{2}}$ $\lt \,4.7\times {10}^{9}$ ${M}_{\odot }$ (3σ). Using the average stellar mass of our undetected sample of log${}_{10}{M}_{* }/{M}_{\odot }\sim 11.5$, this puts a 3σ upper limit on the molecular gas fraction of 1.6%. The stacked sample has an average specific star formation rate of 6 × 10⁻¹¹ yr⁻¹ (likely an upper limit, as explained in Section 2.1.1) and the properties derived from the stack are summarized in Table 2.

Although this is the first systematic study using CO to measure molecular gas in quiescent galaxies at $z\gt 1$, the recent observation of average far-IR properties of 977 quiescent galaxies at $z\gt 1$ found significant dust continuum emission, implying a relatively large molecular gas content (${f}_{{{\rm{H}}}_{2}}$ ∼ 16% when converted to Chabrier IMF; Gobat et al. 2018). The individual measurements of molecular gas in our quiescent sample range from ${f}_{{{\rm{H}}}_{2}}$ ≲ 2%–6%, and are inconsistent with the average ${f}_{{{\rm{H}}}_{2}}$ measurement by Gobat et al. (2018). While a primary uncertainty in our ${f}_{{{\rm{H}}}_{2}}$ measurement is the assumed value of ${\alpha }_{\mathrm{CO}}$, extreme values only observed in low mass and low metallicity systems (${\alpha }_{\mathrm{CO}}$ $\,\gtrsim \,15;$ Narayanan et al. 2012b; Bolatto et al. 2013) would be required to bring our measurements into agreement.

A direct comparison to the Gobat et al. result is difficult owing in part to our differing methodologies, each of which is subject to its own systematic uncertainties. And, as with any photometric selection of quiescent galaxies, there is always some risk of contamination from dusty star-forming galaxies. The contamination may enter the stack either through misidentification because of the age-dust degeneracy of colors (even if only a few bright objects), or due to neighboring dusty galaxies given the low spatial resolution of the far-IR data (∼10''–30''). Neighbors can contaminate either through poor source subtraction, or as hidden dusty galaxies that do not appear in optical/near-IR selection but may remain within the far-IR photometric beam (e.g., Simpson et al. 2017; Schreiber et al. 2018; Williams et al. 2019). Further, both our sample and dusty star-forming galaxies are massive and may be strongly clustered (e.g., Hickox et al. 2012, although see also Williams et al. 2011). In this section we ignore any such possible contamination, and discuss several physical explanations for this disagreement.

First, the relatively large ${f}_{{{\rm{H}}}_{2}}$ observed by Gobat et al. (2018) could be reflecting a heterogeneity of molecular gas properties among the passive galaxy population at $z\gt 1$, as is observed at $z\lt 1$. Our sample represent some of the most massive and oldest passive galaxies known at $1\lt z\lt 1.5$, while the Gobat et al. (2018) sample is dominated by objects less massive than our sample ($\langle {\mathrm{log}}_{10}{M}_{* }\rangle =10.8$). Perhaps lower mass and/or younger additions to the red sequence still have molecular gas leftover, contributing to the far-IR emission observed on average. However, we note that because the stack is average, and our measurements are $\gt 10\times$ lower ${f}_{{{\rm{H}}}_{2}}$, a heterogenous sample would imply even larger ${f}_{{{\rm{H}}}_{2}}\gt 16 \%$ in any sample of gas-rich quiescent galaxies. More surveys that span a larger range of parameter space for individual quiescent galaxies (e.g., lower mass) are required to investigate this explanation further (K. Whitaker et al. 2021, in preparation; J. Caliendo et al. 2021, in preparation).

Alternatively, the calibration to convert the far-IR emission into a measurement of ${M}_{{{\rm{H}}}_{2}}$ might not be universal. These conversions are typically based on assumed dust to gas ratios and/or dust temperatures, calibrated using primarily star-forming galaxies (e.g., Magdis et al. 2012; Scoville et al. 2016). In theory this relies on an intrinsic relationship between dust and gas content that has been shown to accurately describe star-forming galaxies (Kaasinen et al. 2019), and for the most part, also holds for quiescent galaxies in the local universe, albeit with large scatter (e.g., Lianou et al. 2016). In principle, dust traces both atomic (H i) and molecular (H₂) gas phases, and so this could still hold if the H i/H₂ ratio is high in quiescent galaxies, while the dust to H₂ ratio is very low. We note that Spilker et al. (2018) stacked the 2mm dust continuum emission to compare to ${M}_{{{\rm{H}}}_{2}}$ measured from CO, finding consistent values between ${M}_{{{\rm{H}}}_{2}}$ observed via CO and dust, lending support for the idea that dust to H₂ conversions hold for massive quiescent galaxies at high redshift. However, Gobat et al. (2018) make the simplifying assumption that all gas traced by dust is molecular, although H i/H₂ mass ratios in local quiescent galaxies can be large (Zhang et al. 2019) and diverse (Welch et al. 2010; Boselli et al. 2014; Young et al. 2014; Calette et al. 2018). It is therefore a possibility that the significant dust emission detected by Gobat et al. (2018) is not in conflict with our low ${f}_{{{\rm{H}}}_{2}}$ measurements, and instead is primarily tracing atomic H i rather than H₂.

Nevertheless, other factors may affect these conversions, warranting further exploration. For example the dust to gas ratio can also vary with metallicity, as explored in simulations, although the extent to which this disrupts scaling relations is not clear (i.e., gas/dust may plateau above solar metallicity, applicable to most massive galaxies; Privon et al. 2018; Li et al. 2019). Future samples of quiescent galaxies with observations of both CO and dust continuum emission would reveal if the dust to ${M}_{{{\rm{H}}}_{2}}$ calibrations apply across galaxy populations at high redshift, as done locally (Smith et al. 2012). The comparison between our work presented here and that presented in Gobat et al. (2018) thus highlights several avenues of future investigation to understand the intrinsic scatter in molecular gas properties of quiescent galaxies, which will help understand the diversity of pathways that passive galaxies may take to quiescence.

Accretion is now considered to be a primary driver of galaxy growth in the early universe (for a review see Tacconi et al. 2020). While observations support this picture, it remains unclear what disrupts the growth in massive galaxies that become quiescent. Explanations include the destruction or expulsion of gas due to feedback, the suppression of gas accretion (e.g., by virial shocks once log₁₀ ${M}_{\mathrm{halo}}$/${M}_{\odot }$ > 12), or the suppression of gas collapse due to the development of a stellar bulge. Our observations of molecular gas in quenched galaxies can help discriminate between the different processes. In particular, we explore here the timescales for gas expulsion or consumption that are consistent with the low gas fractions we measure. Unfortunately with mostly upper limits to ${M}_{{{\rm{H}}}_{2}}$, and likely only upper limits to the SFR, our data set precludes a robust measurement of (current) depletion times (${t}_{\mathrm{dep}}$ = ${M}_{{{\rm{H}}}_{2}}$/ SFR) and we instead explore the allowable range of ${t}_{\mathrm{dep}}$ given low ${f}_{{{\rm{H}}}_{2}}$ and old mass-weighted stellar age.

4.2.1. Closed-box Toy Model: Constant t_dep

4.2.2. Closed-box Toy Model: Varying t_dep

Further insight is possible by comparing to analytical "bathtub" models, where the gas content of galaxies is an equilibrium of gas inflow, outflow, and consumption by star formation (e.g., Finlator & Davé 2008; Bouché et al. 2010; Davé et al. 2012; Lilly et al. 2013; Peng & Maiolino 2014; Rathaus & Sternberg 2016; Tacchella et al. 2020). This self-regulation, to first order, appears to describe the behavior of ${f}_{{{\rm{H}}}_{2}}$ across star-forming galaxy populations remarkably well (Tacconi et al. 2020). However, as halos grow above ${M}_{\mathrm{halo}}$ $\gt \,{10}^{12}$ ${M}_{\odot }$, the accretion of baryons is slowed down due to shock heating at the virial radius (e.g., Dekel & Birnboim 2006). More massive halos reach this critical mass at higher redshifts, spending a longer fraction of cosmic time without accreting new fuel for star formation.

In this section, we compare our observed gas fractions to that predicted using the simple analytic equilibrium model for ${f}_{{{\rm{H}}}_{2}}$(z, ${M}_{\mathrm{halo}}$) outlined in Davé et al. (2012). For a given halo mass and formation redshift, gas in the galaxy is computed from cosmological accretion as a function of ${M}_{\mathrm{halo}}$ (Dekel et al. 2009), simple stellar and preventative feedback prescriptions that remove gas or keep it hot in the halo, and consumption from the star formation rate. Although the gas fraction in this model is not a self-consistent model for gas evolution because it is computed from the star formation rate with an assumed star formation efficiency, this comparison nonetheless is a simple intuitive tool to qualitatively compare the relative impact of competing processes in galaxies that affect the gas fraction evolution.

In Figure 6 we show a series of these models in comparison to our observations. We show ${f}_{{{\rm{H}}}_{2}}$ for galaxies in halos that reach masses at z = 0 of ${M}_{\mathrm{halo}}$ = 10¹¹ (black), 10¹² (magenta), 10¹³ (orange), and 10¹⁴ (yellow) ${M}_{\odot }$ published in Davé et al. (2012). Observed galaxies are color coded by their inferred halo mass at the redshift of observation, using the stellar mass to halo mass relation of Behroozi et al. (2010) as implemented in halotools (Hearin et al. 2017, assuming no scatter, therefore the uncertainties are likely large).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This model predicts that only halos that reach 10¹⁴ ${M}_{\odot }$ by z = 0 halt accretion early enough in cosmic time to allow gas consumption to reach the low ${f}_{{{\rm{H}}}_{2}}$ we observe in our sample. A halo with 10¹⁴ ${M}_{\odot }$ at z = 0 reaches this critical halo mass of 10¹² ${M}_{\odot }$ at $z\sim 4$, and exceeds the quenching threshold (which evolves slightly with z) around $z\sim 3$ (Dekel et al. 2009; Davé et al. 2012). For ${M}_{\mathrm{halo}}$ ≲ 10¹³ ${M}_{\odot }$ there is not enough time to consume the gas already accreted, and other effects would be required (e.g., gas destruction from feedback) to match our low gas fractions. For M_halo $\gt \,{10}^{13}$ we also plot additional mathematical forms to describe the outflow term in the equilibrium model (variations to the stellar feedback prescription, which vary the star formation efficiency). The dotted line indicates a mass loading factor that lowers efficiency at low masses, and the dashed line indicates an additional dependence on metallicity, that decreases gas consumption at low metallicity. These variations mostly impact growth and gas fraction at low galaxy mass and improve agreement with observations at low masses, but for our case the differences are small and do not impact this result.

Based on these models we speculate that a plausible explanation of our observations is that our galaxies reside in massive halos (10¹⁴ ${M}_{\odot }$ by z = 0) that grew above the critical mass of 10¹² ${M}_{\odot }$ slowing gas accretion early enough in cosmic time ($z\sim 4$) to reach low gas fractions by $z\sim 1.5$. This scenario is qualitatively similar to the idea of cosmological starvation explored in Feldmann & Mayer (2015).

Estimated halo masses for the ALMA sample are consistent with this picture. The stellar mass to halo mass relation predicts typical ${\mathrm{log}}_{10}{M}_{\mathrm{halo}}/{M}_{\odot }$ for our sample of 13.5–14 at their respective redshifts (Behroozi et al. 2010), in general agreement with inferred halo masses from clustering of quiescent galaxies at $z\gt 1$ (e.g., Ji et al. 2018). Furthermore, the relative number density of our ALMA sample from integrating the observed stellar mass function (∼10⁻⁵ Mpc⁻³; Tomczak et al. 2014) is similar to that of log₁₀ ${M}_{\mathrm{halo}}$/${M}_{\odot }$ > 13.5 at our typical redshift $z\sim 1.4$ (halos that will reach 10¹⁴ ${M}_{\odot }$ by z = 0; calculated using the halo mass function calculator hmf published by Murray et al. 2013, assuming the halo mass function of Behroozi et al. 2013). Were our sample too numerous compared to the requisite mass halos, it would require some fraction of lower mass halos have their gas destroyed more rapidly than implied by the equilibrium model (e.g., via stronger AGN feedback). We note that these ballpark estimates are uncertain owing to scatter in the stellar mass to halo mass relation as well as uncertainties in linking progenitor populations through cumulative number density evolution (Torrey et al. 2017; Wellons & Torrey 2017).

Unfortunately, the simplicity of this analytical model and the significant intrinsic scatter in the stellar mass to halo mass relation precludes a rigorous test of the idea that reaching high halo mass and stopping accretion at early times is the primary driver of low gas fractions. We can only speculate here that this could be a contributing factor. With recent improvements in cosmological simulations, they may provide more realistic and self-consistent comparisons to observables like ${f}_{{{\rm{H}}}_{2}}$. We explore these comparisons in the next section.

Historically, cosmological simulations have been challenged to match massive galaxies in their abundances over cosmic time, as well as to prevent continued star formation in massive quiescent galaxies (for a review see Somerville & Davé 2015). Recent advances in feedback prescriptions have enabled progress on both of these fronts (Vogelsberger et al. 2014; Schaye et al. 2015; Davé et al. 2019), and now face a new challenge to match the ISM properties such as the cold gas reservoirs we study here (e.g., Narayanan et al. 2012a; Lagos et al. 2014, 2015a, 2015b). Analysis of recent cosmological hydrodynamical simulations indicate that modern implementations of feedback prescriptions for massive galaxies are able to qualitatively reproduce the global scaling relations for star-forming galaxies across cosmic time (e.g., Scoville et al. 2017; Tacconi et al. 2018; Liu et al. 2019) as well as the low ${f}_{{{\rm{H}}}_{2}}$ that are observed in massive and quiescent galaxies by $z\sim 0$ (e.g., Young et al. 2011; Davis et al. 2016; Saintonge et al. 2017). With our new observations of ${f}_{{{\rm{H}}}_{2}}$ presented here we can now extend these comparisons to massive quiescent galaxies at $z\sim 1.5$.

We compare to the predictions for molecular gas reservoirs in the simba simulation (Davé et al. 2019). simba quenches galaxies primarily via its implementation of jet AGN feedback, in which ∼10⁴ km s⁻¹ jets are ejected bipolarly from low-Eddington ratio black holes. The jets are explicitly decoupled from the ISM, thus presumably the quenching owes to heating and/or removal of halo gas. simba's X-ray feedback is important for removing H₂ from the central regions ($\lt 0.5{R}_{e};$ Appleby et al. 2020), which may also contribute to lowering the global molecular content, in general agreement with evidence for inside out quenching observed in molecular gas reservoirs (Spilker et al. 2019).

We select quiescent galaxies from a snapshot at $z\sim 1.5$ to match our ALMA target selection criteria: log₁₀ ${M}_{* }$/${M}_{\odot }$ > 11.3 and log₁₀ sSFR $\lt \,-10$ yr⁻¹. The comparison of the ${f}_{{{\rm{H}}}_{2}}$ in simba galaxies compared to our ALMA observations can be seen in Figure 7. Remarkably, simba predicts low ${f}_{{{\rm{H}}}_{2}}$ in quiescent galaxies that are consistent with our observational limits. Our ALMA limits on ${f}_{{{\rm{H}}}_{2}}$ lie at the upper envelope of ${f}_{{{\rm{H}}}_{2}}$ predicted for simba galaxies of similar mass and sSFR, with the majority of simba galaxies containing ${f}_{{{\rm{H}}}_{2}}$ < 3%. 90% of simba galaxies similar to our sample reside in halos with log₁₀ ${M}_{\mathrm{halo}}$/${M}_{\odot }$ > 13, and likely truncated accretion of new gas at earlier times. Better observational constraints on ${t}_{\mathrm{dep}}$, the time evolution of gas reservoirs, and the precision of stellar age diagnostics would be required to link this success directly to the destruction from feedback model, and/or the truncation of new gas accretion as explored in the previous section. simba produces a comparable population of "slow quenchers" and "fast quenchers" (Rodríguez Montero et al. 2019) at these redshifts, and in the future we will examine whether the galaxies consistent with our ALMA limits are preferentially in either category, and measure the associated gas depletion times.

Also in Figure 7 we show the scaling relations based on star-forming galaxies across redshifts and quiescent galaxies at $z\sim 0$ (Scoville et al. 2017; Tacconi et al. 2018; Liu et al. 2019). The shaded regions correspond to the scaling relations at $z=1.5$ for log₁₀ ${M}_{* }$/${M}_{\odot }$ = 10.8 (upper bound set by the average mass of the sample studied in Gobat et al. 2018) and log₁₀ ${M}_{* }$/${M}_{\odot }$ = 11.6 (lower bound set by the mass of our most massive galaxy). Both our ALMA limits, as well as the simba predictions, lie well below scaling relations for ${f}_{{{\rm{H}}}_{2}}({M}_{* },z,\mathrm{sSFR})$. This is consistent with the results of Section 4.2, indicating that the simulations also disagree with extrapolations of current scaling relations. Improvements to future scaling relations should include data from surveys such as this one, in the poorly explored parameter space of high redshift and low sSFR.