The Pan-STARRS1 $\mathbf{z>5.6}$ Quasar Survey: III. The $\mathbf{z\approx6}$ Quasar Luminosity Function

We present the $z\!\approx\!6$ type-1 quasar luminosity function (QLF) based on the Pan-STARRS1 (PS1) quasar survey. The PS1 sample includes 125 quasars at $z\approx5.7-6.2$ with $-28\lesssim M_{1450}\lesssim-25$. Complemented by 48 fainter quasars from the SHELLQs survey, we evaluate the $z\approx6$ QLF over $-28\lesssim M_{1450}\lesssim-22$. Adopting a double power law with an exponential evolution of the quasar density ($\Phi(z)\propto10^{k(z-6)}$; $k=-0.7$), we use a maximum likelihood method to model our data. We find a break magnitude of $M^*=-26.38_{-0.60}^{+0.79}\,\text{mag}$, a faint end slope of $\alpha=-1.70_{-0.19}^{+0.29}$, and a steep bright end slope with $\beta=-3.84_{-1.21}^{+0.63}$. Based on our new QLF model we determine the quasar comoving spatial density at $z\!\approx\!6$ to be $n( M_{1450}<-26)=1.16_{-0.12}^{+0.13}\,\text{cGpc}^{-3}$. In comparison with the literature, we find the quasar density to evolve with a constant value of $k\approx-0.7$ from $z\approx7$ to $z\approx4$. Additionally, we derive an ionizing emissivity of $\epsilon_{912}(z=6) =7.23_{-1.02}^{+1.65}\times 10^{22}\,\text{erg}\,\text{s}^{-1}\text{Hz}^{-1}\text{cMpc}^{-3}$ based on the QLF measurement. Given standard assumptions and the recent measurement of the mean free path of Becker et al. (2021) at $z\approx6$ we calculate an HI photoionizing rate of $\Gamma_{\text{HI}}(z{=}6) \approx 6 \times10^{-16}\,\text{s}^{-1}$, strongly disfavoring a dominant role of quasars in hydrogen reionization.


INTRODUCTION
Quasars are rapidly accreting supermassive black holes (SMBHs) at galaxy centers, which shine as the most luminous non-transient light sources in the Universe. At low redshift tight correlations between the SMBH mass and its host galaxy's central properties raised attention on the role of AGNs in galaxy evolution (see Kormendy & Ho 2013, for a review). Specifically, feedback during bright quasar phases has been identified as a prominent avenue to establish this relationship (e.g., Di Matteo et al. 2005). In this context, understanding the evolution of quasars has received growing attention in recent years, especially their evolution in the early Universe. Following the discovery of the first z 6 quasars  it was quickly realized that SMBHs with masses of M BH ≈ 10 9 M already exist less than 1 Gyr after the Big Bang, placing constraints on their formation. These have been significantly tightened by the discovery of the most distant quasars at z ≈ 7.5 known today (Bañados et al. 2018b;Yang et al. 2020;Wang et al. 2021). While single sources highlight the open questions with regard to SMBH formation (see Inayoshi et al. 2020, for a review), understanding the full demographics at high redshifts will play a key role in addressing them with the quasar luminosity function (QLF), the main observational statistic to characterize their population.
At lower redshifts (e.g., Boyle et al. 1988Boyle et al. , 2000Pei 1995), the QLF is most effectively described by a broken double power law, which has been widely adopted at arXiv:2212.04179v2 [astro-ph.GA] 9 Feb 2023 higher redshifts as a good representation of most quasar samples (e.g., Richards et al. 2006;Shen & Kelly 2012;Ross et al. 2013). In this context the QLF is described by a break magnitude, a normalization, and two powerlaw slopes. At z ≈ 6 the most comprehensive measurement of the QLF was presented in Matsuoka et al. (2018) using a combined sample built from the Subaru High-z Exploration of Low-Luminosity Quasars (SHEL-LQs; Matsuoka et al. 2016) project and previous QLF analyses Willott et al. 2010) with a total of 112 sources, covering a redshift range of 5.7 < z < 6.5 and luminosities of −30 M 1450 /mag −22. Placed in context with the lower redshift QLF literature (e.g., Richards et al. 2006;Croom et al. 2009;Glikman et al. 2011;Shen & Kelly 2012;Ross et al. 2013;Akiyama et al. 2018;Schindler et al. 2018Schindler et al. , 2019Boutsia et al. 2021;Pan et al. 2022) the recent results underline the (exponential) increase (Schmidt et al. 1995;Fan et al. 2001a) in quasar activity from z ≈ 6 to its peak at z = 2 − 3 Kulkarni et al. 2019;Shen et al. 2020). The highest redshift constraint on the QLF at z ≈ 6.7 (Wang et al. 2019a) indicates an even more rapid decline of quasar activity at z > 6.5 with consequences for upcoming quasars surveys at z > 8 (e.g., based on the Euclid mission wide survey; Euclid Collaboration et al. (2019); Scaramella et al. (2021)).
Quasars, or more generally active galactic nuclei (AGN), and star formation are the major sources of ultra-violet (UV) radiation that drive the reionization of intergalactic hydrogen. The role of quasars in this process, as inferred from QLF number counts, has been a matter of debate in the literature. Type-1 UV QLFs at z ≈ 5 (e.g., McGreer et al. 2013McGreer et al. , 2018Yang et al. 2016;Kim et al. 2020;Shin et al. 2020), at z ≈ 6 (e.g., Jiang et al. 2008;Kim et al. 2015;Jiang et al. 2016;Matsuoka et al. 2018), at z ≈ 6.5 (Wang et al. 2019a), and from the redshift compilation of Kulkarni et al. (2019) provide substantial evidence for a subdominant contribution of quasars compared to star formation to reionization at z 5. These results are supported by the analysis of the bolometric QLF based on multi-wavelength data sets (Shen et al. 2020). On the other hand, Giallongo et al. (2015), Giallongo et al. (2019), and Grazian et al. (2020) find high quasar number densities for lower luminosity sources, −22.5 ≤ M 1450 ≤ −18.5. These studies are based on multi-wavelength selected sources from the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS) GOODS-South, GOODS-North, and EGS fields. The analysis in Giallongo et al. (2015) and Giallongo et al. (2019) is largely based on photometric candidates, whereas the study of Grazian et al. (2020) has spectroscopy for their two sources.
Based on these number densities at the faint-end the authors argue that quasars could be the dominant source of ionizing photons at z ≈ 4 − 6, which is supported by the analysis of Grazian et al. (2022) based on the QUBRICS quasar survey (Calderone et al. 2019;Boutsia et al. 2020). However, a range of independent deep X-ray studies report significantly lower number densities for the faint quasar population (e.g., Weigel et al. 2015;Cappelluti et al. 2016;Vito et al. 2016;Ricci et al. 2017;Parsa et al. 2018), challenging the results of Giallongo et al. (2015) and Giallongo et al. (2019).
In this work we present a new measurement of the type-1 UV QLF at z ≈ 6 based on the selection strategy and discoveries from the Pan-STARRS distant quasar survey ). Our quasar sample includes 125 sources at z ≈ 5.7 − 6.2 within a luminosity range of −28 M 1450 −25, more than doubling the number counts of previous samples in this range (Jiang et al. 2016, SDSS). Combining the new sample with lower luminosity sources from SHEL-LQs (Matsuoka et al. 2018), we present the most precise measurement of the type-1 UV QLF at these redshifts to date.
In Section 2 we review the quasar selection of the PS1 distant quasar survey and present the new quasar sample used in this work. Section 3 discusses the resulting quasar selection function and completeness. We present the QLF in Section 4 and discuss the implications with regard to quasar evolution and reionization in Section 5. Finally, we summarize this work in Section 6. Interested readers can find the mathematical framework for our QLF analysis described in detail in Appendix A, whereas Appendix B expands on the discussion of our quasar model used for the completeness calculation. In this work we adopt a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7. All magnitudes are reported in the AB photometric system.

DATA
The foundation of the quasar searches this work builds upon is the 3π Steradian Survey (Chambers et al. 2016) that was carried out by the Panoramic Survey Telescope and Rapid Response System Telescope #1 (PS1, Kaiser et al. 2002Kaiser et al. , 2010. From 2009 to 2015 the PS1 3π survey imaged the sky above a declination of −30 • in the five filter bands g P1 , r P1 , i P1 , z P1 , y P1 . The full data releases of the PS1 3π survey are hosted by the Barbara A. Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute 1 .

The PS1 3π Steradian Survey PV2 catalog
The PS1 3π distant quasar survey  based their selection on the internal prereleases of the stacked PS1 3π photometry, in particular the second internal data release PV2. According to Bañados et al. (2016, Section 2.1) the PV2 5σ median limiting magnitudes are (g P1 , r P1 , i P1 , z P1 , y P1 ) = (23.2, 23.0, 22.7, 22.1, 21.1). The quasar selection uses the stacked PSF magnitude signal to noise ratio (SNR x ) and the (3σ) limiting magnitude (m lim, x ) of a given band x as selection criteria. These properties are derived from quantities in the PV2 catalog. We briefly describe them here to provide context for Section 3.1.2, where we discuss the modeling of simulated quasar photometry to assess the selection function completeness. The SNR x is calculated from the stacked PSF magnitude 1σ uncertainty (σ m,x ) of the same filter band x: The 3σ limiting magnitude m lim, x is derived from the band zero point (zp x ) and the 1σ uncertainty on the stacked PSF fit instrumental flux (PSF INST FLUX SIG): m lim,x = −2.5 × log 10 (3 × PSF INST FLUX SIG x ) + zp x . (2) The PSF INST FLUX SIG is a PV2 catalog property from the internal data release, but not provided in the public data release on MAST. The band x zero point, zp x , only depends on the stacked exposure time (EXPTIME x ) in that filter band: zp x = 25 + 2.5 × log 10 (EXPTIME x ) .
All stacked magnitudes are corrected for Galactic extinction using the Schlegel et al. (1998) dust map with the corrections of Schlafly & Finkbeiner (2011). The quasar photometric selection is conducted on the dereddened stacked magnitudes.

Quasar candidate selection
We follow the quasar selection criteria of Bañados et al. (2016, Section 2.1.1), focusing on the search of quasars at 5.7 z 6.2. The selection criteria are applied to the PS1 PV2 catalog generated from PS1 image stacks.
For completeness, we briefly describe the selection here. First, we exclude sources which have been flagged as suspicious (see Bañados et al. 2014, Table 6) by the Imaging Processing Pipeline (Magnier et al. 2020a,b). Additionally, we require 85% of the normalized pointspread-function (PSF) flux in the i P1 , z P1 , and y P1 bands to be located in unmasked pixels (PSF QF > 0.85). This quality cut leans towards a more complete selection by including some lower-quality measurements (Magnier et al. 2020b). We will refer to these requirements on the PV2 catalog as our "photometric quality selection" hereafter.
The Milky Way plane has traditionally been avoided by quasar surveys using criteria based on Galactic latitudes (e.g., Fan 1999;Jiang et al. 2016). The main reason is that the high source density and stronger Galactic extinction leads to unreliable photometry for extragalactic background sources. Following Bañados et al. (2014) we impose a Galactic latitude limit of |b| > 20 • and additionally select only sources with modest degrees of Galactic reddening as determined from the Schlegel et al. (1998) dust map cross-matched to the PS1 PV2 source catalog, our "extinction selection" criterion: Additionally, we exclude all sources around M31 (7 • < R.A. < 14 • ; 37 • < Decl. < 43 • ) as their inclusion results in a large number of candidates that are most likely stars associated with M31. Figure 1 shows a Mollweide projection of the sky with the PS1 quasar survey coverage shown in gray. Confirmed high redshift quasars selected as described in this section are highlighted as orange diamonds. A description of the resulting survey area is provided in Section 3.3, where we discuss the different contributions of the selection criteria to the survey selection function. Apart from M-, L-, and T-dwarf stars (or brown dwarfs), the main contaminants for high-redshift quasars are low redshift galaxies, which mostly appear extended under the PS1 observing conditions. The PS1 3π website 2 lists median seeing conditions of 1.11 , 1.07 , and 1.02 for the i P1 , z P1 , and y P1 filter bands, respectively. To reject extended sources we adopt the "morphology criterion" as discussed in  Section 2.1). We keep sources whose absolute difference between the aperture and PSF magnitudes, |f ext |, is below a value of 0.3 in either the PS1 z -or y-band: Selected PS1 quasars Figure 1. Mollweide projection of the PS1 quasar survey area (20803 deg 2 ) considered for the QLF analysis. Right ascencsion and Declination are noted in hours and degrees. Regions included in the quasar selection based on the PS1 PV2 catalog are shaded grey using a HEALPix tesselation with Lvl 10 (12,582,912 HEALPix cells over the entire sky, see Table 2). The M31 mask is included (RA 0h 42m 44s, Dec +41 • 16' 9 ). High redshift quasars (Section 2.2) are shown as filled orange diamonds.  tested this criterion against known spectroscopic stars and galaxies in SDSS (DR12; Alam et al. 2015) and quasars (SDSS DR10 quasar catalog Pâris et al. 2014). This criterion removes the majority (92%) of galaxies, while retaining 92% of stars and 97% of quasars. The focus of this quasar luminosity function analysis is the redshift range 5.7 z 6.2. At these redshifts quasars can be efficiently differentiated from brown dwarfs by applying color criteria to the PS1 i -, z -, and y-band stacked magnitudes . We summarize the "photometric selection" criteria discussed in Bañados et al. (2014) for the 5.7 z 6.2 range below.
These criteria are applied to the stacked, dereddened PS1 magnitudes. After all the above selection criteria are applied to the PS1 PV2 photometry, we perform forced photom-etry on the PS1 stacked and single-epoch images (see Section 2.2 and 2.3 in Bañados et al. 2014). We remove all sources where the forced photometry is inconsistent with the reported values in the PS1 PV2 catalog. Effectively this removes 80% of the candidates mainly due to discrepancies in the i P1 -band , their Section 2.2).
After the steps above, which are all automatic, this yields a total of 1032 candidates. The PS1 photometry of the z -band and y-band has been taken close in time to each other. Therefore, we visually inspect their stacked and single-epoch images to exclude bright spurious sources (e.g., moving objects) that appear in only one single-epoch z -band y-band pair. EB and JTS visually inspected all 1032 candidates independently and assigned a rank from 1 (good photometry) to 4 (inconsistent/erroneous photometry). Sources that are clearly detected in the z-band, have consistent single epoch measurements, and a i − z color in agreement with the catalog are given rank 1. Sources with rank 2 usually have a fainter detection in the z-band with their stacked y-band measurement close to the the limit of our SNR requirement, such that the single epoch images are harder to assess. Sources with rank 3 usually have some issues with the data. For example, a clear zband detection is only visible in 1 or 2 epochs. Sources that are given rank 4 show clear data artifacts. The final candidates catalog includes sources for which the summed rank is 2, 3, or 4. This results in a total of 640 quasar candidates for follow-up observations, 73 sources with summed rank 2, 117 sources with summed rank 3 and 450 sources with summed rank 4. All sources in the final candidate catalog received an individual rank of either 1 or 2 from either EB or JTS. We exclude a total of 202 sources with summed visual rank 5 and 190 sources with an even higher summed visual rank (> 5) from further follow-up.
The full quasar selection procedure can be summarized in seven individual steps: Each of these selection steps has an impact on the selection function of the survey, which we discuss in Section 3.3.

The PS1 high redshift quasar sample
The detailed follow-up strategy along with descriptions of the photometric and spectroscopic data is presented in Bañados et al. (2014; . Of the 640 good quasar candidates, 41 are published quasars in the literature and 274 sources have been photometrically or spectroscopically followed up. The photometric follow-up rejected 79 candidates. These sources were ruled out as good candidates due to their red y P1 − J color (y P1 − J > 1) or if the followup photometry did not meet the main selection criteria (Equations 6-12). Our confirmation spectroscopy identified 84 quasars among the 195 observed sources. Due to limited telescope time 325 quasar candidates have not been followed up, yet. We display the full identification statistics in Figure 2. The upper panel shows histograms of the confirmed quasars (N CQ (y P1 )), the rejected candidates (N RC (y P1 ); photometry and spectroscopy) and the total number of candidates. We estimate the identification efficiency, Eff ID (y P1 ), as the ratio between confirmed candidates to the total number of identified candidates, We express the identification completeness, S ID (y P1 ), as the number of observed quasars divided by the number of expected quasars given our efficiency: To calculate a continuous distribution of the completeness and efficiency as a function of y P1 , we use a Gaussian kernel density estimate. Kernel density estimation (KDE) provides a non-parametric representation of the probability density function of a random variable. In comparison to binning it is smooth and independent of the end points of the bins. For the full sample we determined the best bandwidth via cross-validation to be around 0.1. We will use this value for all KDE estimates related with the identification completeness. In the follow-up campaigns sources were prioritized by visual rank. This leads to a bias in the completeness between samples of different rank (see Figure 17 in the Appendix). To mitigate the bias we calculate the expected number of quasars for each rank individually, accounting for the different efficiencies of each sample. We sum the number of expected quasars for each of the three samples (ranks 2,3, and 4) and then calculate and adjusted identification completeness: We show the efficiency and the adjusted identification completeness in the bottom panel of Figure 2. The selection efficiency of the full sample is highest at the bright end of the candidate distribution and then declines towards fainter magnitudes. The low number statistics at the faint end (faintest bin the top panel) result in an upturn of the efficiency at the faint end. Our sample has been followed up with very high identification completeness as the bright end (y P1 < 20.0). Towards fainter magnitudes the completeness declines with two visible minima at y P1 = 20.75 and y P1 = 21.3, where it reaches only 50% and 30%, respectively. At the faint end, y P1 > 21.3, our follow-up becomes more complete again and the identification completeness increases. The full quasar sample used in our QLF analysis consists of the 125 confirmed PS1-selected quasars. Their sky distribution is shown in Figure 1 and we provide a complete list in the Appendix (  The analysis of the QLF requires us to quantify our quasar selection function (see Section 2.2), including a correction for the incomplete spectroscopic follow-up (Section 2.3). In order to realistically evaluate the photometric selection for type-1 quasars we simulate quasar photometry taking into account the properties of the PS1 PV2 catalog (e.g., inhomogeneous depth). In particular, to apply the photometric selection criteria from Equations 6-12, we need to produce observed (errorprone) magnitudes, the signal-to-noise ratio and the limiting magnitudes. We describe the photometric modeling in Section 3.1 and continue to derive K-correction terms based on these models to estimate absolute magnitudes with the QLF quasar sample in Section 3.2. The selection function is then evaluated in Section 3.3.
3.1. Modeling PS1 quasar photometry 3.1.1. Modeling quasar photometry with simqso To build a sample of simulated quasar photometry we are using a forked version of the python package simqso 3 (McGreer et al. 2021), originally presented in McGreer et al. (2013McGreer et al. ( , 2018, which contains updates to the emission line defaults and includes photometric error models for additional surveys. The code constructs artificial quasar spectra from a parametric model of different type-1 quasar spectral components. The parameters are informed from the general knowledge of quasar spectra at all redshifts and follow the assumption that quasar spectral energy distributions do not evolve with redshift (e.g., Jiang et al. 2006;Shen et al. 2019;Yang et al. 2021). The quasar model is built from a powerlaw continuum, quasar emission lines, an iron pseudocontinuum component, a dust component, and a component modeling absorption of neutral hydrogen in the intergalactic medium. Using simqso we have chosen the parameters of the different components to produce a reliable high-redshift quasar model specifically developed for this project.
We construct the continuum from a set of broken power laws (f ν = A × ν α ) that are designed to reproduce the continuum emission of the Selsing et al. (2016) quasar template and the broad band photometry of known SDSS quasars (Schneider et al. 2010;Pâris et al. 2017). Break points and mean power law slopes are listed in Table 1. For each spectrum the mean power law slopes are sampled by drawing values from a Gaussian distribution around the mean slope (α) with a dispersion of σ = 0.3. We are modeling a total of 65 individual lines  ity, i.e., the Baldwin effect (Baldwin 1977), are taken into account. To approximate the emission from the large range of iron transitions seen in quasar spectra we included a composite of iron templates. This composite was last updated in simqso by Yang et al. (2016) and uses the Vestergaard & Wilkes (2001) iron template at 1250Å-2200Å, the Tsuzuki et al. (2006) iron template at 2200Å-3500Å, and the Boroson & Green (1992) iron template at 3500Å-7500Å. Following Lyu & Rieke (2017), we add three blackbody components with temperatures of 1800 K, 880 K, and 285 K to model dust emission. We have adjusted the equivalent width scaling of some emission lines, the amplitude scaling of some regions of the iron template, and the amplitudes of the blackbody dust emission components in order to fully reproduce the Selsing et al. (2016) quasar template and the optical to infrared mean quasar colors from the combined sample of the SDSS DR7 (Schneider et al. 2010) and DR12 (Pâris et al. 2017) quasar catalogs. To model absorption due to neutral hydrogen in the intergalactic medium, we use the Lyman-α forest model of McGreer et al. (2013). We use this model to construct 10,000 different quasar sight lines up to z = 7 that are randomly paired with simulated quasars to produce the absorption signatures blueward of the quasar's Lyman-α line.
The simulation that we use as a basis for the Kcorrection to absolute magnitudes and the selection function analysis consists of a uniform grid in absolute magnitude at 1450Å, −30 ≤ M 1450 ≤ −23, and redshift, 5.0 ≤ z ≤ 7.0, with 56 intervals along the axis of absolute magnitude and 40 intervals along the redshift axis. Each cell is uniformly populated with 88 quasars for a total of 197,120 quasars. The resulting quasar spectra are then multiplied with the PS1 filter bands to produce synthetic magnitudes.
Recent studies (Yang et al. 2021;Bischetti et al. 2022) report an increased BAL quasar fraction at z 6 compared to lower redshifts. Similar to previous completeness calculations for z > 6 QLF measurements Wang et al. 2019a;Matsuoka et al. 2018), our quasar model does not take BAL quasars into account. If our quasar selection would be biased against BAL quasars, it would be necessary to account for their different spectral properties in the completeness calculation. Around 40% (6) of the 14 BAL quasars in the Bischetti et al. (2022) sample have weak absorption (balnicity index BI < 1000 kms −1 ), which does not affect the quasar colors beyond the scatter included in our quasar model. On the other hand, strong broad C iv and Si iv absorption troughs (BI > 5000 kms −1 ) will reduce the flux between the C iv and the Lymanα line significantly. At z > 5.8 the Si iv line will be in the y PS1 -band and the C iv line will lie redward of that. Hence, the strongest BAL features will diminish the y PS1 -band flux, resulting in a bluer z PS1 − y PS1 color. As a consequence, our color selection requiring i PS1 − z PS1 > 2.0 and z PS1 − y PS1 < 0.5 also includes strong BAL quasars naturally without the need to add  them to our quasar model. This inherently high completeness to BAL quasars in our color selection may play a role in the higher BAL fractions discovered at z > 6 ( Yang et al. 2021;Bischetti et al. 2022). There is evidence that quasars with weak emission lines (or weak-line quasars; WLQ) could make up ∼ 10% of the population at z 6 (e.g., Shen et al. 2019), which is higher than what is observed at lower redshifts (Diamond-Stanic et al. 2009). The weaker Lyα and N v emission lines can result in bluer i PS1 − z PS1 color and a redder z PS1 − y PS1 color for z > 5.6 quasars. In the recent literature on z > 6 QLFs (e.g., Jiang et al. 2016;Wang et al. 2019a;Matsuoka et al. 2018) a larger fraction of WLQs has not been taken into account, except by Matsuoka et al. (2018) who down-scaled the Lyα line of their model spectra to be more inclusive of quasars with weaker emission lines.
The spread in quasar spectral properties covered within the random sampling of simqso already covers a significant fraction of emission line variation. Nonetheless, quasars with very weak lines (e.g., the blue line in Fig. 10 of ) are underrepresented in our quasar model (and in all models used for previous QLFs) and may be missed at z > 6.0 due to their redder z PS1 − y PS1 (see Figure 16 left). This effect introduces a small systematic uncertainty of < 10% to our QLF analysis, implying that we might be missing < 10% of quasars at 6.0 z 6.2 due to their weak emission lines.
One of the main assumptions for our determination of the selection function is the simqso quasar model itself.
To test the robustness of this model we compare it to the quasar model presented by Temple et al. (2021, here-after qsogen). Euclid Collaboration et al. (2019) used qsogen to predict the quasar yields in the Euclid survey. Both simqso and qsogen are designed to reproduce the median colors of quasars at z < 5. The z P1 − y P1 quasar color of both our simqso model and qsogen agrees well with each other. The random sampling of continuum and emission line properties in the simqso model fully encompasses the variation in z P1 − y P1 color introduced by the spread in emission line properties parametrized by the qsogen emline type parameter. The main difference is the treatment of absorption due to the neutral intergalactic medium, which leads to small differences in the z P1 − y P1 color at z > 5.8 and larger differences in the i P1 − z P1 color over z = 5.2 − 6.6. In comparison to PS1 quasar colors the simqso model colors provide a more faithful representation of the empirical data points, whereas the default qsogen model results in a bluer i P1 − z P1 color. We refer interested readers to a more detailed discussion in Appendix B.

Simulating PS1 photometric uncertainties
The simqso package allowed us to generate synthetic quasar photometry. However, to fully assess the selection function, we need to take the photometric errors into account. For surveys with an approximately homogeneous depth one can derive simple magnitude error relations to sample the photometric error from. However, as shown in Figure 3 the PV2 photometric z-band 3σ limiting magnitude, (Equation 2), has strong variations depending on the sky position. This is also true for all other PS1 bands (e.g., see Figure 19 in Appendix C) and the depth variations do not necessarily correlate with each other across the different filter bands. As a consequence we cannot define a single magnitude error relation to apply to all simulated quasar photometry. Instead we adopt a sampling approach where we randomly associate observed catalog properties with simulated quasars that allow us to calculate the error properties necessary to evaluate the photometric selection criteria on our simulated sample. Our method is similar in spirit to the approach discussed in Section 3.1 of Yang et al. (2016), which was also adopted for PS1 in .
The methodology allows us to calculate the stacked PSF magnitude uncertainty (σ m ), which also provides us with the signal-to-noise ratio (see Equation 1) and the 3σ limiting magnitude (m lim ) for each simulated quasar. Using σ m we can then construct PSF magnitudes from the synthetic photometry and evaluate the photometric selection function from the simulated quasar sample, including photometric uncertainties.
We begin with the PV2 catalog after we applied the quality flags, the extinction selection and the morphology selection. This guarantees that we are using properties of point sources most similar to our quasar candidates that also have good quality photometry within our chosen footprint. In order to sample these properties we use the Hierarchical Equal Area isoLatitute Pixelation (HEALPix 4 ; Górski et al. 2005) as implemented in the healpy (Zonca et al. 2019) python package to associate each source in our clean catalog with a small area on the sky, i.e., a HEALPix cell. For this purpose we subdivided the sky into a map of 196,608 (Lvl=7) cells, of which 100,766 are filled with at least one source. We then sample source positions as uniformly as possible across the survey footprint and retrieve their filter band zero points and their Galactic reddening value, E(B − V ), from the Schlegel et al. (1998) dust map for each of the 197,120 simulated quasars. Each of the simulated quasars is then randomly associated with the noise properties of a real source in the catalog.
We have used the full cleaned PV2 catalog to investigate the relations between σ m , the zp, and the stacked PSF magnitude (CAL PSF MAG). Figure 4 (left) shows a 2D histogram of the median σ m,z as a function of zp z and CAL PSF MAG z . The figure shows the expected dependence of σ m,z on CAL PSF MAG z and zp z , which is further highlighted by the slightly diagonal lines of constant SNR. We note that within each of the 2D bins a range of σ m,z values exists as depicted in Figure 4 (right). In order to approximate the stacked PSF magnitude error for a simulated quasar we first use the sampled E(B−V ) to redden the synthetic magnitude. The reddened magnitude and the sampled zp then determine the bin in the σ m (zp, CAL PSF MAG) 2D histogram (e.g., Figure 4, left) from which we randomly draw a σ m value given the 1D σ m distribution in that particular bin (e.g., Figure 4,right). This is done for all PS1 filter bands. In rare instances where the combination of the sampled zp and simulated CAL PSF MAG does not exist in the cleaned PV2 catalog, we adopt the maximum value of the median σ m in the 2D histogram. This can be the case for synthetic magnitudes in dropout bands (e.g. the PS1 g-band) that exceed the observed range of values. From σ m we calculate the SNR for each given band using Equation 1. Once we have associated a σ m for all synthetic quasars and all bands, we perturb the reddened synthetic magnitude by drawing from a Gaussian distribution with the magnitude as the mean and the associated σ m as the 1σ uncertainty. In a last step we subtract the reddening from the perturbed magnitudes to retrieve the appropriate dereddened magnitudes used for the photometric selection function evaluation. In order to assess the 3σ limiting magnitude m lim (see Equation 2) we not only require knowledge of zp, but also of the stacked PSF fit instrumental flux uncertainty, PSF INST FLUX SIG, at the source position. Similar to σ m , PSF INST FLUX SIG can also be mapped as a function of zp and CAL PSF MAG as shown in Figure 5. The figure shows the median z-band PSF INST FLUX SIG per bin and depicts that it is also dependent on both zp and CAL PSF MAG. In order to associate a PSF INST FLUX SIG 5.6 5.8 6.0 6.2 6.4 Redshift z Figure 6. Luminosity dependent K-correction factor derived from the simulated quasar sample to correct the observed dereddened PS1 stacked PSF y-band magnitude, yP1, to obtain the absolute monochromatic magnitude measured at rest-frame 1450Å, M1450. The colored lines show the Kcorrection factor for different luminosities as a function of redshift. The luminosity dependence is most pronounced at z ≈ 5.5 − 5.8, where the C iv emission line falls into the PS1 y-band.
value with a simulated quasar we proceed just as we have done for the magnitude error above. Then we calculate the limiting magnitudes for all simulated quasars according to Equation 2.

K-correction and absolute magnitudes
In order to evaluate the QLF as a function of absolute monochromatic magnitude measured at rest-frame 1450Å (M 1450 ), we need to be able to calculate M 1450 from the observed dereddened PS1 stacked PSF magnitudes for our quasar sample. The conversion to restframe M 1450 requires a term that accounts for the changing filter response as the quasar spectrum is being redshifted, the K-correction. We compute a K-correction based on our simulated quasar photometry as a function of M 1450 and quasar redshift. As the PS1 z -band cuts out at around 9300Å, the rest-frame wavelength of 1450Å falls into the PS1 y-band above z ≈ 5.4. Therefore, we derive relations between M 1450 , z, and the dereddened y-band magnitude y P1 from our grid of simulated quasars. We then retrieve the y P1 Kcorrection factor by interpolation from this relation. Figure 6 shows the luminosity dependent K-correction factor as a function of redshift for different quasar luminosities. The median absolute deviation (MAD) of the K-correction factor over the entire grid, a measure of its dispersion due to variety in different quasar spectra simulated by the code, is K MAD = 0.07 mag. We use this K-correction to calculate M 1450 from y P1 for all quasars in our QLF sample. The M 1450 values for each individual quasar are listed in Table 7.
Ideally, one would like to measure the flux at 1450Å directly from the quasar spectra to derive M 1450 . Unfortunately, not all quasars in our sample have the necessary spectral coverage and therefore we conclude that our photometry-based M 1450 determination provides a more general approach and allows us to treat all quasars in our sample uniformly.

The selection function
In order to understand the completeness of our quasar selection as a function of absolute magnitude M 1450 and redshift z we need to quantify the selection function following the strategy laid out in Section 2.2. The full survey selection function S(q) can be written as the product of the selection functions for all independent selection steps. The parameter vector q contains all the catalog properties, which are used in the different selections.
We have separated the selection function for the extinction cut and area exclusion S Ω , the selection function based on morphology S morph (Equation 5), the photometric selection S phot (Equations 6-12), the quality flag criteria S qf , and the selection introduced by our incomplete follow-up S ID . Before we present the full survey selection function, we will discuss their individual contributions.

The survey area estimate and Galactic extinction
The Galactic extinction E(B − V )(α, δ) is a function of right ascension α and declination δ of the sources. We apply a selection criterion, Equation 4, to construct our sample, which effectively excludes sky regions with high Galactic extinction. Thus, the extinction selection function only reduces the available survey area irrespective of redshift and luminosity of the sources: We already noted above (Section 2.2) that in addition to the extinction cut, regions close to the Galactic plane (|b| < 20 deg), as well as the area around M31, are excluded for the final quasar selection as well. The complete area selection function is then a product of the three criteria above: To estimate the area included by our selection criteria we again utilize a HEALPix tesselation of the sky. Our area estimate is based on 328,479,702 sources of the PS1 PV2 catalog, which passed the photometric quality criteria, the morphology selection and were included in the area defined by the extinction selection and the Galactic plane and M31 exclusion regions. Using HEALPix we divide the sky into a grid of curvilinear, equal-sized quadrilaterals. We calculate the number of HEALPix cells which are populated by at least one source and sum up the area of all of those cells for a total area estimate. We vary the HEALPix cell sizes to understand the resolution effects on the total estimated area. At the lowest resolution, the sky is represented by 12 HEALpix cells. For the next resolution level these cells are each divided into four sub-cells. Hence, the total number of HEALpix cells N pix depends on the resolution level lvl following N pix = 12 × 4 lvl . At too low resolution levels the survey area is overestimated as the coarse cells cannot capture the fine structure of the Galactic extinction map. As we proceed to finer resolutions the number of HEALpix cells approaches the number of sources for the area estimate. At even higher resolution the size of the HEALpix cells gets smaller than the area between adjacent sources in the survey and we are effectively undersampling the area. As a result the survey area estimate decreases rapidly. Table 2 shows the results from the HEALpix analysis. The columns are the HEALpix level Lvl, the total number of pixels N pix per level, the number of filled pixels N pix, filled per level, the area per pixel Ω pix and the total filled area per level Ω total , the effective survey area. The total area estimate drops significantly above lvl = 11, indicating that the HEALPix cell density is approaching the source density and we are beginning to overresolve the area. Thus, we adopt lvl = 10 for our fiducial area estimate and use the differences to the adjacent levels to reflect the uncertainties on our estimate. This results in a survey area of 20803.38 +61.75 −54.45 deg 2 . This calculated survey area now implicitly takes the extinction selection function into account. We consider the relative uncertainty on the survey area of ≈ 0.3% to be negligible and it is not propagated further.

Morphology selection function
Submillimeter observations indicate that the host galaxies z 5.7 quasars are often compact with effective (half-light) radii of R e ≈ 1.11 kpc . They are effectively unresolved by the PS1 photometry. The selection criterion in Equation 5 is designed to reject extended contaminants from the selection. The quasar completeness was quantified in Bañados et al. (2016, their Section 2.1 and their Figure 3) to be 97%. With the significant detections we require in the PS1 z -and y-band (Equations 6 and 8), it is reasonable to assume that this value is independent of apparent magnitude and quasar redshift. The host galaxies of quasars become only prevalent at much fainter magnitudes (Matsuoka et al. 2016, their Figure 6) beyond the PS1 detection limit. Therefore, we adopt the value of 97% for the quasar selection completeness: Excluding extended sources introduces a bias to our sample against strongly lensed quasars. Multiple images or the foreground lens galaxy could make the lensed quasar appear extended in the imaging.

Photometric selection function
The criteria described in Equations 6-12 form the core of the PS1 photometry 5.7 z 6.2 quasar selection. The photometric selection function is based on the magnitudes and SNRs in different PS1 filter bands and the limiting i -band magnitude, . Interpolated photometric selection function evaluated on a redshift (z) and absolute magnitude (M1450) grid of simulated quasars. Black contours are drawn at a completeness of 20%, 50%, 80% and 90%. The high completeness at 5.7 z 6.2 directly reflects the color selection criteria listed in Section 2.2. The faint end is limited by the depth of the PS1 survey and the signal-to-noise ratio criteria in different filter bands.
We apply the selection criteria to our grid of simulated synthetic quasar photometry (Sections 3.1.1 and 3.1.2) to evaluate its impact as a function of redshift and absolute magnitude. We present the resulting photometric selection function in Figure 7. High redshift quasars are commonly selected by the strong flux break at the Lymanα line, where blueward emission is absorbed by neutral gas in the IGM. The i P1 − −z P1 color criterion (Equation 9) selects quasars with a Lymanα break above redshifts of z ≈ 5.6 and is responsible for the rise of the selection function at this redshift. The z P1 − y P1 color criterion (Equation 9) imposes a certain level of continuity between the two filter bands and thus actively deselects quasars beyond z ≈ 6.3. The width of the redshift transition regions depends on the diversity in quasar spectral properties and the variations of absorbing neutral hydrogen along the line of sight. The signalto-noise ratio criteria (Equations 6 and 8) limit the selection in apparent magnitude and results in a redshift dependent absolute magnitude limit. The inhomogeneity of the PS1 depth and the intrinsic scatter of the flux measurements result in a slow decrease with increasing absolute magnitude rather than a sharp break.

Photometric quality selection
We assess the selection function for these photometric quality criteria (see Section 2.2) using two empirical samples of quasars matched to the PS1 PV2 catalog.
As we want to address the completeness of our method, we first build a sample of z ≈ 5.7, which were discovered from other surveys but lie within the PS1 footprint. After applying the signal-to-noise ratio requirements on the photometry (SNR zP1 > 10, SNR yP1 > 5) to ensure an appropriate comparison, we retained 76 sources. We applied the quality criteria and find 72 sources to be included in the photometric quality selection, resulting in a completeness of ∼ 95%.
To verify the completeness with a larger sample we perform an additional test on quasars selected from the SDSS DR16 quasar catalog (Lyke et al. 2020). We require the SDSS quasars to have low extinction (E(B − V ) < 0.3) and follow our candidate signal-tonoise ratio requirements (SNR zP1 > 10, SNR yP1 > 5). Additionally, we only select sources within the redshift range of 2 < z < 5. We exclude quasars at z > 5 from the DR16 quasar catalog as the majority of these sources are objects misclassified by the automatic pipeline. Catalog quasars at z < 2 have been preemptively excluded to avoid complications that arise when the host galaxy starts to be resolved. At that point the values of PSF QF and the quality flags potentially differ significantly from pure point sources. The full sample contains 179,945 SDSS quasars, of which 162,770 are retained when applying the photometric quality criteria, resulting in an average completeness of ∼ 90%. We further investigate whether the photometric quality selection function is dependent on the dereddened PS1 z -band magnitude. We thus calculate the completeness in 34 magnitude bins between z P1 = 18 to 21.5 with a bin width of ∆z P1 = 0.1. The result shows a dependence of the completeness with z P1 (see Figure 8). The quality selection function dependency on the yband magnitude reflects the correlation between some quality flags (e.g., POORFIT, MOMENT SN, see Bañados et al. (2014), Table 6) and the lower SNR of fainter sources. We model the binned measurements with a hyperbolic tangent: and use the LMFIT python package (Newville et al. 2014) to retrieve the best-fit parameters via the Levenberg-Marquardt algorithm. The values are a = 0.052 ± 0.003, b = 0.904±0.002, η = −1.049±0.127, and φ = −20.161± 0.063. The best-fit model function is shown in Figure 8 as the solid orange line. We adopt this parameterization of the photometric quality completeness for the calculation of the full selection function. For each point in the space of absolute magnitude M 1450 and redshift z, we evaluate the photometric quality selection S(z P1 ) by calculating the apparent magnitude z P 1 (M 1450 , z) at the point using the adopted cosmology and a z P 1 K-correction factor determined analogously to the y P 1 K-correction (Section 3.2). The photometric quality selection in the space of redshift and absolute magnitude is shown in the left panel of Figure 9. The overall completeness is high, > 85%.

Identification completeness
We further have to take into account that not all of our 640 quasar candidates have been identified through follow-up observations or the literature. We note that spectroscopic follow-up observations were prioritized by the candidate's visual rank. However, no other selection relevant information e.g., the coordinate range, played a factor in this process. To account for this we use the weighted identification completeness as a function of y P1 as shown in the bottom panel of Figure 2. Similar to the photometric quality completeness we evaluate the kernel density estimate of the identification completeness Comp ID (y P1 ) in the space of redshift (z) and absolute magnitude (M 1450 ) by mapping these properties to the dereddened PS1 y-band magnitude y P1 (M 1450 , z) using the adopted cosmology and the y P 1 K-correction factor. The right panel of Figure 9 shows the map of the spectroscopic completeness as a function of redshift and absolute magnitude. The low identification completeness (< 80%) at apparent magnitudes y P1 = 20.5 − 21.5 of our follow-up campaigns (Figure 2) is reflected here at absolute magnitudes M 1450 ≈ −24.9 to M 1450 ≈ −26.2. The completeness rises towards the faint end, equivalent to the behavior in Figure 2, as we have identified many of the faintest candidates.

The full survey selection function
In order to obtain the full survey selection function, we combine the individual selection functions multiplicatively according to Equation 16 and present the result in Figure 10. The shape of the selection function is dominated by the photometric selection (see Figure 7). Our incomplete follow-up identification decreases the completeness at the faint end (M 1450 = −26.5 to −25) compared to the pure photometric selection. The impact of the photometric quality selection and morphology selection is more subtle as it decreases the completeness over the full absolute magnitude and redshift range. We use this selection function to correct for our completeness in the measurements of the quasar luminosity function. Figure 10 highlights a few quasars that lie in regions of very low completeness and it is worthwhile to briefly understand why they passed our selection strategy. The lowest redshift quasar of our sample is PSO224.65067+10.21379 at a spectroscopic redshift of z ≈ 5.4. This source shows strong broad absorption in N v that effectively removes all Lyα flux, mimicking a Lyα break at z ≥ 5.7. At the high redshift end our sample includes SDSS J1030+0524 at z = 6.308. This source as well as some fainter sources at z ≈ 6.2 (PSO184.33893+01.52846 and PSO334.01815-05.00488) have especially strong Lyα flux. Due to the resulting blue z P1 − y P1 color they are included in our selection (Equation 9). Finally, our sample includes a number of faint quasars (M 1450 −25.5) with selection probabilities below 20%. This indeed indicates that we expect more than 5 times as many quasars in this parameter region in full agreement with previous QLF results Matsuoka et al. 2018).

Completeness
Full selection function (interpolated) Figure 10. The full PS1 distant quasar survey selection function for the quasar selection focused at 5.7 z 6.2. We show the 125 confirmed quasars that make up the QLF sample as blue diamonds. Black contours are drawn at a completeness of 10%, 20%, 50% and 80%.
In this section we present our QLF measurement at z ≈ 6. As the PS1 distant quasar sample only constrains the bright end of the QLF (M 1450 < −25), we combine our sources with 48 quasars from the SHELLQs quasar sample presented in Matsuoka et al. (2018). We adopt the quasar properties based on their Table 1 and use their selection function (see their Figure 9; the electronic data provided by Yoshiki Matsuoka) to correct the SHELLQs sample for incompleteness. We discuss the PS1 QLF binned as a function of absolute magnitude in Section 4.1. Then, in Section 4.2 we fit a double power law to the combined sample (PS1 + SHELLQs) using a maximum likelihood approach.

The binned quasar luminosity function
We determine the binned quasar luminosity function over the redshift interval 5.65 z 6.25 in bins of absolute magnitude at 1450Å. This redshift range includes a total of 121 of the total 125 quasars in our sample. For the calculation of the binned QLF we implemented the 1/V a method (Schmidt 1968;Avni & Bahcall 1980) with the modifications outlined in Page & Carrera (2000). We correct the quasar number counts using the completeness from the full quasar selection function de-scribed above (Section 3.3.6). Table 3 summarizes the results. The columns describe the absolute magnitude bin, the median absolute magnitude med(M 1450 ), the median redshift med(z), the number of quasars per bin, and the uncorrected and corrected binned QLF values (Φ) with uncertainties reflecting the confidence interval for a Poisson distribution that corresponds to 1σ in Gaussian statistics. 5 The results of our binned PS1 (and SHELLQs) QLF is depicted as (thin) filled orange diamonds in Figure 11 compared to other binned QLF measurements at z ∼ 6 from the literature. We compare our results to the binned QLFs determined by Willott et al. (2010) (41 quasars, 5.74 < z < 6.42, light grey hexagons), Jiang et al. (2016) (47 quasars, 5.7 < z ≤ 6.4, dark grey squares), and Matsuoka et al. (2018) (112 quasars, 5.7 ≤ z ≤ 6.5, blue circles). With a total of 121 quasars, the binned PS1 QLF agrees well with the literature data at z ≈ 6. We note that our sample covers a narrower and slightly lower redshift range than the previous work in the literature. As a consequence we miss the bright quasar J0100+2802 at z = 6.30 (Wu et al. 2015 This demonstrates the consistency of our methods and their implementation, but inadvertently results in the blue circles and orange diamonds to completely overlap in Figure 11. While we use only the SHELLQs data, Matsuoka et al. (2018) derive the binned QLF from a combination of samples, which explains the differences seen between the thin diamonds and the blue circles.

Maximum Likelihood Estimation of the quasar luminosity function
The measurement of the binned QLF, while agnostic to its underlying shape, is dependent on the choice of binning, both in luminosity and redshift. With our choice of only one redshift bin, the analysis in the previous section could not account for any redshift evolution of the QLF. Alternatively, we can assume a parametric model for the QLF, including redshift evolution, and constrain the model QLF Φ(M, z|Θ QLF ) and its parameters Θ QLF by Markov Chain Monte Carlo (MCMC) sampling from the probability of the model QLF given the observed quasar sample N (M, z), P QLF ≡ P (N (M, z)|Φ(M, z|Θ QLF )). We follow Mar- 5 We approximately calculate this confidence interval equivalent to Equations 4 and 5 in Gehrels (1986).
shall et al. (1983) to derive the logarithmic probability ln(P QLF ). We present the full derivation of the logarithmic probability for a luminosity function model Φ(M, z|Θ QLF ) with a selection function S(q) in Appendix A. The logarithmic probability ln(P QLF ) can then be approximated as (see also Equation A16) where Λ(M, z) is the quasar incidence rate as given by Equation A6. Given a large enough dynamic range in luminosity, the QLF at low-redshift is well approximated by a broken double power law (DPL) (e.g., Boyle et al. 2000), defined by the normalization Φ , the break magnitude M and the two power law slopes α and β. By convention, α is most commonly chosen as the faint-end slope, with β then describing the bright-end slope. Generally, all four parameters could evolve with redshift (see e.g., Kulkarni et al. 2019). Our sample only spans a narrow redshift interval, therefore we only adopt a redshift evolution for the normalization in the form of log(Φ (z)) = log(Φ (z = 6)) + k × (z − 6) .
The parameter k describes the exponential evolution of the quasar density with redshift. We implement the luminosity function model and the calculation of the logarithmic probability and then use the python package emcee  Matsuoka et al. (2018), where the same value for k was used. The covariance matrix of the fit parameters is shown Figure 18. We note that the SHELLQs quasar sample covers a larger redshift range than ours. This is taken into account in our maximum likelihood formulation given our assumption on the redshift evolution.  Table 4 along with their uncertainties (16 to 84 percentile range). The first two columns of the table specify the model and data used in the fit. Table 5 summarizes selected QLF studies from the literature for comparison. The figure highlights two characteristics of our new QLF DPL fit. The bright-end slope is significantly steeper than previous measurements and the overall number densities are lower with the orange curve lying beneath the blue and grey curves for the majority of the magnitude range.
To test the robustness of our fiducial QLF fit results we explore four variations on the model. The first variation mimics our fiducial model but allows k to vary during the fit. The results are in good agreement with our fiducial model. However, the best fit value for k is −0.20 ± 0.2 in tension with our assumption of k = −0.7. We will continue to discuss this inconsistency below. In a second variation we fit a single power law (SPL) with k = −0.7 to the PS1 quasar sample without the additional SHELLQs data. The PS1 data alone can be well described by this single power law with a slope of α = −2.81 and a normalization of log Φ * (z = 6) = −8.85. Figure 12 shows the SPL fit (yellow) against our DPL fiducial model (grey). The third variation uses our fiducial model, but only includes PS1 quasars with M 1450 ≤ −26 in addition to the SHELLQs sample. It is designed to exclude the region of lower spectroscopic completeness (see right panel in Figure 9) to understand its influence on the fit results. All fit parameters change slightly within the uncertainties of the fiducial model. Comparing this variation (teal) to the fiducial model (grey) in Figure 12 underlines that the differences over the constrained magnitude range are negligible. We conclude that the magnitude range of low spectroscopic completeness has a minor influence on our fiducial results. Due to the restricted redshift range of the PS1 quasar sample we miss two bright quasars included in Jiang et al. (2016) and Matsuoka et al. (2018). These are J1148+5251 at z = 6.42 with M 1450 = −27.8 and J0100+2802 at 6.30 with M 1450 = −29.1. As a test we artificially add these two sources with a nominal redshift of z = 6 to our quasar sample and fit it again with our fiducial model. We denote this variation as bright end test in Figure 12 and Table,4. The artificial inclusion of these two bright sources significantly changes all QLF parameters, which is reflected in the different shape of this variation (orange) to our fiducial model (grey) in Figure 12. Most notably the bright end slope changes to β = −3.12 and the break is now 0.8 mag fainter. This test highlights the strong influence individual sources at the extreme bright end can have on the QLF measurement. With the conclusions from the different fit variations in mind, we now discuss our fiducial fit results in context with the current literature. Our QLF model favors a break magnitude of M 1450 = −26.38. This value is brighter than previous work at z ≈ 6 Jiang et al. 2016;Matsuoka et al. 2018), but also fainter than work at z ≈ 5 (Yang et al. 2016;McGreer et al. 2018). It is also significantly fainter than the study of Kulkarni et al. (2019), which model the QLF from z = 0 to 6 and fit a break magnitude of M 1450 ≈ −29 at z = 6, effectively constraining only the faint end slope with the data at z = 6. The differences to the results of Jiang et al. (2016) and Matsuoka et al. (2018) are at least in part due to the inclusion of bright quasars at z > 6.25 in these samples as our bright-end test (Table 4, last row) highlights.
Allowing for k to vary or excluding the fainter PS1 sample data with the high spectroscopic incompleteness from our fit (see second and fourth row in Table 4) does not change our best-fit break magnitude significantly, given its uncertainties. Rather, it seems to be a robust result given the combined PS1 + SHELLQs quasar sample (also see Figure 18).
The break magnitude M 1450 and density normalization Φ (z = 6) values are highly covariant in a broken double power law ( Figure 18). As a natural consequence of our brighter break magnitude we measure a lower value of Φ (z = 6) than the previous studies at z = 6. Taking this covariance into account, our best-fit QLF model agrees well with previous determinations especially.
Our fiducial best fit returns a bright-end slope with a relatively value of β = −3.84, significantly steeper than the literature data at z = 6 (β ≈ −2.8 Willott et al. 2010;Jiang et al. 2016;Matsuoka et al. 2018). The exception is the study of Kulkarni et al. (2019). Follow-  ing their global QLF fit, they find a very steep bright end slope of β = −5.05 +0.76 −1.18 . However, their work does not yet include the faint SHELLQs quasars at z 6. In consequence, they derive an extremely bright bestfit break magnitude at z = 6, M * ≈ −29, such that the data they use only constrains the faint end slope at this redshift. We emphasize that the PS1 quasar sample that determines β covers a narrower redshift range. Therefore, our samples do not include some very luminous quasars at z > 6.25, e.g., J1148+J1148+5251 or J0100+2802. Artificially including these sources in our sample with an assumed redshift of z = 6 changes the best-fit model results significantly (Table 4, last row). We conclude that these few sources at very bright magnitudes are the main driver between differences between our fiducial measurements and the studies of Jiang et al. (2016) and Matsuoka et al. (2018). On the other hand, a steep bright end slope is not uncommon. In fact, at z ≤ 5 a range of studies (e.g., Richards et al. 2006;Mc-Greer et al. 2013;Yang et al. 2016) find a bright end slope of β −3, with some studies reporting an even steeper slope of β ≈ −4 (Schindler et al. 2019;Boutsia et al. 2021;Pan et al. 2022). Viewed in this context, our results at z ≈ 6 indicate that the bright end slope is generally steep (β ≈ −4) and does not evolve significantly with redshift.
The faint end slope measurement of our maximum likelihood fit is largely determined by the SHELLQs quasar sample. Our best fit value of α = −1.70 lies between the value of α = −1.23 measured by Matsuoka et al. (2018) and the previous determination of α = −1.9 by Jiang et al. (2016). The addition of the SHELLQs sample to our fit explains the flatter slope compared to Jiang et al. (2016) as data at these faint luminosities was not available at the time. The 2σ differences to the result of Matsuoka et al. (2018) is mostly driven by the influence of bright sources (Table 4, last row) on the QLF results, which are present in their sample but not included in ours. These sources significantly affect the resulting break magnitude, which is covariant with the faint end slope (see Figure 18). Additionally, the QLF of Matsuoka et al. (2018) includes data from Willott et al. (2010) not present in our analysis.
Our fiducial QLF fit assumes an exponential density evolution with k ≈ −0.7. As a test we remove this assumption and allow k to vary, resulting in a best-fit value of k = −0.20 ± 0.2. This value is in tension with our assumption, but does this mean that our assumption of k ≈ −0.7 is not justified? We have based the assumption on literature data, and the following discussion on the quasar density redshift evolution in Section 5.1 strongly supports our assumed value of k for the fiducial fit. Hence, we conclude that the redshift range, which is limited by our sample selection, is not large enough to probe the quasar density evolution sufficiently. The bright-end quasar density n(M 1450 < − 26) has been known to increase from z ≈ 5 to z ≈ 2 (Schmidt et al. 1995). Fan et al. (2001a) modeled the quasar density up to z ≈ 6 with n(M 1450 < − 26) ∝ 10 k(z−z ref ) , finding a value of k = −0.47. McGreer et al. (2013) and Jiang et al. (2016) reported and even steeper increase of the luminous quasar density with k = −0.7, the factor assumed in our maximum likelihood analysis above. From z ≈ 7 to z ≈ 6 the quasar density is reported to increase even more steeply with k = −0.78 ). This finding is supported by quasar searches from the VIKING survey (B. Venemans, private communication) that go beyond the first discoveries (Ven-  Figure 13. Redshift evolution of the QLF from z ≈ 3 − 7 colored by redshift. Filled symbols show the binned QLF at these redshifts from different studies in the literature (Ross et al. 2013;McGreer et al. 2018;Schindler et al. 2019;). This comparison visualizes the strong decline in quasar number densities across redshift z ≈ 3 to 7.
emans et al. 2013). Going backwards in cosmic time, this seemingly accelerating decrease in quasar density has important consequences for the predicted number of discoverable quasars at even higher redshifts, z > 8. In Figure 13 we To determine the quasar density at the bright end, we integrate our best-fit QLF at z = 5.85 down to M 1450 = −26. This results in a value of n(M 1450 < −26) = 1.16 +0.13 −0.12 cGpc −3 , in line with the other measurements at z ≈ 6 Jiang et al. 2016;Matsuoka et al. 2016). We chose a redshift of z = 5.85 close to the median redshift of the PS1 quasar sample (see Table 3), which determines this magnitude range of the QLF. We show our result in comparison with values from the literature in Figure 14. We include a range of studies that provide a significant re-evaluation of the QLF at z = 3 − 5 (Akiyama et al. 2018;Schindler et al. 2018;Giallongo et al. 2019;Schindler et al. 2019;Boutsia et al. 2021;Onken et al. 2022) compared to the first results from SDSS Ross et al. 2013;Shen & Kelly 2012). We also show the redshift evolution from integrating the quasar density from the QLFs of Richards et al. (2006, z = 0.3−5) and Kulkarni et al. (2019, z = 0.3 − 6) as grey dotted and dot-dashed lines. We note that Niida et al. (2020) and Kim et al. (2020) also provide updated measurements on the z ≈ 5  Furthermore, we include the exponential density evolution with k = −0.78 (Wang et al. 2019a) and k = −0.7 anchored on our value as grey dashed and solid lines in Figure 14. Pan et al. (2022) find a single value of k = −0.7 to describe the density evolution across z = 3.5 to z = 5. Based on our results we argue that this evolution continues to z ≈ 7, when excluding the discrepant data from Boutsia et al. (2021) and Giallongo et al. (2019). There is evidence that at z < 4 the density evolution flattens (k > −0.7) before the turnover point at z ≈ 2 − 2.5 Kulkarni et al. 2019) as also discussed in Onken et al. (2022). In light of the recent literature and given the systematic uncertainties inherent in QLF estimates due to differing models for completeness correction, we conclude that the bright end density evolution at z > 4 can be well described by an exponential decline with k = −0.7. Comparing the work of Matsuoka et al. (2018) and our new estimate of the z = 6 quasar density with the value from Wang et al. (2019a), we do not find evidence for a more rapid decrease of the quasar density at z > 6.5 as originally reported from the comparison with the Jiang et al. (2016) quasar density in Wang et al. (2019a).

Forecasting high redshift quasar detections
We explore the impact of our new quasar luminosity function measurement on future high redshift quasar detections. For this purpose we use the Euclid mission (expected launch 2023) as our main example. The Euclid wide-area survey will deliver Y -, J -, and H -band photometry down to a 5σ limiting magnitude of 24.0 mag (for point-like sources) over ∼ 15, 000 deg 2 . We use simqso to simulate quasars and their photometry in the Euclid bands at 7 ≤ z < 10 according to our QLF over the wide-area survey footprint. The Lyman-α break at 1215Å enters Euclid's J -band at z 8.9. Therefore, we require only H < 24.0 for a detection in the Euclid wide-area survey. We show the resulting detection number counts from our best-fit QLF model in comparison with the QLF models of Jiang et al. (2016), Matsuoka et al. (2018), and  in Table 6. The numbers drop to 5-8 at the highest redshift bin, where our new QLF measurement provides a more optimistic forecast. This simple forecast does not claim to be a comprehensive prediction of the quasar yields based on quasar selection strategies as, for example, presented in Euclid Collaboration et al. (2019). In this context our predicted detection numbers should be regarded as an upper limit to the number of quasars that could be discovered with Euclid depending on the selection strategy. Our aim here is to simply illustrate how our new QLF measurement impacts our expectations for quasar discoveries. As the Euclid Collaboration et al. (2019)

Quasar contribution to hydrogen reionization
Based on our new measurement of the QLF at z ∼ 6 we calculate the quasar contribution to the HI photoionization rate of the UV background. Following the literature (e.g., Haardt & Madau 1996, 2012; Faucher-Giguère 2020) the HI photoionization rate is where σ HI (ν) is the frequency dependent HI photoionization cross section and n ν (ν, z) is the number density of ionizing photons per unit frequency at redshift z. The lower boundary of the integral ν 912 corresponds to the frequency at a wavelength of 912Å. At z = 6 we can assume that the optical depth of ionizing photons is smaller than unity, τ eff ≤ 1, allowing us to adopt the 'local source' approximation (e.g., Zuo & Phinney 1993;Madau et al. 1999), simplifying n ν (ν, z) to In the equation above l(ν, z) is the mean free path of ionizing photons and (ν, z) is the comoving emissivity of ionizing sources. We closely follow Shen et al. (2020) in adopting the frequency dependence on the mean free path based on the results of Faucher-Giguère et al. (2008), l(ν, z) = l(ν 912 , z)(ν/ν 912 ) 3(β−1) , where the power law index of the intergalactic HI column density distribution is assumed to be β = 1.5 (Madau et al. 1999). Furthermore, we also assume a power law shape for the extreme UV quasar continuum with an index of α UV = 1.7 (Lusso et al. 2015), Assuming a frequency dependence of σ HI ∝ ν −3 we analytically integrate Equation 27, which yields Γ HI (z) ≈ (1 + z) 3 3 + α UV − 3(β − 1) 912 (z)l(ν 912 , z)σ HI (ν 912 ) .  . We compare our data to various results in the literature Akiyama et al. 2018;Matsuoka et al. 2018;McGreer et al. 2018;Parsa et al. 2018;Kulkarni et al. 2019;Wang et al. 2019a). We include the individual data points from Kulkarni et al. (2019), that were not affected by systematic errors as discussed by the authors. Horizontal error bars on indvidual data points indicate the redshift ranges of the different QLF samples. The blue solid line and the blue shaded area show the derived posterior median emissivity evolution model and 1σ uncertainties of Kulkarni et al. (2019). We also display the models by Haardt & Madau (2012) and Madau & Haardt (2015). Our derived quasar emissivity falls well below other measurements at z = 6 with the exception of the work from Matsuoka et al. (2018) with which we share the SHELLQs quasar sample at the faint end.
For the HI photoionization cross-section we use a value of σ HI (ν 912 ) = 6.35 × 10 −18 cm 2 (Verner et al. 1996;Becker et al. 2015a). Using our new measurement of the QLF we first calculate the ionizing emissivity of quasars at 1450Å, Here we assume that the escape fraction of ionizing photons from the type-1 quasar population measured by the QLF is unity. We adopt an upper integration boundary (faint limit) of M 1450 = −18 for comparison with the recent literature (Kulkarni et al. 2019;Matsuoka et al. 2018;. Assuming a power-law SED for the quasars in the extreme UV (Lusso et al. 2015), We estimate the ionizing emissivity at 912Å as Based on our fiducial double power law fit to the combined PS1+SHELLQs quasar sample, we calculate a value of 912 (z = 6) = 7.23 +1.65 −1.02 × 10 22 ergs −1 Hz −1 cMpc −3 . The errors reflect the statistical fit uncertainty on the QLF as corresponding to the 16% to 84% percentile range. Figure 15 shows our result in comparison with values from the recent literature. Our quasar sample is dominated by the SHELLQs quasar sample at the faint end, which strongly affects the ionizing emissivity. In comparison with other studies of the type-1 UV QLF at z = 6, this fact largely explains the disagreement with the values of Jiang et al. (2016), Parsa et al. (2018) and Kulkarni et al. (2019) and the agreement with the work of Matsuoka et al. (2018). Furthermore, our best-fit QLF model has a steeper brightend slope than all previous measurements, reducing the integrated emissivity of the luminous quasar contribution as well.
In order to calculate the photoionization rate based on our QLF measurement we need to adopt a value for the mean free path of ionizing photons at z ≈ 6. Becker et al. (2021) recently measured the mean free path of ionizing photons and find a value of l(ν 912 , z = 6) = 0.75 +0.65 −0.45 pMpc at z = 6, which falls below extrapolations from lower redshift. We adopt this value, noting that these measurements are in agreement with the independently calculated lower limits reported by Bosman (2021).
With these assumptions we calculate a quasar photoionization rate of Γ HI (z=6) = 5.86 +5.08 −3.51 +1.33 −0.83 × 10 −16 s −1 . The first errors reflect the 1σ uncertainties of the mean free path (Becker et al. 2021), and the second errors are due to the statistical 1σ uncertainty in the QLF double power law fit. The total photoionization rate at z ≈ 6 has been measured to be Γ HI (z = 6) ≈ 10 −13 s −1 based on quasar near-zone sizes (Wyithe & Bolton 2011;Calverley et al. 2011) and the mean transmitted Lyα flux from quasar spectra (D'Aloisio et al. 2018;Davies et al. 2018). As a combination of our low emissivity values based on our new QLF measurement and the short mean free path of Becker et al. (2021), the quasar contribution to the photoionization rate is roughly two orders of magnitude lower than the total value. This result strongly disfavors a dominant contribution of quasars to cosmic hydrogen reionization at high redshifts, in line with other recent studies of the high redshift QLF Parsa et al. 2018;Matsuoka et al. 2018;McGreer et al. 2018;Kulkarni et al. 2019;Wang et al. 2019a;Jiang et al. 2022).

Quasar lensing and the QLF
None of the 125 quasars in our sample is known to be gravitationally lensed by a foreground galaxy. This is a direct consequence of our selection criteria (Section 2.2. Our morphology selection criterion aims at selection point sources and thus naturally excludes lensed quasars, which appear extended if they consist of multiple sources images or for which a foreground lens galaxy is detected. In addition, the required color criteria are designed to select quasars using the Lyman-α break. In the case that a foreground lens galaxy contaminates the bluer bands, the criteria bias against the selection of such sources. Following the serendipitous discovery of the highestredshift lensed quasar, J043947.08+163415.7, at z = 6.51 , Pacucci & Loeb (2019) reevaluated the theoretical consequences of this discovery and conclude that a large fraction of quasars at z > 6 are missed by current surveys. For a bright end slope of β = −3.6, the authors expect about half of the z > 6 population to be lensed. The steep bright-end slope of β = −3.84 of our new QLF measurement would lead to an even larger lensed fraction of the z > 6 quasar population. In consequence, our bright-end slope measurement would strongly suggest that we are missing a large fraction of lensed quasars. Given that only one lensed quasar at z > 6 has been discovered , our bright-end slope measurement and the resulting lensed fraction according to Pacucci & Loeb (2019) is in strong tension with observations. The recent study by Yue et al. (2022) revisits the predicted fraction of high-redshift lensed quasars. By adopting recent galaxy velocity dispersion functions that affect the lensing optical depth, they conclude that the lensed fraction for bright quasars at z ∼ 6 can reach 2% to 6% depending on the QLF bright-end slope. Following their Figure 5 and adopting M lim − M * ≈ 1 mag for our the PS1 quasar selection, results in a lensed fraction of 1% for β ≈ −4. In line with observations , this result suggests that even with our steep bright-end slope we would have found only ∼ 1 lensed quasar with our selection, hadn't it been biased against these sources.

CONCLUSIONS
In this work we present the most precise measurement of the z ≈ 6 QLF at −28 M 1450 −22 to date, based on the combined sample of 121 quasars from the Pan-STARRS 1 z > 5.6 quasar survey (PS1) and 48 quasars from SHELLQs. We determine the full PS1 quasar survey completeness taking into account the different components of the PS1 quasar selection strategy and the state of the spectroscopic observations. We use a maximum likelihood approach (see Appendix A) sampled via MCMC in order to fit a double power law QLF model to the quasar data. Our fiducial model (Table 4, first row) is determined on the combined quasar sample and assumes an exponential evolution of the quasar density with k = −0.7. The four best fit parameters are log(Φ * (z = 6)/Mpc −3 mag −1 ) = −8.75 +0.47 −0.41 , M * = −26.38 +0.79 −0.60 mag, α = −1.70 +0.29 −0.19 , and β = −3.84 +0.63 −1.21 . The combination of the PS1 and SHEL-LQs quasar sample constrains the break magnitude to be ∼ 1 mag brighter and the bright-end slope to be significantly steeper than previous studies at this redshift Jiang et al. 2016;Matsuoka et al. 2018).
Using our fiducial QLF model we calculate the brightend quasar density, n(M 1450 < −26, z = 5.85) = 1.16 +0.13 −0.12 cGpc −3 , and put it in perspective with its redshift evolution at z ≈ 4 − 7. We find that an exponential density evolution model with an exponent of k = −0.7, as assumed in our QLF fit, describes the literature data over this redshift range well without the need for an accelerating decline of the quasar density at z > 6.5 as proposed by .
With our fiducial QLF model we derive the ionizing emissivity of the quasar population and their contribution to cosmic hydrogen reionization. Using standard assumption, we calculate the ionizing emissivity to be 912 (z = 6) = 7.23 +1.65 −1.02 × 10 22 ergs −1 Hz −1 cMpc −3 . This result is lower than some previous results (e.g., Jiang et al. 2016), but shows good agreement with Matsuoka et al. (2018), the most recent estimate of the z ≈ 6 QLF. Adopting the mean free path of Becker et al. (2021), the only measurement at z ≈ 6, we estimate a HI quasar photoionization rate two order of magnitudes below estimates of its total value, strongly disfavoring quasars as a dominant driver of hydrogen reionization at z ≈ 6.
The authors would like to thank the anonymous referee for their insightful comments that helped to improve the manuscript. We would like to further thank Sarah E. I. Bosman, Joseph F. Hennawi, and Romain Meyer for fruitful discussions and comments on the manuscript. JTS and RN acknowledge funding from the European Research Council (ERC) Advanced Grant program under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 885301  To derive the QLF it is imperative to take into account selection effects that are inherent in the parent catalog data and are imposed by the selection criteria for the quasar discovery survey (see Section 2.2). The selection function S(q) describes the probability of a source with attributes q to be within the given sample. In this section we explicitly derive the mathematical formulation of our MCMC maximum likelihood analysis of the QLF.
A.1. The quasar incidence as predicted by the QLF We begin our discussion by closely following Rix et al. (2021, their Equation 1) to describe the expected catalog incidence dΛ(q) of quasars in our sample through the selection function S(q) multiplied with a model family for quasars M(q|Θ mod ), parameterized by Θ mod : In our case the QLF Φ forms the basis of the model family. In its most general form the QLF can be written as a function of the luminosity (in our case the absolute magnitude at 1450Å), the three dimensional position of quasars x = (r, θ, φ) (in spherical coordinates) and the QLF parameters Θ QLF . The distance to quasars r is not a direct observable. Therefore, it is much more practical (and common) to formulate the quasar luminosity function as a function of redshift z. Furthermore, we have good reason to assume that our universe is isotropic on large scales. In this case the QLF is independent of sky position (θ, φ). However, the selection function may depend on the sky position and thus we separate the sky dependence from the QLF using the unit normal vectorn(θ, φ).
M(q|Θ mod ) = Φ(M 1450 , z, θ, φ|Θ QLF ) = Φ(M 1450 , z|Θ QLF )n(θ, φ) For clarity of the mathematical expressions we omit the subscript to the absolute magnitude at 1450Å, M 1450 , in the following and simply denote it with M . We now substitute Equation A2 into Equation A1 and integrate both sides over volume (dV ) and absolute magnitude (dM ) to retain the total expected number of quasars as observed given the model and the selection function: Due to the redshift dependency of the QLF executing the volume integral requires a cosmological model with its own range of parameters Θ Cos . We will now rewrite the volume integral in terms of the differential comoving solid volume element (dV / dz/ dΩ), a standard quantity in any cosmological model: in Equation A3 allows us to separate the volume integration into integrals over redshift and solid angle: For surveys of inhomogeneous depth it can be often difficult to find an analytic expression for the sky position dependence of the selection function. In this work we take the inhomogeneity into account when modeling the observed quasar properties by sampling from the depth distribution. Therefore, we continue as one would with a survey of homogeneous depth and drop the sky position dependence. Now the integral over they survey footprint simply yields the total footprint area Ω.
Given a model for the quasar luminosity function Φ(M, z|Θ QLF ), a cosmological model Θ Cos and a model of the observed properties given absolute magnitude and redshift q(M, z), we can now calculate the total expected number of observed quasars as a function of absolute magnitude and redshift.

A.2. Formulating the likelihood function
We derive the likelihood function for the QLF analysis following Marshall et al. (1983) and Fan et al. (2001b). The probability of detecting n lm quasars given the QLF Φ(M, z|Θ QLF ) in an absolute magnitude bin (∆M ) l and redshift bin (∆z) m can be written in terms of the Poisson distribution function using the incidence rate Λ lm (M, z) of Equation A6: The probability of finding N (M, z) quasars in the entire survey, as characterized by the full absolute magnitude and redshift range, can then be written as the product of the probabilities over all absolute magnitude and redshift bins.
If the absolute magnitude and redshift bins are infitesimally small, then either n lm = 1 or n lm = 0 quasars can be found in each bin. We split the product into two terms for these two cases, simplifying the equation. We can furthermore, rearrange the terms to arrive at the final version of the probability P ( We have already discussed how to express the second term in this equation via the incidence rate in Equation A8. We now formulate the logarithmic probability equivalent to the logarithmic likelihood: To evaluate this equation further we will take a look at the incidence rate Λ lm,j (M, z) for a single quasar j. The quasar j has an absolute magnitude M j and redshift z j in the bin centers of width (∆M ) l and (∆z) m . With these boundary conditions we can write Equation A6 for a single quasar as In the limit of infinitesimal bin sizes around M j and z j the integrals can be trivially evaluated and we can drop the indices l and m: Starting from the first term of the right hand side of Equation A11 we first rewrite the sum over all bins l and m for which n lm = 1 as a sum over all N (M, z) quasars in the data set. Then we apply the natural logarithm to Equation A13 and equivalent to Marshall et al. (1983) drop all terms independent of the QLF and the selection function, which would only add constant values to the logarithmic probability.
The second term in Equation A11 normalizes the logarithmic probability by summing over the full interval in absolute magnitude ∆M and redshift ∆z for which we aim to evaluate the QLF. By choosing the appropriate integration boundaries for the QLF evaluation in Equation A6, we can simplify the expression to To test the robustness of our simqso quasar model we compare its synthetic photometry to the recently published qsogen model by Temple et al. (2021). The latter has been used in the prediction of quasar yields for the Euclid mission (Euclid Collaboration et al. 2019). At the core of qsogen is a parametric quasar spectral model, which has been fit to optical to near-infrared colors of a subsample of SDSS DR16 quasars (Lyke et al. 2020). Equivalent to the simqso quasar model, the qsogen quasar models is explicitly designed to reproduce the median quasar colors of real sources. A particular novel feature of qsogen are the emission line quasar templates, that allow to capture the diversity in quasar line strength and velocity shifts and can be modified via the emline type argument, which has a default value of −0.9936.
The quasar selection criteria on which our new QLF measurements is based rely on color cuts in the reddest PS1 filter bands, i P1 , z P1 , y P1 . We use qsogen to calculate the magnitudes in these bands and then compare the resulting i P1 − z P1 , and z P1 − y P1 colors between the qsogen and our quasar model based on simqso. To capture the quasar diversity in terms of emission line properties with qsogen, we produce three models: the default model uses the default value of emline type. The strong emission line model uses emline type= 1 and the weak emission line model uses emline type= −1. We compare the resulting colors as a function of redshift in Figure 16. The z P1 − y P1 color of our quasar model and the qsogen models is consistent over z = 5.2 − 5.8. The diversity in quasar colors due to the different qsogen emission line models is well captured by the random sampling of continuum and emission line properties used in simqso (Section 3.1.1). Small differences the z P1 − y P1 color appear at z > 5.8, which increase towards z = 6.6. We attribute the difference to the different prescriptions of absorption of neutral hydrogen due to the IGM, which become more pronounced at z > 5.8 as the Lyman-α line moves to the red side of the z P1 filter band. The selection of z > 5.7 quasars is based on a color dropout selection requiring i P1 − z P1 > 2. The right panel of Figure 16 shows i P1 − z P1 as a function of redshift and highlights the dropout criterion with the black solid line. Quasars with colors above the black line would be selected. With our simqso quasar model we start to select quasars at z > 5.6 (blue solid line and blue shaded region). The downturn of the i P1 − z P1 color at z 6.4 is due to the Lyman-α line moving out of the z P1 filter band. As a consequence little to no quasar flux is measured in both the i P1 -and z P1 band. Our simqso 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6  Figure 10). At z 6.3 the color is biased blue by the non-evolving IGM transmission in the templates. Right panel: Synthetic quasar iP1 − zP1 color as a function of redshift. The symbols are the same with the addition of the black solid line that depicts the iP1 − zP1 > 2 color criterion used in this work.

Redshift
quasar model mostly agrees with the qsogen models at z 5.6. However, once the Lyman-α line leaves the i P1 band (z 5.6) color differences become evident. We argue that these are mainly caused by the different prescriptions of IGM absorption in the spectral models. While simqso uses the stochastic Lyman-α forest model of McGreer et al. (2013) (see Section 3.1.1), qsogen uses a prescription based on Becker & Bolton (2013). Compared to our simqso model, the qsogen spectra have more flux blueward of Lyman-α, resulting in a generally bluer color. When compared to quasars selected in this work (grey points) the simqso IGM model provides a better representation of the empirical data points 6 . The qsogen models show significantly bluer i P1 − z P1 , especially at z > 5.6. Introducing a Lyman limit system at the redshift of the quasar (i.e. at rest-frame wavelength of 1215Å) to the qsogen default model, effectively reduces all blueward flux to zero. At z 6.0 this approach is unphysical as the flux blueward of Lyman-α can be transmitted in ionized patches of the IGM. However, it is valid assumption at higher redshifts (see e.g. at z > 7, Euclid Collaboration et al. 2019). We show this qsogen Lyman-limit model (LL) as the green line in Figure 16. The qsogen LL model has a significantly redder i P1 − z P1 color at z < 5.7, which diverges further at z ≈ 5.7. Compared to the PS1 quasar photometry at z ≈ 5.6 − 6.4 it is not a good representation of the data. However, the resulting z P1 − y P1 color of our quasar model and qsogen is now fully consistent over the entire redshift range. This agreement in the z P1 − y P1 color highlights that, barring the differences in the prescription of IGM absorption, both our simqso quasar model and the qsogen model both produce consistent median colors as expected from their design goals. . In these three panels we show the observed survey samples, their selection efficiency and identification completeness for different visual ranks. Comparing the samples highlights that the majority of remaining candidates have the worst visual rank (vis sum= 4) and the lowest selection efficiency. Our identification completeness (Equation 15) takes these differences into account.

C. SUPPLEMENTAL FIGURES
We provide a number of figures that provide supplemental information for a few sections of the paper. Figure 17 shows the spectroscopic identification completeness for quasar samples with visual ranks 2,3, and 4 separately. For context of the covariance of the different QLF parameters in our ML fit we provide the full covariance matrix in Figure 18. Analogous to Figure 3, we also show the PS1 y-band limiting magnitude as a function of survey area in Figure 19.