Galaxy and black hole coevolution in dark matter haloes not captured by cosmological simulations

Abstract

Star formation in galaxies is governed by internal and environmental processes, yet their relative roles are not well understood. In particular, uncertainties in measurements of active galactic nuclei (AGN) host galaxies, combined with modeling limitations, obfuscate the impact of supermassive black hole feedback across environments and over time. Here we address this with a comprehensive analysis of ~60,000 nearby AGNs (z < 0.15 and new environment and halo-mass measurements for ~500,000 AGN and non-AGN host galaxies. This benchmark enables unified comparisons with three prominent cosmological simulations--SIMBA, TNG, and EAGLE--and reveals major, contrasting shortcomings. Simulations fail to reproduce observed trends linking star formation, quiescence, AGN luminosity, stellar mass, and halo mass. While simulations qualitatively capture that AGNs are more common in low-mass halos than in rich groups or clusters, detailed host demographics diverge strongly from observations. Partial agreement exists in the stellar mass distribution within large-scale structures, yet all simulations overproduce quenched low-mass satellites in massive halos, while misrepresenting quenched fractions of massive central galaxies and those in low-density environments, which are sensitive to feedback implementation. Improved AGN physics and modeling of multi-phase gas cooling and flows are required to capture the observed interplay between black holes, galaxies, and halos.

Figures (29)

• Figure 1: Galaxy stellar mass functions. Stellar mass functions (SMFs) of galaxies from the Sloan Digital Sky Survey (SDSS), the Galaxy And Mass Assembly survey (GAMA), and three cosmological simulations (TNG, SIMBA, and EAGLE) at $z \sim 0.1$ are shown. Panel (a) presents the SMF for all galaxies, while panels (b)–(d) show SMFs for galaxy subpopulations ranked by decreasing star-formation rate: high star-forming galaxies (b), low star-forming galaxies (c), and quiescent galaxies (d). Error bars denote standard deviations estimated assuming Poisson statistics. The simulations are calibrated to match the total SMF in panel (a).
• Figure 1: Stellar mass functions of the SIMBA, EAGLE, and TNG simulation variants compared with GAMA. Panels are grouped by simulation (snapshots at $z \approx 0.1$) and, for each simulation, show the stellar mass functions of the global, high-SFR, low-SFR, and quiescent populations, defined by $\Delta \langle \log(\mathrm{SFR})\rangle_{\mathrm{MS}}$ relative to the star-forming main sequence (SFMS). In all panels, the violet curves show the median stellar mass function and the 16th–84th percentile uncertainty range derived from fits to GAMA galaxies at $z < 0.12$. Error bars on the simulation curves indicate standard deviations assuming Poisson statistics. SIMBA (a–d): Black, red, green, and blue curves correspond to SIMBA variants with X-ray feedback disabled; both jet-mode and X-ray AGN feedback disabled; all AGN feedback disabled; and all AGN and supernova feedback processes disabled, respectively. EAGLE (e–h): Red, teal, green, and blue curves correspond to the fiducial EAGLE-50 simulation; the EAGLE-50 variant with increased AGN heating temperature (TAGN $= 10^{9}$ K); the EAGLE-50 simulation with AGN feedback disabled; and the higher-resolution EAGLE-25 simulation, respectively. TNG (i–l): Blue, green, and red curves show the three volume variants of the TNG simulations.
• Figure 1: Star formation rates (SFRs) and stellar masses of observed and simulated galaxies. Contour levels correspond to 5, 15, 25, 50, 75, 85, and 95 percent of the data. Panel (a) shows the completeness-corrected SDSS galaxy distribution (navy) above $M_\star \gtrsim 10^9\, M_\odot$, alongside the GAMA distribution (violet). Panels (b-d): Distributions from the TNG, SIMBA, and EAGLE simulations, respectively. The solid blue and red lines in each panel indicate the mean loci of SDSS star-forming and quiescent galaxies to facilitate comparison. To approximate observational uncertainties, SFR values for quiescent galaxies (QGs) in the simulations have been randomly jittered. The dotted black lines mark the peak location of QGs in SDSS. Note the horizontal offset for SFGs and the vertical shift for massive QGs in the simulations. SIMBA's SFRs for SFGs (main sequence) are significantly offset from the observations, clustering around $M_\star \approx 10^{10}\, M_\odot$ where many galaxies are transitioning from SFGs to QGs. Additionally, SIMBA’s massive QGs peak at lower $M_\star$ than SDSS and TNG, potentially due to differences in black hole mass thresholds for kinetic feedback activation. EAGLE underpredicts the number of massive QGs. These trends are further quantified with the mass function (Fig. 1 and Supplementary Fig. \ref{['fig:SMF_SFG_GV']}).
• Figure 2: Comparison of stellar mass overdensity distributions across multiple scales. Panels (a)–(d) show the stellar mass overdensity distributions measured on scales of 0.5, 1, 2, and 8 Mpc/$h$, respectively, for galaxies at $z<0.12$ in SDSS and GAMA, compared with predictions from the TNG, SIMBA, and EAGLE simulations.
• Figure 2: Quiescent satellite fractions in SIMBA and EAGLE feedback variants. Panels (a)–(d) show the SIMBA50 variants: full physics with all feedback enabled; both jet-mode and X-ray AGN feedback disabled; all AGN feedback disabled; and all AGN and supernova feedback disabled, respectively. Panels (e)–(h) show the EAGLE50 variants: full physics; increased AGN heating temperature (TAGN = $10^{9}$ K); AGN feedback disabled; and the larger-volume EAGLE100 simulation, respectively.