{'fraction_non_alphanumeric': 0.02960614279711197, 'fraction_numerical': 0.013848034534133018, 'mean_word_length': 4.517441247760565, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constraining black hole-galaxy scaling relations from the large-scale clustering of Active Galactic Nuclei and implied mean radiative efficiency

A supermassive black hole has been found at the centre of nearly every galaxy observed with sufficient sensitivity. The masses of these black holes are observed to increase with either the total mass or the mean (random) velocity of the stars in their host galaxies. The origin of these correlations remains elusive. Observational systematics and biases severely limit our knowledge of the local demography of supermassive black holes thus preventing accurate model comparisons and progress in this field. Here we show that the large-scale spatial distribution of local active galactic nuclei (AGN), believed to be accreting supermassive black holes, can constrain the shape and normalization of the black hole-stellar mass relation thus bypassing resolution-related observational biases. In turn, our results can set more stringent constraints on the so-called"radiative efficiency", a fundamental parameter describing the inner physics of supermassive black holes that is closely linked to their spin, geometry, and ability to release energy. The mean value of the radiative efficiency can be estimated by comparing the average total luminous output of AGN with the relic mass density locked up in quiescent supermassive black holes at galaxy centres today. For currently accepted values of the AGN obscured fractions and bolometric corrections, our newest estimates of the local supermassive black hole mass density favour mean radiative efficiencies of ~10-20\%, suggesting that the vast majority of supermassive black holes are spinning moderately to rapidly. With large-scale AGN surveys coming online, our novel methodology will enable even tighter constraints on the fundamental parameters that regulate the growth of supermassive black holes.

A supermassive black hole has been found at the centre of nearly every galaxy observed with sufficient sensitivity. The masses of these black holes are observed to increase with either the total mass or the mean (random) velocity of the stars in their host galaxies. The origin of these correlations remains elusive. Observational systematics and biases severely limit our knowledge of the local demography of supermassive black holes thus preventing accurate model comparisons and progress in this field. Here we show that the large-scale spatial distribution of local active galactic nuclei (AGN), believed to be accreting supermassive black holes, can constrain the shape and normalization of the black hole-stellar mass relation thus bypassing resolution-related observational biases. In turn, our results can set more stringent constraints on the so-called "radiative efficiency", ε, a fundamental parameter describing the inner physics of supermassive black holes that is closely linked to their spin, geometry, and ability to release energy. The mean value of ε can be estimated by comparing the average total luminous output of AGN with the relic mass density locked up in quiescent supermassive black holes at galaxy centres today. For currently accepted values of the AGN obscured fractions and bolometric corrections, our newest estimates of the local supermassive black hole mass density favour mean radiative efficiencies of ε ∼ 10 − 20%, suggesting that the vast majority of supermassive black holes are spinning moderately to rapidly. With large-scale AGN surveys coming online, our novel methodology will enable even tighter constraints on the fundamental parameters that regulate the growth of supermassive black holes.
In the traditional picture, quiescent supermassive black holes at the centres of local galaxies today are the relics of a previous phase of active gas accretion in which they shone as Active Galactic Nuclei (AGN) 1 . The brightness of an AGN depends on how efficiently the gas orbiting the black hole can radiate. In principle, the mean radiative efficiency ε (and hence black hole spin 2 ) can be constrained by equating 3 the mass accreted by all AGN over all times with the mass density in supermassive black holes today [4][5][6][7] . The accreted mass scales roughly with the inverse of the mean radiative efficiency, because a lower radiative efficiency requires a higher mass accretion rate to produce the observed AGN luminosity. In this scenario, the same total AGN emissivity can be matched either by a small black hole mass density today but a large radiative efficiency, or vice-versa. This has fueled interest in estimating the local black hole mass density or "mass function", i.e. the number density per unit comoving volume and bin of black hole mass.
Quiescent black holes are difficult to observe, so large samples are not available. Therefore, the black hole mass function is usually estimated in two steps. First, one uses the small sample that is available to calibrate how a black hole's mass correlates with other properties of its host galaxy, such as its stellar mass M star or the velocity dispersion of its stars σ. Since M star and σ can be estimated in much larger samples, one uses these "scaling relations" to transform the comoving number density of M star or σ into that for black hole mass 4 .
Unfortunately, despite its elegance, the method described above suffers from systematic uncertainties that prevent a robust estimate of ε, the mean radiative efficiency of supermassive black holes. The main difficulty is to estimate the local scaling relations reliably 8 . Over the past decade, a number of groups have put forward a variety of rather different relations. For instance, there are several claims that black holes hosted in galaxies belonging to different morphological classes follow different scaling relations 9,10 . In addition, active and quiescent black holes in the local Universe appear to follow somewhat differently normalized scaling relations, with the degree of the offset depending on the type of galaxies, AGN sample, or scaling relation considered [11][12][13] . Furthermore, limitations on instrumental capabilities or other observational effects may systematically bias the observed scaling relations away from the intrinsic ones. This is a concern especially in quiescent early-type galaxies with dynamically-measured black hole masses 8,14,15 . In fact, recent work 13 suggests that such a bias may explain much of the reported difference between active and quiescent scaling relations: the intrinsic relations may be much more similar than they appear.
In what follows, we present a methodology that uses AGN clustering to set new and valuable independent constraints on the overall shape of the local black hole scaling relations. As we describe below, simultaneously matching the large-scale AGN clustering and the local black hole mass density provides more robust constraints on the mean black hole radiative efficiency ε.

The connection between AGN clustering, scaling relations and radiative efficiency
Within the framework of a cold dark matter Universe, more massive host dark matter haloes appear progressively more strongly clustered, i.e., their spatial distribution shows more marked departures from an underlying random distribution 16 . Galaxies and black holes residing in more massive dark matter haloes are thus naturally expected to appear more clustered. Our methodology builds on this basic notion of galaxy clustering and it can be briefly outlined as follows: • At any redshift of interest z, we first create large catalogues of host dark matter haloes from the halo mass function n(M halo ). In practice, in a given (large) volume V we select the halo masses of mass M halo for which the cumulative halo mass function is an integer, i.e. V × n(> M halo ) = N , with N ∈ I. • Finally, the match between the local mass density and accreted mass density determined from integrated AGN number counts, will yield more secure constraints on the mean radiative efficiency of supermassive black holes.

Results
The local M bh -M star relation and implied M bh -M halo relation We start by discussing the M bh -M star relation between black hole mass and host galaxy (total) stellar mass, which is perhaps the most actively debated correlation between black holes and galaxies ( Figure 1) yet widely adopted as a reference in black hole-galaxy co-evolution cosmological models 19 8 is reported by a solid red line with its scatter marked by the yellow region. Hereafter, we will refer to the de-biased M bh -M star relation put forward by Shankar et al. 8 (solid red line in the left panel of Figure 1, with a scatter of 0.4 dex, marked by the yellow band) as the "intrinsic" or "unbiased" scaling relation.
The very first feature to highlight in the left panel of Figure 1 is the systematic discrepancy between the fits to the local quiescent samples of supermassive black holes and Shankar et al.'s intrinsic M bh -M star relation. Shankar et al. 8 devised targeted Monte Carlo simulations to show that local quiescent (mainly early-type) galaxies having dynamically measured black hole mass, present larger velocity dispersions at fixed stellar mass with respect to galaxies in the Sloan Digital Sky Survey (SDSS). In fact, given the strong and fundamental 8,13,14,28 dependence M bh ∝ σ 4−5 , higher velocity dispersions would on average select higher mass black holes at fixed host galaxy stellar mass. Thus, the local quiescent samples of (mainly early-type) galaxies with dynamically measured black holes would then follow biased scaling relations, with fictitiously higher normalizations than the bulk of the underlying population of supermassive black holes. To a large extent, Shankar et al. 8 showed that a possible reason behind the higher velocity dispersions in the local sample of quiescent early-type galaxies could be ascribed to a bias arising from the requirement that the central black hole's sphere of influence must be resolved to measure black hole masses Figure 1: Overview of local scaling relations between black hole mass and host (total) stellar mass. Left: Correlations between central black hole mass and host galaxy total stellar mass in the local Universe. The black long-dashed line with its scatter (grey band) is the relation inferred from the (whole) quiescent sample by Savorgnan et al. 20 , while the purple, dotted line refers to only their early-type galaxies (ellipticals and lenticulars). The triple dot-dashed orange line is the fits to local quiescent samples of early-type galaxies with dynamical measures of black holes by Sahu et al. 21  with spatially resolved kinematics.
As already noted by Reines & Volonteri, the relations defined by the AGN samples are about an order of magnitude lower than the ones defined by the quiescent early-type samples, in line with several other previous claims 11,29,30 . In fact, Shankar et al. (2019) 13 have highlighted the fact that the vast majority of local AGN samples follow relations that are noticeably lower in normalization than those of quiescent galaxies and closer to the intrinsic relation (yellow region). They argue that AGN do not suffer from angular resolution-selection effects, so they may more faithfully trace the intrinsic black hole-galaxy scaling relations.
In summary, different host galaxy morphology, observational biases and selection effects all tend to result in different black hole mass-host galaxy scaling relations. To make progress, we need an independent constraint on these scaling relations, which we obtain as follows.
It is well known that galaxies with large stellar masses are more strongly clustered 16 . Therefore, it should be possible to set constraints on the normalization and steepness of the M bh -M star relation by including independent galaxy clustering measurements. The right panel of Figure 1 reports the z = 0.1 implied M bh -M halo relation of central galaxies obtained by assigning to each dark matter halo extracted from the parent halo mass function 31 , a galaxy stellar mass from the mean M star -M halo relation (and its scatter), and then a black hole mass from the M bh -M star relation (with scatter) for the quiescent Savorgnan et al. 20 8 . The latter M bh -M star relation clearly predicts a correlation between black hole mass and host halo mass that is more consistent with the variety of data extracted from low-redshift AGN large-scale optical surveys 23-25 and reported in the right panel of Figure 1. In what follows we always adopt as a reference for the "biased/observed" M bh -M star relation the one inferred from the Savorgnan et al. sample. As inferred from the left panel of Figure  Independent constraints on the M bh -M star relation from AGN large-scale clustering To further improve on the constraints that can be derived from AGN samples, we now use the large-scale AGN clustering measurements performed on the ROSAT All-Sky Survey and SDSS data sets 24 in the redshift range 0.16 < z < 0.36. AGN clustering is usually quantified by the two-point correlation function: the excess number of AGN pairs in the data compared to a random distribution of the same number density, as a function of separation. At scales larger than a few Mpc, we are in the linear regime of structure formation and the ratio of the two-point correlation function to that expected for the dark matter is approximately independent of pair separation. The square root of this ratio is called the "large-scale bias factor" b, and it is not to be confused with the bias in the scaling relations we were discussing about in the previous Sections! The large-scale bias factor b differs for different galaxy samples. The bias factor in fact progressively increases for larger halo masses 16 . Due to the monotonic dependence between stellar mass and host halo mass in the M star -M halo relation, especially at the masses of interest in this work (intrinsic scatter of 0.15 dex), samples with larger M star , will on average appear more strongly clustered and characterized by a larger bias factor b 16 . Since black hole mass is also monotonically related to host halo mass (albeit with a larger scatter than in the M star -M halo relation), we would thus expect the b factor to increase with black hole mass.
The two green squares with error bars in Figure 2 (same in both panels) show the large-scale bias factors for two samples of low-and high-mass active supermassive black holes observed in X-ray band. Similarly, the two magenta triangles show the bias factor b in a low-and high-mass sample of optically-selected AGN. The inferred bias factors b appear to be rather independent of wavelength, and indeed larger for the samples with larger black hole mass M bh .
The curves in the left-hand panel of Figure 2 show the expected large-scale bias assuming black hole masses are assigned via the observed/biased (long-dashed black, Equation 1) or intrinsic/unbiased (solid red) M bh -M star scaling relations 8 (Equations 2 and 3), with the same underlying mean M star -M halo relation (see the "Methods" Section). The predicted bias parameter b based on the observed/biased scaling relations is too low and is ruled out at high confidence ( 3σ), whilst that based on the intrinsic/unbiased scaling relations provides a good match to the data. The agreement is remarkable, given that there are no free parameters to fit. In the Methods Section we provide a detailed comparison with the measured cross-correlation of AGN-SDSS sources and careful estimates of the significance of the discrepancies between models and data. The dashed green curve shows that the Reines & Volonteri 12 M bh -M star scaling relation also yields a good match to the clustering measurements, implying that the AGN clustering measurements are still consistent with M bh -M star relations that are even lower in normalization than Shankar et al.'s intrinsic relation, further strengthening the key results of this work. For the Reines & Volonteri relation we assume a nominal intrinsic scatter smaller of 0.3 dex, of the order of what inferred by the Authors, but slightly smaller than the ones characterizing the scaling relations of dynamicallymeasured inactive black holes 13 .
It has been argued that host galaxy morphology may play a substantial role in creating the systematic offset between black hole scaling relations of active and quiescent galaxies. In fact the former active black holes could be mostly hosted in later-type galaxies, which are expected to allegedly be hosting lower black hole masses at fixed stellar mass, while the latter quiescent samples with larger black hole masses tend to inhabit earlier-type galaxies 12 . However, the Davis et al. 10 steep relation of dynamically-measured quiescent black holes in late-type galaxies (blue triple dot-dashed line, with an intrinsic scatter of 0.7 dex), tends to still fall substantially below the bias data at high M bh . The latter result clearly supports the view of a bias in the observed M bh -M star relation which is independent of galaxy morphology or AGN type. We also expect the Krumpe et al. AGN clustering measurements, based on large, serendipitous samples, to be representative of a variety of host galaxy morphologies and not just late types. We finally note that all of the data in Figure 2 mostly refer to Type 1 AGN. Observationally, Type 2/obscured AGN have anyhow been measured to cluster at a comparable or even higher level, at least at z < 1 32, 33 , which would further strengthen our results in favour of lower normalizations in the M bh -M star intrinsic relation.
In summary, the low-redshift AGN clustering measurements considered here favour considerably lower normalizations of the M bh -M star relation. Scalings in M bh -M star with high normalizations such as the observed/biased one characterizing quiescent galaxies (black long-dashed line), tend to map black holes to lower stellar masses and thus to lower host halo masses (see Figure 1 For completeness, the right-hand panel of Figure 2 shows the bias parameter b(M bh ) expected from the M bh -σ scaling relation. Measuring the clustering strength implied by the M bh -σ relation represents an additional independent test to the overall reliability of the black hole scaling relations. This test is also particularly valuable as it circumvents the use of M star by directly applying abundance matching between the velocity dispersion-based local black hole mass function and halo (plus subhalo) mass functions (see the "Methods" Section for full details). Moreover, it has been claimed that the M bh -σ relation is less biased than the M bh -M star relation, implying that the former should be closer to the "intrinsic" relation, in terms of normalization, slope, dispersion around the mean 8,13,14 . It is thus expected that the related (large-scale) clustering properties should be less sensitive to observational biases. Indeed, when adopting reasonable scatter (∼ 0.3 dex) around the mean M bh -M halo relation 7, 34 , when adopting the observed (dashed black) and intrinsic (solid red) M bh -σ relations, in both cases we find mean b-M bh relations broadly consistent with the data, with the former only slightly disfavoured at the highest black hole masses.
Comparing local and accreted mass functions: Constraining the mean radiative efficiency of supermassive black holes We first compute the local mass function of supermassive black holes from the SDSS galaxy sample adopted by Shankar et al. 8 (full details given in the Methods Section). The 1σ error in number densities at fixed black hole mass are shown in Figure 3 for the unbiased and biased scaling relations, respectively. The left panel of Figure 3 points to a large difference between black hole mass functions computed from the observed and intrinsic M bh -M star scaling relations (green and cyan areas, respectively; the vertical grey band marks the region in which black hole mass function estimates become progressively less secure). Of uttermost impor-tance, the black hole mass function derived from the observed M bh -σ relation (dashed red line) still presents noticeable departures from the one derived from the observed M bh -M star relation (green region). As proven in Figure 5, this behaviour is not a consequence of neglecting bulge corrections but it reflects an internal inconsistency in the samples of local quiescent black holes 8,13,14 . On the other hand, the black hole mass function derived from the intrinsic M bh -σ relation 8 (long-dashed blue line in Figure 3) is fully consistent with the black hole mass function derived from the intrinsic M bh -M star relation (cyan region; Equation 2).
In the right panel of Figure 3 we then compare the local black hole mass functions with the accreted mass functions derived from AGN number counts and continuity equation techniques. We show results for the accreted black hole mass function computed in two ways, via the method outlined in Shankar et al. 6 (long-dashed red lines) and as proposed by Aversa et al. 7 (triple dot-dashed purple lines). For both models we adopt the same input X-ray AGN luminosity function 36 with bolometric correction from Yang et al. 37 , and a mean radiative efficiency of ε = 0.15. In all accretion models we include obscured AGN with column densities up to log N H /cm −2 25, as observed in a number of studies 38 . The two sets of lines for each model bracket the systematic uncertainty in including or excluding obscured AGN with higher column densities, 25 < log N H /cm −2) < 26, following the N H column density distributions provided by Ueda et al. 36 . Despite the substantially different input Eddington ratio distributions, both accretion models predict similar accreted black hole mass functions that line up well with the local black hole mass function derived from the intrinsic M bh -M star relation (cyan area). Broadly matching the local black hole mass function derived from the observed M bh -M star relation would instead require at least one order of magnitude lower radiative efficiencies (dotted red lines and yellow area). It is interesting to note the good match between local and accreted mass function even at the largest mass scales. As discussed in previous works 6,7 , frequent black hole mergers tend to generate a relatively long high-mass tail in the accreted black hole mass function. The latter would be in tension with the nearly-exponential fall off of the local black hole mass function, unless an appropriately fine-tuned radiative efficiency progressively increasing with black hole mass is adopted 6 .
In integral form, Figure 4 shows the comparison between the local black hole and accreted mass densities. The former is obtained from direct integration of the local black hole mass function derived from the M bh -M star relation (long-dashed black box, corresponding to integration of the cyan region in Figure 3; also included, for completeness, the mean local black hole mass density implied by the intrinsic M bh -σ relation, pink square). The latter accreted mass density is retrieved from integration of the AGN bolometric luminosity function over luminosity and cosmic time (see Methods Section). The three coloured bands in each panel of Figure 4 mark three different values of the radiative efficiency, as labelled in each panel. Each coloured band brackets the systematic uncertainty induced by including or neglecting the Compton-thick AGN with 25 < log N H /cm −2 < 26. The top row includes model outputs that make use of the Yang et al. 37 and Marconi et al. 5 bolometric corrections, respectively. The right panels also allow for an additional loss of rest-mass energy via, e.g., jets and/or winds, parameterized by the kinetic efficiency which we set 39 to an average value of ε kin = 0.15. It is interesting to note that Yang et al. 37 also computed  6 and Aversa et al. 7 accretion models (dotted red and purple triple dot-dashed lines, respectively; the coloured areas in between the lines, yellow and pink respectively, represent the uncertainty region from excluding/including Compton-thick AGN with 25 < log N H /cm −2 < 26). Irrespective of the input Eddington ratio distributions, both accretion models, which use a radiative efficiency of ε = 0.15, are in good agreement with the unbiased local black hole mass function (cyan region). The grey areas in both panels mark the regions in which local measurements are less robust and cannot be constrained by clustering data. Figure 4: Comparing local and accreted integrated mass densities. "Relic" black hole mass density ρ SMBH inferred from the local unbiased 8 M bh -M star and M bh -σ relations (black dashed box and solid pink square, respectively), compared with the accreted mass density ρ AGN in AGN assuming different values of the radiative efficiency ε, bolometric correction BC, and kinetic efficiency ε kin , as labelled. The coloured region for each model shows the uncertainty from excluding/including Compton-thick AGN with 25 < log N H /cm −2) < 26. For the range of currently accepted values of ε kin , BC and obscured fraction, the data indicate moderate-to-high radiative efficiencies (ε 0.1).
the redshift-dependent M bh -M star relation from integrating the AGN X-ray emissivity at different epochs along (putative) stellar mass accretion histories. By assuming a constant mean radiative efficiency of ε=0.1 to convert from AGN bolometric luminosities to black hole accretion rates, they claimed a local M bh -M star relation in close agreement to the unbiased one by Shankar et al. 8 , and definitely steeper than the one from Kormendy & Ho 9 .
It is clear from Figure 4 that, irrespective of the exact assumed value of the obscured fraction, bolometric correction or kinetic efficiency, all accretion models always require relatively high radiative efficiencies ε ∼ 0.10 − 0.20 to match the local black hole mass density. This yields a mean radiative efficiency a factor of ∼ 2 − 3 higher than previous estimates solely based on observed/biased black hole scaling relations 5,36 . Our results would support black hole spin parameters 2 a 0.7, which are consistent with a number of direct measurements performed via X-ray reflection analysis 40 or spectral energy distribution fits 41, 42 . Hints for relatively high mean radiative efficiencies are also present in some previous works 43, 44 . However those studies adopted biased estimates of the local black hole mass density, sometimes relatively crude estimates of AGN emissivities, and without the vital constraints from AGN clustering. In this work we have ultimately shown that the AGN large-scale clustering as a function of black hole mass represents a fundamental discriminator for black hole-galaxy scaling relations. Upcoming large-scale AGN surveys such as eROSITA 45 , probing millions of AGN, will generate exquisite measurements of AGN clustering, set stringent constraints of the scaling relations down to the lowest accessible black hole masses, and in turn tighter constraints on the radiative and kinetic efficiencies of accreting supermassive black holes.
with an intrinsic scatter of 0.5 ± 0.1 dex.
The zero points of the correlations from Davis et al. 10 and Sahu et al. 21 , their Equations 3 and 11, respectively, have been modified setting the parameter v in their fits equal to v = 0.64, to allow for conversion to the mass-to-light ratios of Bell et al. 46 with a Chabrier 47 stellar initial mass function adopted as a reference in Shankar et al. (2016) 8 and in this work. The Reines & Volonteri fit 12 and the Kormendy & Ho 9 sample, which is also taken from Table 3 with a mass-dependent intrinsic scatter in (the logarithm of) black hole mass at fixed stellar mass (yellow region) well reproduced by 49 Both Equations 2 and 3 are good fits 8, 13 to galaxies with stellar mass above a few times log M star /M > 10. As recently discussed by Shankar et al. (2019) 13 , Equation 2 was obtained by imposing that black holes closely follow, via the M bh ∝ σ 4−5 relation, the velocity dispersion-stellar mass of SDSS early-type galaxies. Any significant contribution from later-type galaxies would tend, if anything, to produce slightly lower normalizations of the unbiased M bh -M star , thus further increasing the mismatch with the scaling relation of quiescent galaxies.
Abundance matching technique. The M star -M halo relation is obtained from abundance matching between the stellar mass and the halo plus subhalo mass functions 17,50 . We here adopt the Tinker et al. (2008) halo mass function 31 and the unevolved, surviving subhalo mass function by Jiang and van den Bosch 51 (see their Section 4.5). The latter is the relevant statistical description of subhaloes that survived mergers down to a given redshift, and thus the most suited quantity to include in accurate abundance matching algorithms. However, we verified that adopting different, more general fits of the subhalo mass function 52 would yield very similar results.
We assume throughout an intrinsic scatter of 0.15 dex in stellar mass at fixed halo mass 17 . We include scatter following Equation 5 in Aversa et al. 7 , though very similar results are found by following the methodology presented in Shankar et al. (2017) 17 . We adopt the Sérsic+Exponential stellar mass function from Bernardi et al. 2013 53 , which is based on the same photometry and mass-to-light ratios used by Shankar et al. (2016) and in this work. We assume no evolution in the stellar mass function up to z = 0.25, which is a very good approximation according to the latest estimates 54,55 .
The M star -M halo relation is now more securely established at least for relatively high stellar masses M star 3 × 10 10 M and at relatively low redshifts 17,18,53,[56][57][58] . We then extract large host halo catalogues from the halo plus subhalo mass function and assign galaxies to haloes via the M star -M halo relation. We note that all central and satellite galaxies are here assigned from the M star -M halo relation at the single redshift of z = 0.25, which is an extremely good approximation to a full multi-epoch abundance matching routine, especially for relatively massive galaxies (M star 3 × 10 10 , M ), as recently demonstrated by Grylls et al. 18 . In this work we always assume, unless otherwise noted, a standard reference cosmology of Ω m = 0.3, h = 0.7, σ 8 = 0.8, n s = 1, as adopted in the observational papers we compare with.
Computing local and accreted mass functions. The local black hole mass function is computed from the SDSS galaxy sample adopted by Shankar et al. 8 . In a Monte Carlo fashion, to each galaxy a supermassive black hole is assigned following a given input M bh -M star relation. To allow for moderate systematics in the galaxy stellar mass estimates due to fine details in the light profiles and/or mass-to-light ratios, at each iteration we randomly add ±0.1 dex to the stellar masses of each host galaxy, and then assign to it a supermassive black hole from the M bh -M star relation inclusive of its scatter. At each iteration, the black hole mass function is then computed from the V max weights assigned to each SDSS galaxy. The final 1σ error in number density at fixed black hole mass is then given by the dispersion in number densities resulting from 100 realizations. Supplementary Table 1 reports the results of the black hole mass function and its 1σ uncertainty at fixed black hole mass derived from the M bh -M star relation. The black hole mass function implied by the M bh -σ relation is fully consistent with the one given in Supplementary Table 1 within the uncertainties, as shown in Figure 3.
In Supplementary Figure 5 we show the results of applying the Monte Carlo methodology described above to only elliptical galaxies, which we select from the SDSS sample imposing a Bayesian classification probability 59 P [E] > 0.5 (our results do not depend on the exact value of this threshold). The left panel shows the comparison between the black hole mass functions implied by the observed M bh -M star and M bh -σ relations of early-type galaxies 13 . A clear mismatch of the order of ∼ 2−3 in integrated black hole mass density, is evident between the black hole mass functions, proving that the difference in black hole mass estimators is not a simple byproduct of including or not bulge corrections in the calculation. On the other hand, the right panel of Figure 5 shows the good match between the black hole mass functions implied by the unbiased/intrinsic M bh -M star and M bh -σ relations of early-type galaxies.
To compute the accreted (or "relic") black hole mass function we instead solve the continuity equation [60][61][62] ∂n bh ∂t where Ṁ bh is the mean accretion rate (averaged over the active and inactive populations) of the black holes of mass M bh and number density n bh (M bh , t) at time t.
In the main text we show results for the accreted black hole mass function computed in two ways, via the method outlined in Shankar et al. 6 and as proposed by Aversa et al. 7 . In the former, Equation 4 is solved by assuming 63 an input Eddington ratio distribution (lognormal) and constantly peaked around an Eddington ratio log λ = −0.6, with λ ∝ L/M bh ) of active black holes convolved with the active mass function to yield the observed AGN luminosity function at any cosmic epoch. In the latter, following the seminal work by Yu & Lu 64 , a light curve is instead assumed in input and the active and total black hole mass functions are then self-consistently derived by direct analytic time integration. Aversa et al. also adopted a slim and thin-disk approximation, in which the radiative efficiency is assumed to depend on the current Eddington ratio, in which the radiative efficiency can be close to 0.3 during super-Eddington accretion, but rapidly approaching the thin-disc constant value at Eddington/sub-Eddington regimes. Given that the latter is the most typical accretion mode for AGN, the Aversa et al. model reduces to a standard, constant radiative efficiency throughout most of the evolution time.
The normalization and shape of the implied duty cycle, i.e., the ratio between the active and total black hole mass functions, depend on the input Eddington ratio distributions P (λ), with λ ∝ L/M bh . As extensively discussed by Shankar et al. 6 , assuming, for example, a constant Eddington ratio distribution, gradually generates a very low duty cycle at fixed black hole mass at late epochs, which is largely inconsistent with present data 65 . On the other hand, an Eddington ratio distribution that drops towards lower values of λ(z) at lower redshifts, can provide orders of magnitude higher duty cycles, as the same black hole mass bins are mapped to lower luminosity and more abundant AGN. We checked that increasing the radiative efficiency from ε ∼ 0.05 to ε ∼ 0.2 does not change these behaviours in the output duty cycle.
Accreted mass densities. The classical 3, 4, 66 argument based on integrated mass densities, and inclusive of both radiative and kinetic efficiencies 39, 67, 68 , can be cast as follows in which the integrated (or "relic") mass density ρ AGN from accretion over all AGN with number densities Φ(L, z), bolometric luminosity L and redshift z, is compared to the local black hole mass density ρ SMBH . The integral in Equation 5 is computed over the luminosity range 42 < log L/erg s −1 < 48 and from z = 6. The exact limits of integration do not significantly change the results. The values of bolometric luminosity L in input into Equation 5, implicitly assume a bolometric correction k bol (L) that is luminosity dependent, folded in the computation of the bolometric luminosity function Φ(L, z) derived from Ueda et al.'s X-ray luminosity function.
Computing large-scale bias To compare models to clustering data, we first perform abundance matching at the average redshift of z = 0.25. We then extract large host halo catalogues from the halo (plus subhalo) mass function and assign galaxies to haloes via the M star -M halo relation as described above.
On the general assumption that black hole scaling relations do not depend on environment/host halo mass, to each mock galaxy in the catalogue we assign a black hole mass adopting an input M bh -M star relation inclusive of scatter. In principle, the mean large-scale bias for black holes with mass in the range M bh and M bh + dM bh could be straightforwardly estimated as where the sum runs over the total number N bin of (central and satellite) parent haloes hosting central black holes with mass in the range M bh and M bh + dM bh .
However, Equation 6 neglects the probabilities for central and satellite black holes to be active. In the case of AGN in fact, not all galaxies in a given host halo mass bin necessarily contribute to the same clustering signal, and Equation 6 should be modified to include the duty cycles, or probabilities, U cen (M bh ) and U sat (M bh ) of, respectively, central and satellite black holes of mass M bh to be active at a given cosmic epoch. We here follow the same convention as in Shankar et al. (2013) 6 and denote the duty cycles at fixed black hole mass M bh at a given redshift z as U (M bh , z). As in this work we only focus on the specific redshift of z = 0.25, we will drop the redshift dependence in the duty cycles from here onwards. Generalising Equation 6 to include the central/satellite probabilities of being active we obtain where both U cen (M bh ) and Throughout this work we adopt as a reference the total duty cycle U (M bh ) in Equation 7 as derived from the continuity equation model by Shankar et al. (2013) 6 , with constant lognormal Eddington ratio distribution. The exact normalization and shape of the total duty cycle U (M bh ) depends on the minimum luminosity threshold considered and input Eddington ratio distributions P (λ), as anticipated in the previous Section . As extensively discussed by, e.g., Shankar et al. 6 , continuity equation models generally tend to generate duty cycles decreasing with increasing black hole mass, as also retrieved by direct data modelling by other groups 6 . The latter trend would then imply that active lower mass black holes would have a larger weight on the bias (as in Figure 10 discussed below). However, when computing the mean bias in narrow bins of black hole mass (as in Figure 2), the exact shape or normalization of the assumed total duty cycle U (M bh ) are irrelevant. It is in fact clear from Equation 7 that, in a relatively small bin of black hole mass, what contributes to the mean bias is not the total duty cycle U (M bh ), rather the relative probabilities for central and satellite black holes to be active, i.e., the ratio Q bh (M bh ) = U sat (M bh )/U cen (M bh ).
Although satellite black holes are usually relatively less abundant than central black holes of similar mass, especially at progressively larger masses of interest to this work, if the ratio Q bh 1, satellite black holes could still noticeably contribute to the total mean bias as they inhabit very massive host dark matter haloes. Several works in the literature 69 implying where f BH sat (M bh ) = N sat (M bh )/[N sat (M bh ) + N cen (M bh )] is the total fraction of satellite black holes with mass within M bh and M bh +dM bh . Based on our mocks, the total fraction of satellites for the high-mass black holes considered by Krumpe et al. 24 , is f BH sat (M bh ) ∼ 0.1, which corresponds to Q bh (M bh ) ∼ 2 when adopting the best-fit f AGN sat = 0.18 from Leauthaud et al. 71 . We will always adopt the latter value of Q bh = 2 independent of black hole mass as our reference value throughout this work, though we will also show results with Q bh = 1. Indeed, Q bh could take on even lower values. Other groups measured a f AGN sat consistent with zero 69 . More recently Man et al. 73 found the SDSS completeness-corrected fraction of local AGN at fixed stellar mass is independent of environment, with the fraction of AGN in satellite galaxies being comparable, if not lower, than that in central galaxies (see their Figure 4). The Man et al. results would imply Q bh 1 largely independent of host stellar mass, as assumed here.
When computing the mean bias in Figure 2 we make use of Equation 7 with Q bh = 2 and the halo bias b h is taken from Tinker et al. (2005) 74 , in line with the Halo Occupation Distribution models adopted by Krumpe et al. 24 . We note that Tinker et al. make use of FoF-based halo masses. We thus correct our reference halo masses, defined as 200 times the critical density, by an average factor 75 of 0.966 before applying the Tinker et al. analytic model, though this has a very minor effect on the results. For completeness, we proceed below to a more thorough comparison with the Krumpe et al. 24 full cross-correlation function between RASS (SDSS) AGN and the Luminous Red Galaxies (LRGs), as shown with filled black squares in Figure 7 for the SDSS Data Release 4 (top panels) and Data Release 7 (bottom panels). The (projected) cross-correlation function at large scales is analytically computed at any redshift of interest z as where b AGN and b LRG are the median large scale biases of AGN and LRGs, respectively, and w DM (r P , z) is the projected linear matter correlation function derived from the linear matter power spectrum of Smith et al. (2003) 76 , though we verified that other power spectra yield nearly identical results.
The bias of LRGs is computed by fitting the predicted linear auto-correlation function Fixing b LRG = 2.25 in Equation 10, we find, when integrating over a scale of 3 − 30 Mpc and including the full covariance in the χ 2 fitting, b AGN = 1.43 ± 0.11 and b AGN = 1.56 ± 0.10 for DR7 and DR4, respectively. In our mocks we then select a black hole mass distribution consistent with the one adopted by Krumpe et al., i.e., with a lognormal distribution peaked at log M bh /M = 8.5 and with a dispersion of 0.3 dex (the "high M BH sample" in their Figure 5). We then compute the mean bias as given in Equation 7 (we stress that medians would, if anything, tend to produce lower biases) assuming Q bh = 1 (left panels) and Q bh = 2 (right panels), which would correspond to a fraction f sat ∼ 0.1−0.2, following Equation 9. When computing the full cross-correlation function we use the more recent analytic model for the bias put forward by Tinker et al. (2010) 79 , which provides a closer match to the halo autocorrelation function from the MultiDark simulation 80 Figure 2. We note that all our assumptions so far have been as conservative as possible, aimed at minimizing the tension in the large-scale clustering between the data and the model based on the observed M bh -M star relation. Adopting even lower Q bh values, and/or more up-to-date observed M bh -M star relations 10 , and/or empirically-based duty cycles 81 , would all tend to yield mean values of the predicted b AGN close to unity, implying > 3 − 4σ departures from the mean bias directly retrieved from the projected correlation function. The implied cross-correlation functions from the intrinsic and observed models, found by setting b LRG = 2.25 into Equation 10, are reported in the Supplementary Figure 7 with solid red and long-dashed black lines, respectively.
Systematics in mean radiative efficiency and obscured fractions This work puts forward a novel methodology that makes use of AGN clustering measurements to set firmer constraints on the shape of the local M bh -M star relation and, in turn, on the mean radiative efficiency of AGN. Some key points should be emphasized. When comparing to a given rendition of the local M bh -M star relation, the accretion models are actually able to constrain the ratio between bolometric correction and radiative efficiency k bol (1 − ε)/ε (see Equation 5).
We showed in Figure 4 that, when adopting the latest renditions of the (luminosity-dependent) bolometric corrections by Yang et al. 37, 82 , we require mean radiative efficiencies of the order of ε 0.1 − 0.2 to match the local black hole mass density implied by the intrinsic/unbiased M bh -M star scaling relation 8 . Up to one order of magnitude lower mean radiative efficiencies ε ∼ 0.02, would instead be necessary to broadly align with the local black hole mass density inferred from the observed M bh -M star relation (yellow region in the right panel of Figure 3).
Ueda et al. 36 performed very similar calculations to ours and with the same AGN luminosity function. When adopting as a local reference a black hole mass function derived from the Kormendy & Ho M bh -M star scaling relation, comparable to the dot-dashed black lines in Figure 3 (filled red circles in their Figure 23), they derived a best-fit radiative efficiency of ε = 0.053. Correcting their estimate of the mean radiative efficiency based on a different bolometric correction 83 , which is a factor of ∼ 2 − 3 higher than our reference one, would imply extremely low values of ε 0.03, fully consistent with our "biased" estimate of the radiative efficiency. The previous results by Ueda et al. and other works would thus be in tension with estimates of independent spectral estimates of the spin of supermassive black holes 40 , which suggest radiative efficiencies of the order of ε 0.1 or above.
Higher values of the bolometric corrections would in fact require proportionally higher radiative efficiencies to generate the same relic black hole mass density (Equation 5). Indeed, Zhang et al. 84 14), though with noticeably less massive black holes (their Figure 15), implying an overall lower integrated black hole mass density falling in between the models presented in the top and bottom panels of Figure 4. A mean radiative efficiency decreasing with cosmic time from close to maximal (ε ∼ 0.2−0.3) to minimal values (ε ∼ 0.05) at low redshift, was originally considered by Shankar et al. as a possibility to reconcile the preliminary analysis of high-redshift clustering data of luminous quasars, which already suggested relatively high radiative efficiencies 87 , with the "biased" local black hole mass densities. A progressively decreasing radiative efficiency at late times is however disfavoured by our current results as it would generate relic black hole mass densities too large with respect to our updated local black hole mass densities.
Another important point to mention, which can be of relevance for all statistical supermassive black hole accretion models, is the possible presence of a selection effect on black hole spin distribution affecting flux-limited AGN surveys 88 . Assuming an intrinsically broad spin distribution at, say, fixed black hole mass, flux-limit effects would then favour the detection of the more luminous, higher-spin/higher-ε sources. This in turn would imply that the average intrinsic radiative efficiency could be lower at fixed black hole or host galaxy stellar mass. On the other hand we note that even assuming a minimum of ε = 0.05, setting 89 L = f 0.1ṁ lM bh , with f = ε/(1 − ε), l ∼ 1.3 × 10 38 erg s −1 M −1 ,ṁ =Ṁ bh /Ṁ Edd , andṀ Edd = L Edd /0.1c 2 , we would get a bolometric luminosity of L 10 43 erg s −1 , for M bh 10 6 M accreting at an average 6ṁ ∼ 0.1. This limit in bolometric luminosity would correspond to L X 10 42 erg s −1 , for typical X-ray bolometric corrections at these luminosities, which is still within the reach of the surveys adopted in this work, at least at z 2 (see left panel of Figure 3 in Ueda et al. 36 ).
Last but not least, a fundamental source of systematic uncertainty in the calculation of the mean radiative efficiency derived from Equation 5, originates from the fraction of obscured AGN included in the AGN luminosity function. Equation 5 in fact implies specifying a complete census of AGN, and thus the AGN luminosity function must be corrected for the fraction of obscured sources f obsc (L, z) missed even in hard X-ray AGN surveys. In this work we followed Ueda et al. 36 , who suggest a fraction of obscured Compton-thick AGN with column densities in the range 24 < log N H /cm −2 < 26 equal to 50% the number densities of AGN with 20 < log N H /cm −2 < 24. The Ueda et al. estimate is in good agreement with the more recent extrapolated estimates from NuSTAR observations of heavily obscured Swift/BAT AGN by Georgantopoulos and Akylas 90 , who predict, from analyzing a fraction of the BAT sample, up to 50% of Compton-thick AGN with respect to the total population, within typical AGN spectral parameters. On the other hand, from cutting-edge comprehensive AGN population synthesis models, Ananna et al. 91 have recently claimed a fraction of Compton-thick AGN with 24 < log N H /cm −2 < 26 a factor of ∼ 2 − 3 higher than what estimated by Ueda et al. (their Figure 10). Increasing the number densities of AGN increases the normalization of the AGN luminosity function Φ(L, z) in Equation 5, implying a higher relic mass density at fixed bolometric correction, thus proportionally higher radiative efficiencies to match the same local black hole mass density.
We note that in hierarchical models of structure formation 75 , a substantial population of won-dering black holes may end up co-existing with the central black hole within the same host galaxy. Accounting for this additional, at present unknown, population of wondering black holes in the estimate of the local black hole mass density would imply a decrease in the mean radiative efficiency if the wondering black holes were themselves active and recorded in the AGN luminosity function. More generally, orbiting and/or ejected black holes could affect the applicability of Equation 5. Cosmological hydrodynamic simulations 92 showed, however, that the correction induced by the integrated contribution of black holes today not locked up at the centre of galaxies, should not amount to more than ∼ 11%, on average. All in all, our inferred mean radiative efficiencies of ε 0.1 − 0.2 can be safely considered as lower limits at fixed kinetic efficiency ε kin and local black hole mass density, as our adopted bolometric corrections and obscured fractions are comparable to, or if anything lower than, other estimates put forward in the literature. Figure 5: Predicted local black hole mass functions of elliptical galaxies. Same format as Figure 3 but only for Elliptical, bulge-dominated galaxies. Lef/right panel show the comparison between the observed/intrinsic M bh -M star and M bh -σ relations of early-type galaxies, as labelled.

Author contributions statement
All authors have contributed in different ways to the work presented in this paper. FS performed the full set up of the AGN mocks, analysis of the results, and writing up of the manuscript. VA independently checked all the results on AGN clustering and largely contributed to text editing and to the referee reports. MB was one of the core Authors in the Shankar et al. (2016) paper on the intrinsic black hole scaling relations and massively contributed to the editing of the manuscript. CM devised accurate determinations of the correlation functions in the MultiDark simulation. AL calculated the accreted black hole mass functions following the models presented in Aversa et al. (2015). NM contributed to the AGN accretion models. PJG and LZ contributed to the galaxy mocks, independently tested some of the key results in the paper, and provided comments. JM performed preliminary large-scale bias estimates of AGN at different luminosities in low redshift galaxies. MK made available a number of data sets on galaxy and AGN clustering inclusive of full covariance matrices. RDB performed independent calculations of some of the AGN mocks. FR contributed to the characterization of the scaling relations in Type 1 and 2 AGN and to the editing of the paper. FLF provided key insights in the use of different AGN bolometric corrections. RKS was one of the core Authors in the Shankar et al. (2016) paper on the intrinsic black hole scaling relations and provided support on the theoretical side and editing of the paper.

Additional information
Competing Interests The authors declare that they have no competing financial interests.
Correspondence Correspondence and requests for materials should be addressed to Francesco Shankar (email: F.Shankar@soton.ac.uk).