Cold Gas in Massive Galaxies as A Critical Test of Black Hole Feedback Models

Black hole feedback has been widely implemented as the key recipe to quench star formation in massive galaxies in modern semi-analytic models and hydrodynamical simulations. As the theoretical details surrounding the accretion and feedback of black holes continue to be refined, various feedback models have been implemented across simulations, with notable differences in their outcomes. Yet, most of these simulations have successfully reproduced some observations, such as stellar mass function and star formation rate density in the local Universe. We use the recent observation on the change of neutral hydrogen gas mass (including both ${\rm H_2}$ and ${\rm HI}$) with star formation rate of massive central disc galaxies as a critical constraint of black hole feedback models across several simulations. We find that the predictions of IllustrisTNG agree with the observations much better than the other models tested in this work. This favors IllustrisTNG's treatment of active galactic nuclei - where kinetic winds are driven by black holes at low accretion rates - as more plausible amongst those we test. In turn, this also indirectly supports the idea that the massive central disc galaxy population in the local Universe was likely quenched by AGN feedback.


INTRODUCTION
Baryons cool and form stars within the potential well of dark matter halos (Rees & Ostriker 1977;White & Rees 1978;White & Frenk 1991). However, only 5%-25% of baryons within dark matter halos are converted into stars by z = 0 efficiently in most galaxies (Wechsler & Tinker 2018;Dutton et al. 2010;Zheng, Coil & Zehavi 2007;Moster, Naab & White 2018;Wang & Jing 2010;Yang, Mo & van den Bosch 2009;Yang et al. 2012;Kravtsov, Vikhlinin & Meshcheryakov 2018;Behroozi et al. 2019). Stellar feedback, such as supernova (SN) explosions and stellar winds, can heat up and even eject the gas in low-mass systems, reducing their star formation (SF) efficiency (Larson 1974;White & Rees 1978;Silk 2003;Springel & Hernquist 2003). However, to suppress the SF efficiency in massive systems where the potential well is deeper, SN feedback alone is not enough (Benson et al. 2003). Active galactic nuclear (AGN) feedback is invoked as a necessary mechanism to reduce further SF activity in massive systems in galaxy formation and evolution models (Silk & Rees 1998;Croton et al. 2006;Bower et al. 2006;Henriques et al. 2015;Yuan & Narayan 2014;Yuan et al. 2015Yuan et al. , 2018, see Somerville &Davé 2015 andOstriker 2017 for a review on the current status of galaxy formation and evolution models).
AGN feedback is launched in quite different ways in different theoretical models of cosmological simulation and semi-analytic models (SAMs). In general, the way that feedback is launched depends on the BH accretion rate: the high accretion rate mode and low accretion rate mode. In SAMs (Croton et al. 2006;Bower et al. 2006;Guo et al. 2011;Henriques et al. 2015), BHs continue to accrete gas from the circumgalactic medium, the induced feedback injects thermal energy into the gas in DM halo to effectively reduce the hot-gas cooling rate. In Illustris (Vogelsberger et al. 2013;Torrey et al. 2014) and IllustrisTNG (Weinberger et al. 2017;Pillepich et al. 2018b), when the BH accretion rate is higher than a certain threshold, the thermal energy is inserted into the surrounding gas within galaxies; when the BH accretion rate is lower, a hot bubble is injected (in Illustris) or a certain amount of kinetic energy is added to the gas surrounding the BH (so called "kinetic feedback", in IllustrisTNG). In EAGLE, a single heating mode from BH is launched, which nevertheless mimics the relatively quiescent 'radio mode' and vigorous 'quasar mode' when the BH accretion rate is much smaller than or similar to the Eddington rate (Schaye et al. 2015;Crain et al. 2015). Despite the non-negligible differences among these implementations, they are all adjusted to reproduce the z = 0 stellar mass function and the stellar mass-BH mass scaling relation to various levels of agreement. More observations, especially those that are not implemented for model calibration, are needed to constrain the AGN feedback models in cosmological simulations and SAMs.
One promising constraint comes from investigating the gas content in galaxies. AGN feedback can act either negatively or positively (Cresci & Maiolino 2018). Negative AGN feedback is believed to suppress the SF in massive galaxies in two ways: through prevention (Cresci et al. 2015) and ejection. Preventative feedback is when the hot gas can not cool efficiently due to the heating from AGN, while ejective feedback is when cold gas is pushed away from the galaxy due to the wind produced by the accreting BH (Zinger et al. 2020). The role of AGN feedback being preventive or ejective is not necessarily directly tied to how the feedback is implemented in practice. For example, the kinetic mode in IllustrisTNG can be also preventive, as kinetic energy can also transform into thermal energy via shocks (Weinberger et al. 2017); likewise, the thermal injection in EAGLE can produce gas ouflows (i.e. can be ejective too) due to the induced pressure gradients. By contrast, and of lesser importance, positive AGN feedback is when additional star formation is induced by the compression of molecular gas in a galaxy's disc (Silk 2013) or directly in the outflowing cold gas (Ishibashi & Fabian 2012;Zubovas et al. 2013;Maiolino et al. 2017;Gallagher et al. 2019). Understanding how AGN feedback impacts galaxies' gas and star formation properties can help us to distinguish/constrain various implementations in these models (Terrazas et al. 2020).
On average, passive galaxies unsurprisingly have less cold gas than star-forming ones (Saintonge et al. 2016;Tacconi et al. 2018;Catinella et al. 2018). However, when massive central disc galaxies (10 10.6 < M /M < 10 11 ) are selected, i.e. focusing on internal quenching mechanisms and excluding external environmental ones, Zhang et al. (2019) (Z19 hereafter) found that, as SFR decreases, their H i gas mass remains surprisingly constant, but both H 2 gas mass and H 2 star formation efficiency decrease. Zhang et al. (2021) ((Z21 hereafter)) further show the change of gas content is also accompanied by the rapid increase of stellar concentration index, bulge-to-total mass ratio, central velocity dispersion and AGN frequency (which in Z19 are mostly LINERs). These altogether suggest more massive black holes, and possibly stronger AGN feedback, in quenching or quenched galaxies. These results are robust against different SFR estimators, including the Hα/D 4000 -based SFR (Brinchmann et al. 2004) with aperture correction by performing spectral energy distribution (SED) fitting to the photometry outside the fiber, and SFRs obtained from SED fitting of UV, optical, and mid-IR bands (Salim et al. 2016;Salim, Boquien & Lee 2018), see more detailed discussion in Z21. Since cold gas is the fuel of star formation and can provide direct evidence of how quenching may happen in galaxies, this observational result of massive central disc galaxies can be used as a new constraint for current galaxy formation and evolution simulations/models, especially for AGN feedback. In particular, since H 2 is mainly distributed in the inner stellar disc and H i is usually located at larger radii, the observed H i and H 2 versus SFR relations as in Z19 potentially put a strong observational constraint on the strength of AGN feedback and on its "inside-out" nature, as it needs to clear out most of the H 2 gas in the inner disc while retaining much of the H i gas in the outer disc.
In this work, we focus mainly on the prediction for the neutral hydrogen gas-SFR relation of massive central discs produced by IllustrisTNG, EAGLE, L-Galaxies, and Illustris. We show that IllustrisTNG agrees with the observed gas-SFR relation better than the others. We then study the gas distribution and quenching mechanism of the central discs in IllustrisTNG in more detail, as it may give useful clues on the quenching process in the real Universe.

The IllustrisTNG simulation
The IllustrisTNG project (Springel et al. 2018;Marinacci et al. 2018;Naiman et al. 2018;Pillepich et al. 2018a;Nelson et al. 2018) is a suite of cosmological hydrodynamical simulations in a ΛCDM Universe run with the moving-mesh code arepo (Springel 2010). It contains simulations of three different volumes, TNG50, TNG100, and TNG300 with respective box lengths of 35 h −1 Mpc, 75 h −1 Mpc, and 205 h −1 Mpc. In this work, we use the publicly available TNG100 data as a good balance between resolution and volume (Nelson et al. 2019b) 1 . The numerical resolution of TNG100 is also closest to the test runs on which the free parameters of the subgrid physics were calibrated, which results in the best match with observational constraints (e.g. stel-lar mass function etc.). The IllustrisTNG galaxy formation and evolution model includes gas cooling and heating, star formation, stellar evolution and chemical enrichment, SN feedback, BH growth, AGN feedback, and cosmic magnetic field (see Weinberger et al. 2017;Pillepich et al. 2018b for more detailed information on the models).
Galaxies in TNG100 are identified using the Friendsof-Friends (FOF) and SUBFIND algorithm (Davis et al. 1985;Springel et al. 2001). The neutral gas fraction of non-star-forming gas cells is calculated selfconsistently in the simulation by computing the cooling rate and photo-ionization rate due to the UV background, while star-forming gas cells are modelled as a two-phase medium following Springel & Hernquist (2003), where the hot phase gas is assumed fully ionized and the cold phase gas fully neutral (see appendix A1 of Stevens et al. 2019b for more information). In each gas cell, the molecular hydrogen fraction can be obtained by post-processing using empirical, simulationbased, or theoretical prescriptions. In this work, we use the molecular hydrogen gas fraction calculated using five different models (Diemer et al. 2018): the empirical model where the molecular gas fraction is found to be correlated with the mid-plane pressure (Leroy et al. 2008, L08 hereafter), the high-resolution chemicalnetwork-inclusive simulations that produced calibrated relations between molecular gas fraction and surface density, metallicity, and UV field (Gnedin & Kravtsov 2011, GK11 hereafter;Gnedin & Draine 2014, GD14 hereafter), and analytical equilibrium models of molecular clouds with detailed calculations of molecular hydrogen creation and destruction (Krumholz 2013, K13 hereafter;Sternberg et al. 2014, S14 hereafter). Galaxies are rotated into a face-on position using the angular momentum of gas and stars, and the H i and H 2 fractions are computed using the projected quantities (such as surface density, SFR surface density etc.) That is, we use the 'map' methods of Diemer et al. (2018). We refer the readers to Diemer et al. (2018Diemer et al. ( , 2019; Stevens et al. (2019bStevens et al. ( ,a, 2021 for more details on the decomposition methods and comparison with observation (also see Davé et al. 2020 for an additional comparison with other hydrodynamic simulations).
In this work, M for each galaxy is calculated within the radius of 2R , where R is the stellar half-mass radius. To approximate the scale of the Arecibo beam width (3.4 at 21 cm, e.g. Jones et al. 2018) for the relevant galaxies in the ALFALFA survey, we only count the H i gas within a fixed radius of 70 kpc (i.e. the corresponding physical scale of the beam witdth at z = 0.037) as in Diemer et al. (2019) and Bahé et al. (2016); but see Stevens et al. (2019b) for a more precise method. We also test our results with fixed radii of 50 kpc and 80 kpc and see no significant variations of our results. We also count the H 2 gas within the radius of 70 kpc and test the results using the radius of 20 kpc (which roughly corresponds to the IRAM aperture for xCOLD GASS survey). Our results stay unchanged since H 2 is mainly located in the central region of galaxies. To make a direct comparison with the observed trend of cold gas content versus SFR as in Z19 and Z21, we apply the same selection criteria of massive central disc galaxies with 10 10.6 < M /M < 10 11 . We use the κ rot parameter to distinguish disc galaxies from elliptical galaxies, defined as the ratio of the ordered rotation versus the total kinetic energy using stellar particles (Sales et al. 2012;Rodriguez-Gomez et al. 2015). If κ rot > 0.5, the galaxy is rotation-dominated, disc-like; if κ rot < 0.5, the galaxy is dispersion-dominated, elliptical-like. Figure A1 in Diemer et al. (2019) compared the elliptical galaxy fraction in TNG100 based on κ rot with the observed fraction given by Calette et al. (2018), finding a good match between simulations and observations in the stellar mass range studied in this work. In Section 3 and Appendix A we give more discussion on the morphology selection. With the above selection, our final sample in TNG100 includes 674 centrals, of which 340 are disc-like. The SFRs are obtained by averaging the star formation over the past 200 Myr within two different apertures: everything that are gravitation-ally self-bounded within the SUBFIND identified subhalo, SFR(all grav.), and that within 2R , SFR(< 2R ). We take SFR(all grav.) as our default choice. In Appendix B, we also test our results with SFR calculated using various apertures and time-scales, and find that the general trends remain similar.
Due to the finite mass resolution in the simulation, galaxies' SFRs suffer a resolution limit (see Donnari et al. 2019 for a detailed discussion). For SFRs averaged in the past 200 Myr, this SFR resolution limit is 10 −2.46 M yr −1 . For galaxies with SFRs below this limit, we reassign each of them a random SFR between 10 −2.8 -10 −2.5 M yr −1 so that they show up in the logscale scatter plots.
2.2. L-Galaxies, EAGLE, and Illustris-1 L-Galaxies 2020 is the newest version of the Munich semi-analytic galaxy formation model (Henriques et al. 2020) 2 . This model is built upon the subhalo merger trees from the Millennium and Millennium-II simulations scaled to first-year Planck cosmology. The evo-lution of baryonic components are traced in the model, including hot gas atmosphere, cold interstellar gas, a gas reservoir ejected by winds, stars in the bulge, disc, and intracluster light, and central supermassive BHs. Elemental abundances (including H) are traced in a galactic chemical enrichment scheme introduced by Yates et al. (2013). In this work, we select the central galaxies with 10 10.6 < M /M < 10 11 and disc-to-total stellar mass ratios (D/T) > 0.5 from the publicly available online database 3 . The amount of hydrogen in cold gas component gives us the neutral hydrogen mass.
The EAGLE simulation (Schaye et al. 2015;Crain et al. 2015;McAlpine et al. 2016) 4 is a cosmological hydrodynamical simulation run with a smoothed particle hydrodynamics' (SPH) solver gadget3 with Planck cosmology. We use the largest-volume run, Ref-L100N1504, with box length of 100 Mpc. Halos and galaxies are identified using the FOF and SUBFIND algorithms. Stochastic, isotropic heating of gas particles is implemented for SF feedback (∆T SF = 10 7.5 K) and BH feedback (∆T BH = 10 8.5 K), where ∆T SF and ∆T BH are the temperature jump of gas particles receiving feedback energy. The neutral hydrogen fractions of gas particles are estimated in post-processing (Lagos et al. 2015;Bahé et al. 2016;Crain et al. 2017), using the fitting function of Rahmati et al. (2013), which is calibrated using TRAPHIC radiative transfer simulations (Pawlik & Schaye 2008). Davé et al. (2020) compared the atomic and molecular hydrogen content of galaxies from Illus-trisTNG, EAGLE, and SIMBA (Davé et al. 2019), we refer the readers to there for more information. Similar to TNG100, the massive central discs are selected using κ rot > 0.5. The neutral hydrogen mass are calculated within a radius of 60 kpc, however, the gas content stays roughly unchanged even with a different radius choice of 40 kpc.
The Illustris simulation is the precursor simulation of IllustrisTNG Sijacki et al. 2015;Genel et al. 2014;Nelson et al. 2015) 5 . In this work, we use Illustris-1, which has the same initial condition, volume, and resolution as TNG100. The main difference of Illustris, compared to IllustrisTNG, is the treatment of BH accretion and feedback (as discussed in Section 1) and of galactic winds. Illustris also lacks magnetohydrodynamics. The central discs in Illustris-1 are selected in the same way as TNG100. The neutral hydrogen gas mass of a galaxy is also summed up over the gas cells within the radius of 70 kpc. An non-trivial caveat in the above models is that none of them attempts to actually resolve the neutral phases of the ISM. Thus the predicted neutral hydrogen fraction is not a real prediction based on the ISM chemistry, but a modeling related to the effective equation of state function, where star formation is initiated once the gas density n > 0.13cm −3 .

RESULTS
3.1. Total neutral hydrogen in simulations and real galaxies Figure 1 shows the neutral hydrogen (H i + H 2 ) gas mass versus SFR for massive central disc galaxies (10 10.6 < M /M < 10 11 ) in TNG100, EAGLE, Illustris, L-Galaxies, and observations (Z19). In those four models, notably different black hole feedback models were implemented.
First, all simulations except for EAGLE reproduce well the observed neutral hydrogen gas content for the galaxies with largest SFRs, while EAGLE predicts about 0.3 dex less gas than observed. With SFR decreasing, only TNG100 roughly reproduces the Z19 observational trend (which is also shown by the results in Figure 2), although the scatter around the median rela-tion is not very small. By contrast, galaxies in EAGLE, L-Galaxies, and Illustris are too gas-poor at fixed SFR. When 10 −0.5 < SFR/(M yr −1 ) < 10 0.5 , the neutral hydrogen gas decreases too fast in EAGLE, L-Galaxies, and Illustris, compared to observations and TNG100. For SFR < 10 −0.5 M yr −1 , the neutral hydrogen gas keeps decreasing rapidly in EAGLE, Illustris, and L-Galaxies, and becomes significantly below the observed values (with either the Hα/D 4000 or UV+IR SFR estimator). We also explored an additional semi-analytic model, Shark (Lagos et al. 2018), and found it behaves similarly to Illustris and EAGLE (not shown here to avoid overcrowding of lines). Since the change in gas as SFR decreases in simulations is mainly driven by the implemented AGN feedback model for galaxies in this stellar mass range, the results in Figure 1 suggest that the AGN feedback during quenching might be too efficient in heating/expelling gas from the galaxies in EAGLE, L-Galaxies, and Illustris.
There is a concern on the Z19 observation result that if different SFR estimator (UV+IR instead of H α /D 4000 ) is used, the lowest SFRs of central disks are around log SFR=−0.5 whereas they go down to ∼ −1 in the former case. So, these galaxies appear having more or less progressed towards full quenching, depending on the adopted SFR diagnostics (see Figure 1 in Cortese et al. 2020, hereafter C20). The effect of using alternative SFR estimator has been discussed in detail in Z21. The general trends of gas content and other galaxy properties with respect to SFR hold for both SFR estimators. It should also be noted that "disc" galaxies are defined in a very different way in Z19&Z21 (visual classification by Galaxy Zoo) and C20 (a threshold in B/T). These two selections result in very different populations. As shown in Z21, the visually defined disc galaxies with low SFR often have massive bulges, and many of them will be classified as ellipticals using B/T; while most S0s have been classified as ellipticals (or uncertains) in Galaxy Zoo, but will be classified as discs by B/T. This different disc selection leads to the fact that the UV+IR derived SFRs in Z19 do not extend to the SFRs as low as in Cortese et al. (2020).
One issue in observation is that the SFRs derived from the SED fitting of UV, optical and mid-IR bands (Salim et al. 2016;Salim, Boquien & Lee 2018) used in Z19 and Z21, and the UV+IR SFRs used in Cortese et al.
(2020) are not corrected for AGN contamination or by hot old stars. Z21 shows almost 100% of the massive discs with the lowest SFRs (by either SFR estimator) host LINERs. Thus the UV+IR SFRs are likely to be the upper limits and true SFRs could reach lower SFRs. One could subtract off an AGN component by doing the SED decomposition, but this is a very uncertain procedure. For these LINERs, Z19 and Z21 used the SFRs derived from D 4000 , while such derivation is calibrated with non-AGN galaxies (Brinchmann et al. 2004). Apparently, SFRs from both estimators contain uncertainties at the low SFR end (dominated by LIN-ERs). Therefore, in our comparison analysis, we use the results derived from both SFR estimators. As shown in Figure 1, the approximate agreement between TNG100 and either SFR estimator holds.

H i and H 2 contents (in TNG100)
3.2.1. MH i and MH 2 versus SFR Given that TNG100 agrees the best with the observed neutral hydrogen mass-SFR relation, in Figure 2 we show the H i and H 2 gas mass as a function of SFR for these massive central disc galaxies in TNG100, for the five different H i/H 2 decomposition models as described in Section 2. In each panel, we show the results for SFRs derived within the entire subhalo (square) and 2R (cross) to check the possible influence of varying apertures. It is evident from Figure 2 that above the SFR resolution limits (the vertical dashed line in each panel), all five of these H i/H 2 decomposition models show qualitatively similar trends. That is, during the quenching of these massive central disc galaxies, as their SFR decreases, their H 2 gas mass and star formation efficiency drop but their H i gas mass remains about constant. These trends are in general agreement with that reported by Z19. The results from S14 and L08 models match the observations remarkably well. Nevertheless, given the large systematic uncertainties in the post-processing (Diemer et al. 2018), we avoid stressing the performance of any individual model, but show the variation of all five models in the last panel in Figure 2. It is important to notice that the H 2 mass-SFR relation is more affected by the choices of aperture size. If the SFR(< 2R ) is used, in some H i/H 2 decomposition models, the H 2 gas stays roughly unchanged with decreasing SFR. This indicates the existence of a population of galaxies with a good amount of SF outside 2R in IllustrisTNG. The observed aperture corrected H α /D 4000 SFR and UV+IR SFR represent more the SFR over the whole galaxies including the outer region, thus we take SFR(all grav.) as our default SFR choice. The relations of H i and H 2 mass with SFR agree with the observations within model uncertainties. We further test our results with SFRs calculated using varying apertures and time-scales, which is shown in Appendix B.
We note that for galaxies with SFR below the resolution limit, their median H i gas mass is still large, with median M H i larger than 10 8.6 M when SFR(< 2R ) is used and larger than 10 9.6 M when SFR(all grav.) is used, but in general is lower than those with SFR above the limit in all five models. As shown in the figure, observations do not extend to such low SFR. Hence the existence of these fully quenched massive central disc galaxies with a good amount but lower H i gas mass becomes a prediction that can be observationally tested by future deeper H i surveys.
Although morphology (i.e. disk galaxies) is measured in a different way in the simulations compared to observations. In Appendix A, we show that the cold gas-SFR relation in TNG100 and EAGLE is almost independent of morphology, implying that the morphology parameter is irrelevant in interpreting the trends in Figures 1  and 2 when compared with observations. However, the independence (or weak dependence) of the H i mass on morphology clearly contradicts observations. Using the same sample selection, Z19 shows that massive central disc galaxies have an average H i detection fraction of > 90%, while the H i detection fraction of massive central elliptical galaxies is only < 20%. AGN feedback by itself does not necessarily directly change the morphology. Therefore, the disagreement between observation and simulation may be due to the gas having not been properly processed during morphological transformation in the simulations, such as when mergers and interactions occur.
Star formation efficiency ≡ SFR/M gas is a widely used observable quantity to explore star formation and quenching in galaxies.
(of the total neutral hydrogen gas) in EAGLE and Illustris is approximately constant at 0.5 Gyr −1 , meaning quenching in massive central galaxies in EAGLE and Illustris is primarily driven by the decrease of cold gas. Similarly, in L-Galaxies, is approximately constant for galaxies with SFR> 10 −0.5 M yr −1 and then drops rapidly for the most quenched galaxies. Hence quenching in L-Galaxies is first driven by the decrease of the cold gas (with a constant ), and then followed by the decrease of . In TNG100, as SFR drops by ∼ 3 dex, total neutral hydrogen gas mass drops by ∼ 0.5 dex and hence drops by ∼ 2.5 dex. Therefore, quenching in TNG100 is mainly driven by a decrease in the SF efficiency of the total neutral hydrogen gas. It should be noted that this conclusion relates to the total neutral hydrogen gas, and hence does not contradict the result regarding H 2 gas shown in Figure 2 (and also in Dou et al. 2021b,a) that quenching is driven by the decrease of both H 2 gas mass and H 2 star formation efficiency H2 . It is straightforward to show that H i+H2 = H2 /(1 + M H i /M H2 ).
(1) The decrease of H2 is smaller than that of H i+H2 , and the exact amount depends on M H i /M H2 . Meanwhile, we stress that in TNG100 the instantaneous SFR in the modeling is calculated based on the local gas density, not based on the H 2 gas density (which is calculated in post-processing). The dramatic decrease of the H i+H2 of ∼ 2.5 dex during quenching in TNG100 is caused by the decrease of the gas density likely as a consequence of black hole feedback. Note that in TNG100, when the gas density falls below the minimum threshold (i.e. 0.13 cm −3 ), star formation halts, by construction in the simulation.
We suspect that the large differences in the cold gas content and during quenching are mainly caused by the different black hole feedback models implemented in IllustrisTNG (kinetic winds), EAGLE (thermal feed-back), Illustris (thermal bubble model) and L-Galaxies (preventive feedback). In many semi-analytic models, such as L-Galaxies, when the average heating rate from BH feedback exceeds the gas cooling rate in the halo, cold gas accretion stops and star formation is quenched. The decrease of the cold gas amount during quenching in central discs in L-Galaxies is caused by the heating from the BH, which reduces the cooling rate of the gas. In EAGLE, the BH feedback heats up the ambient gas by a temperature increment of 10 8.5 K; this thermal injection creates a pressure gradient that can produce an outflow. However, this feedback mode might be overly aggressive in heat up cold gas (Davé et al. 2020;Mitchell et al. 2020), resulting in a deficiency of the neutral hydrogen gas. The overall energy injected by BH feedback is similar in IllustrisTNG and Illustris, but the imple-mentation is very different (Pillepich et al. 2018b). The thermal bubble model (Sijacki et al. 2007) in Illustris expels very large amounts of gas from the halo in massive galaxies. The improved kinetic winds (Weinberger et al. 2017) in IllustrisTNG still remove large quantities of gas from galaxies during quenching (see related discussions in Stevens et al. 2019bStevens et al. , 2021Terrazas et al. 2020), but it retains more gas on average, as shown in Figure 1 and Figure 3. We will discuss the effect of BH feedback in IllustrisTNG on gas properties more later.

Gas distribution profiles
Given the fact that TNG100 approximately reproduces the average cold gas mass-SFR relation as shown in Figures 1 and 2, we further explore in TNG100 the cold gas distributions in galaxies with higher SFR(all grav.)> 10 −0.5 M yr −1 and lower SFR limit <SFR(all grav.)< 10 −0.5 M yr −1 , to understand the cause of these relations. In the analysis here we exclude the galaxies with the SFR below the resolution limit. Current observations in cold gas do not extend to such low SFRs, hence the gas properties in this regime remain a prediction from IllustrisTNG to be tested by future deeper surveys. Figure 3 shows the median surface density profile and the median cumulative radial distributions of the total neutral hydrogen gas (H i + H 2 ), H i, and H 2 for galaxies with higher and lower SFRs in TNG100. The two vertical dashed lines mark radii of 10 kpc and 70 kpc to guide the eye. Comparing the total neutral hydrogen gas distribution of galaxies with higher and lower SFR, the difference between the blue and red dashed line is ∼ 0.9 dex (∼ 0.7 dex) in the inner regions at ∼ 10 kpc, and is 0.2 dex (∼ 0.5 dex) in the outer regions at ∼ 70 kpc for the surface density and cumulative mass. The H 2 surface density profile and cumulative mass profiles of galaxies with SFR(all grav.)> 10 −0.5 M yr −1 is similar to that in H i, in particular in the inner regions within 30 kpc. By contrast, the H 2 gas mass of galaxies with SFR(all grav.)> 10 −0.5 M yr −1 is significantly lower than their H i at all radii. The average H i gas within 70 kpc is ∼ 10 10 M for lower SFR galaxies and ∼ 10 9.65 M for higher FSR galaxies. This difference (∼ 0.35 dex) is smaller than the difference in the total neutral hydrogen gas, which is ∼ 0.5 dex at this radius. This makes sense in light of the flat M H i -SFR relation shown in Figure 2.
The significant decrease of total neutral hydrogen in the central regions of lower SFR galaxies compared to that of the higher SFR galaxies is likely due to the BHdriven kinetic feedback in IllustrisTNG (Terrazas et al. 2020). The kinetic winds in IllustrisTNG are imple-mented by giving the gas immediately surrounding the black hole a momentum kick in a random direction away from the black hole (Weinberger et al. 2017). The kinetic wind this generates is isotropic at the injection scales; it impacts the gas not only in the bi-conical regions, but also in the disc. This kinetic wind quenches galaxy SF by both pushing gas out of galaxies (i.e. ejective) and also heating up the gas (i.e. preventive), more so in the central regions (Zinger et al. 2020;Terrazas et al. 2020;Truong et al. 2020). Since H 2 is mainly located in galaxies' centers, it is the most affected gas phase. While central H i is also decreased, the majority of H i-which lies in the outskirts of galaxies-is less affected. The deficiency of cold gas in the central regions of the lower SFR galaxies is in agreement with Z19's observations that the central regions of the disc galaxies with lower SFR have Hα/Hβ close to the intrinsic value of 3.1, which indicates that there is little dust and gas there. Hence Figure 3 is in fact an important prediction from TNG100, which can be observationally tested with spatially-resolved H i and H 2 observations for both star-forming galaxies and those in the process of being quenched.
Interestingly, at any given radius, the H 2 gas mass difference between lower SFR and higher SFR galaxies is significantly larger than their H i gas mass difference, indicating that the quenching mechanism, i.e. BH feedback, affects H 2 more (at the same radius). One important result here is that even in the inner regions of most lowest SFR galaxies, there is a lower but still significant amount of cold gas. Therefore, the quenching (i.e. decreasing of SFR) is mainly caused by the fact that most of this neutral hydrogen gas is in the form of non-starforming H i, not H 2 , which leads to the very low of the total neutral hydrogen gas as discussed above. As shown in Diemer et al. (2018), for all of the five H i/H 2 decomposition models, the fraction of neutral hydrogen that is H 2 depends steeply on the gas density in a certain range (see also Morselli et al. 2020). The removal of the highdensity gas due to feedback can have a stronger effect on H 2 than H i, leading to the low H 2 -to-H i ratio, as shown here in Figure 3 and also in Figure 2. In Appendix C, we show the stellar light and total gas distribution for two typical central disc galaxies in TNG100, one with SFR(all grav.)> 10 −0.5 M yr −1 and one with SFR(all grav.)< 10 −0.5 M yr −1 , to give a more direct view on the above discussion.

DISCUSSION AND SUMMARY
Recent observations indicate the ubiquitous existence of large regularly-rotating H i discs in massive central disc galaxies, both in those that are star-forming and The surface density (face-on) and cumulative radial mass distribution of neutral hydrogen (dashed lines), H i (shaded area), and H2 (netted shaded region) for higher SFR (blue) and lower SFR (red) central disc galaxies in TNG100, as indicated by the label. The shaded regions indicate the range of median profiles of five H i/H2 decomposition models. The vertical gray dashed lines correspond to radii of 10 kpc and 70 kpc, and are plotted to guide the eye.
in those in the process of being quenched (Z19). This result has yet been used to calibrate the recipes in the hydrodynamical simulations and semi-analytic models, and can therefore serve as a new observational test of these simulations, in particular of their AGN feedback models, as the existence of an H i disc at high stellar masses requires the strength and modes of AGN feedback to be just right. If feedback is too violent, H i discs will not survive during quenching; if it is too weak, it will not be able to deplete the H 2 gas in the galaxy and to quench the star formation.
Among the simulations we tested, only TNG100 of the IllustrisTNG project appears to approximately reproduce the trend in cold gas vs SFR as reported in Z19 and some H 2 /H i decomposition method may even reproduce the observed H 2 and H i trends, separately. This lends some support to the AGN feedback implementation in TNG100, though other factors may also be at play. Given the coarse resolution of cosmological hdyro-simulation relative to the high-resolution simulation of individual galaxy that can resolve Bondi radius (Yuan et al. 2015(Yuan et al. , 2018 and the sub-grid nature of the the AGN feedback models implemented, it is difficult to tell too much details of the AGN feeback physics. However, reproducing the observed gas content of galaxies of various properties will be critical for future development in cosmological simulations. Also, the fact that, none of the models tested in this work reproduce the morphology-HI gas mass relation challenges the current models. In IllustrisTNG, the kinetic winds drive outflows that push gas out of galaxies ) and increase the gas entropy (Zinger et al. 2020), more so in the central regions, while the H i gas at larger radii ( 30kpc) is less affected. This explains the relatively flat M H i -SFR relation shown in Figure 2. Meanwhile, the kinetic feedback also leads to a sharp decrease in gas density, which then leads to a sharp decrease in the fraction of molecular hydrogen computed in the decomposition models (cf. Stevens et al. 2021). This explains the rapid decrease of H 2 gas mass and SFR during quenching. The simultaneous role of being ejective and preventive for AGN feedback in IllustrisTNG also implies that cold gas accretion during quenching is strangulated. These results are also in qualitative agreement with those derived from the observed SFR/M -M H2 /M -SFR/M H2 scaling relations (Dou et al. 2021b), in that to quench massive galaxies, strangulation plus additional H 2 gas removal (with a mass-loading factor of about unity) are required (Dou et al. 2021a).
It should be noted that the existence of the observed regularly rotating H i disk around the massive central disk galaxies that are undergoing quenching process as in Z19 & Z21, can be explained by AGN feedback as above. Alternatively, it can be caused by angular momentum quenching as proposed in Peng & Renzini (2020) and Renzini (2020) that the inflowing gas with excess angular momentum can settle on a stable outer neutral hydrogen disk for a long timescale in the absence of perturbation. This is supported by recent studies in simulations which found that a sufficiently high CGM angular momentum is an important factor in keeping a galaxy quenched (Lu et al. 2021;Wang et al. 2022). Although AGN feedback must be important to quench massive galaxies in simulations (e.g. Su et al. 2019) and can well reproduce the observed cold gas content during quenching as shown in this paper, it may not be the case in the real Universe. We will further investigate both AGN and non-AGN solutions to the quenching of massive galaxies in our future work.

A. COLD GAS-SFR RELATION FOR CENTRALS IN TNG100 AND EAGLE
A major concern for the study of the neutral hydrogen gas-SFR relation is the definition of the galaxy morphology; how we choose what constitutes a disc can affect our results. There are various morphology metrics used for both observations and simulations that allow one to distinguish between discs and ellipticals. However, exploring which one is better is not the task of this work. Here we briefly assess the H i gas-SFR relation for ellipticals in TNG100 and EAGLE. As shown in Figure A.1, the H i mass for ellipticals is slightly lower but not significantly different from that of discs in TNG100. This is mainly caused by the higher-than-observed H i content of TNG100 ellipticals, as pointed out by Diemer et al. (2019). The independence of the cold gas content at given SFR with morphology is also shown in EAGLE, as seen in Figure A.2. These results indicate that the neutral hydrogen gas mass-SFR relation shown in the main text is not significantly affected even if the morphological selection of galaxies in the simulations is not perfect.

B. SFR OVER DIFFERENT TIME-SCALES AND APERTURES
In observations, different SFR indicators, such as UV continuum, Hα emission/D 4000 , IR continuum and 1.4GHz emission, are sensitive to the SFR over different timescales, varying from less than ∼ 10 Myr for Hα to ∼ 100 Myr for UV+IR (Speagle et al. 2014). In simulations, SFRs can also be calculated using different time-scales and apertures. To test the impact of time-scale and aperture choices, in Figure B.1, we plot the gas-SFR relation of the L08 model as in Figure 2, but using SFRs derived within different timescales (from instantaneous to 200 Myr) and apertures (< 2R , < 30 kpc, all grav.) (Donnari et al. 2019). There are differences and variations of the results in different panels, but the general trends for both H i and H 2 remain similar, and are in general agreement with observations. We also tested the results derived with other H i/H 2 decomposition models, and find that the general trends remain unchanged against SFRs calculated using different apertures and time-scales. We want to point out that while the time-scale and aperture choices for the simulation are not particularly important, the SFR indicators used for the observations EAGLE, κ rot > 0.5 EAGLE, κ rot < 0.5 Figure A.2. Neutral hydrogen gas mass versus SFR for central discs and ellipticals in EAGLE. The left-and right-hand panels respectively show the morphology divisions using disc-to-total mass ratio (D/T) and rotation support parameter (κrot). The blue and oranges lines are median relations for discs and ellipticals separately, together with the 16th-84th percentile ranges shown by the shaded areas.
is still very important, as the modeling uncertainties in observations are more complicated than just time-scales and apertures (for example, see the difference between H α /D 4000 and UV+IR SFRs as highlighted in Z19).

C. GAS DENSITY IN EXAMPLE SF AND QUENCHING GALAXIES
In Figure C.1, we show the stellar light and total gas distribution for a typical higher SFR central disc galaxy and a lower SFR central disc in TNG100. The lower SFR galaxy also shows a regular rotating gas disc, similar to the higher SFR one. This is also in good agreement with the results in Z19, where they find that most of the quenching central disc galaxies show characteristically symmetric double-horn H i profiles, indicating regularly rotating H i discs with little significant kinematic perturbations, similar to the star-forming ones. The difference between the two galaxies is more directly shown in the right panel. Comparing to the one with higher SFR, in the lower SFR galaxy, there are fewer regions with high gas density (hence less gas is in the form of H 2 ). The dashed line corresponds to density threshold of 0.13 cm −3 (i.e. n H ∼ 0.1 cm −3 ) set up for star formation. Note the gas cell volumes in quenching galaxy is much larger than the SF one, as a result, the surface density profile of these two galaxies are not significantly different.