The Bolometric Quasar Luminosity Function at z = 0-7

Introduction

Luminous quasars and active galactic nuclei (AGN) in general We use the phrase ”quasar” across the paper. We are not just referring to the optically bright and unobscured systems but the entire AGN population. are observable manifestations of accreting supermassive black holes (SMBHs) at galaxy centers. Gas accreted onto the SMBH forms an accretion disk from which thermal emission is generated through dissipative processes (Shakura & Sunyaev 1973; Rees 1984, e.g.,). Due to their high radiative efficiency, such objects can be extremely luminous and are detected at $z>7$ (Mortlock et al. 2011; Venemans et al. 2015; Bañados et al. 2018). The evolution of quasars is crucial to understand the formation and evolution of SMBHs in the Universe. Apart from that, quasars are one of the most important radiation sources in the Universe. They are luminous in almost all accessible bands and their radiation has a significant impact in the Universe. For example, quasar emission is important for the build-up of cosmic infrared (IR) and X-ray radiation backgrounds. Quasar emission in the extreme ultraviolet (UV) is believed to dominate the reionization of helium in the Universe and may have a non-negligible contribution to the reionization of hydrogen, although star-forming galaxies dominate hydrogen reionization in most current models (Faucher-Giguère et al. 2008a; Faucher-Giguère et al. 2008b; Faucher-Giguère et al. 2009; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2015; Haardt & Salvaterra 2015; Giallongo et al. 2015; Onoue et al. 2017; Parsa et al. 2018, e.g.,). Furthermore, observations have demonstrated that galaxies and SMBHs co-evolve (see reviews of Alexander & Hickox 2012; Fabian 2012; Kormendy & Ho 2013; Heckman & Best 2014, and references therein). For example, the masses of the SMBHs are correlated with the masses, luminosities and velocity dispersions of their host galaxy spheroids (Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Gültekin et al. 2009, e.g.,). AGN are also widely believed to impact star formation in their host galaxies via a "feedback" mechanism that helps quench galaxies (Sanders & Mirabel 1996; Springel et al. 2005; Bower et al. 2006; Croton et al. 2006; Sijacki et al. 2007; Somerville et al. 2008; Hopkins et al. 2008; Feruglio et al. 2010; Fabian 2012; Cicone et al. 2014, e.g.,) and solve the classical "cooling flow" problem (Cowie & Binney 1977; Fabian & Nulsen 1977; Fabian et al. 1984; Tabor & Binney 1993; Fabian 1994; Croton et al. 2006, e.g.,). Therefore, studying the evolution of quasar populations along cosmic time is of great importance in cosmology and galaxy formation.

The quasar luminosity function (QLF), which is the comoving number density of quasars as a function of luminosity, is perhaps the most important observational signature of quasar populations. The study of the QLF goes back decades in the rest-frame optical/UV (Schmidt 1968; Schmidt & Green 1983; Koo & Kron 1988; Boyle et al. 1988; Hartwick & Schade 1990; Hewett et al. 1993; Warren et al. 1994; Schmidt et al. 1995a; Kennefick et al. 1995; Pei 1995; Boyle et al. 2000; Fan et al. 2001c; Fan et al. 2004; Richards et al. 2006a; Croom et al. 2009; Willott et al. 2010; Glikman et al. 2011; Ross et al. 2013; McGreer et al. 2013; Kashikawa et al. 2015; Jiang et al. 2016, e.g.,), soft X-ray (Maccacaro et al. 1991; Boyle et al. 1993; Jones et al. 1997; Page et al. 1997; Miyaji et al. 2000; Hasinger et al. 2005, e.g.,), hard X-ray (Ueda et al. 2003; La Franca et al. 2005; Barger et al. 2005; Silverman et al. 2008; Ebrero et al. 2009; Yencho et al. 2009; Aird et al. 2010; Ueda et al. 2014; Aird et al. 2015a, e.g.,) and IR (Brown et al. 2006; Matute et al. 2006; Assef et al. 2011; Lacy et al. 2015, e.g.,). These studies have conclusively shown that the observed QLF exhibits a strong redshift evolution. This is not simply an evolution in the normalization (number density) but also in the slope of the QLF. For instance, the number density of low luminosity AGN peaks at lower redshift than that of bright quasars indicating the "cosmic downsizing" of AGN (Cowie et al. 1996; Barger et al. 2005; Hasinger et al. 2005, e.g.,). AGN feedback that shuts down the supply of gas for accretion may be responsible for this phenomenon. Both optical and X-ray studies have argued that the faint-end slope of the QLF gets steeper from $z=2$ to $z=0$ (Aird et al. 2015a; Kulkarni et al. 2018, e.g.,). These investigations of the QLF have also found that both the typical spectral shape (Wilkes et al. 1994; Green et al. 1995; Vignali et al. 2003; Strateva et al. 2005; Richards et al. 2006b; Steffen et al. 2006; Just et al. 2007; Lusso et al. 2010; Kashikawa et al. 2015; Lusso & Risaliti 2016, e.g.,) and the obscuring column density distribution of quasars (Hill et al. 1996; Simpson et al. 1999; Willott et al. 2000; Steffen et al. 2003; Ueda et al. 2003; Grimes et al. 2004; Sazonov & Revnivtsev 2004; Barger et al. 2005; Hao et al. 2005; Ueda et al. 2014, e.g.,) have a dependence on quasar luminosity. For example, fainter quasars tend to be more obscured and their emission is more dominated by the X-rays.

In the last decade, the redshift frontier of the observations of quasars have been pushed up to $z>7$ (Mortlock et al. 2011; Bañados et al. 2018; Wang et al. 2018b) and about $40$ quasars are now known at $z\gtrsim 6.5$ (Willott et al. 2010; Venemans et al. 2013; Venemans et al. 2015; Jiang et al. 2016; Reed et al. 2017; Mazzucchelli et al. 2017; Matsuoka et al. 2018; Ross & Cross 2019, e.g.,). These quasars reveal the early growth of SMBHs and also pinpoint the locations for the assembly of massive galaxies in the early Universe. The absorption spectra of these high redshift quasars are important to study the reionization history of the Universe (Miralda-Escudé 1998; Madau & Rees 2000; Fan et al. 2002; Fan et al. 2006, e.g.,). However, due to the rapid decline in the quasar number density at high redshift, detecting quasars and constraining the QLF is currently very difficult at $z\gtrsim 6$. The next generation deep, wide-field infrared surveys will help push the detection of quasars to $z\simeq 9-10$ and deeper optical/UV surveys will provide better constraints on the faint end of the QLF.

Interpreting the observational findings, however, is complicated by the fact that observations in a single band are always subject to selection effects, host galaxy contamination and reddening and obscuration all in a complicated, wavelength-dependent manner. Although quasars are intrinsically very luminous in the optical/UV, dust extinction along some viewing angles (Antonucci 1993; Urry & Padovani 1995, e.g.,) can make quasars much more difficult to detect. Heavily obscured AGN can easily be contaminated with the UV stellar light from their host galaxies (Hickox & Alexander 2018, e.g., see review of). Even in the X-ray, which is much less affected by dust, the Compton-thick (CTK) AGN, which account for $20\%-50\%$ (Burlon et al. 2011; Ricci et al. 2015, e.g.,) of the total AGN population, are still severely blocked and current observations remain largely incomplete. In the mid-IR, due to the strong absorption in the terrestrial atmosphere, observations are more limited and also can be contaminated by the hot dust emission in star forming galaxies. In far-IR to millimeter wavelengths ($30\micron-10\,{\rm mm}$), the majority of AGN are contaminated by emission from dust heated by star formation in host galaxies, which limits the effectiveness of AGN identification. Furthermore, measurements of the QLF based on a single survey are limited in their luminosity coverage and volume probed and are subjected to various biases and uncertainties in completeness corrections.

Given these limitations, what physical models for AGN demographics, SMBH growth and AGN feedback, really require is the bolometric QLF over all redshifts. The bolometric quasar luminosity is the quantity tightly related to the accretion rate of the SMBH and is the ideal quantity to study the physical evolution of quasars. Hopkins et al. 2007 developed a bolometric QLF model that simultaneously fitted the accessible measurements at the time, in different bands. The model has been widely used but has several important shortcomings: First, the model was poorly constrained at $z\gtrsim 3$ due to limited available data at the time and has been shown to deviate significantly from recent observations. Second, the integrated bolometric luminosity at the bright end predicted by this model actually diverges when extrapolating to high redshift ($z\sim 7-8$). Third, the number density normalization of the QLF was assumed to be a constant over redshifts, which does not agree with newer observations at high redshift.

In this paper, we provide a new model for the bolometric QLF at $z=0-7$ constrained by emerging observations of the QLF in the optical, UV, IR and X-ray in the last decade. The paper is organized as follows: In Section 2, we introduce our observational data compilation. In Section 3, we introduce our model linking the observed QLFs with the bolometric QLF. The model includes new bolometric and extinction corrections. In Section 4, we perform a fit to the data and constrain the bolomeric QLF. In Section 5, the evolution of the bolometric QLF is analyzed. In Section 6, we present several predictions from our best-fit bolometric QLF model and demonstrate its consistency with observations from independent channels.

We employ the following cosmological parameters: $\Omega_{\rm m}=0.30$, $\Omega_{\Lambda}=0.70$, $H_{0}=100h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}=70\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. The code of all the analysis in this paper along with the observational data compiled are publicly available (see Appendix C for details).

Observational Data Sets

In this section, we briefly introduce the observations compiled in this work and emphasize the corrections adopted. A full list of the observations compiled is shown in Table 5. We note that some observations used overlapping quasar samples in their binned estimations and are thus not fully independent. We do not include older observations if all the quasar samples used there were covered by later work. For all the observational data, we correct all relevant quantities (distances, luminosities, volumes) to be consistent with our adopted cosmological parameters.

We define "optical" wavelengths as $2500\text{\AA}\leq\lambda\leq 1\micron$ and "UV" wavelengths as $600\text{\AA}\leq\lambda\leq 2500\text{\AA}$ The quasar SED at rest-frame $50\text{\AA}\leq\lambda\leq 600\text{\AA}$ is almost inaccessible in optical/UV observations due to strong extinction at these wavelengths. In the construction of our SED model in Section 3.1, we directly connect the $600\text{\AA}$ flux with the X-ray SED.. We unify the luminosities measured in rest-frame optical (UV) wavelengths in observations to the B band (UV) luminosity defined in Table 1. The optical/UV QLF observations compiled in this work are largely based on the observations listed in Hopkins et al. 2007; Giallongo et al. 2012; Manti et al. 2017 and Kulkarni et al. 2018 (along with their QLF data shared online https://github.com/gkulkarni/QLF/blob/master/Data/allqlfs.dat). The observational compilation from Kulkarni et al. 2018 includes: Bongiorno et al. 2007; Siana et al. 2008; Jiang et al. 2009; Willott et al. 2010; Glikman et al. 2011; Masters et al. 2012; Palanque-Delabrouille et al. 2013; Ross et al. 2013; McGreer et al. 2013; Kashikawa et al. 2015. In the Kulkarni et al. 2018 compilation, the poisson errors in several works were recomputed using the Gehrels 1986 formula. The K-corrections have been unified to that in Lusso et al. 2015, which is based on the stacked spectra of $53$ quasars observed at $z\sim 2.4$. In fact, the uncertainty in K-corrections owing to different spectral assumptions was estimated to be within $0.2\,{\rm mag}$ (Lusso et al. 2015) which is smaller than the uncertainties of the binned estimation itself. The uncertainties in conversion factors between luminosities of different rest-frame bands were also estimated to be smaller than other sources of errors. Other specific corrections have been made in the Kulkarni et al. 2018 compilation are: (1) bins with severe incompleteness from Ross et al. 2013 were discarded; (2) binned estimations in Palanque-Delabrouille et al. 2013 at $z>2.6$ were discarded since Lyman-alpha forest enters g band for those redshifts; (3) data from McGreer et al. 2013 was restricted to $M_{\rm 1450}>-26.73$ to avoid overlapping with Yang et al. 2016; (4) for Willott et al. 2010 and Kashikawa et al. 2015, the redshift intervals were recomputed using consistent completeness estimations.

Outside the Kulkarni et al. 2018 compilation, we include measurements from Fontanot et al. 2007; Croom et al. 2009; Shen & Kelly 2012; Jiang et al. 2016; Palanque-Delabrouille et al. 2016; Yang et al. 2016; Akiyama et al. 2018; Matsuoka et al. 2018; McGreer et al. 2018; Wang et al. 2018a; Yang et al. 2018. Observed optical band luminosities are all converted to UV luminosity either with corrections made in these papers or with the formula in Ross et al. 2013 if no corrections had already been made. Matsuoka et al. 2018 and Yang et al. 2018 have binned estimations that correspond to only one object in the bin, which were interpreted as upper limits there. However, their Poisson error estimations were not correct and we recalculate the Poisson errors using the table in Gehrels 1986. After the correction, these data points have proper upper and lower limits and can be included into our standard fitting procedure. For the observations compiled in Hopkins et al. 2007 (Kennefick et al. 1995; Schmidt et al. 1995a; Fan et al. 2001a; Fan et al. 2001b; Fan et al. 2003; Fan et al. 2004; Wolf et al. 2003; Cristiani et al. 2004; Croom et al. 2004; Hunt et al. 2004; Richards et al. 2005; Richards et al. 2006b; Siana et al. 2006), we include only those whose quasar samples are not completely covered by the more recent work discussed above. The details of all the observations compiled in this paper are listed in Table 5.

2 X-ray

We define "X-ray" wavelengths as $\lambda\leq 50\text{\AA}$ ($E\gtrsim 0.25\,{\rm keV}$) which covers the typical soft X-ray and hard X-ray bands defined in Table 1. In the X-ray, in addition to the observations compiled in Hopkins et al. 2007 (Miyaji et al. 2000; Miyaji et al. 2001; Ueda et al. 2003; Sazonov & Revnivtsev 2004; Barger et al. 2005; La Franca et al. 2005; Hasinger et al. 2005; Nandra et al. 2005; Silverman et al. 2005), we include new observational data from Ebrero et al. 2009; Aird et al. 2008; Silverman et al. 2008; Yencho et al. 2009; Aird et al. 2010; Fiore et al. 2012; Ueda et al. 2014; Aird et al. 2015a; Aird et al. 2015b; Miyaji et al. 2015; Khorunzhev et al. 2018. Among them, Aird et al. 2008 is an update based on Nandra et al. 2005 and Silverman et al. 2008 is an extension to Silverman et al. 2005. Aird et al. 2015a and Ueda et al. 2014 derived binned estimation of the hard X-ray luminosity functions separately based on soft or hard X-ray selected samples. We include both of them in our compilation. Aird et al. 2015b is an observation of the $10-40\,{\rm keV}$ X-ray luminosity function. The luminosities are converted to the hard X-ray luminosities with our SED model which will be discussed in the following section. Some observational works (Ebrero et al. 2009; Ueda et al. 2014; Aird et al. 2015b; Miyaji et al. 2015) have done their own "absorption" corrections and presented the "de-absorbed" compton thin QLFs. This would potentially generate double-counting of the extinction effects since we also intend to do extinction corrections in our model. We address this by reintroducing the extinction effect (only in the compton thin regime) for these data points using our extinction model which will be discussed in Section 3.2.

3 Infrared (IR)

We define "IR" wavelengths as $\lambda\geq 1\micron$. We unify the luminosities measured in rest-frame IR wavelengths to the mid-IR luminosity defined in Table 1. In the IR, in addition to the observations compiled in Hopkins et al. 2007 (Brown et al. 2006; Matute et al. 2006), we include new observations from Assef et al. 2011 and Lacy et al. 2015. The luminosities are converted to the mid-IR ($15\micron$) luminosity with our SED model. These observations have extended the redshift coverage of the IR QLF up to $z=5.8$. However, there is still an apparent deficiency in IR observations compared with other wavelengths. Deep and large field IR surveys are an urgent need in the study of the QLF at high redshift. Though the total number of IR data points are limited and thus they have low statistical significance in the fit of the bolometric QLF, they do provide an independent check for our bolometric QLF model.

Model

In this section, we construct the mean SED model for quasars. With the mean SED, we will calculate the bolometric corrections for the rest-frame B band, UV, soft & hard X-ray and mid-IR, respectively.

In the optical/UV, we start with the SED template in Krawczyk et al. 2013, which was based on $108184$ luminous broad-lined quasars observed at $0.064<z<5.46$. Among these sources, $11468$ showing sign of dust reddening ($\Delta(g-i)>0.3$) had been discarded by Krawczyk et al. 2013 in deriving the mean SED template. Therefore, this SED template can be considered not strongly affected by reddening and obscuration. The extinction corrections on the quasar luminosities will be considered separately in the next section. This SED template starts at $\sim 30\micron$ and truncates at $912\text{\AA}$. We extend the SED to the extreme UV (here defined as $\lambda<912\text{\AA}$) using the power-law model $f_{\nu}=\nu^{\alpha_{\nu}}$ with index $\alpha_{\nu}=-1.70$ reported by Lusso et al. 2015. We truncate this extension at $600\text{\AA}$ where Lusso et al. 2015’s measurement ended and directly connect the flux at $600\text{\AA}$ with the X-ray template which will be discussed then.

Historically, the optical/UV SED was often modelled as a power-law $f_{\nu}=\nu^{\alpha_{\nu}}$. In the UV, Vanden Berk et al. 2001 found that the $1300\text{\AA}$ to $5000\text{\AA}$ continuum roughly has a power-law index $\alpha_{\nu}=-0.44\pm 0.10$. Telfer et al. 2002 found $\alpha_{\nu}=-0.69\pm 0.06$ at $1200\text{\AA}\lesssim\lambda\leq 2200\text{\AA}$. Shull et al. 2012 found $\alpha_{\nu}=-0.68\pm 0.14$ at $1200\text{\AA}\leq\lambda\leq 2000\text{\AA}$. Lusso et al. 2015 found $\alpha_{\nu}=-0.61\pm 0.01$ at $912\text{\AA}\leq\lambda\leq 2500\text{\AA}$. The differences between Vanden Berk et al. 2001 and other updated measurements arise from different continuum regions used to measure the slope. In the extreme UV, Telfer et al. 2002 found $\alpha_{\nu}=-1.76\pm 0.12$ at $500\text{\AA}\lesssim\lambda\leq 1200\text{\AA}$. Scott et al. 2004 found $\alpha_{\nu}=-0.56^{+0.38}_{-0.28}$ at $630\text{\AA}\lesssim\lambda\leq 1155\text{\AA}$. Lusso et al. 2015 found $\alpha_{\nu}=-1.70\pm 0.61$ at $\sim 600\text{\AA}\leq\lambda\leq 912\text{\AA}$. The update of break point from $\sim 1200\text{\AA}$ to $\sim 912\text{\AA}$ mainly attributes to more careful correction on IGM absorption (Lusso et al. 2015). We do not consider the potential redshift/luminosity dependence of the break point, since it has almost no influence on the bolometric corrections. In Figure 1, we show that our optical/UV SED template is generally consistent with the most recent power-law models.

1.2 IR

In the IR, we adopt the SED template in Krawczyk et al. 2013. We extend the template in the long wavelength end to $100\micron$ using the Richards et al. 2006b SED which behaves almost the same as the Krawczyk et al. 2013 SED at $\lambda>10\micron$. We note that this IR SED has already included dust emission. No additional dust emission model will be required.

1.3 X-ray

The X-ray SED template is generated with a cut-off power-law model $f(E)\sim E^{1-\Gamma}\exp(-E/E_{\rm c})$ with the photon index $\Gamma=1.9$ and the cut-off energy $E_{\rm c}=300\,{\rm keV}$ (Dadina 2008; Ueda et al. 2014; Aird et al. 2015a, e.g.,). An additional reflection component is added using the pexrav model (Magdziarz & Zdziarski 1995) assuming the reflection relative strength $R=1$, the inclination angle $i=60^{\circ}$ and solar abundances. Then, we have to properly normalize the X-ray SED relative to the optical SED. Previous studies have reported a correlation between $L_{\nu}(2\,{\rm keV})$ and $L_{\nu}(2500$Å) (the unit of $L_{\nu}$ is $\,{\rm erg}\,{\rm s}^{-1}\,{\rm Hz}^{-1}$):

where $\beta$ is found to be $0.7-0.8$ suggesting a non-linear correlation between the X-ray and optical luminosities. Defining $\alpha_{\rm ox}$ as:

where $A=0.384\,(1-\beta)$ and $C^{\prime}=0.384C$. These prefactors have been measured through observations. However, since there is scatter in this relation, treating $L_{\nu}(2500\text{\AA})$ or $L_{\nu}(2\,{\rm keV})$ as the independent variable will lead to different results if quasars are not perfectly selected in observations. The bisector of the two fitted relation treating either $L_{\nu}(2500\text{\AA})$ or $L_{\nu}(2\,{\rm keV})$ as the independent variable is usually adopted. For example, Steffen et al. 2006 measured $\beta=0.721\pm 0.011$ and $C=4.531\pm 0.688$; Just et al. 2007 measured $\beta=0.709\pm 0.010$ and $C=4.822\pm 0.627$; Lusso et al. 2010 measured $\beta=0.760\pm 0.022$ and $C=3.508\pm 0.641$. Young et al. 2010; Xu 2011; Lusso & Risaliti 2016 found consistent results with previous works though they treated $L_{\nu}(2500\text{\AA})$ as the independent variable. Dependence of $\alpha_{\rm ox}$ on redshift had been reported in Bechtold et al. 2003, but was not confirmed in the following studies. Given these observational results, we conclude that the relation constrained by Steffen et al. 2006, which was adopted in Hopkins et al. 2007, is still consistent with updated observations. We continue to use the parameters measured by Steffen et al. 2006 though varying the parameter choices does not have a significant influence on the bolometric corrections. The X-ray SED is then scaled with the $\alpha_{\rm ox}$ with respect to the optical SED.

1.4 Bolometric corrections

The direct product of our quasar SED model is the bolometric correction, defined as the ratio between the bolometric luminosity, $L_{\rm bol}$, and the observed luminosity in a certain band, $L_{\rm band}$. The definitions of the luminosities in the bands are presented in Table 1. The bolometric luminosity is defined as the integrated luminosity from $30\micron$ to $500\,{\rm keV}$, which represents all the energy budget generated by the accretion of the SMBH. Some studies included the emission beyond $30\micron$ in the bolometric luminosity. But we find that extending the long-wavelength bound to $100\micron$ will only lead to $<0.02\,{\rm dex}$ difference in the bolometric luminosity. Some studies (Marconi et al. 2004; Krawczyk et al. 2013, e.g.,) have discussed that the reprocessed emission in the IR and $>2\,{\rm keV}$ X-ray should be excluded in determining the bolometric luminosity to avoid potential double-counting of quasars’ intrinsic emission. We have tested that using $1\micron$ to $2\,{\rm keV}$ as the range for integration will systematically decrease the bolometric luminosity by $\sim 0.2\,{\rm dex}$.

However, quasars do not have a single universal SED. There are real variations in the spectral shape, which translate to scatters in the bolometric corrections and influence the observed QLFs in the bands. To evaluate this, we first create an ensemble of SEDs. The configuration of these SEDs are similar to our fiducial SED: in the IR, we adopt our fiducial SED; in the optical/UV, for simplicity, we adopt a broken power-law with the break point at $912\text{\AA}$, with a fixed slope $-1.70$ at $\lambda<912\text{\AA}$ and a free slope $\alpha_{\rm opt}$ at $\lambda>912\text{\AA}$; in the X-ray, we adopt our fiducial X-ray SED model but with a free photon index $\Gamma$; the optical/UV and X-ray SEDs are connected with a free $\alpha_{\rm ox}$. We generate an ensemble of $10^{5}$ SEDs with randomly sampled $L_{\nu}(2500\text{\AA})$, $\alpha_{\rm opt}$, $\Gamma$ and $\alpha_{\rm ox}$. In sampling $\alpha_{\rm opt}$, $\Gamma$ and $\alpha_{\rm ox}$, we adopt a normal distribution around median value with a constant scatter. We adopt $\overline{\Gamma}\pm\sigma_{\Gamma}=1.9\pm 0.2$ (Ueda et al. 2014; Aird et al. 2015a, e.g.,), $\overline{\alpha_{\rm opt}}\pm\sigma_{\rm opt}=-0.44\pm 0.125$ (Vanden Berk et al. 2001; Richards et al. 2003), $\sigma_{\rm ox}\simeq 0.1$ (Steffen et al. 2006; Lusso et al. 2010, e.g.,). The bolometric luminosity and bolometric corrections for each realization of the SED are calculated. Then we divide the SEDs based on their bolometric luminosities into $30$ uniformly log-spaced bins from $10^{38}$ to $10^{48}\,\,{\rm erg}\,{\rm s}^{-1}$. We evaluate the standard deviation of the bolometric correction of each band in each bolometric luminosity bin, shown in the bottom panel of Figure 2. Double plateaus show up at the bright and faint ends where a certain band is dominant or negligible in the bolometric luminosity. In Hopkins et al. 2007, the dispersion of the bolometric corrections were fitted with: $\sigma_{\rm corr}(L_{\rm bol})=\sigma_{1}(L_{\rm bol}/10^{9}{\,\rm L_{\odot}})^{\beta}+\sigma_{2}$. However, we find this formula no longer appropriate to fit our results, so we fit the dispersion with an error function:

which naturally exhibits a double plateau shape. The best-fit parameters are listed in Table 1. The fitted relations are also shown in the bottom panel of Figure 2. These results indicate a $\sim 0.1\,{\rm dex}$ uncorrelated dispersion in quasar SEDs that is consistent with observations.

In the top panel of Figure 2, we show the bolometric corrections as a function of bolometric luminosity for all bands along with their dispersions shown with shaded regions. The bolometric corrections are generally similar to the Hopkins et al. 2007 model except for the differences at the faint end driven by the updates in the X-ray SED. Following Hopkins et al. 2007, we fit the dependence of the bolometric corrections on bolometric luminosity with a double power-law:

The best-fit parameters are listed in Table 1.

We note that the derivation of the optical/UV and X-ray luminosities using these bolometric corrections has not considered extinction yet. The observed luminosities will be further affected by extinction, which will be discussed in the following section.

2 Dust and gas extinction

The absorption and scattering of surrounding gas and dust further modifies the intrinsic emission of quasars. Neutral hydrogen photoelectric absorption is crucial to the extinction in the X-ray while dust is crucial to the extinction in the optical/UV. Here, we first introduce the neutral hydrogen column density ($N_{\rm H}$) distribution model which determines the extinction in the X-ray. Then, $N_{\rm H}$ is converted to the column density of dust assuming a dust-to-gas ratio. The dust abundance determines the extinction in the optical/UV.

In Hopkins et al. 2007, where the constant $N_{\rm H}$ model was shown to fail, the $N_{\rm H}$ distribution model from Ueda et al. 2003 was adopted as the fiducial model. Here, we update the $N_{\rm H}$ distribution with the results from Ueda et al. 2014, which was based on measurements of $N_{\rm H}$ and the intrinsic hard X-ray luminosity for each individual object in their sample. The model provides the probability distribution of $N_{\rm H}$, $f(L_{\rm X},z;N_{\rm H})$, at a given intrinsic hard X-ray luminosity (denoted as $L_{\rm X}$) and a redshift. $f(L_{\rm X},z;N_{\rm H})$ is normalized in the compton thin (CTN, $\log{N_{\rm H}}\leq 24$) regime:

where the unit of $N_{\rm H}$ is assumed to be $\,{\rm cm}^{-2}$ and the lower limit of $\log{N_{\rm H}}=20$ is a dummy value introduced for convenience and Ueda et al. 2014 has assigned $\log{N_{\rm H}}=20$ for all the quasars with $\log{N_{\rm H}}<20$.

$f(L_{\rm X},z;N_{\rm H})$ is characterized by three parameters: $\psi(L_{\rm X},z)$, the fraction of absorbed quasars ($22\leq\log{N_{\rm H}}\leq 24$) in total CTN quasars; $f_{\rm CTK}$, the fraction of compton thick (CTK, $\log{N_{\rm H}}\geq 24$) quasars relative to the fraction of absorbed CTN quasars; $\epsilon$, the ratio of the quasars with $23\leq\log{N_{\rm H}}\leq 24$ to those with $22\leq\log{N_{\rm H}}\leq 23$. This $N_{\rm H}$ distribution can then be written as (Ueda et al. 2014):

when $\psi(L_{\rm X},z)<\dfrac{1+\epsilon}{3+\epsilon}$ and:

when $\psi(L_{\rm X},z)\geq\dfrac{1+\epsilon}{3+\epsilon}$. The model assumes $\epsilon=1.7,\,f_{\rm CTK}=1$ and:

where $\psi_{\rm min}=0.2$, $\psi_{\rm max}=0.84$, $\psi_{\rm 43.75}(z)$ depends on redshift as:

The model describes a negative dependence of the absorbed quasar fraction on the intrinsic quasar hard X-ray luminosity as well as redshift at $z<2$.

Given this $N_{\rm H}$ distribution model, both the absorbed and the CTK quasar fractions decrease at higher hard X-ray luminosities and increase at higher redshift with a plateau at $z\geq 2$. The studies of the QLF and extinction properties in the X-ray have many variations in the data used, fitting methods, assumptions of the spectrum form and $N_{\rm H}$ distribution function form, treatments of redshift uncertainties and sources without counterparts. Therefore, it is worth comparing our fiducial extinction model with models determined in other works. In the top and middle panels of Figure 3, we compare the predictions on the absorbed quasar fraction and the CTK quasar fraction from this model with other observational constraints (Ueda et al. 2003; Burlon et al. 2011; Brightman & Ueda 2012; Merloni et al. 2014; Aird et al. 2015a; Buchner et al. 2015; Ricci et al. 2015; Del Moro et al. 2016; Georgakakis et al. 2017; Masini et al. 2018; Lanzuisi et al. 2018). In the comparison, we do not show the hard X-ray luminosity from Ricci et al. 2015 and Masini et al. 2018 since these observations were in harder X-ray bands and the $2-10\,{\rm keV}$ X-ray luminosity was not available. The absorbed fraction $F_{\rm abs}$ in the top panel of Figure 3 is defined as the fraction of absorbed quasars relative to total CTN quasars. The compton thick fraction $F_{\rm CTK}$ in the middle panel of Figure 3 is defined as the fraction of CTK quasars relative to all quasars. We find a good agreement with other observations in the absorbed quasar fraction which monotonically increases towards higher redshift. Our fiducial model (the Ueda et al. 2014 model) is in agreement with the Buchner et al. 2015 and the Aird et al. 2015a models. Besides, we also find a good consistency in the CTK quasar fraction with most of the observations, except for Aird et al. 2015a which determined the $N_{\rm H}$ distribution by reconciling the hard X-ray luminosity function of soft X-ray and hard X-ray selected quasars. Compared with the Buchner et al. 2015 model, the Ueda et al. 2014 model is consistent with it except for mild differences at $z<2$. We note that some recent studies (Masini et al. 2018; Georgantopoulos & Akylas 2019) using NuSTAR, which is more sensitive in the hard X-ray, found very small lower bounds of $F_{\rm CTK}$, $\sim 10-20\%$. Assuming that the CTK quasars are completely absent in observations, the uncertainty in $F_{\rm CTK}$ can result in $\log{\big((1-F_{\rm CTK}^{\rm min})/(1-F_{\rm CTK}^{\rm max})\big)}\sim 0.2\,{\rm dex}$ uncertainty in the binned estimations of the bolometric QLFs. In the bottom panel of Figure 3, we show the $N_{\rm H}$ distribution at $\log{L_{\rm X}}=43.5$, $z=0.05$ comparing different models (Ueda et al. 2003; Gilli et al. 2007; Treister et al. 2009; Ueda et al. 2014; Aird et al. 2015a).

Given $N_{\rm H}$, we calculate the extinction using the photoelectric absorption cross section in Morrison & McCammon 1983 and the non-relativistic Compton scattering cross section. To determine the dust abundance, a dust-to-gas ratio is required. In Hopkins et al. 2007, a constant dust-to-gas ratio was assumed to convert ${\rm N}_{\rm H}$ to dust column density and a SMC-like extinction curve from Pei 1992 was adopted. However, in this work, we find that these assumptions along with our fiducial ${\rm N}_{\rm H}$ distribution model result in a systematic inconsistency between UV, B band and X-ray observations. The UV and B band luminosities are under-predicted and the phenomenon is more severe in the UV than in the B band in a luminosity and redshift-dependent manner. This indicates that the extinction in the optical/UV is over-predicted by the model with the constant dust-to-gas ratio and the SMC-like extinction curve. Observations have revealed that the mass-metallicity relation of galaxies has a redshift evolution (Zahid et al. 2013, e.g.,) with the gas-phase metallicity of typical quasar host galaxies dropping $\sim 0.5\,{\rm dex}$ from $z=0$ to $z=2$. Similar evolution was also seen in numerical simulations (Ma et al. 2016, e.g.,). Assuming that the dust-to-metal ratio remains a constant, the decrement in the gas-phase metallicity of quasar host galaxies will lead to a decrement in the dust-to-gas ratio at higher redshift. In addition, some observations have suggested that the extinction curve of AGN might be shallower than the commonly assumed SMC-like extinction curve (Maiolino et al. 2001; Gaskell et al. 2004; Gaskell & Benker 2007; Czerny et al. 2004, e.g.,). Given the observational updates, we choose to adopt a redshift-dependent dust-to-gas ratio which scales as the gas-phase metallicity given by the fit in Ma et al. 2016. The value of the dust-to-gas ratio in the local Universe still follows Hopkins et al. 2007 with $(A_{\rm B}/N_{\rm H})=8.47\times 10^{-22}\,{\rm cm}^{2}$. We adopt the Milky Way-like extinction curve in Pei 1992 which is shallower than the SMC-like curve. Although the extinction curve of quasars does not exhibit the $2175\text{\AA}$ bump feature as found in the Milky Way, our results are not affected by this since none of the bands we study in this paper are close to $2175\text{\AA}$.

We note that the extinction in the X-ray would also be affected by the decrement of the gas-phase metallicity. In addition, although the metallicities of quasar host galaxies decrease with redshift, the metallicities of broad line regions do not evolve as strongly. The relative contributions of host galaxies and near quasar obscuration are still largely unknown. Here, our choice in the dust-to-metal ratio empirically prefers the scenario that near quasar obscuration is more important in the X-ray and obscuration in host galaxies contributes more to the extinction in the optical/UV. Our choice is motivated by making the X-ray and optical/UV observations more consistent with each other at all redshifts. Similar argument applies to our choice of the extinction curve. The shallow extinction curves found in some studies are still under debate and our choice here is only for empirical needs.

The extinction model and the bolometric corrections introduced in this and previous sections allow us to link the bolometric QLF with the observed QLF in a certain band, and resolve the discrepancies described above. We note that for all the subsequent analysis in the paper, unless otherwise specified, the QLFs presented include both the obscured and unobscured AGN and the observed QLFs presented take account of dust and gas extinction described in this section.

Bolometric Quasar Luminosity Function

We first study the bolometric QLF at a certain redshift. Following the standard practice, we parameterize the bolometric QLF with a double power-law:

where $\phi_{\ast}$ is the comoving number density normalization, $L_{\ast}$ is the break luminosity, $\gamma_{1}$ and $\gamma_{2}$ are the faint-end and bright-end slopes respectively. We note that the conventions for double power-law are sometimes different. In optical/UV studies, the double power-law is usually defined as:

where $\phi^{\prime}_{\ast}$ and $\phi^{\prime\prime}_{\ast}$ are the comoving number density normalizations with different units, $M_{\ast}$ is the break magnitude, $\alpha$ and $\beta$ are the faint-end and bright-end slopes respectively. In our notation, it gives $\alpha=-(\gamma_{1}+1),\,\beta=-(\gamma_{2}+1),\,\phi^{\prime}_{\ast}=\phi_{\ast}/\ln{10}$ and $\phi^{\prime\prime}_{\ast}=0.4\phi_{\ast}$.

For a given bolometric QLF, we can convolve it with the bolometric corrections and extinction corrections discussed in Section 3 to get the predicted observed QLF in a certain band at the redshift we study. We fit the parameters of the bolometric QLF to match the prediction with the observational binned estimations in all bands at the redshift. We select binned estimations of the QLF from our observation compilation listed in Table 5. A data set is selected if the redshift bin of that observation covers the redshift we study. Since the statistical mean redshift of the quasar samples in the binned estimations in observations does not necessarily perfectly match the redshift we study, we correct the binned estimations with a model-dependent method (referred to as "number density correction" in this paper). To be specific, for each data set in the UV, we first use the UV QLF model constrained by Kulkarni et al. 2018 (the Model 2 of the paper) to calculate the "expected" number densities at the redshift we study and at the luminosities where the data points are located. Then, we calculate the mean of the logarithm of the "expected" number densities, representing a mean level of quasar number density. Since the observed quasar samples may center on a slightly different redshift, it is likely that the observed data points exhibit a systematic shift from this "expected" mean level of number density. So we rescale the observed data points to have the "expected" mean value at the redshift we study. We also perform this correction to the X-ray data points with the X-ray QLF model constrained by Miyaji et al. 2015 and to the IR data points with the IR QLF models constrained therein. We note that this correction is model-dependent but the models we choose are representative and have the widest redshift coverage in their bands. They are in good agreement with the observations in their bands. In most of the cases, this correction step improves the clustering of data points from different investigations and reduces the potential bias in redshift estimations of observations. Combining all corrected data points, we can derive the best-fit parameters of the bolometric QLF. The best-fit parameters at some selected redshifts are listed in Table 3. The best-fits at all selected redshifts are shown in Figure 4 with gray points. In the following, we will refer to these fits as the local "free" fits (see Table 2 for details), since none of the parameters are fixed during fitting.

Since the parameters of the double power-law bolometric QLF have significant degeneracy, which manifests as large covariance in fitting, the best-fit parameters exhibit large coherent fluctuations at some redshifts. The degeneracy prevents us from finding the optimal functional form to describe the redshift evolution of the parameters. To improve the fits, we fix the number density normalization to depend linearly on redshift which is quite clear even in the "free" fits. The linear relation is determined by the best-fits at $z=0.4-3.0$. We then redo the fitting at redshifts outside $z=0.4-3.0$ with $\phi_{\ast}(z)$ fixed. Apart from that, we find that the bolometric QLF at $z\geq 5.8$ behaves as a single power-law at least in the regime covered by existing observations. Thus we reduce the fitting formula to a single power-law by restricting the faint and bright-end slope to be the same at these redshifts. The fitting procedure with these updates is referred to as the local "polished" fits (see Table 2 for details). The "polished" best-fits are also shown in Figure 4 with blue points. Based on the local "polished" fits, the bright-end slope and break luminosity evolution clearly have a double power-law shape, similar to what was seen in Hopkins et al. 2007, and the faint-end slope has a polynomial-like dependence on redshift.

2 Parameterized evolution model of the bolometric QLF

In this section, we aim to describe the evolution of the bolometric QLF with simple formulae and to perform a global fit on all the observational data at all redshifts. Following the discussion in the previous section, we describe the QLF as a double power-law with parameters that evolve with redshift as:

where $T_{\rm n}$ is the n-th order Chebyshev polynomial and $z_{\rm ref}$ is chosen to be $2$. The evolution of the bolometric QLF is therefore controlled by 11 parameters: {$a_{0}$, $a_{1}$, $a_{2}$}; {$b_{0}$, $b_{1}$, $b_{2}$}; {$c_{0}$, $c_{1}$, $c_{2}$}; {$d_{0}$, $d_{1}$}. This parameterization is adequate to describe the evolution of the bolometric QLF parameters. We have tried to extend the parameterization with higher order polynomials and find their contributions are negligible.

In the next step, we perform a global fit (referred to as the global fit A, see Table 2 for details) on all the observational data from the compilation at all redshifts simultaneously. To do this, we adopt a Monte Carlo Markov Chain (MCMC) method using the emcee https://emcee.readthedocs.io/en/stable/ package (Foreman-Mackey et al. 2013). Given a proposed parameter set of the evolution model, we calculate the resulting observed QLF in bands and compare that with observational data. For a redshift bin of a given data set, the predicted observed QLF is calculated at the center of the redshift bin. The observational data points are also rescaled to the center of the redshift bin with the number density correction discussed in the previous section. The likelihood function is then calculated in a standard way:

where $\log{\phi_{\rm mod}}$ and $\log{\phi_{\rm obs}}$ are the predicted and observed number density respectively, $\sigma_{\rm n}$ is the uncertainty of the measurement and $W(z)$ is a weighting function introduced to balance the statistical power of high and low redshift data. (Otherwise, the fact that there is more data at low redshifts would skew the fits, sacrificing large discrepancies at high redshifts for marginal improvements at low redshifts.) The summation is taken over all the observational data points at all redshifts. We choose $W(z)=1$ when $z<3$, $W(z)=\Big(\dfrac{1+z}{1+3}\Big)^{2}$ when $3\leq z<4$ and $W(z)=\Big(\dfrac{1+z}{1+4}\Big)^{3}\Big(\dfrac{1+4}{1+3}\Big)^{2}$ when $z\geq 4$. This weighting function makes the weights of data points roughly the same at $z=2-6$ and helps achieve a converged and decent fit on high redshift data. We adopt uniform priors for all the parameters involved, so that the posterior probability function is the same as the likelihood function given above. The global best-fit parameters of this evolution model are listed in Table 4 and this best-fit model will be referred to as the global fit A. In Figure 5, we show the best-fit bolometric QLFs at $6$ selected redshifts compared with the observational data converted onto the bolometric plane with the bolometric corrections and the $N_{\rm obs}/N_{\rm mod}$ method (moving data points across different QLF planes by fixing the ratio between observed and model-predicted number densities). In general, the global fit A does comparably well to the local best-fit at each redshift in matching the observational data. The best-fit bolometric QLFs of the global fit A are qualitatively different from the Hopkins et al. 2007 model. The bright end of the QLF is steeper at $z\gtrsim 2$. The faint end of the QLF is steeper at $z\gtrsim 3$ and becomes progressively steeper at higher redshifts. We achieve a better agreement with observations than the Hopkins et al. 2007 model at $z\gtrsim 3$. The evolution of the double power-law parameters of the bolometric QLF determined by the global fit A is also shown in Figure 4 with purple lines. In the top left (right) panel of Figure 4, we indicate with yellow dashed line (shaded region) the regime where integrated luminosity towards infinite low (high) luminosity will diverge. Compared with Hopkins et al. 2007, extrapolating our new model to $z>7$ will not lead to any divergence at the bright end. However, the integrated luminosity at the faint end will diverge at $z\gtrsim 6$ due to the steep faint-end slope constrained in the global fit A. In the bottom left panel of Figure 4, the colormap shows smoothed distribution of the observational data points converted on to the bolometric plane with the bolometric corrections. Darker colors indicate regions with more data points. At $z\gtrsim 5$, the void of data points approaches the break luminosity, indicating that the fits at those redshifts are affected by limited data points at the faint end. For example, in the extreme case that there are no data points fainter than the break luminosity, the fitted faint-end slope will simply be equal to the bright-end slope. The steepening of the faint-end slope at high redshift we find in both the local fits and the global fit A may be seriously affected by this. So, in parallel to the global fit A introduced above, we perform another independent global fit (referred to as the global fit B, see Table 2 for details) assuming a different evolution model for the faint-end slope. We adopt the function form used in Hopkins et al. 2007:

where $z_{\rm ref}$ is again chosen to be $2$. Different from the function form used in the global fit A, the faint-end slope here is restricted to evolve monotonically as a function of redshift. This will by construction prevent the steepening of the faint-end slope at high redshift. Other double-power-law parameters have the same evolution model as the global fit A. We perform exactly the same fitting procedure for this model and the results are also shown in Figure 4 with pink lines. The best-fit bright-end slopes, break luminosities and number density normalizations of the global fit B are similar to those of the global fit A. However, the faint-end slope in this model remains shallow at $z\gtrsim 5$ on contrary to the global fit A. The best-fit bolometric QLFs of the global fit B are also presented in Figure 5 compared with observational binned estimations. The global fit B is consistent with observations equally well as the global fit A while behaves qualitatively differently at the faint end at $z\gtrsim 5$. Future observations are required to test these two models.

The evolution models of the bolometric QLF we described above are constrained by observational data at $0<z<7$. Making predictions beyond the redshift frontier certainly requires extrapolations of the model. For the global fit A, the best-fit faint-end slope becomes the same as the bright-end slope at $z\sim 7$ and the double power-law bolometric QLF tends to behave like a single power-law approaching $z\sim 6-7$. Therefore, extrapolating to $z>7$, we postulate that the bolometric QLF simply has a single power-law shape. The evolution of the single power-law slope follows the extrapolation of the evolution of the bright-end slope at $z<7$. On the other hand, for the global fit B, the faint-end slope remains shallow at $z>7$ and we can simply extrapolate the evolution of the parameters to get luminosity functions at $z>7$. We note all these extrapolations involve assumptions on the shape of the QLF in the regime where no observational evidence is available. There are serious uncertainties there.

3 Tensions in the UV QLF at z=4−6z=4-6

The measurements of the UV QLF presented in Giallongo et al. 2015, followed by the updates in Giallongo et al. 2019, indicated a high number density of faint AGN at $z\sim 4-6$. This has motivated conjectures on whether quasars alone can be responsible for the reionization of hydrogen at $z>6$. However, other recent observations (Akiyama et al. 2018; Matsuoka et al. 2018; McGreer et al. 2018, e.g.,) have presented measurements that are in conflict with the Giallongo et al. 2015 results (as illustrated in Figure 4 in Giallongo et al. 2019). These tensions serve as a reminder that the potential uncertainties associated with the selection of quasars and host galaxy contamination at high redshift are still substantial.

In the fiducial analysis of this paper, we do not include the Giallongo et al. 2015 data in our fits. In order to check the robustness of our QLF constraints in the UV, we investigate the tensions in the UV QLF at $z\sim 4-6$ in Figure 6. We show the UV QLF determinations with various approaches at $z=4.2,4.8,5.8$ including the Giallongo et al. 2015 measurement, the compiled observational binned estimations in the UV and X-ray and the prediction from the global fits A and B. The redshifts are chosen to be close to the centers of the redshift bins in Giallongo et al. 2015. The Giallongo et al. 2015 data points (the orange crosses at the faint end) are clearly in tension with other observations in the intermediate luminosity range. The X-ray data points are moved onto the UV QLF plane with the $N_{\rm obs}/N_{\rm mod}$ method. They also disfavor the high number density of faint quasars measured by some UV observations. We show the UV QLFs constrained by Kulkarni et al. 2018 in green dashed lines. The overall normalization of the Kulkarni et al. 2018 QLFs is consistent with that of the X-ray data, despite a somewhat steeper evolved faint-end slope. Both of the global fits A and B achieve a better agreement with multi-band observational data than the Kulkarni et al. 2018 model, provided that the observational data themselves are internally consistent. At the faint end where observation data is limited, the global fits A and B behave differently. The global fit B predicts shallow faint-end slope at $z\gtrsim 5$ while the global fit A predicts progressively steeper faint-end slope. Given current available observational constraints, we are not able to tell which model is more accurate.

Evolution of the bolometric QLF

In this section, we explore the evolution of the bolometric QLF in detail and investigate the physical interpretation of the evolution based on our global fits discussed in Section 4.

In Figure 7, we compare the bolometric QLFs at different redshifts predicted by the global fits A (solid lines) and B (dashed lines). We divide the evolution of the bolometric QLF into two phases, the early phase at $z\gtrsim 2-3$ and late phase at $z\lesssim 2-3$. In the early phase, the bolometric QLF rises up monotonically following the hierarchical build-up of structures in the Universe. For the global fit A, approaching lower redshift, the relative abundance of faint quasars decreases accompanied by the increased abundance of brighter populations, forming a sharper "break" in the QLF. As a consequence of this change in the relative abundance, the faint-end slope becomes shallower and the bright-end slope becomes steeper. For the global fit B, the relative abundance of faint and bright quasars remains stable towards lower redshift, accompanied by the growth of the break luminosity. In both fits, the evolution at the bright end ($L_{\rm bol}\gtrsim 48$) is milder than that in the intermediate luminosity range. In the late phase, the bolometric QLF stops rising up. Instead, the bolometric QLF shows a systematic and continuous horizontal shift towards the low luminosity regime. The faint end has almost no evolution in this phase. This indicates processes other than the hierarchical build-up of structures dominating the evolution of the quasar population at late times. AGN feedback is potentially responsible for this evolutionary pattern. AGN feedback is believed to shut down the supply of cold gas to galaxy centers and thus could systematically decrease the bolometric quasar luminosities. Surprisingly, at $z\lesssim 0.5$, the bright end stops evolving and the bright-end slope becomes slightly shallower again. We note that the global fits A and B give similar evolutionary pattern in the late phase. Across the entire evolution history of the QLF, the evolution at the bright end of the bolometric QLF is apparently milder compared to other luminosity regimes. This suggests potential regulation on the abundance of the most luminous quasars. In Figure 7, we also present the bolometric QLF extrapolated to $z=8,10$. The extrapolations are done as introduced in Section 4.2. The rapidly dropping number density normalization makes the detection of quasars progressively difficult at these redshifts.

In Figure 8, we show the evolution of the cumulative number density of quasars in different luminosity bins in different bands. We show the predictions from the global fits A and B with purple and pink lines, respectively. The cumulative number densities predicted from the two models overlap in the bright luminosity bins. The global fit B predicts lower number density than the global fit A in the faintest UV/X-ray luminosity bin at $z\gtrsim 3$. Apparently, the number density of faint quasars peaks at lower redshift than that of bright ones, consistent with the observed "cosmic downsizing" trend (Cowie et al. 1996; Barger et al. 2005; Hasinger et al. 2005, e.g.,) of AGN at $z\lesssim 2-3$. Compared with the Hopkins et al. 2007 model which is shown in red dashed lines, our models agree well at $z\lesssim 2$ but differences show up at high redshift where new data from the past decade modifies the predictions. Since we predict steeper bright-end slopes than the Hopkins et al. 2007 model at $z\gtrsim 2$, it is not surprising that we predict lower number density at $2<z<6$ in the most luminous bin of the bolometric luminosity. The lower number density normalization we predict at high redshift gives rise to the lower cumulative number density in the UV, when integrated down to the faint end, compared with the Hopkins et al. 2007 model. In the faintest bin of the UV luminosity, at $z\gtrsim 6$, the prediction of the global fit A in the cumulative number density does not drop as fast as the Hopkins et al. 2007 model and the global fit B primarily because it predicts steeper faint-end slopes at those redshifts. In the UV, we also compare our prediction with the results in Kulkarni et al. 2018 which is an optical/UV-only study. In the bright luminosity bins, we are consistent with their estimations. However, at the faint end, we predict much lower number density of quasars at $z\gtrsim 2$ primarily driven by the much less steep faint-end slope we constrain. We note that the estimations of the number density in Kulkarni et al. 2018 did not reach $M_{\rm UV}\sim-21/-18$, so their predictions on the cumulative number density depends on the extrapolation of their measurements at brighter parts ($M_{\rm UV}\sim-23$) of the QLF. The steeper faint-end slope they constrained results in the higher cumulative number density in their prediction at $z\gtrsim 2$. The steep faint-end slope of UV QLF constrained in Kulkarni et al. 2018 is potentially affected by the paucity of X-ray observations in their study, which provide better constraints at the faint end than present UV observations. Crucially, our models do not have the unphysical upturn at $z>6$ in the Kulkarni et al. 2018 model.

In the top left panel of Figure 9, we compare the faint-end slope of our best-fit bolometric QLF with the faint-end slope of the rest-frame UV luminosity function of galaxies observed at $z=0-8$. The predictions from the global fits A and B are shown in purple and pink lines, respectively. For the galaxy UV luminosity function (GUVLF), constraints on the faint-end slope come from: observations (Duncan et al. 2014; Bowler et al. 2015; Bouwens et al. 2015; Parsa et al. 2016; Finkelstein 2016; Mehta et al. 2017; Atek et al. 2015; Atek et al. 2018; Ishigaki et al. 2018) and theoretical studies (Jaacks et al. 2012; Tacchella et al. 2013; Mason et al. 2015; Wilkins et al. 2017; Tacchella et al. 2018; Yung et al. 2018). The faint-end slopes of the bolometric QLF and the GUVLF are roughly the same at $z\lesssim 2$. Towards higher redshift, the faint end of the GUVLF starts to become steeper at $z=2-3$ where the QLF is still flat. For the global fit A, the faint-end slope of the QLF soon catches up that of the GUVLF and become even steeper at $z\gtrsim 5$. For the global fit B, the faint-end slope of the QLF remains shallow at high redshift. Again, these differences are caused by the paucity of observations at the faint end of the QLF.

In the other three panels of Figure 9, we compare the UV QLF with the GUVLF at $z=2,4,6$. Both the binned estimations and the best-fit luminosity function models are shown. The binned estimations of the GUVLF include: the compilation from Finkelstein 2016 at $z=4-10$, Alavi et al. 2014; Mehta et al. 2017 at $z=2$, Parsa et al. 2016 at $z=2,4$, van der Burg et al. 2010 at $z=4$, Bouwens et al. 2017; Atek et al. 2018 at $z=6$. We use the best-fit Schechter function in Finkelstein 2016 for the blue curves in the figure. The point where the UV QLF and the GUVLF cross each other becomes progressively higher from $z=6$ to $z=2$ which indicates enhanced significance of quasars at late times. But at all redshifts, the GUVLF appears to strongly dominate over the faint quasar UVLF at $M_{\rm 1450}\gtrsim-23$, wherever data exists. This is true even in models predicting high number density of faint quasars. It is also clear that the global fits A and B only show discrepancies in the regime where no observation is available. For the shaded regions, the two vertical boundaries show the single-visit and final detection limits of the Legacy Survey of Space and Time (LSST Science Collaboration et al. 2009, LSST, ) which will be conducted with the Simonyi Survey Telescope at the Vera Rubin Observatory. The horizontal boundary shows a reference number density corresponding to one object in the field-of-view of LSST ($\sim 20000\,{\rm deg}^{2}$) with a survey depth $\Delta z=1$. As illustrated in the figure, LSST will expand our knowledge by observing faint quasars at high redshift. This will be particularly important in resolving the knee of the QLF at $z\gtrsim 6$, if it exists, and more reliably determine the faint end slope of the QLF at high redshift. Meanwhile, LSST will boost the statistics of both galaxies and quasars at the bright end.

Implications and Predictions

In this section, we will make predictions based on our global best-fit bolometric QLF models. The predictions involve quasars’ cumulative emissivity in the UV and their contribution to hydrogen ionization, the cosmic X-ray radiation background spectrum, the evolution of the cosmic SMBH mass density and the local SMBH mass functions. Comparing the predictions with observations in these independent channels tests the validity of our bolometric QLF model. We note that, for all the predictions made in this section, the global fits A and B give indistinguishable predictions. The major difference between the global fits A and B is the faint-end slope which does not have strong impact on cumulative luminosities of all quasars, unless the faint-end slope is extremely steep. Therefore, in the following sections, we will only show the result of the global fit A and refer to it as "the global fit" for simplicity.

Faint galaxies have long been considered the dominant source of ionizing photons for the reionization of hydrogen in the Universe (Kuhlen & Faucher-Giguère 2012; Robertson et al. 2015, e.g.,). However, some observations of high-redshift quasars (Giallongo et al. 2015; Giallongo et al. 2019, e.g.,) have inferred much higher number density of quasars at the faint end than other measurements. This suggets the idea that faint quasars could potentially account for the reionization photons (Haardt & Salvaterra 2015). In this section, we quantify the quasar contribution to the photoionization of intergalactic hydrogen using the bolometric QLF derived in this paper.

Following standard modeling of UV background (Haardt & Madau 1996; Haardt & Madau 2012; Faucher-Giguère et al. 2009; Faucher-Giguère 2020; Khaire & Srianand 2019, UVB; e.g.,), the $\rm HI$ photoionization rate is:

where $\sigma_{\rm HI}(\nu)$ is the $\rm HI$ photoionization cross section and $n_{\nu}(\nu,z)$ is the number density of ionizing photons per unit frequency at redshift $z$. In principle, ionizing photons emitted at all $z^{\prime}>z$ should contribute to the ionizing background $n_{\nu}(\nu,z)$:

where $\epsilon_{\nu}(\nu_{\rm em},z^{\prime})$ is the comoving emissivity of $\rm HI$ Lyman continuum sources at redshift $z^{\prime}>z$ at emitting frequency $\nu_{\rm em}=\nu(1+z^{\prime})/(1+z)$ and $\tau_{\rm eff}(z,z^{\prime},\nu)$ is the effective optical depth of photons at $z$ emitted at $z^{\prime}$. First, to simplify the calculation, we adopt the "local source" approximation (Schirber & Bullock 2003; Faucher-Giguère et al. 2008b; Hopkins et al. 2007, e.g.,), which assumes that only ionizing sources with optical depth $\tau_{\rm eff}\leq 1$ contribute to the ionizing background (we will relax this assumption below). Then approximately, Equation 18 is reduced to:

where $\Delta l(\nu,z)$ is the mean free path of ionizing photons defined by $\tau_{\rm eff}(\Delta l)=1$. Based on the results in Faucher-Giguère et al. 2008b, the frequency dependence of the mean free path can be described as $\Delta l(\nu,z)=\Delta l(\nu_{912},z)\,(\nu/\nu_{912})^{3(\beta-1)}$, where the $\beta$ is the power-law index of the intergalactic HI column density distribution. For our local source approximation, we assume that the HI column distribution can be approximated by a single power-law index $\beta\approx 1.5$ (Madau et al. 1999, e.g.,). Assuming a power-law shape for the extreme UV quasar continuum, we have:

Since $\sigma_{\rm HI}(\nu)\propto\nu^{-3}$, the $\sigma_{\rm HI}(\nu)\,c\,n_{\nu}(\nu,z)$ term in Equation 17 will be proportional to $\nu^{-(4+\alpha_{\rm UV}-3(\beta-1))}$. Then integrating Equation 17 gives $\Gamma_{\rm HI}(z)=\dfrac{\sigma_{\rm HI}(\nu_{912})\,c\,n_{\nu}(\nu_{912},z)\,\nu_{912}}{3+\alpha_{\rm UV}-3(\beta-1)}$. Plugging in Equation 19, we finally obtain:

where we adopt $\alpha_{\rm UV}$ = 1.7 (Lusso et al. 2015) and $\Delta l^{912}_{z=3.5}=50\,{\rm Mpc}$ with a power-law index $\eta=4.44$ for the redshift dependence of $\Delta l$ (Songaila & Cowie 2010). Here, we only consider the contribution from quasars. The emissivity at Lyman limit $\epsilon_{912}(z)$ can be linked with the UV emissivity $\epsilon_{1450}(z)$ of quasars as:

assuming a power-law shape of the UV continuum with index $-0.61$ (Lusso et al. 2015), which is in good agreement with our SED model. We note that here we have assumed the escape fraction $f_{\rm esc}=100\%$ for the ionizing photon produced by quasars. It is common to adopt $100\%$ escape fractions for optically-bright quasars. However, some fraction of quasars have known dust and gas obscuration that would severely limit the escape of ionizing photons. So, the results we derive here should be interpreted as an upper limit of quasars’ contribution to ionization. To derive the comoving UV emissivity of quasars, we integrate luminosity over the UV QLF predicted by our global best-fit models:

where $L^{0}_{\nu}$ is the zero-point luminosity of the AB magnitude system, $M_{\rm min}$ and $M_{\rm max}$ are the magnitude bounds for integration. We adopt $M_{\rm min}=-18$ and $M_{\rm max}=-35$.

In the left panel of Figure 10, we present the predicted Lyman limit comoving emissivity $\epsilon_{\rm 912}$ versus redshift. At low redshifts, our prediction is close to the results of Hopkins et al. 2007 and Kulkarni et al. 2018. At high redshifts, our prediction agrees well with the Haardt & Madau 2012 model and is much lower than the Kulkarni et al. 2018 prediction due to the less steep faint-end slope we constrain. The prediction is in agreement with observational estimations in narrow redshift bins from Masters et al. 2012; Palanque-Delabrouille et al. 2016; Akiyama et al. 2018. We predict lower emissivity compared to the estimations of McGreer et al. 2018; Parsa et al. 2018. We fit the redshift dependence of the emissivity with a five-parameter functional form (Haardt & Madau 2012):

In the right panel of Figure 10, we present our prediction for the hydrogen photoionization rate contributed by quasars. We find that the prediction using the local source approximation severely over-predicts the hydrogen photoionization rate at $z\lesssim 2$ where the mean free path of ionizing photons grows comparable to (and eventually larger than) the Hubble radius, so that the local source approximation fails significantly. Therefore, we perform a full UVB calculation using the method described in Faucher-Giguère 2020. The result is also shown in the right panel of Figure 10. The prediction from this full UVB calculation almost overlaps with the prediction with the local source approximation at $z\gtrsim 4$, despite slight differences. The slight differences are due to more physics incorporated in the full UVB calculation that make the UVB spectrum (filtered by IGM absorption and including recombination emission) different from the simple power-law that we have assumed above (Faucher-Giguère 2020, see). We compare our predicted hydrogen photoionization rates (from quasars only) with observational inferences of the total rates from Wyithe & Bolton 2011; Calverley et al. 2011; Becker & Bolton 2013; Gaikwad et al. 2017; D’Aloisio et al. 2018. The predicted hydrogen photoionization rate contributed by quasars is an order of magnitude lower than the measured total rate at $z\sim 6$ and only becomes close to the total rate at $z\lesssim 3$. The results indicate that quasars are subdominant to the hydrogen reionization at $z\gtrsim 6$, but they start to dominate the ionization budget at $z\lesssim 3$. Interestingly, that the hydrogen photoionization rates predicted using our new bolometric QLF are quite similar to the results of Hopkins et al. 2007, which used a different bolometric QLF and adopted a different mean free path model. We have assumed that all ionizing photons produced by the quasar can escape the host galaxy even for faintest quasars. Given this favorable assumption, the predicted contribution of quasars to the hydrogen reionization is still subdominant. Similar conclusion has been reached in Ricci et al. 2017 who adopted a different approach. We have tested that, even including the Giallongo et al. 2015 data in the fit and neglecting all the data points that are incompatible with it, quasars can only have a maximum of $\sim 50\%$ contribution to the ionization budget at $z\sim 5.8$, under the assumption that the escape fraction $f_{\rm esc}=100\%$ even for quasars much fainter than typical star-forming or Seyfert galaxies.

2 Cosmic X-ray background

Since quasars dominate the radiation budget in the X-ray in the Universe, the cosmic X-ray radiation background (CXB) serves as an important channel to cross check our model of the bolometric QLF. The observation of the CXB does not require spatially resolving and identifying quasars and thus can even probe the contribution from faint-end quasars at any redshift.

In general, to get the cosmic radiation background contributed by quasars, we integrate the spectrum of quasars at $z=0-7$ as:

where $\nu_{\rm em}=(1+z)\nu$ and $\dfrac{{\rm d}V}{{\rm d}\Omega{\rm d}z}(z)$ is the differential comoving volume element at $z$. $\epsilon_{\nu}(\nu_{\rm em},z)$ is derived by integrating over the luminosity function of the emission at $\nu_{\rm em}$ predicted by our best-fit model.

In practice, we have found that simply adopting the X-ray SED template with the median photon index $\Gamma=1.9$ leads to an under-prediction for the CXB. Considering that the photon index has a significant scatter, $\sim 0.2$, the stacked SED of quasars should have a very different shape from a simple cut-off power-law. Therefore, in making predictions on the CXB, we adopt the stacked SED of $1000$ sampled SEDs with a normal distribution of photon indexes with median value $1.9$ and scatter $0.2$. In Figure 11, we show the predicted CXB spectrum and compare it with the measurements from Gendreau et al. 1995; Gruber et al. 1999; Churazov et al. 2007; Ajello et al. 2008; Moretti et al. 2009; Cappelluti et al. 2017. For simplicity, we have assumed the galaxies’ contribution to the CXB to be a constant $2\,{\rm keV}^{2}{\rm s}^{-1}{\rm sr}^{-1}\,{\rm keV}^{-1}$. We find our predicted CXB spectrum agrees well with observations at high energy end while it is roughly $\sim 0.05\,{\rm dex}$ lower at $E\lesssim 20\,{\rm keV}$. Imperfectness in the extinction model may be responsible for this though it is hard to argue the source of this level of inconsistency. The Hopkins et al. 2007 model systematically over-predicts the CXB spectrum. We also show separately the contribution to CXB from CTK, absorbed CTN and unabsorbed CTN AGN. The absorbed CTN AGN are the major sources of the CXB in the high energy regime while the unabsorbed CTN AGN overtake at $E\lesssim 3\,{\rm keV}$. The CTK AGN are subdominant to the CXB.

3 Growth history of SMBHs

The bolometric quasar luminosity is connected with the accretion of the SMBH that powers quasar activities. Thus, based on our bolometric QLF model, constraints can be put on the growth history of SMBHs in the Universe. Here we focus on the evolution of the cosmic SMBH mass density and the SMBH mass function.

Assuming a constant averaged radiative efficiency $\epsilon_{\rm r}\simeq 0.1$ for the SMBH accretion, the bolometric quasar luminosity can be related to the accretion rate of the SMBH as:

Therefore, the integrated luminosity density can be translated to the rate of change in the total SMBH mass density as:

where we adopt $\log{L_{\rm min}}=43,\log{L_{\rm max}}=48$ here. Starting from an initial redshift for SMBH growth $z_{\rm i}$ and integrating over redshift, we derive the evolution of $\rho_{\rm BH}$. In Figure 12, we show the redshift evolution of the SMBH mass density with $z_{\rm i}=10,7,4$ in the red, blue and green lines. Shaded regions show the uncertainties when increasing or decreasing $\epsilon_{\rm r}$ by $2$ times. The build-up of the SMBH mass density is completely dominated by the accretion at $z<4$. Compared with local constraints (Shankar et al. 2004; Marconi et al. 2004; Graham & Driver 2007; Yu & Lu 2008; Shankar et al. 2009), we predict slightly higher SMBH mass density at $z=0$. There are several uncertainties that could impact the comparison made here. These local constraints were calculated by translating galaxy central spheroid properties to the mass of SMBH. New calibrations of the scaling relations between the mass of SMBH and galaxy spheroid properties (Kormendy & Ho 2013; McConnell & Ma 2013, e.g.,) have generally found higher intercepts and steeper slopes than the old calibrations. Besides, as discussed in Section 3.1, variations on the definition of the bolometric luminosity could also lead to systematic shift in the estimated radiation energy budget of SMBHs. Both of these two factors could drive the local constraints and our predictions to be more consistent with each other. However, on the other hand, the selection biases in observed scaling relations could result in an over-estimation of the local SMBH mass density (Shankar et al. 2016, e.g.,). In that case, the discrepancy of our result with local estimations indicates a higher averaged radiative efficiency than the assumed value $0.1$.

3.2 SMBH mass function

The mass function is one of the most important statistical properties of the SMBH population. In the local Universe, SMBH mass can be determined by various properties of galaxy spheroids, e.g. the velocity dispersion, the bulge mass. Both quiescent and active SMBHs’ masses can be estimated in this way. At high redshift, SMBH masses are measured based on direct radiation from the vicinity of active SMBHs. Alternatively, the SMBH mass can be related to the bolometric quasar luminosity with the Eddington ratio. Assuming an Eddington ratio distribution, one can convert the bolometric QLF to the SMBH mass function. Technically, there are two ways to achieve this:

convolve the bolometric QLF with the measured relation between Eddington ratio and bolometric quasar luminosity. This method is referred to as "convolution".

assuming an Eddington ratio distribution, fit the parameterized SMBH mass function based on the bolometric QLF. This method is referred to as "deconvolution".

For the first approach, we adopt the scaling relation (Nobuta et al. 2012):

where $\lambda_{\rm Edd}$ is the Eddington ratio. The relation was measured based on X-ray selected AGN at $z\sim 1.4$ and was demonstrated (Nobuta et al. 2012) to be consistent with what had been found in the SDSS DR5 broad-line AGN (Shen et al. 2009). We also consider the $\sim 0.4\,{\rm dex}$ scatter of this relation (Nobuta et al. 2012). Convolving the bolometric QLF with this relation, we can derive the SMBH mass function for active SMBHs. We further multiply the fraction of unabsorbed CTN AGN $F\sim 0.38$, estimated at the knee of the local X-ray QLF with our fiducial extinction model, to get the SMBH mass function of Type-1 AGN. We present the predicted SMBH mass function of Type-1 AGN in Figure 13 with the blue dashed line which is in good agreement with the observation (Kelly & Shen 2013). In order to further deduce the total SMBH (including quiescent ones) mass function, we need to correct for the fraction of AGN that are in the active phase, $f_{\rm duty}$. We find that in order to match the observational constrained total SMBH mass function in the local Universe (Vika et al. 2009; Shankar et al. 2009; Marconi et al. 2004), $f_{\rm duty}$ should take the value $\sim 0.03$. After multiplying $1/f_{\rm duty}$ to the predicted SMBH mass function of active SMBHs, we derive the total SMBH mass function shown with the blue solid line in Figure 13.

For the second approach, we assume a two component Eddington ratio distribution function (ERDF) for AGN (Tucci & Volonteri 2017):

The first component takes a Schechter function format and describes the ERDF of Type-2 AGN. The prefactor $A$ is set to normalize the total probability of this component to be $1-F$. We choose $\lambda_{1}=1.5$ and $\alpha=-0.6$ which were found in agreement with observations on low redshift Type-2 AGN (Hopkins & Hernquist 2009; Kauffmann & Heckman 2009; Aird et al. 2012). The second component takes a log-normal format and describes the ERDF of Type-1 AGN of which the parameters were determined by fitting the shape of the ERDFs from Kelly & Shen 2013 in different redshift bins and interpolating the results with a linear function (Tucci & Volonteri 2017):

We note that a consensus on the shape of the ERDF has not been reached. However, the potential influence of the ERDF assumptions should be limited (see the Appendix of Weigel et al. 2017) for our purpose here. We parameterize the total SMBH mass function as a double power-law function. For a proposed total SMBH mass function, multiplying $f_{\rm duty}=0.03$ where we found through the other method, we can derive the SMBH mass function of the active SMBHs with parameters left for fitting. We can convolve this active SMBH mass function with the assumed ERDF to derive the resulting bolometric QLF. By comparing the result with our bolometric QLF model, we derive the best-fit parameter choice for the SMBH mass function. In Figure 13, we present constraints on the local SMBH mass function from two different methods and compare it with observations of the total SMBH mass function (Marconi et al. 2004; Shankar et al. 2009; Vika et al. 2009) and observations of the Type-1 AGN mass function (Kelly & Shen 2013). The constraints from this work are in decent agreement with all the observations in the range $10^{7}$ to $10^{9.5}{\,\rm M_{\odot}}$. The "convolution" method does better at the massive end while the "deconvolution" method does better at the low mass end.

We limit our prediction to the local SMBH mass function, since the uncertainties in the ERDF, the active fraction and the absorbed fraction grow much larger at high redshift. A more comprehensive model of the SMBH population and constraints on the evolution of the SMBH mass function will be explored in future works.

Summary and Conclusions

In this paper, we update the constraints on the bolometric QLF at $z=0-7$ and make various predictions based on this model. Our technique follows the method of Hopkins et al. 2007 but with an updated quasar mean SED model and bolometric and extinction corrections. We have also extended the observational compilation in Hopkins et al. 2007 with new binned estimations of the QLF from the recent decade. These new observations allow more robust determination of the bolometric QLF at $z\gtrsim 3$. Our findings on the bolometric QLF can be summarized as:

We obtain two global best-fit models A and B with different assumptions on the evolution of the faint-end slope at high redshift. As shown in Figure 4 and Figure 5, comparing with the Hopkins et al. 2007 model, we find the bright-end slope steeper at $z\gtrsim 2$ in both the global fits A and B. In the global fit A, the faint-end slope is steeper than the Hopkins et al. 2007 model at $z\gtrsim 3$ and becomes progressively steeper at higher redshift. In the global fit B, where we adopt a monotonically evolved faint-end slope, the faint-end slope remains shallow at high redshift and is close to the prediction of the Hopkins et al. 2007 model. The uncertainties on the faint-end slope arise from the paucity of measurements of the faint-end QLF at high redshift. Apart from that, we have fixed some extrapolation problems of the Hopkins et al. 2007 model. The integrated luminosity of bright-end quasars would not blow up at $z\gtrsim 7$ and the number density normalization exhibits a more natural evolution towards higher redshift.

We investigate the current tension in the UV QLF at $z\simeq 4-6$ shown in Figure 6. We find that the high number density of faint quasars found in Giallongo et al. 2015 is disfavored when compared with current available X-ray observations. Our QLF models achieve a better agreement with the X-ray data at the faint end than the previous QLF models based on optical/UV observations only.

The evolution of the bolometric luminosity function can be interpreted as two phases separted at $z\simeq 2-3$, illustrated in Figure 7. In the early phase, the bolometric QLF rises up monotonically following the hierarchical build-up of structures in the Universe. In the late phase, the bolometric QLF shows a systematic and continuous horizontal shift towards the low luminosity regime. AGN feedback is potentially responsible for this evolutionary pattern. Surprisingly, in both phases, the evolution at the bright end ($L_{\rm bol}\gtrsim 48$) of the bolometric QLF is apparently milder compared to other luminosity regimes. This suggests potential regulation on the abundance of the most luminous quasars.

We have made predictions with this new model on the hydrogen photoionization rate contributed by quasars, the CXB spectrum, the evolution of the cosmic SMBH mass density and the local SMBH mass function. We find a general consistency with observations in these channels and our findings can be summarized as:

We find that quasars are subdominant to the hydrogen photoionization rate during the epoch of reionization at $z\gtrsim 6$. They start to dominate the UV background at $z\lesssim 3$.

The predicted CXB spectrum shown in Figure 11 agrees well with observations in the high energy regime while lies slightly lower than observations at $E\lesssim 20\,{\rm keV}$.

We predict the evolution of the SMBH mass density at $z=0-7$ shown in Figure 12. We find that the prediction is consistent with local observations and the evolution is dominated by the growth of SMBHs at $z<4$.

We make predictions on the local total SMBH mass function and the Type-1 AGN mass function shown in Figure 13. We explore two different methods, a "convolution" method and a "deconvolution" method. Both of them can generate consistent results with observations.

The new bolometric QLF model constrained in this paper can simultaneously match the multi-band observations on QLF over a wide redshift range up to $z\sim 7$. The model reveals an evolutionary pattern of the bolometric QLF at high redshift that is qualitatively different from the Hopkins et al. 2007 model. The predictions from the new model is in consistent with observations in various channels. We demonstrate the new bolometric QLF model as a solid basis for future studies of high redshift quasar populations and their cosmological impacts.

Acknowledgements

Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research grant #1715847 and CAREER grant #1455342. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; and by a Cottrell Scholar Award and Scialog Award #26968 from the Research Corporation for Science Advancement. NPR acknowledges support from the STFC and the Ernest Rutherford Fellowship scheme. GTR was supported in part by NASA-ADAP grant NNX17AF04G. DMA thanks the Science and Technology Facilities Council (STFC) for support through grant code ST/P000541/1. Numerical calculations were run on the Caltech computer cluster ’Wheeler’.

References

Appendix A Compiled observations

In Table 5, we list the observational papers compiled in this work along with the details of their observations, including the survey/fields, the band, the luminosity/redshift coverage and the number of quasar samples adopted.

Appendix B Posterior distribution in the global fit

In Figure 14, we show the posterior distribution of the four double-power-law parameters at $z=5$ in our global fit A (see Table 2 for details). The global fit A is originally done in a $11$ dimension parameter space of the QLF evolution model. Here, we project the posterior distribution onto the $4$ dimension parameter space of the double power-law function at $z=5$.

Appendix C Code and data

The code and data used in this work are publicly available at https://bitbucket.org/ShenXuejian/quasarlf/src/master/. It includes the compiled observational datasets of the QLF, the mean SED model, the pipeline for bolometric and extinction corrections, the global best-fit bolometric QLF models and all other code for the analysis done in this paper.