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

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 $\Lambda$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}\,\rm Mpc$, $75\,h^{-1}\,\rm Mpc$, and $205\,h^{-1}\,\rm Mpc$. In this work, we use the publicly available TNG100 data as a good balance between resolution and volume (Nelson et al., 2019b)¹¹1https://www.tng-project.org. 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. stellar 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 Friends-of-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 self-consistently 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, simulation-based, 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 chemical-network-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 $\rm H\,{\textsc{i}}$ and ${\rm 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. (2018, 2019); Stevens et al. (2019b, 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_{\star}$ for each galaxy is calculated within the radius of $2R_{\star}$, where $R_{\star}$ is the stellar half-mass radius. To approximate the scale of the Arecibo beam width ($3.4^{\prime}$ at $21$ cm, e.g. Jones et al. 2018) for the relevant galaxies in the ALFALFA survey, we only count the $\rm H\,{\textsc{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 $\rm 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 $\rm 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_{\star}/{\rm M}_{\odot}<10^{11}$. We use the $\kappa_{\rm 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 $\kappa_{\rm rot}>0.5$, the galaxy is rotation-dominated, disc-like; if $\kappa_{\rm 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 $\kappa_{\rm 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\,{\rm Myr}$ within two different apertures: everything that are gravitation-ally self-bounded within the SUBFIND identified subhalo, SFR(all grav.), and that within $2R_{\star}$, SFR($<2R_{\star}$). 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\,{\rm Myr}$, this SFR resolution limit is $10^{-2.46}\,{\rm M}_{\odot}\,{\rm yr^{-1}}$. For galaxies with SFRs below this limit, we reassign each of them a random SFR between $10^{-2.8}$–$10^{-2.5}\,{\rm M}_{\odot}\,{\rm yr^{-1}}$ so that they show up in the log-scale scatter plots.

L-Galaxies 2020 is the newest version of the Munich semi-analytic galaxy formation model (Henriques et al., 2020)²²2https://lgalaxiespublicrelease.github.io/. This model is built upon the subhalo merger trees from the Millennium and Millennium-II simulations scaled to first-year Planck cosmology. The evolution 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_{\star}/{\rm M}_{\odot}<10^{11}$ and disc-to-total stellar mass ratios (D/T) $>0.5$ from the publicly available online database³³3http://www.mpa-garching.mpg.de/millennium. 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)⁴⁴4http://eagle.strw.leidenuniv.nl 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 ($\Delta T_{\rm SF}=10^{7.5}$K) and BH feedback ($\Delta T_{\rm BH}=10^{8.5}$K), where $\Delta T_{\rm SF}$ and $\Delta T_{\rm 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 IllustrisTNG, 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 $\kappa_{\rm 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 (Vogelsberger et al., 2014; Sijacki et al., 2015; Genel et al., 2014; Nelson et al., 2015)⁵⁵5https://www.illustris-project.org. 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.13$cm^-3.

Figure 1 shows the neutral hydrogen (${\rm H\,{\textsc{i}}+\rm H_{2}}$) gas mass versus SFR for massive central disc galaxies ($10^{10.6}<M_{\star}/{\rm M}_{\odot}<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 relation is not very small. By contrast, galaxies in EAGLE, L-Galaxies, and Illustris are too gas-poor at fixed SFR. When $10^{-0.5}<{\rm SFR}/({\rm M}_{\odot}{\rm yr^{-1}})<10^{0.5}$, the neutral hydrogen gas decreases too fast in EAGLE, L-Galaxies, and Illustris, compared to observations and TNG100. For ${\rm SFR<10^{-0.5}\,{\rm M}_{\odot}\,{\rm 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$\alpha$/$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 ${\rm H}_{\alpha}/D_{\rm 4000}$) is used, the lowest SFRs of central disks are around log SFR=$-0.5$ whereas they go down to $\sim$ $-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_{\rm 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 LINERs). 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.