Subaru High-$z$ Exploration of Low-Luminosity Quasars (SHELLQs). X. Discovery of 35 Quasars and Luminous Galaxies at $5.7\leq$ $z$ $\leq 7.0$

Quasars in the high-$z$ universe¹¹1 Throughout this paper, “high-$z$” denotes $z>5.7$, where quasars are observed as $i$-band dropouts in the Sloan Digital Sky Survey (SDSS) filter system (Fukugita et al., 1996). The term “$X$-band dropouts” or “$X$-dropouts” refers to objects which are much fainter (and often invisible) in the $X$ and bluer bands than in the redder bands. have been used as an unique probe into early cosmic history. The progress of cosmic reionization can be measured from H I absorption imprinted on the rest-frame ultraviolet (UV) spectrum of a high-$z$ quasar; this absorption is very sensitive to the neutral fraction of the foreground intergalactic medium (IGM; Gunn & Peterson, 1965; Fan et al., 2006a). The luminosity and mass functions of quasars reflect the seeding and growth mechanisms of supermassive black holes (SMBHs), which can be studied through comparison with theoretical models (e.g., Volonteri, 2012; Ferrara et al., 2014; Madau et al., 2014). Measurements of quasar host galaxies and surrounding environments tell us about the earliest mass assembly of galaxies, which may happen in the highest-density peaks of the underlying dark matter distribution (e.g., Goto et al., 2009; Decarli et al., 2017).

For the past several years, we have been carrying out a project to search for high-$z$ quasars, named “Subaru High-$z$ Exploration of Low-Luminosity Quasars (SHELLQs)”, with the Subaru 8.2-m telescope. Matsuoka et al. (2016, 2018a, 2018b, 2019) have already reported the discovery of 65 high-$z$ quasars, along with 24 high-$z$ luminous galaxies, 6 [O III] emitters at $z\sim 0.8$, and 43 Galactic cool dwarfs (low-mass stars and brown dwarfs). These objects include many low-luminosity quasars with redshifts up to $z=7.07$, objects that have been difficult to find in past shallower surveys. Combined with samples of luminous quasars from the SDSS (Jiang et al., 2016) and the Canada-France High-$z$ Quasar Survey (Willott et al., 2010), the SHELLQs sample has allowed us to establish the quasar luminosity function (LF) at $z=6$ over an unprecedented range of magnitude ($-22<M_{1450}<-30$ mag; Matsuoka et al., 2018c). We found that the LF has a break at $M_{1450}\sim-25$ mag and flattens significantly toward the faint end. With this LF shape, we predict that quasars are responsible for less than 10% of the ionizing photons that are necessary to keep the IGM fully ionized at $z=6$.

We are also carrying out multi-wavelength follow-up observations of the discovered quasars. Near-infrared (IR) spectra were obtained for 6 quasars, which revealed that the SHELLQs quasars have a wide range of accretion properties, from sub-Eddington accretion onto a massive ($M_{\rm BH}>10^{9}M_{\odot}$) black hole to Eddington accretion onto a less massive ($M_{\rm BH}<10^{8}M_{\odot}$) black hole (Onoue et al., 2019). Atacama Large Millimeter / submillimeter Array (ALMA) observations were conducted toward 7 quasars, and we found that the host galaxies are forming stars at a rate below that of luminous quasar hosts (e.g., Decarli et al., 2018), and that the low-luminosity quasars are more or less on the local relation between black hole mass and host mass (Izumi et al., 2018, 2019). We are also exploiting archival Wide-field Infrared Survey Explorer data to look for red quasars in the SHELLQs sample, which may represent a younger phase of evolution than normal quasars (N. Kato et al., in prep.).

This paper is the tenth in a series of SHELLQs publications, and reports spectroscopic identification of an additional 67 objects. The paper is structured as follows. We give brief descriptions of the data and methods used for candidate selection in §2. Spectroscopic follow-up observations are presented in §3. Results and discussion appear in §4. We adopt the cosmological parameters $H_{0}$ = 70 km s^-1 Mpc^-1, $\Omega_{\rm M}$ = 0.3, and $\Omega_{\rm\Lambda}$ = 0.7. All magnitudes in the optical and near-IR bands are presented in the AB system (Oke & Gunn, 1983), and are corrected for Galactic extinction (Schlegel et al., 1998). We use the point spread function (PSF) magnitude ($m_{\rm AB}$) and the CModel magnitude ($m_{\rm CModel,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). In what follows, we refer to $z$-band magnitudes with the AB subscript (“$z_{\rm AB}$”), while redshift $z$ appears without a subscript.

The SHELLQs project is based on the multi-band imaging data produced by the Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP) survey (Aihara et al., 2018). HSC is a wide-field camera on Subaru (Miyazaki et al., 2018). It has a nearly circular field of view of 1^∘.5 diameter, covered by 116 2K $\times$ 4K Hamamatsu fully depleted CCDs. The pixel scale is 0″.17. The HSC-SSP survey has three layers with different combinations of area and depth. The Wide layer is observing 1400 deg² mostly along the celestial equator, with 5$\sigma$ point-source depths of ($g_{\rm AB}$, $r_{\rm AB}$, $i_{\rm AB}$, $z_{\rm AB}$, $y_{\rm AB}$) = (26.5, 26.1, 25.9, 25.1, 24.4) mag measured in 2″.0 aperture. 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 $\sim$3 minutes each. The Deep and the UltraDeep layers are observing smaller areas (27 and 3.5 deg²) down to deeper limiting magnitudes ($r_{\rm AB}$ = 27.1 and 27.7 mag, respectively). The survey data are reduced with the dedicated pipeline hscPipe (Bosch et al., 2018) derived from the Large Synoptic Survey Telescope software pipeline (Jurić et al., 2017).

The procedure of our quasar candidate selection has been described in the previous papers (Matsuoka et al., 2016, 2018a, 2018b), so here we repeat the essence only briefly. The selection starts from the HSC-SSP survey data in the three layers, over the area observed in the $i$, $z$, and $y$ bands for $i$-dropout selection or in the $z$ and $y$ bands for $z$-dropout selection. We include data that haven’t yet reached the planned full depth. All sources meeting the following criteria, and that are not flagged as having suspect detection or photometry, are selected:

or

Equations (1) and (2) select $i$-dropout and $z$-dropout candidates, respectively. Here ($\sigma_{i}$, $\sigma_{z}$, $\sigma_{y}$) refer to the errors of the PSF magnitudes ($i_{\rm AB}$, $z_{\rm AB}$, $y_{\rm AB}$), as measured by hscPipe. The $y$-band limiting magnitude in Equation (2) is deeper than the limit used in our previous papers ($y_{\rm AB}<24.0$), as we are shifting our focus to establishing the quasar LF at $z\sim 7$. The difference between the PSF and CModel magnitudes ($m_{\rm AB}-m_{\rm CModel,AB}$) is used to exclude extended sources. Sources with more than $3\sigma$ detection in the $g$ or $r$ band (if these bands are available) are removed as likely low-$z$ interlopers. We also remove those sources whose images appear to be moving objects, transients, or artifacts, spotted by automatic or eye inspection (see Matsuoka et al., 2016). We always keep ambiguous cases in our sample for spectroscopic follow-up, so that we do not discard any real quasars in the selection.

The candidates selected from the HSC data are matched, within 1″.0 separation, to the public near-IR catalogs from the United Kingdom Infrared Telescope Infrared Deep Sky Survey (UKIDSS; Lawrence et al., 2007) data release (DR) 10 and 11, the Visible and Infrared Survey Telescope for Astronomy (VISTA) Kilo-degree Infrared Galaxy survey (VIKING) DR 4 and 5, and the VISTA Deep Extragalactic Observations Survey (Jarvis et al., 2013) DR 5. In practice, only a small fraction of the final quasar candidates have near-IR photometry, due to limited sky coverage and/or depth of the above surveys. Using all available magnitudes in the $i$, $z$, $y$, $J$, $H$, and $K$ bands, we calculate the Bayesian quasar probability ($P_{\rm Q}^{\rm B}$) for each candidate. The calculation is based on models of spectral energy distributions and surface densities of high-$z$ quasars and contaminating brown dwarfs, as a function of magnitude; galaxy models are not included in the algorithm at present (see Matsuoka et al., 2016). We keep those sources with $P_{\rm Q}^{\rm B}>0.1$ in the sample of candidates, while removing sources with lower quasar probabilities.

We used the latest HSC-SSP (internal) DR, which includes observations carried out between 2014 March and 2018 January. The data cover roughly 900 deg², when we limit to the fields where at least (1, 2, 2) exposures were taken in the ($i$, $z$, $y$) bands for $i$-dropout selection, or at least (1, 2) exposures were taken in the ($z$, $y$) bands for $z$-dropout selection. We found $\sim$500 final candidates over this area, and put highest priority on (i) all the $z$-dropouts and (ii) the $i$-band dropouts with $i_{\rm AB}-z_{\rm AB}>2.0$, $z_{\rm AB}<$ 24.0 mag, and $y$-band detection. We have almost completed follow-up spectroscopy of these $\sim$200 high-priority candidates, as reported in the past and present papers, and continue to observe the remaining low-priority candidates.

We carried out follow-up spectroscopy of 67 candidates from 2018 March through 2019 July, using the Optical System for Imaging and low-intermediate-Resolution Integrated Spectroscopy (OSIRIS; Cepa et al., 2000) mounted on the 10.4-m Gran Telescopio Canarias (GTC), and the Faint Object Camera and Spectrograph (FOCAS; Kashikawa et al., 2002) mounted on Subaru. We prioritized observations in such a way that the targets with brighter magnitudes and higher $P_{\rm Q}^{\rm B}$ were observed at earlier opportunities. Table 1 is a journal of the spectroscopic observations.

GTC is a 10.4-m telescope located at the Observatorio del Roque de los Muchachos in La Palma, Spain. We used OSIRIS with the R2500I grism and 1″.0-wide longslit, which provides spectral coverage from $\lambda_{\rm obs}$ = 0.74 to 1.0 $\mu$m with a resolution $R\sim 1500$. The observations (Program IDs: GTC3-18A, GTC8-18B, and GTC32-19A) were carried out in queue mode on dark and gray nights, under spectroscopic sky conditions and with seeing 0″.6 – 1″.3.

The Subaru observations were carried out as part of two Subaru Intensive Programs (Program IDs: S16B-071I and S18B-011I). We used FOCAS in the multi-object spectrograph (MOS) mode with the VPH900 grism and SO58 order-sorting filter. The widths of the slitlets were set to 1″.0. This configuration provides spectral coverage from $\lambda_{\rm obs}$ = 0.75 to 1.05 $\mu$m with a resolution $R\sim 1200$. All the observations were carried out on gray nights, which were occasionally affected by cirrus, with seeing 0″.4 – 1″.0.

All the data obtained with GTC and Subaru were reduced using the Image Reduction and Analysis Facility (IRAF). Bias correction, flat fielding with dome flats, sky subtraction, and 1d extraction were performed in the standard way. The wavelength was calibrated with reference to sky emission lines. The flux calibration was tied to white dwarfs (Feige 34, Feige 110, G191-B2B, GD 153, Ross 640) or a B-type standard star (HILT 600), observed as standard stars within a few days of the target observations. We corrected for slit losses by scaling the spectra to match the HSC magnitudes in the $z$ and $y$ bands for the $i$- and $z$-band dropouts, respectively.

We present the reduced spectra in Figures 1 – 7. Based on these spectra, we identified 28 high-$z$ quasars, 7 high-$z$ galaxies, 10 strong [O III] emitters at $z\sim 0.8$, and 22 cool dwarfs, as detailed below. The photometric properties of the observed candidates are listed in Table 1. Table 2 lists the objects detected in the near-IR bands.

Figures 1 – 3 present 28 new quasars we identified at $5.7\leq z\leq 7.0$. Their spectroscopic properties are presented in the first section of Table 3. The seven quasars at $z>6.6$ are observed as $z$-dropouts on the HSC images, while the other objects are $i$-dropouts. The quasars have broad Ly$\alpha$ lines, blue rest-UV continua, and/or sharp continuum breaks just blueward of Ly$\alpha$, the properties which are characteristic of high-$z$ quasars. Following the previous papers, we classified the objects with very luminous ($L>10^{43}$ erg s^-1) and narrow (full width at half maximum [FWHM] $<$ 500 km s^-1) Ly$\alpha$ emission as possible quasars (see the discussion below). These objects tend to be found at the faintest magnitudes of our selection (Matsuoka et al., 2018b), and have been missed in past shallower surveys.

Redshifts of the discovered quasars were determined from the Ly$\alpha$ lines, assuming that the observed line peaks correspond to the intrinsic Ly$\alpha$ wavelength (1216 Å in the rest frame). This assumption is not always correct, due to the strong H I absorption by the neutral IGM. More accurate redshifts require observations of other emission lines, such as Mg II $\lambda$2800 observed in the near-IR or [C II] 158 $\mu$m accessible with ALMA. When there is no clear Ly$\alpha$ lines, we obtained rough estimates of redshift from the wavelengths of the onset of the Gunn & Peterson (1965) trough. Therefore the redshifts presented here are only approximate, with the uncertainties up to ${\Delta}z\sim 0.1$.

Absolute magnitudes ($M_{1450}$) and Ly$\alpha$ line properties of the quasars were measured as follows. For every object, we defined a continuum window at wavelengths relatively free from strong sky emission lines, and extrapolated the measured continuum flux to estimate $M_{1450}$. A power-law continuum model with a slope $\alpha=-1.5$ ($F_{\lambda}\propto\lambda^{\alpha}$; e.g., Vanden Berk et al. (2001)) was assumed. Since the continuum windows fall in the range of $\lambda_{\rm rest}$ = 1220 – 1350 Å, which are close to $\lambda_{\rm rest}$ = 1450 Å, these measurements are not sensitive to the exact value of $\alpha$. The Ly$\alpha$ properties (luminosity, FWHM, and rest-frame equivalent width [EW]) of a quasar with relatively weak continuum emission, such as $J084456.62+022640.5$, was measured with a local continuum defined on the red side of Ly$\alpha$. For the remaining objects with strong continuum, we measured the properties of the broad Ly$\alpha$ + N V $\lambda$1240 complex, with a local continuum defined by the above power-law model. The resultant line properties are summarized in Table 3.

$J235646.33+001747.3$ has a Ly$\alpha$ line peak at 9740 Å, which corresponds to a redshift $z=7.01$. But this redshift may be an overestimate, as the intrinsic Ly$\alpha$ peak may be at a shorter wavelength and absorbed by the IGM. The flux spike at $\sim$9800 Å in the $J011257.84+011042.4$ spectrum is likely due to residual of a strong sky emission line. $J151248.71+442217.5$ has a relatively red continuum, and indeed the HSC magnitudes from the latest DR ($i-z=1.5$ and $z-y=0.5$) indicates a Bayesian quasar probability of $P_{\rm Q}^{\rm B}=0.000$ (see Table 1); this object was selected from an older DR, which happened to give a higher $P_{\rm Q}^{\rm B}$.

Only three quasars are detected in the near-IR bands (see Table 2) and, interestingly, two of them ($J113753.64+004509.7$ and $J142307.21-022519.0$) apparently lack strong Ly$\alpha$ in emission. We re-examined all the 75 broad-line quasars we discovered so far, and found that out of the eight quasars with near-IR detection, five are such weak-line quasars. This is significantly higher than the average fraction of weak-line quasars among all 75 quasars ($\sim$20 %; Y. Matsuoka et al., in prep.), or among more luminous quasars in the literature ($\sim$10 % Diamond-Stanic et al., 2009; Bañados et al., 2016). This may suggest a link between observed Ly$\alpha$ weakness and the continuum shape at longer wavelength. We defer further analysis of this topic to a future paper.

Figure 4 presents 7 high-$z$ objects with no or weak Ly$\alpha$ emission line, which are most likely galaxies at $z\sim 6$. Redshifts of these galaxies were estimated from the observed wavelengths of Ly$\alpha$, the interstellar absorption lines of Si II $\lambda$1260, Si II $\lambda$1304, C II $\lambda$1335, and/or the onset of the Gunn & Peterson (1965) trough. Due to the limited S/N of the spectra, the estimated redshifts should be regarded as only approximate (${\Delta}z\lesssim 0.1$). The absolute magnitudes and Ly$\alpha$ properties were measured in the same way as for quasars, except that we assumed a continuum slope of $\beta=-2.0$ ($F_{\lambda}\propto\lambda^{\beta}$; Stanway et al., 2005).

Our quasar survey explores the luminosities where quasars and galaxies have comparable number densities (Matsuoka et al., 2018c), and hence contamination of galaxies is inevitable. Given the limited wavelength coverage of our spectra ($\lambda_{\rm rest}\sim 1200-1500$ Å), high-$z$ objects without broad Ly$\alpha$ emission are difficult to classify unambiguously into quasars or galaxies. Therefore the above quasar/galaxy classification is not perfect, and it may change in the future, with new data providing additional S/N or wavelength coverage.

In addition to the above high-$z$ objects, we found 10 [O III] emitters at $z\sim 0.8$, as displayed in Figure 5. Their strong [O III] lines contribute significantly to the HSC $z$-band magnitude, and thus mimic colors of $z\sim 6$ quasars. We measured the properties of the H$\gamma$, H$\beta$, [O III] $\lambda$4959 and $\lambda$5007 emission lines, as listed in Table 3. Since these galaxies have very weak continuum, we estimated the continuum levels by summing up all available pixels after masking the above emission lines. As we discussed in Matsuoka et al. (2018a), their extremely high [O III] $\lambda$5007/H$\beta$ ratios may indicate that these are galaxies with sub-solar metallicity and high ionization state of the interstellar medium (e.g., Kewley et al., 2016), and/or contribution from an active galactic nucleus (AGN).

The remaining 22 objects presented in Figures 6 – 7 were found to be Galactic cool dwarfs (low-mass stars and brown dwarfs). Their rough spectral classes were estimated by fitting the spectral standard templates of M4- to T8-type dwarfs, taken from the SpeX Prism Spectral Library (Burgasser, 2014; Skrzypek et al., 2015), to the observed spectra at $\lambda=7500-10200$ Å. The results are summarized in Table 4 and plotted in the figures. Due to the low S/N and limited wavelength coverage of the spectra, the classifications presented here are rather uncertain, and should be regarded as only approximate.

Figure 8 presents shape indicators of the HSC objects with spectroscopic identification from our survey. The horizontal axis uses second-order adaptive moments in two image directions ($\mu_{11}$ and $\mu_{22}$; Hirata & Seljak, 2003), which should be equal to those of PSF ($\mu_{11}^{\rm PSF}$ and $\mu_{22}^{\rm PSF}$) for ideal point sources. The figure shows that quasars with narrow Ly$\alpha$ lines are slightly more extended than the other quasars, which indicates contribution of stellar emission (see the discussion further below). We could eliminate some of the galaxy contamination by a stricter cut of point source selection (e.g., $m_{\rm PSF}-m_{\rm CModel}<0.10$) than currently used, but that would also reject those potentially important population of narrow-line objects. On the other hand, the figure suggests that the adoptive moments ($\mu_{11}$ and $\mu_{22}$) may be less affected by catastrophic measurement errors, and hence may be a better shape indicator, than the magnitude difference ($m_{\rm PSF}-m_{\rm CModel}$); this option will be considered in future candidate selection.

Figure 9 presents the composite spectra of all 75 broad-line quasars, 18 narrow-line quasars, and 31 galaxies discovered by our survey. These spectra were generated by converting the individual spectra to the rest frame and normalizing to $M_{1450}=-23$ mag, and then median-stacking. Thus the individual objects have equal weights in the stacking, regardless of the brightness or spectral S/N. The overall shape of the quasar spectrum is similar to that of the local quasar composite (Vanden Berk et al., 2001), except for the narrower Ly$\alpha$ line and the IGM absorption. Our composite spectrum is also similar in shape to a composite spectrum of high-$z$ luminous quasars from the Panoramic Survey Telescope & Rapid Response System 1 (Pan-STARRS1; Chambers et al., 2016) survey (Bañados et al., 2016), again except for the narrower Ly$\alpha$ line, which may reflect the Baldwin (1977) effect. While a clear sign of the quasar near-zone effect is visible just blueward of Ly$\alpha$, it should be noted that the redshifts of most of our quasars have been determined with Ly$\alpha$, and are thus not very accurate. A detailed analysis of emission and absorption profiles around Ly$\alpha$, including the effect of IGM damping-wing absorption, must wait for accurate measurements of systemic redshifts via other emission lines.

The composite spectrum of the narrow-line quasars is characterized by the strong and asymmetric Ly$\alpha$ line. The rest-frame EW of this line is 29 $\pm$ 2 Å in the composite, while it ranges up to $\sim$500 Å in the individual spectra. The spectrum shows a P Cygni-like profile around 1240 Å, which is likely due to N V $\lambda$1240 emission line and an associated mini broad absorption line system. This spectral feature, along with the very luminous Ly$\alpha$ line usually associated with AGN ($L_{\rm Ly\alpha}>10^{43}$ erg s^-1; Konno et al., 2016) and the absence of interstellar absorption in the continuum, strongly suggest that these objects are narrow-line quasars (see also the discussion in Matsuoka et al., 2018a). Compared with narrow-line quasar candidates in SDSS (Alexandroff et al., 2013), these high-$z$ objects have apparently narrower Ly$\alpha$, indeed many are spectrally unresolved (most of the narrow-line quasars were observed with Subaru/FOCAS, which has a resolving power of $\sim$250 km s^-1 with our observing mode). This may suggest a significant contribution from the host galaxies, thus these objects may be composites of quasars (or AGNs) and star-forming galaxies.

We note that a color selection of high-$z$ quasars (including a more sophisticated Bayesian selection we used) is generally more sensitive to objects with stronger Ly$\alpha$ lines. Thus the composite spectra may not represent the whole populations of broad- and narrow-line quasars that reside in the high-$z$ universe. On the other hand, we confirmed in a previous work (LF paper; Matsuoka et al., 2018c) that our quasar selection is fairly complete, as long as high-$z$ quasars have a similar distribution of intrinsic spectral shapes to low-$z$ SDSS quasars.

Finally, the galaxy spectrum in Figure 9 shows strong absorption lines, with rest-frame EWs of $1.7\pm 0.8$ Å, $1.8\pm 0.7$ Å, and $2.2\pm 1.6$ Å for Si II $\lambda$1260, Si II $\lambda$1304, and C II $\lambda$1335, respectively. These are broadly consistent with the values measured in a composite spectrum of $z\sim 3$ Lyman break galaxies (Shapley et al., 2003). The galaxies discovered by our survey are very luminous, with $M_{1450}\simeq-22$ to $-25$ mag, well above the break magnitude of the galaxy LF at $z=6$ ($M_{1450}\simeq-21$ mag; Ono et al., 2018).

We also present a composite spectrum of the 16 [O III] emitters at $z\sim 0.8$ in Figure 10. This spectrum was created by converting the individual spectra to the rest frame and normalizing to the continuum flux of $10^{18}$ erg s^-1 cm^-2 Å^-1, and then median-stacking. In addition to H$\gamma$, H$\beta$, [O III] $\lambda$4959 and $\lambda$5007, we detected [O III] $\lambda$4363 and He I $\lambda$4473 emission; these two lines are weak and undetected in the individual spectra. The rest-frame EWs of the six emission lines are 53 $\pm$ 3 Å (H$\gamma$), 35 $\pm$ 3 Å ([O III] $\lambda$4363), 21 $\pm$ 3 Å (He I $\lambda$4473), 165 $\pm$ 7 Å (H$\beta$), 264 $\pm$ 10 Å ([O III] $\lambda$4959), 1030 $\pm$ 40 Å ([O III] $\lambda$5007). All lines are spectrally unresolved.

The HSC-SSP survey has completed observations on more than 80 % of the planned 300 nights, and we are making steady progress on our high-$z$ quasar survey. We plan to continue our follow-up spectroscopy in the next few years, and will report new discoveries, along with measurements of the properties of individual objects with multi-wavelength observations. The quasar sample thus established will provide a more accurate quasar LF at $z\geq 6$ than is currently available (e.g., Matsuoka et al., 2018c; Wang et al., 2018). In a few years, we also aim to start a large SSP survey with the Prime Focus Spectrograph, a new wide-field multi-object spectrograph for the Subaru Telescope under development (Takada et al., 2014). This will enable us to take spectra of a significant number of HSC objects, including quasar candidates at all redshifts.