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

TL;DR

This study uses semi-analytic galaxy formation models to investigate biases in supermassive black hole mass estimates, revealing inconsistencies and the importance of selection effects in matching observed scaling relations.

Contribution

It demonstrates that current models struggle to simultaneously reproduce black hole scaling relations and residual correlations, highlighting potential issues in AGN feedback modeling.

Findings

Selection effects align model predictions with observations.

Weak residual correlations challenge current AGN feedback models.

Models show internal inconsistencies in black hole-galaxy relations.

Abstract

Recent work has confirmed that the masses of supermassive black holes, estimated from scaling relations with global properties such as the stellar masses of their host galaxies, may be biased high. Much of this may be caused by the requirement that the gravitational sphere of influence of the black hole must be resolved for the black-hole mass to be reliably estimated. We revisit this issue by using a comprehensive galaxy evolution semi-analytic model, which self-consistently evolves supermassive black holes from high-redshift seeds via gas accretion and mergers, and also includes AGN feedback. Once tuned to reproduce the (mean) correlation of black-hole mass with velocity dispersion, the model is unable to also account for the correlation with stellar mass. This behaviour 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 local scaling relations are more effective than the scaling relations themselves at constraining AGN feedback models. In fact, we find that our semi-analytic 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 $z\sim 0$ outputs of the hydrodynamic cosmological simulation Horizon-AGN, which includes state-of-the-art treatments of AGN feedback. This suggests that current AGN feedback models may be too weak or are simply not capturing the effect of the black hole on the stellar velocity dispersion.

Equations

$$r_{\mathrm{infl}}\equiv k\frac{G{M}_{\mathrm{bh}}}{{\sigma }^2},\frac{r_{\mathrm{infl}}}{d_{\mathrm{Ang}}}>0.1^{\prime \prime }$$

$$\Delta (Y|X)\equiv \log Y-⟨\log Y|\log X⟩$$

Full text

††pagerange: Selection bias in dynamically-measured supermassive black hole samples: Scaling relations and correlations between residuals in semi-analytic galaxy formation models– References††pubyear: 2017

1 Introduction

Supermassive black holes, with masses $M_{\rm 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; Graham, 2016; McConnell & Ma, 2013, for reviews). A surprising finding that has puzzled astrophysicists for the last forty years 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 a thousand 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 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; Barausse, 2012; Dubois et al., 2012a, 2016; Fontanot et al., 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 $\sim 80\%$ 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ález et al., 2011; Shankar et al., 2013; 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; Di Matteo et al., 2012; Barausse, 2012; Dubois et al., 2012b).

Recently, Shankar et al. (2016, see also , , and ) showed that the local sample of galaxies with dynamical mass measurements of supermassive black holes is biased. Local galaxies hosting supermassive black holes, irrespective of their exact 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 sample 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 gravitational 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 dispersion 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 normalisation), 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 et al., 2009). Shankar et al. (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 the present 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 context of a comprehensive semi-analytic model of galaxy formation and evolution (Barausse, 2012). This model also evolves supermassive black holes self-consistently from high-redshift “seeds”, and accounts for black-hole mergers and for the feedback from active galactic nuclei (AGNs). After briefly reviewing the model in Section 2, we discuss (Sections 3.1 and 3.2) the slope, normalisation 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 summarise our conclusions in Section 5.

2 Model

The full description of the semi-analytic model adopted as a reference 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 et al. (2015a, b). Here, we briefly summarise 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 generated via extended Press-Schechter algorithms (e.g., Press & Schechter, 1974; Parkinson et al., 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 initially seeded with either light black holes of $M_{\rm seed}\sim 200\,M_{\odot}$ (to be interpreted e.g. as the remnants of PopIII stars), or with heavy black holes of mass $M_{\rm 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_{\rm 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 observations, focusing on the normalisation of the scaling relations and 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 normalisation of the scaling relations and the role of selection bias

The left panel of Figure 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_{\rm 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_{\rm star}$ are from Bernardi et al. (2013). They are the product of luminosity $L$ and mass-to-light ratio $M_{\rm star}/L$; the $L$ values are from Meert et al. (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_{\rm b}/{r_{\rm b}}}[1+(V_{\rm b}/\sigma)^{2}]$, where $M_{\rm b}$ is the bulge dynamical mass, $r_{\rm 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, Figure 2, lower panel), and the ratio $V_{\rm b}/\sigma$ accounts for the contribution of the bulge rotation and is modeled 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 panel of Figure 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 $2\times 10^{10}\lesssim M_{\rm 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_{\rm star}$ relation shown in the left 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. Figure 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 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% 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,

[TABLE]

($d_{\rm 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_{\rm 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 Figure 2 suggest that the model predicts intrinsic scaling relations that lie slightly below the data, especially for $M_{\rm bh}$-$M_{\rm star}$ (top right), but which are broadly consistent with the intrinsic relations suggested by Shankar et al. (2016, long dashed red lines). The match to the $M_{\rm bh}$-$\sigma$ relation improves (slightly) if we assign velocity dispersions by using the SDSS $\sigma$-$M_{\rm star}$ relation, e.g. via the analytic fits provided by Sesana et al. (2016, dashed magenta line, left panels).111Note that it makes sense to assign $\sigma$ from the model-predicted $M_{\rm star}$ via the observed $\sigma$-$M_{\rm star}$ relation, rather than $M_{\rm 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. At the same time, the model substantially underpredicts the $M_{\rm bh}$-$M_{\rm 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_{\rm bh}$-$M_{\rm star}$ relation. This is in line with the results of Shankar et al. (2016). Calibrating the model to match the observed $M_{\rm bh}$-$M_{\rm star}$ relation would instead overestimate the observed $M_{\rm 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).222Volonteri 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_{\rm bh}$-$\sigma$ relation. However, their predictions for $M_{\rm bh}$ are larger than the observations even at large $\sigma$ (c.f. their Figure 7), while they are in good agreement with the $M_{\rm bh}$-$M_{\rm star}$ relation.

When the selection effect on the sphere of influence of the black hole is applied to the model (middle and bottom panels of Figure 2), the median normalisations of the predicted scaling relations increase (especially for the $M_{\rm bh}$-$M_{\rm 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 Figure 2. This also induces an overall flattening of the scaling relations, which is again more obvious in the $M_{\rm bh}$-$M_{\rm star}$ plane: selection hardly matters for the most massive galaxies, but it causes a factor $\lesssim 10$ increase in the median observed $M_{\rm bh}$ at lower masses. Selection-biased models are flatter in the $M_{\rm bh}$-$\sigma$ plane as well.333As 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 in the sample of Savorgnan et al. (2016). Note however that other samples, such as that of Kormendy et al. (2011), contain black holes with masses up to $\sim 10^{8}M_{\odot}$ at $\sigma\sim 10^{2}$ km/s, which is in better agreement with our model.

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.

Figure 2 also clearly shows that in addition to changing the normalisation 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 Figure 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 Figure 2. The intrinsic distributions of the model-predicted galaxies (top panels) are narrower than in the full samples (c.f. top panels of Figure 2), and the median scaling relations have higher normalisations, 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_{\rm bh}$-$\sigma$ relation is offset slightly above the data, while the predictions for the intrinsic $M_{\rm bh}$-$M_{\rm 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 Figure 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 normalisation, slightly flatter, and slightly tighter.

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

It is important to emphasise that, whatever the exact black-hole seeding recipe adopted in the model, the predicted local (intrinsic) scaling relations reported in Figure 2 for all galaxy types show very large dispersions, especially for galaxy stellar masses $M_{\rm 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_{\rm star}\gtrsim 5\times 10^{11}\,M_{\odot}$ having black holes with $M_{\rm 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, Figure 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% confidence region (i.e. the 5th and 95th percentiles) and the 99% confidence region (i.e. the $0.5$th and $99.5$th) of the (model-predicted) $M_{\rm bh}$-$M_{\rm 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/{\rm 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% 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 et al. (2007) and Shankar et al. (2009).

First, let us note, as emphasised by Shankar et al. (2016), that a lower limit of $L\gtrsim 10^{42}\,{\rm erg\,s^{-1}}$ should still allow black holes down to a mass of $M_{\rm 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 ) 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_{\rm 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 average, below the model median predictions, which were tuned to reproduce the data on inactive local galaxies with a significantly higher normalisation. While the predicted 90% and 99% 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 toward low-luminosity active systems, or our model should be normalised to lower values by a factor $\gtrsim 3$.

The large scatter in our model means that, if the normalisation 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 normalisation of the $M_{\rm bh}$-$M_{\rm star}$ relation is by no means straightforward. While it could be achieved by decreasing black-hole accretion, this would imply a proportional reduction in 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 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_{\rm bh}$-$\sigma$ and $M_{\rm bh}$-$M_{\rm star}$ correlations alone do not provide insight about whether $\sigma$ is more important than $M_{\rm star}$ in determining $M_{\rm bh}$. This is because the $M_{\rm bh}$-$\sigma$ and $M_{\rm bh}$-$M_{\rm star}$ correlations do not encode complete information about the joint distribution of $M_{\rm bh}$, $\sigma$ and $M_{\rm star}$. Correlations between the residuals encode some of this extra information.

To this purpose, the left and right hand panels of Figure 5 show $\Delta(M_{\rm bh}|M_{\rm star})$ vs $\Delta(\sigma|M_{\rm star})$ and $\Delta(M_{\rm bh}|\sigma)$ vs $\Delta(M_{\rm star}|\sigma)$, where

[TABLE]

is the residual in the $Y$ variable (at fixed $X$) from the log-log-linear fit of $Y(X)$ vs $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) dataset: 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 $\sim 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 Figure 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 realise 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 summarised as $M_{\rm bh}\propto\sigma^{\beta}M_{\rm star}^{\alpha}\propto\sigma^{\beta+\alpha\,\gamma}$, where $\gamma$ comes from $M_{\rm 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 Figure 5 (solid blue lines) yield $\beta\sim 0.5$ and $\alpha\sim 0.6$, one obtains a total dependence of $M_{\rm bh}\propto\sigma^{3}$, consistent with Figure 3 (middle left panel). Also note that a slope $\gamma\approx 4$ for the $M_{\rm 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_{\rm star}\gtrsim 2\times 10^{10}M_{\odot}$ (c.f. also Figure 7 below).

For completeness, Figure 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_{\rm bh}\propto\sigma^{4.5}M_{\rm star}^{0.5}$ and include the selection bias on $r_{\rm 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. Figure 1, right panel). By slightly decreasing this efficiency (e.g. by a factor $\sim 3$), the agreement with the observed stellar mass function in Figure 1 slightly worsens in the high-mass end, while the $M_{\rm bh}$-$\sigma$ relation steepens444Note that, in our model, increasing the AGN feedback efficiency decreases the slope of the $M_{\rm 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. ($M_{\rm bh}\propto\sigma^{3.4}$ after including the effect of the selection bias), but not enough to fully match the observations (for which slopes $\sim 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_{\rm star}$, and $r=0.3$ when fixing $\sigma$). Second, we have checked that these correlations are not improved by considering a 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 the SDSS $\sigma$-$M_{\rm star}$ relation (e.g., by the fits of Sesana et al., 2016) to compute velocity dispersions. Third, the model’s residuals shown in Figure 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_{\rm 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 dispersion 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$, 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 (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 Figures 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_{\rm star}$. Also note that in Figure 6 we have not applied any selection bias (unlike in the case of the semi-analytic model in Figure 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 Figures 5 and 6 clearly show that correlations between the residuals of the $M_{\rm bh}$-$\sigma$ and $M_{\rm bh}$-$M_{\rm 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 normalisation of the black-hole scaling relations, at face value our model clearly fails to reproduce the observed $M_{\rm bh}$-$\sigma$ and $M_{\rm bh}$-$M_{\rm 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 onto 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_{\rm bh}$-$M_{\rm star}$ relation (top right panel of Figure 2). However, this would then spoil the good match to the velocity dispersion-stellar mass relation. Similarly, one could increase the accretion onto 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_{\rm bh}$-$M_{\rm star}$ relation. This however would proportionally increase the $M_{\rm bh}$-$\sigma$ relation above the data, when velocity dispersions are chosen to faithfully track those from the observed SDSS $\sigma-M_{\rm star}$ relation (magenta lines in the left panels of Figure 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_{\rm 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_{\rm 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_{\rm bh}$-$\sigma$ relation, inclusive of observational bias (e.g., middle and bottom left panels of Figure 2), is approximately $M_{\rm 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_{\rm star}\gtrsim 5\times 10^{11}\,M_{\odot}$. In this respect, our model does not support the conjecture put forward by Batcheldor et al. (2007), according to which the $M_{\rm 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 Figure 11), in which a very broad distribution of black holes extending to the lowest masses in all types of galaxies is highly disfavoured. Nevertheless, Figure 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_{\rm star}\lesssim 3\times 10^{11}\,M_{\odot}$, extending down to black holes of mass $M_{\rm 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 Figures 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_{\rm bh}$-$\sigma$ relation and those in the $M_{\rm star}$-$\sigma$ relation, at fixed velocity dispersion. However, the model predicts almost no correlation between the residuals of the $M_{\rm bh}$-$M_{\rm star}$ relations and those of the $\sigma$-$M_{\rm 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_{\rm 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 onto the stellar and dark-matter dynamics (Peirani et al., 2016). Nonetheless, 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:

At $M_{\rm 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_{\rm bh}\lesssim 10^{7}\,M_{\odot}$ are predicted in such massive galaxies. However, for galaxies having $M_{\rm star}\lesssim 3\times 10^{11}\,M_{\odot}$, the distribution of $M_{\rm bh}$ at fixed $\sigma$ or $M_{\rm star}$ is broad. 2. 2.

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

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_{\rm bh}$-$\sigma$ and $M_{\rm bh}$-$M_{\rm star}$ relations. 4. 4.

Galaxy-evolution models (our semi-analytic one as well as the Horizon-AGN simulation) predict almost no correlation between the residuals of the $M_{\rm bh}$-$M_{\rm star}$ relations and those of the $\sigma$-$M_{\rm 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).

Acknowledgments

We thank M. Volonteri for insightful conversations. EB acknowledges support from the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904 and from the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche. 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.