# Selection bias in dynamically-measured supermassive black hole samples: Scaling relations and correlations between residuals in semi-analytic galaxy formation models

Enrico Barausse, Francesco Shankar, Mariangela Bernardi, Yohan Dubois, Ravi K. Sheth

► To cite this version:

Enrico Barausse, Francesco Shankar, Mariangela Bernardi, Yohan Dubois, Ravi K. Sheth. Selection bias in dynamically-measured supermassive black hole samples: Scaling relations and correlations between residuals in semi-analytic galaxy formation models. Monthly Notices of the Royal Astronomical Society, 2017, 468 (4), pp.4782-4791. 10.1093/mnras/stx799. hal-01555149

HAL Id: hal-01555149

https://hal.science/hal-01555149v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Selection bias in dynamically measured supermassive black hole samples: scaling relations and correlations between residuals in semi-analytic galaxy formation models

Enrico Barausse,<fnref n="1" />★ Francesco Shankar,<fnref n="2" /> Mariangela Bernardi,<fnref n="3" /> Yohan Dubois<fnref n="1" /> and Ravi K. Sheth<fnref n="3" />

Accepted 2017 March 27. Received 2017 March 23; in original form 2017 February 7

## ABSTRACT

Recent work has confirmed that the scaling relations between the masses of supermassive black holes and host-galaxy properties such as stellar masses and velocity dispersions may be biased high. Much of this may be caused by the requirement that the black hole sphere of influence must be resolved for the black hole mass to be reliably estimated. We revisit this issue with a comprehensive galaxy evolution semi-analytic model. Once tuned to reproduce the (mean) correlation of black hole mass with velocity dispersion, the model cannot account for the correlation with stellar mass. This is independent of the model's parameters, thus suggesting an internal inconsistency in the data. The predicted distributions, especially at the low-mass end, are also much broader than observed. However, if selection effects are included, the model's predictions tend to align with the observations. We also demonstrate that the correlations between the *residuals* of the scaling relations are more effective than the relations themselves at constraining models for the feedback of active galactic nuclei (AGNs). In fact, we find that our model, while in apparent broad agreement with the scaling relations when accounting for selection biases, yields very weak correlations between their residuals at fixed stellar mass, in stark contrast with observations. This problem persists when changing the AGN feedback strength, and is also present in the hydrodynamic cosmological simulation Horizon-AGN, which includes state-of-the-art treatments of AGN feedback. This suggests that current AGN feedback models are too weak or simply not capturing the effect of the black hole on the stellar velocity dispersion.

**Key words:** black hole physics – galaxies: fundamental parameters – galaxies: nuclei – quasars: supermassive black holes – galaxies: structure.

## 1 INTRODUCTION

Supermassive black holes, with masses $M_{bh} \sim 10^6–10^9 M_{\odot}$, have been identified at the centres of all local galaxies observed with high enough sensitivity (see e.g. Ferrarese & Ford 2005; Shankar 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Graham 2016, for reviews). A surprising finding that has puzzled astrophysicists for the last 40 yr or so is that the masses of these black holes appear to be tightly linked to the global properties of their hosts, such as stellar mass or velocity dispersion, defined on scales up to 1000 times the sphere of influence of the central black hole. The origin of these

correlations is still hotly debated, though there is general agreement that understanding this origin will shed light on the more general and still unsolved problem of the formation of galaxies.

Supermassive black holes are thought to have formed in a highly star-forming, gas-rich phase at early cosmological epochs. Central ‘seed’ black holes are thought to gradually grow via mainly gas accretion, eventually becoming massive enough to shine as quasars or Seyfert galaxies and trigger powerful winds and/or jets capable of removing gas and quenching star formation in the host galaxy. This feedback from active black holes has become a key ingredient in almost all galaxy evolution models (e.g. Granato et al. 2004; Cirasuolo et al. 2005; Vittorini, Shankar & Cavaliere 2005; Croton et al. 2006; Hopkins et al. 2006; Lapi et al. 2006; Shankar et al. 2006; Monaco, Fontanot & Taffoni 2007; Guo et al. 2011; Barausse 2012; Dubois

<footnote n="1">¹Institut d'Astrophysique de Paris, UMR 7095, CNRS and UPMC, Sorbonne Universités, 98bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="2">²Department of Physics and Astronomy, University of Southampton, Highfield SO17 1BJ, UK</footnote>

<footnote n="3">³Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd St, Philadelphia, PA 19104, USA</footnote>

<footnote>* E-mail: barausse@iap.fr</footnote>
et al. 2012a, 2016; Fontanot, Monaco & Shankar 2015; Bower
et al. 2016). At later times, both the host galaxy and its black hole
may further increase their mass (and size) via mergers with other
galaxies/black holes, which could contribute up to ~80 per cent of
their final mass (e.g. De Lucia & Blaizot 2007; Malbon et al. 2007;
Oser et al. 2010, 2012; Shankar et al. 2010; Shankar et al. 2013;
González et al. 2011; Dubois et al. 2013, 2016; Rodriguez-Gomez
et al. 2016; Welker et al. 2017). Additional mechanisms, besides
mergers, can also contribute to the growth of the stellar bulge and
feeding of the central black hole, most notably disc instabilities (e.g.
Bower et al. 2006; Bournaud et al. 2011; Barausse 2012; Di Matteo
et al. 2012; Dubois et al. 2012b).

Recently, Shankar et al. (2016, see also Bernardi et al. 2007,
Gültekin et al. 2011, Morabito & Dai 2012 and van den Bosch
et al. 2015) showed that the local sample of galaxies with dynami-
cal mass measurements of supermassive black holes is biased. Local
galaxies hosting supermassive black holes, irrespective of their ex-
act morphological type or the aperture within which the velocity
dispersion is measured, typically present velocity dispersions that
are substantially larger than those of a very large and unbiased sam-
ple from the Sloan Digital Sky Survey (SDSS) with similar stellar
masses. One of the main reasons for this bias can be traced back to
the observationally imposed requirement that the black hole gravi-
tational sphere of influence must be resolved for the black hole mass
to be reliably estimated. Via dedicated Monte Carlo simulations and
accurate analysis of the residuals around the observed black hole
scaling relations, Shankar et al. (2016) found the velocity disper-
sion to be a more fundamental quantity than stellar mass or effective
radius. Indeed, the observed black hole scaling relation involving
the stellar mass was found to be much more biased than the one
involving velocity dispersion (up to an order of magnitude in nor-
malization), and its apparent tightness could be entirely ascribed to
a selection effect.

Shankar et al. (2016) also suggested that a selection bias more
prominent in stellar mass than in velocity dispersion may explain
several discrepancies often reported in the literature, i.e. the fact that
the observed relation between black hole and stellar mass predicts a
local black hole mass density two to three times higher than inferred
from the relation between black hole mass and velocity dispersion
(e.g. Graham et al. 2007; Tundo et al. 2007; Shankar, Weinberg &
Miralda-Escudé 2009). Shankar, Bernardi & Sheth (2017) further
extended the comparison between the set of local galaxies with
dynamically measured black hole masses and SDSS galaxies. They
found evidence that even the correlation between black hole mass
and Sérsic index, recently claimed to be even tighter than the one
with velocity dispersion (Savorgnan 2016), is severely biased, with
the correlation to velocity dispersion remaining more fundamental.
The bias in the local scaling relations could also have profound
implications for the background of gravitational waves expected
from binary supermassive black holes, which could be a factor of
a few lower than what current pulsar timing arrays can effectively
detect (Sesana et al. 2016).

The aim of this work is to revisit the local scaling relations
between the masses of supermassive black holes and host-galaxy
properties, namely velocity dispersion and stellar mass, in the con-
text of a comprehensive semi-analytic model of galaxy formation
and evolution (Barausse 2012). This model also evolves supermas-
sive black holes self-consistently from high-redshift ‘seeds’, and
accounts for black hole mergers and for the feedback from ac-
tive galactic nuclei (AGNs). After briefly reviewing the model in
Section 2, we discuss (Sections 3.1 and 3.2) the slope, normaliza-
tion and scatter of the black hole scaling relations with and without

the aforementioned selection effect on the resolvability of the black
hole sphere of influence. In Section 3.3, we study the correlations
between the residuals from fitted scaling relations and show that
they are useful for constraining theoretical models such as ours as
well as the hydrodynamic cosmological simulation Horizon-AGN
(Dubois et al. 2014, 2016; Volonteri et al. 2016). Finally, we will
discuss our results in Section 4 and summarize our conclusions in
Section 5.

## 2 MODEL

The full description of the semi-analytic model adopted as a ref-
erence in this work can be found in Barausse (2012), with later
updates of the prescriptions for the black hole spin and nuclear star
cluster evolution described, respectively, in Sesana et al. (2014) and
Antonini, Barausse & Silk (2015a,b). Here, we briefly summarize
the key points about the growth of the central supermassive black
holes, which are the main focus of this paper.

The model is built on top of dark-matter merger trees gen-
erated via extended Press–Schechter algorithms (e.g. Press &
Schechter 1974; Parkinson, Cole & Helly 2008) tuned to reproduce
the results of *N*-body simulations (Parkinson et al. 2008). Galaxies
form in each halo via the interplay and balance of gas cooling, star
formation and (supernova) feedback. Dark matter haloes are also ini-
tially seeded with either *light* black holes of $M_{\text{seed}} \sim 200 M_{\odot}$ (to be
interpreted e.g. as the remnants of PopIII stars), or with *heavy* black
holes of mass $M_{\text{seed}} \sim 10^5 M_{\odot}$, which may arise for instance from
protogalactic disc instabilities. The seeding of haloes is assumed
to happen at early epochs $z > 15$, with halo occupation fractions
depending on the specific seeding model (see Barausse 2012; Klein
et al. 2016, for details).

In our model, seed black holes initially grow via (mainly) gas accretion from a gas reservoir, which is in turn assumed to form at a rate proportional to the bulge star formation rate (e.g. Granato et al. 2004; Lapi et al. 2006). As a result, the feeding of this reservoir and the ensuing black hole accretion events typically happen after star formation bursts triggered by major galactic mergers and disc instabilities. In both their radiatively efficient (‘quasar’) and inefficient (‘radio’) accretion modes, the black holes also exert a feedback on the host galaxies, thus reducing their (hot and cold) gas content and quenching star formation. As discussed by a number of groups (Granato et al. 2004; Cirasuolo et al. 2005), AGN feedback prescriptions such as these tend to induce a correlation between black hole mass and velocity dispersion of the bulge component. Also accounted by the model is the black hole growth via black hole mergers, following the coalescence of the host galaxies. This mechanism becomes particularly important for high black hole masses at recent epochs.

The model is calibrated against a set of observables, such as the
local stellar and black hole mass functions, the local gas fraction,
the star formation history, the AGN luminosity function, the local
morphological fractions and the correlations between black holes
and galaxies and between black holes and nuclear star clusters (c.f.
Barausse 2012; Sesana et al. 2014; Antonini et al. 2015a,b). In more
detail, as we will show in the following (c.f. Fig. 2), the model’s
default calibration attempts to match the observed $M_{bh}-\sigma$ relation
without accounting for any observational bias (on morphological
type or on the resolvability of the black hole influence sphere).

## 3 RESULTS

We will now compare the predictions of our model with observa-
tions, focusing on the normalization of the scaling relations and
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: velocity dispersion as a function of total stellar mass for the SDSS galaxies with $P(E+S0) > 0.8$ (long-dashed red line), with its $1\sigma$ dispersion (grey), compared with the prediction of the light-seed model for bulge-dominated/elliptical galaxies, i.e. galaxies with $B/T > 0.7$ (solid line, with dotted lines marking the 70 per cent confidence region). The model's predictions for heavy seeds are very similar. Right: stellar mass function of the SDSS galaxies based on Sérsic plus Exponential fits to the observed surface brightness (shaded area, accounting for the uncertainties due to the stellar population modelling, fitting and assumptions about dust in the galaxies, c.f. Bernardi et al. 2013; Bernardi et al. 2017). Solid and dotted lines show the prediction of the model with light seeds and its (Poissonian) $1\sigma$ uncertainties. The predictions for heavy seeds are very similar.</figcaption></figure>

on the role played by selection biases; the dispersion around the scaling relations; and the correlations between the residuals of the data from the scaling relations.

## 3.1 The normalization of the scaling relations and the role of selection bias

The left-hand panel of Fig. 1 shows the relation between velocity dispersion and stellar mass for early-type galaxies in the SDSS. Here, 'early-type' means that the probability of being elliptical or lenticular, $p(E+S0)$, according to the automatic morphological classification of Huertas-Company et al. (2011), exceeds 0.8. We restrict to this specific SDSS subsample as velocity dispersions in late-type galaxies are not spatially resolved, though the correlation does not depend on the exact cut in $p(E+S0)$. For consistency with the data of Savorgnan et al. (2016) to which we will compare, we follow Shankar et al. (2016) and correct the velocity dispersions $\sigma_{HL}$, as in Cappellari et al. (2006), to a common aperture of 0.595 kpc (i.e. the one adopted by the Hyperleda data base; Paturel et al. 2003). Henceforth, unless stated otherwise, we will always define velocity dispersions $\sigma$ at the aperture of Hyperleda. Stellar masses $M_{star}$ are from Bernardi et al. (2013). They are the product of luminosity $L$ and mass-to-light ratio $M_{star}/L$; the $L$ values are from Meert, Vikram & Bernardi (2015), based on Sérsic+Exponential fits to the light profiles.

The black solid line marks the median velocity dispersion-stellar mass relation as predicted by the model (for bulge-dominated/elliptical galaxies only) and black dotted lines show the 15th and 85th percentiles of the predicted distribution (at fixed stellar mass). Central velocity dispersions in the model are computed as $\sigma = A\sqrt{GM_b/r_b}[1 + (V_b/\sigma)^2]$, where $M_b$ is the bulge dynamical mass, $r_b$ is the scale radius of the Hernquist profile (which the model adopts to describe the bulge, see Barausse 2012), $A \approx 0.4$ accounts for the anisotropy of the distribution function of the bulge stellar population (c.f. Baes & Dejonghe 2002, fig. 2, lower panel) and the ratio $V_b/\sigma$ accounts for the contribution of the bulge rotation and is modelled based on observations (c.f. Sesana et al. 2014, for details). As can be seen, the predicted correlation is similar to the observed one, although slightly flatter.

For completeness, the right-hand panel of Fig. 1 compares the stellar mass function predicted by the model with the observed one

(Bernardi et al. 2013; Bernardi et al. 2017). While the model lies slightly above the data at the highest masses, it lies below over the range of $2 \times 10^{10} \lesssim M_{star}/M_{\odot} \lesssim 10^{11}$. This is not a major issue in the present context, since the model is consistent with the empirical galaxy scaling relations, and most notably with the $\sigma-M_{star}$ relation shown in the left-hand panel.

Having checked that our model reproduces the dynamical scaling relations of early-type galaxies, we now study the scaling relations with the central supermassive black hole. Fig. 2 compares the model (with no morphological selection) with the whole sample (i.e. spirals as well as ellipticals and lenticulars) of Savorgnan et al. (2016, blue diamonds). The left- and right-hand panels show the scaling of black hole mass with velocity dispersion and total stellar mass, respectively. In each panel, the solid and dotted lines show the median and the region containing 70 per cent of the model objects in a given bin. The top panels show the full black hole sample, while the middle and bottom panels only show the subset for which the sphere of influence exceeds the typical (HST) resolution limit,

$$r_{\text{infl}} \equiv k \frac{GM_{\text{BH}}}{\sigma^2}, \quad \frac{r_{\text{infl}}}{d_{\text{Ang}}} > 0.1'' \quad (1)$$

$(d_{Ang}$ being the angular diameter distance). We use the parameter $k$ to take into account different galaxy mass profiles: $k \sim 4$ for the Hernquist profiles assumed by the model (see Barausse 2012, for details), but $k \sim 10$ (or even larger) is possible if a core is present. On the other hand, strong lensing and accurate dynamical modelling have shown that the mass profiles of intermediate-mass, early-type galaxies are consistent with nearly isothermal profiles down to (at least) tenths of the effective radius (e.g. Cappellari et al. 2015, and references therein). These have $k \sim 1$. To bracket these uncertainties, we show results for both $k = 10$ and 1.

To match the data as closely as possible, we draw the angular diameter distances $d_{Ang}$ from an empirical probability distribution function, which we construct from the distances in the Savorgnan et al. (2016) sample, and which peaks around $\sim 15-20$ Mpc. Using a distribution that is uniform in comoving volume, as done in Shankar et al. (2016), yields comparable results.

The top panels of Fig. 2 suggest that the model predicts intrinsic scaling relations that lie slightly below the data, especially for $M_{bh}-M_{star}$ (top right), but that are broadly consistent with the intrinsic relations suggested by Shankar et al. (2016, long dashed red lines).
<figure><img src="image_4.png" /></figure>

**Figure 2.** Black hole mass as a function of velocity dispersion (left-hand panels) and total stellar mass (right-hand panels) as predicted by our model with light seeds (the results for heavy seeds are qualitatively similar). Results are shown for the full outputs of the model (labelled as 'Intrinsic', top panels), and for the subsample of black holes with a gravitational sphere of influence above 0.1 arcsec (labelled as 'Observed', middle and bottom panels; see the text for details). The solid and dotted blue lines in all panels represent the medians and 70 per cent confidence region of the distributions. The dashed magenta lines in the left-hand panels are the median model predictions when assigning velocity dispersions to galaxies from the observed SDSS $\sigma_{\text{HL}}-M_{\text{star}}$ relation (long-dashed red line in the left-hand panel of Fig. 1). Long-dashed red lines are the intrinsic scaling relations as inferred from Monte Carlo simulations by Shankar et al. (2016). Blue diamonds are data collected and updated by Savorgnan et al. (2016) on local galaxies with dynamical measurements of supermassive black holes. Note that observational biases can increase the normalization and reduce the scatter in the intrinsic scaling relations.

The match to the $M_{\text{bh}}-\sigma$ relation improves (slightly) if we assign velocity dispersions by using the SDSS $\sigma-M_{\text{star}}$ relation, e.g. via the analytic fits provided by Sesana et al. (2016, dashed magenta line, left-hand panels).<fnref n="1" /> At the same time, the model substantially under-predicts the $M_{\text{bh}}-M_{\text{star}}$ relation by up to an order of magnitude (top right panel), suggesting that there is some internal inconsistency in the data. In other words, a model like ours, tuned to match the local velocity dispersion-stellar mass and black hole mass-velocity dispersion relations, tends to severely underpredict the $M_{\text{bh}}-M_{\text{star}}$ relation.

This is in line with the results of Shankar et al. (2016). Calibrating the model to match the observed $M_{\text{bh}}-M_{\text{star}}$ relation would instead overestimate the observed $M_{\text{bh}}-\sigma$ relation. Such an overestimate has indeed been seen in two (very different) cosmological hydrodynamic simulations (Sijacki et al. 2015; Volonteri et al. 2016).<fnref n="2" />

When the selection effect on the sphere of influence of the black hole is applied to the model (the middle and bottom panels of Fig. 2),

<footnote n="1">¹ Note that it makes sense to assign $\sigma$ from the model-predicted $M_{\text{star}}$ via the observed $\sigma-M_{\text{star}}$ relation, rather than $M_{\text{star}}$ from the model-predicted $\sigma$, because masses are more ‘primitive’ quantities than velocity dispersions for a semi-analytic galaxy formation model such as ours.</footnote>

<footnote n="2">² Volonteri et al. (2016) mention resolution, which is at best 1 kpc in their simulation, as one reason why their results do not match the $M_{\text{bh}}-\sigma$ relation. However, their predictions for $M_{\text{bh}}$ are larger than the observations even at large $\sigma$ (c.f. their fig. 7), while they are in good agreement with the $M_{\text{bh}}-M_{\text{star}}$ relation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Fig. 2, but for galaxies with a bulge-to-total ratio $B/T > 0.7$ (and for light seeds; the results are similar for heavy seeds). The results are broadly similar to Fig. 2, though the dispersion of the model’s intrinsic predictions is substantially lower than for the full galaxy sample. Note that the observational data set is also restricted to early-type galaxies only, to match the sample of simulated galaxies.</figcaption></figure>

the median normalizations of the predicted scaling relations increase (especially for the $M_{bh}-M_{star}$ relation), because a substantial fraction of low-mass black holes are excluded. This is because, for a given angular aperture, equation (1) preferentially removes objects with the smallest gravitational spheres of influence; these tend to be the lowest mass black holes. Therefore, this effect tends to select the ‘upper end’ (in black hole mass) of the intrinsic distributions shown in the top panels of Fig. 2. This also induces an overall flattening of the scaling relations, which is again more obvious in the $M_{bh}-M_{star}$ plane: selection hardly matters for the most massive galaxies, but it causes a factor of $\lesssim 10$ increase in the median observed $M_{bh}$ at lower masses. Selection-biased models are flatter in the $M_{bh}-\sigma$ plane as well.<fnref n="3" />

We conclude that to agree with resolution-biased observed scaling relations, models must predict intrinsic scaling relations that are significantly steeper than the observed relations. We have tried to obtain steeper intrinsic scaling relations in our model by changing the AGN feedback, but we have found this to be insufficient to achieve better agreement with the data at low masses and velocity dispersions. In fact, the results and conclusions of this paper are robust to changes in the AGN feedback strength as well as to changes in the black hole accretion prescriptions.

Fig. 2 also clearly shows that in addition to changing the normalization and slope, selection dramatically decreases the dispersion around the median relations, especially at lower masses. Note that these ‘corrections’ do not depend on the exact choice of $k$, since they are present for both $k = 10$ and $k = 1$.

Similar comments apply to Fig. 3, which compares model-predicted galaxies with bulge-to-total ratios of $B/T > 0.7$ with the E+S0 galaxies of Savorgnan et al. (2016), in the same format as Fig. 2. The intrinsic distributions of the model-predicted galaxies (top panels) are narrower than in the full samples (c.f. the top panels

<footnote n="3">³ As a result of this flatter slope, the model’s prediction (after applying the selection bias) lies above the (few) data with $\sigma \sim 10^2$ km s$^{-1}$ in the sample of Savorgnan et al. (2016). Note however that other samples, such as that of Kormendy, Bender & Cornell (2011), contain black holes with masses up to $\sim 10^8$ M$_\odot$ at $\sigma \sim 10^2$ km s$^{-1}$, which is in better agreement with our model.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Median black hole mass as a function of total host-galaxy stellar mass (solid lines) as predicted by the model, for the black holes with bolometric luminosity $log L/\text{erg s}^{-1} > 42$. The dotted (long-dashed) lines mark the 90 per cent (99 per cent) confidence regions at given stellar mass, as predicted by our model. Blue squares are data from Reines & Volonteri (2015). The results are for the light-seed model, but the heavy-seed one gives very similar results.</figcaption></figure>

of Fig. 2), and the median scaling relations have higher normalizations, in better agreement with the data and in line with the findings of Barausse (2012). Although the selection effect is smaller for this specific subsample of galaxies, the model-predicted intrinsic $M_{bh}-\sigma$ relation is offset slightly above the data, while the predictions for the intrinsic $M_{bh}-M_{\text{star}}$ relation are slightly below the data. The effect of including the selection bias on the resolvability of the black hole sphere of influence (middle and bottom panels) is less important than in Fig. 3, because the small systems for which the sphere of influence is not resolvable tend to live in late-type galaxies in our semi-analytic model. Nevertheless, the selection bias still tends to make the correlations higher in normalization, slightly flatter and slightly tighter.

## 3.2 The dispersion around the scaling relations: a comparison with observations

It is important to emphasize that whatever the exact black hole seeding recipe adopted in the model, the predicted local (intrinsic) scaling relations reported in Fig. 2 for all galaxy types show very large dispersions, especially for galaxy stellar masses $M_{\text{star}} \lesssim 3 \times 10^{11} M_{\odot}$, with a broad distribution up to two orders of magnitude or more in black hole mass at fixed stellar mass or velocity dispersion. This is because the model tends to retain a significant number of low-mass black holes, which did not accrete much gas along their history (because they live in spirals or satellite galaxies) and which therefore remain closer to their seed masses (e.g. Volonteri et al. 2005; Barausse 2012). However, for more massive galaxies the dispersion is smaller, with almost no galaxies with $M_{\text{star}} \gtrsim 5 \times 10^{11} M_{\odot}$ having black holes with $M_{bh} \lesssim 10^7 M_{\odot}$.

One interesting way to probe the existence of very low mass black holes in relatively low-mass galaxies could be to compare with the scaling relations in active galaxies, which are not limited by the spatial resolution issues that heavily affect dynamical measurements of black holes. To this purpose, Fig. 4 compares the predictions of our model with the recent sample of 262 broad-line AGNs collected by Reines & Volonteri (2015). In more detail, the observations are represented by blue squares, while the lines

represent the median, the 90 per cent confidence region (i.e. the 5th and 95th percentiles) and the 99 per cent confidence region (i.e. the 0.5th and 99.5th) of the (model-predicted) $M_{bh}-M_{\text{star}}$ relation, by assuming a light black hole seed scenario (the heavy-seed scenario gives very similar results) and considering only systems with bolometric luminosity $log (L/\text{erg s}^{-1}) > 42$ (roughly the minimum luminosity probed by Reines & Volonteri 2015). The model's distribution of active black holes has been built by randomly drawing Eddington ratios from a Schechter distribution that extends up to the Eddington limit, in agreement with a number of observations (Kauffmann & Heckman 2009; Aird et al. 2012; Bongiorno et al. 2012; Schulze et al. 2015; Jones et al. 2016). We have verified that our predicted luminosity function at $z=0$, computed by assuming an average duty cycle of active black holes of 10 per cent consistent with the results from local surveys (e.g. Goulding & Alexander 2009; Shankar et al. 2013; Pardo et al. 2016, and references therein), agrees with the (obscuration-corrected) bolometric luminosity functions of Hopkins, Richards & Hernquist (2007) and Shankar et al. (2009).

First, let us note, as emphasized by Shankar et al. (2016), that a lower limit of $L \gtrsim 10^{42} \text{ erg s}^{-1}$ should still allow black holes down to a mass of $M_{bh} \sim 10^4 M_{\odot}$ to be detected, at least if a non-negligible fraction of these black holes are still accreting at the Eddington limit. Such low-mass black holes do not seem to exist in the Reines & Volonteri (2015, see also Baldassare et al. 2015) sample (and in our model, at least in sufficiently large numbers and with high enough Eddington ratios to warrant detection). Even assuming lower virial factors $f_{\text{vir}}$ than those adopted by Reines & Volonteri (2015) in deriving black hole masses from their measured FWHMs, as suggested by some groups (e.g. Shankar et al. 2016; Yong, Webster & King 2016, and references therein), would not alter these conclusions.

Secondly, it is clear that the observational sample lies, on average, below the model median predictions, which were tuned to reproduce the data on inactive local galaxies with a significantly higher normalization. While the predicted 90 per cent and 99 per cent confidence regions for the active black hole population encompass the data of Reines & Volonteri (2015), the model also predicts the existence of a large number of active higher mass black holes (above the median), which are not observed. Therefore, either the sample of Reines & Volonteri (2015) is biased towards low-luminosity active systems, or our model should be normalized to lower values by a factor of $\gtrsim 3$.

The large scatter in our model means that if the normalization is decreased, then our model would predict a large tail of very low mass black holes. Since these are not observed, this would have important consequences for constraining models of the seeds of the supermassive black hole population. On the other hand, decreasing the normalization of the $M_{bh}-M_{\text{star}}$ relation is by no means straightforward. While it could be achieved by decreasing black hole accretion, this would imply a proportional reduction in the AGN luminosity, unless a higher radiative efficiency and/or duty cycle are also assumed. The present calibration of our model already predicts rather high radiative efficiencies/black hole spins (Sesana et al. 2014). Duty cycles of active black holes could instead be constrained by comparing with independent AGN clustering measurements (Gatti et al. 2016, and references therein). We plan to explore some of these important interrelated issues in future work.

## 3.3 Constraints from correlations between residuals

We now compare our model's predictions for the *residuals* of the black hole-galaxy scaling relations to the data. Such correlations are
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Correlation between the residuals of the $M_{bh}-M_{star}$ relations and those of the $\sigma-M_{star}$ relation, at fixed stellar mass (left-hand panel), and between the residuals of the $M_{bh}-\sigma$ relation and those of the $M_{star}-\sigma$ relation, at fixed velocity dispersion (right-hand panel). Red circles and green triangles show the ellipticals and lenticulars in the Savorgnan et al. (2016) sample. The solid blue and long-dashed magenta lines mark the best fits (for the model's predictions and the observations, respectively), while the dotted lines show the $1\sigma$ uncertainty on the slope. Also reported are the best fits and the Pearson correlation coefficient $r$. The grey bands show the residuals extracted from the Monte Carlo simulations from Shankar et al. (2016), with selection bias on the black hole gravitational sphere of influence. Note that the correlation of the residuals at fixed stellar mass (left-hand panel) is very strong in the data, but essentially absent in the model. These results are for the light-seed model (the heavy-seed one yields similar results), selecting only early-type galaxies ($B/T > 0.7$) for which the black hole sphere of influence is resolvable.</figcaption></figure>

an efficient way of going beyond pairwise correlations between the variables themselves (Bernardi et al. 2005; Sheth & Bernardi 2012). For example, measurements of the $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ correlations alone do not provide insight about whether $\sigma$ is more important than $M_{star}$ in determining $M_{bh}$. This is because the $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ correlations do not encode complete information about the joint distribution of $M_{bh}$, $\sigma$ and $M_{star}$. Correlations between the residuals encode some of this extra information.

To this purpose, the left- and right-hand panels of Fig. 5 show $\Delta(M_{bh}|M_{star})$ versus $\Delta(\sigma|M_{star})$ and $\Delta(M_{bh}|\sigma)$ versus $\Delta(M_{star}|\sigma)$, where

$$ \Delta(Y|X) \equiv \log Y - \langle\log Y|\log X\rangle \quad (2) $$

is the residual in the $Y$ variable (at fixed $X$) from the log-log-linear fit of $Y(X)$ versus $X$, i.e. $\langle\log Y|\log X\rangle$. The magenta long-dashed and dotted lines in each panel show the best fit and $1\sigma$ uncertainties on the correlations between residuals in the Savorgnan et al. (2016) data set: red circles and green triangles represent ellipticals and lenticulars. We obtained the magenta lines by running 200 iterations following the steps outlined in Shankar et al. (2017), which include errors in both variables. At each iteration, we eliminate three random objects from the original sample. From the full ensemble of realizations, we measure the mean slope and its $1\sigma$ uncertainty.

The blue solid and dotted lines show a similar analysis in our semi-analytic model. However, in this case, we randomly produce 30 mock samples of ~75 galaxies each, with $B/T > 0.7$ and resolvable black hole spheres of influence. From the full ensemble of mock realizations, we then extract the mean slope and Pearson coefficient.

The correlations in Fig. 5 show that, in the data, the velocity dispersion is more strongly correlated with the black hole mass, with a mean Pearson coefficient of $r = 0.72$, than stellar mass, for which the Pearson coefficient is $r = 0.56$. The model instead predicts just the opposite, with almost zero correlation with velocity dispersion (mean $r = 0.05$), but with a correlation with stellar mass consistent with the data, though still rather weak (mean $r = 0.28$). It is nevertheless important to realize that even an intrinsically weak

correlation with velocity dispersion at fixed stellar mass does not necessarily imply that the total correlation with velocity dispersion is small. In fact, following Appendix B in Shankar et al. (2017), the total dependence of the black hole mass on velocity dispersion can be summarized as $M_{bh} \propto \sigma^{\beta} M_{star}^{\alpha} \propto \sigma^{\beta+\alpha}\gamma$, where $\gamma$ comes from $M_{star} \propto \sigma^{\gamma}$ (where we have ignored any other explicit dependence on, e.g. Sérsic index or galaxy size). Since the model predicts $\gamma \approx 4$, and the residuals in Fig. 5 (solid blue lines) yield $\beta \sim 0.5$ and $\alpha \sim 0.6$, one obtains a total dependence of $M_{bh} \propto \sigma^3$, consistent with Fig. 3 (middle left panel). Also note that a slope $\gamma \approx 4$ for the $M_{star}-\sigma$ relation is in itself already in tension with the ($V_{max}$-corrected) SDSS observations, which rather prefer a slope $\gamma \approx 2.2-2.5$, at least for $M_{star} \gtrsim 2 \times 10^{10} M_{\odot}$ (c.f. also Fig. 7).

For completeness, Fig. 5 also shows results from the Monte Carlo simulations performed by Shankar et al. (2016, grey bands), which assume an intrinsic correlation of the type $M_{bh} \propto \sigma^{4.5} M_{star}^{0.5}$ and include the selection bias on $r_{infl}$. These models reproduce the observed trends.

Several comments are in order. First, our results are robust against the efficiency of the AGN feedback. Indeed, in our reference model, we fixed the AGN feedback efficiency to obtain a stellar mass function as close as possible to the observations at the high-mass end (c.f. Fig. 1, right-hand panel). By slightly decreasing this efficiency (e.g. by a factor of ~3), the agreement with the observed stellar mass function in Fig. 1 slightly worsens in the high-mass end, while the $M_{bh}-\sigma$ relation steepens<fnref n="4" /> ($M_{bh} \propto \sigma^{3.4}$ after including the effect of the selection bias), but not enough to fully match the observations (for which slopes ~4.5–5 are usually quoted in the literature, e.g. Graham et al. 2011; Kormendy & Ho 2013). More importantly, the correlations of the residuals remain almost unchanged ($r = 0.1$ when fixing $M_{star}$, and $r = 0.3$ when fixing $\sigma$). Secondly, we have checked that these correlations are not improved by considering a

<footnote n="4">⁴ Note that in our model, increasing the AGN feedback efficiency decreases the slope of the $M_{bh}-\sigma$ relation, mainly because AGN feedback is more effective at clearing the gas in more massive galactic hosts, thus inhibiting the growth of especially the more massive black holes.</footnote>
<figure><img src="image_8.png" /></figure>

**Figure 6.** Same as Fig. 5, but for the hydrodynamic cosmological simulation Horizon-AGN (Dubois et al. 2014, 2016). However, unlike in Fig. 5, no selection bias on the resolvability of the black hole sphere of influence has been applied. Doing so leads to slightly weaker correlations and slightly lower slopes for the fits to the residuals. Note that the velocity dispersions in the simulation are measured within the effective radius of the galaxy, and are not corrected to the Hyperleda aperture.

<figure><img src="image_9.png" /></figure>

**Figure 7.** Total stellar mass as a function of velocity dispersion for the SDSS galaxies with $P(E-SO) > 0.8$ (shaded region) and for the Horizon-AGN simulation (the solid line represents the mean and the dotted lines represent the 70 per cent confidence region; Dubois et al. 2014, 2016). As in the observational sample, we select early-type galaxies from the simulations by only considering systems with $V_c/\sigma < 0.7$, where $V_c$ is the rotational velocity. Applying a different cut does not change this figure significantly. Here, for both the SDSS and Horizon-AGN data, $\sigma_e$ is computed within the effective radius. Given the resolution of Horizon-AGN, this quantity is more reliably estimated than the more central velocity dispersion $\sigma_{HL}$ adopted in the previous figures.

different choice of $k$ in equation (1); by assuming no bias on the resolvability of the black hole sphere of influence; by considering all galaxies rather than just bulge-dominated ones; by using the bulge mass instead of the total stellar mass or by using the SDSS $\sigma-M_{star}$ relation (e.g. by the fits of Sesana et al. 2016) to compute velocity dispersions. Third, the model's residuals shown in Fig. 5 do not account for measurement errors (i.e. we compute the residuals for the model's exact predictions, without folding in any observational uncertainties). We have verified that including these errors can yield a steeper slope $\beta \sim 2$ for the residuals at fixed $M_{star}$, but this does not strengthen the correlations of the residuals shown above (essentially because the error on $\beta$ also grows when $\beta$ grows).

To check if our results are due to the particular implementation of black hole accretion and AGN feedback in our semi-analytic model, and/or lack of sufficient 'coupling' between velocity dis-

persion and AGN feedback, we have performed a similar analysis of the Horizon-AGN simulation (Dubois et al. 2014, 2016). This is a hydrodynamic cosmological simulation of a box with size 100 Mpc $h^{-1}$, run with the adaptive mesh refinement code RAMSES (Teyssier 2002), with $10^9$ dark-matter particles and a minimum mesh size of 1 kpc. The simulation includes gas cooling, star formation, feedback from stars and AGNs (see Dubois et al. 2014 for the details of the numerical modelling and Volonteri et al. 2016 for the discussion of the correlations between galaxies and black holes in Horizon-AGN). Galaxies are extracted with a galaxy finder running on star particles. Since the simulation assumes a Salpeter (1955) IMF, we have reduced galaxy stellar masses by 0.25 dex (e.g. Bernardi et al. 2010) to match the Chabrier IMF (Chabrier 2003) adopted in this work. The velocity dispersion is measured from its components along each direction of the cylindrical coordinates oriented along the galaxy's spin axis (i.e. $\sigma_r$, $\sigma_t$ and $\sigma_z$ for the radial-, tangential- and z-component, respectively), thus $\sigma^2 = (\sigma_r^2 + \sigma_t^2 + \sigma_z^2)/3$. The velocity dispersion of each galaxy is measured using only star particles within the effective radius of the galaxy. The results are shown in Figs 6 and 7. They are in qualitative agreement with those obtained with our semi-analytic model. In more detail, Horizon-AGN also predicts $\gamma \approx 4$, in tension with the data, and a very weak correlation between the residuals when fixing $M_{star}$. Also note that in Fig. 6, we have not applied any selection bias (unlike in the case of the semi-analytic model in Fig. 5). Restricting to systems with resolvable spheres of influence actually makes the correlations in the residuals of the Horizon-AGN simulation even weaker.

We plan to consider different (and possibly stronger) models of AGN feedback in future work. For now, the discrepancies highlighted in Figs 5 and 6 clearly show that correlations between the residuals of the $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ scaling relations are more powerful than the scaling relations themselves at constraining models for the co-evolution of black holes and their host galaxies.

# 4 DISCUSSION

Concerning the normalization of the black hole scaling relations, at face value our model clearly fails to reproduce the observed $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ relations at the same time. Without invoking any selection effect, one could attempt to improve the simultaneous match
to both the observed scaling relations by fine-tuning some of the key parameters in the model, namely the one controlling gas accretion on to the black hole and/or the energetic feedback from the central active nucleus. For example, increasing the efficiency of AGN feedback could in principle decrease the stellar masses of the host galaxies at fixed black hole mass, thus possibly improving the match to the $M_{bh}-M_{star}$ relation (the top right panel of Fig. 2). However, this would then spoil the good match to the velocity dispersion-stellar mass relation. Similarly, one could increase the accretion on to the central black hole at fixed star formation rate and final stellar mass of the host galaxy, with the aim of improving the match to the $M_{bh}-M_{star}$ relation. This however would proportionally increase the $M_{bh}-\sigma$ relation above the data, when velocity dispersions are chosen to faithfully track those from the observed SDSS $\sigma-M_{star}$ relation (magenta lines in the left-hand panels of Fig. 2). Concerning the dispersion around the scaling relations, our semi-analytic model, which self-consistently evolves black holes from seeds at high redshifts, naturally predicts very broad distributions in black hole mass at fixed velocity dispersion or stellar mass, at relatively low masses $M_{star} \lesssim 3 \times 10^{11} M_{\odot}$.

Overall, these effects are in line with the results of the Monte Carlo simulations presented in Shankar et al. (2016), though the intrinsic scatters assumed there, especially in the $M_{bh}-\sigma$ relation, were always below 0.3 dex. This was needed to avoid too flat slopes in the ‘observed’, biased relations. Indeed, here our predicted slope of the $M_{bh}-\sigma$ relation, inclusive of observational bias (e.g. middle and bottom left panels of Fig. 2), is approximately $M_{bh} \propto \sigma^\beta$ with slope $\beta \sim 3-3.5$, significantly flatter than the slopes $\beta \gtrsim 4.5-5$ usually quoted in the literature (e.g. Graham et al. 2011; Kormendy & Ho 2013).

Conversely, a key point of our present work is that, irrespective of the chosen masses for the seed black holes, the model predicts relatively tight scaling relations at high stellar masses $M_{star} \gtrsim 5 \times 10^{11} M_{\odot}$. In this respect, our model does not support the conjecture put forward by Batchelor et al. (2007), according to which the $M_{bh}-\sigma$ relation is only an upper limit of a more or less uniform distribution of black holes. This is also in line with the Monte Carlo tests performed by Shankar et al. (2016, their fig. 11), in which a very broad distribution of black holes extending to the lowest masses in all types of galaxies is highly disfavoured. Nevertheless, Fig. 4 shows that the present data on active galaxies may still be consistent with our model, once the proper flux limits and Eddington ratio distributions are accounted for. Therefore, an intrinsically broad distribution for the black hole mass in relatively small galaxies with mass $M_{star} \lesssim 3 \times 10^{11} M_{\odot}$, extending down to black holes of mass $M_{bh} \gtrsim 10^5 M_{\odot}$ or so, cannot be excluded by present data, though it is highly disfavoured in more massive galaxies. Models with such broad distributions, however, still tend to produce scaling relations flatter than observed, once the bias on the black hole sphere of influence and on galaxy morphology is folded in the model predictions, as shown in Figs 2 and 3.

Finally, we stress again that our semi-analytic model fails to reproduce the observations when it comes to correlations between the *residuals* of the scaling relations. In more detail, the model is consistent with the data for the correlation between the residuals in the $M_{bh}-\sigma$ relation and those in the $M_{star}-\sigma$ relation, at fixed velocity dispersion. However, the model predicts almost no correlation between the residuals of the $M_{bh}-M_{star}$ relations and those of the $\sigma-M_{star}$ relation, at fixed stellar mass, while the data hint at a rather strong correlation. We have verified that these results are robust against changing the parameters of the model (namely the AGN feedback efficiency and the parameter regulating black hole

accretion). Moreover, we have shown that the same weak correlation between the residuals at fixed stellar mass is also obtained in the hydrodynamic cosmological simulation Horizon-AGN (Dubois et al. 2014, 2016), which includes thermal quasar-mode feedback and jet-structured radio-mode AGN feedback (i.e. processes that are expected to induce a stronger coupling between black hole mass and velocity dispersion). Another noteworthy point is that both our model and Horizon-AGN tend to produce slopes for the $\sigma-M_{star}$ relation that are significantly steeper than the data, although in Horizon-AGN numerical resolution effects may bias the measurements of the velocity dispersion in lower mass galaxies (see discussion in Dubois et al. 2016).

Our interpretation is that the weak correlation of the residuals at fixed stellar mass may be at least partly due to current AGN feedback models being possibly too weak to capture the full effect of the black hole on the stellar velocity dispersion. While this is expected in a semi-analytic model such as ours, where the AGN feedback is typically assumed to simply eject gas from the nuclear region (Granato et al. 2004; Barausse 2012), it is more surprising in the Horizon-AGN simulation, which captures the back-reaction of the ejected gas on to the stellar and dark-matter dynamics (Peirani et al. 2016). None the less, in hydrodynamic cosmological simulations such as Horizon-AGN, the spatial resolution is 1 kpc at best and may limit our capability to properly capture the interaction of AGN winds with gas and its impact on the dynamics within galactic bulges. Moreover, the Horizon-AGN simulation currently employs a rather crude model of AGN feedback, mostly based on simple gas heating, to mimic the so-called quasar-mode feedback. More realistic models of AGN feedback, also inclusive of momentum-driven winds and stronger coupling with the surrounding interstellar medium (e.g. Bieri et al. 2017), might possibly be more effective at improving the comparison with the black hole-galaxy scaling relations and their residuals. We should also mention that another possibility that has been put forward in the literature is that AGN feedback may *not* be the main cause of the scaling relations, which might be ascribed instead to a common gas supply for the galaxy and the black hole, regulated by gravitational torques (Anglés-Alcázar et al. 2017).

## 5 CONCLUSIONS

We have compared the predictions of a comprehensive semi-analytic model of galaxy formation, which self-consistently evolves supermassive black holes from high-redshift seeds by accounting for gas accretion, mergers and AGN feedback, with the observed scaling relations of the masses of supermassive black holes with stellar mass and velocity dispersion. Our main conclusions are as follows:

(i) At $M_{star} \gtrsim 5 \times 10^{11} M_{\odot}$, the dispersion in black hole mass at fixed stellar mass is $\lesssim 1$ dex – very few black holes with masses $M_{bh} \lesssim 10^7 M_{\odot}$ are predicted in such massive galaxies. However, for galaxies having $M_{star} \lesssim 3 \times 10^{11} M_{\odot}$, the distribution of $M_{bh}$ at fixed $\sigma$ or $M_{star}$ is broad.

(ii) Observational selection effects, associated with resolving the black hole sphere of influence and/or with selecting bulgedominated/elliptical galaxies, tighten the $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ scaling relations, bringing them into better agreement with the observations.

(iii) No evident variation in AGN feedback and/or black hole accretion efficiencies can provide a *simultaneous* match to both scaling relations. This supports previous work suggesting an internal inconsistency between the observed $M_{bh}-\sigma$ and $M_{bh}-M_{star}$ relations.

(iv) Galaxy-evolution models (our semi-analytic one as well as the Horizon-AGN simulation) predict almost *no* correlation
between the residuals of the $M_{bh}-M_{star}$ relations and those of the $\sigma-M_{star}$ relation, at fixed stellar mass. Since the data hint at a rather strong correlation, this calls for revamped AGN feedback recipes in the next generation of cosmological galaxy-evolution models, or for a re-assessment of the importance of gravitational torques in regulating the black hole-galaxy co-evolution (Anglés-Alcázar et al. 2017).

## ACKNOWLEDGEMENTS

We thank M. Volonteri for insightful conversations. EB acknowledges support from the H2020-MSCA-RISE-2015 grant no. StronGrHEP-690904. This work has made use of the Horizon Cluster, hosted by the Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running this cluster so smoothly. This work was partially supported by the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche.

## REFERENCES

Aird J. et al., 2012, ApJ, 746, 90

Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840

Antonini F., Barausse E., Silk J., 2015a, ApJ, 806, L8

Antonini F., Barausse E., Silk J., 2015b, ApJ, 812, 72

Baes M., Dejonghe H., 2002, A&A, 393, 485

Baldassare V. F., Reines A. E., Gallo E., Greene J. E., 2015, ApJ, 809, L14

Barausse E., 2012, MNRAS, 423, 2533

Batcheldor D., Marconi A., Merritt D., Axon D. J., 2007, ApJ, 663, L85

Bernardi M., Sheth R. K., Nichol R. C., Schneider D. P., Brinkmann J., 2005, AJ, 129, 61

Bernardi M., Sheth R. K., Tundo E., Hyde J. B., 2007, ApJ, 660, 267

Bernardi M., Shankar F., Hyde J. B., Mei S., Marulli F., Sheth R. K., 2010, MNRAS, 404, 2087

Bernardi M., Meert A., Sheth R. K., Vikram V., Huertas-Company M., Mei S., Shankar F., 2013, MNRAS, 436, 697

Bernardi M. et al., 2017, MNRAS, preprint (arXiv:1604.01036)

Bieri R. et al., 2017, MNRAS, 464, 1854

Bongiorno A. et al., 2012, MNRAS, 427, 3103

Bournaud F. et al., 2011, ApJ, 730, 4

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Bower R., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2016, preprint (arXiv:1607.07445)

Cappellari M. et al., 2006, MNRAS, 366, 1126

Cappellari M. et al., 2015, ApJ, 804, L21

Chabrier G., 2003, PASP, 115, 763

Cirasuolo M., Shankar F., Granato G. L., De Zotti G., Danese L., 2005, ApJ, 629, 816

Croton D. J. et al., 2006, MNRAS, 365, 11

De Lucia G., Blaizot J., 2007, MNRAS, 375, 2

Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C., Lopez J., Springel V., 2012, ApJ, 745, L29

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012a, MNRAS, 420, 2662

Dubois Y. et al., 2012b, MNRAS, 423, 3616

Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y. et al., 2016, MNRAS, 463, 3948

Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523

Fontanot F., Monaco P., Shankar F., 2015, MNRAS, 453, 4112

Gatti M. et al., 2016, MNRAS, 456, 1073

González J. E., Lacey C. G., Baugh C. M., Frenk C. S., 2011, MNRAS, 413, 749

Goulding A., Alexander D., 2009, preprint (arXiv:0906.0772)

Graham A. W., 2016, Galactic Bulges, 418, 263

Graham A. W., Driver S. P., Allen P. D., Liske J., 2007, MNRAS, 378, 198

Graham A. W., Onken C. A., Athanassoula E., Combes F., 2011, MNRAS, 412, 2211

Granato G. L., De Zotti G., Silva L., Bressan A., Danese L., 2004, ApJ, 600, 580

Gültekin K., Tremaine S., Loeb A., Richstone D. O., 2011, ApJ, 738, 17

Guo Q. et al., 2011, MNRAS, 413, 101

Hopkins P. F., Hernquist L., Cox T. J., DiMatteo T., Robertson B., Springel V., 2006, ApJS, 163, 1

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

Huertas-Company M., Aguerri J. A. L., Bernardi M., Mei S., Sánchez Almeida J., 2011, A&A, 525, A157

Jones M. L., Hickox R. C., Black C. S., Hainline K. N., DiPompeo M. A., Goulding A. D., 2016, ApJ, 826, 12

Kauffmann G., Heckman T. M., 2009, MNRAS, 397, 135

Klein A. et al., 2016, Phys. Rev. D, 93, 024003

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Kormendy J., Bender R., Cornell M. E., 2011, Nature, 469, 374

Lapi A., Shankar F., Mao J., Granato G. L., Silva L., De Zotti G., Danese L., 2006, ApJ, 650, 42

McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184

Malbon R. K., Baugh C. M., Frenk C. S., Lacey C. G., 2007, MNRAS, 382, 1394

Meert A., Vikram V., Bernardi M., 2015, MNRAS, 446, 3943

Monaco P., Fontanot F., Taffoni G., 2007, MNRAS, 375, 1189

Morabito L. K., Dai X., 2012, ApJ, 757, 172

Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63

Pardo K. et al., 2016, preprint (arXiv:1603.01622)

Parkinson H., Cole S., Helly J., 2008, MNRAS, 383, 557

Paturel G., Petit C., Prugniel P., Theureau G., Rousseau J., Brouty M., Dubois P., Cambrésy L., 2003, A&A, 412, 45

Peirani S. et al., 2016, preprint (arXiv:1611.09922)

Press W. H., Schechter P., 1974, ApJ, 187, 425

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Rodriguez-Gomez V. et al., 2016, preprint (arXiv:e-prints)

Salpeter E. E., 1955, ApJ, 121, 161

Savorgnan G. A. D., 2016, ApJ, 821, 88

Savorgnan G. A. D., Graham A. W., Marconi A., Sani E., 2016, ApJ, 817, 21

Schulze A. et al., 2015, MNRAS, 447, 2085

Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104

Sesana A., Shankar F., Bernardi M., Sheth R. K., 2016, MNRAS, 463, L6

Shankar F., 2009, New Astron. Rev., 53, 57

Shankar F., Lapi A., Salucci P., De Zotti G., Danese L., 2006, ApJ, 643, 14

Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20

Shankar F. et al., 2010, MNRAS, 403, 117

Shankar F., Marulli F., Bernardi M., Mei S., Meert A., Vikram V., 2013, MNRAS, 428, 109

Shankar F. et al., 2016, MNRAS, 460, 3119

Shankar F., Bernardi M., Sheth R. K., 2017, MNRAS, preprint (arXiv:1701.01732)

Sheth R. K., Bernardi M., 2012, MNRAS, 422, 1825

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Teyssier R., 2002, A&A, 385, 337

Tundo E., Bernardi M., Hyde J. B., Sheth R. K., Pizzella A., 2007, ApJ, 663, 53

van den Bosch R. C. E., Gebhardt K., Gültekin K., Yıldırım A., Walsh J. L., 2015, ApJS, 218, 10

Vittorini V., Shankar F., Cavaliere A., 2005, MNRAS, 363, 1376

Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620, 69

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Welker C. et al., 2017, MNRAS, 465, 1241

Yong S. Y., Webster R. L., King A. L., 2016, preprint (arXiv:1602.04672)

This paper has been typeset from a TEX/LATEX file prepared by the author.