SILVERRUSH. IV. Ly$\alpha$ Luminosity Functions at $z = 5.7$ and $6.6$ Studied with $\sim$ 1,300 LAEs on the $14-21$ deg$^2$ Sky

TL;DR

This study derives precise Lyα luminosity functions at redshifts 5.7 and 6.6 using a large sample of LAEs from Subaru/HSC data, providing insights into galaxy evolution and reionization.

Contribution

It presents the first Lyα luminosity functions at these redshifts with unprecedented accuracy, addressing previous uncertainties and systematic effects.

Findings

Lyα LFs are consistent with previous studies.

Faint-end slope is very steep (~ -2.5).

Estimated neutral hydrogen fraction at z=6.6 is 0.3 ± 0.2.

Abstract

We present the Ly$\alpha$ luminosity functions (LFs) at $z=$5.7 and 6.6 derived from a new large sample of 1,266 Ly$\alpha$ emitters (LAEs) identified in total areas of 14 and 21 deg$^2$, respectively, based on the early narrowband data of the Subaru/Hyper Suprime-Cam (HSC) survey. Together with careful Monte-Carlo simulations that account for the incompleteness of the LAE selection and the flux estimate systematics in the narrowband imaging, we have determined the Ly$\alpha$ LFs with the unprecedentedly small statistical and systematic uncertainties in a wide Ly$\alpha$ luminosity range of $10^{42.8-43.8}$ erg s$^{-1}$. We obtain the best-fit Schechter parameters of $L^{*}_{{\rm Lya}}=1.6^{+2.2}_{-0.6} (1.7^{+0.3}_{-0.7}) \times10^{43}$ erg s$^{-1}$, $\phi^{*}_{{\rm Lya}}=0.85^{+1.87}_{-0.77}\ (0.47^{+1.44}_{-0.44})\times10^{-4}$ Mpc$^{-3}$, and $\alpha=-2.6^{+0.6}_{-0.4}\ (-2.5^{+0.5}_{-0.5})$ at $z=5.7$ ($6.6$). We confirm that our best-estimate Ly$\alpha$ LFs are consistent with the majority of the previous studies, but find that our Ly$\alpha$ LFs do not agree with the high number densities of LAEs recently claimed by Matthee/Santos et al.'s studies that may overcorrect the incompleteness and the flux systematics. Our Ly$\alpha$ LFs at $z=5.7$ and $6.6$ show an indication that the faint-end slope is very steep ($\alpha \simeq -2.5$), although it is also possible that the bright-end LF results are enhanced by systematic effects such as the contribution from AGNs, blended merging galaxies, and/or large ionized bubbles around bright LAEs. Comparing our Ly$\alpha$ LF measurements with four independent reionization models, we estimate the neutral hydrogen fraction of the IGM to be $x_{\rm HI}=0.3\pm0.2$ at $z=$6.6 that is consistent with the small Thomson scattering optical depth obtained by Planck 2016.

Tables

Table 1. Table 1: Summary of HSC/ NB816 and NB921 Data

Field	Area (NB816)	Area (NB921)	$g$a	$r$a	$i$a	NB816a	$z$a	NB921a	$y$a
	(deg2)	(deg2)	(ABmag)	(ABmag)	(ABmag)	(ABmag)	(ABmag)	(ABmag)	(ABmag)
UD-COSMOS	1.97	2.05	26.9	26.6	26.2	25.7	25.8	25.6	25.1
UD-SXDS	1.93	2.02	26.9	26.4	26.3	25.5	25.6	25.5	24.9
D-COSMOS	—	5.31	26.5	26.1	26.0	—	25.5	25.3	24.7
D-DEEP2-3	4.37	5.76	26.6	26.2	25.9	25.2	25.2	24.9	24.5
D-ELAIS-N1	5.56	6.08	26.7	26.0	25.7	25.3	25.0	25.3	24.1
Total	13.8	21.2	—	—	—	—	—	—	—

Table 2. Table 2: Photometric Sample of z = 5.7 𝑧 5.7 z=5.7 and 6.6 6.6 6.6 LAEs

Field	LAE All samplea	LAE Ly$\alpha$ LF sampleb
The $z=5.7$ LAE sample
UD-COSMOS	201	201
UD-SXDS	224	224
D-DEEP2-3	423	423
D-ELAIS-N1	229	229
Total	1077	1077
The $z=6.6$ LAE sample
UD-COSMOS	338	50
UD-SXDS	58	21
D-COSMOS	244	48
D-DEEP2-3	164	38
D-ELAIS-N1	349	32
Total	1153	189

Table 3. Table 3: Ly α 𝛼 \alpha LFs at z = 5.7 𝑧 5.7 z=5.7 and 6.6 6.6 6.6 from this work

$\log L\left(\mathrm{Ly}\alpha \right)$a	$\log \varphi$b
(erg s-1)	([$\Delta \log L\left(\mathrm{Ly}\alpha \right)$]-1 Mpc-3)
$z=5.7$
$42.95$	$-{3.478}_{-0.193}^{+0.038}$
$43.05$	$-{3.735}_{-0.199}^{+0.044}$
$43.15$	$-{3.953}_{-0.198}^{+0.043}$
$43.25$	$-{4.163}_{-0.210}^{+0.055}$
$43.35$	$-{4.427}_{-0.231}^{+0.076}$
$43.45$	$-{4.970}_{-0.308}^{+0.147}$
$43.55$	$-{5.170}_{-0.355}^{+0.187}$
$43.65$	$-{5.318}_{-0.401}^{+0.224}$
$43.75$	$-{5.717}_{-0.606}^{+0.365}$
$z=6.6$
$43.15$	$-{4.194}_{-0.317}^{+0.154}$
$43.25$	$-{4.407}_{-0.258}^{+0.101}$
$43.35$	$-{4.748}_{-0.243}^{+0.087}$
$43.45$	$-{5.132}_{-0.300}^{+0.140}$
$43.55$	$-{5.433}_{-0.374}^{+0.203}$
$43.65$	$-{5.609}_{-0.438}^{+0.253}$
$43.75$	$-{6.212}_{-0.917}^{+0.519}$
$43.85$	$-{6.226}_{-0.917}^{+0.519}$

Table 4. Table 4: Best-fit Schechter Parameters and Ly α 𝛼 \alpha Luminosity Densities

Redshift	$L^{\ast }$	${\varphi }^{\ast }$	$\alpha$	${\rho }_{\mathrm{obs}}^{\mathrm{Ly}\alpha }$a
	(${10}^{43}$ erg s-1)	(${10}^{-4}$ Mpc-3)		(${10}^{39}$ erg s-1 Mpc-3)
Classical method for $\log L\left(\mathrm{Ly}\alpha \right)$ [erg s-1] $=42.4-43.5$b
5.7	${1.07}_{-0.38}^{+0.77}$	${2.46}_{-1.86}^{+3.48}$	$-{2.26}_{-0.44}^{+0.76}$	$3.39$
6.6	${0.82}_{-0.30}^{+0.86}$	${2.83}_{-2.38}^{+3.52}$	$-{1.86}_{-0.67}^{+0.79}$	$1.96$
Classical method for $\log L\left(\mathrm{Ly}\alpha \right)$ [erg s-1] $=42.4-44.0$c
5.7	${1.64}_{-0.62}^{+2.16}$	${0.849}_{-0.771}^{+1.87}$	$-{2.56}_{-0.45}^{+0.53}$	$3.49$
6.6	${1.66}_{-0.69}^{+0.30}$	${0.467}_{-0.442}^{+1.44}$	$-{2.49}_{-0.50}^{+0.50}$	$1.82$
End-to-end Monte Carlo simulations
5.7	$2.0$	$0.63$	$-2.6$ (fix)	$3.5$
6.6	$1.3$	$0.63$	$-2.5$ (fix)	$1.7$

Equations

$$\begin{array}[]{l}i-\textit{NB816}\geq 1.2,\\ g>g_{3\sigma},\\ \mathrm{and}\ [(r\leq r_{3\sigma}\ \mathrm{and}\ r-i\geq 1.0)\ \mathrm{or}\ (r>r_{3\sigma})]\end{array}$$

$$\begin{array}[]{l}z-\textit{NB921}\geq 1.0,\\ g>g_{3\sigma},\\ r>r_{3\sigma},\\ \mathrm{and}\ [(z\leq z_{3\sigma}\ \mathrm{and}\ i-z\geq 1.0)\ \mathrm{or}\ (z>z_{3\sigma})]\end{array}$$

$$\varphi (L_{\mathrm{Ly}\alpha })d{L}_{\mathrm{Ly}\alpha }={\varphi }_{\mathrm{Ly}\alpha }^{\ast }{\left(\frac{L_{\mathrm{Ly}\alpha }}{L_{\mathrm{Ly}\alpha }^{\ast }}\right)}^{\alpha }\exp \left(-\frac{L_{\mathrm{Ly}\alpha }}{L_{\mathrm{Ly}\alpha }^{\ast }}\right)d\left(\frac{L_{\mathrm{Ly}\alpha }}{L_{\mathrm{Ly}\alpha }^{\ast }}\right),$$

$${\rho }^{\mathrm{Ly}\alpha }=\kappa T_{\mathrm{Ly}\alpha ,z}^{\mathrm{IGM}}{f}_{\mathrm{Ly}\alpha }^{\mathrm{esc}}{\rho }^{\mathrm{UV}},$$

$$\frac{T_{\mathrm{Ly}\alpha ,z=6.6}^{\mathrm{IGM}}}{T_{\mathrm{Ly}\alpha ,z=5.7}^{\mathrm{IGM}}}=\frac{{\kappa }_{z=5.7}}{{\kappa }_{z=6.6}}\frac{f_{\mathrm{Ly}\alpha ,z=5.7}^{\mathrm{esc}}}{f_{\mathrm{Ly}\alpha ,z=6.6}^{\mathrm{esc}}}\frac{{\rho }_{z=6.6}^{\mathrm{Ly}\alpha }/{\rho }_{z=5.7}^{\mathrm{Ly}\alpha }}{{\rho }_{z=6.6}^{\mathrm{UV}}/{\rho }_{z=5.7}^{\mathrm{UV}}}.$$

Full text

1 Introduction

Ly$\alpha$ emission lines are one of the key properties of galaxies for exploring a high-$z$ universe. Ly$\alpha$ emitters (LAEs), which generally have a spectrum of a luminous Ly$\alpha$ line and a faint ultraviolet (UV) continuum, have been found at a wide redshift range of $z=0-8$ by several approaches including narrowband surveys (e.g., [Cowie & Hu (1998)]; [Hu et al. (1998)]; [Rhoads et al. (2000)]; [Steidel et al. (2000)]; [Malhotra & Rhoads (2002)]; [Ajiki et al. (2002)]; [Ouchi et al. (2003)]; [Hayashino et al. (2004)]; [Matsuda et al. (2004)]; [Taniguchi et al. (2005)]; [Iye et al. (2006)]; [Kashikawa et al. (2006)]; [Shimasaku et al. (2006)]; [Gronwall et al. (2007)]; [Murayama et al. (2007)]; [Guaita et al. (2010)]; [Shibuya et al. (2012)]; [Yamada et al. (2012)]; [Konno et al. (2014)]) and spectroscopic observations (e.g., [Deharveng et al. (2008)]; [Adams et al. (2011)]; [Finkelstein et al. (2013)]; [Schenker et al. (2014)]; [Cassata et al. (2015)]; [Oesch et al. (2015)]; [Zitrin et al. (2015)]; [Song et al. (2016)]; [Stark et al. (2017)]). From these observations, it has been revealed that LAEs are in an early phase of galaxy evolution, i.e., LAEs are young, less massive, less dusty, and in highly ionized state (e.g., [Ono et al. (2010b)]; [Ono et al. (2010a)]; [Nakajima & Ouchi (2014)]; [Kusakabe et al. (2015)]; [Inoue et al. (2016)]).

Ly$\alpha$ luminosity functions (LFs) and their evolution can be a probe for the early evolution of galaxies and cosmic reionization (e.g., [Haiman & Spaans (1999)]; [McQuinn et al. (2007)]; [Mao et al. (2007)]; [Kobayashi et al. (2007)]; [Mesinger & Furlanetto (2008)]; [Dayal et al. (2011)]). Previous studies have found that Ly$\alpha$ LFs increase from $z\sim 0$ to $z\sim 3$, show a moderate plateau between $z\sim 3$ to $z\sim 6$, and decrease toward $z\gtrsim 6$ (e.g., [Deharveng et al. (2008)]; [Ouchi et al. (2008)]; [Kashikawa et al. (2011)]). The evolution of Ly$\alpha$ LFs is different from that of UV LFs, which increases from $z\sim 0$ to $z\sim 2$, and turns to the decrease beyond $z\gtrsim 3$ (e.g., [Schiminovich et al. (2005)]; [Reddy & Steidel (2009)]; [Bouwens et al. (2015b)]; see also Figure 7 of [Konno et al. (2016)]). The difference of the evolutionary trend between Ly$\alpha$ and UV LFs would be related to the escaping process of Ly$\alpha$ photons not only from the Hi ISM of a galaxy, but also from the Hi intergalactic medium (IGM). The Ly$\alpha$ escape fraction, $f^{\mathrm{Ly}\alpha}_{\mathrm{esc}}$, which is defined by the ratio of the star formation rate densities (SFRDs) estimated from observed Ly$\alpha$ luminosity densities (LDs) to those estimated from intrinsic UV LDs, largely increases from $z\sim 0$ to $z\sim 6$ by two orders of magnitudes, and turns to the decrease beyond $z\gtrsim 6$ (e.g., [Hayes et al. (2011)]). The rapid evolution of the Ly$\alpha$ escape fraction from $z\sim 6$ to $z\sim 0$ would be explained by the combination of the Ly$\alpha$ attenuation by dust and the Ly$\alpha$ resonance scattering effect by Hi in ISM. In the case that the ISM Hi density of a galaxy is large, the path lengths of Ly$\alpha$ photons become longer due to the resonant scattering, and these Ly$\alpha$ photons are subject to the attenuation by dust. Konno et al. (2016) have used simple expanding shell models, which compute the Ly$\alpha$ radiative transfer by Monte Carlo simulations (MCLya; Verhamme et al. (2006); Schaerer et al. (2011)), and have suggested that the large increase of Ly$\alpha$ escape fraction at $z=0-6$ can be reproduced by the combination of the Hi column density decrease (by two orders of magnitude) and the average dust extinction values. The decrease of the Ly$\alpha$ LFs at $z\gtrsim 6$ is related to the cosmic reionization, because the Ly$\alpha$ damping wing of Hi in IGM attenuates Ly$\alpha$ photons from a galaxy. Previous studies have found that Ly$\alpha$ LFs at $z\sim 7$ significantly decrease from those at $z\sim 6$ (e.g., Kashikawa et al. (2006); Ouchi et al. (2010); Hu et al. (2010); Santos et al. (2016)), and especially at $z\gtrsim 7$, Ly$\alpha$ LFs decrease rapidly (e.g., Konno et al. (2014)). The neutral hydrogen fraction of IGM, $x_{\mathrm{HI}}$, can be estimated by the Ly$\alpha$ LD evolution subtracting the galaxy evolution effect. Ouchi et al. (2010) have constrained $x_{\mathrm{HI}}=0.2\pm 0.2$ at $z=6.6$ from the Ly$\alpha$ LF evolution at $z=5.7-6.6$ (see also Malhotra & Rhoads (2004); Kashikawa et al. (2006)). Similarly, the neutral hydrogen fractions at $z\gtrsim 7$ have also been estimated from the Ly$\alpha$ LF evolution (Ota et al. (2010); Konno et al. (2014); Ota et al. (2017)). These $x_{\mathrm{HI}}$ estimates could constrain the history of cosmic reionization by the comparison with the Thomson scattering optical depth of cosmic microwave background (CMB).

Recently, a large number of wide-field narrowband imaging surveys have been conducted not only to spread the Ly$\alpha$ luminosity ranges of Ly$\alpha$ LFs, but also to reveal physical properties for luminous LAEs. At $z\sim 2-3$, luminous LAEs are known to have counterparts in multiwavelength data (e.g., X-ray and radio) and/or extended Ly$\alpha$ haloes (e.g., Steidel et al. (2000); Ouchi et al. (2008); Cantalupo et al. (2014); Cai et al. (2017)). A recent study, for example, has confirmed that there are excesses found in Ly$\alpha$ LFs at $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $\gtrsim 43.4$, and the excesses are made by (faint) AGNs based on multiwavelength imaging data (Konno et al., 2016). Interestingly, such luminous LAEs have also been discovered at a higher redshift of $z\sim 6.6$ (e.g., Himiko by Ouchi et al. (2009), CR7 and MASOSA by Sobral et al. (2015), and COLA1 by Hu et al. (2016); see also IOK-1 by Iye et al. (2006)). A number of observational and theoretical studies have aimed to uncover the physical origins of these bright LAEs (e.g., Ouchi et al. (2013) and Zabl et al. (2015) for Himiko; Bowler et al. (2017b), Pacucci et al. (2017), and Shibuya et al. (2017b) for CR7).

In this paper, we present the Ly$\alpha$ LFs at $z=5.7$ and $6.6$ based on the Subaru/Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP; Aihara et al. (2017b)). Because the field of view of HSC is about seven times wider than that of Subaru/Suprime-Cam, HSC can identify a large number of high-$z$ LAEs with a wide range of Ly$\alpha$ luminosity more efficiently than Suprime-Cam. In our HSC SSP survey, a total of $\sim 13.8$ deg2 and $\sim 21.2$ deg2 sky areas are covered by NB816 and NB921 observations, respectively (see also Section 2.1, Ouchi et al. (2017) and Shibuya et al. (2017a) for details). These wide field HSC NB data sets allow us to determine the Ly$\alpha$ LFs at $z=5.7$ and $6.6$ with unprecedented accuracy. By examining the evolution of these Ly$\alpha$ LFs at $z=5.7-6.6$, we can constrain the $x_{\mathrm{HI}}$ value at $z=6.6$. Moreover, based on these HSC SSP data, we can push the Ly$\alpha$ luminosity range toward brighter luminosity, and investigate the abundance of luminous high-$z$ LAEs. We describe a summary of our HSC surveys and the sample construction for $z=5.7$ and $6.6$ LAEs in Section 2. We derive the Ly$\alpha$ LFs at these redshifts, and compare the Ly$\alpha$ LFs with those of previous studies in Section 3. We examine the Ly$\alpha$ LF evolution at $z=5.7-6.6$, and discuss cosmic reionization in Section 4. This paper is placed in a series of papers from twin programs studying high-$z$ objects based on the HSC SSP data products. One program is our high-$z$ LAE studies named Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH). This program provides the clustering measurements of $z=5.7$ and $6.6$ LAEs (Ouchi et al., 2017), the photometric and spectroscopic properties of LAEs at these redshifts (Shibuya et al., 2017a, b), the systematic survey for LAE overdense region (R. Higuchi et al. in preparation), and our Ly$\alpha$ LF studies. The other program is the high-$z$ dropout galaxy study, Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH; Ono et al. (2017); Harikane et al. (2017); Toshikawa et al. (2017)). Throughout this paper, we use magnitudes in the AB system (Oke, 1974). We adopt $\Lambda$CDM cosmology with a parameter set of ($h$, $\Omega_{m}$, $\Omega_{\Lambda}$, $\sigma_{8}$) = ($0.7$, $0.3$, $0.7$, $0.8$), which is consistent with the nine-year WMAP and the latest Planck results (Hinshaw et al., 2013; Planck Collaboration et al., 2016a).

2 Observations and Sample Selection

2.1 Hyper Suprime-Cam Imaging Observations and

Data Reduction

In our sample construction for $z=5.7$ and $6.6$ LAEs, we use narrowband (NB816, NB921) imaging data as well as broadband ($g,r,i,z,y$) imaging data, which are taken with Subaru/HSC (Miyazaki et al. (2012); see also Miyazaki et al. (2017); Furusawa et al. (2017); Kawanomoto et al. (2017); Komiyama et al. (2017)). The narrowband filters, NB816 and NB921, have central wavelengths of 8170 Å and 9210 Å, respectively, and FWHMs of 131 Å and 120 Å to identify LAEs in the redshift range of $z=5.67-5.77$ and $z=6.52-6.63$, respectively. We show the response curves of the narrowband filters as well as the broadband filters in Figure 1. These narrowband and broadband images are obtained in our ongoing HSC legacy survey under the Subaru Strategic Program (SSP; PI: S.Miyazaki, see also Aihara et al. (2017b)). The HSC SSP has been allocated 300 nights over 5 years, and started in March 2014. The HSC SSP survey has three layers with different sets of area and depth: the Wide, Deep, and UltraDeep layers. These layers will cover the sky area of $\sim 1400$ deg2, $\sim 30$ deg2, and $\sim 4$ deg2 with the $5\sigma$ limiting magnitudes (in r band) of $\sim 26$ mag, $\sim 27$ mag, and $\sim 28$ mag, respectively. While the broadband images are taken in all the three layers, the NB816 and NB921 images are obtained only in the Deep and UltraDeep layers. We use early datasets of the HSC SSP survey taken from March 2014 to April 2016 (S16A), where all additional data taken in January to April 2016 have been merged with the data of Public Data Release 1 (Aihara et al., 2017a). With the NB816 filter, the HSC SSP survey has observed two blank fields in the Deep layer, the D-DEEP2-3 ($23^{\mathrm{h}}30^{\mathrm{m}}00^{\mathrm{s}}$, $+00^{\mathrm{d}}00^{\prime}00\farcs 0$) and D-ELAIS-N1 ($16^{\mathrm{h}}10^{\mathrm{m}}00^{\mathrm{s}}$, $+54^{\mathrm{d}}00^{\prime}00\farcs 0$) fields, and two blank fields in the UltraDeep layer, the UD-COSMOS ($10^{\mathrm{h}}00^{\mathrm{m}}29^{\mathrm{s}}$, $+02^{\mathrm{d}}12^{\prime}21\farcs 0$) and UD-SXDS ($02^{\mathrm{h}}18^{\mathrm{m}}00^{\mathrm{s}}$, $-05^{\mathrm{d}}00^{\prime}00\farcs 0$) fields. For the NB921 filter, a blank field of the D-COSMOS ($10^{\mathrm{h}}00^{\mathrm{m}}29^{\mathrm{s}}$, $+02^{\mathrm{d}}12^{\prime}21\farcs 0$) field in the Deep layer has also been observed as well as the four fields described above. Each field in the Deep layer is covered by three or four pointing positions of HSC, while in the UltraDeep layer, each field is covered by one pointing position of HSC. The details of our HSC SSP survey is listed in Table 1.

The HSC data are reduced by the HSC SSP survey team with hscPipe (Bosch et al., 2017), which is based on the Large Synoptic Survey Telescope (LSST) pipeline (Ivezic et al., 2008; Axelrod et al., 2010; Jurić et al., 2015). This HSC pipeline performs CCD-by-CCD reduction, calibrates astrometry, mosaic-stacking, and photometric zeropoints, and generates catalogs for sources detected and photometrically measured in the stacked images. The photometric and astrometric calibrations are based on the data from the Panoramic Survey Telescope and Rapid Response System 1 imaging survey (Pan-STARRS1; Schlafly et al. (2012); Tonry et al. (2012); Magnier et al. (2013)). In the stacked images, regions contaminated with diffraction spikes and halos of bright stars are masked by using the mask extension outputs of the HSC pipeline (Coupon et al., 2017). After the masking, the total effective survey areas in the S16A data are $13.8$ deg2 and $21.2$ deg2 for NB816 and NB921, respectively. These survey areas are $70-87$ times larger than those of the Subaru Deep Field studies (Shimasaku et al. (2006); Kashikawa et al. (2011)), $14-21$ times larger than those of the Subaru/XMM-Newton Deep Survey (Ouchi et al. (2008); Ouchi et al. (2010)), and $2-5$ times larger than those of other subsequent studies with Subaru/Suprime-Cam (Matthee et al. (2015); Santos et al. (2016)). Under the assumption of a simple top-hat selection function for LAEs whose redshift distribution is defined by the FWHM of a narrowband filter, these survey areas correspond to comoving volumes of $\simeq 1.16\times 10^{7}$ Mpc3 and $\simeq 1.91\times 10^{7}$ Mpc3 for $z=5.7$ and $6.6$ LAEs, respectively. The narrowband images reach the $5\sigma$ limiting magnitudes in a $1\farcs 5$-diameter circular aperture of $24.9-25.3$ mag in the Deep layer, and $25.5-25.7$ mag in the UltraDeep layer. Note that the PSF sizes of the HSC images are typically $<0\farcs 8$, which is sufficiently smaller than the aperture diameter of $1\farcs 5$ (see Aihara et al. (2017a) for details). We summarize the $5\sigma$ limiting magnitudes of the NB816 and NB921 images in Table 1. For the total magnitudes, we use cmodel magnitudes. The cmodel magnitude is derived from a linear combination of exponential and de Vaucouleurs profile fits to the light profile of each object (Bosch et al., 2017). We make use of the cmodel magnitudes for color measurements, because the HSC data used in this study are reduced with no smoothing to equalize the PSFs and fixed aperture photometry does not provide good measurements of object colors (Aihara et al., 2017a). The total magnitudes and colors are corrected for Galactic extinction (Schlegel et al., 1998).

2.2 Photometric Samples of $z=5.7$ and $6.6$ LAEs

LAE samples at $z=5.7$ and $6.6$ are constructed based on narrowband color excess by Ly$\alpha$ emission, $i-\textit{NB816}$ and $z-\textit{NB921}$, respectively, and no detection of blue continuum fluxes. We first select objects with magnitudes brighter than the $5\sigma$ limit in NB816 or NB921 from the HSC SSP database. We then apply similar selection criteria to those of Ouchi et al. (2008) and Ouchi et al. (2010):

[TABLE]

for $z=5.7$ LAEs, and

[TABLE]

for $z=6.6$ LAEs, where ($g_{3\sigma}$, $r_{3\sigma}$, $z_{3\sigma}$) are the $3\sigma$ limiting magnitudes of ($g$, $r$, $z$) bands. Note that the criterion in the former parentheses of the third criterion in Equation (1) and the fourth criterion in Equation (2) are used to select bright objects whose SED is consistent with a Lyman break due to intergalactic absorption. In addition to the color selection criteria, we use the countinputs parameter, which represents the number of exposures for each object in each band. We apply countinputs $\geq 3$ for the narrowband images. We also remove objects affected by bad pixels, proximity to bright stars, or poor photometric measurement by using the following flags: flags_pixel_edge, flags_pixel_interpolated_center, flags_pixel_saturated_center, flags_pixel_cr_center, and flags_pixel_bad. After the visual inspection for the rejection of spurious sources and cosmic rays, we identify 1,081 and 1,273 LAE candidates at $z=5.7$ and $6.6$, respectively (Shibuya et al., 2017a). The samples of these LAE candidates are referred to as the ‘LAE All’ samples. The LAE All samples are $\sim 2-6$ times larger than photometric samples in previous studies (e.g., Ouchi et al. (2008); Ouchi et al. (2010); Matthee et al. (2015); Santos et al. (2016)). This sample is used for clustering analyses in our companion paper (Ouchi et al., 2017). The details of the sample construction including the color-magnitude diagrams of NB$-$BB vs. NB are presented in Shibuya et al. (2017a).

In this Ly$\alpha$ LF study, we create subsamples of the LAE All samples to directly compare our results with previous work. The only difference between the subsamples and the LAE All samples is the $z-$NB921 color criterion for $z=6.6$ LAEs. The color selection criterion for $z=5.7$ LAEs (i.e., $\textit{i}-\textit{NB816}>1.2$ in Equation 1) corresponds to the rest-frame Ly$\alpha$ equivalent width (EW), EW0, of $\mathrm{EW}_{0}\gtrsim 10$Å in the case of a flat UV continuum (i.e., $f_{\nu}=$ const.) with IGM attenuation (Madau, 1995). This EW limit is similar to those of previous studies ($\mathrm{EW}_{0}\gtrsim 10-30$Å; e.g., Shimasaku et al. (2006); Ouchi et al. (2008); Santos et al. (2016)). Thus, the $z=5.7$ LAE sample of the LAE All samples can be used for comparison with the previous Ly$\alpha$ LF results. On the other hand, the color criterion of $z-\textit{NB921}>1.0$ in Equation (2) for $z=6.6$ LAEs corresponds to the EW0 limit significantly lower than those of previous studies using Subaru/Suprime-Cam (e.g., Ouchi et al. (2010); Matthee et al. (2015)). This is because the relative wavelength position of NB921 to $z^{\prime}$ (or $z$) band filter is different between Suprime-Cam and HSC (Figure 1). Specifically, the central wavelength of the HSC $z$-band filter is $\simeq 160$Å shorter than that of the Suprime-Cam $z^{\prime}$-band filter. For consistency of comparison, we adopt a more stringent color criterion of $z-\textit{NB921}>1.8$. This criterion corresponds to EW${}_{0}>14$Å ($f_{\nu}=$ const.), which is the same as that used in Ouchi et al. (2010). We refer to these $z=5.7$ and $6.6$ LAE samples as the ‘LAE Ly$\alpha$ LF’ samples. We use the LAE Ly$\alpha$ LF samples to derive surface number densities and color distributions (Section 3.3), and Ly$\alpha$ LFs at $z=5.7$ and $6.6$ (Section 3.4). The numbers of our LAE candidates at $z=5.7$ and $6.6$ are summarized in Table 2. Note that the number of $z=5.7$ LAEs found in D-DEEP2-3 is about two times larger than that in D-ELAIS-N1, although the area of D-DEEP2-3 is about 1.3 times smaller than that of D-ELAIS-N1 and the depths of the NB816 data for these two fields are comparable. This is probably because the seeing of the NB816 data for D-DEEP2-3 is better than that for D-ELAIS-N1. This is also the case for the difference of the numbers of $z=6.6$ LAEs between UD-COSMOS and UD-SXDS.

3 Ly$\alpha$ Luminosity Functions

3.1 Detection Completeness

We estimate detection completeness as a function of the NB816 and NB921 magnitude by Monte Carlo simulations with the SynPipe software (Huang et al., 2017; Murata et al., 2017). Using the SynPipe software, we distribute $\sim$ 18,000 pseudo LAEs with various magnitudes in NB816 and NB921 images. These pseudo LAEs have a Sérsic profile with the Sérsic index of $n=1.5$, and the half-light radius of $r_{\mathrm{e}}\sim 0.9$ kpc, which corresponds to $0.15$ and $0.17$ arcsec for $z=5.7$ and $6.6$ sources, respectively. These Sérsic index and half-light radius values are similar to the average ones of $z\sim 6$ LBGs with $L_{\mathrm{UV}}=0.3-1L^{*}_{z=3}$ (Shibuya et al., 2015). We then perform source detection and photometry with hscPipe, and calculate the detection completeness. We define the detection completeness in a magnitude bin as the fraction of the numbers of the detected pseudo LAEs to all of the input pseudo LAEs in the magnitude bin. Figure 2 shows the detection completeness of the NB816 and NB921 images for the D-DEEP2-3 field. We find that the detection completeness is typically $\gtrsim 80$% for bright objects with $\textit{NB}\lesssim 24.5$ mag, and $\sim 40$% at the $5\sigma$ limiting magnitudes of these narrowband images. We correct for the detection completeness to derive the surface number densities and the Ly$\alpha$ LFs of LAEs in Sections 3.3 and 3.4. For the D-DEEP2-3 field, we use the detection completeness shown in Figure 2, and for the other fields, we shift it along the magnitude considering the limiting magnitudes of the narrowband images.

3.2 Contamination

In our companion paper Shibuya et al. (2017b), we estimate the contamination fractions in our $z=5.7$ and $6.6$ LAE samples based on $81$ LAE candidates whose spectroscopic redshifts are obtained by our past and present programs with Subaru/Faint Object Camera and Spectrograph (FOCAS; Kashikawa et al. (2002)), Magellan/Low Dispersion Survey Spectrograph 3 (LDSS3), and Magellan/Inamori Magellan Areal Camera and Spectrograph (IMACS; Dressler et al. (2011)). We find that $28$ ($53$) LAE candidates at $z=6.6$ ($z=5.7$) have been spectroscopically observed and $4$ out of the $28$ ($4$ out of the $53$) LAE candidates are found to be low-$z$ interlopers. Based on these results, the contamination fraction, $f_{\mathrm{cont}}$, is estimated to be $f_{\mathrm{cont}}=4/28\simeq 14$% ($4/53\simeq 8$%) for the $z=6.6$ ($z=5.7$) LAE sample. We also estimate the contamination fractions for bright LAE candidates with NB $<24$ mag. We have spectroscopically observed $18$ bright LAE candidates. Out of the $18$ candidates, $13$ sources are confirmed as LAEs and the other $5$ objects are strong [Oiii] emitters at low $z$. Based on our spectroscopy results, the contamination rates for the bright $z=6.6$ and $z=5.7$ LAE samples are $f_{\mathrm{cont}}\simeq 33$% ($=4/12$) and $\simeq 17$% ($=1/6$), respectively. Although the contamination rates appear to depend on NB magnitude, the estimated values are in the range of around $0-30$% and have large uncertainties due to the small number of our spectroscopically confirmed sources at this early stage of our program. In this study, we take into account this systematic uncertainty by increasing the lower $1\sigma$ confidence intervals of the Ly$\alpha$ LFs by $30$% (see Section 3.4). Note that our estimated $f_{\mathrm{cont}}$ values are similar to those obtained in Ouchi et al. (2008), Ouchi et al. (2010), and Kashikawa et al. (2011) ($f_{\mathrm{cont}}=0-30$%), who have conducted the Subaru/Suprime-Cam imaging survey for LAEs at $z=5.7$ and $6.6$.

3.3 Surface Number Densities and Color Distributions

Figures 3 and 4 represent the LAE surface number densities at $z=5.7$ and $6.6$, respectively, derived with our HSC SSP survey data. We obtain the surface number densities by dividing the number counts of LAEs by our survey areas (Section 2.1). These surface number densities are corrected for the detection completeness (Section 3.1). The $1\sigma$ error bars of the surface number densities are calculated based on the Poisson statistics (Gehrels, 1986), because the number counts of LAEs are small in some bright-end bins and their errors are not well represented by the square root values of the number counts. We use the Poisson single-sided limit values in the columns of “0.8413” in Tables 1 and 2 of Gehrels (1986) for the $1\sigma$ upper and lower confidence intervals, respectively. Note that the surface number densities decrease at faint magnitude bins due to the color-selection incompleteness. For comparison, we show the surface densities at $z=5.7$ and $6.6$ of Ouchi et al. (2008) and Ouchi et al. (2010) in Figures 3 and 4, respectively. These previous studies have conducted deep narrowband imaging surveys for LAEs in the SXDS field, which is the sky region overlapping the UD-SXDS field in our HSC SSP survey. In these figures, we find that our surface densities are broadly consistent with those of Ouchi et al. (2008) and Ouchi et al. (2010).

Figures 5 and 6 show the color distributions of $\textit{i}-\textit{NB816}$ and $z-\textit{NB921}$ for $z=5.7$ and $6.6$ LAEs, respectively. Magnitudes with a detection significance below $2\sigma$ are replaced with the $2\sigma$ limiting magnitudes. Based on Figures 3$-$6, we estimate the best-fit Schechter functions and Ly$\alpha$ EW0 distributions by the Monte Carlo simulations in Section 3.5.

3.4 Ly$\alpha$ Luminosity Functions at $z=5.7$ and $6.6$

We present Ly$\alpha$ LFs at $z=5.7$ and $6.6$ based on our HSC Ly$\alpha$ LF samples constructed in Section 2.2. We derive the Ly$\alpha$ LFs in the same manner as Ouchi et al. (2008) and Ouchi et al. (2010). We calculate the Ly$\alpha$ EW0 values of $z=5.7$ (6.6) LAEs from the magnitudes of NB816 (NB921) and $z$ band, and estimate the Ly$\alpha$ luminosities of LAEs from these EW0 values and the total magnitudes of NB816 (NB921), under the assumption that the spectrum of LAEs has a Ly$\alpha$ line and a flat UV continuum (i.e., $f_{\nu}=$ constant) with the IGM absorption of Madau (1995), following the methods described in Shimasaku et al. (2006), Ouchi et al. (2010), and Konno et al. (2014). Ly$\alpha$ luminosities are calculated, assuming that Ly$\alpha$ emission is placed at the central wavelength of the narrowbands. The uncertainties of the Ly$\alpha$ luminosities are calculated based on the uncertainties of the NB and $z$ band magnitudes. We obtain the volume number density of LAEs in each Ly$\alpha$ luminosity bin by dividing the number of observed LAEs in each bin by our survey volume (Section 2.1). We correct these number densities for the detection completeness estimated in Section 3.1. The $1\sigma$ uncertainties of the Lya LF measurements are calculated based on Poisson statistics (Gehrels, 1986). Note that we do not include the field-to-field variance in the uncertainties of our Ly$\alpha$ LFs, because the survey areas for $z=5.7$ and $6.6$ LAEs are very large (see Section 2.1). This procedure of Ly$\alpha$ LF derivation is known as the classical method.

We first show our derived Ly$\alpha$ LFs at $z=5.7$ and $z=6.6$ with the classical method in Figure 7. To check field-to-field variations, we present the $z=5.7$ and $z=6.6$ Ly$\alpha$ LF results for the four and five fields in the top and bottom panels, respectively, as well as the results averaged over these fields. We find that our results for these separate fields are consistent with each other, although they have relatively large uncertainties.

In Figure 8, we show our Ly$\alpha$ LF at $z=5.7$ derived with the classical method and previous results. The filled circles represent our $z=5.7$ Ly$\alpha$ LF, which is derived from the HSC SSP data. Our Ly$\alpha$ LF covers a Ly$\alpha$ luminosity range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.9-43.8$. The wide area of the HSC SSP survey allows us to probe this brighter luminosity range than those of previous studies (e.g., Shimasaku et al. (2006); Ouchi et al. (2008); Hu et al. (2010)). We take into account the contamination fractions in our samples (Section 3.2) in the calculations of the Ly$\alpha$ LF uncertainties by increasing the lower $1\sigma$ confidence intervals by $30$%. Similarly, in Figure 9, we show our $z=6.6$ Ly$\alpha$ LF from the HSC SSP data derived with the classical method. The uncertainties from the $f_{\mathrm{cont}}$ value (Section 3.2) are considered. Our $z=6.6$ Ly$\alpha$ LF covers a bright Ly$\alpha$ luminosity range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=43.0-43.8$ thanks to the wide area of the HSC SSP survey. Table 3 shows the values of our Ly$\alpha$ LFs at $z=5.7$ and $z=6.6$.

We fit a Schechter function (Schechter, 1976) to our $z=5.7$ and $z=6.6$ Ly$\alpha$ LFs by minimum $\chi^{2}$ fitting. The Schechter function is defined by

[TABLE]

where $L^{*}_{\mathrm{Ly}\alpha}$ is the characteristic Ly$\alpha$ luminosity, $\phi^{*}_{\mathrm{Ly}\alpha}$ is the normalization, and $\alpha$ is the faint-end slope. We consider two cases. In one case, we use our Ly$\alpha$ LF measurements at $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $<43.5$, where AGN contamination is not significant in lower-$z$ LAE studies (Ouchi et al., 2008; Konno et al., 2016).111 As mentioned in Section 4.1, Shibuya et al. (2017b) have found no clear signature of AGNs for several bright LAEs with $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $>43.5$. Their bright LAEs show narrow Ly$\alpha$ line widths of $<400$ km s*-1* and no clear detection of UV lines such as N v and C iv. However, their investigation is based on the rest-frame UV spectroscopic observations and they cannot rule out the possibility that the bright LAEs host an AGN with faint highly-ionized UV lines (e.g., Hall et al. (2004); Martínez-Sansigre et al. (2006)). In this paper, we present the Ly$\alpha$ LF fitting results for the two cases where we include and exclude the bright-end bins of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $>43.5$ for a conservative discussion.

In the other case, we include the bright-end LF results at $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $\geq 43.5$. In both of these cases, we also use the faint-end Ly$\alpha$ LFs of Ouchi et al. (2008) and Ouchi et al. (2010) for $z=5.7$ and $z=6.6$, respectively. This is because the faint-end Ly$\alpha$ LFs of these studies cover faint Ly$\alpha$ luminosity ranges that we do not reach. Specifically, we include the $z=5.7$ Ly$\alpha$ LF data points of Ouchi et al. (2008) in the range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.4-42.9$ and the $z=6.6$ Ly$\alpha$ LF data points of Ouchi et al. (2010) in the range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.4-43.0$, both of which are not overlapped with the luminosity ranges of our derived LFs. The best-fit Schechter function parameters are listed in Table 4 and the best-fit Schechter functions are shown in Figures 8 and 9 (black thin curve and dashed curve).

The classical method is accurate if the narrowband filter has an ideal boxcar transmission shape. However, the actual narrowband filter transmission shapes are close to a triangle, which causes mainly the following two systematic uncertainties in Ly$\alpha$ LF estimates by the classical method. (I) A Ly$\alpha$ flux of a LAE at a given narrowband magnitude depends on the redshift of the LAE. (II) The minimum EW0 value that corresponds to a given $\textit{BB}-\textit{NB}$ color criterion changes with redshift. These two systematic effects are closely related to each other. Moreover, there are many other systematic uncertainties including the survey volume definitions. We evaluate such systematic uncertainties in our HSC Ly$\alpha$ LFs by carrying out end-to-end Monte Carlo simulations that are conducted in Shimasaku et al. (2006) and Ouchi et al. (2008). We generate a mock catalog of LAEs with a given set of Schechter function parameters ($\phi^{*}_{\mathrm{Ly}\alpha}$, $L^{*}_{\mathrm{Ly}\alpha}$, $\alpha$) and a standard deviation ($\sigma$) of a Gaussian Ly$\alpha$ EW0 probability distribution. LAEs in the mock catalog are uniformly distributed in a comoving volume over the redshift range that a narrowband covers, and their narrowband and broadband magnitudes are measured. We then select LAEs using the same criteria as used for our LAE selections from the actual HSC data. Finally, we derive the surface number densities and color distributions of the selected LAEs, and compare these results with the actual ones (see Shimasaku et al. (2006) and Ouchi et al. (2008) for more details of the simulations). In this comparison, we use the surface number densities and color distributions that are obtained for the $z=5.7$ ($z=6.6$) LAEs in the four (five) fields separately to take into account the different relative depths of these fields. Free parameters in our end-to-end Monte Carlo simulations are $L^{*}_{\mathrm{Ly}\alpha}$ and $\phi^{*}_{\mathrm{Ly}\alpha}$ of the Schechter funtions and $\sigma$ of Gaussian Ly$\alpha$ EW0 probability distributions. The faint-end slope $\alpha$ is fixed at $\alpha=-2.6$ for $z=5.7$ and $\alpha=-2.5$ for $z=6.6$, which are the same as those obtained with the classical method for the Ly$\alpha$ LF measurements in the range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.4-44.0$. Comparing the surface number densities (Figure 3) and color distributions (Figure 5) from the real data with those from the Monte Carlo simulations, we search for the best-fitting set of the three parameters that minimizes $\chi^{2}$. The best-fit Schechter parameters are summarized in Table 4 and examples of the fitting results are shown in Figure 10.

We show the best-fit functions from the Monte Carlo simulations for our Ly$\alpha$ LFs at $z=5.7$ and $6.6$ in Figures 8 and 9, respectively. We find that the best-fit Schechter functions from the simulations are consistent with our HSC Ly$\alpha$ LFs derived by the classical method. Similar conclusions are obtained by Shimasaku et al. (2006) and Ouchi et al. (2008), who have derived the Ly$\alpha$ LFs at $z\sim 3-6$ with Subaru/Suprime-Cam. We confirm that the classical method for the Ly$\alpha$ LF calculations gives a good approximation to the true Ly$\alpha$ LF even in the case of our HSC SSP data. The top panel of Figure 8 compares the luminosities from the classical method ($L_{\rm c}$) and from the simulations ($L_{\rm s}$) at the same number densities as a function of $L_{\rm c}$. We find that the difference between these two luminosities is only $\lesssim$ 0.1 dex. Similarly, the middle panel of Figure 8 shows the ratios of the number densities derived from the classical method to those from the simulations. We find that this ratio is also nearly equal to unity, where the departures of the classical-method data points from the simulation results are smaller than the statistical $\sim 1\sigma$ uncertainties shown with the error bars. Moreover, we also find that the classical-method data points are not always underestimated (Figures 8 vs. 9 ). We thus think that the large correction factors beyond our statistical errors should not be applied to our data points of the classical method, which rather give additional systematics.

As shown in Figures 8 and 9, the best-fit Schechter functions can explain the Ly$\alpha$ LF measurements in the wide luminosity range. If this is true, the faint-end slopes of Ly$\alpha$ LFs are very steep. The best-fit faint-end slope values are $\alpha=-2.5--2.6$ (Table 4), which may indicate that the faint-end slopes of Ly$\alpha$ LFs are steeper than those of the UV LFs at similar redshifts (e.g., Bouwens et al. (2015b)). Note that our best-fit faint-end slopes are steeper than that obtained in previous work on the $z=5.7$ Ly$\alpha$ LF (Dressler et al., 2015).

It should be noted that, if we compare our Ly$\alpha$ LF measurements with the best-fit Schechter function results obtained from the classical method where we consider only the fainter Ly$\alpha$ luminosity range of $\log L_{\rm Ly\alpha}$ [erg s*-1*] $\gtrsim 42.5-43.5$, we find that there is a significant bright-end excess of the $z=5.7$ and $z=6.6$ Ly$\alpha$ LF measurements at $\log L_{\rm Ly\alpha}$ [erg s*-1*] $\gtrsim 43.5$. Based on the deviation of the bright-end data points from the best-fit Schechter function, the significance value of the bright-end excesses is $\simeq 3\sigma$ ($2.6\sigma$ for $z=5.7$ and $3.2\sigma$ for $z=6.6$). For $z=6.6$, similar results are also claimed by some previous studies (e.g., Matthee et al. (2015); Santos et al. (2016); Castellano et al. (2016); Bagley et al. (2017); Zheng et al. (2017)). Although our results suggest that the LF fittings including the bright-end LF results may reveal the true shapes of the Ly$\alpha$ LFs, it is also possible that the bright-end LF results are enhanced by some systematic effects. We discuss possible origins of the bright-end excesses in Section 4.1.

3.5 Comparison with Previous Studies

In this section, we compare our Ly$\alpha$ LFs at $z=5.7$ and $6.6$ with those obtained by previous studies. As shown in Figures 8 and 9, our Ly$\alpha$ LFs are generally consistent with those of the previous results. However, our Ly$\alpha$ LF results do not agree with the high number densities of LAEs recently claimed by Matthee et al. (2015) and Santos et al. (2016). The reason of this discrepancy is unclear. This study and most of the previous studies have derived the Ly$\alpha$ LFs by the classical method and/or by using Monte Carlo simulations that take account of the two systematic uncertainties (I) and (II) in Ly$\alpha$ LF estimates (Section 3.4). Matthee et al. (2015) and Santos et al. (2016) also appear to have considered these two uncertainties; they have applied filter profile correction for Ly$\alpha$ flux estimates and taken into account the incompleteness of the NB-excess color selection. One possible explanation for the discrepancy is that their corrections are redundant, and that the correction factors are overestimated. In fact, in our end-to-end Monte Carlo simulations, we have adopted a Schechter functional form for Ly$\alpha$ LFs and a Gaussian for Ly$\alpha$ EW0 probability distributions, and have determined their best-fit functions simultaneously based on $\chi^{2}$ fitting to the observed surface number densities and the $\textit{BB}-\textit{NB}$ color distributions (Section 3.4). In other words, the two systematic uncertainties are considered at the same time in our simulations. This is because these two systematic effects are closely related to each other. On the other hand, it seems that Matthee et al. (2015) have estimated the effects of the two uncertainties separately in their Sections 4.1 and 4.3 (See also Santos et al. (2016)), which might cause overcorrections due to the redundancy. Another possibility is the difference of the Ly$\alpha$ EW0 distributions. In our simulations, we have adopted a Gaussian Ly$\alpha$ EW0 probability distribution (e.g., Shimasaku et al. (2006); Gronwall et al. (2007); Ouchi et al. (2008)). On the other hand, Matthee et al. (2015) do not describe what functional form is used for the Ly$\alpha$ EW0 distribution in their calculations of the filter profile correction estimates and the color selection incompleteness estimates (see also Santos et al. (2016)). For example, if they assume an EW0 value that is significantly smaller than the typical value for LAEs, they would obtain too large correction factors and thus too large Ly$\alpha$ LF measurements.

4 Discussion

4.1 Systematic Effects in the Ly$\alpha$ LF Measurements

As shown in Section 3.4, our best-fit Schechter functions derived with the end-to-end Monte Carlo simulations as well as the ones derived with the classical method for the Ly$\alpha$ luminosity range of $\log L_{\rm Ly\alpha}$ [erg s*-1*] $\gtrsim 42.5-44.0$ are fitted to the Ly$\alpha$ LF measurements well both at the bright end and fainter magnitude bins. However, the best-fit values of the faint-end slope $\alpha$ are very steep, compared to the shallower slopes of the UV LFs at similar redshifts (e.g., Bouwens et al. (2015b)). Although our results may imply that the wide luminosity range of our Ly$\alpha$ LFs allow us to reveal the true shapes of the Ly$\alpha$ LFs, it is also possible that the bright-end measurements have some systematic effects. There are four possibilities for such systematics. One possibility is the contribution of AGNs, which is the same as the origin of the bright-end excess at $z\sim 2-3$ (e.g., Konno et al. (2016)). Another possibility is the formation of large ionized bubbles in the IGM around bright LAEs during the epoch of reionization (EoR; e.g., Santos et al. (2016); Bagley et al. (2017); Zheng et al. (2017)). The possibility of the gravitational lensing effect also needs to be considered (e.g., Wyithe et al. (2011); Takahashi et al. (2011); Mason et al. (2015)). The other possibility is that merger systems which are blended at ground-based resolution appear as very bright LAEs (e.g., Bowler et al. (2017a)).

Firstly, we discuss the possibility of AGNs. Although the number densities of AGNs rapidly decrease from $z\sim 3$ toward higher redshift (e.g., Haardt & Madau (2012)), some previous studies suggest the existence of (faint) AGNs at $z\sim 6-7$ (e.g., Willott et al. (2010); Mortlock et al. (2011); Kashikawa et al. (2015); Giallongo et al. (2015); Jiang et al. (2016); Bowler et al. (2017b); Parsa et al. (2017) ), which may systematically enhance the bright end of our $z=5.7$ and $6.6$ Ly$\alpha$ LFs. To evaluate this possibility quantitatively, we compare the number densities of faint AGNs presented in the literature with those of bright-end LAEs with $\log L_{\mathrm{Ly}\alpha}$ [erg s*-1*] $>43.5$. The numbers of bright-end LAEs at $z=5.7$ and $6.6$ are 10 and 13, respectively. Dividing the numbers of bright-end LAEs by the survey volumes (Section 2.1), we obtain their number densities of $8.6\times 10^{-7}$ Mpc*-3* and $6.8\times 10^{-7}$ Mpc*-3* at $z=5.7$ and $6.6$, respectively. Since the UV magnitudes of the bright-end LAEs are $M_{\mathrm{UV}}\gtrsim-21$ mag, we compare their number densities with extrapolations of the previous QSO UV LF results for brighter magnitudes (e.g., Willott et al. (2010); Kashikawa et al. (2015); Jiang et al. (2016)). We find that the number densities of bright-end LAEs are consistent with the QSO UV LF results at $z\sim 6$, which indicates that bright-end LAEs with $\log L_{\mathrm{Ly}\alpha}$ [erg s*-1*] $\gtrsim 43.5$ at $z=5.7$ and $6.6$ could be AGNs. It should be noted that our recent deep near-infrared spectroscopic follow-up observations for several bright-end LAEs at $z=5.7$ and $6.6$ reveal no clear signature of AGNs such as a broad Ly$\alpha$ emission line and strong highly-ionized metal lines, e.g., Nv and Civ (Shibuya et al., 2017b). Although these spectroscopy results imply that the observed bright-end LAEs are unlikely to host an AGN, the number of spectroscopically observed bright-end LAEs is still small. To further examine the possibility of AGNs, we will continue to carry out deep follow-up near-infrared spectroscopy.

Secondly, we discuss the possibility of large ionized bubbles. During the EoR, Ly$\alpha$ photons can easily escape into the IGM in the case that the galaxy is surrounded by an ionized bubble which is large enough to allow the Ly$\alpha$ photons to redshift out of resonant scattering before entering the IGM at the edge of the ionized bubble (e.g., Matthee et al. (2015); Bagley et al. (2017)). In this case, it is expected that bright-end LAEs are preferentially observed, which can enhance the number densities of LAEs at the bright end. In other words, the $z=6.6$ bright-end LF may be enhanced by the effect of large ionized bubbles to some extent, although this effect is unlikely to happen at $z=5.7$, where the IGM is already highly ionized (e.g., Fan et al. (2006)). We further consider this possibility speculatively. By using the analytic models of Furlanetto et al. (2006) (See also Furlanetto & Oh (2005)), we quantify the typical size of ionized bubbles around LAEs at $z=6.6$. We use their results of the relations between the globally averaged ionized fraction of the IGM and the typical size of ionized bubbles, where overlaps of ionized bubbles are considered. As we will describe in Section 4.3, we estimate the neutral hydrogen fraction at $z=6.6$ to be $x_{\mathrm{HI}}=0.3\pm 0.2$ from the evolution of the Ly$\alpha$ LFs at $z=5.7-6.6$. Based on the $x_{\mathrm{HI}}$ value and the top panel of Figure 1 of Furlanetto et al. (2006), we obtain the typical size of ionized bubbles at $z=6.6$ of $\sim 15$ comoving Mpc. If the bright-end excess at $z=6.6$ is caused by large ionized bubbles, the sizes of ionized bubbles around bright-end LAEs would be larger than $\sim 15$ comoving Mpc. To estimate the sizes of ionized bubbles around bright-end LAEs, we use the following formula for the Strömgren radius $R_{\mathrm{S}}$ of an ionized bubble around a source at $z=6.6$ by Haiman (2002): $R_{\mathrm{S}}=0.8\times(\mathrm{SFR}/10\ M_{\odot}\ \mathrm{yr}^{-1})^{1/3}(t_{*}/100\ \mathrm{Myr})^{1/3}[(1+z_{*})/7.56]^{-1}$ proper Mpc. In this equation, Haiman (2002) considers an ionizing source at a given redshift $z_{*}$ with a constant SFR and a Salpeter IMF (the $0.1-120M_{\odot}$ mass range), assuming that the source produces ionizing photons during the lifetime ($t_{*}$). From this equation and the UV magnitudes of the bright-end LAEs at $z=6.6$ (i.e., $M_{\mathrm{UV}}\gtrsim-21$ mag), we calculate the size of the ionized bubbles of $R_{\mathrm{S}}\lesssim 7$ comoving Mpc.222 The SFRs can be estimated from UV luminosities with the following equation: $\mathrm{SFR}\ (M_{\odot}\ \mathrm{yr}^{-1})=L_{\mathrm{UV}}\ (\mathrm{erg}\ \mathrm{s}^{-1}\ \mathrm{Hz}^{-1})/(8\times 10^{27})$ (Madau et al., 1998). From this equation, the SFR corresponding to $M_{\mathrm{UV}}=-21$ is $13.6\ M_{\odot}\ \mathrm{yr}^{-1}$. We estimate the ionized bubble size under the assumption that that these bright LAEs have a constant SFR of $13.6\ M_{\odot}\ \mathrm{yr}^{-1}$, and emit ionizing photons during their age of $100$ Myr.

This size is smaller than that estimated from the analytic model of Furlanetto et al. (2006) ($\sim 15$ comoving Mpc). This result implies that, if the bright end of the Ly$\alpha$ LF at $z=6.6$ is enhanced by large ionized bubbles, ionizing sources that are different from the bright LAEs would be clustered around bright LAEs and form large ionized regions by overlapping their ionized bubbles.

Thirdly, we discuss the possibility of the gravitational lensing effect. The lensing effect by foreground massive galaxies boosts apparent magnitudes of LAEs, which can make a bright-end excess of LFs (Wyithe et al. (2011); Takahashi et al. (2011); Mason et al. (2015); Barone-Nugent et al. (2015)). To investigate whether the bright-end LAEs are affected by the gravitational lensing, we identify foreground sources around them which can act as lenses. We check a catalog of massive galaxy clusters that have been found by using the Cluster finding Algorithm based on Multi-band Identification of Red-sequence gAlaxies (CAMIRA; Oguri (2014); Oguri et al. (2017)). In addition, we check the positions of massive ($M_{\rm star}>10^{10.3}M_{\odot}$) red galaxies with photometric redshift of $z_{\rm photo}=0.05-1.05$ (M. Oguri et al. in preparation). However, we find that out of the 23 bright-end LAEs only two have a nearby foreground galaxy on the sky, which may produce modest lensing magnifications of $\mu\sim 1.2-1.7$. Thus, we conclude that the impact of the gravitational lensing on the shapes of the Ly$\alpha$ LFs is small.

Finally, we discuss the possibility of blended merging galaxies. Recently, Bowler et al. (2017a) have found that multi-component systems account for more than $40$% of their bright $z\sim 7$ galaxies based on the analyses of their Hubble images. In fact, our bright-end LAEs include well-studied Himiko and CR7, whose morphologies in the Hubble WFC3 images show possible signatures of galaxy mergers (Ouchi et al. (2013); Sobral et al. (2015)). At least we confirm that the light profiles of our bright-end LAEs in the HSC images are mostly consistent with point sources (Shibuya et al., 2017a). However, the relatively coarse ground-based resolution cannot rule out the possibility that they are merging systems. To examine this possibility, we plan to investigate the morphologies of bright-end LAEs with higher resolution images taken with Hubble.

In summary, the bright end of our Ly$\alpha$ LFs could be systematically enhanced by the contribution of AGNs and/or blended merging galaxies. It may also be possible that large ionized bubbles contribute to the bright end at $z=6.6$ if ionizing sources are clustered around bright-end LAEs. To further investigate the remaining possibilities, follow-up observations are needed.

4.2 Evolution of Ly$\alpha$ LF at $z=5.7-6.6$

We investigate the evolution of the Ly$\alpha$ LF at $z=5.7-6.6$. In Figure 11, we show our Ly$\alpha$ LFs at $z=5.7$ and $6.6$, which are obtained from the $13.8$ deg2 and $21.2$ deg2 sky area of the HSC SSP survey. Here, we show the best-fit Schechter functions for the LF data points in the luminosity range of $\log L_{\rm Ly\alpha}$ [erg s*-1*] $\gtrsim 42.4-44.0$ derived with the classical method, which are good approximations to the true LFs (Section 3.4). We also present the Ly$\alpha$ LF at $z=7.3$ derived by Konno et al. (2014) in this figure, who have conducted the ultradeep $z=7.3$ LAE survey with Subaru/Suprime-Cam. Ouchi et al. (2008) and Ouchi et al. (2010) have derived the Ly$\alpha$ LFs at $z=5.7$ and $6.6$ based on their $\sim 1$ deg2 narrowband imaging data taken with Subaru/Suprime-Cam, and have found the decrease of the Ly$\alpha$ LF from $z=5.7$ to $6.6$. The same results have been obtained by other previous studies (e.g., Kashikawa et al. (2006); Hu et al. (2010); Kashikawa et al. (2011); Santos et al. (2016)). We find such evolution from our Ly$\alpha$ LFs at $z=5.7$ and $6.6$ in Figure 11. To evaluate this evolution at $z=5.7-6.6$ quantitatively, we investigate the error distribution of Schechter parameters. Figure 12 presents the error contours of the Schechter parameters, $L^{*}_{\mathrm{Ly}\alpha}$ and $\phi^{*}_{\mathrm{Ly}\alpha}$, of our $z=5.7$ and $6.6$ Ly$\alpha$ LFs shown with the blue and red ovals, respectively. We also show the error contours for the Ly$\alpha$ LF at $z=7.3$ of Konno et al. (2014). From this figure, the Schechter parameters of the $z=6.6$ Ly$\alpha$ LF are different from those of the $z=5.7$ Ly$\alpha$ LF, and the Ly$\alpha$ LF decreases from $z=5.7$ to $6.6$ at the $>90\%$ confidence level. Note that the evolution of the Ly$\alpha$ LFs that we derive is similar to the one reported by Santos et al. (2016), although our best-fit $L^{*}_{\mathrm{Ly}\alpha}$ values are smaller than theirs. The decreasing trend of the Ly$\alpha$ LFs with increasing redshift obtained in this study is also consistent with those of Ouchi et al. (2010), who have investigated the evolution of LFs in the faint Ly$\alpha$ range ($\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $\lesssim 43$) as shown in Figure 11. It should be noted that the best-fit Ly$\alpha$ LF parameters of $\phi^{*}_{\mathrm{Ly}\alpha}$ and $L^{*}_{\mathrm{Ly}\alpha}$ presented in Figure 12 appear to be shifted from those of Ouchi et al. (2010). This is caused by the difference of the faint-end slope $\alpha$ values. In our Schechter function fitting with the classical method, the slope $\alpha$ is treated as a free parameter and the best-fit value is about $-2.5$. On the other hand, in Ouchi et al. (2010) the $\alpha$ value has been fixed at $-1.5$.

4.3 Estimation of $x_{\mathrm{HI}}$ at $z=6.6$

We estimate the neutral hydrogen fraction, $x_{\mathrm{HI}}$, at $z=6.6$ based on our Ly$\alpha$ LFs at $z=5.7$ and $6.6$ in the same manner as Ouchi et al. (2010) and Konno et al. (2014). We first calculate $T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=6.6}/T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=5.7}$, where $T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z}$ is a Ly$\alpha$ transmission through the IGM at a redshift $z$. The observed Ly$\alpha$ LD, $\rho^{\mathrm{Ly}\alpha}$, can be obtained from

[TABLE]

where $\kappa$ is the conversion factor from UV to Ly$\alpha$ fluxes, $f^{\mathrm{esc}}_{\mathrm{Ly}\alpha}$ is the Ly$\alpha$ escape fraction through the ISM of a galaxy, and $\rho^{\mathrm{UV}}$ is the intrinsic UV LD. Based on the equation, we can estimate the Ly$\alpha$ transmission fraction $T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=6.6}/T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=5.7}$ by

[TABLE]

To calculate $\rho^{\mathrm{Ly}\alpha,\mathrm{tot}}_{z=6.6}/\rho^{\mathrm{Ly}\alpha,\mathrm{tot}}_{z=5.7}$, we use the Ly$\alpha$ LD results in Section 3.4. We adopt the Ly$\alpha$ LDs derived for the Ly$\alpha$ LF measurements in the luminosity range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.4-44.0$, to take account of the contribution from bright-end LAEs as well as from the fainter ones. Based on the UV LF measurements of Bouwens et al. (2015b), $\rho^{\mathrm{UV}}_{z=6.6}/\rho^{\mathrm{UV}}_{z=5.7}=0.74\pm 0.10$ is obtained. Under the assumption of $\kappa_{z=5.7}/\kappa_{z=6.6}=1$ and $f^{\mathrm{esc}}_{\mathrm{Ly}\alpha,z=5.7}/f^{\mathrm{esc}}_{\mathrm{Ly}\alpha,z=6.6}=1$, we obtain $T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=6.6}/T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha,z=5.7}=0.70\pm 0.15$. from Equation (5).

We obtain constraints on $x_{\mathrm{HI}}$ based on comparisons of our results with theoretical models. Santos (2004) have calculated the IGM Ly$\alpha$ transmission fraction as a function of $x_{\mathrm{HI}}$ in two cases of galactic outflow: the Ly$\alpha$ velocity shifts of 0 and 360 km s*-1* from the systemic velocity. It is noted from recent studies that the average velocity shift of Ly$\alpha$ emission is $\sim 200$ km s*-1* for LAEs at $z\sim 2$ (e.g., Hashimoto et al. (2013); Shibuya et al. (2014)). Based on Figure 25 of Santos (2004), our Ly$\alpha$ transmission fraction result is consistent with $x_{\mathrm{HI}}$ $\sim 0.0-0.2$ considering the two cases. Next, we compare our Ly$\alpha$ LF result with the theoretical results of McQuinn et al. (2007), who have derived $z=6.6$ Ly$\alpha$ LFs for various $x_{\mathrm{HI}}$ values based on their radiative transfer simulations. From Figure 4 of McQuinn et al. (2007), we obtain constraints of $x_{\mathrm{HI}}\sim 0.3-0.5$. Finally, we compare our result with a combination of two theoretical models. Dijkstra et al. (2007b) have derived expected Ly$\alpha$ transmission fractions of the IGM as a function of the typical size of ionized bubbles (see also Dijkstra et al. (2007a)). The relation between the typical size of ionized bubbles and $x_{\mathrm{HI}}$ has been calculated by Furlanetto et al. (2006) based on their analytic model. A comparison of our Ly$\alpha$ transmission fraction result with these two models (Figure 6 of Dijkstra et al. (2007b) and the top panel of Figure 1 of Furlanetto et al. (2006)) yields $x_{\mathrm{HI}}\sim 0.1-0.3$. Based on the results described above, we conclude the neutral hydrogen fraction is estimated to be $x_{\mathrm{HI}}=0.1-0.5$, i.e., $x_{\mathrm{HI}}=0.3\pm 0.2$ at $z=6.6$, where the variance of the theoretical model predictions as well as the uncertainties in our Ly$\alpha$ transmission fraction estimates are considered.

Figure 13 shows our $x_{\mathrm{HI}}$ estimate at $z=6.6$ and those taken from the previous studies. The previous results of the $z\gtrsim 7$ Ly$\alpha$ LFs imply $x_{\mathrm{HI}}=0.3-0.8$ at $z=7.3$ (Konno et al., 2014) and $x_{\mathrm{HI}}<0.63$ at $z=7.0$ (Ota et al., 2010). The studies of Ly$\alpha$ emitting fractions indicate $x_{\mathrm{HI}}\gtrsim 0.5$ at $z\sim 7$ (e.g., Pentericci et al. (2011); Schenker et al. (2012); Ono et al. (2012); Treu et al. (2012); Caruana et al. (2012); Caruana et al. (2014); Pentericci et al. (2014); Schenker et al. (2014)). The Ly$\alpha$ damping wing absorption measurements of QSOs suggest $x_{\mathrm{HI}}\gtrsim 0.1$ at $z=7.1$ (Mortlock et al., 2011; Bolton et al., 2011).

As already pointed out in our previous work (Konno et al., 2014), the decrease of the Ly$\alpha$ LF from $z=6.6$ to $7.3$ is larger than that from $z=5.7$ to $6.6$. In Figure 13, this accelerated evolution could be also found, although the uncertainties are large. The Ly$\alpha$ LF evolves from $z=6.6$ to $7.3$ at the $>90$% confidence level, while the difference of $x_{\mathrm{HI}}$ between $z=6.6$ and $7.3$ is only within $1\sigma$. This is because, in our $x_{\mathrm{HI}}$ estimates, we take into account the uncertainties of the UV LFs and the various theoretical model results as well as the uncertainties of the Ly$\alpha$ LFs (see Konno et al. (2014) for details).

Here, we investigate whether the $x_{\mathrm{HI}}$ evolution obtained by our and previous studies can explain the Thomson scattering optical depth, $\tau_{\mathrm{el}}$, value obtained from the latest Planck 2016 data. Because one needs to know $\tau_{\mathrm{el}}$ from a given $x_{\mathrm{HI}}$ evolution, we use the semi-analytic models of Choudhury et al. (2008). They have derived $x_{\mathrm{HI}}$ and $\tau_{\mathrm{el}}$ evolutions by considering three models which differ the minimum halo masses for reionization sources to cover typical scenarios of the cosmic reionization history. These three models are referred to as models A, B, and C corresponding to the minimum halo masses of $\sim 10^{9}$, $\sim 10^{8}$, and $\sim 5\times 10^{5}\ M_{\odot}$, respectively, at $z=6$. We present the $x_{\mathrm{HI}}$ evolutions of the three models in Figure 13, and their $\tau_{\mathrm{el}}$ evolutions in Figure 14. The gray (hatched) region in Figure 14 shows the $1\sigma$ range of $\tau_{\mathrm{el}}$ obtained by Planck (WMAP). The latest results from the Planck observations indicate that the Thomson scattering optical depth is $\tau_{\mathrm{el}}=0.058\pm 0.009$ (Planck Collaboration et al., 2016b), which is significantly lower than the one obtained from the WMAP data. In Figure 13, the models A and B are consistent with our $x_{\mathrm{HI}}$ estimates at $z=6.6$ and $7.3$, and also explain the Thomson scattering optical depth obtained by the latest Planck 2016 data in Figure 14. The model C can barely explain our $x_{\mathrm{HI}}$ value at $z=7.3$, but is placed above the $\tau_{\mathrm{el}}$ of Planck beyond the $1\sigma$ error (Figure 14). Thus, these results show that the cosmic reionization history such as the models A and B can explain both the $x_{\mathrm{HI}}$ estimates and the Planck 2016 $\tau_{\mathrm{el}}$ value simultaneously. Similar conclusions are reached by Robertson et al. (2015) and Bouwens et al. (2015a), who have discussed the UV LF evolution of reionization sources that is independent from our Ly$\alpha$ LF study.

5 Summay

We have derived the Ly$\alpha$ LFs at $z=5.7$ and $6.6$ based on the first-year narrowband and broadband imaging data products obtained by the HSC SSP survey. Our major results are listed below:

Our HSC narrowband images for $z=5.7$ and $6.6$ LAEs have the effective areas of $\sim 13.8$ deg2 and $\sim 21.2$ deg2, respectively. The $5\sigma$ limiting magnitudes of the narrowband images are $\sim 25.0$ mag and $\sim 25.5$ mag in the Deep and UltraDeep layers, respectively. Using these narrowband images, we have identified, in total, $\sim$ 2,000 LAEs at $z=5.7$ and $6.6$ with a bright Ly$\alpha$ luminosity range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $\simeq 42.9-43.8$. Our HSC LAE sample is $\sim 2-6$ times larger than those of previous studies of $z\sim 6-7$ LAEs. 2. 2.

Based on the LAE samples, we have derived the Ly$\alpha$ LFs at $z=5.7$ and $6.6$. We have obtained the best-fit Schechter parameters of $L^{*}_{\mathrm{Ly}\alpha}=1.6^{+2.2}_{-0.6}\times 10^{43}\ \mathrm{erg}\ \mathrm{s}^{-1}$, $\phi^{*}_{\mathrm{Ly}\alpha}=0.85^{+1.87}_{-0.77}\times 10^{-4}\ \mathrm{Mpc}^{-3}$, and $\alpha=-2.6^{+0.6}_{-0.4}$ for the $z=5.7$ Ly$\alpha$ LF, and $L^{*}_{\mathrm{Ly}\alpha}=1.7^{+0.3}_{-0.7}\times 10^{43}\ \mathrm{erg}\ \mathrm{s}^{-1}$, $\phi^{*}_{\mathrm{Ly}\alpha}=0.47^{+1.44}_{-0.44}\times 10^{-4}\ \mathrm{Mpc}^{-3}$, and $\alpha=-2.5^{+0.5}_{-0.5}$ for the $z=6.6$ Ly$\alpha$ LF, if we consider the Ly$\alpha$ luminosity range of $\log L(\mathrm{Ly}\alpha)$ [erg s*-1*] $=42.4-44.0$. 3. 3.

Our Ly$\alpha$ LFs at $z=5.7$ and $z=6.6$ show a very steep faint-end slope, although there is a possibility that the bright-end measurements are enhanced by some systematic effects such as the contribution from AGNs, blended merging galaxies, and/or large ionized bubbles around bright LAEs. 4. 4.

We have confirmed the decrease of the Ly$\alpha$ LF from $z=5.7$ to $6.6$. This evolution is caused by the Ly$\alpha$ damping wing absorption of neutral hydrogen in the IGM. Based on the decrease of the Ly$\alpha$ LF at $z=5.7-6.6$, we have estimated the IGM neutral hydrogen fraction of $x_{\mathrm{HI}}=0.3\pm 0.2$ at $z=6.6$. The $x_{\mathrm{HI}}$ evolution obtained from our and previous studies can explain the Thomson scattering optical depth measurement of the latest Planck 2016.

{ack}

We thank Mamoru Doi, Kentaro Motohara, Toshitaka Kajino, and Masafumi Ishigaki for useful discussion and comments. We appreciate Masayuki Umemura and Masao Mori, who provided the fund for the narrowband filters.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.

Based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center, NAOJ.

A.K. acknowledges support from the Japan Society for the Promotion of Science (JSPS) through the JSPS Research Fellowship for Young Scientists. This work is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (15H02064) Grant-in-Aid for Scientific Research (A) through Japan Society for the Promotion of Science. N.K. acknowledges supports from the JSPS grant 15H03645.