Stochastic Processes as the Origin of the Double Power-law Shape of the Quasar Luminosity Function

Keven Ren; Michele Trenti; Tiziana Di Matteo

doi:10.3847/1538-4357/ab86ab

1. Introduction

Quasars or quasi-stellar objects (QSOs) are highly luminous objects that are a class of active galactic nuclei (AGNs) powered by the accretion of matter from a disk onto the central supermassive black hole (SMBH). While quasars have been observed at all redshifts up to z ∼ 7.5 (Bañados et al. 2017), inferences of SMBH masses from luminous high redshift quasars (z > 6) can reach up to ∼10⁹M_⊙. The progenitors of these extreme quasars, and the mechanism behind their immense rapid growth remain a topic of open active research and a stringent theoretical hurdle to overcome. At face value, the relative rarity of the luminous quasars in question (comoving densities of the order ∼10⁻⁹ Mpc⁻³) would suggest that these objects reside at regions of extreme overdensities and contained within the most massive halos (Springel et al. 2005). Recent work, however, suggests this may not be the case (Fanidakis et al. 2013; Aversa et al. 2015; di Matteo et al. 2017) and that instead the first quasars may live in relatively more common halos. Thus, studying the connection between the early quasars to their host halos can provide insight to the origin of the first quasars and their subsequent evolution over cosmic time.

The role of AGN in galaxy formation has been a topic of considerable research interest in the past decade. The remarkably tight local M–σ relation hints toward a deep connection between a host galaxy and its respective black hole (Gebhardt et al. 2000; Kormendy & Ho 2013) and by extension, a connection between galaxies and their AGNs. The nature and the implications of this association have been explored both in observations (Fabian 2012; Page et al. 2012; Cicone et al. 2014) and theory/simulations (di Matteo et al. 2003; Croton et al. 2006; Lucia & Blaizot 2007; Sijacki et al. 2007; Jahnke & Macciò 2011; Conroy & White 2012, EAGLE: Schaye et al. 2014, Meraxes: Qin et al. 2017, SAGE/RSAGE: Croton et al. 2016; Raouf et al. 2017, Horizon-AGN: Beckmann et al. 2017, BLUETIDES: Feng et al. 2015; Ni et al. 2018). In this context, the luminosity functions of galaxies and AGNs are key observables that are easily measured owing to modern wide-area surveys such as the 2dF Galaxy Redshift Survey (2dFGRS: Colless et al. 2001), the Sloan Digital Sky Survey (SDSS, Richards et al. 2002), and the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (Grogin et al. 2011; Koekemoer et al. 2011). On one hand, the galaxy luminosity function and that of its underlying dark-matter halo mass function (HMF) are well described by Schechter-like shapes (a power law with an exponential drop-off at the bright/massive end) to all redshifts (Arnouts et al. 2005; Alavi et al. 2013; Schmidt et al. 2014; Weisz et al. 2014; Bouwens et al. 2015; Finkelstein et al. 2015), while the quasar luminosity function (QLF) has been canonically described at all redshifts by a double power law, or where faint end data is not available, by a single power law at its bright end (Croom et al. 2009; Schneider et al. 2010; Ross et al. 2013; Akiyama et al. 2017; Kulkarni et al. 2019; Shen et al. 2020). Understanding why and how this functional shape transition happens from the distribution of halos/galaxies to luminous quasars would refine our understanding of the paradigm around halo–galaxy–AGN interactions. Furthermore, it would not only contribute to elucidate the subsequent evolution of the QLF, but also potentially enable additional constraints for simulations/theories on the formation and growth of galaxies and SMBHs.

In this work, we investigate a physical origin for the difference in the shapes of the QLF and the HMF. Specifically, we look at the role stochasticity (or scatter) plays in determining the shape of the QLF. Then, we investigate constraints on feedback processes, AGN/galaxy quenching and local clustering that are set by the amount of scatter in the quasar luminosity versus halo mass relation. Scatter has been a core component surrounding discussions of clustering around quasars based on measurements of luminosity (Shen et al. 2007; Wyithe & Loeb 2009; Shankar et al. 2010), and more recently on feedback processes of galaxies in the largest halos (Ren et al. 2019). In principle, stochasticity in the quasar luminosity to halo mass relation facilitates the probability of an overluminous quasar to reside in a common lower mass halo instead of a very rare, massive halo. This effect is compounded at the massive-end of the HMF, because its Schechter-like shape makes lower mass halos exponentially more abundant than the rarest counterparts. Additionally, this also implies that adding scatter in the quasar luminosity versus halo mass relation increases the abundance of the brightest objects. This effect can thus be potentially used to constrain feedback processes if the observed bright end is sufficiently well constrained. The broadening of a distribution due to stochasticity is already well established for a number of observables, e.g., stellar mass and halo mass (Behroozi et al. 2010), galaxy luminosity and halo mass (Cooray & Milosavljević 2005; Ren et al. 2019), black hole mass and stellar mass (Hirschmann et al. 2010), black hole mass to stellar velocity dispersion (Volonteri & Stark 2011), stellar mass and AGN luminosity (Veale et al. 2014).

The inclusion of stochasticity in the QLF can be explicitly developed through the conditional luminosity function (CLF) approach that is routinely used in galaxy luminosity modeling (Yang et al. 2003; Cooray & Milosavljević 2005) just by constructing the median quasar luminosity versus halo mass relation in a way that is consistent with the observed QLF given the quasar duty cycle and the amount of log-normal scatter. The CLF is a powerful tool that links the quasar luminosity directly to its host dark-matter halo, and describes the distribution of luminosities inside a halo of given mass. Qualitatively, we expect scatter to suppress the massive end of the median quasar luminosity to halo mass relation (see Ren et al. 2019 for a galaxy analog).

In this work, we run such analysis for quasars at z ∼ 4, motivated from the availability of a new robust determination of the LF thanks to the recent survey from the Subaru Hyper Suprime Camera (Akiyama et al. 2017). In addition, this epoch operates as an ideal bridge to make statements on the environments that quasars inhabit at high redshifts (z > 6), as well as infer consequences in the galaxy-AGN evolution down to lower redshift.

This paper is structured as follows. In Section 2 we describe the method used to derive the median quasar luminosity versus halo mass relation given scatter and quasar duty cycle. In Section 3 we present the median quasar luminosity versus halo mass relations and we discuss the broader implication of this broadening under physical context of feedback, quenching and local clustering. Section 4 summarizes our findings and concludes. Throughout this paper, we use WMAP-7 cosmology with parameters, Ω_m = 0.272, Ω_b = 0.0455, Ω_Λ = 0.728, h = 0.704, σ₈ = 0.81, and n_s = 0.967 (Komatsu et al. 2011). We use the Jenkins et al. (2001) HMF. Magnitudes are given in the AB system (Oke & Gunn 1983).

2. Modeling

The typical approach to construct a relation between quasar luminosity and halo mass is through abundance matching which cumulatively matches the number densities of quasar luminosities with those of halos masses (Vale & Ostriker 2004). This technique is inherently a one-to-one matching process that does not include stochasticity expected to originate from more fundamental processes such as those driving SMBH growth. By taking stochasticity into account, we can expect a qualitatively different picture compared to a deterministic abundance matching. Specifically, the most luminous quasars are likely to be outliers in terms of accretion rate, i.e., overluminous for their host halo, instead of being hosted within the most massive halos (see Figure 1). Increasing the magnitude of scatter boosts the probability of the less massive halos to host an overluminous source, as halo abundance is strongly dependent on mass. In fact, compared to the galaxy (UV) luminosity function, it becomes even more critical to consider the effects of stochasticity in the QLF where observations typically extend well beyond the L > L_* regime. In this work, we largely follow the steps used in Allen et al. (2018) and introduce scatter into our modeling by defining a CLF for quasars. The method is summarized in the following subsections.

**Figure 1.** Leftmost three columns of the panels: Simulated z = 4 quasar magnitudes as a function of host halo mass as derived from our conditional luminosity function model, assuming a (1.5 Gpc)³ comoving volume, and different values for log-normal scatter Σ and quasar duty cycle _DC in each panel. The median of the data points is shown as a solid colored line. Arrows on the right edge of the bottom row of plots show quasars with the indicated luminosity inside halos beyond plot limits. The rightmost column has panels showing the resulting quasar luminosity functions from binning the model points as histograms (colored), compared to the input calibration function (solid black).
Download figure:
Standard image High-resolution image

epsilon — **Figure 1.** Leftmost three columns of the panels: Simulated z = 4 quasar magnitudes as a function of host halo mass as derived from our conditional luminosity function model, assuming a (1.5 Gpc)³ comoving volume, and different values for log-normal scatter Σ and quasar duty cycle _DC in each panel. The median of the data points is shown as a solid colored line. Arrows on the right edge of the bottom row of plots show quasars with the indicated luminosity inside halos beyond plot limits. The rightmost column has panels showing the resulting quasar luminosity functions from binning the model points as histograms (colored), compared to the input calibration function (solid black).
Download figure:
Standard image High-resolution image

2.1. Distribution of Quasar Luminosities

We explicitly model the stochasticity with a CLF approach. Here, the CLF ${\rm{\Phi }}({M}_{\mathrm{UV}}| {M}_{{\rm{h}}})$ can be interpreted as the probability distribution for quasar magnitudes, M_UV, given a halo mass M_h,

$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}({M}_{\mathrm{UV}}| {M}_{{\rm{h}}})=(1-{\epsilon }_{\mathrm{DC}})\delta (L=0)\\ \quad +\,\displaystyle \frac{{\epsilon }_{\mathrm{DC}}}{\sqrt{2\pi }(2.5\times {\rm{\Sigma }})}\exp \left(\displaystyle \frac{-{\left[{M}_{\mathrm{UV}}-{M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }},{\epsilon }_{\mathrm{DC}})\right]}^{2}}{2{\left(2.5\times {\rm{\Sigma }}\right)}^{2}}\right),\end{array}\end{eqnarray} \tag{ 1 }$

where ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ is the median quasar magnitude at our given halo mass, Σ is the width of the dispersion in dex, _DC is the halo-mass independent quasar duty cycle defined as the constant fraction of SMBHs that have active quasar-mode radiation and δ(x) is the Dirac delta function. The form of the CLF is log-normal and is justified by two reasons: (1) Observationally, the dispersion in the M_BH–σ relation is well fit with a log-normal (Gültekin et al. 2009). Having M_BH as a proxy for quasar luminosity and σ for halo mass could plausibly suggest a similar form of intrinsic scatter. (2) The luminosity of a quasar can be considered as a product of many processes, which can tend toward a log-normal dispersion by the central limit theorem. Hence, to first order, the scatter in this relation can be thought to be log-normal. The usual QLF, ϕ(M_UV) can then be derived by the equation,

$\begin{eqnarray}&&\phi ({M}_{\mathrm{UV}})={\int }_{0}^{\infty }\displaystyle \frac{{dn}}{{{dM}}_{{\rm{h}}}}{\rm{\Phi }}({M}_{\mathrm{UV}}| {M}_{{\rm{h}}}){{dM}}_{{\rm{h}}},\end{eqnarray} \tag{ 2 }$

where $\tfrac{{dn}}{{{dM}}_{{\rm{h}}}}$ is the HMF. The median halo mass versus quasar magnitude relation, ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ , given Σ and _DC can then be obtained by deconvolving Equation (2). The free parameters in this model are _DC, Σ, and the input QLF, ϕ. We adopt the z ∼ 4 QLF data set as reported by Akiyama et al. (2017). The form of the typical double power-law QLF is given by,

$\begin{eqnarray}&&{\phi }_{\mathrm{obs}}({M}_{\mathrm{UV}})=\displaystyle \frac{{\phi }^{* }}{{10}^{0.4(\alpha +1)({M}_{\mathrm{UV}}-{M}^{* })}+{10}^{0.4(\beta +1)({M}_{\mathrm{UV}}-{M}^{* })}},\end{eqnarray} \tag{ 3 }$

with the Akiyama et al. (2017) parameters we have ϕ* = 2.66 × 10⁻⁷ Mpc⁻³ mag⁻¹ as the LF normalization factor, M* = −25.36 as the characteristic break magnitude, α = −1.3 as the faint end slope, and β = −3.11 as the bright end slope. We apply two different deconvolution methods outlined in Allen et al. (2018) and Ren et al. (2019) to determine ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ through a least-squares fit.

2.2. Deconvolution Method

To prepare Figure 1 we follow the deconvolution method outlined in Allen et al. (2018), further discussed in Behroozi et al. (2010). The steps for this iterative process can be summarized as follows:

1.
We derive the QLF using Equation (2), given inputs Σ and _DC, and assuming the Jenkins et al. (2001) HMF. We start from an initial guess ${M}_{\mathrm{UV},{\rm{c}}}^{{\prime} }={M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }}=0,{\epsilon }_{\mathrm{DC}})$ , which can be derived from direct abundance matching, and construct ${\phi }_{M}({M}_{\mathrm{UV}}^{{\prime} })$ from Equation (2).
2.
We apply abundance matching between ${\phi }_{M}({M}_{\mathrm{UV}}^{{\prime} })$ and the input calibration QLF to derive a correction to the relation between ${M}_{\mathrm{UV}}^{{\prime} }$ and M_UV.
3.
We transform our median quasar magnitude versus halo mass relation, ${M}_{\mathrm{UV},{\rm{c}}}^{{\prime} }({M}_{{\rm{h}}},{\rm{\Sigma }},{\epsilon }_{\mathrm{DC}})$ according to the relation derived in Step (2).
4.
Steps (1) to (3) are iteratively repeated.
5.
The iteration process is terminated when the squared residual difference between successive iteration steps at a fixed number density, ϕ = 6 × 10⁻¹¹ Mpc⁻³ first changes by ${({\rm{\Delta }}{M}_{\mathrm{UV}})}^{2}\leqslant 0.01$ . This choice of the number density corresponds to the value of the fit in the brightest available magnitude data point of Akiyama et al. (2017).

This deconvolution process approximately derives ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ , which best matches the initial calibration QLF in Equation (2) through progressive improvements, with the caveat that numerical instabilities may arise for a large number of iterations.

3. Results and Discussion

3.1. Impact of Scatter on the Bright End of the QLF

Figure 1 shows one Monte Carlo realization for the distribution of quasar luminosities after sampling from the CLF and HMF, under different values of Σ. The number of sampled quasar points corresponds to an equivalent cosmological volume of (∼1.5 Gpc)³. The figure clearly demonstrates that our treatment of the median halo mass versus quasar magnitude ( ${M}_{\mathrm{UV},{\rm{c}}}$ ) relation successfully returns the z ∼ 4 QLF (within errors) after assuming values of Σ = 0, 0.3, and 0.5. We include choices for the quasar duty cycle at values of _DC = 0.01, 1.

In our modeling where Σ = 0, i.e., using the typical deterministic abundance matching approach, the resulting function ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }}=0)$ becomes increasingly steep at high halo masses. This is due to the need to abundance-match the exponential drop-off of the HMF at high masses with the power law bright-end of the QLF, and implies that hosts are required to be increasingly efficient at producing quasar luminosity at higher halo masses. We show this behavior in the green curve of Figure 2, highlighting a positive slope in the ratio of the median quasar luminosity to halo mass. In fact, a similar feature is also present for Σ = 0.3 (red curve, Figure 2), also yielding a positive slope, albeit the trend is less extreme. We note that the small fluctuations in this curve are just a consequence of the numerical deconvolution technique. For low values of Σ, the brightest quasars are generally hosted by the most massive halos. In this context, the ratio of median quasar luminosity per unit halo mass over halo mass still shows positive. Identifying a physical mechanism for such a requirement can be difficult in the presence of feedback processes that are regulating SMBH (and galaxy) growth at the high mass end in particular.

**Figure 2.** Quasar luminosity to halo mass ratio as a function of halo mass, $\mathrm{log}(\tfrac{L/{L}_{\odot }}{{M}_{{\rm{h}}}/{M}_{\odot }})$ for different values of scatter Σ (colored lines) and duty cycle _DC (solid vs. dashed lines) at $z\sim 4$ . For comparison, the same relation is also shown for galaxies (black solid line, assuming no scatter and unity duty cycle). The oscillations in the curves for Σ > 0 are from numerical instabilities in the Allen et al. (2018) deconvolution process.
Download figure:
Standard image High-resolution image

**Figure 2.** Quasar luminosity to halo mass ratio as a function of halo mass, $\mathrm{log}(\tfrac{L/{L}_{\odot }}{{M}_{{\rm{h}}}/{M}_{\odot }})$ for different values of scatter Σ (colored lines) and duty cycle _DC (solid vs. dashed lines) at $z\sim 4$ . For comparison, the same relation is also shown for galaxies (black solid line, assuming no scatter and unity duty cycle). The oscillations in the curves for Σ > 0 are from numerical instabilities in the Allen et al. (2018) deconvolution process.
Download figure:
Standard image High-resolution image

However, Figure 1 highlights that the shape of ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ changes when larger Σ values are considered. In the case of Σ = 0.5, the median halo mass for a luminous quasar is reduced. This effect is compounded by the exponential increase of the number density of less massive halos. Therefore, the more typical halos around the characteristic halo mass value ( ${M}_{{\rm{h}}}\lesssim {10}^{12.4}{M}_{\odot }$ ) begin to play a larger role in shaping the bright end of the QLF for sufficiently large Σ (Σ = 0.5 panel of Figure 1). As a result, the high-mass end of ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }}=0.5)$ flattens out in order for the CLF to successfully fit the observed QLF. Under these circumstances, the shape of the resulting ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ naturally resembles the functional shape of the galaxy luminosity versus halo mass relation (e.g., see Ren et al. 2019 Figure 1), and emphasizes a deeper fundamental connection between the galaxies and quasars. In this instance, ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ no longer implies the existence of an increasingly efficient mechanism for quasar accretion as the halo mass grows. Two conclusions can be drawn from the Σ = 0.5 panel: (1) the growth of BHs in the most massive halos is suppressed, and (2) the model suggests that the bright-end of the QLF is populated by objects that are outliers with very high accretion rate compared to the average value for their host halo mass.

Furthermore, the flattening induced in ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ by large Σ values effectively implies that the bright-end of the QLF becomes insensitive to changes in the value of the quasar duty cycle at high halo masses, as the objects are increasingly hosted by relatively common lower-mass halos. This is a generic feature of high Σ models and is independent of _DC. Thus, models with high values of Σ are justified in assuming a mass-independent duty cycle to describe the bright end of the QLF. The weak dependence on halo mass for _DC has been a key feature for several models in describing luminosities of quasars (e.g., Shankar et al. 2010; Conroy & White 2012).

The nature of Σ is not restricted to the intrinsic stochasticity in the fundamental processes that build the quasar itself. The dependence on environment also adds an effective "scatter" that is encapsulated in the Σ parameter. Thus a high value of Σ in ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ could also indicate that powering a quasar accretion mode is contingent on having ideal environmental conditions that funnel a sufficient amount of gas to the center of the host galaxy and thus promote SMBH growth. Furthermore, the natural variation in accretion properties as a function of environment is expected to be more considerable for a single AGN compared to variation of luminosity in galaxies, where scatter in star formation efficiency from individual star-forming regions is minimized from summing over a population of both star clusters and stars. The significant dependence of quasar properties on the local environment is consistent with results of high redshift (z > 7) hydrodynamical simulations. For example, di Matteo et al. (2017) identify the conditions of the local tidal field as instrumental to early black hole growth relative to the large-scale overdensity of the host halo.

3.2. Halo Mass Threshold for Feedback

As quasars in less massive hosts contribute significantly to the brightest end of the QLF for sufficient Σ, under these conditions and to zeroth order, the QLF becomes insensitive to the exact shape of the high-end for both the HMF and ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ (see the galaxy analog in Ren et al. 2019). From this, we can approximate the numerical deconvolution process to derive ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }}\gt 0)$ by scaling the ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}},{\rm{\Sigma }}=0)$ relation by a constant factor and by concurrently substituting the high-mass end past some characteristic halo mass ${M}_{{\rm{h}}}^{{\rm{c}}}$ with a constant luminosity value, i.e., ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}}\gt {M}_{{\rm{h}}}^{{\rm{c}}})\,={M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}}^{{\rm{c}}})$ . We find the best values for the scaling factor and ${M}_{{\rm{h}}}^{{\rm{c}}}$ by least-squares fits based on observed data points M < M*. The differences between the deconvolution method described in Allen et al. (2018) and this approximate method are marginal in terms of the resulting QLF (a detailed analysis of the differences between these two calculations are highlighted in Ren et al. 2019 for galaxies). However, one benefit from the approximate deconvolution method is that ${M}_{{\rm{h}}}^{{\rm{c}}}$ is a well-defined parameter that can be interpreted as a pseudo-scale for feedback.

In Figure 3, we show the extent of our modeled QLFs with this simple approximation method over a large range of scatter, 0.3 < Σ < 0.7 and look at two cases, using as inputs either the Akiyama et al. (2017) or Kulkarni et al. (2019) fits for the z ∼ 4 QLF. The modeled QLFs are broadly consistent with the data points compiled by the observations for all values of Σ considered here. In Figure 4, we determine the characteristic feedback threshold ${M}_{{\rm{h}}}^{{\rm{c}}}$ , in each of these cases given _DC and Σ as fixed parameters. We note that both Σ and _DC are fully degenerate with each other for the purpose of determining ${M}_{{\rm{h}}}^{{\rm{c}}}$ . Increasing Σ or decreasing _DC both lead to a lower ${M}_{{\rm{h}}}^{{\rm{c}}}$ . Still, while it is challenging to disentangle the two parameters, a characterization of the inherently constrained parameter space of _DC reveals useful information on the association between halos, galaxies and their AGNs.

**Figure 3.** Modeled QLFs at z ∼ 4 for a range of scatter, 0.3 < Σ < 0.7 (shaded area), assuming different observed QLFs as inputs (black data points). The left panel corresponds to an input using Akiyama et al. (2017) QLF and the right panel is for the Kulkarni et al. (2019) QLF. The dashed lines are the respective best-fit double power law to the data points (see Equation (3)). The parameters for the fits are as follows: (left panel, Akiyama et al. 2017) $\mathrm{log}({\phi }^{* })=-6.58$ , M* = −25.36, α = −1.3, and β = 3.11; (right panel, Kulkarni et al. 2019) log(ϕ^*) = −7.99, M* = −27.28, α = −2.11, and β = −4.64.
Download figure:
Standard image High-resolution image

**Figure 4.** Distribution of values for the critical halo mass threshold ${M}_{{\rm{h}}}^{c}$ as a function of input scatter Σ and duty cycle _DC. The left panel corresponds to models using the Akiyama et al. (2017) QLF, while the right panel is for the Kulkarni et al. (2019) QLF. The solid black lines are contours for select characteristic halo masses (M_⊙ log-scale).
Download figure:
Standard image High-resolution image

It is evident from Figure 4 that ${M}_{{\rm{h}}}^{{\rm{c}}}$ depends on the choice of the calibration QLF, and specifically on the robust determination of the QLF for the population around the characteristic M ∼ M* quasars. We investigate the magnitude of this effect and compare the distribution of ${M}_{{\rm{h}}}^{{\rm{c}}}$ thresholds across both cases. Quantitatively, the comparison of the two panels shows a small-to-modest difference, $\mathrm{Max}\ ({\rm{\Delta }}\mathrm{log}{M}_{{\rm{h}}}^{{\rm{c}}})\sim 0.3$ , in the distribution of ${M}_{{\rm{h}}}^{{\rm{c}}}$ at fixed Σ and _DC, with higher Σ and lower _DC values having the largest variation in ${M}_{{\rm{h}}}^{{\rm{c}}}$ between our choices of the calibration QLF.

A point of interest for this figure is that a number of studies have noted luminous quasars to preferentially reside in ∼10¹²M_⊙ hosts, coinciding with the halo mass range that has maximal specific star formation efficiency (e.g., Conroy & White 2012 and references therein). This idea is appealing as it would suggest a single origin for the joint regulation of the galaxy and black hole growth. However, in our model, ${M}_{{\rm{h}}}^{{\rm{c}}}\sim {10}^{12}{M}_{\odot }$ requires _DC ≲ 10⁻² for either choices of the QLF. Thus, if larger duty cycles are present at high z, the two processes would appear to be distinct and affected by feedback operating at different halo-mass scales. Furthermore, the situation can be more complex as we have assumed a monotonic relation between halo and quasar luminosity. In fact the numerical deconvolution method of Behroozi et al. (2010) is limited to finding the solution under this condition. It is also entirely possible that ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ experiences a turnover after SMBH growth becomes self-regulated at the highest masses. That would have minimal impact on the QLF for high Σ values as those rare sources in high-mass halos would be mixed with the population of more common halos, but it would make it increasingly difficult reconcile a self-consistent single halo-mass scale for quasar and galaxy feedback.

Deriving constraints on the average duty cycle _DC at high redshifts remain an open problem, as it critically depends on a variety of unsolved processes, including but not limited to the formation of seeds and the dominant mode of early SMBH growth. Clustering studies provide one avenue of partially constraining _DC; however, results are mixed due to the wide range of environments occupied by quasars, hence resulting in duty cycles between ${10}^{-3}\gtrsim {\epsilon }_{\mathrm{DC}}\gtrsim 6\times {10}^{-1}$ for z ∼ 4 quasars (Shen et al. 2007; He et al. 2017). On the other hand, Aversa et al. (2015) infer an average _DC ∼ 1 at z > 3 using a physically motivated light-curve parameterization to derive the BH mass function from the AGN luminosity function, which corresponds to a mass threshold ${M}_{{\rm{h}}}^{{\rm{c}}}\gt {10}^{12.9}{M}_{\odot }$ . Additionally, we can draw a comparison using ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ from the hydrodynamical suite MassiveBlack-II (MBII, Khandai et al. 2015) with _DC ∼ 1, finds a scatter of Σ ∼ 0.55 for M_h ∼ 10^11.8 M_⊙ halos. The MBII simulation volume of (100 Mpc h⁻¹)³ is capable of hosting halos up to M_h ∼ 10^12.7 M_⊙, but only expects ∼100 halos with M > 10¹²M_⊙. Thus, despite an insufficient simulation volume to fully capture the flattening in M_UV,c(M_h), the MBII analysis is consistent with a halo mass threshold of ${10}^{12.9}{M}_{\odot }\lesssim {M}_{{\rm{h}}}^{{\rm{c}}}\lesssim {10}^{13.2}{M}_{\odot }$ (see the Appendix for further details on modeling consistency with MB-II). From a physical interpretation perspective, assuming that AGN feedback can act independently in quenching galaxy and SMBH growth (Cielo et al. 2018) would naturally lead to a separation of the characteristic halo masses where star formation and black hole growth become affected.

3.3. Consequences on Luminosity-based Clustering Measurements

Rare quasars detected at the epoch of cosmic dawn are understood to be powered by SMBH of masses around 10⁹M_⊙, hence there is a general expectation that these quasars trace extreme overdensities and reside within massive hosts. However, attempts at observational confirmation based on quasar clustering have historically reported mixed results, in particular, at high redshift when quasar counts are increasingly sparse. For example, at z ∼ 4 high biases have been reported using quasar–quasar correlation function measurements (Shen et al. 2007; Onoue et al. 2017), while other studies have found little to no evidence of significant clustering based on quasar-galaxy cross correlation (Fukugita et al. 2004; He et al. 2017). Likewise, a number of studies (both simulations and observations) around z ∼ 5–6 quasars have shown that they belonged to a wide range of environments: overdensities (Stiavelli et al. 2005; Romano-Diaz et al. 2011; Husband et al. 2013; Costa et al. 2014; Morselli et al. 2014; García-Vergara et al. 2017) or average/underdensities (Kim et al. 2009; Mazzucchelli et al. 2017; Champagne et al. 2018; Ota et al. 2018). Compared to bright galaxy samples, that have been widely used to infer luminosity dependence of galaxy-halo properties (see Zheng et al. 2009; Trenti et al. 2012; Harikane et al. 2017), luminous quasars are rarer. Hence, small number count stochasticity intrinsically limits the ability to draw robust inference from current observations.

From Figure 1, it is evident that both _DC and Σ impact local clustering by reducing the average mass of the quasar host. To demonstrate this explicitly, we derive the distribution of the linear bias factor, b for our quasar host halos using:

$\begin{eqnarray}&&p(b| {M}_{\mathrm{UV}})=\displaystyle \frac{\tfrac{{dn}}{{{dM}}_{{\rm{h}}}}({M}_{{\rm{h}}}(b)){\rm{\Phi }}({M}_{{\rm{h}}}| {M}_{\mathrm{UV}})}{\phi ({M}_{\mathrm{UV}})},\end{eqnarray} \tag{ 4 }$

where M_h(b) is taken here as the inversion of the analytical Sheth & Tormen (1999) bias relation, and ${\rm{\Phi }}({M}_{{\rm{h}}}| {M}_{\mathrm{UV}})$ is the inverse CLF for the distribution of halo masses given some quasar magnitude, M_UV. In Figure 5 we show the range of linear bias for quasars in 2 groups: quasars that populate the extreme bright end (M_UV = −28) and the faint end (M_UV = −22). The sharp peak feature is a consequence of imposing a flat cutoff in ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ and we do not expect this to impact the results in any significant qualitative way. It is clear that Σ and _DC dictate the variety of environments for quasars more luminous than the characteristic magnitude. For the most luminous quasars, the dispersion in the bias is predominantly dependent on the threshold value, ${M}_{{\rm{h}}}^{{\rm{c}}}$ and the spread in M_h at a quasar magnitude to a lesser extent. This is because the bias factor has a strong nonlinear dependence on M_h. The nonlinearity also induces a positive skew that scales with Σ in the bias distribution, suggesting an uneven distribution of environments around the mode value of the bias. In contrast, the distribution of the linear bias around fainter quasars is seen to be substantially insensitive to Σ, but not to _DC. This highlights the opportunity to use the local clustering around fainter quasars as a probe to constrain the quasar duty cycle. Indeed, such a task is easily within reach using next-generation facilities such as the James Webb Space Telescope, which will be able to probe with both imaging and spectroscopy the fainter companions around high-z quasar halos. One caveat to note is that Figure 1 has been obtained using the median quasar magnitude versus halo mass relation inferred from the Akiyama et al. (2017) QLF. Using the relation derived from the Kulkarni et al. (2019) QLF would both shift to more positive values and broaden the linear bias distribution. Therefore, the analysis is dependent on a precise determination of the QLF.

**Figure 5.** Probability distribution of the quasar bias value for different model input parameters Σ (colored lines) and _DC (top vs. bottom panel). The halo mass corresponding to the bias value on the horizontal axis is shown on top of the upper panel. These distributions have been derived using the median quasar magnitude versus halo mass relations in our model with the Akiyama et al. (2017) QLF as inputs.
Download figure:
Standard image High-resolution image

In addition, it is important to stress that large values of Σ dilute signatures of luminosity dependent clustering, assuming a quasar duty cycle largely independent of halo mass. We check using the high/low luminosity bins of He et al. (2017), taking halo biases for quasars of M_UV ∼ −25.5 and M_UV ∼ −23.5 we find that our entire range of Σ > 0.3 models is consistent with no luminosity dependent clustering between the two bins. This is generally in agreement with the conclusion of He et al. (2017); however, we still find a weak luminosity dependence for clustering if we extend the baseline of our luminosity bins. In fact, we predict that the bias of the brightest quasars (M_UV ∼ −28) can still be quantified as different from that of the population of faint quasars (M_UV ∼ −22) to a confidence of 97.5% provided that Σ ∼ 0.5.

Future wide-area surveys are required to create a representative sample of brighter M_UV < −25.5 quasars in order to conclusively establish the luminosity dependence of quasar clustering.

4. General Remarks and Conclusion

The broad power-law bright end of the QLF relative to the exponential drop-off in the host HMF suggests there could be significant stochasticity, Σ, in the quasar magnitude versus halo mass relation, ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ . In this work, we use a CLF approach to derive ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ from the observed z ∼ 4 QLF assuming values for our free parameters scatter Σ and the quasar duty cycle _DC, and investigate how these parameters shape ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ . In addition to a full deconvolution study, we also construct an approximate best-fit ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ , with a functional form characterized by a constant value above a critical halo mass threshold at ${M}_{{\rm{h}}}^{{\rm{c}}}$ , i.e., ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}}\gt {M}_{{\rm{h}}}^{{\rm{c}}})={M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}}^{{\rm{c}}})$ . In this framework, ${M}_{{\rm{h}}}^{{\rm{c}}}$ can be interpreted as the critical value beyond which feedback is required to significantly regulate black hole growth/quasar radiation to avoid overproducing luminous objects that are not observed from the QLF. Finally, we investigated how this threshold depends on model assumptions and parameters, discussing physical interpretations of our model, with the following key results:

1.
We show that Σ induces a flattening in ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ to account for the abundance of lower mass quasar hosts with extreme accretion rates populating the bright end of the QLF (Figure 1). The flattening effect from stochasticity has been previously explored for a variety of relations (see the following examples: stellar mass, Behroozi et al. 2010; galaxy luminosity, Ren et al. 2019; BH mass/AGN luminosity, Aversa et al. 2015). We find that values of Σ < 0.3 lead to a rising quasar magnitude-halo mass ratio for more massive halos. This can be difficult to reconcile with ideas of self-regulating black hole growth. In contrast, Σ ≳ 0.5 implies a flattened ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ for massive halos, and is indicative of a turnover in quasar efficiency. This suggests both that the median black hole growth is regulated at the massive halo end and that there is an increasing likelihood that the most luminous quasars are extreme outliers in accretion efficiency hosted in relatively common medium-mass dark-matter halos.
2.
Following this, we note that constant duty cycle is a good approximation for modeling the bright end of the QLF assuming significant scatter (Σ ≳ 0.3). Since halo abundance is strongly dependent on mass, the abundance of lower mass halos with extreme accretion rates dominates over the influence of any changes in duty cycle for larger quasar hosts.
3.
The characteristic mass threshold for feedback (defined as the halo mass, ${M}_{{\rm{h}}}^{{\rm{c}}}$ where ${M}_{\mathrm{UV},{\rm{c}}}({M}_{{\rm{h}}})$ flattens) is strongly dependent on Σ and _DC (Figure 4). ${M}_{{\rm{h}}}^{{\rm{c}}}$ is relatively insensitive to variations in observed QLF determinations, with Max( ${\rm{\Delta }}\mathrm{log}{M}_{{\rm{h}}}^{{\rm{c}}})\sim 0.3$ between the use of the Akiyama et al. (2017) and of the Kulkarni et al. (2019) QLF as modeling inputs.
4.
Σ and _DC are strongly degenerate. An increase in Σ has essentially the same impact in decreasing ${M}_{{\rm{h}}}^{{\rm{c}}}$ as a reduction in _DC. Disentangling this degeneracy would require additional measurements, such as local clustering using cross-correlations between quasars and galaxies around faint (M_UV ∼ −22) quasars (Figure 5). These complementary observations would need to reach about a factor 10× fainter than the quasar (−20 ≲ M_UV ≲ −19). This limit is already within imaging capabilities of current facilities (e.g., Wide Field Camera 3 instrument on the Hubble Space Telescope), and sufficiently bright for spectroscopy with the upcoming James Webb Space Telescope.
5.
Matching the halo mass for optimal efficiency in both quasar radiation emission and stellar formation in galaxies, ∼10¹²M_⊙, requires _DC < 10⁻². Observations of clustering and hydrodynamical simulations infer a range from 10⁻³ < _DC < 1. A high _DC ∼ 1 would return ${M}_{{\rm{h}}}^{{\rm{c}}}\gtrsim {10}^{12.9}{M}_{\odot }$ for quasars hosts, which could suggest that AGN feedback quenches galaxies and SMBH growth independently of each other (Cielo et al. 2018).
6.
We calculate the distribution of halo bias around bright (M_UV = −28) and faint (M_UV = −22) quasar hosts (Figure 5). Σ increases the spread of the biases, while a rising _DC increases the median halo bias. Σ > 0.3 weakens the luminosity dependence for clustering such that there is effectively no luminosity dependence for clustering (He et al. 2017).

The framework developed in this manuscript allows us to speculate on the evolution of quasar demographics across z. For example, if we make the assumption that the magnitude of scatter, Σ and the halo mass threshold, ${M}_{{\rm{h}}}^{{\rm{c}}}$ remain relatively unchanged across redshift, then we would generally expect the bright end slope of the QLF to become shallower at higher z from the decreased abundance of $M\gt {M}_{{\rm{h}}}^{{\rm{c}}}$ halos. However, the change in Σ over z is not well constrained observationally owing to a number of factors, such as the intrinsic rarity of the brightest objects and the challenges in observing the more typical M ∼ M* objects. A recent work by Marshall et al. (2020) investigating the primary growth mechanisms of SMBHs finds that merger-driven SMBH growth is subdominant compared to instability-driven growth at z > 2. The smaller contribution of mergers on the mass history of SMBHs can point to Σ having a weaker dependence on z. The assumption that ${M}_{{\rm{h}}}^{{\rm{c}}}$ is independent over z can be justified on a theoretical basis from the relative z independence in halo-mass where we expect radio-mode feedback to become nonnegligible (Croton et al. 2006). Additionally, the mass-independence from this mode of feedback is also partially supported through sophisticated empirical modeling of stellar-mass to halo-mass ratio (SHMR) demonstrating a weak evolution in the peak of the SHMR relation across redshifts (Tacchella et al. 2018; Behroozi et al. 2019).

However, it should be noted that the current extent of available evidence through observations does not yield a compelling case on the direction the bright end slope should evolve with redshift. In one case, the QLFs from Akiyama et al. (2017) and Matsuoka et al. (2018) show that the bright end slope evolves to be more shallow at higher z (at least from z > 4), while empirical modeling by Kulkarni et al. (2019) suggests the contrary, that the bright end slope should steepen toward higher z. It is worthwhile to note that the bolometric QLFs recently recompiled by Shen et al. (2020) suggest that the bright end slope becomes shallower at higher z (from z ∼ 2), which is broadly consistent with the picture provided here. However, one surprising element is that Shen et al. (2020) finds the bright end slope to also become shallower at z < 2. A conclusive observational picture of the evolution in the QLF across a large range of magnitudes is therefore essential in order to answer such questions. This task is both challenging and time consuming, requiring the need to leverage both: (1) wide-field surveys to capture abundance of the rarest and brightest of quasars, and (2) deep imaging plus spectroscopy together with extensive modeling to deconvolve the contribution of quasar light within its host galaxy for the characteristic M ∼ M* quasars. In this context, simple but effective models that physically capture the dominant contribution in the evolution of the QLF can provide effective tools to assess expectations in preparation of future surveys.

This research was conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. K.R. is additionally supported through the Research Training Program Scholarship from the Australian Government. T.D.M. acknowledges funding from NSF ACI-1614853, NSF AST-1517593, NSF AST-1616168, NASA ATP 80NSSC18K1015, and NASA ATP 17-0123.

Appendix: Comparison of Modeling Methods and Simulations

In Section 3.2, we determine a range for the halo-mass threshold in MassiveBlack-II (MB-II). In Figure A1, we show the quasar luminosity to halo mass ratios from individual z = 4 quasars in MB-II. It is clear that MB-II lacks the volume to capture any possible turnover in M_h. The output displayed is generated from the deconvolution process outlined in Section 2.2. The input uses Σ = 0.55 corresponding to the scatter in MB-II quasars inside M_h ∼ 10^11.8 M_⊙ halos and assumes the Kulkarni et al. (2019) QLF, which is a closer fit to faint quasars present in the MB-II QLF. The modeling shown here indicates a turnover at M_h ∼ 10^13.3 M_⊙, slightly higher than the one inferred through our approximate modeling at M_h ∼ 10^13.2 M_⊙. As discussed in Section 3.2, the key parameter between the choice of the QLF that impacts ${M}_{{\rm{h}}}^{{\rm{c}}}$ is the position and normalization of the characteristic magnitude break, M^*, hence only a small difference in ${M}_{{\rm{h}}}^{{\rm{c}}}$ is expected when using either Akiyama et al. (2017) and Kulkarni et al. (2019) QLFs.

Stochastic Processes as the Origin of the Double Power-law Shape of the Quasar Luminosity Function

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction