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

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₁₄₅₀ = −30 to −24 mag. The absolute magnitudes (M₁₄₅₀) 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², 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 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)

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₁₄₅₀ = −27 to −24 mag. The remaining quasar, with M₁₄₅₀ = −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², 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₁₄₅₀) 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₁₄₅₀ by extrapolating the continuum spectrum redward of Lyα, assuming a power-law shape f_λ ∝ λ^−1.5. The resultant M₁₄₅₀ 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.

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 15 diameter, covered by 116 2K × 4K fully depleted Hamamatsu CCDs, with a pixel scale of 017. 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² 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 20 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²) down to deeper limiting magnitudes (r_AB = 27.1 and 27.7 mag, respectively). Data reduction was performed with the dedicated pipeline hscPipe (Bosch et al. 2018). We use the point-spread function (PSF) magnitude (${m}_{\mathrm{PSF},\mathrm{AB}}$, or simply m_AB) and the CModel magnitude (${m}_{\mathrm{CModel},\mathrm{AB}}$), which are measured by fitting the PSF models and two-component, PSF-convolved 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:

$$\begin{eqnarray}\begin{array}{r}{z}_{\mathrm{AB}}\lt 24.5\ \&\ {\sigma }_{z}\lt 0.155\ \&\ {i}_{\mathrm{AB}}-{z}_{\mathrm{AB}}\gt 2.0\\ \&\ {z}_{\mathrm{AB}}-{z}_{\mathrm{CModel},\mathrm{AB}}\lt 0.15\\ \&\ {\mathtt{merge}}.{\mathtt{peak}}.({\mathtt{g}},{\mathtt{r}},{\mathtt{z}},{\mathtt{y}})=({\mathtt{f}},{\mathtt{f}},{\mathtt{t}},{\mathtt{t}})\\ \&\ ({\mathtt{z}},{\mathtt{y}}).{\mathtt{inputcount}}.{\mathtt{value}}\geqslant (2,2)\\ \&\ ({\mathtt{i}},{\mathtt{z}},{\mathtt{y}}).{\mathtt{pixelflags}}.{\mathtt{edge}}=({\mathtt{f}},{\mathtt{f}},{\mathtt{f}})\\ \&\ ({\mathtt{i}},{\mathtt{z}},{\mathtt{y}}).{\mathtt{pixelflags}}.{\mathtt{saturatedcenter}}\\ =\,({\mathtt{f}},{\mathtt{f}},{\mathtt{f}})\\ \&\ ({\mathtt{i}},{\mathtt{z}},{\mathtt{y}}).{\mathtt{pixelflags}}.{\mathtt{crcenter}}=({\mathtt{f}},{\mathtt{f}},{\mathtt{f}})\\ \&\ ({\mathtt{i}},{\mathtt{z}},{\mathtt{y}}).{\mathtt{pixelflags}}.{\mathtt{bad}}=({\mathtt{f}},{\mathtt{f}},{\mathtt{f}})\\ \&\ ({\mathtt{i}},{\mathtt{z}},{\mathtt{y}}).{\mathtt{pixelflags}}.{\mathtt{bright}}.{\mathtt{objectcenter}}\\ =\,({\mathtt{f}},{\mathtt{f}},{\mathtt{f}})\end{array}\end{eqnarray} \tag{ 1 }$$

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 z- and 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 10, 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}_{{\rm{Q}}}^{{\rm{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 ${P}_{{\rm{Q}}}^{{\rm{B}}}\gt 0.1$ 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}_{\mathrm{AB}}^{\mathrm{splim}}\simeq 24.0$ mag. The actual ${z}_{\mathrm{AB}}^{\mathrm{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}_{\mathrm{AB}}\leqslant {z}_{\mathrm{AB}}^{\mathrm{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}_{{\rm{AB}}}^{{\rm{splim}}}$, outside the above redshift range, or fail to meet one or more of the criteria listed in Equation (1). The absolute magnitudes (M₁₄₅₀) were estimated in the same way as used for the SDSS quasars (see above).

Table 3.  SHELLQs Survey Fields

Name	R.A. Range	Decl. Range	Area	${z}_{\mathrm{AB}}^{\mathrm{splim}}$	${N}_{\mathrm{obj}}$
	(deg)	(deg)	(deg2)	(mag)
XMM	28–41	−7–+3	83.7	24.1	5
GAMA09H	127–155	−3–+6	165.1	23.8	8
WIDE12H	173–200	−3–+3	106.5	23.8	10
GAMA15H	205–227	−3–+3	100.7	24.0	8
VVDS	330–357	−2–+7	124.7	24.2	13
HECTOMAP	220–252	+42–+45	65.4	24.0	4
Total	⋯	⋯	646.1	⋯	48

Note. The field names refer to the distinct areas covered in the HSC-SSP survey to date; see Aihara et al. (2018a) for details. ${z}_{\mathrm{AB}}^{\mathrm{splim}}$ and ${N}_{\mathrm{obj}}$ represent the spectroscopic limiting magnitude and the number of quasars included in the present complete sample, respectively.

Download table as:  ASCIITypeset image

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⁻², 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.

Source detection in the HSC data processing pipeline (hscPipe; Bosch et al. 2018) is performed on PSF-convolved images by finding pixels with flux >5σ above the background sky. Here, σ is the root-mean-square (rms) of the local background fluctuations. For a point source, this thresholding is approximately equivalent to ${m}_{\mathrm{AB}}\lt {m}_{\mathrm{AB}}^{5\sigma }$, where ${m}_{\mathrm{AB}}^{5\sigma }$ represents the PSF limiting magnitude at which S/N = 5 (see Bosch et al. 2018 for a description of the theory). The HSC database stores ${m}_{\mathrm{AB}}^{5\sigma }$ measurements for each patch (12' × 12') in the survey. As shown in Figure 2, all but a small fraction of the survey patches have ${z}_{\mathrm{AB}}^{5\sigma }\gt 24$ 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}_{\mathrm{AB}}^{\mathrm{splim}}=23.8\mbox{--}24.2$ mag.

Zoom In Zoom Out Reset image size

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.²⁹ 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):

$$\begin{eqnarray}&&f({m}_{\mathrm{AB}})=\displaystyle \frac{{f}_{\max }\,-\,{f}_{\min }}{2}(\tanh [\alpha ({m}_{\mathrm{AB}}^{50}-{m}_{\mathrm{AB}})]+1)+{f}_{\min },\end{eqnarray} \tag{ 2 }$$

where f_max, f_min, α, and ${m}_{\mathrm{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}_{\mathrm{AB}}^{50}$ values with the 5σ limiting magnitudes (${m}_{\mathrm{AB}}^{5\sigma }$) described above. These two quantities agree very well with each other, as expected given that the hscPipe detection threshold is approximately equivalent to ${m}_{\mathrm{AB}}\lt {m}_{\mathrm{AB}}^{5\sigma }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ${f}_{\det }({z}_{\mathrm{AB}},p)$ and ${f}_{\det }({y}_{\mathrm{AB}},p)$ were defined using Equation (2). We retrieved ${z}_{\mathrm{AB}}^{5\sigma }$ and ${y}_{\mathrm{AB}}^{5\sigma }$ from the survey database and used them as surrogates for ${z}_{\mathrm{AB}}^{50}$ and ${y}_{\mathrm{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.

The SHELLQs algorithm uses the criterion

$$\begin{eqnarray}&&{z}_{\mathrm{AB}}-{z}_{\mathrm{CModel},\mathrm{AB}}\lt 0.15\end{eqnarray} \tag{ 3 }$$

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}_{\mathrm{AB}}^{50}$) = (1.00, 0.72, 0.76, 24.5). This best-fit function, ${f}_{\mathrm{ps}}({z}_{\mathrm{AB}})$, is used to simulate the selection completeness of point sources in the following.

Zoom In Zoom Out Reset image size

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, 2015): (i) the typical host galaxy luminosity ranges from M_UV = −18 to −21 mag (corresponding to ${z}_{\mathrm{CModel},\mathrm{AB}}\sim 25.5\mbox{--}28.5$ 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 (Harikane et al. 2018; Ono et al. 2018). We used 231 LBGs with $24.0\lt {z}_{\mathrm{CModel},\mathrm{AB}}\lt 25.0$ mag, where AGN contamination to the sample is small (Ono et al. 2018). For each LBG, we randomly assigned M_UV from −18 to −21 mag, assumed a flat UV spectral slope (β = −2.0; Stanway et al. 2005), and calculated the corresponding CModel flux (${f}_{\mathrm{CModel}}^{\mathrm{sim}}$) at z = 6. The PSF flux was calculated as ${f}_{\mathrm{PSF}}^{\mathrm{sim}}={f}_{\mathrm{CModel}}^{\mathrm{sim}}\times ({f}_{\mathrm{PSF}}^{\mathrm{obs}}/{f}_{\mathrm{CModel}}^{\mathrm{obs}})$, where ${f}_{\mathrm{PSF}}^{\mathrm{obs}}/{f}_{\mathrm{CModel}}^{\mathrm{obs}}$ is the ratio between the PSF and CModel fluxes observed for the individual LBGs. Then we added various AGN fluxes (${f}^{\mathrm{AGN}}={f}_{\mathrm{PSF}}^{\mathrm{AGN}}={f}_{\mathrm{CModel}}^{\mathrm{AGN}}$) artificially and calculated the fraction of the simulated objects that satisfy Equation (3) and are thus "unresolved:"

$$\begin{eqnarray}&&-2.5\mathrm{log}\left(\displaystyle \frac{{f}_{\mathrm{PSF}}^{\mathrm{sim}}+{f}^{\mathrm{AGN}}}{{f}_{\mathrm{CModel}}^{\mathrm{sim}}+{f}^{\mathrm{AGN}}}\right)\lt 0.15.\end{eqnarray} \tag{ 4 }$$

We found that the unresolved fraction is 100% at AGN magnitudes z_AB < 25.0 mag and decreases to 90% at 26.0 mag. 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}_{\mathrm{AB}}-{z}_{\mathrm{CModel},\mathrm{AB}}\lt 0.15$ 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.

As we wrote previously, we failed to recover one of the eight previously known quasars in our survey footprint, due to the i-band 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³⁰ ${f}_{\mathrm{fgd}}({f}_{\nu })\,=3.3\ {e}^{-5\sqrt{{f}_{\nu ,29}}}+0.0014$ (${f}_{\nu ,29}={f}_{\nu }\times {10}^{29}$ erg s⁻¹ cm⁻² Hz⁻¹) truncated at ${f}_{\nu ,29}=5.8$ (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 ${f}_{\mathrm{fgd}}({f}_{\nu })$ in the following simulations.

Zoom In Zoom Out Reset image size

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.

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-radio-selected 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 color-selected 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.

Zoom In Zoom Out Reset image size

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₁₄₅₀ ∼ −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 rest-frame 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₁₄₅₀ = −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 (${m}_{\mathrm{AB}}^{5\sigma };$ see above). We simulated the foreground flux contamination using the PDF ${f}_{\mathrm{fgd}}({f}_{\nu })$, 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.

Zoom In Zoom Out Reset image size

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 ${f}_{\det }({z}_{\mathrm{AB}},p)\times {f}_{\det }({y}_{\mathrm{AB}},p)\times {f}_{\mathrm{ps}}({z}_{\mathrm{AB}})$ of being selected. This accounts for the field variance of the detection completeness. We further selected those meeting the following conditions:

$$\begin{eqnarray}&&{z}_{\mathrm{AB}}\lt {z}_{\mathrm{AB}}^{\mathrm{slim}}\ \mathrm{and}\ {\sigma }_{z}\lt 0.155\ \mathrm{and}\ {i}_{\mathrm{AB}}-{z}_{\mathrm{AB}}\gt 2.0.\end{eqnarray} \tag{ 5 }$$

Finally, we calculated Bayesian quasar probabilities (${P}_{{\rm{Q}}}^{{\rm{B}}}$) for the selected sources, using the method described in Matsuoka et al. (2016), and counted the number of sources with ${P}_{{\rm{Q}}}^{{\rm{B}}}\gt 0.1$. The total completeness, ${f}_{\mathrm{comp}}(z,{M}_{1450})$, is given by the ratio between the output and input numbers of random sources, calculated in bins of z and M₁₄₅₀. There are roughly 400 simulated quasars in each bin with sizes Δz = 0.01 and ${\rm{\Delta }}{M}_{1450}=0.05$.

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₁₄₅₀ < −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}_{{\rm{Q}}}^{{\rm{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₁₄₅₀ > −22.5 mag simply fail to meet the condition ${z}_{\mathrm{AB}}\lt {z}_{\mathrm{AB}}^{\mathrm{slim}}$.

Zoom In Zoom Out Reset image size

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.