© The Authors 2023

1. Introduction

Understanding the processes governing the growth of supermassive black holes (BHs) constitutes a major step toward the building up of a theoretical framework for galaxy formation. In all current ab initio galaxy formation models, this growth results from the interplay between the growth of dark matter (DM) halos, the collapse and cooling of gas in the DM potential wells, the dynamical processes affecting the distribution of gas (merging, interaction, and gravitational instabilities), and the counteracting effects of outflows (e.g., Granato et al. 2004; Menci et al. 2005, 2014; Hopkins et al. 2006; Lapi et al. 2006; Shankar et al. 2006; Guo et al. 2011; Fanidakis et al. 2012; Habouzit et al. 2021, and references therein). The latter are due both to the supernovae explosions following star formation (see, e.g., White & Rees 1978; Dekel & Silk 1986) and to the radiation field of AGN powered by the process of accretion onto the central BHs.

A putative coevolution between central BHs and their host galaxies and dark matter (DM) halos is supported by the observational evidence of a series of correlations between the BH mass M_BH and several host galaxy properties (see, e.g., Ferrarese & Ford 2005; Kormendy & Ho 2013; Graham 2016, for reviews). In particular, an evident correlation between the mass of the BH M_BH and the stellar mass of the host galaxy bulge M_*, b, with an intrinsic scatter of ≲0.4 dex in BH mass at fixed bulge mass, has been reported several times (e.g., Magorrian et al. 1998; Marconi & Hunt 2003; Kormendy & Ho 2013; McConnell & Ma 2013; Läsker et al. 2014; Saglia et al. 2016; Davis et al. 2018; De Nicola et al. 2019; see, e.g., Graham 2016 for a review). A possibly even tighter correlation between M_BH and stellar velocity dispersion σ, with an intrinsic scatter as low as ≲0.3 dex (e.g., Gebhardt et al. 2000, but see more recent estimates by, e.g., Sahu et al. 2019a) has often been interpreted as evidence of the key role played by AGN in self-regulating both BH growth and star formation in the host galaxy via some wind/jet-driven feedback mechanism (see reviews by, e.g., Shankar 2009; Alexander & Hickox 2012; Somerville & Davé 2015), although galaxy-galaxy mergers could still preserve its tightness but might modify its shape (e.g., Graham 2023, and references therein). Analyses of the pairwise residual correlations applied to the BH-galaxy scaling relations have also revealed that stellar velocity dispersion σ appears more strongly correlated to the central BH mass M_BH than other galactic variables (e.g., Bernardi et al. 2007; Hopkins et al. 2007; Shankar et al. 2016, 2017, 2019; Iannella et al. 2021). This trend is expected in AGN feedback-driven coevolution models (e.g., Silk & Rees 1998; King 2003).

At even larger scales, central BHs continue to show evidence for correlations between the BH mass and the total dynamical mass of the host galaxy (Bandara et al. 2009), with the number of globular clusters (Burkert & Tremaine 2010), the galactic rotation velocity v_rot (Davis et al. 2019a,b; Robinson et al. 2021), and the mass of the entire host DM halo mass (e.g., Ferrarese 2002; Baes et al. 2003; Pizzella et al. 2005; Volonteri et al. 2011). A clear link between the BH mass and host dark matter halos is also evident from the clustering strength of AGN, at least at z ≲ 0.5 (e.g., Shankar et al. 2020a,b; Allevato et al. 2021; Powell et al. 2022). A correlation between M_BH and M_h has also been measured by Marasco et al. (2021), who considered both early- and late-type systems for which M_h was determined either from globular cluster dynamics or from spatially resolved rotation curves. Correlations between BHs and the larger-scale systems such as the host halos might also be secondary correlations that are ultimately induced by the underlying correlation between galaxy mass/velocity dispersion and host halo mass (e.g., Ferrarese 2002).

It is the aim of this paper to explore the predicted normalization, shape, scatter, and even pairwise residuals of the main scaling relations between BH mass and host galaxy properties, namely stellar mass, stellar velocity dispersion, and the mass of the host dark matter halo, with the goal of determining how these correlations arise from the complex interplay between DM halo growth, the inflow of gas, galaxy mergers, and the stellar-AGN feedback mechanisms at play in galaxy evolution models. Shedding light on the origin of the BH-host scaling relations has constituted a major goal of galaxy formation models in past decades, and all current cosmological models now probably include more or less refined recipes for self-consistent growth of massive BHs at the centers of galaxies, as also mentioned above (e.g., Cirasuolo et al. 2005; Fontanot et al. 2015, 2020; Habouzit et al. 2021, and references therein). However, this work differs from previous attempts in two main respects: (i) we adapted a well-tested and realistic kinetic blast-wave AGN feedback recipe to a comprehensive semi-analytic model, and (ii) we self-consistently explore the possible impact of selection biases while analyzing the shapes and residuals in the BH-galaxy scaling relations. We further discuss these points below.

A critical aspect in realistic models of the coevolution of BHs and their hosts is constituted by the modeling of AGN feedback. While the role played by the growth of DM halos and the ensuing inflow and cooling of gas in the building of massive BHs is relatively well established in galaxy formation models and simulations (see, e.g., Somerville & Davé 2015; Naab & Ostriker 2017 for a review), large uncertainties still affect the implementation of reliable descriptions of AGN feedback, which constitutes a key self-regulation mechanism in the process of BH growth (see, e.g., Alexander & Hickox 2012 for a review). Cosmological simulations adopt various different phenomenological approximations to represent AGN feedback, often including multiple feedback modes associated with rapid (close to Eddington) or slow (very sub-Eddington) accretion onto the central BH (see Somerville & Davé 2015, and references therein for reviews). The large dispersion in the predicted AGN luminosity functions from different simulations presumably largely reflects the main uncertainties in modeling these processes. On the other hand, semi-analytic models often adopt a parametric phenomenological description of AGN feedback and compute the amount of ejected gas on the basis of parametric laws based on energy-balance arguments. They do not consider the complex two-dimensional structure of the outflows. While an approach based on a physical model for the expansion of AGN-driven shocks in the galactic gas has been adopted in Menci et al. (2008), this was still based on a isotropic model for shock expansion. However, to act as an effective negative feedback on star formation and BH accretion, AGN winds must couple to and affect cold galactic gas that is rotationally supported and settled in a disk. The high density of the gas distribution in the direction parallel to the plane of the disk strongly inhibits the shock expansion in this direction, while the gas is preferentially ejected perpendicular to the disk, resulting in an overall fraction of ejected interstellar medium lower than in one-dimensional (isotropic) models (see, e.g., Faucher-Giguere & Quataert 2012; Hartwig et al. 2018). Moreover, all treatments of AGN feedback adopted so far in galaxy formation models have not been tested in detail against the observed properties of AGN winds, which have been increasingly frequently studied in the recent years. Since the works by Cicone et al. (2014) and Fiore et al. (2017), samples with more than one hundred outflow measurements have been assembled, with detected massive winds at different scales (subparsec to kiloparsec) and with different molecular/ion compositions. In a previous paper (Menci et al. 2019), we computed the two-dimensional expansion of outflows driven by AGN in galactic disks as a function of the global properties of the host galaxy and of the luminosity of the central AGN. We derived the expansion rate, the mass outflow rate, and the density and temperature of the shocked shell in the case of an exponential profile for the disk gas for different expansion directions θ with respect to the plane of the disk. We expressed our model results in terms of the global properties of the host galaxies and compared our predictions to a large sample of 19 outflows (mostly molecular, except for one object) in galaxies with a measured AGN luminosity and gas mass, and with the estimated total mass. This allowed us to test the model through a detailed one-by-one comparison with the data. The two-dimensional structure that we obtained for the outflows is characterized by a shock expansion that follows the paths of least resistance with an elongated shock front in the direction perpendicular to the disk. The higher outflow velocities attained in the direction perpendicular to the disk easily exceed the escape velocity at the virial radius on a short timescale, while the slower expansion of the shock in the plane of the disk can prevent the escape of gas in this direction within the lifetime of the AGN. The overall ejected gas fraction differs substantially from that obtained in spherical models, both in the total value and in its dependence on the galaxy properties. Here we implement this two-dimensional description in our semi-analytic model of galaxy formation. For each model galaxy in our Monte Carlo realizations, we associate the outflow solutions corresponding to its properties (gas mass, total DM mass, and AGN luminosity) and derive the opening angle generated by the outflows and the total ejected gas mass. This improved and observationally tested description of AGN feedback constitutes a novel upgrade of our galaxy formation model that allows us to compute the scaling relations between the BH mass and the different galaxy properties on more physically motivated grounds. In addition, we have improved our model by implementing the computation of the size R_e of the bulge component of the model galaxies, and then by deriving the bulge stellar velocity dispersion σ. This is computed by solving the spherical Jeans equation for each model galaxy. This allows us not only to predict (among other observables) the M_BH − σ relation, but also to derive the BH sphere of influence for each model galaxy and hence to investigate the effects of angular resolution-related selection effects.

We also explore the impact of possible observational biases here that may affect the local sample of dynamically measured supermassive BHs, and thus in turn affect the comparison between theoretical predictions and the observed correlations because the latter may differ from the intrinsic relations. Shankar et al. (2016, 2017, 2019), following in the footsteps of other groups (e.g., Batcheldor 2010, Morabito & Dai 2012), have put forward the hypothesis that local samples of quiescent galaxies with dynamically measured M_BH may be affected by an angular resolution-related selection that might bias the observed scaling relations between M_BH and host galaxy properties away from the intrinsic relations. Insufficient resolution might prevent reliable BH mass estimates whenever the spatial resolution r_crit is higher than the BH sphere of influence r_infl ≡ GM_BH/σ². Specifically, the condition r_infl ≥ r_crit may lead to a sample that is biased toward higher masses of dynamically measured M_BH at fixed galaxy stellar mass or stellar velocity dispersion. The above authors investigated the effect of this bias by assuming different intrinsic relations as an input for dedicated Monte Carlo simulations, concluding that the M_BH − M_* relation is more strongly biased than the M_BH − σ relation. Because this bias only affects dynamical measurements of M_BH, this might explain the difference in the scaling relations obtained for samples of active AGN (Ho & Kim 2014; Martin-Navarro & Mezcua 2018; van den Bosch 2016; Greene et al. 2016; Busch et al. 2014, Reines & Volonteri 2015; Bentz & Manne-Nicholas 2018). In these samples, M_BH is not determined from spatially resolved dynamical measurements, but rather from the assumed virial motions of the gas in the broad line region (BLR) orbiting in the vicinity of the central BH. However, these measurements are affected by the uncertainty related to the kinematics, geometry, and inclination of the BLR clouds (e.g., Ho & Kim 2014, and references therein), usually parameterized in terms of the free parameter f_vir that relates the measured velocity of the clouds ΔV and its size r to the BH mass M_BH = f_vir r (ΔV)²/G (but see also, e.g., Marculewicz & Nikolajuk 2020 on how to adapt the virial factor f_vir based on the spectral properties of the quasars). Although the degree and extent of a resolution bias in the current sample of dynamically measured supermassive BHs is still a matter of debate (e.g., Sahu et al. 2023), it is worth exploring its possible impact on the predicted scaling relations in a full and self-consistent galaxy and supermassive BH coevolution model (see, e.g., Barausse et al. 2017).

The paper is organized as follows. In Sect. 2 we first summarize the basic features of the semi-analytic galaxy formation model. The implementation of the above two-dimensional description of AGN outflows in the semi-analytic model is presented in Sect. 2.1, and in Sect. 2.2 we describe how we compute the bulge stellar velocity dispersion for each model galaxy. In Sect. 3 we briefly recall how the computed stellar velocity dispersion allows us to study angular resolution-related selection effect affecting the comparison between the observed and the predicted scaling relation between the dynamically measured BH mass and the galaxy properties. In Sect. 4 we present our results for the local scaling relations between the BH mass and the bulge and the global properties of the host galaxies, and we present our predictions for the redshift evolution of the M_BH − M_* relation. Section 5 is devoted to the discussion and conclusions.

2. Method

We based this work on the semi-analytic model (SAM) described in previous papers (see Menci et al. 2005, 2014, 2016). However, the specific version of the SAM adopted here differs from the one presented in the above papers in that it implements a new detailed description of AGN feedback, as we discuss in detail below (Sect. 2.1), and it is extended to describe the galaxy stellar velocity dispersion (Sect. 2.2).

The backbone of the model is based on the merging trees of the DM halos, which are generated through a Monte Carlo procedure, with merging probabilities given by the extended Press & Schechter formalism (see Bond et al. 1991; Lacey & Cole 1993), assuming a CDM power spectrum of perturbations in a concordance cosmology with density parameter Ω₀ = 0.3, a baryon density parameter Ω_b = 0.04, a dark energy density parameter Ω_Λ = 0.7, and a Hubble constant h = 0.7.

The subhalos included into the main halo can coalesce with the central galaxy after the orbital decay due to dynamical friction, or merge with other satellite subhalos, as described in our previous papers (Menci et al. 2005, 2014, 2016). Gas cools inside the halos due to radiative processes and settles into a rotationally supported disk with mass M_gas, with a scale length r_d and a circular velocity v_d related to the DM circular velocity V_c, as given in Mo et al. (1998). The stars are converted from the gas through three channels: (1) quiescent star formation with long timescales: ∼1 Gyr, (2) starbursts following galaxy interactions with timescales ≲100 Myr, according to BH feeding, (3) the loss of angular momentum triggered by the internal disk instabilities causing the gas inflows to the center, which stimulates star formation (as well as BH accretion). The stellar feedback is also considered by calculating the energy released by the supernovae associated with the total star formation, which returns a fraction of the disk gas into a hot phase. A Salpeter IMF is adopted in the SAM simulation.

We assumed a BH seed M_seed = 100 M_⊙ (Madau & Rees 2001) to be initially present in all galaxy progenitors at the initial redshift z = 15. This constitutes an approximate way of rendering the effect of the collapse of PopIII stars. However, the detailed value of M_seed has a negligible impact on the final BH masses as long as they remain in the range M_seed = 50 − 500 M_⊙.

The BH accretion is triggered by both interactions disk instabilities

(1) triggered by interactions. Each model galaxy (with a tidal radius r_t and a circular velocity V_c) in our Monte Carlo simulation may interact inside a host DM halo with a circular velocity V at a rate ${\tau }_r^{-1}=n_T\mathrm{\Sigma }(r_{\text{t}},V_{\text{c}},V)V_{\text{rel}}(V)$. Here V_rel is the average relative velocity of subhalos inside a host DM halo with a circular velocity V, and n_T is their number density in the common DM halo. The interaction rate determines the probability for encounters, either fly-by or merging, through the corresponding cross sections Σ given in Menci et al. (2014). The fraction f = (1/2)|Δj/j| of gas that is destabilized in each interaction corresponds to the relative loss Δj of orbital angular momentum j. In the above equation, the pre-factor accounts for the probability 1/2 of inflow rather than outflow related to the sign of Δj. Both the destabilized gas fraction f and the cross section Σ depend on DM the mass of the considered galaxy M_h and on the DM mass $M_h^{\prime }$ of the partner galaxy in the interaction that we extracted in our Monte Carlo procedure. Both major mergers (M_h ∼ $M_h^{\prime }$) and minor mergers (M_h ≫ $M_h^{\prime }$) are considered.

(2) Induced by disk instabilities. We assumed that disk instabilities arise in galaxies whose disk mass exceeds $M_{\text{crit}}=v_{\text{max}}^2{r}_d/Gϵ$ with ϵ = 0.75, where v_max is the maximum circular velocity associated with each halo (Efstathiou et al. 1982). Here, the parameter ϵ ≈ 0.5 − 0.75 is calibrated in the simulations. Although more refined treatments might be implemented (see, e.g. Romeo et al. 2023), this criterion is commonly adopted in semi-analytic models: Hirschmann et al. (2012) adopted a value ϵ = 0.75, and a similar value was adopted by Fanidakis et al. (2012). In order to investigate the maximum affect of DI on the statistical evolution of AGN, we also adopted the value ϵ = 0.75. This criterion strongly suppresses the probability for disk instabilities to occur not only in massive gas-poor galaxies, but also in dwarf galaxies that are characterized by low values of the gas-to-DM mass ratios. The instabilities induce loss of angular momentum, which results in strong inflows that we computed following the description in Hopkins & Quataert (2011), recast and extended as in Menci et al. (2014).

The BH accretion rates associated with the two feeding modes above yield an AGN emission with a bolometric luminosity L_AGN = η Ṁ_BH c², where we assumed the accretion efficiency to take the value η = 0.1 (Yu & Tremaine 2002; Marconi et al. 2004; Shankar et al. 2020a,b). The radiation field of the AGN pushes the cold gas disk outward, determining the expansion of a shock front (Silk & Rees 1998; Cavaliere et al. 2002; King 2003; Lapi et al. 2005; Granato et al. 2004; Silk & Nusser 2010; King et al. 2011; Zubovas & King 2012; Faucher-Giguere & Quataert 2012; King & Pounds 2015). This results in the expulsion of a fraction of disk gas (AGN feedback), which is also described in the model. However, we substantially upgraded this description, as shown in detail in the next subsection.

This model was upgraded in two fundamental respects: First, the AGN feedback is now described on the basis of a complete two-dimensional model that describes the expansion of the blast wave associated with the outflow in a disk geometry, and the outflow mass is computed from this model by integrating the mass outflow rate over the directions where the outflows escape from the galactic disk. Second, we implemented a description of the bulge size and of the velocity dispersion of stars in the bulge. We describe the two upgrades in turn.

2.1. Implementing a two-dimensional physical description of outflows in galactic disks

In a previous paper (Menci et al. 2019), we provided a compact two-dimensional description for the expansion of AGN-driven shocks in realistic galactic disks with exponential gas density profiles in a disk geometry. The gas distribution was assumed to be axisymmetric (with the X-axis aligned with the plane of the disk and the Y coordinate aligned with the rotation axis), with a density profile ρ(r) = ρ₀ exp(−r/r_d) depending only on the galactocentric distance r and on the scale length r_d, but with a cutoff in the Y direction at a distance corresponding to the disk scale height h, that is, ρ = 0 for Y ≥ h. The normalization ρ₀ was such as to obtain the disk gas mass M_gas when integrated out to large radii (i.e. r → ∞), while the scale height was assumed to be constant with radius for a given galaxy and to increase with the galaxy circular velocity according to the observed average relation h = 0.45 (V_c/100 km s⁻¹ −0.14 kpc (see van der Kruit & Freeman 2011 and references therein). Details are given in Menci et al. (2019).

The description accounts for the balance between the pressure term (fueling the expansion) that acts on the surface element corresponding to the solid angle in the considered direction and the counteracting gravitational term that is determined by the total mass M (contributed by the DM and by the central BH) within the shock radius. It includes the cooling of both the shell and gas inside the bubble due to the relevant atomic processes (inverse Compton and free-free emission; see Richings & Faucher-Giguere 2018).

We derived solutions for the outflow velocity V_S, θ, the mass outflow rate Ṁ_S,θ, and the shock position R_S, θ as a function of time in different directions defined by the angle θ with respect to the plane of the disk. These depend on three fundamental quantities related to host galaxy: the AGN bolometric luminosity L_AGN (fueling the expansion), the mass of cold gas in the disk M_gas (the mass that is pushed by the wind), and the total circular velocity V_c (determining the depth of the potential wells that counteract the expansion).

The solutions are plotted in Fig. 1 for a reference galaxy with a DM mass M_h = 10¹² M_⊙, a gas mass M_gas = 10¹⁰ M_⊙, and an AGN bolometric luminosity L = 10⁴⁵ erg s⁻¹. The X coordinate represents the distance from the center in the direction parallel to the plane of the disk, while the Y coordinate corresponds to the distance in the (vertical) direction perpendicular to the disk.

Fig. 1. Properties of shock expansion for our reference galaxy. In the top panel we show the positions of the shock at the different times represented by different colors and displayed in the right bar. The central and bottom panels show the velocity map and mass outflow rate map, respectively. The values corresponding to the colored contours are displayed in the bars. The X and Y coordinates correspond to the distance from the galaxy center in the directions parallel and perpendicular to the plane of the disk, respectively.

The shock expansion radius follows the paths of least resistance (see the left panel of Fig. 1), yielding an elongated shock front in the vertical direction. Inspection of Fig. 1 shows that while in the direction perpendicular to the disk, the outflows reaches a distance of 20 kpc in approximately 10⁷ yr, it takes about 10⁸ yrs to reach the same distance in the plane of the disk. This has important implications for studies of AGN feedback in galaxy formation models. or an AGN life time ∼10⁸ yrs, for instance, this would result in null gas expulsion along the plane of the disk.

Along the plane of the disk, the velocity V_S, θ rapidly decreases with increasing radius (central panel), while in the vertical direction, the shock decelerates until it reaches the disk vertical boundary h, but it rapidly accelerates afterward due to the drop in the gas density outside the disk. The opposite is true for the mass outflow rate (right panel), which instead grows appreciably only along the plane of the disk, where the higher densities allow values of Ṁ_{S,θ = 0} ∼ 10³ M_⊙ yr⁻¹.

These shock quantities depend on the global properties of the host galaxy and of the central AGN. For a given set of galaxy properties (cold gas mass M_gas in the disk, circular velocity of the DM halo V_c, and AGN luminosity L_AGN), the predictions of the model (mass outflow rates and outflow velocities in any direction) can be compared with observation. Menci et al. (2019) tested the model against a compilation that was new then (Fiore et al. 2017) of observed outflows in 19 galaxies with different measured gas and dynamical mass, allowing for a detailed one-by-one comparison with the model predictions. Specifically, for each observed object, the measured AGN and host galaxy properties (AGN luminosity, gas mass, and circular velocity) were used to obtain the input quantities for the model. The resulting velocity and mass outflow rates computed from the model were then compared with the observed corresponding wind properties. For each considered galaxy and at the observed outflow radii, the model yielded values of the velocity and mass outflow rates (averaged over the directions) that agreed well with the observations. The predicted densities of the shocked shell were consistent with the observed molecular emission of the outflows in the vast majority of cases. The agreement we obtained for a wide range of host galaxy gas masses (10⁹ M_⊙ ≲ M_gas ≲ 10¹² M_⊙) and AGN bolometric luminosities (10⁴³ erg s⁻¹ ≲ L_AGN ≲ 10⁴⁷erg s⁻¹) provides a quantitative systematic test for the modeling of AGN-driven outflows in galactic disks. This provides a solid basis for a reliable implementation of AGN feedback in galaxy formation models.

Implementing this modeling of AGN outflows in the SAM constitutes a challenging task. Directly solving the equation for the expansion of outflows in each of the simulated galaxies generated by our SAM would be too demanding in terms of computation time. Thus, we proceeded as follows.

We considered a grid of values for the three input quantities of the outflow model. Specifically, we considered 20 equally spaced logarithmic values for the AGN luminosity in the range log (L_AGN/erg s⁻¹) = 42 − 47 for a gas mass in the range log (M_gas/M_⊙) = 8 − 11.5 and for a circular velocity of the host galaxy log (V_c/100 km/s) = 0.5 − 4.

For each combination (L_AGN, M_gas, V_c), we ran our model for the expansion of the outflows and computed R_S, θ, V_S, θ, and the mass outflow rate Ṁ_S,θ as a function of time in different directions defined by the angle θ. We then computed the total outflow mass M_S by integrating Ṁ_S,θ over time and over all the directions in which the outflow escapes, that is, the directions in which the velocity V_S, θ accelerates after it reaches the disk boundary (see Fig. 1) to reach velocities well beyond the escape velocity from the galactic halo. This provided us with tabulated values of M_S (and of the critical opening angle θ_c corresponding to escaping outflows) for each combination (L_AGN, M_gas, V_c). The expelled gas fraction f_expell = M_S/M_gas and the opening angle θ_c as a function of L_AGN and M_gas are shown in Fig. 2 for a reference value of the circular velocity V_c = 200 km s⁻¹.

Fig. 2. Properties of the gas outflows. The color code shows the dependence of the fraction of cold gas that is expelled by AGN outflows (fexpell, top panel) and of the solid angle θc swept out by AGN outflows (bottom panel) on the AGN bolometric luminosity LAGN and of the cold gas mass Mcold. A DM circular velocity Vc = 200 km s−1 is assumed for both panels.

To compute the effect of AGN feedback on each galaxy of our SAM, we interpolated these tabulated values of M_S and θ_c to determine the escaped gas mass and opening angle corresponding to the combination (L_AGN, M_gas, V_c) associated with each model galaxy in the SAM. The distribution of the average values of the expelled gas fraction f_expell as a function of the stellar mass M_* and star formation rate Ṁ_∗ is shown in Fig. 3.

Fig. 3. Average fraction of gas expelled by AGN outflows in galaxies with different stellar mass M* and star formation rate Ṁ∗ shown by the color code.

The distribution shows a bimodal shape at low values of M_*. HIgh values of f_expell require either low values of M_gas (and hence low values of Ṁ_∗), and/or high values of L_AGN (attained for a high gas mass M_gas, yielding high values of Ṁ_∗ on average). For high values of M_*, the outflows can overcome the effect of the deep galaxy potential wells (yielding large expelled fractions f_expell ≈ 1) only when powered by a large AGN energy injection L_AGN, which can only be attained for high values of M_gas and hence of Ṁ_∗. The inclusion of a kinetically driven AGN feedback can thus trigger powerful outflows in massive galaxies that can push large portions of the gas in the intergalactic medium. This expelled gas is not reaccreted at later times. We verified that this feedback model yields stellar galaxy luminosity functions that are consistent with observations, although in this respect, the isotropic feedback model previously adopted in our SAM performed equally well.

We verified that the new treatment of AGN feedback described above yields AGN luminosity functions that are consistent with existing observations. This is shown in Fig. 4, where the predicted AGN bolometric luminosity function is compared with the compilation by Shen et al. (2020) in different redshift bins extending up to z = 6. The slight overprediction at z ≳ 3 is expected. The predicted luminosity functions include Compton-thick AGN, which constitute a relevant fraction of the AGN population, with an appreciable uncertainty range ∼20 − 50% at z ≳ 3 (Shen et al. 2020), which is consistent with the slight overestimate of the observed abundance in the same redshift range. In addition, large observational uncertainties still affect the measurement of the faint end of the AGN luminosity function at z ≳ 3 (see, e.g., Shankar & Mathur 2007; Ricci et al. 2017; Shen et al. 2020). Based on results by recent surveys, the number of quasars at redshifts z ≳ 3 is being constantly revised upward. Several authors have obtained abundances of low-luminosity AGN that are up to 40 % higher than the estimates reported in the analysis by Shen et al. (2020) at both faint (M_UV ≈ −23.5, Boutsia et al. 2018; Giallongo et al. 2019) and bright (M_UV ≈ −27, Boutsia et al. 2021; Grazian et al. 2020) magnitudes.

Fig. 4. Predicted AGN bolometric luminosity functions (solid lines) are compared to the compilation of data in Shen et al. (2020, shaded regions) in six redshifts bins shown in the labels.

After we tested the new treatment of AGN feedback in the SAM, the AGN luminosity functions were consistent with existing observations over a wide range of redshift. We then proceeded with the extension of the SAM to include a description of the stellar velocity dispersion, as detailed below.

2.2. Computing Jeans-based stellar velocity dispersions

We here further extended our SAM with respect to previous renditions of the model to include the most recent calculations of the half-mass radii and stellar velocity dispersions. The former were presented in Zanisi et al. (2020). Half-mass galactic radii in progenitor disks are computed from angular momentum conservation with the host halo spin. When galaxies merge, a stellar bulge is formed with the half-mass radius R_h, which is computed via energy conservation between the gravitational potential energies of the progenitors and the descendant (e.g., Cole et al. 2000). We tested that the sizes computed in our SAM match the observed correlation with stellar mass for the redshift range covered by observations for star-forming and quiescent galaxies. To perform the comparison, we adopted the approach in Shankar et al. (2014). Assuming that light traces mass, the above authors converted R_h into projected half-light radii R_e using the tabulated factors from Prugniel & Simien (1997), that is, R_e ≈ 2S(n) R_h, with the scaling factors S(n) dependent on the Sérsic index n. Following Shankar et al. (2014), for we adopted a fixed value n = 4 for this. As noted by the above authors, while it is possible to predict a Sérsic index a priori from the models in principle (e.g., Hopkins et al. 2009), this relies on several additional assumptions on the exact profile and the evolution with time of the dissipational and dissipationless components, so that the true advantage with respect to simply empirically assigning a constant Sérsic index is modest. By the same token, we note that none of our main results would change if we were to assign stellar mass-dependent Sérsic indices and/or effective radii directly extracted from the observed distributions to our mock galaxies.

The predicted correlation between the size R_e and the stellar mass is compared with observations in Fig. 5 for different redshift bins and for both quiescent and star-forming galaxies. Following the approach in Genel et al. (2018), we classified our SAM-generated galaxies as early- or late-type through their specific star formation SSFR = Ṁ_∗/M_∗. Galaxies are defined as quenched when their SSFR is at least 1 dex below the ridge of the star formation main sequence. This is defined, for simplicity, as log(SSFR/Gyr⁻¹) =  − 0.94, −0.85, −0.35, and 0.05 for z = 0, 0.1, 1, and 2, respectively. On the observational side, van der Wel et al. (2014) classified early- and late-type galaxies on the basis of their star formation properties, defined in terms of cuts in the (U–V)–(V–J) rest-frame colors, while Bernardi et al. (2014) adopted a morphological classification based on a Bayesian automated procedure. We assumed that our model half-light radius R_e was computed assuming that light traces mass. On the other hand, the different observations are made in different bands. While the measurements of Bernardi et al. (2014) were performed directly in the Sloan Digital Sky Survey (SDSS) r band, those from van der Wel et al. (2014) were made in the Hubble Space Telescope Wide Field Camera 3 (WFC3) F125W filter, which is centered on an observed frame wavelength of 1.25 μm. In order to match the SDSS-based studies, we converted the measurements in van der Wel et al. (2014) into the SDSS r band using the gradient d log R_e/dlogλ that they directly measured between the F814W (814 nm) and F125W filters.

Fig. 5. Predicted distribution of the effective radius Re shown as a function of the stellar mass M* at redshifts z = 0.25, z = 1.51, and z = 2.25 (from left to right) for early-type (top row) and late-type galaxies (bottom row). The predictions are compared with the data from van der Wel et al. (2014; shaded areas). We also show the size-stellar mass relation measured at z = 0.1 by Bernardi et al. (2014) and Mowla et al. (2019) as dashed lines and dots, respectively.

From the galaxy size computed as shown above, we derived the corresponding velocity dispersion by solving the spherical Jeans equation, following Desmond & Wechsler (2017). For the bulge component, this reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}\rho (r){\sigma }_r^2}{\mathrm{d}r}+\frac{2\beta (r)}{r}\rho (r){\sigma }_r^2=-G\rho (r)\frac{M_{\ast ,b}(r)}{r^2},\end{array}$$(1)

where $\beta (r)=1-{\sigma }_{\theta }^2/{\sigma }_r^2$ is the orbit anisotropy profile relating the radial velocity dispersion σ_r to the tangential component σ_θ, ρ(r) is the mass density profile of the bulge, and M_*, b is the bulge stellar mass.

The line-of-sight velocity dispersion solving Eq. (1) is then given by Mamon & Łokas (2005; see also Desmond & Wechsler 2017) in the form

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2(R)=\frac{2G}{\mathrm{\Sigma }(R)}{\displaystyle {\int }_R^{\infty }}K(r/R)\rho (r),M_{\ast ,b}(r)\frac{\mathrm{d}r}{r},\end{array}$$(2)

where R is the projected distance from the center, and Σ(R) is the projected density profile. The function K is defined as a combination of gamma functions Γ and incomplete β functions B,

$$\begin{array}{c}\hfill K(u)\equiv 1/2u^{2\beta -1}[(3/2-\beta )\sqrt{\pi }(\mathrm{\Gamma }(\beta -1/2)/\mathrm{\Gamma }(\beta ))\end{array}$$(3)

$$\begin{array}{cc}& +\beta B(1/u^2,\beta +1/2,1/2)-B(1/u^2,\beta -1/2,1/2)],\hfill \end{array}$$(4)

where u ≡ r/R, and β is the assumed orbit anisotropy.

For each model galaxy, we assumed for Σ(R) a Sersic form Σ(R) = Σ₀ g(R) with

$$\begin{array}{c}\hfill g(R/R_{\mathrm{e}})=e^{-b_n[(\frac{R}{R_{\mathrm{e}}}{)}^{1/n}-1]},\end{array}$$(5)

where b_n ≈ 2n − 1/3 + 0.009876/n. This was normalized to yield the bulge mass M_*, b when integrated over R, so that ${\mathrm{\Sigma }}_0=(1/2\pi Q)M_{\ast ,b}/R_{\mathrm{e}}^2$, where $Q\equiv {\int }_0^{\infty }g(\mathit{y})\mathit{y}\mathrm{d}\mathit{y}$.

The corresponding deprojected mass density profile is ρ(r) = ρ₀ f(r/R_e), where

$$\begin{array}{c}\hfill f(r/R_{\mathrm{e}})=(\frac{r}{R_{\mathrm{e}}}{)}^{-p_n}{e}^{-b_n(\frac{r}{R_{\mathrm{e}}}{)}^{1/n}},\end{array}$$(6)

where p_n = 1 − 0.6097/n + 0.00563/n². Again, the normalization ρ₀ was computed so as to obtain the total mass M_*, b when the above mass density profile was integrated over the volume. Thus, ${\rho }_0=(1/4\pi F)M_{\ast ,b}/R_{\mathrm{e}}^3$, where $F\equiv {\int }_0^{\infty }g(\mathit{y}){\mathit{y}}^2\mathrm{d}\mathit{y}$.

Finally, the bulge mass density profile is $M_{\ast ,b}(r)=4\pi {\int }_0^r\rho (\mathit{y}){\mathit{y}}^2\mathrm{d}\mathit{y}$. Inserting the above form for the density profile ρ(r), we obtain

$$\begin{array}{c}\hfill M_{\ast ,b}(r/R_{\mathrm{e}})=(M_{\ast ,b}/F){\displaystyle {\int }_0^{r/R_{\mathrm{e}}}}{y}^{2}f(y)\mathrm{d}y\equiv (M_{\ast ,b}/F)h(r/R_{\mathrm{e}}).\end{array}$$(7)

Inserting Eqs. (4)–(6) in Eq. (2), we obtain

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2(x)=\frac{G{M}_{\ast ,b}}{R_{\mathrm{e}}}\frac{1}{FQ}{\displaystyle {\int }_x^{\infty }}K(y/x)f(y)h(y)\mathrm{d}y/y,\end{array}$$(8)

where x ≡ R/R_e.

In the following, we consider the average velocity dispersion

$$\begin{array}{c}\hfill {\sigma }^2=\frac{{\int }_0^{\infty }\mathrm{\Sigma }(x){\sigma }_{\mathrm{los}}^{2}x\mathrm{d}x}{{\int }_0^{\infty }\mathrm{\Sigma }(x)x\mathrm{d}x}.\end{array}$$(9)

When a finite-aperture radius R_ap is considered, the upper limits in the above integrals should be replaced with x_ap ≡ R_ap/R_e. Inserting the expression in Eq. (7), we obtain

$$\begin{array}{c}\hfill {\sigma }^2=\frac{G{M}_{\ast ,b}}{R_{\mathrm{e}}}W,\end{array}$$(10)

with the shape factor W defined as

$$\begin{array}{c}\hfill W(\beta ,n)\equiv \frac{1}{FQ}\frac{{\int }_0^{\infty }g(x)x\mathrm{d}x{\int }_x^{\infty }K(y/x)f(y)h(y)\mathrm{d}y/y}{{\int }_0^{\infty }g(x)x\mathrm{d}x}.\end{array}$$(11)

The velocity dispersion is thus related to the virial value GM_*, b/R_e through the shape factor W(β, n), which depends on the structural properties of the bulge (anisotropy parameter β and Sérsic index n).

For each simulated galaxy in our SAM, the effective bulge mass content M_*, b and radius R_e allowed us to derive the average bulge velocity dispersion from Eqs. (10) and (11) for a given choice of β and n. For the first, we extracted a random value from a Gaussian distribution with an average value ⟨β⟩ = 0.3 and dispersion 0.35. This choice has been shown to provide a good fit to observational data (Desmond & Wechsler 2017). For the Sérsic index n, we adopted a fixed vale n = 4 based on the considerations presented above. To test that the treatment above was consistent with available data, we compared (Fig. 6) the velocity dispersion of passive model galaxies (Ṁ_∗/M_∗ ≤ 0.1 Gyr⁻¹) at low redhsift (z ≤ 0.1) to existing data for early-type galaxies in the SDSS database (Bernardi et al. 2011), selected on the basis of morphological indicators (Hyde & Bernardi 2009). The observational correlation does not change appreciably when passive galaxies are selected on the basis of spectral or color indicators (Belli et al. 2014; Zahid et al. 2016). These galaxies are characterized by high bulge-to-disk stellar mass ratios B/T ≳ 0.8, for which the total galaxy velocity dispersion is dominated by the bulge component computed in Eq. (8)–(11). For consistency with the data by Hyde & Bernardi (2009), we adopted x_ap = 1/8 when computing Eq. (8), (this leads to a small offset −0.03 in log(σ) when compared to the global average corresponding to x_ap → ∞). The excellent agreement shown in Fig. 6 provides a reliable basis for computing the scaling relations between the BHs and the properties of the host galaxies, as well as for including the effects of observational selection biases in the analysis of model results, as we discuss below.

Fig. 6. Predicted relation between the velocity dispersion σ and the stellar mass M* for passive galaxies (see text). The color code represents the logarithm of the fraction of model galaxies with the considered σ for each bin of M*. The shaded band shows the relation for early-type galaxies selected on the basis of morphological indicators after Hyde & Bernardi (2009). The dashed black and red lines show the average relation derived by Zahid et al. (2016) for passive SDSS galaxies and for the Smithsonian Hectospec Lensing Survey (SHELS) sample, selected on the basis of spectral properties (Dn4000 > 1.5, where Dn4000 is the flux ratio of two spectral windows adjacent to the 4000 Å break (see Balogh et al. 1999).

3. Selection bias

Shankar et al. (2016) showed that a selection bias may affect the observe scaling relations when dynamical measurements of the BH mass are considered. These mass estimates are only possible when the sphere of influence of the BH,

$$\begin{array}{c}\hfill r_{\mathrm{infl}}\equiv G{M}_{\mathit{bh}}/{\sigma }^2,\end{array}$$(12)

has been resolved in the observations (e.g., Peebles 1972; Gültekin et al. 2011; Graham & Scott 2013).

This means that the observed correlations between dynamically measured BH masses and the galaxy properties automatically select objects for which

$$\begin{array}{c}\hfill {\theta }_{\mathrm{infl}}\equiv r_{\mathrm{infl}}/d_{\mathrm{ang}}\ge {\theta }_{\mathrm{crit}},\end{array}$$(13)

where d_ang is the angular distance of the object, and θ_crit is the resolution limit of the instrument (e.g., θ_crit ≈ 0.1 for observations with the Hubble Space Telescope, HST).

Shankar et al. (2016) estimated the impact of this observational bias on the comparison between the predicted and the observed scaling relations under different assumptions on the intrinsic scaling relations.

Our improved semi-analytic model allowed us to investigate the impact of this observational bias on the basis of the ab initio model for galaxy formation presented in the previous section, tested against a wide number of observations. Our computation of the velocity dispersion σ in terms of the properties of the host galaxy enables the prediction of the BH sphere of influence r_inf for each simulated galaxy of our model.

In the following, we exploit the model predictions of r_inf to investigate the impact of the selection bias affecting the dynamical BH mass measurements on the comparison between the predicted and the observed BH scaling relations.

4. Results

Before performing a detailed comparison of our predictions for the BH scaling relations with available data, we first show the effects of the model improvements presented in the previous section, that is, the inclusion of the two-dimensional angle-dependent AGN feedback model presented in Sect. 2.1 and the σ dependent selection bias described in Sect. 3. With this aim, we compare in Fig. 7 the M_BH − M_*, b relation that we obtained by adopting our two-dimensional feedback model (upper panels) with that obtained using the isotropic AGN feedback model previously adopted in our SAM (lower panels). To show the effect of the selection bias, the results from models including the bias (right columns) are compared to those without this effect (left columns).

Fig. 7. Effect of feedback treatment on the MBH − M*, b relation. The relation predicted for galaxies at z ≤ 0.1 by our new two-dimensional model for AGN feedback (top panels) presented in Sect. 2.1 is compared with the predictions obtained with our previous isotropic treatment of feedback (bottom panels). For both AGN feedback models, we show the results with and without the selection effect discussed in Sect. 2.1, as indicated by the labels. In all panels, the color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M*, b bin. The standard deviations ϵ for the new and previous treatment of AGN feedback is shown as a function of bulge stellar mass in the rightmost histograms. The horizontal strip represents the range of values obtained from different observational works: the upper bound is the intrinsic scatter ϵ derived by De Nicola et al. (2019), and the lower bound is the value of ϵ derived by the least-square fit analysis by Kormendy & Ho (2013) using individual errors in log σ, adding individual errors in log MBH to the intrinsic scatter in quadrature, and iterating the intrinsic scatter until the reduced χ2 = 1.

The comparison between the first and the second column in Fig. 7 shows that a true relation between the BH mass M_BH and the bulge mass is indeed expected in the isotropic and in the two-dimensional, angle-dependent AGN feedback models, regardless of whether the selection effect resolution of the BH sphere of influence is considered. The main effect of including the selection bias is that the distribution of the model galaxies is extended toward higher BH masses (especially for high M_*, b), which agrees better with the observations.

It is interesting to note that the slope of the scaling relations in the case of two-dimensional, angle-dependent AGN feedback is almost unchanged for an isotropic AGN feedback. We therefore conclude that the detailed form of AGN feedback is not the main cause for the scaling relations, which instead appear mostly to be driven in our model by the dependence of the BH mass on the processes connected to the growth of the stellar content and of the host galaxy, which in turn imply a correlation between M_BH and σ. Interestingly, Marsden et al. (2022) have recently shown that a steep and nearly constant M_BH − σ relation naturally arises from the coupling of the observed L_X − M_* relation and a Jeans-based stellar velocity dispersion (similarly to what we carried out here), without any additional fine-tuning.

However, the adoption of the two-dimensional, angle-dependent description of AGN feedback has an important effect on the scatter of the relations, as shown quantitatively in the rightmost panels of Fig. 7. The adoption of our new two-dimensional model for AGN feedback reduces the intrinsic scatter from values ϵ ≲ 0.5 − 0.6 (obtained in the isotropic feedback case) to values ϵ ≲ 0.3. This is due to the following reasons: In the new two-dimensional model for feedback, the blast wave expansions stalls along the direction of the disk, and the radius within which the expansion stops strongly depends on the gas density of the disk and the AGN luminosity. This means that the opening angle (and hence the fraction of expelled gas) is larger when the gas density is low (because less energy has to be spent to push the gas outward) and when the AGN luminosity is high (because more energy is available to push the blast wave outward), as shown in Fig. 2. Both quantities depend on the merging histories and are related because the AGN luminosity L_AGN depends on the available cold gas reservoir M_gas. The high efficiency of feedback in galaxies with particularly low M_gas (for a given L_AGN) or in those with particularly large L_AGN (for a given M_gas) inhibits the BH growth in all the host galaxies that are outliers with respect to the average relation between M_gas and L_AGN. This results in a smaller scatter.

The M_BH − M_*, b relation predicted by our improved model is compared with the data in Fig. 8. The observational data with which we compared the relation include BH estimates for both quiescent and active BHs.

Fig. 8. Comparison of our predicted MBH − M*, b relation with the data for all galaxies at z ≤ 0.1. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M*, b bin. Panel A: Predicted relation between the BH mass MBH and the bulge stellar mass M*, b. The data points are from Kormendy & Ho (2013; black squares and triangles for elliptical galaxies and bulges, respectively), Savorgnan & Graham (2016; red dots), and Ho & Kim (2014; orange filled squares). Panel B: Same as panel A, but including the selection bias discussed in Sect. 3, i.e., considering only model galaxies satisfying the condition in Eq. (12) with θcrit = 0.1. Panel C: Same as panel A, but considering only model galaxies hosting AGN with bolometric luminosity LAGN ≥ 1044 erg s−1. Bottom right panel: Rms value of the scatter of the distributions shown in panel A (brown histogram), panel B (blue histogram), and panel C (orange histogram) for different bins of M*, b. The horizontal strip corresponds to the range of measured intrinsic scatter by different authors. The lower value ϵ = 0.28 was obtained by Kormendy & Ho (2013) and the upper value ϵ = 0.4 was obtained by Savorgnan et al. (2016).

First, we considered the BH mass measurements collected from spatially resolved estimates available from the literature for bulges or bulge-dominated galaxies, that is, the data set including bulges and elliptical galaxies in Kormendy & Ho (2013), and the extended data sets in Savorgnan & Graham (2016) and De Nicola et al. (2019). The first data set consists of 66 galaxies with dynamical estimates of their black hole masses as reported by Graham & Scott (2013) or Rusli et al. (2013). Using 3.6 μm (Spitzer satellite) images, Savorgnan & Graham (2016) modeled the one-dimensional surface brightness profile of each of the 66 considered galaxies and estimated the structural parameters of their spheroidal component by simultaneously fitting a Sérsic function. Following Shankar et al. (2016), we retained only the galaxies with secure black hole mass measurements from the original Savorgnan & Graham (2016) sample and removed the sources that wer eclassified as ongoing mergers. This limited the final sample to 48 galaxies, 37 of which are early-type galaxies (ellipticals or lenticulars) that we considered in our comparison. De Nicola et al. (2019) derived a sample of 83 galaxies based on the compilation of Saglia et al. (2016) whose BH masses were derived from spatially resolved kinematics.

We also compared our predicted M_BH − M_*, b relation with the observed relation for active galaxies whose BH masses were derived from the (presumed) virial motions of the broad line region (BLR) gas cloud orbiting in the vicinity of the central compact object: M_BH = f_vir r(ΔV)²/G. Here, r is the radius of the BLR, which is derived from reverberation mapping (e.g., Blandford & McKee 1982; Peterson 1993), or reverberation-based methods that use the radius-luminosity relation (e.g., Bentz et al. 2006). The characteristic velocity ΔV was derived from the width of the emission lines (a common line is Hβ). As motions in the BLR are not perfectly Keplerian, the parameter f_vir is included in the equation to account for the uncertainties in kinematics, geometry, and inclination of the clouds. This was calibrated by comparing reverberation-mapped AGNs with measured bulge stellar velocity dispersions against the M_BH − σ relation of inactive galaxies. We considered the compilation by Ho & Kim (2014), who derived a virial coefficient f_vir = 6.3  ±  1.5 for classical bulges and ellipticals.

The treatment of the AGN feedback we adopted in our semi-analytic model yields BH scaling relations that match the observations when the observational bias related to the resolution of the BH sphere of influence is considered. Again, the inclusion of this bias extends the distribution of model galaxies toward higher BH masses (especially for high M_*, b), which agrees better with the observations. Shankar et al. (2016) discussed that the degree of the impact of the bias between the normalization of the intrinsic and observed BH-galaxy scaling relations depends on the type of the intrinsic relation that is considered. Assuming an intrinsic relation scaling M_BH ∝ σ^α with α ≳ 4 would lead to a jump between intrinsic and observed relations that is larger than assuming, for example, $M_{\text{BH}}\propto M_{*}^{\beta }$ with β ≈ 1 as the intrinsic relation. Our current SAM predicts an average difference of ≲2 between intrinsic and observed scaling relations, which is what would be expected if the BH mass depends more weakly on the velocity dispersion and/or is more strongly correlated with stellar mass than with the velocity dispersion (models 2 and 3 in Shankar et al. 2016), as also highlighted in the residuals below.

We also note that the predicted distribution of model galaxies with an active AGN (L_AGN ≥ 10⁴⁴ erg s⁻¹) does not appear to differ substantially from the distribution of the whole population. However, some of the observational points related to active galaxies lie out of the strip representing the model predictions. While a possible explanation is that some of the outflows in these galaxies are spherical, it is possible that the apparent discrepancy in the normalizations of the M_BH − M_*, b scaling relations characterizing local AGN data and local inactive BH samples with dynamical mass measurement is due to the presence of some additional biases in the BH and/or the stellar mass estimates between the two BH samples, so that the AGN data represent a biased sample of the whole population (see, e.g., Shankar et al. 2019; Farrah et al. 2023, and references therein). We further investigate this point below.

The scatter we obtain agrees well with the observations, for which most observations yield values in the range 0.3 ≲ ϵ ≲ 0.4 (see the references in the caption), although some authors find somewhat higher values (see Sahu et al. 2019a,b). Tehe novel treatment of AGN feedback we considerd reduces the scatter of these relations with respect to the adoption of an isotropic average feedback. This is shown in Fig. 9, where we compare the M_BH − M_*, b relation (including the selection bias described in Sect. 3) obtained with the two-dimensional model presented in Sect. 2.1 with our previous isotropic treatment of AGN feedback described in Menci et al. (2014).

Fig. 9. Comparison of our predicted MBH − σ relation with data for all galaxies at z ≤ 0.1. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log σ bin. Panel A: Predicted relation between the BH mass MBH and the bulge velocity dispersion σ brightness-weighted averaged over the galaxy. The data points are from Kormendy & Ho (2013; black squares and triangles for elliptical galaxies and bulges, respectively), De Nicola et al. (2019; brown squares), Savorgnan & Graham (2016; red dots), and Ho & Kim (2014; orange filled squares). Panel B: Same as panel A, but including in the model the selection bias discussed in Sect. 3, i.e., considering only model galaxies that satisfy the condition in Eq. (12) with θcrit = 0.1. Panel C: Same as panel A, but considering only model galaxies hosting AGN with bolometric luminosity LAGN ≥ 1044 erg/s. Bottom right panel: Histograms represent the rms values of the scatter for the distributions shown in the top left (brown), top right (blue), and bottom left panel (orange). The horizontal strip represents the range of values obtained from different observational works: the upper bound is the intrinsic scatter ϵ derived by De Nicola et al. (2019), and the lower bound is the value of ϵ derived by the least-square fit analysis by Kormendy & Ho (2013) using individual errors in log σ, adding individual errors in log MBH to the intrinsic scatter in quadrature, and iterating the intrinsic scatter until the reduced χ2 = 1.

We present in Fig. 9 the predicted relation between BH mass M_BH and the bulge velocity dispersion σ for galaxies at z ≤ 0.1. In addition to the intrinsic relation, we also show the relations that we predict when the selection bias presented in the previous section is considered (i.e., for model galaxies that satisfy condition in Eq. (12)) and when only active galaxies (i.e., galaxies hosting an AGN with luminosity L_bol ≥ 10⁴⁴ erg s⁻¹) are considered. Although in the plot, σ refers to the bulge velocity dispersion, we verified that for galaxies with B/T ≥ 0.3, the model predictions do not change substantially when the global σ is computed from the total galaxy size and total stellar mass M_* (in this case, we used a stellar mass dependent Sérsic index as given in Terzic & Graham 2005).

In the comparison, we assumed that stellar mass traces light because the mass-weighted averaged velocity dispersion in Eq. (8) derived for model galaxies is compared with observed values derived from a brightness-weighted average (see, e.g., De Nicola et al. 2019).

The comparison between the observational data and the model predictions provides interesting insight on some crucial points.

– The comparison between panels A and B shows that a true relation between the BH mass M_BH and the bulge mass is indeed expected, regardless of the inclusion of selection effects. In the model, the observed correlation does not consist of the upper envelope of the scatter of a distribution of BH masses that extends to lower BH masses, as suggested by some authors (see, e.g., Batcheldor 2010), which is also in line with the findings by Shankar et al. (2016).

– However, when we compare the results in panels A and B, the distribution of the model galaxies in the M_BH − σ plane agrees better with the observational data when the selection effect discussed in Sect. 3 is considered. When this bias is included, the distribution of the model galaxies extends toward higher BH masses, which agrees better with observations. However, this extension is less pronounced than the extension that affects the predicted M_BH − M_*, b relation. This effect has been noticed by Shankar et al. (2016) and led these authors to conclude that this suggests that the M_BH − σ relation is the true fundamental relation, while the M_BH − M_*, b relation is mostly a consequence of the M_BH − σ and σ − M_*, b relation. We investigate this point in detail below.

– When the selection bias is considered, not only the bulk of the population of model galaxies is consistent the observed relation, but the scatter in the distribution also agrees excellently with the distribution of data points.

– It is extremely interesting to note that the SAM predicts a M_BH − σ relation of similar normalization and dispersion with respect to the collection of available data sets in the local Universe, but with a flatter slope. The predicted M_BH − σ relation has a slope of ∼2.5 when selection effects are included (panel B), and somewhat lower x (∼2.2) in the intrinsic relation. This result is in line with what has been suggested in other theoretical and numerical works. Cavaliere & Vittorini (2002; see also Vittorini et al. 2005) highlighted the fact that if the quasar-feedback condition (e.g., Silk & Rees 1998; King 2003) is not constantly applied throughout the evolution of BHs, the resulting M_BH − σ relation would have a slope closer to ∼3 rather than 4 − 5. This also agrees with some recent (Caglar et al. 2023) and previous (Woo et al. 2013, Park et al. 2012a) observational findings. On the numerical side, Li et al. (2020) recently showed that the TNG100 simulation also produces a M_BH − σ relation with a slope of ∼3 (their Fig. 1, bottom panels) that is in general above the data, where the latter effect could also at least in part be ascribed to selection effects in the local observational data sets (e.g., Shankar et al. 2016, Li et al. 2020). Sijacki et al. (2015) and Thomas et al. (2019) found in the Illustris and SIMBA simulations, respectively, a slope and normalization for the M_BH − σ relation that is consistent with observations. However, a M_BH ∝ σ^α, with α ≳ 4 − 5, in some simulations/models might also in part be an artifact of a steeper M_* − σ relation in the same models, as discussed in Barausse et al. (2017), as in the case of the Horizon-AGN simulation (Dubois et al. 2016).

– Analogously to what we predict for the M_BH − M_*, b relation, the predicted distribution of model galaxies with an active AGN (L_AGN ≥ 10⁴⁴ erg s⁻¹) does not seem to differ substantially from the distribution of the whole population.

In summary, the predicted distribution of galaxies in the M_BH − σ plane covers the regions allowed by available data. However, our SAM, in line with some other current models discussed above, does not fully match the steep slope of the M_BH − σ relation, as inferred from current local data. This might be interpreted in two ways: (1) either the local sample is incomplete, and/or (2) current cosmological models for the coevolution of BHs and galaxies are still missing some ingredients, either in the AGN feedback process itself and/or in other events during the growth of the BH (e.g., Cavaliere et al. 2002). In our model, the only viable way to increase the slope in the M_BH − σ relation would be to slightly increase the anisotropy parameter with stellar mass from 0.3 to 0.35, but this would also steepen the σ − M_* relation, which would worsen the match with observations. A trade-off between the inclusion of the dark matter component and stellar anisotropy might improve the fit to both observables simultaneously. The main lesson we learn from previous studies and the one conducted in this work is that at face value, AGN kinetic feedback by itself may not be sufficient to entirely reproduce the current observational data on local BHs. This is further corroborated by the analysis of the residuals that we present below.

To go beyond the study of pairwise correlations between M_BH, M_*, b and σ and gain better insight into the joint distribution of these quantities, we now focus on our predictions for the residuals of the BH-galaxy scaling relations. For any pair of quantities X and Y, these are defined as Δ(Y|X)≡log Y − ⟨log Y|log X⟩. As proposed by various authors (see, e.g., Shankar et al. 2016), the correlations between the residuals provide additional information about the joint distribution of M_BH, σ and M_*, and can provide a way of determining which variable is more important in determining the M_BH − M_* and M_BH − σ relations. In the ideal case where M_BH was fundamentally determined by σ alone, for instance, the residuals from correlations with σ would be uncorrelated, while residuals from the M_BH − M_* relation would correlate with residuals from the σ − M_* relation.

In the left and the right panels of Fig. 10 we show Δ(M_BH|M_*) vs Δ(σ|M_*) and Δ(M_BH|σ) versus Δ(M_*|σ) predicted for bulge-dominated model galaxies (B/T ratio ≥0.8), and compare these relations with those derived for the sample of Savorgnan et al. (2016) of early-type galaxies, for which residuals were computed in Shankar et al. (2016) and Barausse et al. (2017).

Fig. 10. Predictions for the residuals of the BH-galaxy scaling relations. In the left panel, we show the correlation between the residuals of the MBH − M*, b relations and those of the σ − M* relation at fixed stellar mass. The correlation between the residuals of the MBH − σ relation and those of the M* − σ relation at fixed velocity dispersion are shown in the right panel. Data points refer to the sample of Savorgnan et al. (2016). The contours refer to model galaxies with a bulge-to-total stellar mass ratio B/T > 0.8 for which the black hole sphere of influence can be resolved.

The predicted residuals occupy the same regions as the observational data points, showing that the model captures the complex intertwining between the considered quantities. As expected, the predicted behavior of the residuals is complex, and none of the considered quantities can be considered as really fundamental in our SAM. In the framework of cosmological galaxy evolution models such as our SAM, the BH mass stems from the complex interplay between the growth history of the halo mass and the physics of baryons, which cannot be easily reduced to a single fundamental dependence on a single quantity.

Testing the consistency of this picture with existing data is not straightforward because the model does not provide definite correlations between the residuals, but rather distributions in the Δ(M_BH|M_*)−Δ(σ|M_*) and Δ(M_BH|σ)−Δ(M_*|σ|) planes. Thus, we performed a two-dimensional Komlogorov-Smirnov (KS) test to assess the consistency of the predicted distributions with the observational points. To do this, we followed the approach in Fasan & Franceschini (1987) using the extended two-dimensional KS distance D defined in Peacock (1983). We ran a Monte Carlo simulation to generate synthetic data sets from the predicted distributions, each set with the same number of points N_data as the real data set. For each Monte Carlo realization, we also generated N_data observational data points distributed according to the measured data points and error bars. We then computed D for each synthetic data set using the approach mentioned above and counted the fraction of the time for which these synthetic Ds exceeded the D from the real data. The resulting fraction provides the significance level Q_KS for the consistency between the data distribution and the model distribution. For the Δ(M_BH|σ)−Δ(M_*|σ|) relation, we obtain a significance level Q_KS = 0.24, indicating that the former data distribution is compatible with the model predictions. For Δ(M_BH|M_*)−Δ(σ|M_*), we obtain Q_KS = 0.1. This value leaves open the possibility that either the sample of local BHs is still incomplete and/or that some element is still missing in the cosmological models that generates the dependence between M_BH and σ (see also the results by Barausse et al. 2017, who claimed the absence of the correlation Δ(M_BH|M_*) and Δ(σ|M_*) in both hydrodynamic simulations and SAMs).

In this context, we expect BH masses to be correlated with different DM and baryonic properties of galaxies, all emerging from the above complex interplay. To investigate this point in closer detail, we computed the relation between the BH mass and the global properties of galaxies, such as the total stellar mass content and the DM mass. We stress that these relations are expected in all cosmological models of galaxy formation because all baryonic processes (gas cooling and accretion, inflows, star formation, and feedback) are ultimately related to the depth and growth history of the gravitational potential wells. In this framework, the growth of the stellar, BH, and DM components are intertwined so that M_BH, M_*, and M_h are expected to be related to each other. The connection among these quantities is investigated in Fig. 11 and is compared with recent observational measurements by Marasco et al. (2021). These authors identified a sample of 55 nearby galaxies with dynamically measured M_BH ≥ 10⁶ M_⊙ for which dynamical measurements of M_h were available, based either on globular cluster dynamics (for galaxies of earlier Hubble types) or on spatially resolved rotation curves (for galaxies of later Hubble types). The correlation we find between the BH mass M_BH and the total stellar mass content M_* agrees excellently with observations in both strengths and scatter, the latter being characterized by an intrinsic value ϵ ≈ 0.45. As expected, we find an extremely well-defined correlation between M_BH and M_h, in agreement with the results in Marasco et al. (2021) and with previous authors (Ferrarese 2002; Pizzella et al. 2005; Volonteri et al. 2011; Davis et al. 2019a,b; Shankar et al. 2020a,b; Smith et al. 2021) who used the v_rot as a proxy for the DM mass M_h. Although our model includes the merging of BHs as a key process, our predicted M_BH − M_* relation does not show any sign of flattening at high masses. This feature has been suggested to be indicative of merging events as the leading process in the building up of the relation for early-type massive galaxies (see Graham & Sahu 2023). This suggest that the gas physics is relevant in the build-up of the relation even for massive galaxies.

Fig. 11. Predicted relations between MBH and the stellar and DM properties of the host galaxies. The MBH − M* relation (panel A), the MBH − Mh relation (panel B), the f* − MBH relation (panel C), and the f* − Mh relation (panel D) are compared with the data from Marasco et al. (2021; orange points). In panel A we also report data from Sahu et al. (2019a,b; black dots). The gray contours show the predictions of our model in the absence of AGN feedback. In all panels and for each abscissa value, the color code corresponds to the logarithm of the fraction of galaxies with different values of the quantities in the y-axis.

The tightness of the relation is directly related to the presence of AGN feedback. In a model without AGN feedback, the relation between M_BH and M_* is expected to be much broader, as shown by the gray contours in panel A of Fig. 11.

Thus, the observed tight correlation between M_BH and M_* results from the interplay between the growth of DM halos (and the ensuing collapse and cooling of gas in the DM potential wells, and the related dynamical processes affecting the distribution of gas) and the counteracting effects of AGN outflows, which keeps the relation tight.

Because all baryonic processes entering the BH growth are ultimately related to the growth of the DM halos, we expect a tight correlation between the BH mass M_BH and the DM halo mass of the host M_H. This is shown in panel B and is compared with data. The relations resulting from the model appear to be tighter than the data distribution. However, the large error bars of the data points currently do not allow us to validate or falsify the model prediction. Here we note that recent observations of central galaxies in galaxy clusters show an extremely tight relation between the BH and the halo mass (Gasparri et al. 2019), tighter indeed than the M_BH − M* relations. This supports our prediction of a strong correlation between M_BH and M_h, as found in ab initio models and cosmological simulations (see Davies et al. 2023 for recent results, and references therein).

However, as stressed by the authors above, the complex physics of baryons in the growing DM potential wells results in a number of correlations. Like all cosmological models of galaxy formation, our SAM predicts a strict relation between the depth of the gravitational potential wells and the star formation processes (panel D). In particular, the star formation efficiency f_* ≡ M_*/f_b M_h (here f_b is the Universal baryon fraction) is predicted to have a maximum for DM masses of about M_h ≈ 10¹² M_⊙. For lower values of M_h, star formation is suppressed by the strong effects of stellar feedback, which effectively expels a relevant fraction of gas from the shallow potential wells, while in massive halos with M_h ≳ 10¹² M_⊙, the inefficiency of gas cooling and strong AGN winds combine to quench star formation, leading to a declining f_* (see also Shankar et al. 2006). Thus, matching the observed behavior (panel D) provides a key test for our implemented 2-dimensional AGN feedback model. It is important to stress that such a model provides a combined correlation between the quenching properties of AGN feedback (the quantity f_*) and its function of self-regulation in the growth of BHs, as shown in panel C. The consistency between the predicted and the observed relation between M_BH and f_* provides an important consistency check for our implemented AGN feedback model because it provides a simultaneous quantitative description of both the star formation quenching and the regulation of the BH growth that is consistent with available measurements.

This correlation between the AGN feedback and the star formation properties of galaxies is expected to impact over the distribution of BHs with different degree of activity in the M_BH − M_* plane. This is investigated in detail in Fig. 12, where we show the M_BH − M* relation for different BH Eddington ratios λ, and compare them with data for active AGNs with BH masses measured through either reverberation mapping or maser techniques (van den Bosch 2016), so that they are not affected by uncertainties related to f_vir. As noted by previous authors (e.g., Reines & Volonteri 2015; Shankar et al. 2019), these objects appear to lie below the relation for the global BH population. Interestingly, although the model predicts BH masses of active galaxies somewhat higher than observed, it predicts that objects with higher Eddington ratios lie systematically along the lower envelope of the global M_BH − M* relation. The physical interpretation of this trend is that extremely high accretion rates only take place in objects in which large gas reservoirs are still available and have not been converted into stars and high BH masses at higher redshifts. Although the distribution is consistent with available data, we note that a detailed verification of our prediction would require a larger collection of observational points, which would enable us to sample the distribution of Eddington ratios with sufficient statistical significance. The apparent discrepancy/offset in the normalizations of the M_BH − M* scaling relations characterizing local AGN data and local inactive BH samples with dynamical mass measurements is a matter of intense debate (e.g., Reines & Volonteri 2015; Shankar et al. 2019; and references therein), and is not yet solved. Although it is beyond the scope of the current work to discuss this difference in depth, we note that the dispersion predicted by our SAM is more contained than the one observed in the combined data sets, possibly suggesting the presence of some additional biases in the BH and/or the stellar mass estimates between the two BH samples.

Fig. 12. Average Eddington ratio of BHs in the MBH − M* plane (color-coded). The blue data are taken from Savorgnan & Graham (2016), and the red squares and triangles are the data obtained from reverberation mapping and maser techniques, respectively. These are taken from van den Bosch (2016). To the error bars reported by the latter authors (referring to measurement errors in the K-band magnitude), we have added (in quadrature) a systematic error estimate of 0.3 dex corresponding to the uncertainties in the conversion from k-band magnitudes to stellar masses.

Finally, we investigated the predictions of our updated model on the evolution of the M_BH − M_* relation. Specifically, we compared our results with the recent measurements at 1.2 ≤ z ≤ 1.7 based on 32 X-ray selected broad-line AGN performed by Ding et al. (2020). By applying modern tools to decompose the HST images including available ACS data, the authors measured the host galaxy stellar mass through the two-dimensional model fitting, while the BH mass was determined using the broad Hα line, detected in the near-infrared with the Subaru Fiber Multi-Object Spectrograph.

To ensure a fair comparison with these observations, we reproduced the selection biases when considering the model galaxies. Specifically, we adopted the same selection window as in the observations λ_Edd ≥ −0.5 − 0.11 * (log(M_BH/M_⊙)−7.5, based on the Eddington ratio λ_Edd (see Fig. 1 of Ding et al. 2020). In addition, following Ding et al. (2020), we injected random noise to the sample of model galaxies to mimic the scatter in the data due to measurement errors, that is, δM_BH = 0.4 dex and δM_* = 0.17 dex, so as to allow for a quantitative comparison between the observed and predicted scatter in the M_BH − M_* relations.

The result of the comparison is shown in Fig. 13, where we show the prediction for both the whole sample of model galaxies and for the subsample that was selected to mimic the observational measurement errors and selection effects. When the latter are considered, the distribution of model galaxies matches the observations in both central value and scatter. The former shows a small positive evolution ΔM_BH ≈ 0.5 with respect to the local relation, which is consistent with results from N-body simulation (see, e.g., Sijacki et al. 2015; Khandai et al. 2015; De Graf et al. 2015). On the other hand, the scatter takes values ranging from ϵ ≈ 0.3 to ϵ ≈ 0.38, which excellently agrees with the observed value ϵ ≈ 0.35.

Fig. 13. Model predictions for the evolution of the MBH − M* relation. Left panel: comparison of the observed (orange dots) and predicted MBH − M* relation for galaxies in the redshift range 1.2 ≤ z ≤ 1.7. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M* bin. The data are taken from Ding et al. (2020). The total predicted distribution of model galaxies in the MBH − M* plane is shown by the blue contours, and the red region corresponds to the subsample of model galaxies that satisfies the same selection criteria as the observations. Right panel: distribution of the residuals for the observations (orange histogram), the model galaxies (blue histogram), and the model galaxies resulting in simulations with the isotropic feedback model implemented in previous versions of our SAM (green histogram, taken from Ding et al. 2020). Model galaxies are selected so as to mimic the observational measurement errors and selection effects (see text).

The latter constitutes a remarkable result that can entirely be ascribed to our improved description of AGN feedback. The same comparison with observations was performed in Ding et al. (2020) based on the SAM with the simplified isotropic description of AGN feedback. The authors reported a much larger scatter of ϵ ≈ 0.7. The distribution of the residual based on M_* (i.e., along the x-axis) Δ(M_*) is shown on the right of Fig. 12, where the observed distribution is compared with the predictions from our SAM with the new 2-dimensional treatment of feedback (Sect. 2.1), and with those referring to our SAM with the previous isotropic treatment of AGN feedback, showing that the latter are characterized by a much wider distribution of residuals. Thus, our comparative test based on one numerical model and altering the AGN feedback prescription enlightens the role of feedback in regulating the scatter of the observed relations not only at low redshift (Fig. 9), but also at higher redshifts z ≈ 1.5.

5. Summary and conclusions

The true role of AGN feedback in shaping galaxies is still a matter of intense debate both observationally and theoretically. To solve this fundamental problem, we implemented a new physical treatment of AGN-driven winds in our semi-analytic model of galaxy formation. This is based on a two-dimensional description for the expansion of AGN-driven shocks in realistic galactic gaseous disks, which has been successfully tested on a wide set of observed outflows (Menci et al. 2019). We associated solutions for the outflow expansion and the mass outflow rates in different directions with respect to the plane of the disk with each galaxy in our model (Fig. 1), depending on the AGN luminosity and on the circular velocity and gas content of the considered galaxy (Figs. 2–3). In addition, we also updated our semi-analytic model by self-consistently computing the stellar velocity dispersion via detailed Jeans modeling for each galaxy in our SAM.

We applied our SAM to the study of the scaling relations between the mass M_BH of central black holes and several properties of galactic bulges, such as the stellar mass and the velocity dispersion of the bulge, and some global properties of the model galaxies, such as the total stellar content, the dark matter mass, and the star formation efficiency. Our main results can be summarized as follows:

– The implementation in the SAM of the observational limited-resolution effects allowed us to confirm that a true correlation between the BH mass M_BH and galaxy properties (bulge mass and velocity dispersion, total stellar mass, and DM mass) does exist (Figs. 8–9), and that the observed relations do not consist of the upper envelope of the scatter of a distribution of BH masses that extends to lower BH masses, as suggested by some authors (see, e.g., Batcheldor 2010). In this context, the model results are consistent with the description of BH growth that emerges from hydrodynamical N-body simulations, such as the Magneticum Pathfinder SPH simulations (Steinborn et al. 2015), the evolution and assembly of galaxies and their environments suite of SPH simulations (Schaye et al. 2015), the Illustris moving mesh simulation (Vogelsberger et al. 2014; Sijacki et al. 2015; Li et al. 2020), and the SIMBA simulation (Thomas et al. 2019).

– When the above limited-resolution selection bias is included, the model predictions yield a normalization of the relations that slightly increased (by a factor ≈1.5 − 2) toward higher BH masses, as expected, bringing them in better agreement with observations in both normalization and scatter (Figs. 7–8).

– Our predicted BH-galaxy scaling relations are flatter than the observed ones, in particular the M_BH − σ relation, even though we adopted an improved kinetic AGN feedback model (Fig. 8). Our current work thus does not support the view that AGN feedback is the main driver of the origin of a steep M_BH − σ relation.

– Compared to the isotropic treatment of AGN feedback implemented in the previous version of our SAM, the effect of a 2-dimensional treatment of AGN outflows in realistic galactic disks is that the scatter of the observed relations is reduced (Fig. 7). Thus, the reduction (from ϵ ≈ 0.55 to ϵ ≈ 0.3 in the local Universe) is larger at high redshifts (from ϵ ≈ 0.75 to ϵ ≈ 0.3). The resulting scatter of the M_BH − M_* relation is approximately constant out to z ≈ 1.5.

– We predict a small positive evolution ΔM_BH ≈ 0.5 in the normalization of the average M_BH − M_* relation at z ≈ 1.5 with respect to the local relation (Fig. 13), which is consistent with results from N-body simulation (see, e.g., Sijacki et al. 2015; Khandai et al. 2015; De Graf et al. 2015).

– In our SAM, none of the considered galactic properties can be considered as really fundamental in determining the observed scaling relations, not even for the stellar velocity dispersion. In the framework of cosmological galaxy evolution models such as our SAM, the BH mass stems from the complex interplay between the growth history of halo mass and the physics of baryons, which cannot be reduced to a single fundamental dependence on a single quantity. To test this prediction, we analyzed the residuals of the scaling relations (Fig. 10) and performed a Kolmogorov-Smirnov test to probe the consistency of their distribution with the distribution of data points. For the Δ(M_BH|σ)−Δ(M_*|σ|) relation, we obtain a significance level Q_KS = 0.24, indicating that the former data distribution is compatible with the model predictions, while for Δ(M_BH|M_*)−Δ(σ|M_*), we obtain Q_KS = 0.1. This value leaves open the possibility that either the sample of local BHs is still incomplete, and/or that some element is still missing in the cosmological models that generates the dependence between M_BH and σ.

– The model predicts a correlation between the BH mass M_BH and the total DM mass M_h of the host galaxy, which is consistent with available data (Fig. 11, panel B), with a small scatter that is comparable to the scatter of the correlations with M_* or σ. As found in several N-body simulations (see, e.g., Davies et al. 2023), the main trigger for BH growth is constituted by galaxy merging. When disk instabilities are switched off in our model, the correlations we find between the BH mass and the galaxy properties are not affected. However, we do not find breaks in the M_BH.

– The implemented AGN feedback model provides a quenching of star formation in massive halos (Fig. 11, panel D) that is consistent with current observations.

– Our treatment of the AGN feedback allowed our SAM not only to yield correlations between the BH mass and both M_* and M_h separately, but also to provide relations with the efficiency of star formation in different DM halos f_* ≡ M_*/f_b M_h (Fig. 11, panel C). Thus, it provides a simultaneous quantitative description of the star formation quenching and the regulation of the BH growth that is consistent with available measurements.

The implementation of the two-dimensional treatment of AGN feedback in our model opens the way to several additional studies. For instance, the physical treatment of outflows in a realistic two-dimensional geometry allows us to compute the escape fraction of ionizing photons associated to each galaxy of the model. This will be tightly related to the opening angle (Fig. 2) computed for each combination of galaxy quantities (M_gas, L_AGN and M_h). We plan to investigate this issue in a future paper and to address the contribution of AGN to the reionization of the Universe (see Giallongo et al. 2012) based on a more realistic and reliable description of AGN outflows.

Acknowledgments

We thank the referee for helpful comments. We acknowledge grants from MIUR (PRIN MIUR contract 2017PH3WAT). F.S. acknowledges partial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 860744.

References

All Figures

Fig. 1. Properties of shock expansion for our reference galaxy. In the top panel we show the positions of the shock at the different times represented by different colors and displayed in the right bar. The central and bottom panels show the velocity map and mass outflow rate map, respectively. The values corresponding to the colored contours are displayed in the bars. The X and Y coordinates correspond to the distance from the galaxy center in the directions parallel and perpendicular to the plane of the disk, respectively.

In the text

	Fig. 2. Properties of the gas outflows. The color code shows the dependence of the fraction of cold gas that is expelled by AGN outflows (fexpell, top panel) and of the solid angle θc swept out by AGN outflows (bottom panel) on the AGN bolometric luminosity LAGN and of the cold gas mass Mcold. A DM circular velocity Vc = 200 km s−1 is assumed for both panels.
In the text

	Fig. 3. Average fraction of gas expelled by AGN outflows in galaxies with different stellar mass M* and star formation rate Ṁ∗ shown by the color code.
In the text

	Fig. 4. Predicted AGN bolometric luminosity functions (solid lines) are compared to the compilation of data in Shen et al. (2020, shaded regions) in six redshifts bins shown in the labels.
In the text

Fig. 5. Predicted distribution of the effective radius Re shown as a function of the stellar mass M* at redshifts z = 0.25, z = 1.51, and z = 2.25 (from left to right) for early-type (top row) and late-type galaxies (bottom row). The predictions are compared with the data from van der Wel et al. (2014; shaded areas). We also show the size-stellar mass relation measured at z = 0.1 by Bernardi et al. (2014) and Mowla et al. (2019) as dashed lines and dots, respectively.

In the text

Fig. 6. Predicted relation between the velocity dispersion σ and the stellar mass M* for passive galaxies (see text). The color code represents the logarithm of the fraction of model galaxies with the considered σ for each bin of M*. The shaded band shows the relation for early-type galaxies selected on the basis of morphological indicators after Hyde & Bernardi (2009). The dashed black and red lines show the average relation derived by Zahid et al. (2016) for passive SDSS galaxies and for the Smithsonian Hectospec Lensing Survey (SHELS) sample, selected on the basis of spectral properties (Dn4000 > 1.5, where Dn4000 is the flux ratio of two spectral windows adjacent to the 4000 Å break (see Balogh et al. 1999).

In the text

Fig. 7. Effect of feedback treatment on the MBH − M*, b relation. The relation predicted for galaxies at z ≤ 0.1 by our new two-dimensional model for AGN feedback (top panels) presented in Sect. 2.1 is compared with the predictions obtained with our previous isotropic treatment of feedback (bottom panels). For both AGN feedback models, we show the results with and without the selection effect discussed in Sect. 2.1, as indicated by the labels. In all panels, the color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M*, b bin. The standard deviations ϵ for the new and previous treatment of AGN feedback is shown as a function of bulge stellar mass in the rightmost histograms. The horizontal strip represents the range of values obtained from different observational works: the upper bound is the intrinsic scatter ϵ derived by De Nicola et al. (2019), and the lower bound is the value of ϵ derived by the least-square fit analysis by Kormendy & Ho (2013) using individual errors in log σ, adding individual errors in log MBH to the intrinsic scatter in quadrature, and iterating the intrinsic scatter until the reduced χ2 = 1.

In the text

Fig. 8. Comparison of our predicted MBH − M*, b relation with the data for all galaxies at z ≤ 0.1. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M*, b bin. Panel A: Predicted relation between the BH mass MBH and the bulge stellar mass M*, b. The data points are from Kormendy & Ho (2013; black squares and triangles for elliptical galaxies and bulges, respectively), Savorgnan & Graham (2016; red dots), and Ho & Kim (2014; orange filled squares). Panel B: Same as panel A, but including the selection bias discussed in Sect. 3, i.e., considering only model galaxies satisfying the condition in Eq. (12) with θcrit = 0.1. Panel C: Same as panel A, but considering only model galaxies hosting AGN with bolometric luminosity LAGN ≥ 1044 erg s−1. Bottom right panel: Rms value of the scatter of the distributions shown in panel A (brown histogram), panel B (blue histogram), and panel C (orange histogram) for different bins of M*, b. The horizontal strip corresponds to the range of measured intrinsic scatter by different authors. The lower value ϵ = 0.28 was obtained by Kormendy & Ho (2013) and the upper value ϵ = 0.4 was obtained by Savorgnan et al. (2016).

In the text

Fig. 9. Comparison of our predicted MBH − σ relation with data for all galaxies at z ≤ 0.1. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log σ bin. Panel A: Predicted relation between the BH mass MBH and the bulge velocity dispersion σ brightness-weighted averaged over the galaxy. The data points are from Kormendy & Ho (2013; black squares and triangles for elliptical galaxies and bulges, respectively), De Nicola et al. (2019; brown squares), Savorgnan & Graham (2016; red dots), and Ho & Kim (2014; orange filled squares). Panel B: Same as panel A, but including in the model the selection bias discussed in Sect. 3, i.e., considering only model galaxies that satisfy the condition in Eq. (12) with θcrit = 0.1. Panel C: Same as panel A, but considering only model galaxies hosting AGN with bolometric luminosity LAGN ≥ 1044 erg/s. Bottom right panel: Histograms represent the rms values of the scatter for the distributions shown in the top left (brown), top right (blue), and bottom left panel (orange). The horizontal strip represents the range of values obtained from different observational works: the upper bound is the intrinsic scatter ϵ derived by De Nicola et al. (2019), and the lower bound is the value of ϵ derived by the least-square fit analysis by Kormendy & Ho (2013) using individual errors in log σ, adding individual errors in log MBH to the intrinsic scatter in quadrature, and iterating the intrinsic scatter until the reduced χ2 = 1.

In the text

Fig. 10. Predictions for the residuals of the BH-galaxy scaling relations. In the left panel, we show the correlation between the residuals of the MBH − M*, b relations and those of the σ − M* relation at fixed stellar mass. The correlation between the residuals of the MBH − σ relation and those of the M* − σ relation at fixed velocity dispersion are shown in the right panel. Data points refer to the sample of Savorgnan et al. (2016). The contours refer to model galaxies with a bulge-to-total stellar mass ratio B/T > 0.8 for which the black hole sphere of influence can be resolved.

In the text

Fig. 11. Predicted relations between MBH and the stellar and DM properties of the host galaxies. The MBH − M* relation (panel A), the MBH − Mh relation (panel B), the f* − MBH relation (panel C), and the f* − Mh relation (panel D) are compared with the data from Marasco et al. (2021; orange points). In panel A we also report data from Sahu et al. (2019a,b; black dots). The gray contours show the predictions of our model in the absence of AGN feedback. In all panels and for each abscissa value, the color code corresponds to the logarithm of the fraction of galaxies with different values of the quantities in the y-axis.

In the text

Fig. 12. Average Eddington ratio of BHs in the MBH − M* plane (color-coded). The blue data are taken from Savorgnan & Graham (2016), and the red squares and triangles are the data obtained from reverberation mapping and maser techniques, respectively. These are taken from van den Bosch (2016). To the error bars reported by the latter authors (referring to measurement errors in the K-band magnitude), we have added (in quadrature) a systematic error estimate of 0.3 dex corresponding to the uncertainties in the conversion from k-band magnitudes to stellar masses.

In the text

Fig. 13. Model predictions for the evolution of the MBH − M* relation. Left panel: comparison of the observed (orange dots) and predicted MBH − M* relation for galaxies in the redshift range 1.2 ≤ z ≤ 1.7. The color code corresponds to the logarithm of the fraction of galaxies with different MBH in a given log M* bin. The data are taken from Ding et al. (2020). The total predicted distribution of model galaxies in the MBH − M* plane is shown by the blue contours, and the red region corresponds to the subsample of model galaxies that satisfies the same selection criteria as the observations. Right panel: distribution of the residuals for the observations (orange histogram), the model galaxies (blue histogram), and the model galaxies resulting in simulations with the isotropic feedback model implemented in previous versions of our SAM (green histogram, taken from Ding et al. 2020). Model galaxies are selected so as to mimic the observational measurement errors and selection effects (see text).

In the text