From Haloes to Galaxies. III. The Gas Cycle of Local Galaxy Populations

Figure 1 shows the average H₂ molecular gas mass (left panel) and average H ɪ gas mass (right panel) in the SFR–M_* plane. In both panels, the dashed line indicates the position of the star-forming main sequence (MS) defined in Saintonge et al. (2016). The dotted lines indicate ±0.4 dex scatter around the MS. The lower dotted line is the approximate divide between star-forming and quiescent galaxies. The upper dotted line indicates the approximate divide between normal star-forming galaxies and galaxies with elevated SFR (starburst galaxies). It is evident that the color-coded data distribution in panels (a) and (b) are different. This is because panel (a) is produced by galaxies with H₂ detections regardless of their H ɪ detection status, while panel (b) is produced by galaxies with H ɪ detections regardless of their H₂ detection status. These two samples share many common galaxies that have detections in both H ɪ and H₂, but are not the same due to that some galaxies have only H ɪ detections/observations while some have only H₂ detections/observations. We discuss below the reasons that cause this difference in data distribution and the impact to the results.

Comparing the data distribution in panels (a) and (b), the small extra extension to the star-bursting regime in H₂ at log M_* ∼ 10–11 is due to the xCOLD GASS survey containing some extra star-bursting galaxies that have not been included in the xGASS survey. In other words, this is not because these star-bursting galaxies have H₂ detections but no H ɪ detections, but simply because they have not been observed in H ɪ. Meanwhile, the fact that many starbursts/LIRGs have strong 21 cm continuum in the star-bursting regions results in strong H ɪ absorption and hence obvious H ɪ emission may not be easily detected by single-dish telescopes (e.g., van Driel et al. 2001). The more obvious difference between panels (a) and (b) is in the low SFR regime, where only H ɪ data is present. First, this difference here is not due to different samples as in the star-bursting regime discussed above. These galaxy samples are identical, i.e., they have been observed in both H ɪ and H₂. Second, the observation limit in H₂ is similar to that of H ɪ, both go down to a gas fraction of about 2%. This suggests that some galaxies in the low SFR regime are genuine H₂ poor (hence, not detected in CO) but H ɪ rich (hence, detected in H ɪ). Indeed, as discussed in Zhang et al. (2019, 2021), during the quenching of the massive central disk galaxies, their H₂ gas mass drops rapidly by more than 1 dex from MS to the lowest observed sSFR, but their H ɪ gas mass remains surprisingly constant, i.e., the H ɪ/H₂ ratio rapidly increases with decreasing sSFR. If we select galaxies that have detections in both H ɪ and H₂, the results in panels (a) and (b) in the overlapped region of H ɪ and H₂ distributions on the SFR–M_* plane remain the same.

Comparing panels (a) and (b), H₂ gas mass and H ɪ gas mass follow a similar trend of distribution on the SFR–M_* plane. Both of them increase with increasing stellar mass and SFR. At a given stellar mass and SFR, on average, the H ɪ gas mass is higher than the H₂ gas mass, especially for low-mass galaxies that are H ɪ dominated. This indicates the SFE of H₂ gas (SFE${}_{{{\rm{H}}}_{2}}$) is higher than that of the H ɪ gas (SFE_HI). Or equivalently, the gas depletion timescale τ_dep (=1/SFE) of H₂ gas (${\tau }_{\mathrm{dep},{{\rm{H}}}_{2}}$) is shorter than that of the H ɪ gas (τ_dep,Hɪ). It should be noted that SFE_Hɪ can be reformed as SFE_Hɪ = SFE${}_{{{\rm{H}}}_{2}}\times {M}_{{{\rm{H}}}_{2}}/{M}_{{\rm{H}}{\rm{I}}}$. Hence, SFE_Hɪ is simply the product of SFE${}_{{{\rm{H}}}_{2}}$ and the H₂ to H ɪ gas mass ratio.

The typical timescale of H₂ formation in equilibrium is about 10⁷ yr (e.g., Goldsmith et al. 2007; Liszt 2007), which is much shorter than ${\tau }_{{{\rm{dep}},{\rm{H}}}_{2}}$ of typical star-forming galaxies. This hence suggests that a significant fraction of H ɪ gas will not be able to form H₂, in particular in low-mass galaxies, probably due to the low surface density or lack of dust. H₂ formation efficiency is also expected to depend on the gas-phase metallicity and/or the gas pressure (e.g., Blitz & Rosolowsky 2006; Krumholz et al. 2009). In addition, H₂ can also be destroyed by star formation feedback (Morselli et al. 2020).

Various studies have shown that ${\tau }_{\mathrm{dep},{{\rm{H}}}_{2}}$ is tightly correlated with the SFR offset from the MS (e.g., Saintonge et al. 2011, 2012, 2017; Huang & Kauffmann 2014; Genzel et al. 2015; Scoville et al. 2017; Tacconi et al. 2018, 2020, for a review; Eales et al. 2020). In D21, we show the ${\tau }_{\mathrm{dep},{{\rm{H}}}_{2}}$–sSFR relation is part of the FFR. Similar to Figure 4 in D21, Figure 2 shows the average SFE as a function of sSFR for H₂ (upper six panels) and for H ɪ (lower six panels). We split the galaxies into different stellar mass (panels (a) and (g)), B/T (panels (b) and (h)) and effective radius (panels (c) and (i)), centrals and satellites (panels (d) and (j)), and different gas-phase metallicities (panels (e) and (k)). In panels (f) and (l), we also show the distribution of all galaxies in the H₂ and H ɪ samples, respectively, without splitting into subsamples. The vertical light blue and red shades in each panel indicate the position of the star-forming MS and the approximate position of the quenched galaxies at M_* = 10¹⁰ M_⊙, respectively.

Since SFE = sSFR/μ, where μ is the gas to stellar mass ratio, for a given position on the SFE–sSFR plane, the value of μ (or the gas fraction f_gas) is uniquely determined. At a given sSFR, a larger SFE corresponds a lower gas fraction. The two dark diagonal lines in the upper six panels mark the constant ${\mu }_{{{\rm{H}}}_{2}}$ of 1.5% and 2.5%, respectively. These are the detection limits for the xCOLD GASS galaxies at different stellar mass (1.5% for log M_* > 10 and 2.5% for 9 < log M_* < 10). The three dark diagonal lines in the lower six panels mark the constant μ_HI of 2%, 4%, and 10%, respectively. These are the detection limits for the xGASS galaxies at different stellar mass (see Section 2). The dark shades in each panel mark the regions below the detection limit for different M_*.

The sample selection effects in xCOLD GASS and the uncertainties in SFR measurements have been discussed in detail in D21. The average H₂ detection ratio is higher than 80% at log sSFR > −1.6 Gyr⁻¹, where the results should be reliable. For H ɪ, as in Figure A1 and A2, the average H ɪ detection ratio in xGASS is higher than 80% at log sSFR > −1.5 Gyr⁻¹ (where the results should be reliable). For galaxies with log sSFR < −1.5 Gyr⁻¹, the detection ratio for both H₂ and H ɪ decreases rapidly due to the observation limit of the surveys. Hence, results in this regime are for H₂- or H ɪ-detected galaxies, which may not represent the statistics of an unbiased sample. We also note that, as shown in the upper six panels in Figure 2, the general trend of the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation seems to persist below log sSFR < −1.5 Gyr⁻¹, i.e., keep a similar slope and is still independent of other parameters. In the modeling part below, we assume this is true. Future deeper surveys will provide critical observation evidences in the low sSFR regime. If the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation does change its slope and normalization in the low sSFR regime, or becomes dependent on other galaxies properties, the derived gas cycle during quenching should be revised accordingly.

The results shown in Figure 2 confirm the results in D21 that there exists a tight correlation between SFE${}_{{{\rm{H}}}_{2}}$ and sSFR. This relation holds from star-bursting galaxies to those in the process of being quenched, probing most accurately galaxies with log sSFR ≳ −2 Gyr⁻¹. It is independent of all key parameters that we have explored, including stellar mass, B/T, size, environment (in terms of central/satellite), and gas-phase metallicity. However, no such a tight correlation is found for H ɪ gas. The brown dashed line in all panels indicates the best orthogonal distance regression (ODR) fitted relation between SFE${}_{{{\rm{H}}}_{2}}$ and sSFR as in D21:

$$\begin{eqnarray}&&\mathrm{log}\ {\mathrm{SFE}}_{{H}_{2}}=0.50\ \mathrm{log}\ \mathrm{sSFR}\ +0.56,\end{eqnarray} \tag{ 1 }$$

where the 1σ systematic uncertainties in the two fitting parameters are 0.5 ± 0.021 and 0.56 ± 0.027, respectively. The slope of 0.5 is in good agreement with previous studies (e.g., Saintonge et al. 2012; Sargent et al. 2014; Genzel et al. 2015; Tacconi et al. 2018). For the question of whether changes in star formation are primarily driven by the change in SFE or gas fraction, it becomes immediately clear with the equation sSFR = μ × SFE. The answer is given by the logarithmic slope of the ${\mu }_{{{\rm{H}}}_{2}}$-sSFR relation, or equivalently, of the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation. Given sSFR = ${\mu }_{{{\rm{H}}}_{2}}$ × SFE${}_{{{\rm{H}}}_{2}}$, the sum of the two logarithmic slopes should be unity. As shown in D21 and Equation (1), the measured logarithmic slope is ∼0.5 for SFE${}_{{{\rm{H}}}_{2}}$–sSFR, ∼0.6 for ${\mu }_{{{\rm{H}}}_{2}}$–sSFR. Hence, their sum is indeed about unity. The two logarithmic slopes are similar (∼0.5), which means the change in log sSFR is due to similar contribution from the change in log SFE${}_{{{\rm{H}}}_{2}}$ and log ${\mu }_{{{\rm{H}}}_{2}}$, i.e., ${\mu }_{{{\rm{H}}}_{2}}$ and SFE${}_{{{\rm{H}}}_{2}}$ are comparably important in determining the star formation level (in terms of sSFR), for both quenching and star formation enhancement. Since the exact values of the two slopes may depend on sample selection, observation limit, CO conversion factor, and also fitting method, it remains to be determined more precisely with further surveys.

It should be noted that knowing SFR alone is not enough to know the star formation status of the galaxy, i.e., it is star forming or quenching. The SFR of a massive galaxy in quenching may be higher than the SFR of a star-forming low-mass galaxy. In addition to the SFR, one also needs to know the stellar mass (hence, the sSFR). The same argument also applies to ${M}_{{{\rm{H}}}_{2}}$, i.e., a higher ${M}_{{{\rm{H}}}_{2}}$ alone does not necessarily mean the galaxy is a star-forming MS galaxy, while a higher ${M}_{{{\rm{H}}}_{2}}$/M_* (i.e., μ) does, at least in a statistical sense via the observed average μ-sSFR relation. Therefore, the specific relation sSFR =${\mu }_{{{\rm{H}}}_{2}}$ × SFE${}_{{{\rm{H}}}_{2}}$ is more effective in studying the star formation status and its determining factor.

Interestingly, the sSFR depending on both SFE and μ holds also locally as derived from the resolved data, as illustrated in Section 4.2 in Morselli et al. (2020) and the various studies of the ALMA-MaNGA Quenching and Star Formation survey (Ellison et al. 2020a, 2020b). Besides, in the lower panels, all the curves are essentially below the dashed line, which indicates, as mentioned before, on average the SFE${}_{{{\rm{H}}}_{2}}$ is higher than SFE_HI.

The difference between the SFE–sSFR relation for H ɪ and that for H₂ is dramatic. It is well known that the extended Hɪ gas is very sensitive to environmental effects for satellites (Giovanelli & Haynes 1985; Catinella et al. 2013), while for massive central disk galaxies, it remains largely unchanged during quenching (Zhang et al. 2019). Therefore, different processes, in particular quenching processes (e.g., the two separated quenching channels as in (Peng et al. 2010), mass quenching and environment quenching) may have very different effects on the H ɪ gas and hence produce very different SFE_HI–sSFR relations. On the contrary, the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation remains the same for galaxies with vastly different stellar masses, structures, sizes, and in different environments. It should be noted that, this does not mean different processes (e.g., quenching, keeping galaxies evolving on the MS or starburst triggering) have zero effect on the H₂ gas. As we will show in the next subsection, different processes require different H₂ gas supply or removal rate, which controls how the galaxy evolves on the (same) SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation (i.e., positions on the relation, evolution speed, evolving to the low or high sSFR direction).

As discussed in D21, the scatter of the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is 0.2 dex, similar to the combined measurement errors of SFE and sSFR on the orthogonal direction to the fitted line which is 0.21 dex, i.e., the scatter of the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation can be entirely explained by the measurement errors. This means the intrinsic scatter of this relation is extremely small. This is supported by the fact that the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is independent of all key galaxy properties that we have explored in the previous section. This unique feature of the relation makes it simpler and more convenient in modeling the evolution of different galaxy populations than other scaling relations. This is because when galaxy evolving with time, its stellar mass increases due to star formation. If a mass-dependent scaling relation is used, for instance, the integrated SFR–${M}_{{{\rm{H}}}_{2}}$ relation (as shown in Figure 3 in D21), then different SFR–${M}_{{{\rm{H}}}_{2}}$ relations have to be used at different times (due to its strong dependence on stellar mass). Conversely, the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is stellar mass independent; hence, there is no need to consider this effect. The same argument is also applied to other galaxy properties, such as structure, size, and environment, that are all changing with time. For instance, as shown in D21 and Figure 2, the SFR–${M}_{{{\rm{H}}}_{2}}$ relation also depends on B/T and the B/T for a given galaxy is also changing with time. On the contrary, the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is independent of all other galaxy properties that we have explored, and therefore we consider it more fundamental that the SFR–${M}_{{{\rm{H}}}_{2}}$ relation.

Since the sSFR is directly related to the cosmic time t and acts as the cosmic clock (Peng et al. 2010; Peng & Maiolino 2014), the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is in fact an SFE${}_{{{\rm{H}}}_{2}}-\mathrm{sSFR}\sim t$ relation. For simplicity in modeling, we further assume that galaxies with different properties (e.g., stellar masses, metallicities, structures, sizes, and in different environments) are all strictly evolving on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation given by Equation (1). Then we use the gas regulator model to study how different galaxy populations evolve on this scaling relation.

The change in the total cold gas (i.e., H ɪ + H₂) mass of the galaxy per unit of time is given by Equation (8) in Peng & Maiolino (2014). In a similar way, given the mass conservation of H₂, the change in the H₂ gas mass is given by the net H₂ supply (or removal) rate minus the H₂ consumed by star formation, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{M}_{{{\rm{H}}}_{2}}}{{\rm{d}}t} & = & {{\rm{\Phi }}}_{{{\rm{H}}}_{2}}-{{\rm{\Psi }}}_{{{\rm{H}}}_{2}}-(1-R)\times \mathrm{SFR}\\ & = & {\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}\times \mathrm{SFR}-(1-R)\times \mathrm{SFR}\\ & = & ({\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}-(1-R))\times {\mathrm{SFE}}_{{H}_{2}}\times {M}_{{{\rm{H}}}_{2}},\end{array}\end{eqnarray} \tag{ 2 }$$

where R is the fraction of the mass of the newly formed stars as measured by the SFR, which is quickly returned to the interstellar medium and is assumed to be about 0.4. Here, we have also assumed that stars form only from H₂ gas, i.e., (1 – R) × SFR is the total mass permanently consumed by star formation and is also the H₂ mass permanently locked up into long-lived stars. ${{\rm{\Phi }}}_{{{\rm{H}}}_{2}}$ is the H₂ supply rate and ${{\rm{\Psi }}}_{{{\rm{H}}}_{2}}$ is the H₂ removal rate. Hence, ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=({{\rm{\Phi }}}_{{{\rm{H}}}_{2}}-{{\rm{\Psi }}}_{{{\rm{H}}}_{2}}$)/SFR is the net H₂ supply loading factor. It should be noted that we use the term H₂ supply/removal instead of H₂ inflow/outflow. This is because H₂ does not typically accrete onto galaxies or haloes as inflow, except for mergers, but is rather formed inside galaxies from pre-accreted H ɪ. Hence, the H₂ supply here includes H₂ formed inside galaxies, and those brought in by mergers or accreting filaments. The H₂ removal here includes the H₂ outflow driven by star formation or active galactic nuclei (AGN) feedback, and it can also include H₂ destroyed by photodissociation (see, e.g., Blitz & Rosolowsky 2006; McKee & Krumholz 2010; Morselli et al. 2021), or removed by external environment effects (e.g., strong ram pressure stripping).

As noted above, if galaxies with different properties are all assumed to evolve strictly on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation given by Equation (1), given the analytic form of Equation (2), the evolution locus of a galaxy on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation is uniquely controlled by ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$. A positive value of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ means net H₂ supply, and a negative value means net H₂ removal. In particular, when λ_net,eq = (1 − R) ∼ 0.6, ${\rm{d}}{M}_{{{\rm{H}}}_{2}}$/dt = 0, and the model achieves equilibrium at a constant H₂ gas mass. Therefore, we can use the evolution locus of a given galaxy population on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation to accurately determine ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ and the net H₂ gas supply or removal rate.

Figures 3 and 4 show the evolution locus of model galaxies with different stellar mass and ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ (labeled on the top of each panel) on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation (thick brown line). Although the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation itself does not depend on stellar mass, the sSFR of the MS do. For the MS as shown in Figure 1, the sSFR at log M_* ∼ 9, 10, and 11 is −0.7, −1, and −1.7 Gyr⁻¹, respectively. These are indicated by the vertical dashed line, labeled as "z ∼ 0 MS" in each panel of Figure 3 (for log M_* ∼ 10) and Figure 4 (for log M_* ∼ 9 and 11). As shown in Figure 3(b) in Peng et al. (2015), the stellar age difference between the star-forming and quiescent galaxies is almost constant across the entire observable stellar mass range in SDSS and is equal to Δage ∼ 4 Gyr. This is in broad agreement with other quenching timescales measured in Balogh et al. (2016), Fossati et al. (2017), and Guo et al. (2017). Therefore, in the toy model we start the evolution from the MS at z = 0.5 (i.e., 4 Gyr ago), shown by the circle dot labeled as "0 Gyr" in each panel. The observed sSFR of the MS at z = 0.5 is about 0.5 dex higher than the MS at z = 0, derived from sSFR ∼ (1 + z)³ (e.g., Lilly et al. 2013; Speagle et al. 2014; Ilbert et al. 2015).

Zoom In Zoom Out Reset image size

The initial conditions are the initial M_*, redshift and the assumed ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ (labeled in each panel in Figure 3). As mentioned above, the appearance of M_* in the initial conditions is due to the sSFR of the MS is mass dependent at a given redshift. The SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation itself is mass independent. One great feature of this approach is that there is no need to assume any initial ${M}_{{{\rm{H}}}_{2}}$ or SFR. The evolution of the model galaxy is calculated as follows. The sSFR of the MS is given by the observed sSFR_MS(M_*, z). Then SFE${}_{{{\rm{H}}}_{2}}$ is given by the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation in Equation (1). ${\mu }_{{{\rm{H}}}_{2}}$ is then given by sSFR/SFE${}_{{{\rm{H}}}_{2}}$. Then ${\mu }_{{{\rm{H}}}_{2}}$ and M_* together give ${M}_{{{\rm{H}}}_{2}}$. Putting SFE${}_{{{\rm{H}}}_{2}}$, ${M}_{{{\rm{H}}}_{2}}$, and ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ into Equation (2) gives the change of ${M}_{{{\rm{H}}}_{2}}$ during dt. Repeating this iteration gives the evolution loci of the model galaxy under different assumed gas cycle regime (i.e., different ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$) in each panel in Figure 3.

As shown in Figure 1 in Peng & Maiolino (2014), the average accretion rate of a given dark matter halo and the implied baryonic accretion rate change little in the past few gigayears (e.g., from z = 0.5–0). The H₂ supply may also have not changed much; hence, a constant ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ has been assumed here. At z > 1, both halo and baryonic accretion rate change rapidly with time; hence, a time-dependent ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ will need to be used.

As shown in Figure 3(a) and (b), with ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=0$ (no net H₂ supply or removal, like a closed box) or −0.25 (weak net H₂ removal), the galaxy evolves rapidly along the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation toward low SFE${}_{{{\rm{H}}}_{2}}$ and low sSFR. It will be fully quenched and join the quiescent galaxy population (whose SFRs are to be 1 dex below that of the MS galaxies) within 4–5 Gyr, in good agreement with the observed Δage ∼ 4 Gyr. This suggests that a strong H₂ removal is not required to quench the low-redshift MS galaxies. As shown in panel (c), with a stronger H₂ removal ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=1$, the galaxy will be fully quenched in about 2 Gyr. However, as shown in Figure 4 of Peng et al. (2015) (see also in Trussler et al. 2020), with a net gas removal of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=1$, the stellar metallicity enhancement during quenching will be significantly smaller than the observed values and hence argues against strong net outflow during quenching for middle to low-mass galaxies, on average.

As shown in Figure 3(d), with ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=0.58$ (modest net H₂ supply), the galaxy evolves slowly downwards along the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation. After about 4–5 Gyr (the cosmic time difference between z = 0.5 and 0), it arrives at the MS at z = 0. This suggests that to keep the galaxy evolving on the MS, it needs a net H₂ gas supply rate of about 0.58 SFR (i.e., ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=0.58$), which interestingly is about the equilibrium value predicted in the gas regulator model, of λ_net,eq = (1 − R) ∼ 0.6 with R assumed to be 0.4. Therefore, the MS galaxies with log M_* = 10 are indeed evolving around the equilibrium state, in which their H₂ gas mass is about constant with time. Since its stellar mass will continue to grow due to star formation, their sSFR and gas fraction will still decrease gradually with time.

We investigate the effect of the uncertainties in the fitting parameters in Equation (1) on the results by allowing the two fitting parameters to change within the estimated uncertainties. We also test our results using fitting parameters derived from ordinary least squares, which produces a slightly shallower slope of 0.44. These results are shown in Figure B1 in Appendix B. The evolution loci of the model galaxy on different lines are all similar. Hence, the uncertainties in the fitting parameters only have a small effect on our results.

Figure 4 is similar to Figure 3, but for galaxies with log M_* = 9 (upper panels) and log M_* = 11 (lower panels). As mentioned above, the difference due to a different M_* is the sSFR of the MS, which is now −0.7 Gyr⁻¹ for log M_* ∼ 9 and −1.7 Gyr⁻¹ for log M_* ∼ 11 at z ∼ 0. The case for log M_* ∼ 9 is very similar to log M_* ∼ 10 as shown in Figure 3. With no net H₂ supply, the galaxy will be fully quenched in about 4 Gyr. To keep it evolving on the MS, a net H₂ supply of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}\sim 0.6$ is required, which is also equal to the equilibrium value of λ_net,eq = (1 − R), see Equation (2), with R assumed to be 0.4. For a massive MS galaxies with log M_* ∼ 11, the situation is different. To fully quench it in 4 Gyr (i.e., within the observed Δage ∼ 4 Gyr as discussed before), the halt of H₂ supply (through either prior strangulation of predominantly H ɪ gas, and/or prevention of H ɪ-to-H₂ conversion inside galaxies) is not enough, and it requires an additional net H₂ removal with ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}=-1$ (i.e., a mass-loading factor of unity). This means the massive galaxies are more resistant to quenching due to the halt of H₂ supply. This is because massive MS galaxies have a lower sSFR than low-mass MS galaxies and hence a lower SFE or longer τ_dep. While the definition of quench is the same for all galaxies, i.e., to decrease the sSFR by more than 1 dex relative to the MS. This effect is clearly seen in other panels that the decrease of sSFR is largest during the first Gyr of quenching. For the same reason, it is also relatively easier to keep the massive galaxy stay on the MS, which requires a net H₂ supply of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}\sim 0.4$, lower than the λ_net,eq. This means the H₂ gas mass in the massive MS galaxies will slowly decrease with time.

Putting together, to quench middle to low-mass galaxies, the halt of H₂ supply is required. The observed stellar metallicity difference between star-forming and passive galaxies argues against any strong outflow or gas removal in this mass range (Peng et al. 2015; Trussler et al. 2020). Massive galaxies are more resistant to quenching due to the halt of H₂ supply. To quench them, additional H₂ gas removal is required (with a net mass-loading factor about unity), probably driven by AGN feedback in massive galaxies (e.g., Cicone et al. 2014; Förster Schreiber et al. 2019; Veilleux et al. 2020, for a review) or by photodissociation of H₂. Interestingly, this does not contradict the stellar metallicity requirement as in Peng et al. (2015). This is because although gas outflow or gas removal will act to reduce the stellar metallicity enhancement during quenching, the f_gas of the MS massive galaxy on average is already very low and the stellar metallicity enhancement will be very limited even if the galaxy is quenched with zero gas removal. Therefore, counterintuitively, although the low-mass galaxies have a significantly higher f_gas in H ɪ than in massive galaxies (as evident in Figure 2(g)), they are more sensitive to the H₂ gas supply and are easier to be quenched by the halt of H₂ supply (through either prior strangulation of predominantly H ɪ gas due to environmental effect, and/or prevention of H ɪ-to-H₂ conversion inside galaxies).

Zoom In Zoom Out Reset image size

It is also interesting to notice that for the cases of quenching (i.e., with no H₂ supply or net H₂ removal, including Figures 3(a)–(c) and Figures 4(a) and (c)), the quenching speed (dsSFR/dt) is faster in the beginning due to the higher SFE (hence, shorter τ_dep). The sSFR drops more than 1 dex in the first 1 or 2 Gyr (as in Figure 3 and Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While keeping galaxies evolving on the MS, the middle to low-mass galaxies require a net H₂ supply of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}\sim 0.6$, which is the equilibrium state predicted by the gas regulator model of λ_net,eq = (1 − R) ∼ 0.6 with R assumed to be 0.4. Massive galaxies require a slightly lower net H₂ supply of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}\sim 0.4$, which is smaller than the λ_net,eq due to their relatively lower sSFR.

It should be appreciated that in the simple framework of FFR, the value of M_* only matters due to its (weak) dependence on sSFR (at a given z). As a consequence of this, at a given z, the MS galaxies with different M_* locate on slightly different positions on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation. Apart from this, their evolution on the SFE${}_{{{\rm{H}}}_{2}}$–sSFR(t) relation is controlled solely by ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$. A zero or negative ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ (net H₂ removal) in general quenches the galaxy and makes it rapidly evolve to join the quiescent population within 4–5 Gyr; a positive ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ around λ_net,eq (modest net H₂ supply) can keep the galaxy evolve on the MS (with slowly decreasing sSFR), and a large positive ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ (strong net H₂ supply introduced by, e.g., mergers or strong dynamical interactions) can boost the sSFR (e.g., trigger starbursts).

We stress that although the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is very tight and is independent of other key galaxy properties, such as M_*, structure, metallicity, and environment, the value of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ may strongly depends on these properties. For instance, when moving from fields to groups, and to clusters, as environmental quenching becomes more and more effective, the value of ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ is expected to become progressively smaller (i.e., suppressed H ɪ gas inflow and hence less available H ɪ gas to be converted into H₂), then around zero (i.e., strangulation of H ɪ gas inflow and no more H₂ supply), and then negative (i.e., strangulation plus additional gas removal, such as gas stripping in clusters). This means that galaxies in dense regions, on average, will evolve faster along the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation to the low sSFR direction (i.e., being quenched faster). Therefore, the fact that the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is independent of environment does not contradict to the well-known fact that the fraction of the passive galaxies is higher in dense regions (Kauffmann et al. 2004; Baldry et al. 2006; Peng et al. 2010). This argument can also be applied to explain morphological quenching (Martig et al. 2009; Genzel et al. 2014; Gensior et al. 2020), i.e., the SFE${}_{{{\rm{H}}}_{2}}$–sSFR relation is independent of B/T, but the galaxies with a more massive bulge (hence, higher B/T) may have a smaller ${\lambda }_{\mathrm{net},{{\rm{H}}}_{2}}$ and hence more suppressed star formation.

As shown in Figure 2, the SFE_{H I}–sSFR relation shows strong systematic dependence on all key parameters that we have explored (to a lesser extent for B/T), in particular at the star-forming regions at log sSFR > −1.5 (where the sample has a high detection rate >80%). We note that at log sSFR < −1.5, the H ɪ detection rate drops rapidly below 20% (Figure A2). This is evident in Figures 2(g)–(j) that the differences between different curves become smaller at lower sSFR, which are likely caused by the increasing non-detections in H ɪ. The true differences could be much larger, like those at log sSFR > −1.5. Comparing panels (f) and (l), the scatter of the SFE–sSFR relation for H ɪ is significantly larger than that for H₂, in particular at log sSFR > −1.5, where both the H ɪ and H₂ sample have a high detection rate >80%. In general, given the fact that the SFE_{H I}–sSFR relation is very different for galaxies with different M_*, structures, metallicities, and in different environments (i.e., central/satellite), it is impossible to derive a simple H ɪ gas supply/removal rate as for H₂. As discussed before, this reflects that different quenching mechanisms may have very different effects on H I gas and hence produce very different SFE_{H I}–sSFR relations. We will further investigate this in the future work.