A Physical Model for the Quasar Luminosity Function evolution between Cosmic Dawn and High Noon

The model best fit parameters and their associated uncertainties are listed in Table 1. In Fig. 1, we show $L_{c}(M_{h})$ and the resulting $z\geq 4$ QLF for $\varepsilon_{DC}=1,0.1$ and $0.01$. Here, $L_{c}(M_{h})$ at other redshifts were derived using abundance matching between the HMF and the ‘intrinsic’ QLF, i.e. finding $L(M)$ at $z_{f}$ such that

is satisfied. The $\Sigma=0$ case with an appropriately selected $k$ value is shown for comparison purposes in the $\varepsilon_{DC}=0.01$ panels of Fig. 1. We highlight here, the impact of $\Sigma$ in reproducing the bright end of the QLF after redshift evolution. By construction, a single QLF at one redshift can be adequately fit with any $\Sigma$, although with some small deviations around $M_{UV}\sim M^{*}$ for the highest $\Sigma$ values (Ren et al., 2020). We note that low values of $\Sigma$ underestimate the number density of the brightest quasars after redshift evolution and yields a bright end that is substantially steeper than what is observed (as seen by the dot-dashed grey line in Fig. 1). On the contrary, a sufficiently high $\Sigma$ will result in a flatter bright end that is consistent with the observed evolution of the QLF. Furthermore, we find that our modeled QLFs appear to become shallower as a function of increasing $z$. This behavior of the bright end with high $\Sigma$ is better aligned with the redshift evolution of the bolometric QLF (Shen et al., 2020), whom finds shallower slopes for higher redshifts ($z>2.4$), compared to the evolution of the UV QLF with bright end slopes that become steeper at higher redshifts (Kulkarni et al., 2019). However, we emphasize that the functional form of our model QLFs are not double power laws (DPL), but more akin to ‘broadened’ Schechter functions. Any physical interpretation to the steepness of the QLF slopes can be obfuscated after forcing a double power law parameterisation to fit the QLFs, where the bright end slope, $\beta$ is sensitive to the placement of the ‘knee’, $M^{*}$ (see Section 3.5). In this model, the evolution for the shallower bright end is primarily due to the decreased abundance of $M_{h}>M_{h,0}$ halos at higher redshifts, which leads to substantial variation in quasar luminosities within the population of the most massive halos (Ren et al., 2020).

In Table 1 we find that the best-fitting $\Sigma$ is insensitive to the choice of $\varepsilon_{DC}$, whereas $k$ is relatively sensitive to $\varepsilon_{DC}$. With the best fit parameters, the evolved QLFs are consistent with existing QLF observations at $z\geq 4$ for all $\varepsilon_{DC}$, including the $z\sim 6$ Matsuoka et al. (2018) QLF that is outside of our redshift range for calibration, thus providing a validation set for our model. Since the only redshift dependent terms in Equation 3 are the HMF and the halo assembly times, our model therefore suggests that the evolution of the overall quasar density is driven by the HMF, with some possible contribution from $R_{A}(M_{h})$ and its associated power index $k$, depending on $\varepsilon_{DC}$. One interpretation for $R_{A}(M_{h})$ is to associate this quantity as a proxy for the contribution of merger-driven growth to the SMBH, which can be expected to account for some part in the SMBH growth history (Marshall et al., 2020). We see in Equation 5 that the effect of $R_{A}(M_{h})$ is to provide a boost in $L_{c}(M_{h})$. From Fig. 1, we find that this boost can be significant, up to $2$ mag, for high values of $\varepsilon_{DC}$ (e.g. $\varepsilon_{DC}=1$), and marginal, with corresponding $z\geq 4$ $L_{c}(M_{h})$ curves within $1\sigma$ of each other, for smaller values of $\varepsilon_{DC}$ (e.g. $\varepsilon_{DC}=0.01$).

While our model contains a number of parameters, the calibration step effectively fixes all of them except for $\varepsilon_{DC}$. Furthermore, from the left panels of Fig. 1, we see that the model results are only weakly depending on the choice of $\varepsilon_{DC}$. It is noteworthy that a constant $\varepsilon_{DC}$ used in our model is sufficient to fully reproduce all observed QLFs at $z\geq 4$. For all $\varepsilon_{DC}$ considered, the modeled QLFs at $z\leq 7.5$ are within $2\sigma$ from each other, with the most difference at high $z$. This potentially offers the novel possibility of using a future determination of the QLF at $z\gtrsim 7.5$ as an independent constraint for the duty cycle, distinct from derivations of the duty cycle based on clustering measurements. We explore further the effect of $\Sigma$ and $\varepsilon_{DC}$ for measurements of the duty cycle in Section 3.3.

Fig. 1 also allows us to investigate how the characteristic mass scale $M_{h,0}$ in $L_{c}(M_{h})$ evolves as a function of $z$. We find that the evolution of $M_{h,0}$ depends on $\varepsilon_{DC}$, with $\varepsilon_{DC}=0.01$ showing no visible redshift evolution and $M_{h,0}\sim 10^{12.1}M_{\odot}$, while the $\varepsilon_{DC}=1$ case displays a range of $M_{h,0}\sim 10^{12.6}-10^{13}M_{\odot}$ for $4\leq z\leq 10$. A more consistent feature is that of a characteristic luminosity $L_{c,0}$ that is constant for all $\varepsilon_{DC}$ and $z$. If we make the assumption that the brightest quasars are radiating at Eddington luminosity, then our model implies a characteristic SMBH mass. One can then interpret $L_{c,0}$ as an indication that the main mechanism for the self-regulation of SMBH growth in the average quasar is due to the quasar itself, as opposed to the radio mode accretion that is typically expected to quench star formation in massive halos with $M_{h}>10^{12}M_{\odot}$ (Croton et al., 2006). We stress that these feedback modes are not mutually exclusive with each other and that an individual quasar can be affected by either mechanism. Furthermore, we note that even for our maximal $\varepsilon_{DC}=1$ case, the shift in $M_{h,0}$ is relatively small, just $\Delta M_{h,0}=0.4$, thus it does not exclude the possibility that the main mode of self-regulation for the average quasar is still due to radio mode feedback. For smaller values of $\varepsilon_{DC}$, the model begins to completely lose the ability to discriminate between the primary mechanism that self-regulates SMBH growth, whether it is a characteristic mass $M_{h,0}$, luminosity $L_{c,0}$, or a combination thereof.

We investigate the applicability of this model for $z\leq 4$ as well (i.e. below the redshift range at which the parameters were calibrated). In Fig. 2, we show our fiducial model predictions for $z=3.1$ against the observed QLF resulting in an excellent agreement to observations for all $\varepsilon_{DC}$. However, as redshift decreases, a small evolution in $\Sigma$ is necessary to prevent the over-representation of bright quasars. For example, the $z=2.2$ QLF requires a tighter scatter with $\Delta\Sigma\sim-0.09$ relative to the nominal constant $\Sigma$ used for higher $z$ for all values of $\varepsilon_{DC}$. As redshift decreases, $\Delta\Sigma$ grows, e.g. we find $\Delta\Sigma\sim-0.19$ at $z=1.0$. Physically, it is not implausible to expect a tighter scatter in $L_{c}(M_{h})$ at lower redshift. For example, mass averaging over numerous successive merging events could restrict the spread in SMBH masses (Peng, 2007), therefore resulting in a reduced dispersion in $\Sigma$. We overlay the fiducial $\Sigma$ predictions in Fig. 2 as dotted colored lines to highlight the extent of the overproduction of bright quasars if $\Sigma$ was not corrected for in the low-$z$ regime.

The low redshift evolution of the QLF in our model is also impacted by the choice of the power law index parameter $k$. Since $R_{A}(M_{h})<1$ for $z_{f}<z_{c}$, a high $k$ damps the effective growth of the QLF compared to what would be expected simply from the evolution of the HMF. For substantial values of $k$, this may even lead to a decrease of the overall quasar number density going to lower redshifts, as seen in the $z=2$ to $z=1$ QLFs for $\varepsilon_{DC}=1$ in Fig. 2. However, we stress that $\varepsilon_{DC}=1$ is an extreme scenario where the large $k$ can lead to possibly unphysical consequences, such as a significant evolution in the SMBH - host galaxy mass relation which is not expected theoretically nor inferred from observations. In fact, the SMBH - host galaxy mass relation is broadly expected to follow the evolution in $L_{c}(M_{h})$ since neither the stellar efficiency parameter (Behroozi et al., 2012) or the Eddington ratio distribution (Kollmeier et al., 2006) show strong evidence of redshift evolution up to $z\sim 4$. Both observations (Decarli et al., 2010) and simulations (Volonteri et al., 2016; Marshall et al., 2020) support only a mild evolution in the SMBH - host galaxy mass relation. However, our modeling finds that any evolution can be largely suppressed by assuming a lower intrinsic duty cycle $\varepsilon_{DC}\leq 0.1$, which yields a lower $k$. On the other hand, assuming $k<<1$ also leads to difficulties in reproducing the observed QLF at $z=1$ under fiducial model assumptions, due to an over-representation of fainter quasars. In this scenario, we find that our model is able to approximately match the $z=1$ QLF only if we allow the duty cycle $\varepsilon_{DC}$ at $z=1$ to scale by a factor of $0.2(0.5)$ for starting values $\varepsilon_{DC}=0.01(0.1)$. A decrease in the value of $\varepsilon_{DC}$ at $z=1$ is plausible through observations of shorter quasar lifetimes, $t_{Q}\sim 10^{7}$ Myr at $z\lesssim 1$ (Porciani et al., 2004; Padmanabhan et al., 2009), corresponding to a duty cycle of $\sim t_{Q}/t_{H}\sim 10^{-3}$.

It is interesting that the assumption of a redshift-independent duty cycle yields a very good match to observations throughout a large redshift range, well beyond the calibration interval. In fact, as in the $z\geq 4$ case, we find that the evolution in the overall number density for $z\leq 4$ can be determined from the combined evolution of the HMF and $R_{A}(M_{h})$. Thus our model suggests that the ensemble phenomenology of dark matter halos in combination with AGN feedback alone is sufficient to fully describe the evolution of the QLF. There is no requirement for any change in the global ISM conditions (such as a global depletion of gas or a redshift-dependent cosmic dust extinction), at least for $z\geq 2$.

One consistent feature that our model is unable to fully replicate, even after allowing an evolution for $\Sigma$ and $\varepsilon_{DC}$, is the faint end of the $z\leq 2$ QLF. In fact, our model routinely underestimates the distribution of fainter quasars at $z=1$, reaching up to $\sim 1$ order of magnitude difference in counts for faint $M_{UV}=-22$ objects. We attribute this discrepancy to a limitation of our modeling method, specifically to our choice to calibrate the QLF at $z=4$ and $z=5$. In reality, at $z\leq 2$, we may expect a population of AGN with massive SMBHs ($M_{\rm BH}>10^{7}M_{\odot}$) that are quiescently accreting in radio-mode which ultimately contributes to the faint end of the QLF (Hirschmann et al., 2014). The calibration at high $z$ is not sensitive to this secondary population as the number of these objects is expected to be very low at high $z$. Interestingly, the inability to account for this feature is beneficial for our modeling of the QLF across most $z$, but particularly at high $z$, as this allows us to take advantage of using a simple flattening of $L_{c}(M_{h})$ at the high mass end to capture the general qualities in the evolution of the QLF. However, the natural drawback is that our model does not account for the population of AGN inside massive halos that may play a larger role at lower redshifts.

The effect of quasar luminosity scatter $\Sigma$ can impact measurements of the duty cycle derived from quasar-quasar or quasar-galaxy clustering, as $\Sigma$ decreases median halo mass hosting highly luminous quasars (Padmanabhan et al., 2009; Shankar et al., 2010; Ren et al., 2020). Our modeling framework presents the opportunity to quantitatively investigate how our calibrated $\Sigma$ can impact the duty cycle that would be inferred in a mock survey as a function of $\varepsilon_{DC}$ and a survey’s limiting magnitude. We employ the standard definition of the duty cycle (e.g. as in Eftekharzadeh et al. 2015; He et al. 2017) defined as,

where $n_{\rm QSO}$ is the number of visible quasars in a given survey volume, $V$. $\bar{M}_{\rm min}$ is the median halo mass from the set of observable (i.e. sufficiently bright) quasars and $\frac{dn}{dM_{h}}$ is the HMF. The survey volume $V$ is only important in determining the variation of $\varepsilon_{\rm{measured}}$. Thus, instead of a fixed survey area, we use a fixed volume $100$Gpc³ to simply represent an arbitrary wide-area survey. For a mock ‘survey’, we sample the HMF, calculate quasar luminosities through the stochastic relationship $L_{c}(M_{h})$, and solve Equation 6 assuming some limiting magnitude, $M_{\rm lim}$. We investigate two distinct survey arrangements over $2\leq z\leq 10$: (1) First, we use a limiting magnitude $M_{\rm lim}=-22.9$, roughly consistent with the lower luminosity limit from both of the quasar samples used in the $z<3$ SDSS-III BOSS catalog from Eftekharzadeh et al. (2015) and the low-luminosity sample from the $z\sim 4$ HSC wide-layer catalog of He et al. (2017); For (2) we assume a brighter magnitude limit of $M_{\rm lim}=-27$ to be more representative of the quasar samples used in the SDSS DR5 catalogue from Shen et al. (2007) and also the high-luminosity sample of He et al. (2017).

We conduct $10$ mock surveys for each of the scenarios presented above. In Fig. 3, we show the measured duty cycle, $\varepsilon_{\rm{measured}}$ as a function of redshift and $\varepsilon_{DC}$. In all cases, we find that $\varepsilon_{\rm{measured}}$ is relatively independent of $z$. Furthermore, we note that there is an overall good agreement between $\varepsilon_{DC}$ and $\varepsilon_{\rm{measured}}$ on the provision that the limiting magnitude is sufficiently sensitive (blue lines of Fig. 3). This surprising result is contrary to the earlier expectation that the quasar luminosity scatter, $\Sigma$ should reduce the measured clustering and therefore should result in a lower $\varepsilon_{\rm{measured}}$. While we do expect the general qualitative impact of changing $\Sigma$ on clustering to hold, we predict that given the range of $\Sigma$ and the above specific survey parameters, the significantly larger population of fainter quasars with a smaller spread of possible halo masses is sufficient to overcome the dilution of clustering strength from $\Sigma$ to recover the internal duty cycle as the measured duty cycle. In contrast, we see that a shallower magnitude limit $M_{\rm lim}=-27$ results in a reduced $\varepsilon_{\rm{measured}}$ from the impact of $\Sigma$ (orange lines of Fig. 3). In this scenario, $M_{\rm lim}$ is much brighter than the magnitude corresponding to $L_{c}(M_{h,0})$, which is determined by $\Sigma$, thus any detected quasar is biased towards having a halo mass close to $M_{h,0}$. This effectively suppresses the measured duty cycle determined by Equation 6. Thus, in order to obtain a robust measurement of the duty cycle, our modeling implies that the limiting depth of a survey needs to be at least comparable to $L_{c}(M_{h,0})$.

Under the survey conditions of Eftekharzadeh et al. (2015) and He et al. (2017), we find that that our modeling is consistent with their results assuming an internal duty cycle of a few percent, i.e. $\varepsilon_{DC}\sim 0.01$. Furthermore, this modeling prescription offers an explanation for the drastically lower duty cycle ($\varepsilon_{\rm{measured}}\sim 8\times 10^{-4}$) determined from the high luminosity $z\sim 4$ quasar sample of He et al. (2017). As the high luminosity quasar sample used in He et al. (2017) spans a magnitude range between $-24\lesssim M_{UV}\lesssim-28$, we speculate that observed quasars are overwhelmingly biased in residing within lower mass halos. In contrast, Shen et al. (2007) determine a measured duty cycle of $\varepsilon_{\rm{measured}}\sim 0.03-0.6$ from a $z>3.5$ quasar sample which includes objects $M_{UV}<-25.3$ using the band conversion in Ross et al. (2013). The upper end of this value is largely discrepant to the other measurements by Eftekharzadeh et al. (2015) and He et al. (2017) with the difference in measurements being attributed to the treatment of the large scale ($>30$ Mpc) data points. If we were to take the measurements of Shen et al. (2007) at face value for constraining the value $\varepsilon_{DC}$, then this would imply an internal duty cycle that is larger than the measured duty cycle, as the survey’s limiting depth is brighter than $L_{c}(M_{h,0})$, thus exacerbating the existing discrepancy in duty cycle values. Our modeling finds an equivalent lower limit of $\varepsilon_{DC}\sim 0.1$ using the survey parameters of Shen et al. (2007).

In Fig. 4 we use our model to provide a forecast accessible by the proposed wide-field Nancy Grace Roman Telescope (RST) mission ($\sim 2000$ deg²) as well as the total number of quasars expected over the all-sky area ($\sim 40000$ deg²). We find that RST is easily capable of finding $z\sim 8$ quasars assuming an apparent magnitude limit of $m_{UV}=26.5$, with an expected total of $N\gtrsim O(10^{1})$ detected objects for all values of $\varepsilon_{DC}$. Additionally, at the same magnitude depth, we predict that the earliest quasar observable in the all-sky area is $z\sim 10(11)$ for $\varepsilon_{DC}=1(0.01)$.

As standard in the literature, we provide the double power law fits (DPL) for the modeled $z\geq 3$ QLF in Table 2. The DPL fits are represented as the dashed lines in Fig. 1. The DPL function is given by

where $\phi^{*}$ is the normalization, $\alpha$ is the faint end slope, $\beta$ is the bright end slope and $M^{*}$ is the characteristic break magnitude. To calculate the best fit parameters we compute the posterior probability distribution using the Monte Carlo Markov Chain technique (MCMC; Foreman-Mackey et al. 2013) with the standard likelihood, $\mathcal{L}\sim\exp(-\chi^{2})$ over a luminosity range of $M_{UV}\in[-30,-20]$. We assume uniform priors for all parameters and enforce a DPL shape with the constraint such that $\alpha$ and $\beta$ do not overlap prior domains.

We note that the DPL fit parameters $\alpha,\beta$ and $M^{*}$ are fully consistent with each other in $\varepsilon_{DC}$, while the normalization $\phi^{*}$ is the sole parameter to vary with $\varepsilon_{DC}$. The best fit parameters suggest a steepening of the faint end slope $\alpha$ in redshift, following the evolution of the faint end slope of the HMF. In contrast, we find a very mild steepening of the bright end slope $\beta$ in redshift. However, we strongly caution the significance of this finding as the evolutionary trend in $\beta$ also correlates with an evolution of the DPL break point $M^{*}$. Since our modeled QLF is effectively a Schechter shape (from the HMF) convolved with a log-normal kernel, we can expect a fit that sets $M^{*}$ beyond the ‘knee’ of the modeled QLF will result in a steeper $\beta$ ‘fit’. Observational determinations of $(M^{*},\beta)$ seem to be consistent with this correlation (Kulkarni et al., 2019). Therefore, we find that it is difficult to associate the steepening trend of the DPL bright end slope to any physical mechanism, especially in our case where the trend is mild. Furthermore, an inspection of the QLF in Fig. 1 suggests the contrary, i.e. the bright end of the HMF appears to flatten due to the broadening from quasar luminosity scatter $\Sigma$.

With respect to the calibration QLFs, the DPL parameters derived for the $z=4$ QLF show good agreement with Akiyama et al. (2017), except for a slightly underestimated $\beta$. However, the calibrated QLF is still fully consistent with the entire set of observed data points. In contrast, our DPL parameters does not agree with the fit for the $z=5$ QLF by McGreer et al. (2018). The discrepancy in the latter is attributed to the choice of parameterization in fixing $\beta=-4.0$. We find that the DPL fits systematically overproduces UV-bright quasars compared to the results of our full modeling. This is expected as from the Schechter nature of the modeled QLF, because the bright end retains its exponential shape even after accounting for quasar luminosity scatter $\Sigma$. However, we note that the differences between model and fit are only restricted to the brightest objects and thus are still well within current observational uncertainties.