The Infrared Medium-deep Survey. VIII. Quasar Luminosity Function at $z\sim5$

Faint $z\sim5$ quasars with $M_{1450}\sim-23$ mag are known to be the potentially important contributors to the ultraviolet ionizing background in the post-reionization era. However, their number density has not been well determined, making it difficult to assess their role in the early ionization of the intergalactic medium (IGM). In this work, we present the updated results of our $z\sim5$ quasar survey using the Infrared Medium-deep Survey (IMS), a near-infrared imaging survey covering an area of 85 deg$^{2}$. From our spectroscopic observations with the Gemini Multi-Object Spectrograph (GMOS) on the Gemini-South 8 m Telescope, we discovered eight new quasars at $z\sim5$ with $-26.1\leq M_{1450} \leq -23.3$. Combining our IMS faint quasars ($M_{1450}>-27$ mag) with the brighter Sloan Digital Sky Survey (SDSS) quasars ($M_{1450}<-27$ mag), we derive the $z\sim5$ quasar luminosity function (QLF) without any fixed parameters down to the magnitude limit of $M_{1450}=-23$ mag. We find that the faint-end slope of the QLF is very flat ($\alpha=-1.2^{+1.4}_{-0.6}$), with a characteristic luminosity of $M^{*}_{1450}=-25.8^{+1.4}_{-1.1}$ mag. The number density of $z\sim5$ quasars from the QLF gives an ionizing emissivity at 912 $\unicode{x212B}$ of $\epsilon_{912}=(3.7$--$7.1)\times10^{23}$ erg s$^{-1}$ Hz$^{-1}$ Mpc$^{-3}$ and an ionizing photon density of $\dot{n}_{\rm ion}=(3.0$--$5.7)\times10^{49}$ Mpc$^{-3}$ s$^{-1}$. These results imply that quasars are responsible for only 10-20% (up to 50% even in the extreme case) of the photons required to completely ionize the IGM at $z\sim5$, disfavoring the idea that quasars alone could have ionized the IGM at $z\sim5$.


INTRODUCTION
Quasars are known as key objects for understanding the universe along cosmic time, especially at high redshifts (z 5) where galaxies are hard to detect. The * KIAA Fellow Sloan Digital Sky Survey (SDSS) pioneered the identification of quasars at z ∼ 6, providing the first glimpse of the reionization process in the early universe (e.g., Fan et al. 2001aFan et al. , 2006Jiang et al. 2016). As new surveys have explored wider areas and deeper limits, the number of known quasars at high redshifts has steadily increased in the last few decades. Along with the record holder ULAS J1345+0928 at an extreme redshift of z = 7.54 , hundreds of high-redshift quasars are being discovered by various surveys (Willott et al. arXiv:2010.05859v1 [astro-ph.GA] 12 Oct 2020 2010; Mortlock et al. 2011;McGreer et al. 2013McGreer et al. , 2018Venemans et al. 2013Venemans et al. , 2015aBañados et al. 2014Bañados et al. , 2016Kashikawa et al. 2015;Kim et al. 2015Kim et al. , 2019Wu et al. 2015;Jiang et al. 2016;Matsuoka et al. 2016Matsuoka et al. , 2018aMatsuoka et al. , 2019Wang et al. 2016Wang et al. , 2017Wang et al. , 2018Wang et al. , 2019Yang et al. 2016Yang et al. , 2017Yang et al. , 2019aYang et al. ,b, 2020Jeon et al. 2017;Reed et al. 2017;Shin et al. 2020). This large sample of distant quasars enables us to construct quasar luminosity functions (QLFs) at high redshifts, which are important for investigating the contribution of high-redshift quasars to the cosmic ultraviolet (UV) ionizing background and the evolution of quasar populations along the redshift. Since the cosmic reionization at 6 < z < 8.8 (e.g., McGreer et al. 2015;Planck Collaboration et al. 2016;Madau 2017;Davies et al. 2018;Mason et al. 2018Mason et al. , 2019, most of the hydrogen atoms in the intergalactic medium (IGM) have remained ionized until z = 0. It has been debated whether high-redshift quasars that are bright in UV wavelengths can provide the enormous quantity of ionizing photons required to keep hydrogen atoms ionized. The number of ionizing photons from quasars can be calculated from their UV emissivity, proportional to both luminosity and quasar number density. According to recently derived QLFs, the emissivity critically depends on the faint quasars at M 1450 ∼ −23.5 mag (absolute magnitude at 1450Å in the rest frame). Figure 1 shows the UV emissivity from the quasars for different magnitude bins. There are differences between QLFs, but M 1450 ∼ −23.5 mag quasars make a considerable contribution to the total quasar emissivity. However, it has been difficult to locate them in the past owing to observational limitations, such as the shallow imaging depths of large-area surveys or the small survey areas of deep surveys. As a result, past determination of the QLF at M 1450 ∼ −23.5 mag often relied on uncertain extrapolations of the QLF (Giallongo et al. 2015;Yang et al. 2016).
Owing to new deep and wide-area optical/nearinfrared (NIR) surveys, tens of such faint quasars are now being uncovered (e.g., Kashikawa et al. 2015;Kim et al. 2015Kim et al. , 2019Matsuoka et al. 2016Matsuoka et al. , 2018aMatsuoka et al. , 2019Akiyama et al. 2018;McGreer et al. 2018;Shin et al. 2020). However, the contribution of quasars to the IGMionizing UV photons remains controversial. On one hand, the results of optical/NIR wide-area surveys suggest that high-redshift quasars are not the main provider of UV photons, contributing < 50% of the required photons at z ∼ 5 (Yang et al. 2016;McGreer et al. 2018;Shin et al. 2020) and < 10% at z ∼ 6 (Willott et al. 2010;Kashikawa et al. 2015;Kim et al. 2015;Matsuoka et al. 2018b). On the other hand, from very deep surveys covering small areas ( 4 deg 2 ), other studies have re- . Differential contributions to z ∼ 5 quasar emissivity, calculated from the QLFs in the literature: Giallongo et al. (2015;red), Yang et al. (2016;orange), McGreer et al. (2018;teal), and Giallongo et al. (2019;purple). The gray shaded region shows the magnitude range in which quasars make a considerable contribution to total quasar emissivity.
vealed a large number of faint quasars at high redshifts, supporting a significant contribution of quasars to the IGM ionizing background (Glikman et al. 2011;Boutsia et al. 2018;Grazian et al. 2020). In particular, for Xray-selected quasars with M 1450 ≥ −22 mag, Giallongo et al. (2019; hereafter G19) have suggested a substantial contribution of quasars ( 50%, depending on z) to the IGM ionizing background (see also Giallongo et al. 2015). At z ∼ 5, McGreer et al. (2018, hereafter M18), performed a quasar search using the data from the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS; Hudelot et al. 2012). They derived the QLF at 4.7 < z < 5.4, combining the bright SDSS quasar sample of McGreer et al. (2013) and 25 newly discovered CFHTLS quasars with M 1450 −23 mag. Their QLF, in the form of a double power-law function, has a faint-end slope of −1.97 and a low number density, implying that quasars are unlikely to be the main UV ionizing source at z ∼ 5. Due to the lack of faint quasars identified with spectroscopy, however, they fixed the bright-end slope of the QLF to −4.0, which may have given biased interpretations of the quasar population (and corresponding emissivity) with underestimated uncertainties (Kulkarni et al. 2019). To date, only Yang et al. (2016, hereafter Y16), have derived the z ∼ 5 QLF without fixed parameters, including the quasar sample of McGreer et al. (2013). However, most of the quasars they used are M 1450 < −24 mag, which might not be sufficient to de-rive the QLF precisely. In fact, their best-fit parameters show discontinuity along the redshift, compared to the recently reported QLFs at z ∼ 4 (Akiyama et al. 2018) and 6 (Matsuoka et al. 2018b), which are from a wide range of luminosities (−30 M 1450 −22) over a large survey area without any fixed parameters.
Recently, our group has been performing a z ∼ 5 quasar survey (Kim et al. 2019, hereafter K19), with the Infrared Medium-deep Survey (IMS; M. Im et al. 2020, in preparation), for which NIR imaging data were obtained with the Wide Field Camera (WFCam; Casali et al. 2007) on the United Kingdom Infrared Telescope (UKIRT). The data cover an area of ∼ 100 deg 2 over several extragalactic fields and reach 5σ depths of J ∼ 23 mag. We combined the IMS data with the optical data from the CFHTLS (Hudelot et al. 2012) and additional NIR data from the Deep eXtragalactic Survey (DXS) of the UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007). Adopting the medium-band color selections that are known to improve the high-redshift quasar selection efficiency , we spectroscopically identified 13 faint quasars at z ∼ 5 in the past, of which 10 were newly discovered. Increasing the number of spectroscopically identified quasars enabled us to refine and constrain the faint-end slope of the z ∼ 5 QLF. In this study, we have added more samples and derived the QLF at z ∼ 5.
The remainder of this paper is organized as follows. We present the newly discovered z ∼ 5 quasars from our survey in Section 2, including a summary of our survey progress. In Section 3, we describe the derivation of the QLF using the IMS z ∼ 5 quasars. We discuss the implications of the z ∼ 5 QLF regarding the cosmic ionizing background in Section 4. We adopt the cosmological parameters of Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 , which are supported by observations in recent decades (e.g., Im et al. 1997;Planck Collaboration et al. 2018). All magnitudes given in this paper are in the AB system.

NEW IMS QUASARS AT z ∼ 5
Here we summarize the IMS z ∼ 5 quasar survey, described in detail in K19, with the newly discovered quasars at z ∼ 5.

Initial Broadband Selection
The initial photometric dataset is a combination of the optical imaging data from the CFHTLS Wide Survey and the NIR imaging data from IMS and DXS, covering four extragalactic fields: XMM-Large Scale Structure survey region (XMM-LSS), CFHTLS Wide survey second region (CFHTLS-W2), Extended Groth . Medium-band color-color diagram of the quasars observed in all three medium bands (m675/m725/m775). The circles denote 37 quasars that are spectroscopically identified at 4.7 < z < 5.4, while the 32 quasars satisfying our color selection criteria (dotted lines) are highlighted as filled circles. The arrows indicate the upper limits on the magnitude according to the detection limits of their medium-band images. For the five quasars outside the selection boxes, their (z, M1450) values are denoted at the lower left of each symbol. The newly spectroscopically confirmed quasars by us (K19 and this work) are highlighted with squares. The gray diamonds with lines denote the quasar color tracks from z = 4.5 to 5.3.

Strip (EGS)
, and Small Selected Area 22h (SA22). The average 5σ limiting magnitudes for point sources are u = 26.1, g = 26.4, r = 25.9, i = 25.6, z = 24.6, and J = 22.9 mag. Note that the optical photometric system here is based on the SDSS filter system transformed from those of the CFHTLS 1 . The total survey area used for the z ∼ 5 quasar search covers ∼ 85 deg 2 of the sky, calculated from the full overlapping region of the surveys. For the i -band detected sources, we selected z ∼ 5 quasar candidates using the broadband color selection criteria: (1) i < 23, (2) S/N (u ) < 2.5, In the visual inspection stage, seven sources that are unusually elongated or ex-tended were rejected, giving the final 69 z ∼ 5 quasar candidates 2 .

SQUEAN Follow-up Imaging in Medium Bands
To narrow down the number of plausible z ∼ 5 quasar candidates, follow-up observations of these broadbandselected candidates were carried out in medium bands with the SED Camera for Quasars in Early Universe (SQUEAN; Choi et al. 2015;Kim et al. 2016) on the Otto Struve 2.1 m Telescope at McDonald Observatory. At z ∼ 5, the Lyα λ1216 line is located at ∼ 7300Å, so we have obtained m675, m725, and m775 images, where the filter names are chosen by combining "m" for medium-band and the number of the central wavelength of the filter in nanometers. These 50 nm width filters finely sample the redshifted Lyα line to improve quasar identification and determine their redshifts to nearly 1% accuracy ; K19). The total exposure times in medium bands are proportional to the i -band magnitude of each target, as are the resultant imaging depths. These inhomogeneous imaging depths between targets are important when we calculate the medium-band selection functions in Section 3.2. Most of the broadband-selected candidates were observed in the m725 and m775 bands (63/69), while 50 were also observed in the m675 band. We introduced the medium-band color selection criteria of Jeon et al. (2016): (1) m675 − m725 > 1.0 and m675 − m725 > m725 − m775 + 1.5 (4.7 < z < 5.1), (2) m725 − m775 > 1 (5.1 < z < 5.5). Among the 50 medium-band observed candidates, 33 satisfy the criteria. In Figure 2, we plot a m675−m725 vs m725−m775 diagram with the medium-band color criteria shown as dotted lines. For simplicity, we only present the 37 spectroscopically identified 4.7 < z < 5.4 quasars that are also observed in the three medium bands. There are five quasars excluded by our selection (open circles), possibly due to their m675 imaging depths being shallower than expected for a given i -band magnitude as described in Section 3.2, and/or their marginal redshift of z 4.7.

GMOS Spectroscopy
We obtained the optical spectra of the plausible candidates, which satisfied either the above broadband or medium-band color selection criteria, with the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004) on the Gemini-South 8 m Telescope on 2018 Sep 17 2 The total number of broadband-selected candidates is decreased from 70 in K19 because of the miscoded procedure. Four candidates, including two spectroscopically identified non-quasars, are rejected, and three candidates without medium-band observations are included instead. GMOS optical spectra of the newly discovered IMS z ∼ 5 quasars. The black solid lines represent the binned spectra according to spectral resolution, while the red solid lines are the best-fit models. The blue marks are the wavelengths of possible quasar emission lines: Lyβ, Lyα, N V, O I, Si IV, and C IV, from short to long wavelengths. The shaded regions represent the bad columns, hot pixels, CCD gaps, or wavelength range that are not covered by the observational configuration.
(PID: GS-2018B-Q-217), 2019 Jan 2-11, and 2019 Feb 14-24 (PID: GS-2019A-Q-218), under the seeing condition of 1. 0. The observing configurations were set to increase the signal-to-noise ratio (S/N) of the spectrum; we used the Nod & Shuffle mode for sky subtraction, an R150 G5326 grating with a resolution of R ∼ 315, and spectral/spatial binning of 4 × 4. In total, we obtained the spectra of ten candidates observed in at least two medium bands. For spectral data reduction, we followed the procedure in K19 which is summarized here briefly. We used the Gemini IRAF package for the basic reductions: the bias subtraction, flat-fielding, sky-lines subtraction, wavelength calibration with CuAr arc lines, and flux calibration with standard stars. After the 1D extraction with an aperture with a diameter of 1. 0, the overall flux of each spectrum was scaled with the iband magnitude of each target. To maximize S/N, we binned the spectra along the wavelength accounting of the instrumental resolution of ∼ 300, using the inversevariance weighting method (e.g., Kim et al. 2018).
Through these spectroscopic observations, eight z ∼ 5 quasars were newly discovered, the optical spectra of which are shown in Figure 3. All of them show strong Lyα emission lines with sharp breaks, and some also show other possible emission lines like C IV, while there are broad absorption lines (BAL) in the spectrum of IMS J090508−021523.
We list the broadband and medium-band photometry of the new IMS z ∼ 5 quasars in Table 1, and their spectroscopic observations are summarized in Table 2. Note-All magnitudes are given in the AB system. The detection limits (∼ 3σ) are given as the upper limits for objects that were undetected or fainter than the limit.  The other two candidates identified as non-quasars are listed in Appendix A.

Spectral Fitting
We measured the spectral properties of the eight new quasars by fitting their spectra with a refined version of our high-redshift quasar model (K19). The model is based on the composite quasar spectrum of SDSS quasars Vanden Berk et al. (2001). Since no IGM correction has been applied in this spectrum, we replaced the spectrum at λ < 1216Å to that of Lusso et al. (2015).
Note that the two spectra were normalized at 1450Å before merging. For this composite spectrum, the IGM attenuation by neutral hydrogen was applied in the form of a function of redshift (Madau et al. 1996). The model has four parameters: redshift (z), absolute magnitude (M 1450 ), continuum slope (α λ ), and equivalent width (EW) of Lyα λ1216 and N V λ1250 (EW Lyα+NV ). The SDSS composite spectrum is first decomposed into the emission line component and the continuum component by fitting a power-law continuum to it. The slope of the continuum component is allowed to vary to a given α λ by multiplying a factor of (λ rest /1000Å) α λ +1.54 . EW Lyα+NV , estimated from the composite continuumsubtracted flux at 1160 ≤ λ rest (Å) ≤ 1290, is also scaled with (λ rest /1290Å) p , where p is the appropriate power to adjust the EV Lyα+NV value of the model to an arbitrary one. In the original K19 model, only EW of Lyα+N V was allowed to scale. On the other hand, in the refined K19 model, we also adjusted the other emission lines, such as C IV λ1549 and C III] λ1909, by scaling fluxes over continuum with EW Lyα+NV . This is because EWs of these emission lines are known to scale with EW Lyα+NV (Dietrich et al. 2002) and have EW distributions with scatters similar to that of Lyα+N V (0.2-0.3 dex; Diamond-Stanic et al. 2009;Shen et al. 2011).
We found the best-fit model for each quasar spectrum by finding the minimum χ 2 values. Considering the narrow wavelength coverage of our spectra, we fixed α λ to −1.54, same as K19. In Figure 3, the best-fit models for the eight new quasars are shown as the red solid lines, and we listed the best-fit values in Table 3. For IMS J090508−021523, a BAL quasar, the fluxes at BAL wavelengths are excluded from the fitting. Note that this model is also used to determine the selection functions in Section 3.2.
Including our eight newly discovered quasars, 49 z ∼ 5 quasars with i < 23 mag have been spectroscopically identified in our survey area (McGreer et al. 2013;Wang et al. 2016;M18;K19). With these quasars, we constructed the QLF at z ∼ 5. However, for the QLF construction, we excluded four quasars that were not selected by broadband color (see K19) and two quasars outside the redshift range of 4.7 < z < 5.4. Figure 4 shows the sky distributions of the remaining 43 quasars, color-coded by their magnitudes. We used the z and To examine the effect of adding the medium-band color selection, we considered either the broadband selection only (case 1) or a combination of broadband and medium-band selection (case 2). The numbers for each case are 43 and 32 in cases 1 and 2, respectively, and their z and M 1450 values are given in Table 4. The five quasars observed in the three medium bands but outside the selection boxes (the open circles in Figure 2) could have been selected if we obtained deeper images. However, we excluded them from the case 2 sample and considered this issue in Section 3.2

Quasar Selection Function
Using our high-redshift quasar model described in Section 2.4, we calculated the selection efficiency of our color selection criteria. We generated mock spectra of 100,000 quasars randomly distributed in the range of 4.4 < z < 6.0 and −28 < M 1450 < −22, while α λ and EW Lyα+NV were randomly generated by Gaussian distributions with mean and standard deviation values of α λ = −1.6 ± 1.0 (Mazzucchelli et al. 2017) and log EW Lyα+NV = 1.803 ± 0.205 (Diamond-Stanic et al. 2009), respectively, including the Baldwin effect of Lyα (Dietrich et al. 2002). We integrated the mock quasar spectra over the broadband and medium-band transmission curves to obtain the magnitudes used for color selection. In Figure 5, we show the color distributions of these mock quasars. They are in line with not only the IMS quasars but also the known z ∼ 5.5 quasars in the same photometric system (Yang et al. , 2019a. Note that the relatively larger scatter in r − i color appears to be due to the various proximity zone sizes of high-redshift quasars according to their ages (e.g. Eilers et al. 2017).
The imaging survey depths are not always homogeneous, and the varying imaging depths at different locations can affect the selection function. For case 1, we generated broadband depth maps resampled at a pixel scale of 1. 0. For each pixel of the broadband depth maps and z-M 1450 bin with a size of ∆z = 0.05 and ∆M 1450 = 0.1 mag, we calculated the selection completeness, a ratio of the number of quasars satisfying our broadband color selection criteria to the total number of mock quasars in the bin (the average number of mock quasars in each bin is ∼ 50). In this process, we applied additional Gaussian noise to the model magnitudes, considering the image depths and model magnitudes. In Figure 6, the calculated selection functions are shown as gray contours on the left panel. Although our J-band data were inhomogeneous (K19), the differences between the selection functions of the four extragalactic fields of our survey were negligible. We used the combined selection function weighted by the area of each survey field.
For case 2, we calculated a selection function including our medium-band color selection criteria. The imaging depths of the follow-up medium-band observations are not uniform but depend on the i -band magnitudes of the targets, because we gave exposures proportional to the i -band magnitude as much as we could, although the detection limits vary depending on the observing conditions at the time of observation and whether the object was easily detected or not. If the object was detected with a sufficient S/N ( 10), we terminated observation of the target in the corresponding band to save the observing time. Figure 7 shows the detection limits (2.7σ) for a point source in the m675, m725, and m775 images of 32 quasars of the case 2 sample, as a function of their i -band magnitude. We fitted a linear regression model to these imaging depths, shown as the solid lines with confidence intervals calculated with the seaborn Python package 3 . Out of 32 quasars, 22 are not detected in m675 (open circles), implying small fluxes below the redshifted Lyα. However, we used all of them for linear regression fitting, because we consider such undetected cases in determining the selection function of case 2. As shown in Figure 2, for objects without detection, we used their detection limits as their magnitudes for the color selection. In the case of mock quasars, similarly, medium-band magnitudes fainter than the detection limit expected from their i -band magnitudes (i.e., linear regression models) are replaced with the expected values. The color selection with these modified colors agreed with our medium-band color selection process, which could miss some quasars in marginal conditions (e.g., open circles in Figure 2). The right panel of Figure 6 shows the selection function for case 2. Note that Note-The "Case" columns indicate whether the quasars are included in the case 1 and 2 samples (see details in Section 3.1). The spectral properties are from (1)  a These M1450 values are not taken from their spectra but from the fitting with their photometric data, as described in K19.
b This quasar is excluded in our parametric QLF derivation owing to M1450 < −27 mag.
c These values are revised in K19 from those published in M18. there is a valley at z ∼ 5.1, which corresponds to a specific redshift for which the medium-band selection fails (see Figure 2). It has been suggested that the EW Lyα+NV values of z > 5.7 quasars are lower (e.g., log EW Lyα+NV = 1.542 ± 0.391; Bañados et al. 2016) than their lowredshift counterparts, which might be due to the decrease in Lyα transmission at λ > 1216Å by neutral hydrogen (e.g., Davies et al. 2018). At z ∼ 5, the neutral hydrogen fraction is much lower than 0.1 (Fan et al. 2006;McGreer et al. 2015), and thus it has a negligible effect on the redward Lyα fluxes. If we introduce the EW distribution of Bañados et al. (2016), the selection efficiencies of cases 1 and 2 are reduced by about 30% and 50% overall, respectively.

IMS Survey Completeness
We first consider the completeness for case 1. In the process of selecting the 69 broadband-selected candidates, we used the magnitude cut of i < 23 mag, at which the CFHTLS i -band (or i CFHTLS ) images we used for the source detection were almost complete (Hudelot et al. 2012  Medium-band imaging depths (2.7σ) of 32 quasars of case 2 sample along the i-band magnitude. The blue, orange, and red colors denote the cases of m675, m725, and m775, respectively. The filled/open circles represent whether a quasar is detected/undetected, respectively, while the solid lines with shaded regions are the linear regression models with confidence intervals. 50% and 80% completeness limits for a point source of each i CFHTLS -band image. We fitted these limits with an analytic completeness function (Fleming et al. 1995), resulting in 99% completeness at i CFHTLS < 23 mag for all CFHTLS images used. Even if we consider the transformation from i CFHTLS to i (i.e., from CFHTLS to SDSS), the point-source detection at i < 23 mag is almost complete. Therefore, for the broadband photometric completeness, we assume that our detection is complete. The point-source separation process was not included in our selection 4 , but note that we performed the visual inspection as mentioned in Section 2.1. The left panel of Figure 8 shows the number and completeness as a function of the i -band magnitude for case 1. Of the 69 candidates, 51 were spectroscopically observed and there are 43 quasars at 4.7 < z < 5.4 as mentioned above. We used the ratio of the candidates with spectroscopy to the broadband-selected candidates as the spectroscopic completeness.
For case 2, the completeness for both broadband and medium-band photometry should be considered. As with case 1, we also assume complete broadband photometry. In addition, we adopted medium-band photometric completeness. According to our criteria, quasars at 5.1 < z < 5.4 could be selected with the m725 and m775 bands alone, while 63 of the 69 broadband-selected 4 M18 adopted i AUTO − i PSF > −0.15 that showed high completeness of ∼ 98% for i 23 mag sources, where i AUTO is the automatic aperture magnitude (MAG AUTO) from SExtractor (Bertin & Arnouts 1996) and i PSF is the magnitude based on the point spread function (PSF). From our broadband-selected sample, only three candidates did not satisfy the criteria due to adjacent object(s) or subtle elongation. Two of them were identified as high-redshift quasars (IMS J022113−034252 at z = 5.02 and IMS J085028−050607 at z = 5.36), while the other was identified as a non-quasar object (IMS J090540−011038; see Appendix A). The exclusion of these targets does not affect the QLF determination. candidates (91%) were observed in these two bands. However, the selection of quasars at 4.7 < z < 5.1 requires the m675, m725, and m775 bands, in which 50 of 69 candidates (72%) were observed. To avoid complex and ambiguous calculations, the candidates with observations in the three medium bands were considered as the main sample in case 2. Therefore, we used the ratio of the number of candidates with observations in the three medium bands to the total number of broadband-selected candidates as the medium-band photometry completeness in case 2 (green histogram in the bottom right panel of Figure 8). As described in Section 2, 33 of the 50 candidates observed in the three medium bands satisfy our medium-band color selection criteria (yellow histogram in Figure 8). We also obtained spectra for 32 of the candidates, all of which were identified as 4.7 < z < 5.4 quasars. The red histogram in the bottom right panel of Figure 8 shows the spectroscopic completeness of the case 2 sample. To apply these photometric/spectroscopic completenesses to each case, we generated the photometric/spectroscopic completeness functions of (z, M 1450 ) using our high-redshift quasar model described in Section 3.2. For the mock quasars in each bin of the selection functions, we took the number-weighted mean values of their photometric/spectroscopic completenesses. We first calculated the binned QLF at z ∼ 5 using the 1/V a method from Avni, & Bahcall (1980). Given a bin with sizes of ∆M 1450 and ∆z, a specific co-moving volume of our survey V a can be calculated as Binned and parametric QLFs of quasars at 4.7 < z < 5.4 within the IMS coverage. The red and blue points are Φ bin for the IMS z ∼ 5 quasar samples in cases 1 and 2, respectively. The orange points are that of the Y16 sample, which are re-binned for bright quasars (M1450 < −27 mag). The red and blue solid lines with shaded regions show our best-fit results of Φpar with 1σ confidence level in cases 1 and 2, respectively. The dotted and dashed lines show the Φpar with fixed slopes of β = −3.0 and −4.0, respectively, in each case. (1) where p(M 1450 , z) is the total selection probability combining the quasar selection functions (Section 3.2) and the photometric/spectroscopic completeness functions (Section 3.3), and dV /dz is the co-moving volume element of our survey area. The binned QLF Φ bin (M 1450 ) can then be written as

Binned and Parametric QLFs
where N Q is the number of quasars in the bin. The bin size was set to ∆M 1450 = 0.5 mag (see Figure 6). After experimenting with several different bin sizes, we settled on this size, which gave a good compromise between the number of available quasars per bin and the fine sampling of the QLF. Note that the uncertainties of Φ bin are calculated using the root-sum-square method. We have listed the derived Φ bin for cases 1 and 2 in Table 5, which are also denoted by red and blue points in Figure  9, respectively. We also crosschecked our binned QLF with another nonparametric QLF for IMS quasar sample using the C − estimator (Lynden-Bell 1971). The comparison shows  Figure 10. Confidence regions of the best-fit parameters for Φpar in cases 1 (left) and 2 (right). The color-filled regions denote the 1σ (68.3%), 2σ (95.4%), and 3σ (99.7%) confidence levels from darker to brighter. The best-fit values are also marked.  that the QLFs derived from different methods agree with each other (see details in Appendix B). We also derived the parametric QLF at z ∼ 5 (Φ par ), which has a functional form with faint-and bright-end slopes (α and β, respectively; e.g., M18): where Φ * = Φ * z=6 × 10 k(z−6) with a normalization factor at z = 6 (Φ * z=6 ) and k = −0.47 5 (Fan et al. 2001b), and M * 1450 is the break magnitude. Accounting for the lack of bright quasars (M 1450 < −27 mag) in our sample, we used an ancillary SDSS quasar sample (Y16). The IMS and Y16 samples overlap at M 1450 ∼ −27 mag. We simply assigned the two samples to cover different ranges; we used the IMS sample at M 1450 ≥ −27 mag and the Y16 sample at M 1450 < −27 mag. This is because of the potentially underestimated number of Y16 quasars at M 1450 −27 mag, which is lower than that of M18 by a factor of 2-5. For ease of comparison, we re-binned the QLF of the Y16 sample (4.7 < z < 5.4 and M 1450 < −27 mag) and marked them with orange points in Figure 9. Note that we used the completeness functions for the SDSS quasars, provided by Jinyi Yang personally.
We used the maximum likelihood method including selection efficiencies (Marshall et al. 1983) to fit the function with four parameters (Φ * , M * 1450 , α, β) by minimizing S, defined as (4) The first term of this equation is the sum of all the IMS and Y16 samples. The second term stands for the expected number of quasars with given Φ par and p, which is integrated over the ranges of −29.0 < M 1450 < −23.0 and 4.7 < z < 5.4. Therefore, minimizing S is equivalent to finding the appropriate Φ par to satisfy both the observations and expectations. In Figure 9, the bestfit results of Φ par of cases 1 and 2 are shown as red and blue solid lines, respectively, in line with the Φ bin of each case. Note that the uncertainties of the free parameters are calibrated from the ∆S = S − S min distribution, under the assumption of ∆S ∼ χ 2 ν (Lampton et al. 1976), as shown in Figure 10 (Akiyama et al. 2018;Matsuoka et al. 2018b), and this will be discussed in a forthcoming paper (Y. Kim et al. 2020, in preparation). We also obtained the Φ par parameters with fixed bright-end slopes of β = −3.0 and −4.0 (the dotted and dashed lines in Figure 9, respectively). These are consistent with the best-fit results within the 1σ confidence level. All of the derived values are listed in Table 6.
When we compare the estimated QLFs of cases 1 and 2, they are consistent with each other, in both the Φ bin and Φ par cases. The only minor differences between them are observed at M 1450 > −25 mag, especially at the faintest bin. This possible difference at the faintest bin is due to the serendipitous identification of faint quasars, such as IMS J085028−050607 with M 1450 = −23.5 mag at z = 5.36. As shown in Figure  6, the quasar is located at a very low completeness region and the V a of the faintest bin is smaller than the other bins, giving a high number density at the faintest bin after the completeness correction. In fact, in case 2 where this quasar is not included, the Φ bin (−23.25) is more in line with our best-fit QLF. But here, we stress that the faintest bin is also consistent with our best-fit QLF within 1σ confidence level.
Considering the high success rate of our medium-band selection, especially at fainter magnitudes (K19), a future wide-field imaging survey in medium bands will allow us to estimate reliable QLFs at high redshifts even without deep spectroscopic observations. In the following discussion, the best-fit result in case 1 is used as a representative result.
For the IMS quasars only, we also performed testing with a single power-law function: Φ par ∝ 10 −0.4(α+1) . The resultant slopes are α = −1.70 +0.24 −0.26 and −1.67 +0.29 −0.31 in cases 1 and 2, respectively, which are slightly steeper than our best-fit results. These values agree with the recent estimate of α = −2.1 ± 0.7, obtained from a smaller sample of z ∼ 5 faint quasars using a single power-law function (pink line in Figure 11; Shin et al. 2020).
4. DISCUSSIONS 4.1. Reliability of the z ∼ 5 QLF Figure 11 shows various QLFs at z ∼ 5 from the literature (Glikman et al. 2011;Y16;Boutsia et al. 2018;M18;G19;Shen et al. 2020;Shin et al. 2020), in comparison to our best-fit QLFs. To obtain the binned QLF of G19, we averaged their binned QLFs at z = 4.5 and 5.6 (purple crosses). Note that the upper and lower limits of each point indicate the binned QLFs at z = 4.5 and 5.6, respectively. Their parametric QLF (model 2; purple solid line) is also adjusted from z = 4.5 to z = 5 with k = −0.47. Considering that the binned QLFs of Glikman et al. (2011) and Boutsia et al. (2018) are used for the parametric QLF derivation of G19, for ease of comparison, they are represented by purple open circles. Similarly, the z ∼ 4.8 QLF of Shen et al. (2020; gray crosses) is also shifted to z = 5.
While the faint-end slope of our parametric QLF has a flatter value, the binned QLFs in this work are consistent with most QLFs at −25 < M 1450 < −23 except for the Glikman+11, Boutsia+18 Shen+20 Shin+20 Figure 11. QLFs at z ∼ 5. The red filled circles show our binned QLF in case 1, while the red solid line with a shaded region indicates the best-fit parametric QLF with 1σ confidence level. The blue filled circles with the blue solid line are those in case 2. The orange filled circles show the re-binned Y16 QLF consisting of the bright quasars that we used for fitting. The open circles and solid lines are the QLFs from the optical/NIR surveys: Y16 (orange), M18 (teal), and Shin et al. (2020;pink). The X-ray QLFs of G19 (purple) and Shen et al. (2020;gray) are shown as crosses with solid lines, which are shifted to z = 5 (see text). The QLFs of Glikman et al. (2011) and Boutsia et al. (2018) from the optical/NIR surveys, which are used by G19, are represented by purple open circles for ease of comparison. In the case of Shen et al. (2020), we only plotted their X-ray QLFs, while their global QLF (solid line) was derived from not only X-ray but also UV/optical quasars. result of G19 and results from their group in previous years. Possible reasons for the discrepancy have been discussed in detail in Parsa et al. (2018), but here we discuss if the discrepancy is caused by cosmic variance. In this magnitude range, G19 used the z ∼ 4 QLFs (open purple circles) determined from small-area surveys (3.76 deg 2 of Glikman et al. 2011; 1.73 deg 2 of Boutsia et al. 2018). We randomly selected 1.73 deg 2 area in a circular shape from our survey area, and counted the numbers of M 1450 = −23.5 and −24.5 mag quasars (∆M 1450 = 1 mag), which is repeated 100,000 times. Then, we examined the frequency of obtaining the Φ bin values of Boutsia et al. (2018) or higher at the two magnitude bins by chance due to cosmic variance; only 2.2% and 1.4% for M 1450 = −24.5 and −23.5 mag, respectively. If we require the number densities to be high in these two consecutive magnitude bins, then the probability drops to zero, which can be also inferred from Figure 4. We also tested for 5.49 deg 2 area, resulting in the very low probability of 1% even if only one magnitude bin was considered. The result shows that the likelihood of gaining a high quasar number density due to cosmic variance is small, even with an area of 1.73 deg 2 . All quasars that we identified with spectroscopy (K19 and this work) have EW Lyα+NV higher than 15.4Å, the criterion for weak Lyα quasars (Diamond-Stanic et al. 2009). The fraction of weak Lyα quasars is known to be 13.7% at z > 5.7 (Bañados et al. 2016), higher than the fraction of 1.8% from the EW distribution adopted for our calculation of the selection function (Diamond-Stanic et al. 2009). Using the EW distribution of Bañados et al. (2016) for our calculation of the selection function instead (see Section 3.2), we find that our QLF could be underestimated by a factor of ∼ 1.3 at the faint end (both in cases 1 and 2). The exact value of this factor is more likely to be smaller than 1.3, assuming the gradual change in the EW distributions between z ∼ 5 and 6. Future investigation of the EW distributions is desired to improve the QLF estimate.
The G19 QLF at M 1450 > −23 mag mainly relies on X-ray-selected quasars that do not always overlap with those selected in the UV/optical range. One possible explanation for their higher number density than UV/optical quasars is due to the high fraction of obscured X-ray quasars at high redshifts (e.g., 50-80% at 3 < z < 6; Vito et al. 2018), if the obscured fraction applies equally to UV/optical selected quasars. However, recent X-ray QLFs provided by Shen et al. (2020;gray crosses), corrected by their model prediction, also agree with our results, reinforcing the suggestion of the overestimation of G19 even at M 1450 > −22 mag. We note that the Φ bin of M18 has a value lower than those of our and other QLFs at the faintest bin they explored (M 1450 > −23 mag). QLFs can be over-or underestimated easily at the faintest end depending on how well the selection function is constructed. We suggest this as a possible reason for the discrepancy.

Contribution of z ∼ 5 Quasars to Ionizing Background
We followed the method in Kulkarni et al. (2019) to calculate the UV ionizing emissivity of z ∼ 5 quasars at 1450 ( 1450 ) and 912Å ( 912 ), which can be written as Haardt & Madau 12 Kulkarni+19 Shen+20 Haardt & Madau 12 Kulkarni+19 Shen+20 respectively. The thick solid gray lines show the emissivity values required to fully ionize the IGM as a function of redshift, for a given clumping factor C (=2, 3, 5) and fesc = 1 (Madau et al. 1999). The right-hand vertical axis shows the correspondingṅion for a given 912, assuming the UV spectral shape of Lusso et al. (2015) and fesc = 1. The gray squares show the photon density required to maintain the ionization state of IGM, predicted from the transmitted Lyα flux measurements (D'Aloisio et al. 2018).
To evaluate 912 from 1450 , we assume a double powerlaw shape of the quasar UV spectra (f ν ∝ ν −αν , where α ν = 0.61 if λ ≥ 912Å and α ν = 1.70 if λ < 912Å; Lusso et al. 2015). For the best-fit case 1 QLF, 1450 and 912 are 6.64 +2.81 −1.70 (or 4.94 to 9.45) and 5.00 +2.12 −1.28 (or 3.72 to 7.12), respectively, in units of 10 23 erg s −1 Hz −1 Mpc −3 . Figure 12 shows the observed 912 values of high-redshift quasars. Our result is the lowest value at z ∼ 5, an order of magnitude lower than Giallongo et al. (2015) estimates (brown open symbols). Considering large errors, most of the 912 values are consistent with the synthesis model of Haardt & Madau (2012) and the best-fit model of Shen et al. (2020). Notable exceptions to the above trend are the results of Giallongo et al. (2015) and G19. However, the former were derived from X-ray candidates without spectroscopic confirmation, while the latter are also less reliable than others due to the possible overestimation of their QLF as mentioned earlier in this section. The above results, however, could be sensitive to the uncertain extrapolation of the QLF down to faint magnitudes of M 1450 = −23 to −18 mag. If we integrate Equation (5) up to M 1450 = −23 mag, at which our QLF is well constrained, 912 is 4.48 × 10 23 erg s −1 Hz −1 Mpc −3 . This is consistent with 912 integrated up to M 1450 = −18 mag within 1σ error, meaning that the additional contribution of quasars with M 1450 > −23 mag is not significant unless their number density is unexpectedly high.
Despite the uncertainty in our emissivity estimates at the faintest end, we can explore the significance of the contribution of quasars to the UV ionizing radiation. Under the same assumptions, we calculate the number density of ionizing photons, where f esc is an escape fraction assumed to be unity (f esc = 1; e.g., Grazian et al. 2018), and ξ ion is the ionizing photon production efficiency of a quasar, normalized at 1450Å. We obtain ξ ion 6.05 × 10 25 Hz erg −1 with the quasar spectral shape of Lusso et al. (2015). Using the 1450 value from our z ∼ 5 QLF, we obtainṅ ion of (4.02 +1.70 −1.03 ) × 10 49 Mpc −3 s −1 . At z = 5, the required (critical) photon density to keep balance with hydrogen recombination isṅ cri ion ∼ 6.3×10 49 C Mpc −3 s −1 (Madau et al. 1999;Bolton, & Haehnelt 2007; thick gray lines in Figure 12), where C is the H II clumping factor recently predicted as C < 5 at z = 5 (Bolton, & Haehnelt 2007;Shull et al. 2012;Jeeson-Daniel et al. 2014;D'Aloisio et al. 2020). For a plausible clumping factor of C = 3, our result suggests that z ∼ 5 quasars radiate about 21 +9 −5 % of the UV ionizing photons required to balance the ionized state of hydrogen at that time. This fraction can reach 45% if C = 2 andṅ ion is the 1σ upper value of 5.72 × 10 49 Mpc −3 s −1 . In contrast, if we adopt the 1σ lower limit, the quasar contribution is 9%. Therefore, z ∼ 5 quasars alone radiate only up to ∼ 50% of the minimum requirement of UV photons to balance with hydrogen recombination rates (ṅ cri ion of Madau et al. 1999), under the extreme assumptions given above.
An alternative method to estimateṅ cri ion is by using observational inference, such as the Lyα forest. Previous studies calculated the photoionization rate of the UV ionizing background required to match the mean transmitted Lyα flux of high-redshift quasars (Bolton, & Haehnelt 2007;Wyithe & Bolton 2011;D'Aloisio et al. 2018), and the corresponding critical emissivity cri 912 . Theṅ cri ion can then be approximated toṅ cri ion ( cri 912 /h p α ν ) Mpc −3 s −1 (Meiksin 2005;Becker & Bolton 2013), where h p is the Planck constant and α ν = −1.70 (Lusso et al. 2015). In Figure 12, the gray squares indicate the predictedṅ cri ion values with the recently estimated cri 912 values of D' Aloisio et al. (2018). At z = 5, we obtainṅ cri ion = 6.4 +1.1 −2.3 × 10 50 erg s −1 , indicating that the quasar contribution with our QLF is only approximately 6 +8 −2 % (or 4-14%) considering the 1σ errors oḟ n ion andṅ cri ion . Such a low contribution is similar to that at z ∼ 6, where the quasar contribution is considered to be less than 10% (Willott et al. 2010;Kashikawa et al. 2015;Kim et al. 2015;Onoue et al. 2017;Matsuoka et al. 2018b). However, note that the cri 912 determination from the Lyα opacity of high-redshift quasars is sensitive to the mean free path of hydrogen photons, which can also be affected by the main ionizing sources (D'Aloisio et al. 2018).
Unlike at z ∼ 6, there is a discrepancy between thė n cri ion from the two methods at z = 5. This is due to the increases in the number of ionizing sources and the volume filled with ionized hydrogen from z = 6 to 5. Considering this difference, we conclude that quasars are likely to occupy only a minor portion of the total UV ionizing background at z = 5, estimated from the Lyα opacity of high-redshift quasars. However, they can provide nearly half the quantity of photons required to balance with recombination rates. A better understanding of the quasar contribution for z = 5 IGM ionizing photons requires the construction of a QLF at M 1450 > −23 mag where we rely on extrapolation to estimate the QLF shape.
In addition, improving constraints on C and other physical properties of high-redshift quasars (e.g., f esc ) is also required. For instance, our calculation ofṅ ion /ṅ cri ion depends on the assumed UV spectral slope of quasars. The assumed value of α λ = α ν − 2 = −1.39 (Lusso et al. 2015) is based on quasars at z = 2.4 and may not represent UV spectral slopes of z = 5 quasars in this work. For example, Mazzucchelli et al. (2017) find α λ = −1.6 ± 0.1 forz = 6.6 quasars, which is steeper than the Lusso et al. (2015) value. We directly infer the UV slope of our sample using the SED-fitting result of a sub-sample of our z = 5 quasars in K19. Using the values presented in Table 5 of K19, we find α λ = −1.9±0.1, which is the steepest UV slope among other results from various quasar samples (Laor et al. 1997;Vanden Berk et al. 2001;Telfer et al. 2002;Stevans et al. 2014;Selsing et al. 2016). If we adopt this rather extreme slope of −1.9, the exponent in Eq. (6) changes to 0.1. However, the UV emissivity increases only by a factor of 1.27. So, the maximal UV ionizing photon fraction changes from 45% to 57%-but this can be considered as only a modest increase.

SUMMARY
We present the advanced results of our z ∼ 5 quasar survey with IMS (K19). Our findings are as follows.
1. Based on the follow-up spectroscopy carried out with GMOS on the Gemini-South 8 m Telescope, we newly identified eight z ∼ 5 quasars in the IMS survey area. Using our high-redshift quasar model, we measured their redshifts and magnitudes in the ranges of 4.71 ≤ z ≤ 5.15 and −26.2 ≤ M 1450 ≤ −23.3.

2.
Considering the survey completeness with our selection criteria, we derived the binned and parametric z ∼ 5 QLFs. The selection criteria are considered either with only broadband colors (case 1; 43 quasars) or broadband and medium-band colors (case 2; 32 quasars). Including the SDSS bright quasar sample of Y16, the parametric QLFs in both cases are well determined without any fixed parameters.
3. We find a relatively flat faint-end slope of −1.2 compared to previous studies at z ∼ 5, although our 1σ limit allows steeper slopes down to α ∼ −2. We calculated a low ionizing emissivity of 912 = 5.0 × 10 23 erg s −1 Hz −1 Mpc −3 and a number density of UV ionizing photons ofṅ = 4.0 × 10 49 Mpc −3 s −1 . This implies that the quasars do not produce all of the ionizing photons at z ∼ 5.
Our z ∼ 5 QLF that is well defined at M 1450 < −23 mag with a better parameter estimation than previous surveys, but our imaging survey depths limits us from determining the QLF at M 1450 > −23 mag. Hence, our conclusion regarding the contribution of quasars to the IGM-ionizing photons relies on the uncertain extrapolation of our QLF to M 1450 > −23 mag. The above conclusion may change significantly if there is an abrupt increase in the number of quasars at the faintest limit such as a large number of obscured quasars at high redshifts (e.g., Aird et al. 2015). Deeper imaging data are necessary to settle down this issue, which could come from the Legacy Survey of Space and Time (LSST) and future space telescope missions. Reliable identification of these faint quasars may be challenging even with the future 30-m class telescopes. But even in such a challenging regime, we expect that our medium-band technique will be able to identify faint high-redshift quasars if used on 4 to 8-m class telescopes. This paper includes data obtained at the Gemini Observatory, acquired through the Gemini Science Archive, and processed using the Gemini IRAF package, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). This paper includes data obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at Terapix available at the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS. The United Kingdom Infrared Telescope (UKIRT) is supported by NASA and operated under an agreement between the University of Hawaii, the University of Arizona, and Lockheed Martin Advanced Technology Center; operations are enabled through the cooperation of the Joint Astronomy Centre of the Science and Technology Facilities Council of the U.K. This paper includes data gathered with the 6.5 m Magellan telescopes located at Las Campanas Observatory, Chile. This paper includes data taken at The McDonald Observatory of The University of Texas at Austin. We would like to thank Editage (www.editage.co.kr) for English language editing.  Figure 13. Left: IMACS optical spectra of two of the IMS z ∼ 5 quasars that were newly discovered in this work. IMS J021507−051849 was not detected with IMACS. Right: GMOS-S optical spectra of two non-quasars. They have no Lyα break with continuum emissions. The shaded regions represent the bad columns, hot pixels, CCD gaps, or wavelength ranges that are not covered by the observational configuration.  Here, we describe the optical spectra of our quasars/candidates, which are not included in the analysis in this paper. The observations are summarized in Table 7, and their spectra are shown in Figure 13.
Of the eight newly discovered z ∼ 5 quasars, IMS J021507−051849 and IMS J221054+013430 were also observed with the Inamori-Magellan Areal Camera and Spectrograph (IMACS; Dressler et al. 2011) on the Magellan Baade 6.5 m Telescope. The observing configuration was set to follow that used in K19: using a grating of 150 lines mm −1 with a resolution of R ∼ 600, OG570 filters to avoid a zeroth-order overlap, and spectral/spatial binning of 2 × 2. We followed the general steps for spectral reduction, except for the flux calibration with standard stars. Instead, we scaled the fluxes with their i -band magnitudes. These are excluded in the main text because the GMOS with which our fiducial spectra were obtained is more sensitive than IMACS.
We also present the two candidates that were identified as non-quasars with GMOS-S: IMS J090540−011038 and IMS J090540−025524. The observing configurations were as described in Section 2 and in K19. In their 1D and 2D spectra, there is no Lyα break but continuum emission across the whole wavelength range, meaning that they are not high-redshift quasars.

B. NONPARAMETRIC DETERMINATION OF QLF
The QLFs determined in the main text are based on the small number statistics, so the methods we used could affect the exact QLF estimations. For instance, the 1/V max method critically depends on how we bin the sample especially in the case with a handful of quasars. In this section, we describe the nonparametric determination of QLF, using the Lynden-Bell's C − method (Lynden-Bell 1971) that requires no assumption on the distribution for the one-side-truncated data sets. Thus, it is suitable for our quasar samples and the survey limits; the brightest quasar in our sample (IMS J220107+030208; i = 19.1 mag) is much fainter than the saturation limit of CFHTLS (i ∼ 17 mag).
QLF is canonically described in the form of a bivariate function of magnitude and redshift, e.g., Φ(M, z). However, our IMS quasars are in the narrow redshift range of 4.7 < z < 5.4, so what we actually derived in Section 3.4 is similar to the univariate function of magnitude at the specific redshift of z = 5. Hence, in this narrow redshift range, the marginal distribution of the bivariate QLF in the magnitude direction, ψ(M ), can be approximated by ψ(M ) ≈ Φ(M, z). Note that we performed the standard correlation test between the magnitude and redshift of our sample with the Kendall's τ statistic (see details in Efron & Petrosian 1992;Fan et al. 2001b;Schindler et al. 2018Schindler et al. , 2019, giving τ = 0.04 which means they are uncorrelated at ∼ 1σ level. With the Lynden-Bell's C − estimator (Lynden-Bell 1971), the cumulative luminosity function Ψ(M ) = ψ(M )dM can be recovered: for the sample sorted by magnitude (M 1 < ... < M i−1 < M i < ... < M N ). The estimator C − (M i ) indicates the total number of object of which magnitudes are brighter than M i (but the object i itself is omitted). But the problem is that this direct usage of C − (M i ) can be applied to a sample with a sharp boundary for selection, while our selection probability has a smooth boundary in the magnitude direction. Instead of C − (M i ), Fan et al. (2001a) generalized the total weighted number of quasars N i for arbitrary selection functions (see also Schindler et al. 2018Schindler et al. , 2019. Given a selection probability function p(M, z), one can construct the quantity N i :  Figure 14 shows the ψ(M 1450 ) of the IMS quasars from the 100 bootstrap samples. While the C − method can recover the shape of the marginal luminosity function, it does not provide any information related to its normalization. Therefore, the marginal distribution is normalized by the total number of quasars expected, Ψ(M 1450 < −23) = −23 −∞ ψ(M 1450 )dM 1450 . We also plot our binned and parametric QLFs, also normalized by the cumulative luminosity function of the parametric QLF at M 1450 < −23 mag. They are consistent with the marginal distribution ψ(M 1450 ) within 1σ errors. Therefore, we can confidently say that our binning for the binned QLF is reliable, and the fitting for the parametric QLF as well.