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
aa r X i v : . [ a s t r o - ph . GA ] A p r Mon. Not. R. Astron. Soc. , 1– 11 (2017) Printed 17 July 2018 (MN L A TEX style file v2.2)
Selection bias in dynamically-measured supermassive black holesamples: Scaling relations and correlations between residuals insemi-analytic galaxy formation models
Enrico Barausse ⋆ , Francesco Shankar † , Mariangela Bernardi , Yohan Dubois and Ravi K. Sheth Institut d’Astrophysique de Paris, UMR 7095, CNRS & UPMC, Sorbonne Universit´es, 98bis Bd Arago, F-75014 Paris, France Department of Physics and Astronomy, University of Southampton, Highfield, SO17 1BJ, UK Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd St, Philadelphia, PA 19104
ABSTRACT
Recent work has confirmed that the scaling relations between the masses of supermassiveblack holes and host-galaxy properties such as stellar masses and velocity dispersions maybe biased high. Much of this may be caused by the requirement that the black-hole sphereof influence must be resolved for the black-hole mass to be reliably estimated. We revisitthis issue with a comprehensive galaxy evolution semi-analytic model. Once tuned to repro-duce the (mean) correlation of black-hole mass with velocity dispersion, the model cannotaccount 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, especiallyat the low-mass end, are also much broader than observed. However, if selection effects areincluded, the model’s predictions tend to align with the observations. We also demonstratethat the correlations between the residuals of the scaling relations are more effective than therelations themselves at constraining AGN feedback models. In fact, we find that our model,while in apparent broad agreement with the scaling relations when accounting for selectionbiases, yields very weak correlations between their residuals at fixed stellar mass, in starkcontrast with observations. This problem persists when changing the AGN feedback strength,and is also present in the hydrodynamic cosmological simulation Horizon-AGN, which in-cludes state-of-the-art treatments of AGN feedback. This suggests that current AGN feedbackmodels are too weak or simply not capturing the effect of the black hole on the stellar velocitydispersion.
Key words: (galaxies:) quasars: supermassive black holes – galaxies: fundamental parame-ters – galaxies: nuclei – galaxies: structure – black hole physics
Supermassive black holes, with masses M bh ∼ – M ⊙ , havebeen identified at the centres of all local galaxies observed withhigh enough sensitivity (see, e.g., Ferrarese & Ford 2005; Shankar2009; Kormendy & Ho 2013; Graham 2016; McConnell & Ma2013, for reviews). A surprising finding that has puzzled astro-physicists for the last forty years or so is that the masses of theseblack holes appear to be tightly linked to the global properties oftheir hosts, such as stellar mass or velocity dispersion, defined onscales up to a thousand times the sphere of influence of the centralblack hole. The origin of these correlations is still hotly debated,though there is general agreement that understanding this origin ⋆ E-mail: [email protected] † E-mail: [email protected] will shed light on the more general and still unsolved problem ofthe formation of galaxies.Supermassive black holes are thought to have formed in ahighly star-forming, gas-rich phase at early cosmological epochs.Central “seed” black holes are thought to gradually grow viamainly gas accretion, eventually becoming massive enough toshine as quasars or Seyfert galaxies and trigger powerful windsand/or jets capable of removing gas and quenching star forma-tion in the host galaxy. This feedback from active black holeshas become a key ingredient in almost all galaxy evolution mod-els (e.g., Granato et al. 2004; Cirasuolo et al. 2005; Vittorini et al.2005; Croton et al. 2006; Hopkins et al. 2006; Lapi et al. 2006;Shankar et al. 2006; Monaco et al. 2007; Guo et al. 2011; Barausse2012; Dubois et al. 2012a, 2016; Fontanot et al. 2015; Bower et al.2016). At later times, both the host galaxy and its black hole c (cid:13) E. Barausse et al. may further increase their mass (and size) via mergers withother galaxies/black holes, which could contribute up to ∼ of their final mass (e.g., De Lucia & Blaizot 2007; Malbon et al.2007; Oser et al. 2010; Shankar et al. 2010; Oser et al. 2012;Gonz´alez et al. 2011; Shankar et al. 2013; Dubois et al. 2013,2016; Rodriguez-Gomez et al. 2016; Welker et al. 2017). Ad-ditional mechanisms, besides mergers, can also contribute tothe growth of the stellar bulge and feeding of the centralblack hole, most notably disc instabilities (e.g., Bower et al.2006; Bournaud et al. 2011; Di Matteo et al. 2012; Barausse 2012;Dubois et al. 2012b).Recently, Shankar et al. (2016, see also Bernardi et al.2007, G¨ultekin et al. 2011, Morabito & Dai 2012 andvan den Bosch et al. 2015) showed that the local sample ofgalaxies with dynamical mass measurements of supermassiveblack holes is biased. Local galaxies hosting supermassive blackholes, irrespective of their exact morphological type or the aperturewithin which the velocity dispersion is measured, typically presentvelocity dispersions that are substantially larger than those of avery large and unbiased sample from the Sloan Digital Sky Survey(SDSS) with similar stellar masses. One of the main reasonsfor this bias can be traced back to the observationally imposedrequirement that the black-hole gravitational sphere of influencemust be resolved for the black-hole mass to be reliably estimated.Via dedicated Monte Carlo simulations and accurate analysis ofthe residuals around the observed black-hole scaling relations,Shankar et al. (2016) found the velocity dispersion to be a morefundamental quantity than stellar mass or effective radius. Indeed,the observed black-hole scaling relation involving the stellar masswas found to be much more biased than the one involving velocitydispersion (up to an order of magnitude in normalisation), and itsapparent tightness could be entirely ascribed to a selection effect.Shankar et al. (2016) also suggested that a selection bias moreprominent in stellar mass than in velocity dispersion may ex-plain several discrepancies often reported in the literature, i.e. the fact that the observed relation between black-hole and stel-lar mass predicts a local black-hole mass density two to threetimes higher than inferred from the relation between black-holemass and velocity dispersion (e.g., Graham et al. 2007; Tundo et al.2007; Shankar et al. 2009). Shankar et al. (2017) further extendedthe comparison between the set of local galaxies with dynamicallymeasured black-hole masses and SDSS galaxies. They found evi-dence that even the correlation between black-hole mass and S´ersicindex, recently claimed to be even tighter than the one with veloc-ity dispersion (Savorgnan 2016), is severely biased, with the cor-relation to velocity dispersion remaining more fundamental. Thebias in the local scaling relations could also have profound impli-cations for the background of gravitational waves expected frombinary supermassive black holes, which could be a factor of a fewlower than what current pulsar timing arrays can effectively detect(Sesana et al. 2016).The aim of the present work is to revisit the local scal-ing relations between the masses of supermassive black holesand host-galaxy properties, namely velocity dispersion and stel-lar mass, in the context of a comprehensive semi-analytic modelof galaxy formation and evolution (Barausse 2012). This modelalso evolves supermassive black holes self-consistently from high-redshift “seeds”, and accounts for black-hole mergers and for thefeedback from active galactic nuclei (AGNs). After briefly review-ing the model in Section 2, we discuss (Sections 3.1 and 3.2)the slope, normalisation and scatter of the black-hole scaling re-lations with and without the aforementioned selection effect on the resolvability of the black-hole sphere of influence. In Sec-tion 3.3 we study the correlations between the residuals from fit-ted scaling relations, and show that they are useful for constrain-ing theoretical models such as ours as well as the hydrodynamiccosmological simulation Horizon-AGN (Dubois et al. 2014, 2016;Volonteri et al. 2016). Finally, we will discuss our results in Sec-tion 4 and summarise our conclusions in Section 5. The full description of the semi-analytic model adopted as a ref-erence in this work can be found in Barausse (2012), with laterupdates of the prescriptions for the black-hole spin and nuclearstar cluster evolution described respectively in Sesana et al. (2014)and Antonini et al. (2015a,b). Here, we briefly summarise the keypoints 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 treesgenerated via extended Press-Schechter algorithms (e.g.,Press & Schechter 1974; Parkinson et al. 2008) tuned to re-produce the results of N-body simulations (Parkinson et al. 2008).Galaxies form in each halo via the interplay and balance of gascooling, star formation and (supernova) feedback. Dark matterhaloes are also initially seeded with either light black holes of M seed ∼ M ⊙ (to be interpreted e.g. as the remnants of PopIIIstars), or with heavy black holes of mass M seed ∼ M ⊙ , whichmay arise for instance from protogalactic disc instabilities. Theseeding of haloes is assumed to happen at early epochs z > ,with halo occupation fractions depending on the specific seedingmodel (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 toform 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 ofthis reservoir and the ensuing black-hole accretion events typi-cally happen after star formation bursts triggered by major galac-tic mergers and disc instabilities. In both their radiatively effi-cient (“quasar”) and inefficient (“radio”) accretion modes, the blackholes also exert a feedback on the host galaxies, thus reducing their(hot and cold) gas content and quenching star formation. As dis-cussed by a number of groups (Granato et al. 2004; Cirasuolo et al.2005), AGN feedback prescriptions such as these tend to induce acorrelation between black hole mass and velocity dispersion of thebulge component. Also accounted by the model is the black-holegrowth via black-hole mergers, following the coalescence of thehost galaxies. This mechanism becomes particularly important forhigh black-hole masses at recent epochs.The model is calibrated against a set of observables, suchas the local stellar and black-hole mass functions, the local gasfraction, the star-formation history, the AGN luminosity function,the local morphological fractions, and the correlations betweenblack holes and galaxies and between black holes and nuclear starclusters (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 ob-served M bh - σ relation without accounting for any observationalbias (on morphological type or on the resolvability of the black-hole influence sphere). c (cid:13) , 1– 11 MBHs: selection bias and models Figure 1.
Left: Velocity dispersion as a function of total stellar mass for SDSS galaxies with P ( E + S > . (long-dashed red line), with its 1 σ dispersion(grey), compared with the prediction of the light-seed model for bulge-dominated/elliptical galaxies, i.e. galaxies with B/T > . (solid line, with dottedlines marking the 70% confidence region). The model’s predictions for heavy seeds are very similar. Right: Stellar mass function of SDSS galaxies based onS´ersic 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 lightseeds and its (Poissonian) σ uncertainties. The predictions for heavy seeds are very similar. We will now compare the predictions of our model with observa-tions, focusing on the normalisation of the scaling relations and onthe role played by selection biases; the dispersion around the scal-ing relations; and the correlations between the residuals of the datafrom the scaling relations.
The left panel of Figure 1 shows the relation between velocity dis-persion and stellar mass for early-type galaxies in the SDSS. Here‘early-type’ means that the probability of being elliptical or lentic-ular, p (E+S0), according to the automatic morphological classifica-tion of Huertas-Company et al. (2011), exceeds 0.8. We restrict tothis specific SDSS subsample as velocity dispersions in late-typegalaxies are not spatially resolved, though the correlation does notdepend on the exact cut in p (E+S0). For consistency with the dataof Savorgnan et al. (2016) to which we will compare, we followShankar et al. (2016) and correct the velocity dispersions σ HL , asin 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 velocitydispersions σ at the aperture of Hyperleda. Stellar masses M star arefrom 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 et al.(2015), based on S´ersic+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 the15th and 85th percentiles of the predicted distribution (at fixedstellar mass). Central velocity dispersions in the model are com-puted as σ = A p GM b /r b [1 + ( V b /σ ) ] , where M b is the bulgedynamical mass, r b is the scale radius of the Hernquist profile(which the model adopts to describe the bulge, see Barausse 2012), A ≈ . accounts for the anisotropy of the distribution functionof the bulge stellar population (c.f. Baes & Dejonghe 2002, Fig-ure 2, lower panel), and the ratio V b /σ accounts for the contribu-tion of the bulge rotation and is modeled based on observations (c.f. Sesana et al. 2014, for details). As can be seen, the predictedcorrelation is similar to the observed one, although slightly flatter.For completeness, the right panel of Figure 1 compares thestellar mass function predicted by the model with the observed one(Bernardi et al. 2013; Bernardi et al. 2017). While the model liesslightly above the data at the highest masses, it lies below overthe range × . M star /M ⊙ . . This is not a majorissue in the present context, since the model is consistent with theempirical galaxy scaling relations, and most notably with the σ - M star relation shown in the left panel.Having checked that our model reproduces the dynamicalscaling relations of early-type galaxies, we now study the scal-ing relations with the central supermassive black hole. Figure 2compares the model (with no morphological selection) with thewhole sample (i.e., spirals as well as ellipticals and lenticulars) ofSavorgnan et al. (2016, blue diamonds). The left and right panelsshow the scaling of black hole mass with velocity dispersion andtotal stellar mass, respectively. In each panel, the solid and dottedlines show the median and the region containing 70% of the modelobjects in a given bin. The top panels show the full black-hole sam-ple, while the middle and bottom panels only show the subset forwhich the sphere of influence exceeds the typical (HST) resolutionlimit, r infl ≡ k GM bh σ , r infl d Ang > . ′′ (1)( d Ang being the angular-diameter distance). We use the parame-ter k to take into account different galaxy mass profiles: k ∼ for the Hernquist profiles assumed by the model (see Barausse2012, for details), but k ∼ (or even larger) is possible ifa core is present. On the other hand, strong lensing and accu-rate dynamical modelling have shown that the mass profiles ofintermediate-mass, early-type galaxies are consistent with nearlyisothermal profiles down to (at least) tenths of the effective radius(e.g., Cappellari et al. 2015, and references therein). These have k ∼ . To bracket these uncertainties, we show results for both k = 10 and 1.To match the data as closely as possible, we draw the an-gular diameter distances d Ang from an empirical probability dis-tribution function, which we construct from the distances in theSavorgnan et al. (2016) sample, and which peaks around ∼ –20 c (cid:13)000
The left panel of Figure 1 shows the relation between velocity dis-persion and stellar mass for early-type galaxies in the SDSS. Here‘early-type’ means that the probability of being elliptical or lentic-ular, p (E+S0), according to the automatic morphological classifica-tion of Huertas-Company et al. (2011), exceeds 0.8. We restrict tothis specific SDSS subsample as velocity dispersions in late-typegalaxies are not spatially resolved, though the correlation does notdepend on the exact cut in p (E+S0). For consistency with the dataof Savorgnan et al. (2016) to which we will compare, we followShankar et al. (2016) and correct the velocity dispersions σ HL , asin 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 velocitydispersions σ at the aperture of Hyperleda. Stellar masses M star arefrom 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 et al.(2015), based on S´ersic+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 the15th and 85th percentiles of the predicted distribution (at fixedstellar mass). Central velocity dispersions in the model are com-puted as σ = A p GM b /r b [1 + ( V b /σ ) ] , where M b is the bulgedynamical mass, r b is the scale radius of the Hernquist profile(which the model adopts to describe the bulge, see Barausse 2012), A ≈ . accounts for the anisotropy of the distribution functionof the bulge stellar population (c.f. Baes & Dejonghe 2002, Fig-ure 2, lower panel), and the ratio V b /σ accounts for the contribu-tion of the bulge rotation and is modeled based on observations (c.f. Sesana et al. 2014, for details). As can be seen, the predictedcorrelation is similar to the observed one, although slightly flatter.For completeness, the right panel of Figure 1 compares thestellar mass function predicted by the model with the observed one(Bernardi et al. 2013; Bernardi et al. 2017). While the model liesslightly above the data at the highest masses, it lies below overthe range × . M star /M ⊙ . . This is not a majorissue in the present context, since the model is consistent with theempirical galaxy scaling relations, and most notably with the σ - M star relation shown in the left panel.Having checked that our model reproduces the dynamicalscaling relations of early-type galaxies, we now study the scal-ing relations with the central supermassive black hole. Figure 2compares the model (with no morphological selection) with thewhole sample (i.e., spirals as well as ellipticals and lenticulars) ofSavorgnan et al. (2016, blue diamonds). The left and right panelsshow the scaling of black hole mass with velocity dispersion andtotal stellar mass, respectively. In each panel, the solid and dottedlines show the median and the region containing 70% of the modelobjects in a given bin. The top panels show the full black-hole sam-ple, while the middle and bottom panels only show the subset forwhich the sphere of influence exceeds the typical (HST) resolutionlimit, r infl ≡ k GM bh σ , r infl d Ang > . ′′ (1)( d Ang being the angular-diameter distance). We use the parame-ter k to take into account different galaxy mass profiles: k ∼ for the Hernquist profiles assumed by the model (see Barausse2012, for details), but k ∼ (or even larger) is possible ifa core is present. On the other hand, strong lensing and accu-rate dynamical modelling have shown that the mass profiles ofintermediate-mass, early-type galaxies are consistent with nearlyisothermal profiles down to (at least) tenths of the effective radius(e.g., Cappellari et al. 2015, and references therein). These have k ∼ . To bracket these uncertainties, we show results for both k = 10 and 1.To match the data as closely as possible, we draw the an-gular diameter distances d Ang from an empirical probability dis-tribution function, which we construct from the distances in theSavorgnan et al. (2016) sample, and which peaks around ∼ –20 c (cid:13)000 , 1– 11 E. Barausse et al.
Figure 2.
Black-hole mass as a function of velocity dispersion (left panels) and total stellar mass (right panels) as predicted by our model with light seeds (theresults for heavy seeds are qualitatively similar). Results are shown for the full outputs of the model (labelled as “Intrinsic”, top panels), and for the subsampleof black holes with a gravitational sphere of influence above 0.1” (labelled as “Observed”, middle and bottom panels; see text for details). The solid and dottedblue lines in all panels represent the medians and 70% confidence region of the distributions. The dashed magenta lines in the left panels are the medianmodel predictions when assigning velocity dispersions to galaxies from the observed SDSS σ HL − M star relation (long-dashed red line in the left panel ofFigure 1). Long-dashed red lines are the intrinsic scaling relations as inferred from Monte Carlo simulations by Shankar et al. (2016). Blue diamonds are datacollected and updated by Savorgnan et al. (2016) on local galaxies with dynamical measurements of supermassive black holes. Note that observational biasescan increase the normalisation and reduce the scatter in the intrinsic scaling relations. Mpc. Using a distribution that is uniform in comoving volume, asdone in Shankar et al. (2016), yields comparable results.The top panels of Figure 2 suggest that the model predictsintrinsic scaling relations that lie slightly below the data, espe-cially for M bh - M star (top right), but which are broadly consis-tent with the intrinsic relations suggested by Shankar et al. (2016,long dashed red lines). The match to the M bh - σ relation improves(slightly) if we assign velocity dispersions by using the SDSS σ - M star relation, e.g. via the analytic fits provided by Sesana et al.(2016, dashed magenta line, left panels). At the same time, the Note that it makes sense to assign σ from the model-predicted M star viathe observed σ - M star relation, rather than M star from the model-predicted model substantially underpredicts the M bh - M star relation by upto an order of magnitude (top right panel), suggesting that there issome internal inconsistency in the data. In other words, a modellike ours, tuned to match the local velocity dispersion-stellar massand black-hole mass-velocity dispersion relations, tends to severelyunderpredict the M bh - M star relation. This is in line with the re-sults of Shankar et al. (2016). Calibrating the model to match theobserved M bh - M star relation would instead overestimate the ob-served M bh - σ relation. Such an overestimate has indeed been seen σ , because masses are more “primitive” quantities than velocity dispersionsfor a semi-analytic galaxy formation model such as ours.c (cid:13) , 1– 11 MBHs: selection bias and models Figure 3.
Same as Figure 2, but for galaxies with a bulge-to-total ratio
B/T > . (and for light seeds; the results are similar for heavy seeds). The resultsare broadly similar to Figure 2, though the dispersion of the model’s intrinsic predictions is substantially lower than for the full galaxy sample. Note that theobservational dataset is also restricted to early-type galaxies only, to match the sample of simulated galaxies. in two (very different) cosmological hydrodynamic simulations(Sijacki et al. 2015; Volonteri et al. 2016). When the selection effect on the sphere of influence of theblack hole is applied to the model (middle and bottom panels ofFigure 2), the median normalisations of the predicted scaling rela-tions increase (especially for the M bh - M star relation), because asubstantial fraction of low-mass black holes are excluded. This isbecause, for a given angular aperture, Equation 1 preferentially re-moves objects with the smallest gravitational spheres of influence;these tend to be the lowest-mass black holes. Therefore, this effecttends to select the “upper end” (in black-hole mass) of the intrinsicdistributions shown in the top panels of Figure 2. This also induces Volonteri et al. (2016) mention resolution, which is at best 1 kpc in theirsimulation, as one reason why their results do not match the M bh - σ rela-tion. However, their predictions for M bh are larger than the observationseven at large σ (c.f. their Figure 7), while they are in good agreement withthe M bh - M star relation. an overall flattening of the scaling relations, which is again moreobvious in the M bh - M star plane: selection hardly matters for themost massive galaxies, but it causes a factor . increase in themedian observed M bh at lower masses. Selection-biased modelsare flatter in the M bh - σ plane as well. We conclude that, to agree with resolution-biased observedscaling relations, models must predict intrinsic scaling relationsthat are significantly steeper than the observed relations. We havetried to obtain steeper intrinsic scaling relations in our model bychanging the AGN feedback, but we have found this to be insuffi-cient to achieve better agreement with the data at low masses andvelocity dispersions. In fact, the results and conclusions of this pa- As a result of this flatter slope, the model’s prediction (after applying theselection bias) lies above the (few) data with σ ∼ km/s in the sampleof Savorgnan et al. (2016). Note however that other samples, such as that ofKormendy et al. (2011), contain black holes with masses up to ∼ M ⊙ at σ ∼ km/s, which is in better agreement with our model.c (cid:13)000
B/T > . (and for light seeds; the results are similar for heavy seeds). The resultsare broadly similar to Figure 2, though the dispersion of the model’s intrinsic predictions is substantially lower than for the full galaxy sample. Note that theobservational dataset is also restricted to early-type galaxies only, to match the sample of simulated galaxies. in two (very different) cosmological hydrodynamic simulations(Sijacki et al. 2015; Volonteri et al. 2016). When the selection effect on the sphere of influence of theblack hole is applied to the model (middle and bottom panels ofFigure 2), the median normalisations of the predicted scaling rela-tions increase (especially for the M bh - M star relation), because asubstantial fraction of low-mass black holes are excluded. This isbecause, for a given angular aperture, Equation 1 preferentially re-moves objects with the smallest gravitational spheres of influence;these tend to be the lowest-mass black holes. Therefore, this effecttends to select the “upper end” (in black-hole mass) of the intrinsicdistributions shown in the top panels of Figure 2. This also induces Volonteri et al. (2016) mention resolution, which is at best 1 kpc in theirsimulation, as one reason why their results do not match the M bh - σ rela-tion. However, their predictions for M bh are larger than the observationseven at large σ (c.f. their Figure 7), while they are in good agreement withthe M bh - M star relation. an overall flattening of the scaling relations, which is again moreobvious in the M bh - M star plane: selection hardly matters for themost massive galaxies, but it causes a factor . increase in themedian observed M bh at lower masses. Selection-biased modelsare flatter in the M bh - σ plane as well. We conclude that, to agree with resolution-biased observedscaling relations, models must predict intrinsic scaling relationsthat are significantly steeper than the observed relations. We havetried to obtain steeper intrinsic scaling relations in our model bychanging the AGN feedback, but we have found this to be insuffi-cient to achieve better agreement with the data at low masses andvelocity dispersions. In fact, the results and conclusions of this pa- As a result of this flatter slope, the model’s prediction (after applying theselection bias) lies above the (few) data with σ ∼ km/s in the sampleof Savorgnan et al. (2016). Note however that other samples, such as that ofKormendy et al. (2011), contain black holes with masses up to ∼ M ⊙ at σ ∼ km/s, which is in better agreement with our model.c (cid:13)000 , 1– 11 E. Barausse et al.
Figure 4.
Median black-hole mass as a function of total host-galaxy stellarmass (solid lines) as predicted by the model, for the black holes with bolo-metric luminosity log L/ erg s − > . The dotted (long-dashed) linesmark the 90% (99%) confidence regions at given stellar mass, as predictedby our model. Blue squares are data from Reines & Volonteri (2015). Theresults are for the light-seed model, but the heavy-seed one gives very sim-ilar results. per are robust to changes in the AGN feedback strength as well asto changes in the black-hole accretion prescriptions.Figure 2 also clearly shows that in addition to changing thenormalisation and slope, selection dramatically decreases the dis-persion 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 Figure 3, which compares model-predicted galaxies with bulge-to-total ratios of B/T > . withthe E+S0 galaxies of Savorgnan et al. (2016), in the same format asFigure 2. The intrinsic distributions of the model-predicted galaxies(top panels) are narrower than in the full samples (c.f. top panels ofFigure 2), and the median scaling relations have higher normalisa-tions, in better agreement with the data and in line with the findingsof Barausse (2012). Although the selection effect is smaller for thisspecific subsample of galaxies, the model-predicted intrinsic M bh - σ relation is offset slightly above the data, while the predictionsfor the intrinsic M bh - M star relation are slightly below the data.The effect of including the selection bias on the resolvability of theblack hole sphere of influence (middle and bottom panels) is lessimportant than in Figure 3, because the small systems for which thesphere of influence is not resolvable tend to live in late-type galax-ies in our semi-analytic model. Nevertheless, the selection bias stilltends to make the correlations higher in normalisation, slightly flat-ter, and slightly tighter. It is important to emphasise that, whatever the exact black-holeseeding recipe adopted in the model, the predicted local (intrin-sic) scaling relations reported in Figure 2 for all galaxy typesshow very large dispersions, especially for galaxy stellar masses M star . × M ⊙ , with a broad distribution up to two or-ders of magnitude or more in black hole mass at fixed stellar massor velocity dispersion. This is because the model tends to retaina significant number of low-mass black holes, which did not ac-crete much gas along their history (because they live in spirals or satellite galaxies) and which therefore remain closer to their seedmasses (e.g., Volonteri et al. 2005; Barausse 2012). However, formore massive galaxies the dispersion is smaller, with almost nogalaxies with M star & × M ⊙ having black holes with M bh . M ⊙ .One interesting way to probe the existence of very low-massblack holes in relatively low-mass galaxies could be to comparewith the scaling relations in active galaxies, which are not lim-ited by the spatial resolution issues that heavily affect dynamicalmeasurements of black holes. To this purpose, Figure 4 comparesthe 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 linesrepresent the median, the 90% confidence region (i.e. the 5th and95th percentiles) and the 99% confidence region (i.e. the . th and . th) of the (model-predicted) M bh - M star relation, by assum-ing a light black-hole seed scenario (the heavy-seed scenario givesvery similar results) and considering only systems with bolomet-ric luminosity log( L/ erg s − ) > (roughly the minimum lu-minosity probed by Reines & Volonteri 2015). The model’s distri-bution of active black holes has been built by randomly drawingEddington ratios from a Schechter distribution that extends up tothe 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 thatour predicted luminosity function at z = 0 , computed by assumingan average duty cycle of active black holes of 10% consistent withthe 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 luminosityfunctions of Hopkins et al. (2007) and Shankar et al. (2009).First, let us note, as emphasised by Shankar et al. (2016), thata lower limit of L & erg s − should still allow black holesdown to a mass of M bh ∼ M ⊙ to be detected, at least if anon-negligible fraction of these black holes are still accreting at theEddington limit. Such low mass black holes do not seem to exist inthe Reines & Volonteri (2015, see also Baldassare et al. 2015) sam-ple (and in our model, at least in sufficiently large numbers and withhigh enough Eddington ratios to warrant detection). Even assuminglower virial factors f 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 et al.2016, and references therein), would not alter these conclusions.Second, it is clear that the observational sample lies, on av-erage, below the model median predictions, which were tuned toreproduce the data on inactive local galaxies with a significantlyhigher normalisation. While the predicted 90% and 99% confidenceregions for the active black-hole population encompass the data ofReines & Volonteri (2015), the model also predicts the existenceof a large number of active higher-mass black holes (above themedian), which are not observed. Therefore, either the sample ofReines & Volonteri (2015) is biased toward low-luminosity activesystems, or our model should be normalised to lower values by afactor & .The large scatter in our model means that, if the normalisationis decreased, then our model would predict a large tail of very low-mass black holes. Since these are not observed, this would haveimportant consequences for constraining models of the seeds ofthe supermassive black-hole population. On the other hand, de-creasing the normalisation of the M bh - M star relation is by nomeans straightforward. While it could be achieved by decreasingblack-hole accretion, this would imply a proportional reduction in c (cid:13) , 1– 11 MBHs: selection bias and models AGN luminosity, unless a higher radiative efficiency and/or dutycycle are also assumed. The present calibration of our model al-ready predicts rather high radiative efficiencies/black-hole spins(Sesana et al. 2014). Duty cycles of active black holes could in-stead be constrained by comparing with independent AGN clus-tering measurements (Gatti et al. 2016, and references therein). Weplan to explore some of these important interrelated issues in futurework.
We now compare our model’s predictions for the residuals of theblack hole-galaxy scaling relations to the data. Such correlationsare an efficient way of going beyond pairwise correlations betweenthe variables themselves (Bernardi et al. 2005; Sheth & Bernardi2012). For example, measurements of the M bh - σ and M bh - M star correlations alone do not provide insight about whether σ is moreimportant than M star in determining M bh . This is because the M bh - σ and M bh - M star correlations do not encode complete infor-mation about the joint distribution of M bh , σ and M star . Correla-tions between the residuals encode some of this extra information.To this purpose, the left and right hand panels of Figure 5 show ∆( M bh | M star ) vs ∆( σ | M star ) and ∆( M bh | σ ) vs ∆( M star | σ ) ,where ∆( Y | X ) ≡ log Y − h log Y | log X i (2)is the residual in the Y variable (at fixed X ) from the log-log-linearfit of Y ( X ) vs X , i.e. h log Y | log X i . The magenta long-dashedand dotted lines in each panel show the best fit and 1 σ uncertain-ties on the correlations between residuals in the Savorgnan et al.(2016) dataset: red circles and green triangles represent ellipticalsand lenticulars. We obtained the magenta lines by running 200 iter-ations following the steps outlined in Shankar et al. (2017), whichinclude errors in both variables. At each iteration we eliminate threerandom objects from the original sample. From the full ensembleof realizations, we measure the mean slope and its 1 σ uncertainty.The blue solid and dotted lines show a similar analysis in oursemi-analytic model. However, in this case, we randomly produce30 mock samples of ∼ galaxies each, with B/T > . andresolvable black-hole spheres of influence. From the full ensembleof mock realizations, we then extract the mean slope and Pearsoncoefficient.The correlations in Figure 5 show that, in the data, the veloc-ity dispersion is more strongly correlated with the black-hole mass,with a mean Pearson coefficient of r = 0 . , than stellar mass, forwhich the Pearson coefficient is r = 0 . . The model instead pre-dicts just the opposite, with almost zero correlation with velocitydispersion (mean r = 0 . ), but with a correlation with stellar massconsistent with the data, though still rather weak (mean r = 0 . ).It is nevertheless important to realise that even an intrinsically weakcorrelation with velocity dispersion at fixed stellar mass does notnecessarily imply that the total correlation with velocity dispersionis small. In fact, following Appendix B in Shankar et al. (2017), thetotal dependence of the black-hole mass on velocity dispersion canbe summarised as M bh ∝ σ β M α star ∝ σ β + α γ , where γ comesfrom M star ∝ σ γ (where we have ignored any other explicit de-pendence on, e.g., S´ersic index or galaxy size). Since the modelpredicts γ ≈ , and the residuals in Figure 5 (solid blue lines)yield β ∼ . and α ∼ . , one obtains a total dependence of M bh ∝ σ , consistent with Figure 3 (middle left panel). Also notethat a slope γ ≈ for the M star - σ relation is in itself already intension with the ( V max -corrected) SDSS observations, which rather prefer a slope γ ≈ . –2.5, at least for M star & × M ⊙ (c.f.also Figure 7 below).For completeness, Figure 5 also shows results from theMonte Carlo simulations performed by Shankar et al. (2016, greybands), which assume an intrinsic correlation of the type M bh ∝ σ . M . and include the selection bias on r infl . These models re-produce the observed trends.Several comments are in order. First, our results are robustagainst the efficiency of the AGN feedback. Indeed, in our refer-ence model we fixed the AGN feedback efficiency to obtain a stellarmass function as close as possible to the observations at the high-mass end (c.f. Figure 1, right panel). By slightly decreasing thisefficiency (e.g. by a factor ∼ ), the agreement with the observedstellar mass function in Figure 1 slightly worsens in the high-massend, while the M bh - σ relation steepens ( M bh ∝ σ . after includ-ing the effect of the selection bias), but not enough to fully matchthe observations (for which slopes ∼ . − are usually quotedin the literature, e.g. Graham et al. 2011; Kormendy & Ho 2013).More importantly, the correlations of the residuals remain almostunchanged ( r = 0 . when fixing M star , and r = 0 . when fix-ing σ ). Second, we have checked that these correlations are notimproved by considering a different choice of k in Equation 1;by assuming no bias on the resolvability of the black-hole sphereof influence; by considering all galaxies rather than just bulge-dominated ones; by using the bulge mass instead of the total stellarmass; or by using the the SDSS σ - M star relation (e.g., by the fitsof Sesana et al. 2016) to compute velocity dispersions. Third, themodel’s residuals shown in Figure 5 do not account for measure-ment errors (i.e. we compute the residuals for the model’s exactpredictions, without folding in any observational uncertainties). Wehave verified that including these errors can yield a steeper slope β ∼ for the residuals at fixed M star , but this does not strengthenthe correlations of the residuals shown above (essentially becausethe error on β also grows when β grows).To check if our results are due to the particular implementa-tion of black-hole accretion and AGN feedback in our semi-analyticmodel, and/or lack of sufficient “coupling” between velocity dis-persion and AGN feedback, we have performed a similar anal-ysis of the Horizon-AGN simulation (Dubois et al. 2014, 2016).This is a hydrodynamic cosmological simulation of a box with size Mpc /h , run with the adaptive mesh refinement code RAM-SES (Teyssier 2002), with dark-matter particles and a min-imum 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 andblack holes in Horizon-AGN). Galaxies are extracted with a galaxyfinder running on star particles. Since the simulation assumesa Salpeter (Salpeter 1955) IMF, we have reduced galaxy stel-lar masses by 0.25 dex (e.g., Bernardi et al. 2010) to match theChabrier IMF (Chabrier 2003) adopted in this work. The velocitydispersion is measured from its components along each directionof the cylindrical coordinates oriented along the galaxy’s spin axis(i.e. σ r , σ t and σ z for the radial-, tangential-, and z -component re-spectively), thus σ = ( σ r + σ t + σ z ) / . The velocity dispersionof each galaxy is measured using only star particles within the ef- Note that, in our model, increasing the AGN feedback efficiency de-creases the slope of the M bh - σ relation, mainly because AGN feedbackis more effective at clearing the gas in more massive galactic hosts, thusinhibiting the growth of especially the more massive black holes.c (cid:13) , 1– 11 E. Barausse et al.
Figure 5.
Correlation between the residuals of the M bh - M star relations and those of the σ - M star relation, at fixed stellar mass (left panel), and betweenthe residuals of the M bh - σ relation and those of the M star - σ relation, at fixed velocity dispersion (right panel). Red circles and green triangles show theellipticals 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 predictionsand the observations, respectively), while the dotted lines show the 1 σ uncertainty on the slope. Also reported are the best fits and the Pearson correlationcoefficient r . The grey bands show the residuals extracted from the Monte Carlo simulations from Shankar et al. (2016), with selection bias on the black-holegravitational sphere of influence. Note that the correlation of the residuals at fixed stellar mass (left panel) is very strong in the data, but essentially absent inthe model. These results are for the light-seed model (the heavy-seed one yields similar results), selecting only early-type galaxies ( B/T > . ) for whichthe black-hole sphere of influence is resolvable. Figure 6.
Same as Figure 5, but for the hydrodynamic cosmological simulation Horizon-AGN (Dubois et al. 2014, 2016). However, unlike in Figure 5, noselection bias on the resolvability of the black-hole sphere of influence has been applied. Doing so leads to slightly weaker correlations and slightly lowerslopes 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 notcorrected to the Hyperleda aperture. fective radius of the galaxy. The results are shown in Figures 6 and7. They are in qualitative agreement with those obtained with oursemi-analytic model. In more detail, Horizon-AGN also predicts γ ≈ , in tension with the data, and a very weak correlation be-tween the residuals when fixing M star . Also note that in Figure 6we have not applied any selection bias (unlike in the case of thesemi-analytic model in Figure 5). Restricting to systems with re-solvable spheres of influence actually makes the correlations in theresiduals of the Horizon-AGN simulation even weaker.We plan to consider different (and possibly stronger) modelsof AGN feedback in future work. For now, the discrepancies high-lighted in Figures 5 and 6 clearly show that correlations betweenthe residuals of the M bh - σ and M bh - M star scaling relations aremore powerful than the scaling relations themselves at constrainingmodels for the co-evolution of black holes and their host galaxies. Concerning the normalisation of the black-hole scaling relations, atface value our model clearly fails to reproduce the observed M bh - σ and M bh - M star relations at the same time. Without invoking anyselection effect, one could attempt to improve the simultaneousmatch to both the observed scaling relations by fine-tuning someof the key parameters in the model, namely the one controlling gasaccretion onto the black hole and/or the energetic feedback fromthe central active nucleus. For example, increasing the efficiencyof AGN feedback could in principle decrease the stellar masses ofthe host galaxies at fixed black hole mass, thus possibly improv-ing the match to the M bh - M star relation (top right panel of Fig-ure 2). However, this would then spoil the good match to the veloc-ity dispersion-stellar mass relation. Similarly, one could increase c (cid:13) , 1– 11 MBHs: selection bias and models Figure 7.
Total stellar mass as a function of velocity dispersion for SDSSgalaxies with P ( E − S > . (shaded region), and for the Horizon-AGN simulation (the solid line represents the mean, and the dotted linesrepresent the 70% confidence region; Dubois et al. 2014, 2016). As in theobservational sample, we select early-type galaxies from the simulationsby only considering systems with V c /σ < . , where V c is the rotationalvelocity. Applying a different cut does not change this figure significantly.Here, for both the SDSS and Horizon-AGN data, σ e is computed within theeffective radius. Given the resolution of Horizon-AGN, this quantity is morereliably estimated than the more central velocity dispersion σ HL adopted inthe previous figures. the accretion onto the central black hole at fixed star formation rateand final stellar mass of the host galaxy, with the aim of improvingthe match to the M bh - M star relation. This however would propor-tionally increase the M bh - σ relation above the data, when velocitydispersions are chosen to faithfully track those from the observedSDSS σ − M star relation (magenta lines in the left panels of Fig-ure 2). Concerning the dispersion around the scaling relations, oursemi-analytic model, which self-consistently evolves black holesfrom seeds at high redshifts, naturally predicts very broad distribu-tions in black-hole mass at fixed velocity dispersion or stellar mass,at relatively low masses M star . × M ⊙ .Overall, these effects are in line with the results of the MonteCarlo simulations presented in Shankar et al. (2016), though the in-trinsic scatters assumed there, especially in the M bh - σ relation,were always below . dex. This was needed to avoid too flatslopes in the “observed”, biased relations. Indeed, here our pre-dicted slope of the M bh - σ relation, inclusive of observational bias(e.g., middle and bottom left panels of Figure 2), is approximately M bh ∝ σ β with slope β ∼ − . , significantly flatter thanthe slopes β & . − usually quoted in the literature (e.g.,Graham et al. 2011; Kormendy & Ho 2013).Conversely, a key point of our present work is that, irrespec-tive of the chosen masses for the seed black holes, the modelpredicts relatively tight scaling relations at high stellar masses M star & × M ⊙ . In this respect, our model does not supportthe conjecture put forward by Batcheldor et al. (2007), accordingto which the M bh - σ relation is only an upper limit of a more orless uniform distribution of black holes. This is also in line withthe Monte Carlo tests performed by Shankar et al. (2016, their Fig-ure 11), in which a very broad distribution of black holes extendingto the lowest masses in all types of galaxies is highly disfavoured.Nevertheless, Figure 4 shows that the present data on active galax-ies may still be consistent with our model, once the proper fluxlimits and Eddington ratio distributions are accounted for. There-fore, an intrinsically broad distribution for the black-hole mass in relatively small galaxies with mass M star . × M ⊙ , extend-ing down to black holes of mass M bh & M ⊙ or so, cannot beexcluded by present data, though it is highly disfavoured in moremassive galaxies. Models with such broad distributions, however,still tend to produce scaling relations flatter than observed, once thebias on the black-hole sphere of influence and on galaxy morphol-ogy is folded in the model predictions, as shown in Figures 2 and3. Finally, we stress again that our semi-analytic model failsto reproduce the observations when it comes to correlations be-tween the residuals of the scaling relations. In more detail, themodel is consistent with the data for the correlation between theresiduals in the M bh - σ relation and those in the M star - σ rela-tion, at fixed velocity dispersion. However, the model predicts al-most no correlation between the residuals of the M bh - M star re-lations and those of the σ - M star relation, at fixed stellar mass,while the data hint at a rather strong correlation. We have veri-fied that these results are robust against changing the parametersof the model (namely the AGN feedback efficiency and the pa-rameter regulating black-hole accretion). Moreover, we have shownthat the same weak correlation between the residuals at fixed stel-lar mass is also obtained in the hydrodynamic cosmological sim-ulation Horizon-AGN (Dubois et al. 2014, 2016), which includesthermal quasar-mode feedback and jet-structured radio-mode AGNfeedback (i.e., processes that are expected to induce a stronger cou-pling between black-hole mass and velocity dispersion). Anothernoteworthy point is that both our model and Horizon-AGN tendto produce slopes for the σ - M star relation that are significantlysteeper than the data, although in Horizon-AGN numerical reso-lution effects may bias the measurements of the velocity dispersionin lower-mass galaxies (see discussion in Dubois et al. 2016).Our interpretation is that the weak correlation of the residu-als at fixed stellar mass may be at least partly due to current AGNfeedback models being possibly too weak to capture the full effectof the black hole on the stellar velocity dispersion. While this isexpected in a semi-analytic model such as ours, where the AGNfeedback is typically assumed to simply eject gas from the nu-clear region (Granato et al. 2004; Barausse 2012), it is more sur-prising in the Horizon-AGN simulation, which captures the back-reaction of the ejected gas onto the stellar and dark-matter dynam-ics (Peirani et al. 2016). Nonetheless, in hydrodynamic cosmolog-ical simulations such as Horizon-AGN, the spatial resolution is 1kpc at best and may limit our capability to properly capture theinteraction of AGN winds with gas and its impact on the dynam-ics within galactic bulges. Moreover, the Horizon-AGN simulationcurrently employs a rather crude model of AGN feedback, mostlybased on simple gas heating, to mimic the so-called “quasar-mode”feedback. More realistic models of AGN feedback, also inclusive ofmomentum-driven winds and stronger coupling with the surround-ing interstellar medium (e.g., Bieri et al. 2017), might possibly bemore effective at improving the comparison with the black hole-galaxy scaling relations and their residuals. We should also mentionthat another possibility that has been put forward in the literatureis that AGN feedback may not be the main cause of the scalingrelations, which might be ascribed instead to a common gas sup-ply for the galaxy and the black hole, regulated by gravitationaltorques (Angl´es-Alc´azar et al. 2017). c (cid:13) , 1– 11 E. Barausse et al.
We have compared the predictions of a comprehensive semi-analytic model of galaxy formation, which self-consistentlyevolves supermassive black holes from high-redshift seeds by ac-counting for gas accretion, mergers and AGN feedback, with theobserved scaling relations of the masses of supermassive blackholes with stellar mass and velocity dispersion. Our main conclu-sions are:(i) At M star & × M ⊙ , the dispersion in black-hole massat fixed stellar mass is . dex – very few black holes with masses M bh . M ⊙ are predicted in such massive galaxies. However,for galaxies having M star . × M ⊙ , the distribution of M bh at fixed σ or M star is broad.(ii) Observational selection effects, associated with resolvingthe black-hole sphere of influence and/or with selecting bulge-dominated/elliptical galaxies, tighten the M bh - σ and M bh - M star scaling relations, bringing them into better agreement with the ob-servations.(iii) No evident variation in AGN feedback and/or black-holeaccretion efficiencies can provide a simultaneous match to bothscaling relations. This supports previous work suggesting an in-ternal inconsistency between the observed M bh - σ and M bh - M star relations.(iv) Galaxy-evolution models (our semi-analytic one as wellas the Horizon-AGN simulation) predict almost no correlationbetween the residuals of the M bh - M star relations and those ofthe σ - M star relation, at fixed stellar mass. Since the data hintat a rather strong correlation, this calls for revamped AGN feed-back recipes in the next generation of cosmological galaxy-evolution models, or for a re-assessment of the importance ofgravitational torques in regulating the black hole-galaxy co-evolution (Angl´es-Alc´azar et al. 2017). ACKNOWLEDGMENTS
We thank M. Volonteri for insightful conversations. EB ac-knowledges support from the H2020-MSCA-RISE-2015 GrantNo. StronGrHEP-690904 and from the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche. Thiswork has made use of the Horizon Cluster, hosted by the Institutd’Astrophysique de Paris. We thank Stephane Rouberol for runningthis cluster so smoothly.
REFERENCES
Aird J. et al., 2012, ApJ, 746, 90Angl´es-Alc´azar, D., Dav´e, R., Faucher-Gigu`ere, C.-A., ¨Ozel, F.,& Hopkins, P. F. 2017, MNRAS, 464, 2840Antonini F., Barausse E., Silk J., 2015a, ApJ, 812, 72Antonini F., Barausse E., Silk J., 2015b, ApJ, 806, L8Baes M., Dejonghe H., 2002, A&A, 393, 485Baldassare V. F., Reines A. E., Gallo E., Greene J. E., 2015, ApJ,809, L14Barausse E., 2012, MNRAS, 423, 2533Batcheldor D., Marconi A., Merritt D., Axon D. J., 2007, ApJ,663, L85Bernardi, M., Meert, A., Sheth, R. K., et al. 2017, MNRAS, inpress, arXiv:1604.01036 Bernardi M., Meert A., Sheth R. K., Vikram V., Huertas-CompanyM., Mei S., Shankar F., 2013, MNRAS, 436, 697Bernardi M., Shankar F., Hyde J. B., Mei S., Marulli F., ShethR. K., 2010, MNRAS, 404, 2087Bernardi M., Sheth R. K., Nichol R. C., Schneider D. P.,Brinkmann J., 2005, AJ, 129, 61Bernardi M., Sheth R. K., Tundo E., Hyde J. B., 2007, ApJ, 660,267Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854Bongiorno A. et al., 2012, MNRAS, 427, 3103Bournaud F. et al., 2011, ApJ, 730, 4Bower R., Schaye J., Frenk C. S., Theuns T., Schaller M., CrainR. A., McAlpine S., 2016, ArXiv:160707445Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S.,Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645Cappellari M. et al., 2006, MNRAS, 366, 1126Cappellari M. et al., 2015, ApJ, 804, L21Chabrier G., 2003, PASP, 115, 763Cirasuolo M., Shankar F., Granato G. L., De Zotti G., Danese L.,2005, ApJ, 629, 816Croton D. J. et al., 2006, MNRAS, 365, 11De Lucia G., Blaizot J., 2007, MNRAS, 375, 2Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C.,Lopez J., Springel V., 2012, ApJ, 745, L29Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS,420, 2662Dubois, Y., Pichon, C., Haehnelt, M., et al. 2012, MNRAS, 423,3616Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS,433, 3297Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444,1453Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463,3948Ferrarese L., Ford H., 2005, Space Science Reviews, 116, 523Fontanot F., Monaco P., Shankar F., 2015, MNRAS, 453, 4112Gatti, M., Shankar, F., Bouillot, V., et al. 2016, MNRAS, 456,1073Gonz´alez J. E., Lacey C. G., Baugh C. M., Frenk C. S., 2011,MNRAS, 413, 749Goulding A., Alexander D., 2009, ArXiv:0906.0772Graham A. W., 2016, Galactic Bulges, 418, 263Graham A. W., Driver S. P., Allen P. D., Liske J., 2007, MNRAS,378, 198Graham A. W., Onken C. A., Athanassoula E., Combes F., 2011,MNRAS, 412, 2211Granato G. L., De Zotti G., Silva L., Bressan A., Danese L., 2004,ApJ, 600, 580G¨ultekin, K., Tremaine, S., Loeb, A., & Richstone, D. O. 2011,ApJ, 738, 17Guo Q. et al., 2011, MNRAS, 413, 101Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., RobertsonB., Springel V., 2006, ApJS, 163, 1Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731Huertas-Company M., Aguerri J. A. L., Bernardi M., Mei S.,S´anchez Almeida J., 2011, A&A, 525, A157Jones M. L., Hickox R. C., Black C. S., Hainline K. N., DiPompeoM. A., Goulding A. D., 2016, ApJ, 826, 12Kauffmann G., Heckman T. M., 2009, MNRAS, 397, 135Klein A. et al., 2016, Phys. Rev. D, 93, 024003Kormendy, J., Bender, R., & Cornell, M. E. 2011, Nature, 469,374 c (cid:13) , 1– 11 MBHs: selection bias and models Kormendy J., Ho L. C., 2013, ARA&A, 51, 511Lapi A., Shankar F., Mao J., Granato G. L., Silva L., De Zotti G.,Danese L., 2006, ApJ, 650, 42Malbon R. K., Baugh C. M., Frenk C. S., Lacey C. G., 2007, MN-RAS, 382, 1394McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184Meert, A., Vikram, V., & Bernardi, M. 2015, MNRAS, 446, 3943Monaco P., Fontanot F., Taffoni G., 2007, MNRAS, 375, 1189Morabito, L. K., & Dai, X. 2012, ApJ, 757, 172Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A.2010, ApJ, 725, 2312Oser, L., Naab, T., Ostriker, J. P., & Johansson, P. H. 2012, ApJ,744, 63Pardo K. et al., 2016, ArXiv:160301622Parkinson H., Cole S., Helly J., 2008, MNRAS, 383, 557Paturel G., Petit C., Prugniel P., Theureau G., Rousseau J., BroutyM., Dubois P., Cambr´esy L., 2003, A&A, 412, 45Peirani, S., Dubois, Y., Volonteri, M., et al. 2016,arXiv:1611.09922Press W. H., Schechter P., 1974, ApJ, 187, 425Reines A. E., Volonteri M., 2015, ApJ, 813, 82Rodriguez-Gomez V. et al., 2016, ArXiv e-printsSalpeter, E. E. 1955, ApJ, 121, 161Savorgnan G. A. D., 2016, ApJ, 821, 88Savorgnan G. A. D., Graham A. W., Marconi A., Sani E., 2016,ApJ, 817, 21Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS,447, 2085Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794,104Sesana A., Shankar F., Bernardi M., Sheth R. K., 2016, MNRAS,463, L6Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006,ApJ, 643, 14Shankar F., 2009, New Astron. Rev., 53, 57Shankar F., Weinberg D. H., Miralda-Escud´e J., 2009, ApJ, 690,20Shankar, F., Marulli, F., Bernardi, M., et al. 2010, MNRAS, 403,117Shankar F., Marulli F., Bernardi M., Mei S., Meert A., Vikram V.,2013, MNRAS, 428, 109Shankar F. et al., 2016, MNRAS, 460, 3119Shankar, F., Bernardi, M., & Sheth, R. K. 2017, MNRAS, in press,arXiv:1701.01732Sheth R. K., Bernardi M., 2012, MNRAS, 422, 1825Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Sny-der G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575Teyssier, R. 2002, A&A, 385, 337Tundo E., Bernardi M., Hyde J. B., Sheth R. K., Pizzella A., 2007,ApJ, 663, 53van den Bosch R. C. E., Gebhardt K., G¨ultekin K., Yıldırım A.,Walsh J. L., 2015, ApJS, 218, 10Vittorini V., Shankar F., Cavaliere A., 2005, MNRAS, 363, 1376Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620,69Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MN-RAS, 460, 2979Welker, C., Dubois, Y., Devriendt, J., et al. 2017, MNRAS, 465,1241Yong S. Y., Webster R. L., King A. L., 2016, ArXiv:1602.04672 c (cid:13)000