Quasar UV Luminosity Function at $3.5<z<5.0$ from SDSS Deep Imaging Data

The SDSS is an imaging and spectroscopic survey using a dedicated wide-field 2.5 m telescope (Gunn et al., 2006) with five broad bands, $ugriz$, at Apache Point Observatory. An SDSS imaging run consists of six parallel scanlines, and two interleaving runs slightly overlap, leading to duplicate observations in a small area. The imaging survey was along great circles and had two common poles, so regions near the survey poles overlap substantially. In addition, if a run (or part of a run) did not satisfy the SDSS quality criteria, the relevant region was re-observed, yielding duplicate observations in this region. Due to the above survey strategy and geometry, SDSS has a large amount of duplicate observations (referred to as overlap regions in this paper). The total area of the overlap regions is more than one-third of the SDSS footprint. Detailed information about these overlap regions can be found in Jiang et al. (2015, 2016). We selected quasars in part of the overlap regions in this paper. These overlap regions provide a unique dataset that allows us to select quasars fainter than those found from the SDSS single-epoch data.

In addition to the single-epoch main survey, SDSS conducted a deep imaging survey of $\sim 300\ \mathrm{deg}^{2}$ (Stripe 82) on the Celestial Equator in the SGC (Annis et al., 2014; Jiang et al., 2014). Stripe 82 roughly spans $20^{\mathrm{h}}<\mathrm{R.A.}<4^{\mathrm{h}}$ and $-1\fdg 26<\mathrm{Decl.}<1\fdg 26$, and was scanned around $70-90$ times. The combined data are 1.5–2 mag deeper than the single-epoch data.

There are a variety of quasar selection methods using optical imaging data. Early searches of type 1 quasars largely rely on point-like morphology and blue UV continuum colors. This color selection is efficient for quasars at relatively low redshift, as quasars and stars have different loci in color-color diagrams (Fan, 1999). Later, more methods and more sophisticated techniques were developed. For example, transfer learning has been a useful tool (Fu et al., 2021) for sky regions with large dust extinction and large contamination like the Galactic plane. Other methods including the likelihood approach (Kirkpatrick et al., 2011), the Neutral Network approach (Yèche et al., 2010), and the Extreme Deconvolution (Bovy et al., 2011) have been applied to recent quasar surveys such as the SDSS-III BOSS quasar survey (Ross et al., 2012).

Our goal was to select quasars at $3.6<z<5.5$. Searches of higher-redshift quasars in SDSS have been carried out by other programs (e.g., Fan et al., 2006; Jiang et al., 2016). We chose to use the traditional color-color diagrams to select our targets. Specifically, we selected quasars at $3.6<z<4.5$ and $4.5<z<5.5$ using the $gri$ and $riz$ colors, respectively (see Figure 1 and Table 2). At $z>3.6$ ($z>4.5$), the Ly$\alpha$ emission line enters the $r$ ($i$) band, and the Lyman forest absorption makes quasars much fainter in bluer bands. Therefore, the color-color diagrams are efficient for the selection of quasars in these two redshift ranges (Fan et al., 1999; Richards et al., 2002).

Our targets (like other targets for the SDSS BOSS survey) were selected from the SDSS Data Release 7 (DR7) imaging data. We first present our quasar selection in the overlap regions. We use ‘1’ to denote the SDSS primary detection and ‘2’ to denote the SDSS secondary detection. For example, $i_{1}$ ($i_{2}$) is the $i$-band magnitude for the primary (secondary) detection. SDSS primary detections generally have slightly higher signal-to-noise ratios than secondary detections. The selection procedure consists of two major steps. In the first major step, we retrieved a preliminary candidate list from the SDSS Query CasJobs online server. We searched the following area at high Galactic latitude: $\rm 100\arcdeg<R.A.<300\arcdeg$, $\rm Decl.>-5\arcdeg$, and Galactic latitude $b>40\arcdeg$. The object type of both primary and secondary detections is ‘star’, i.e., they were classified as point sources. Although distant quasars are point-like objects in ground-based images, faint point sources can be misclassified as extended objects. We will correct this effect in Section 3. The positional separation between a primary detection and its secondary detection was required to be smaller than $0\farcs 5$, which ensures that they are the same object. We excluded objects with the SDSS processing flags ‘BRIGHT’, ‘EDGE’, ‘SATUR’, and ‘BLENDED’. We then imposed initial color cuts to reduce the number of objects in the preliminary target list. The initial cuts are similar to (but much looser than) the final color cuts addressed in Appendix A. We do not expect to lose real quasars in this step. In addition, objects with previous spectroscopic observations were excluded using ‘specObjID’$=0$ (meaning no spectroscopic observations). This was to reduce the number of targets for follow-up spectroscopy.

After we obtained the preliminary list of targets, we combined the primary and secondary detections. For each object, we first converted its primary and secondary magnitudes to flux, and calculated weighted mean flux. Errors were added in quadrature. We then converted the combined flux and errors to AB magnitudes and errors. For example, $i$ ($i_{\rm{err}}$) denote the combined $i$-band magnitude (error). The combined magnitudes and errors will be used in the following color selection.

Our color selection criteria were primarily based on the criteria for SDSS I and II from Richards et al. (2002). The selection of candidates at $3.6<z<4.5$ ($gri$ candidates) was based on the $r-i$ versus $g-r$ diagram (the left column of Figure 1). In Figure 1, the black dots represent randomly selected point sources from the SDSS Query CasJobs online server, the blue and red dots represent a sample of randomly selected quasars from SDSS DR5 (Schneider et al., 2007), and the solid lines indicate our selection criteria. Note that there were no known quasars fainter than $i=20.2$ mag here. In order to reduce the number of contaminants, we used slightly different criteria in three different magnitude ranges, $19.0<i<20.2$ mag, $20.2<i<20.8$ mag, and $20.8<i<21.3$ mag. All selection criteria are provided in Appendix A.

The selection of quasar candidates at $4.5<z<5.5$ ($riz$ candidates) was based on the $i-z$ versus $r-i$ diagram (the right column of Figure 1). We also used slightly different criteria in the three magnitude ranges and the criteria are shown in Appendix A. These criteria are very similar to those used in the literature (e.g., Richards et al., 2002; McGreer et al., 2013; Wang et al., 2016; Yang et al., 2016).

The target selection in Stripe 82 is straightforward. The combined images and photometric catalogs were available in the archive, and the data were much deeper (Annis et al., 2014). We selected candidates down to $i=21.5$ mag and included objects brighter than $i=19.0$ mag from the SDSS Query CasJobs online server (using ‘Run’$=106$ or 206 and ‘specObjID’$=0$). The overall selection criteria are very similar to Equations A1 and A4. They are shown in Appendix A. The search area is $22^{\mathrm{h}}<\mathrm{R.A.}<4^{\mathrm{h}}$. We did not use the region of R.A. $<22^{\mathrm{h}}$ deg as the Galactic latitude becomes lower.

Our targets were observed by the BOSS spectrograph in SDSS-III (Eisenstein et al., 2011; Dawson et al., 2013). The BOSS main survey covered $\sim$10,000 deg² in the north and south galactic caps and was completed in 2014 (Alam et al., 2015). The BOSS quasar survey mainly focused on quasars at $2.2<z<3.5$. Our program was selected as one of the ancillary programs to fill spare fibers. The SDSS bitmasks used in SDSS targeting can be found in the website¹¹1https://www.sdss.org/dr12/algorithms/ancillary/boss/highz/ of our program. From the selection procedure above, we obtained 4374 quasar candidates, including 3454 candidates from the overlap regions and 920 candidates from Stripe 82. A total of 3406 candidates were spectroscopically observed. The mean fraction of targets with spectroscopic observations reaches about 78% and we will correct this incompleteness in Section 3.3.

From the spectroscopic observations, we obtained 887 quasars. They have been included in the SDSS DR16 quasar catalogue (Lyke et al., 2020). Their redshift distribution is shown in Figure 2. This sample consists of 35 quasars at $z<3$ and 852 quasars at $z>3$. The $z>3$ sample (hereafter new sample) includes 652 quasars in the overlap regions and 200 quasars in Stripe 82.

As we mentioned earlier, we did not observe the objects that had already been spectroscopically observed in SDSS I and II. Some of them also satisfy our target selection criteria. We recovered this quasar sample (hereafter the archival sample) as follows. For quasars in the overlap regions, we changed one criterion (using ‘specObjID’$!=0$) and repeated the selection procedure. For quasars in Stripe 82, we directly used the criteria to match quasars in the DR7 quasar catalogue (hereafter DR7Q, Schneider et al., 2010). We recovered a total of 346 quasars. Our final high-redshift sample is the combination of the archival and new samples, and consists of 1198 quasars at $z>3$ (Table 2).

We measure the continuum properties of the quasars assuming a power-law shape $f_{\lambda}\propto\lambda^{\alpha_{\lambda}}$. We fit this power-law to the spectral regions with little line emission. The resultant slope $\alpha_{\lambda}$ distribution is shown in Figure 3. The mean $\alpha_{\lambda}$ value is about $-1.1$, similar to previous measurements of high-redshift quasars (e.g., Fan et al., 2001c; Schneider et al., 2001), but quite softer than the results from low-redshift works due to the short wavelength coverages for high-redshift quasars (Vanden Berk et al., 2001). Continuum luminosity/magnitude $M_{1450}$ is also calculated in this step. Figure 4 shows the redshift and $M_{1450}$ distributions of the archival sample and new sample. The median value of $M_{1450}$ is around $-25.5$ mag. The new sample is about 1 magnitude deeper on average, so the QLF calculated in this paper will reach a lower luminosity compared to that from the SDSS single-epoch data. Table 3 lists our high-redshift quasar sample.

The area coverage of the SDSS overlap regions is complex. We use the Hierarchical Equal Area isoLatitude Pixelization (HEALPix; Górski et al., 2005) to estimate the effective area of the overlap regions. The basic idea is to pixelize the sky sphere into a mesh of quadrilateral pixels and the effective area is calculated by adding up all pixels that cover our data points.

For a given dataset, the starting resolution level of HEALPix is important for the area calculation. We follow Jiang et al. (2016) and adopt HEALPix Level 10 (i.e., 11.8 square arcminutes per pixel) as the best starting level for the overlap regions (see Figure 5 in Jiang et al., 2016). We classify all pixels into three categories. The first category consists of empty pixels that do not cover any objects, so they do not contribute to the effective coverage. The close neighbors to empty pixels are boundary pixels that will result in the uncertainty of the area calculation. The remaining pixels are all in the third category (hereafter non-boundary pixels). All non-boundary pixels at Level 10 contribute to the total effective area. We then gradually increase the resolution level for the boundary pixels. The resultant new pixels are once again classified into two categories, boundary and non-boundary pixels. The new non-boundary pixels are added to the total area and the new boundary pixels are refined again by increasing the resolution level. This procedure stops when the resolution roughly matches the average surface density of the data points. In this paper, we reach the best resolution at Level 12.

From the above calculation, the total area of the overlap regions is $1292\pm 266\ \rm{deg^{2}}$. The uncertainty of area coverage will be included in the measurement of the QLF. The calculation of the Stripe 82 area is very straightforward, since it is one rectangular piece of the sky. Its coverage is $225\ \rm{deg^{2}}$ with a negligible uncertainty.

In this subsection we estimate our sample incompleteness that is critical to derive QLFs. The first incompleteness is from the fact that the BOSS survey did not observe all our targets. For example, about 82% (71%) of the targets in the $gri$ ($riz$) sample at $i<20.2$ mag were observed in the overlap regions. When we correct this incompleteness, we assume that the quasar fraction in the unobserved candidates is the same as that in the observed sample. The incompleteness slightly varies with the $i$-band magnitude. This variation is considered as the uncertainty of this incompleteness, and the results (1-3%) are negligible.

The second incompleteness arises from the morphological bias. The quasar candidates that we observed are point sources, but faint point sources with low signal-to-noise ratios can be misclassified as extended sources by the SDSS photometric pipeline. We correct this bias for the targets in the overlap regions (i.e., the single-epoch data). The Stripe 82 imaging data are much deeper and we assume that the targets in this region do not suffer from the morphological bias. We will see below that this assumption is reasonable. To estimate the bias for the single-epoch data, we use 27,593 point sources classified in Stripe 82. We divide the data into narrow magnitude bins. For each magnitude bin, we calculate the fraction of the objects that are misclassified as extended sources in the single-epoch data. The resultant fractions are from 0.04 at $19.0<i<19.6$ mag to 0.55 at $21.2<i<21.3$ mag. There is a clear relation between the fraction and brightness. The fraction in the brightest range is nearly zero, suggesting little bias for the targets in Stripe 82. These fractions are considered as sample incompletenesses in individual bins, and will be included when we calculate QLFs. In order to estimate the uncertainty of the incompleteness, we resample the data 1000 times. For each time, we randomly select 500 point sources per bin to estimate the incompleteness and finally regard the standard deviation as uncertainty. The resultant uncertainties (1-2%) are negligible.

The next incompleteness comes from our color selection criteria, i.e., the color cuts introduced in Section 2. This incompleteness is described by a selection function, the probability that a quasar with a given magnitude ($M_{1450}$), redshift ($z$), and intrinsic spectral energy distribution (SED) meets the color selection criteria. We calculate the average selection probability $p_{s}(M_{1450},z)$ by assuming that intrinsic SEDs have certain distributions. Following the procedure in Fan (1999) and McGreer et al. (2013), we use a simulation tool to estimate the selection function $p_{s}(M_{1450},z)$. McGreer et al. (2013) updated the quasar spectral model in the code based on the colors of $\sim$ 60,000 quasars at $2.2<z<3.5$ from SDSS III (Ross et al., 2012). This model consists of a broken power-law continuum, most emission lines, an IGM absorption model, and an Fe emission template. It also accounts for the Baldwin effect. Jiang et al. (2016) and Yang et al. (2016) extended this model to higher redshifts with the assumption that quasar SEDs do not evolve with redshift.

Based on the model of Yang et al. (2016), we generate a sample of simulated quasars with proper photometric errors. We follow the procedure in Jiang et al. (2016) to add photometric errors. We withdraw a large representative sample of point sources in the overlap regions and Stripe 82, and derive error distributions as a function of magnitude in the $u,g,r,i,z$ bands. Finally, we construct a grid of 2,000,000 mock quasars in a redshift range of $3.5<z<5.5$ and a luminosity range of $-27.5<M_{1450}<-23.5$, with step sizes of $\Delta M=0.02$ and $\Delta z=0.02$. Then we calculate the selection function $p_{s}(M_{1450},z)$, the fraction of simulated quasars that meet our selection criteria. Figure 5 shows the selection functions for the overlap regions and Stripe 82. We estimate the uncertainty of the selection function using the same method as we did for the morphological incompleteness, and the result is around 1%. As we will see that the uncertainties of the incompleteness corrections are negligibly small compared to other uncertainties, so they are not included in the following QLF calculations.

We use a traditional $1/V_{\mathrm{a}}$ method (Avni & Bahcall, 1980) to derive the binned differential QLFs. The available volume for a quasar with absolute magnitude $M$ and redshift $z$ in a magnitude bin $\Delta M$ and a redshift bin $\Delta z$ is

where $p(M,z)$ is the final selection function that includes all incompleteness corrections discussed above.

In general, the binned QLF and its statistical uncertainty can be expressed as

where the sum is over all quasars in each bin. When a density approaches the Poisson limit, its uncertainty is corrected using Equation 7 in Gehrels (1986).

We divide our sample into several luminosity bins and redshift bins, and we focus on three redshift ranges, $3.6<z<4.0,4.0<z<4.5$, and $4.5<z<4.9$. The magnitude limits in each redshift range are determined by the faintest and/or brightest quasars and the selection functions. The binned QLF results are listed in Table 4 and also displayed in Figure 6 as the blue circles (Overlap regions) and red circles (Stripe 82). The horizontal locations of the symbols are at the centers of each magnitude bin, and the horizontal bars indicate the magnitude coverage ranges. The binned QLFs calculated for the two datasets are consistent within $1\sigma$.

We combine the two datasets from the overlap regions and Stripe 82 and derive a parametric QLF using the maximum likelihood method (Marshall et al., 1983). This method aims to minimize the function $S$, which is equal to $-2\ln L$, where L is the likelihood function:

where $p(M,z)$ includes all the incompleteness corrections discussed above. The first term is the sum over all observed quasars in the sample. The second term is integrated over the whole magnitude and redshift range of the sample. It represents the total number of expected quasars for a given luminosity function. The confidence intervals are determined from the logarithmic likelihood function using a $\chi^{2}$ distribution of $\Delta S$ ($\mathbf{=S-S_{min}}$) (Lampton et al., 1976).

We choose a double power-law form (Boyle et al., 2000) as the parametric QLF model:

where $\alpha$ and $\beta$ are the faint- and bright-end slopes, $M^{*}$ is the characteristic magnitude (or break magnitude), and $\Phi^{*}$ is the density normalization. We assume that these parameters do not change in small redshift ranges, such as the ranges considered here. We will discuss the QLF evolution in Section 4.2. We perform a grid search to determine the best-fit results and the confidence intervals. The grid resolutions of log$\Phi^{*}$, $M^{*}$, and $\alpha$ are 0.05, 0.05, and 0.1, respectively. There is a strong degeneracy between $M^{*}$ and $\alpha$, so we set a bright limit of $-28.0$ mag for $M^{*}$. The best-fit results are listed in Table 5.

Figure 6 shows the results in three redshift ranges. The open circles denote the data points (some of the faintest bins) that have very low completeness and significantly deviate from the general trend. It is unclear what causes this deviation. This has been frequently seen in previous studies and is likely due to some unknown selection effects. We did not use these data points in the above calculation. Our sample covers a limited range of luminosity, so it is not able to constrain both slopes $\alpha$ and $\beta$. Therefore, we use three fixed values for $\beta$ in each redshift range (see Figure 6) and derive the other three parameters. The best-fitted $\alpha$ values are about $-1.8$ at $3.6<z<4.9$, indicating that the results for different $\beta$ values and different redshift ranges are not significantly different. In addition, most of the best-fitted $M^{*}$ values for three redshift ranges are lower than $-27$ (see Table 5), making the double power-law model degenerating into a single power-law model. These results suggest that our sample alone is not enough to constrain all parameters in the above QLF model. In the next section, we will combine our binned QLFs with some results in the literature.

In Figure 7 we show a collection of previous QLF measurements at $3.6<z<4.9$ (Richards et al., 2006; McGreer et al., 2013; Yang et al., 2016; Boutsia et al., 2018; McGreer et al., 2018; Schindler et al., 2019; Kulkarni et al., 2019; Kim et al., 2020; Boutsia et al., 2021). All data points have been scaled to $z=3.8,4.25$, or $4.7$ using the density evolution model of Schindler et al. (2019) with $\gamma=-0.38$ (e.g., $\mathbf{\Phi(z=3.8)=\Phi(z)*10^{-\gamma(z-3.8)}}$). Richards et al. (2006) constructed a sample of 15,343 quasars at $0<z<5$ in 1,622 $\mathrm{deg}^{2}$ from SDSS DR3, with a small fraction of quasars at $z>3.6$. Our sample is roughly 1.5 mag deeper than this sample. Boutsia et al. (2018) focused on faint quasars at $3.6<z<4.2$ in the COSMOS field. Schindler et al. (2019) built a sample of very luminous quasars at $2.8<z<4.5$ to constrain the bright-end slope $\beta$ of the QLF. Boutsia et al. (2021) identified 58 bright quasars at $3.6<z<4.2$ and the brightest ones reach $M_{1450}=-29.5$ mag. Kulkarni et al. (2019) combined multiple datasets from previous surveys and generated a sample of more than 80,000 quasars, with a small fraction of quasars at $z>3.6$. They used this sample to derive QLFs and study quasar evolution from $z=7.5$ to 0. In the bottom panel of Figure 7, we particularly compared our results with previous QLF measurements at $z\sim 5$ (e.g., McGreer et al., 2013; Yang et al., 2016; McGreer et al., 2018; Kim et al., 2020). We did not include samples with no or very few spectroscopic observations (e.g., Akiyama et al., 2018; Niida et al., 2020).

Figure 7 shows that our results are generally consistent with previous measurements. In the top panel, our luminosity coverage for $z=3.8$ partly overlap with the luminosities covered by Richards et al. (2006), Boutsia et al. (2018), and Kulkarni et al. (2019). In this overlap range, the binned QLFs from different studies roughly agree with each other. Our binned QLFs at $-26.5<M_{1450}<-25.5$ are about 1.5 times the results of Richards et al. (2006). It is unclear whether their results were underestimated or our results were overestimated. In the middle panel for $z=4.25$, our result is well consistent with the previous results from Richards et al. (2006) and Kulkarni et al. (2019) except two high data points from Kulkarni et al. (2019). The bottom panel shows several studies of the QLF at $z=4.7$ and most of these results are consistent with ours within 1 $\sigma$ level. It is worth noting that discrepancies are relatively larger at the faint- and bright- ends where uncertainties are also significantly large.

In order to better constrain the shape of QLFs, we combine our QLF measurements with the results from some previous studies (Richards et al., 2006; McGreer et al., 2013; Yang et al., 2016; Boutsia et al., 2018; McGreer et al., 2018; Schindler et al., 2019; Kulkarni et al., 2019; Kim et al., 2020; Boutsia et al., 2021). We caution that these previous samples are not fully independent, because the same quasars were used in some samples (thought analyzed with different methods). For the previous results, detailed selection functions for individual samples are not publicly available, so we directly use their binned QLFs. The derived parameters from binned data may be subject to some biases (e.g., La Franca & Cristiani, 1997; Page & Carrera, 2000). To reduce the impact from potential biases, we only use our results in the luminosity ranges that our sample covers. In addition, we exclude those points with very low completeness. The points that are not used in the fitting are shown as the open symbols in Figure 7 and the purple line in each panel denotes the best-fit QLF. We fit the observed QLF data points ($\Phi_{\mathrm{obs}}$) using the maximum likelihood estimation. We use the logarithmic likelihood function $\mathcal{L}$:

where $\sigma=\sqrt{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{extra}}^{2}}$, $\sigma_{\mathrm{obs}}$ is the 1$\sigma$ uncertainty of $\Phi_{\mathrm{obs}}$ from the literature, and $\sigma_{\mathrm{extra}}$ is calculated from the magnitude bin size by the error propagation. The details of $\sigma_{\mathrm{extra}}$ and the chosen priors on the parameters are shown in Appendix B. We use the emcee Python package ²²2https://emcee.readthedocs.io/en/stable/ (Foreman-Mackey et al., 2013) for the Markov chain Carlo sampling of the QLF parameters. The best-fit results and uncertainties are estimated based on the 16th, 50th and 84th percentiles of the samples in the marginalized distributions. Figure 7 shows the best-fit QLFs and parameters are listed in Table 5.

The best-fit values at $z=3.8$ are $\alpha=-1.7^{+0.2}_{-0.2}$, $\beta=-4.0^{+0.4}_{-0.5}$, and $M^{*}=-26.7^{+0.4}_{-0.4}$. The $\beta$ and $M^{*}$ values are consistent with the results of Boutsia et al. (2018) ($\beta=-4.025^{+0.575}_{-0.425}$ and $M^{*}=-26.50^{+0.85}_{-0.60}$) with reduced errors. The binned QLF data points of Boutsia et al. (2018) are higher than previous measurements, resulting in a steeper faint-end slope $\alpha$ and a higher density $\Phi^{*}$ than our results. The QLF at $z=4.25$ is also well constrained by our results. As mentioned earlier, the faintest binned QLF calculated by Kulkarni et al. (2019) is very high, so they obtained a steep faint-end slope $\alpha=-2.20^{+0.16}_{-0.14}$. Our slope $\alpha$ is slightly flatter.

For the QLF at $z=4.7$, we adopt the results with $\alpha=-1.7^{+0.6}_{-0.3}$ and $\beta=-3.2^{+0.5}_{-1.2}$ as our best fit. We combined the data from Kim et al. (2020) and McGreer et al. (2018) for the faint end and the data from McGreer et al. (2013) and Yang et al. (2016) for the bright end. In this high-redshift range, the current sample size is still small and thus the above measurements are associated with relatively large uncertainties.

Based on the measurements from the above subsection, we explore quasar evolution at $3.5<z<5$. We use a linear function to describe the evolution of the four QLF parameters in this small redshift range:

where $X\in\{\mathrm{log}\Phi^{*},\ M^{*},\ \alpha,\ \beta\}$. Here we consider four cases:

1. $\circ$

   Case 1: A PDE model where only $\Phi^{*}$ evolves, i.e., we fit 5 parameters: $\mathrm{log}\Phi^{*}_{0},\ M^{*}_{0},\ \alpha_{0},\ \beta_{0},\ k_{\Phi}$.

2. $\circ$

   Case 2: We allow $\Phi^{*}$ and $\beta$ to evolve, i.e., we fit 6 parameters: $\mathrm{log}\Phi^{*}_{0},\ M^{*}_{0},\ \alpha_{0},\ \beta_{0},\ k_{\Phi},\ k_{\beta}$.

3. $\circ$

   Case 3: A LEDE model. We allow $\Phi^{*}$ and $M^{*}$ to evolve, i.e., we fit 6 parameters: $\mathrm{log}\Phi^{*}_{0},\ M^{*}_{0},\ \alpha_{0},\ \beta_{0},\ k_{\Phi},\ k_{M}$.

4. $\circ$

   Case 4: We allow all parameters to evolve, i.e., we fit 8 parameters: $\mathrm{log}\Phi^{*}_{0},\ M^{*}_{0},\ \alpha_{0},\ \beta_{0},\ k_{\Phi},\ k_{M},\ k_{\alpha},\ k_{\beta}$.

We fit these 4 models to the observed QLF data points using the Maximum likelihood estimation introduced earlier. We calculate the reduced $\chi^{2}$ as below,

where $n$ is the number of data points and $\nu$ is the number of the free parameters. The results are listed in Table 5.

The models of Cases 1, 2, and 3 have similar fitting performance in terms of $\chi^{2}_{\nu}$. The model in Case 4, with up to 8 free parameters, has likely overfitted the observed data because its $\chi_{\nu}^{2}$ is much smaller than 1. Besides, the PDE model in Case 1 has the smallest number of free parameters, and is enough to describe the evolution of the QLF at $z\sim 3.5-5$. This is consistent with the conclusion of Kim & Im (2021).

We compare the four models in Figure 8. The Case 1 result (purple solid line) and the Case 4 result (purple dashed line) have similar fitting performance, so it is difficult to find the best-fit model based on current data. These two results are both higher than the result of Kulkarni et al. (2019) in the bright end. This discrepancy can be attributed to the different data that two studies used. For example, Kulkarni et al. (2019) did not use the data from Boutsia et al. (2018) and Schindler et al. (2019). These data may increase the measurement of the bright-end QLF. In addition, the redshift coverage of Kulkarni et al. (2019) is much larger than this work. The data in the range of $z<3.5$ and $z>5$ will affect the result in $3.5<z<5$ when computing the QLF evolution.

Compared with Kim & Im (2021), our QLF slopes are steeper. The measurements of the slopes sensitively depend on the data points at the brightest and faintest ends that usually suffer from large incompleteness and uncertainties. In addition, the combination of different samples introduces extra uncertainties that are often difficult to characterize. A large sample with a full coverage of the both ends is needed to improve the measurement of the QLF.

Finally, we explore the cumulative spatial density evolution of quasars at high redshift. The spatial density of quasars brighter than a given magnitude $M$ is calculated by integrating the QLF,

where we use the PDE model (case 1) of this work as $\Phi(M,z)$. Figure 9 shows the cumulative density as a function of redshift for different magnitude ranges using our model and previous results (Richards et al., 2006; Ross et al., 2013; McGreer et al., 2013; Jiang et al., 2016; Yang et al., 2016; Wang et al., 2019; Kim & Im, 2021). It indicates a rapid density decline from $z\sim 3.5$ to $z\sim 7$, consistent with previous studies. For example, Fan et al. (2001a) fit an exponential decline, $\rho(<M,z)\propto 10^{kz}$, and found $k=-0.47$ at high redshift. McGreer et al. (2013) found that the slope at $4<z<5$ was $k=-0.38\pm 0.07$. Previous results also suggested that the spatial density of quasars drops faster with increasing redshift at $z>3.5$. As shown in Figure 9, the slopes at $5<z<6$ and $6<z<7$ are $k=-0.72\pm 0.11$ and $k=-0.78\pm 0.18$, respectively (Jiang et al., 2016; Wang et al., 2019). The PDE model from Kim & Im (2021) (dashed lines) also supports this scenario as shown in Figure 9.

From our sample, we measured $k=-0.7\pm 0.1$ at $3.5<z<5$, which is slightly steeper than previous measurements for the same redshift range, but is similar to those measured at $z=5\sim 7$. The main reason for the steeper slope is that our binned QLF within $-26.5<M_{1450}<-25.5$ at $z=3.8$ is about 1.5 times the previous results (see grey squares in Figure 7 and blue crosses in Figure 9). However, the cumulative densities between our model and the observed results are consistent within 1$\sigma$. To confirm our results, we need a larger and more complete sample, such as a quasar sample from the Chinese Space Station Telescope wide-area slitless spectroscopic survey (Zhan, 2021). If confirmed, the quasar density at $3.5<z<5$ declines faster than previous measurements, and as fast as the density evolution at $z>5$.