Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs). V. Quasar Luminosity Function and Contribution to Cosmic Reionization at z = 6

We present new measurements of the quasar luminosity function (LF) at z ∼ 6 over an unprecedentedly wide range of the rest-frame ultraviolet luminosity M1450 from −30 to −22 mag. This is the fifth in a series of publications from the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs) project, which exploits the deep multiband imaging data produced by the Hyper Suprime-Cam Subaru Strategic Program survey. The LF was calculated with a complete sample of 110 quasars at 5.7 ≤ z ≤ 6.5, which includes 48 SHELLQs quasars discovered over 650 deg2 and 63 brighter quasars discovered by the Sloan Digital Sky Survey and the Canada–France–Hawaii Quasar Survey (including one overlapping object). This is the largest sample of z ∼ 6 quasars with a well-defined selection function constructed to date, which has allowed us to detect significant flattening of the LF at its faint end. A double power-law function fit to the sample yields a faint-end slope , a bright-end slope , a break magnitude , and a characteristic space density Gpc−3 mag−1. Integrating this best-fit model over the range −18 < M1450 < −30 mag, quasars emit ionizing photons at the rate of s−1 Mpc−3 at z = 6.0. This is less than 10% of the critical rate necessary to keep the intergalactic medium ionized, which indicates that quasars are not a major contributor to cosmic reionization.


Introduction
The first billion years of the universe, corresponding to redshift z>5.7, have been the subject of major observational and theoretical studies in the last few decades. The first generation of stars, galaxies, and supermassive black holes (SMBHs) are thought to have formed during this epoch, and the universe became reionized during that time, most likely due to the ionizing photons from these light sources. A large number of high-z galaxies and galaxy candidates have been identified up to z∼10 and beyond, and the evolution of the galaxy luminosity function (LF) has been intensively studied (e.g., Bouwens et al. 2011Bouwens et al. , 2015McLeod et al. 2016;Oesch et al. 2016Oesch et al. , 2018Ishigaki et al. 2018). Robertson et al. (2015) demonstrated that these high-z galaxies produced sufficient quantities of ionizing photons to dominate the reionization process, based on the Planck measurements of the cosmic microwave background polarization (Planck Collaboration et al. 2016) and an assumed value of the Lyman continuum escape fraction.
The search for high-z quasars 27 has also undergone significant progress in recent years, thanks to the advent of wide-field (1000 deg 2 class) multiband red-sensitive imaging surveys such as SDSS (York et al. 2000), the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS), the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS1; Chambers et al. 2016), and the United Kingdom Infrared Telescope Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007). At the time of writing of this paper, there were 242, 145, 18, and 2 quasars reported in the literature at redshifts beyond z=5.7, 6.0, 6.5, and 7.0, respectively. The two highest-z quasars were found at z=7.09 (Mortlock et al. 2011) and z=7.54 (Bañados et al. 2018). The quasar LF at z=6 has been measured with the complete samples of quasars from SDSS (Jiang et al. 2016) and the Canada-France-Hawaii Quasar Survey (CFHQS; Willott et al. 2010) based on CFHTLS. However, the above measurements were limited mostly to M 1450 <−24 mag, where the LF is approximated by a single power law, with only a single CFHQS quasar known at a fainter magnitude (M 1450 =−22.2 mag). Thus, it has remained unclear whether or not the LF has a break, and what the faint-end slope is if the break exists. This is a critical issue, since the faint-end shape of the LF reflects a more typical mode of SMBH growth than probed by luminous quasars, and it has a direct impact on the estimate of the quasar contribution to cosmic reionization.
In the past few years, there have been several attempts to find low-luminosity quasars at z∼6. Kashikawa et al. (2015) found two quasars (one of which may in fact be a galaxy) with M 1450 ∼−23 mag over 6.5 deg 2 imaged by Suprime-Cam (Miyazaki et al. 2002), a former-generation wide-field camera on the Subaru 8.2 m telescope. The number densities derived from these two (or one) quasars and the faintest CFHQS quasar may point to a flattening of the faint-end LF, but the small sample size hampered accurate measurements of the LF shape. Onoue et al. (2017) took over the analysis of the above Suprime-Cam data, but found no additional quasars, confirming the number density measured by Kashikawa et al. (2015). On the other hand, Giallongo et al. (2015) reported Chandra X-ray detection of five very faint active galactic nuclei (AGNs) at z∼6, with −19M 1450 −21 mag, over 170 arcmin 2 of the Great Observatories Origins Deep Survey (Giavalisco et al. 2004) field. This surprisingly high detection rate could indicate a significant AGN contribution to cosmic reionization. However, their results have been challenged by a number of independent deep X-ray studies, finding much lower number densities of faint AGNs (e.g., Weigel et al. 2015;Cappelluti et al. 2016;Vito et al. 2016;Ricci et al. 2017;Parsa et al. 2018). A high number density of high-z faint AGNs may also be in tension with the epoch of He II reionization inferred from observations (D'Aloisio et al. 2017;Khaire 2017;Mitra et al. 2018).
There have also been extensive efforts to measure the quasar LF at lower redshifts, e.g., at z∼4 (Glikman et al. 2011;Ikeda et al. 2011;Masters et al. 2012) and at z∼5 (Ikeda et al. 2012;McGreer et al. 2013;Yang et al. 2016). Recently, Kulkarni et al. (2018) reanalyzed a large sample of quasars compiled from the above and other papers, and reported very bright break magnitudes ( * < -M 27 1450 mag) with steep faint-end slopes at 4z6. On the other hand, more recent data reaching ∼1 mag fainter than the previous measurements seem to suggest that the LF breaks at fainter magnitudes both at z∼4 (Akiyama et al. 2018) and z∼5 (McGreer et al. 2018; see the discussion in Section 4 of this paper).
This paper presents new measurements of the quasar LF at z∼6, exploiting a complete sample of 110 quasars at 5.7z6.5. The sample includes 48 low-luminosity quasars recently discovered by the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs; Matsuoka et al. 2016) project. SHELLQs rests on the Subaru Strategic Program (SSP) survey (Aihara et al. 2018a) with the Hyper Suprime-Cam (HSC; Miyazaki et al. 2018), a wide-field camera mounted on the Subaru telescope. We are carrying out follow-up spectroscopy of high-z quasar candidates imaged by the HSC and have so far identified 150 candidates over 650 deg 2 , which include 74 high-z quasars, 25 high-z luminous galaxies, 6 [O III] emitters at z∼0.8, and 45 Galactic cool dwarfs (Matsuoka et al. 2016(Matsuoka et al. , 2018a(Matsuoka et al. , 2018bY. Matsuoka et al. 2018, in preparation). We are also carrying out near-infrared (near-IR) spectroscopy and Atacama Large Millimeter/submillimeter Array (ALMA) observations of the discovered objects. The first ALMA results were published in Izumi et al. (2018), and further results are in preparation.
This paper is organized as follows. In Section 2, we describe our quasar sample to establish the LF, drawn from SDSS, CFHQS, and SHELLQs. The completeness of the SHELLQs quasar selection is evaluated in Section 3. The binned and parametric LFs are presented and discussed in Section 4, and the quasar contribution to cosmic reionization is estimated in Section 5. A summary appears in Section 6. We adopt the cosmological parameters H 0 =70 km s −1 Mpc −1 , Ω M =0.3, and Ω Λ =0.7. All magnitudes are presented in the AB system (Oke & Gunn 1983) and are corrected for Galactic extinction (Schlegel et al. 1998). In what follows, we refer to z-band magnitudes with the AB subscript ("z AB "), while the redshift z appears without a subscript.

Quasar Sample
We derive the quasar LF with a complete sample of 110 quasars at 5.7z6.5, as summarized in Table 1 and plotted in Figure 1. These quasars are drawn from SDSS, CFHQS, and SHELLQs, which roughly cover the bright, middle, and faint portions of the magnitude range we probe (−22<M 1450 < −30 mag), respectively. 28 Table 2 lists the number of objects in 27 Throughout this paper, "high-z" denotes z>5.7, where the cosmic age is less than a billion years and objects are observed as i-band dropouts in the Sloan Digital Sky Survey (SDSS) filter system (Fukugita et al. 1996). 28 The present measurements do not include the bright quasars discovered by Pan-STARRS1 (Bañados et al. 2016), whose selection completeness has not been published yet. each M 1450 bin used for the LF calculation and the corresponding survey volumes (V a ; see below).

SDSS
We exploit a complete sample of 47 SDSS quasars at 5.7z6.5, presented in Jiang et al. (2016). Of these, 24 quasars with z AB 20 mag were discovered in the SDSS main survey, using single-epoch imaging data with 54 s exposures. Seventeen quasars (in which seven quasars were also found in the main survey) with 20z AB 20.5 mag were discovered in the SDSS overlap regions, where two or more exposures were taken, due to the scanning strategy and repeated observations of some fields in the main survey. The remaining 13 quasars with z AB 22 mag were discovered in SDSS Stripe 82 on the celestial equator, which was repeatedly scanned 70-90 times (Annis et al. 2014;Jiang et al. 2014). In total, these 47 quasars span the magnitude range from M 1450 =−30 to −24 mag. The absolute magnitudes (M 1450 ) were estimated by extrapolating the continuum spectrum redward of Lyα to rest-frame 1450 Å, assuming a power-law shape f λ ∝λ −1.5 (except for a few quasars, whose observed spectra covered that rest-frame wavelength, or whose near-IR spectra provided estimates of the continuum slope). The effective areas of the main, overlap, and Stripe 82 surveys are 11,240, 4223, and 277 deg 2 , respectively.
The selection completeness was estimated with model quasars, which were created using spectral simulations presented in McGreer et al. (2013). The models were designed to reproduce the observed colors of ∼60,000 quasars at 2.5<z<3.5 in the SDSS Baryon Oscillation Spectroscopic Survey (Ross et al. 2012) and took into account the observed relations between  Jiang et al. (2016) for the SDSS quasars, in Willott et al. (2010) for the CFHQS quasars, and in our previous papers for the SHELLQs quasars. J231546.58-002357.9 was also recovered by CFHQS and SHELLQs, and is hence included in the complete samples of all the three surveys. Five quasars in the SHELLQs sample (J021013.19-045620.8, J022743.29-060530.3, J121721.34+013142.6, J220417.92+011144.8, and J221917.22+010249.0) were originally discovered by other surveys (see Table 1 of Matsuoka et al. 2018a for details), but are not included in the SDSS or CFHQS complete sample.
spectral features and luminosity, such as the Baldwin effect. The effect of IGM absorption was modeled using the prescription of Worseck & Prochaska (2011) extended to higher redshifts with the data from Songaila & Cowie (2010) and was checked against the measurements of Songaila (2004) and Fan et al. (2006). The electronic data of the completeness functions of each of the three surveys were kindly provided by Linhua Jiang (2017, private communication)

CFHQS
We use a complete sample of 17 CFHQS quasars at 5.7z6.5, presented in Willott et al. (2010). Of these, 12 quasars were discovered in the Red-sequence Cluster Survey 2 (RCS-2) field observed with the MegaCam on CFHT, with exposure times of 500 and 360 s in the i-and z-bands, respectively. Four quasars were discovered in the CFHTLS Very Wide (VW) field, imaged for 540 and 420 s in the MegaCam i-and z-bands, respectively. These 16 quasars ("CFHQS-wide quasars" hereafter) span the magnitude range from M 1450 =−27 to −24 mag. The remaining quasar, with M 1450 =−22.2 mag, was discovered in the CFHQS-deep field, which is a combination of the CFHTLS-Deep and the Subaru XMM-Newton Deep Survey (SXDS) fields. The effective areas of the CFHQS-wide (RCS-2 + CFHTLS-VW) and deep (CFHTLS-Deep + SXDS) fields are 494 and 4.47 deg 2 , respectively. The selection completeness was estimated with quasar models created from the observed spectra of 180 SDSS quasars at 3.1<z<3.2. The effect of IGM absorption was incorporated based on the data taken from Songaila (2004). The electronic data of the completeness functions were kindly provided by Chris Willott (2012, private communication).
The absolute magnitudes (M 1450 ) of the CFHQS quasars were originally estimated from the observed J-band fluxes with a template quasar spectrum. For consistency with the measurements in SDSS and SHELLQs, we remeasured their M 1450 by extrapolating the continuum spectrum redward of Lyα, assuming a power-law shape f λ ∝ λ −1.5 . The resultant M 1450 values differ from the original (CFHQS) values by −0.4to+0.2 mag for all but one quasar; the exception is the faintest quasar J021627.81-045534.1, for which the new measurement indicates 0.7 mag fainter continuum luminosity than in the original measurement. This quasar has an unusually strong Lyα line, contributing about 70% of the observed z-band flux (Willott et al. 2009). It has a similar z−J color to other high-z quasars despite the strong contribution of Lyα to the z-band flux, suggesting that the J-band also has significant contribution from strong lines like C IV λ1549. If so, the continuum flux is significantly fainter than the J-band magnitude would indicate.

SHELLQs
We use 48 SHELLQs quasars at 5.7z6.5, discovered from the HSC-SSP Wide survey fields. HSC is a wide-field camera mounted on the Subaru Telescope (Miyazaki et al. 2018). It has a nearly circular field of view of 1°.5 diameter, covered by 116 2K×4K fully depleted Hamamatsu CCDs, with a pixel scale of 0 17. The HSC-SSP survey (Aihara et al. 2018a) has three layers with different combinations of area and depth. The Wide layer is observing 1400 deg 2 in several discrete fields mostly along the celestial equator, with 5σ point-source depths of (g AB , r AB , i AB , z AB , y AB )=(26.5, 26.1, 25.9, 25.1, 24.4) mag measured in 2 0 apertures. The total exposure times range from 10 minutes in the g-and r-bands to 20 minutes in the i-, z-, and y-bands, divided into individual exposures of ∼3 minutes each. The Deep and the UltraDeep layers are observing smaller areas (27 and 3.5 deg 2 ) down to deeper limiting magnitudes (r AB =27.1 and 27.7 mag, respectively). Data reduction was performed with the dedicated pipeline hscPipe . We use the point-spread function (PSF) magnitude (m PSF,AB , or simply m AB ) and the CModel magnitude (m CModel,AB ), which are measured by fitting the PSF models and two-component, PSFconvolved galaxy models to the source profile, respectively (Abazajian et al. 2004;Bosch et al. 2018). We utilize forced photometry, which measures source flux with a consistent aperture in all bands. The aperture is usually defined in the z-band for i-band dropout sources, including high-redshift quasars. A full description of the HSC-SSP survey may be found in Aihara et al. (2018a).
The SHELLQs quasars used in this work were drawn from the HSC-SSP Wide survey fields. While the candidate selection procedure has changed slightly through the course of the survey, we defined a single set of criteria to select the 48 objects. We first queried the "S17A" internal data release (containing all the data taken before 2017 May) of the SSP survey, with the following conditions:  The first line defines the selection limits of magnitude, photometry signal-to-noise ratio (S/N), and color, while the second line rejects apparently extended objects (see Matsuoka et al. 2016, and the following section). The merge.peak flag is true (t) if the source is detected in the specified band and false (f) if not. The quasars in the present complete sample are required to be observed in the i-, z-, and y-bands (but not necessarily in the g-or r-band), and to be detected both in the zand y-bands. The condition on the inputcount.value flag requires that the query is performed on the fields where two or more exposures were taken in each of the z-and y-bands. The last five conditions reject sources on the pixels that are close to the CCD edge, saturated, affected by cosmic rays, registered as bad pixels, or close to bright objects, in any of the i-, z-, or y-bands. The sources selected above were matched, within 1 0, to near-IR sources from the UKIDSS (Lawrence et al. 2007) and Visible and Infrared Survey Telescope for Astronomy (VISTA) Kilo-degree Infrared Galaxy (VIKING) surveys (Edge et al. 2013). We then calculated a Bayesian probability (P Q B ) for each candidate being a quasar rather than a Galactic brown dwarf (BD), based on models for the spectral energy distribution (SED) and surface density as a function of magnitude (see Matsuoka et al. 2016 for details). Our algorithm does not include galaxy models at present. We consider those sources with in the list of candidates for spectroscopy. Only ∼10% of the final SHELLQs quasars have near-IR counterparts in practice, and they would have been selected as candidates with HSC photometry alone; the near-IR photometry is mainly used to reject contaminating BDs, which have much redder optical to near-IR colors than do high-z quasars.
Finally, the candidates went through a screening process using the HSC images. We first used an automatic algorithm with Source Extractor (Bertin & Arnouts 1996) to remove apparently spurious sources (e.g., cosmic rays, transient objects, and CCD artifacts). The algorithm rejects those sources whose photometry (in all the available bands) is not consistent within 5σ error between the stacked and individual pre-stacked images, and those sources whose shapes are too compact, diffuse, or elliptical to be celestial point sources. We checked a portion of the rejected sources and confirmed that no real, stable sources were rejected in this automatic procedure. Indeed, we adopted conservative rejection criteria here, so that any ambiguous cases were passed through to the next stage. The remaining candidates were then screened by eye, which removed additional problematic objects (mostly cosmic rays and transient sources). The automatic procedure rejected >95% of the input candidates, and ∼80% of the remaining candidates were removed by eye.
The final spectroscopic identification is still underway, but has now been completed down to a limiting magnitude of  z 24.0 AB splim mag. The actual z AB splim values vary from field to field, depending on the available telescope time when the individual fields were observable, and are summarized in Table 3. In total, 48 quasars with  z z AB AB splim and spectroscopic redshifts 5.7z6.5 were selected as the complete sample for the present work. The remaining SHELLQs quasars were not in the sample because they are fainter than z AB splim , outside the above redshift range, or fail to meet one or more of the criteria listed in Equation (1). The absolute magnitudes (M 1450 ) were estimated in the same way as used for the SDSS quasars (see above).
The effective survey area was estimated with a random source catalog stored in the HSC-SSP database (Coupon et al. 2018). The random points are placed over entire survey fields, with surface density of 100 arcmin −2 , and each point contains the survey information at the corresponding position (number of exposures, variance of background sky, pixel quality flags, etc.) for each filter. We queried this random catalog with the pixel flag conditions presented in Equation (1). The number of output points was then divided by the input surface density, giving the effective survey area as listed in Table 3.
The SDSS, CFHQS, and SHELLQs samples contain one quasar in common (J231546.58-002357.9). This quasar is treated as an independent object in each of the individual survey volumes in order not to underestimate the number density.  Notes. M 1450 and DM 1450 represent the center and width of each magnitude bin, respectively. The numbers in the parentheses represent the cosmic volumes contained in the individual surveys (V a ; see Equation (6)

SHELLQS Completeness
The SHELLQs quasar selection is known to be fairly complete at bright magnitudes, to which past wide-field surveys (such as SDSS and CFHQS) were sensitive. The HSC-SSP S17A survey footprint contains eight previously known high-z quasars with , and our selection recovered seven of them. The remaining quasar is blended with a foreground galaxy, which boosted the i-band flux of the quasar measured by the HSC pipeline and caused it to be rejected. We evaluate the actual selection completeness in this section.

Source Detection
Source detection in the HSC data processing pipeline mag. The z-band detection completeness is thus expected to be close to 100% for the quasars in our complete sample, which are brighter than = z 23.8 24.2 AB splim mag. We tested the detection completeness in each band with simulations, in which artificial point sources were inserted on random positions of the stacked HSC images, and then recovered with hscPipe. The input source models were created with the PSFs measured at each image position. The same simulations were used in Aihara et al. (2018b) to evaluate the detection completeness of the HSC-SSP Public Data Release 1. 29 These simulations were performed on 180 12′×12′ patches selected randomly from the survey area (the computer time required to run over the entire survey area would have been prohibitively long). The recovery rate of the input sources, as a function of magnitude, is then fitted with a function (Serjeant et al. 2000): where f max , f min , α, and m AB 50 represent the detection completeness at the brightest and faintest magnitudes, the sharpness of the transition between f max and f min , and the magnitude at which the detection completeness is 50%, respectively.
The resultant completeness functions are presented in Figure 3. Overall, they have similar shapes to each other, except for varying depths from patch to patch. It is worth noting that the completeness at the faintest magnitudes ( f min ) is higher than zero, which is due to the chance superposition of   input sources with true sources in the original HSC images used. Figure 4 compares the m AB 50 values with the 5σ limiting magnitudes ( s m AB 5 ) described above. These two quantities agree very well with each other, as expected given that the hscPipe detection threshold is approximately equivalent to < s m m AB AB 5 . Based on the above measurements and simulations, we quantified the detection completeness in the z-and y-bands over the entire survey area, as follows. For each 12′×12′ patch ("p"), the completeness functions were defined using Equation (2). We retrieved s z AB 5 and s y AB 5 from the survey database and used them as surrogates for z AB 50 and y AB 50 in the individual patches. The parameters f max and f min were fixed to 1.0 and 0.0, respectively. Finally, we assumed α=2.4, the median value measured in both the z-and y-bands for the 180 patches in which we ran the simulations (the dispersion in this quantity measured by the median absolute deviation is Δα∼0.4 in both bands). We checked that the present results are not sensitive to the choice of α, since the detection completeness is close to 100% at the present magnitude limit of z AB <24.2 mag.

Point-source Selection
The SHELLQs algorithm uses the criterion to identify point sources from the HSC database. The completeness of this selection was evaluated with a special HSC data set on the COSMOS field, one of the two UltraDeep fields of the SSP survey, for which we have many more exposures than in a Wide field. This data set was created by stacking a portion of the UltraDeep data taken during the best, median, or worst seeing conditions to match the Wide depth. We selected stars on this field with the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) catalog (Leauthaud et al. 2007) and measured the fraction of stars meeting Equation (3). The results are presented in Figure 5.
The completeness of our point-source selection is close to 100% at bright magnitudes and decreases mildly to 90% at z AB ∼24.0 mag. No significant difference was observed between the different seeing conditions at z AB <24 mag. We fitted the above results for the median seeing with Equation (2) and obtained the best-fit parameters ( f max , f min , α, z AB 50 )=(1.00, 0.72, 0.76, 24.5). This best-fit function, ( ) f z ps AB , is used to simulate the selection completeness of point sources in the following.
On the other hand, we found that the effect of resolved host galaxies on our quasar selection is negligible. This was simulated as follows. Since the luminosities of high-z quasar host galaxies are unknown, we assumed the following, based on the low-z results for SDSS quasars with similar nuclear luminosity to the SHELLQs quasars (Matsuoka et al. 2014(Matsuoka et al. , 2015: (i) the typical host galaxy luminosity ranges from M UV =−18 to −21 mag (corresponding tõ z 25.5 28.5 CModel,AB mag at z=6), and (ii) there is no correlation between the nuclear and host galaxy luminosities. The host galaxies were simulated with a sample of Lyman Break Galaxies (LBGs) at z∼6, found from the HSC-SSP Wide data Ono et al. 2018 We found that the unresolved fraction is 100% at AGN magnitudes z AB <25.0 mag and decreases to 90% at 26.0 mag.  The open circles, dots, and crosses represent the best, median, and worst seeing conditions, respectively. The best-fit function (Equation (2)) to the median seeing data is represented by the dashed curve.
We thus conclude that our point-source selection loses only a negligible fraction of quasars due to the resolved host galaxies, at the present magnitude limit of z AB <24.2 mag.
Here we note that compact galaxies could have -< z z 0.15 AB CModel,AB and contaminate our quasar candidates. Indeed, so far we have discovered 25 high-z galaxies in addition to 74 high-z quasars from the HSC candidates. However, the present work uses only spectroscopically confirmed quasars and thus is not affected by galaxy contamination.

Foreground Flux Contamination
As we wrote previously, we failed to recover one of the eight previously known quasars in our survey footprint, due to the iband flux contamination of a foreground galaxy. The forced photometry can overestimate the i-band flux of an i-band dropout object superposed on a foreground source, because the aperture is defined by the object image in a redder band.
In order to simulate this effect, we randomly selected 10,000 points from the HSC-SSP random source catalog in the way that we described in Section 2.3 and measured the i-band flux in an aperture placed at each point. The aperture size was set to twice the seeing FWHM at each position. The probability density distribution (PDF) of the measured fluxes is presented in Figure 6. The distribution around f ν =0 follows a Gaussian distribution, which represents the sky background fluctuation. In addition, the measured distribution has a tail toward higher f ν , which can be approximated by the function 30 = (corresponding to i AB =22.0 mag, above which the measured PDF contains less than 0.5% of the total probability). This tail contains 12% of the total probability, which is the fraction of sources affected by the foreground flux contamination. We use this function n ( ) f f fgd in the following simulations.
The foreground flux contamination is much less significant in the z-and y-bands, in which high-z quasars (meeting Equation (1)) are clearly detected and the hscPipe deblender properly apportions the measured flux. Huang et al. (2018) demonstrated that the HSC flux measurement is accurate within 0.1 mag after deblending for the vast majority of the sources.

Total Completeness
The total completeness of our selection was estimated with quasar models, created from 319 SDSS spectra of luminous (−27M i −30) quasars at z;3. This SDSS sample contains 29 radio-selected quasars, which are not sensitive to incompleteness in the color selection (see, e.g., Worseck & Prochaska 2011). We selected a sample of 29 non-radioselected quasars (i.e., objects selected for SDSS spectroscopy with other targeting criteria) from the remaining 290 objects, matched in luminosity to the radio-selected quasars, and compared the composite spectra of the two samples. This is shown in Figure 7. The composite spectra are almost identical to each other, indicating that the colors of radio-and colorselected quasars are similar and that we introduce no significant bias by using the spectra of all 319 quasars in the simulations that follow. We note that the above radio-selected quasars are still a part of the magnitude-limited SDSS sample and are biased against optically faint populations such as obscured quasars. The present estimate does not include incompleteness due to such quasars that are missing from the SDSS spectroscopic sample.
Each of the above 319 spectra was redshifted to z=5.6-6.6, with Δz=0.01 steps, with appropriate correction for the different amounts of IGM H I absorption between z∼3 and z∼6. The IGM absorption in the original SDSS spectra was removed using the mean IGM effective optical depth (τ eff ) at z3 presented by Songaila (2004). We then added IGM absorption to the redshifted model spectra by assuming the mean and scatter of τ eff taken from Eilers et al. (2018). The absorption started at a wavelength corresponding to 1 proper Mpc from the quasar, to model the effect of quasar proximity zones. The assumed proximity radius is appropriate for the mean luminosity of the SHELLQs quasars (M 1450 ∼−23 mag; Eilers et al. 2017). The damping wing of the IGM absorption was modeled following the prescription in Totani et al. (2006).  At this stage, we found that the mean and the scatter of restframe Lyα equivalent widths (EWs) of the model quasars were 64±16 Å (this includes the effect of IGM absorption and was measured with a subset of model quasars matched in redshift to the observed sample; the scatter was measured with the median absolute deviation), which are larger than those of the observed sample, 38±12 Å. This trend is opposite to the luminosity dependence known as the Baldwin effect and may be in part due to the redshift dependence of quasar SEDs, including a higher fraction of weak-line quasars found at higher redshifts (e.g., Bañados et al. 2016;Shen et al. 2018). We scaled the Lyα line of the model spectra, with the scaling factor chosen randomly from a Gaussian distribution of mean 0.6 and standard deviation 0.2, which roughly reproduces the observed EW distribution. Since the HSC bands cover only a limited portion (rest-frame wavelength 1500 Å) of the high-z quasar spectra redward of Lyα, differences in other emission lines or continuum slopes between the z∼3 SDSS quasars and the SHELLQs quasars would not be very relevant here.
The simulations of our quasar selection were performed with five million points selected from the HSC-SSP random source catalog, using the pixel flag conditions in Equation (1). We randomly assigned one of the above quasar models to each random point and calculated apparent magnitudes, assuming an absolute magnitude drawn from a uniform distribution from M 1450 =−20 to −28 mag. We then added simulated errors to the apparent magnitudes, assuming a Gaussian error distribution with standard deviation (σ) equal to the sky background rms, computed from the 5σ limiting magnitudes of the corresponding patches ( s m ; AB 5 see above). We simulated the foreground flux contamination using the PDF n ( ) f f fgd , derived in Section 3.3.
We then applied additional flux scatter with a Gaussian distribution with standard deviation 0.3 mag in each of the three bands. This was necessary to match the color distributions of the model and observed quasars, while it does not change the derived LF significantly. This additional scatter may account for sources of flux fluctuation other than those explicitly considered above, including photometry errors due to cosmic rays, image artifacts, and imperfect source deblending, the host galaxy contribution, and difference in the intrinsic SED shapes between the above SDSS quasars and the SHELLQs quasars (see, e.g., Niida et al. 2016). The resultant color distributions of the model and observed quasars are presented in Figure 8.
We selected a portion of the above simulated quasars, such that a quasar with simulated magnitudes (z AB , y AB ) on a patch p has a probability´( ps AB of being selected. This accounts for the field variance of the detection completeness. We further selected those meeting the following conditions: Finally, we calculated Bayesian quasar probabilities (P Q B ) for the selected sources, using the method described in Matsuoka et al. (2016), and counted the number of sources with  Figure 9 presents the total completeness derived above. The selection of the present complete sample is most sensitive to 5.9<z<6.5 and M 1450 <−22.5 mag. The completeness drops at z5.9 due to the color cut of i−z>2.0, while it drops more gradually at z6.5 due to the increasing contamination of brown dwarfs (which reduces the quasar probability P Q B ). The figure also shows that several quasars located in the high completeness region are not included in the complete sample. This is caused by various reasons; some quasars are in survey fields that fail to meet the pixel flag conditions (Equation (1)) in the S17A data release, and some quasars have i−z colors just below the threshold of 2.0. The faintest quasars with M 1450 >−22.5 mag simply fail to meet the condition < z z AB AB slim . In the following section, we use the completeness functions of SDSS, CFHQS, and SHELLQs to derive a single LF. These functions were all derived with quasar models tied to spectra of SDSS quasars at z∼3, while the IGM absorption models in SDSS and CFHQS were created from τ eff data older than those we used here for the SHELLQs sample. We tested another IGM absorption model for the SHELLQs sample, with the mean and scatter of the τ eff determined empirically to reproduce the data in Songaila (2004), and found little change in the derived completeness or LF. In addition, while the completeness correction is most important at the faintest luminosity of a given sample, the faintest SDSS/CFHQS quasars have smaller available volumes (V a ; see below and Table 2) and thus smaller weights in LF calculation than do the CFHQS/SHELLQs quasars with similar luminosities and high completeness. Thus, we conclude that no significant bias is introduced by combining the completeness functions of the three surveys.

Luminosity Function
First, we derive the binned LF using the 1/V a method (Avni & Bahcall 1980). The cosmic volume available to discover a quasar, in a magnitude bin ΔM 1450 , is given by where the sum is taken over the quasars in the magnitude bin. This expression ignores the redshift evolution of the LF over the measured range (5.7z6.5); we will take this evolution into account in the parametric LF described below.
Here we combine the three complete samples of quasars from SDSS, CFHQS, and SHELLQs to derive a single binned LF over −22<M 1450 <−30 mag (we use the completeness functions and the survey areas of SDSS and CFHQS described in Sections 2.1 and 2.2). We set the bin size ΔM 1450 = 0.5 mag, except at both ends of the luminosity coverage where the sample size is small. The results of this calculation are listed in Table 4 and presented in Figure 10. The derived LF agrees well with the previous results from SDSS (Jiang et al. 2016) and CFHQS (Willott et al. 2010) at M 1450 <−25 mag, and significantly improves the accuracy at fainter magnitudes. It may be worth mentioning that the number density of the brightest bin measured by Jiang et al. (2016) and in this work do not exactly match, although the two works use a single SDSS quasar in common. This is due to the different choice of bin center and width, which is known to have a significant impact on the binned LF when the sample    size is small. On the other hand, we significantly increased the available survey volume for the faintest bin at M 1450 =−22.00 and found a number density lower than (but consistent within 1σ) the previous measurement by Willott et al. (2010).
Next, we derive the parametric LF, using a commonly used double power-law function: where α and β are the faint-and bright-end slopes, respectively. We fix the redshift evolution term to k=−0.47 (Willott et al. 2010) or k=−0.7 (Jiang et al. 2016); we found that the choice makes little difference in the determination of other parameters (see below). Following the argument in Jiang et al. (2016), we adopt k=−0.7 as our standard value. The parameters * M 1450 and * F give the break magnitude and normalization of the LF, respectively.
We perform a maximum likelihood fit (Marshall et al. 1983) to determine the four free parameters (α, β, * M 1450 , and Φ * ). Specifically, we maximize the likelihood L by minimizing = -S L 2 ln , given by where the sum in the first term is taken over all quasars in the sample. The resultant parametric LF is presented in Figure 10, and the best-fit LF parameters are listed in the first row of Table 5. Figure  We also performed LF calculations with k fixed to −0.47 or allowed to vary as a free parameter, and found that the other LF parameters are not very sensitive to the choice of k. These results are listed in the second and third rows of Table 5. The fitting with the variable k favors relatively flat LF evolution ( = --+ k 0.2 0.1 0.2 ), which may be consistent with the tendency for k to be smaller for lower-luminosity quasars seen in Jiang et al. (2016, their Figure 10). But, given the short redshift baseline of the present sample, we chose to adopt the fixed value k=−0.7 for our standard LF.
Recently, Kulkarni et al. (2018) reported a very bright break magnitude of * = --+ M 29.2 1450 1.9 1.1 mag at z∼6 by reanalyzing the quasar sample constructed by Jiang et al. (2016), Willott et al. (2010), and Kashikawa et al. (2015). However, their data favor a single power-law LF, and thus the break magnitude was forced to be at the bright end of the sample in their LF fitting (Kulkarni et al. 2018). The present work indicates that the LF breaks at a much fainter magnitude, in the luminosity range that has been poorly explored previously.
It may be worth noting that the CFHQS-deep survey discovered one quasar in the M 1450 =−22.00 bin from Figure 11. Confidence regions (light gray: 1σ, gray: 2σ, dark gray: 3σ) of the individual LF parameters. The best-fit values are marked by the crosses. Table 5 Parametric Luminosity Function  Table 2). This is presumably due to statistical fluctuations. Based on the present parametric LF, the expected total number of quasars in the CFHQS-deep survey is roughly one, with the most likely luminosity in the range −25M 1450 −22 mag. In reality, the survey discovered one quasar with M 1450 =−21.5 mag and none at brighter magnitudes, which is consistent with the expectation. On the other hand, the expected number of SHELLQs quasars in the M 1450 =−22.00 bin is roughly one. This is consistent with the actual discovery of no quasars in this bin, given Poisson noise. The SHELLQs complete sample used here includes five objects with narrow Lyα lines (FWHM<500 km s −1 ) at −23.5<M 1450 <−22.5. We classified them as quasars based on their extremely high Lyα luminosities, featureless continuum, and possible mini broad absorption-line system of N V λ1240 seen in their composite spectrum (Matsuoka et al. 2018b). It is possible that they are not in fact type 1 quasars, so for reference, we recalculated the binned LF at −23.5< M 1450 <−22.5, omitting these five objects, and listed the results in the last two rows of Table 4. The parametric LF in this case is reported in the fourth row of Table 5, which shows a modest difference from the standard case.
We also calculated the LF by limiting the sample to the 89 quasars in our complete sample at z>5.9, the redshift range over which CFHQS and SHELLQs are most sensitive (see Figure 1). The resultant parametric LF is listed in the last row of Table 5. The LF in this case has a slightly brighter * M 1450 and steeper α than the standard LF, but the difference is smaller than the fitting uncertainty.  Kashikawa et al. 2015), who had only a few low-luminosity quasars in their samples. The extrapolation of our LF underpredicts the number densities of faint AGNs compared to those reported by Giallongo et al. (2015), while the former is consistent with the more recent measurements by Parsa et al. (2018). On the other hand, we note that the above X-ray measurements are immune to dust obscuration, and that the discrepancy with the rest-UV measurements, if any, could be due to the presence of a large population of obscured AGNs in the high-z universe. Finally, Figure 12 indicates that LBGs (taken from Ono et al. 2018) outnumber quasars at M 1450 >−23 mag. This is consistent with our experience from the SHELLQs survey, which found increasing numbers of LBGs contaminating the quasar candidate sample at z AB >23 mag (Matsuoka et al. 2016(Matsuoka et al. , 2018a(Matsuoka et al. , 2018b.
We compare the present LF with that recently derived at z∼4 (Akiyama et al. 2018) andz∼5 (McGreer et al. 2018) in Figure 13. The overall shape of the binned LF remains relatively similar, while there is a steep decline of the total number density toward higher redshifts. However, the best-fit break magnitudes reported in the above studies differ substantially, i.e.,  (2018), which is significantly steeper than measured at z∼4 (β∼−3.1; Akiyama et al. 2018) or at z∼6 (β∼−2.7; this work). As shown in the middle panel of Figure 11, the brightend slope and the break magnitude are strongly covariant in the parametric LF fitting. We found that the binned LF of McGreer et al. (2018) can also be fitted reasonably well with β=−3.0, as shown in Figure 13    to the break magnitudes at z∼4 and z∼6. The figure also displays the parametric LFs reported by Kulkarni et al. (2018); while these LFs match the data in the luminosity ranges covered by their sample, the LFs seem to overpredict the number densities of fainter quasars presented in the recent studies by Akiyama et al. (2018), McGreer et al. (2018, and this paper. Since the LF is a product of the mass function and the Eddington ratio function of SMBHs, it is not straightforward to interpret the significant flattening observed at M 1450 −25 mag in terms of a unique physical model. It could indicate relatively mass-independent number densities and/or quasar radiation efficiency at low SMBH masses. We will compare our LF with theoretical models in a forthcoming paper. Alternatively, as discussed above, the LF flattening may indicate an increasing fraction of obscured AGNs toward low luminosities, especially in light of the X-ray results in Figure 12. This could be an interesting subject for future deep X-ray observations, such as those that ATHENA (Nandra et al. 2013) will achieve.

Contribution to Cosmic Reionization
There is much debate about the source of photons that are responsible for cosmic reionization, as we discussed in Section 1. Here we derive the total ionizing photon density from quasars per unit time,ṅ ion (s −1 Mpc −3 ), and compare with that necessary to keep the IGM fully ionized. The ionizing photon density can be calculated as where f esc is the photon escape fraction, ò 1450 (erg s −1 Hz −1 Mpc −3 ) is the total photon energy density from quasars at 1450 Å,  (12), we used a broken power-law quasar SED ( f ν ∝ ν −1.70 at λ<912 Å and ∝ν −0.61 at λ>912 Å) presented by Lusso et al. (2015) and integrated from the H I Lyman limit (frequency ν=ν LL ) to the He II Lyman limit (ν=4ν LL ). The implicit assumptions here are that the above SED, created from luminous quasars at z∼2.4, holds for the present high-z quasars, and that all ionizing photons with ν<4ν LL are absorbed by the IGM. The resultant photon density is =  n 10 ion 48.8 0.1 s −1 Mpc −3 at z=6.0 for f esc =1. We would get lowerṅ ion for f esc <1, which may be the case for lowluminosity quasars (Cristiani et al. 2016;Micheva et al. 2017;Grazian et al. 2018). The energy density at 912 Å is estimated to be  =  10 912 22.9 0.1 erg s −1 Hz −1 Mpc −3 , which is close to the value reported by Haardt & Madau (2012) at z=6. The results presented in this section change very little when the faint limit of the integral in Equation (11) is changed to M 1450 =−10 mag, or when the five SHELLQs quasars with narrow Lyα (see Section 4) are excluded.
On the other hand, the evolution of the H II volume-filling factor in the IGM, ( ) wheren H andt rec are the mean hydrogen density and recombination time, respectively (Madau et al. 1999 where C H II represents an effective H II clumping factor (Bolton & Haehnelt 2007). The ionizing photon density we found above, given our LF, is less than 10% ofṅ ion crit for the plausible range of = -C 1.0 5.0 H II (Shull et al. 2012). This means that quasars alone cannot sustain reionization. For reference, we would get = n 10 ion 50.3 s −1 Mpc −3~ṅ ion crit if we assumed no LF break (α=β=−2.73) and integrated Equation (11) from M 1450 =−18 to −30 mag.
Finally, we numerically integrate Equation (13) and track the evolution of Q H II driven solely by quasar radiation. We assume that the IGM was neutral at z=15, and thatṅ ion was constant in time (i.e., it stayed at 10 48.8±0.1 s −1 Mpc −3 ) or evolved as µ -10 z 0.7 (i.e., proportional to the LF normalization found around z = 6) at 5<z<15. We followed Robertson et al. (2015) to estimaten H andt rec . The results of this calculation are presented in Figure 14. For reference, we also plot the Q H II evolution driven by star-forming galaxies, using the star formation rate density at < z 15 presented in Robertson et al. (2015). This figure demonstrates that star-forming galaxies can supply enough high-energy photons to ionize the IGM by z=6, while quasars cannot. We thus conclude that quasars are not a major contributor to reionization. Even if there is a large population of obscured AGNs that are missed by rest-UV surveys (see the discussion in Section 4), they are unlikely to release many ionizing photons, since the ionizing photon escape fraction from these objects would be close tof 0 esc .

Summary
This paper presented new measurements of the quasar LF at z 6, which is now established over an unprecedentedly wide magnitude range from = -M 30 1450 to −22 mag. We collected a complete sample of 110 quasars from the SDSS, the CFHQS, and the SHELLQs surveys. The completeness of the SHELLQs quasar selection was carefully evaluated, and we showed that the selection is most sensitive to quasars with < < z 5.9 6.5 and < - mag. This accounts for <10% of the critical rate necessary to keep the IGM fully ionized at z=6.0. We conclude that quasars are not a major contributor to cosmic reionization.
The HSC-SSP survey is making steady progress toward its goal of observing 1400 deg 2 in the Wide layer. We will continue follow-up spectroscopy to construct a larger complete sample ofz 6 quasars, down to luminosities lower than probed in the present work. We are also starting an intensive effort to explore higher redshifts, with the aim of establishing the quasar LF atz 7. At the same time, we are collecting near-IR spectra to measure the SMBH masses and mass accretion rates, which will be used in combination with the LFs to understand the growth of SMBHs in the early universe. The ALMA follow-up observations are also ongoing, which will provide valuable information on the formation and evolution of the host galaxies.