SILVERRUSH. IV. Lyα luminosity functions at z = 5.7 and 6.6 studied with ∼1300 Lyα emitters on the 14–21 deg2 sky

SILVERRUSH. IV. Lyα luminosity functions at z = 5.7 and 6.6 studied with ∼1300 Lyα emitters... Abstract We present the Lyα luminosity functions (LFs) at z = 5.7 and 6.6 derived from a new large sample of 1266 Lyα emitters (LAEs) identified in total areas of 14 and 21 deg2, respectively, based on the early narrowband data of the Subaru/Hyper Suprime-Cam 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α LFs with unprecedentedly small statistical and systematic uncertainties in a wide Lyα luminosity range of 1042.8–43.8 erg s−1. We obtain best-fit Schechter parameters of $$L^{*}_{\mathrm{Ly}\alpha } = 1.6^{+2.2}_{-0.6} \ (1.7^{+0.3}_{-0.7}) \times 10^{43}\:\mathrm{erg}\:\mathrm{s}^{-1}$$, $$\phi ^{*}_{\mathrm{Ly}\alpha } = 0.85^{+1.87}_{-0.77} \ (0.47^{+1.44}_{-0.44}) \times 10^{-4}\:\mathrm{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α LFs are consistent with the majority of the previous studies, but find that our Lyα 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α LFs at z = 5.7 and 6.6 show an indication that the faint-end slope is very steep (α ≃ −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α LF measurements with four independent reionization models, we estimate the neutral hydrogen fraction of the intergalactic medium to be $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ at z = 6.6, which is consistent with the small Thomson scattering optical depth obtained by Planck 2016. 1 Introduction Lyα emission lines are one of the key properties of galaxies for exploring a high-z universe. Lyα emitters (LAEs), which generally have a spectrum of a luminous Lyα 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 a highly ionized state (e.g., Ono et al. 2010a, 2010b; Nakajima & Ouchi 2014; Kusakabe et al. 2015; Inoue et al. 2016). Lyα 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α LFs increase from z ∼ 0 to z ∼ 3, show a moderate plateau between z ∼ 3 and z ∼ 6, and decrease toward z ≳ 6 (e.g., Deharveng et al. 2008; Ouchi et al. 2008; Kashikawa et al. 2011). The evolution of Lyα LFs is different from that of UV LFs, which increase from z ∼ 0 to z ∼ 2, and turn to decrease beyond z ≳ 3 (e.g., Schiminovich et al. 2005; Reddy & Steidel 2009; Bouwens et al. 2015a; see also figure 7 of Konno et al. 2016). The difference in the evolutionary trend between Lyα and UV LFs could be related to the escaping process of Lyα photons not only from the H i interstellar medium (ISM) of a galaxy, but also from the H i intergalactic medium (IGM). The Lyα 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α luminosity densities (LDs) to those estimated from intrinsic UV LDs, largely increases from z ∼ 0 to z ∼ 6 by two orders of magnitude, and turns to decrease beyond z ≳ 6 (e.g., Hayes et al. 2011). The rapid evolution of the Lyα escape fraction from z ∼ 6 to z ∼ 0 can be explained by the combination of Lyα attenuation by dust and the Lyα resonance scattering effect by H i in the ISM. In the case that the ISM H i density of a galaxy is large, the path lengths of Lyα photons become longer due to the resonant scattering, and these Lyα photons are subject to attenuation by dust. Konno et al. (2016) used simple expanding shell models, which compute the Lyα radiative transfer by Monte Carlo simulations (MCLya: Verhamme et al. 2006; Schaerer et al. 2011), and suggested that the large increase of the Lyα escape fraction at z = 0–6 can be reproduced by the combination of an H i column density decrease (by two orders of magnitude) and the average dust extinction values. The decrease of the Lyα LFs at z ≳ 6 is related to the cosmic reionization, because the Lyα damping wing of H i in the IGM attenuates Lyα photons from a galaxy. Previous studies have found that Lyα LFs at z ∼ 7 significantly decrease from those at z ∼ 6 (e.g., Kashikawa et al. 2006; Ouchi et al. 2010; Hu et al. 2010; Santos et al. 2016), and especially at z ≳ 7, Lyα LFs decrease rapidly (e.g., Konno et al. 2014). The neutral hydrogen fraction of IGM, $$x_\mathrm{H\,{\small I}}$$, can be estimated by the Lyα LD evolution subtracting the galaxy evolution effect. Ouchi et al. (2010) have constrained $$x_\mathrm{H\,{\small I}} = 0.2 \pm 0.2$$ at z = 6.6 from the Lyα LF evolution at z = 5.7–6.6 (see also Malhotra & Rhoads 2004; Kashikawa et al. 2006). Similarly, the neutral hydrogen fractions at z ≳ 7 have also been estimated from the Lyα LF evolution (Ota et al. 2010, 2017; Konno et al. 2014). These $$x_\mathrm{H\,{\small I}}$$ estimates could constrain the history of cosmic reionization by comparison with the Thomson scattering optical depth of the cosmic microwave background (CMB). Recently, a large number of wide-field narrowband imaging surveys have been conducted not only to spread the Lyα luminosity ranges of Lyα LFs, but also to reveal physical properties for luminous LAEs. At z ∼ 2–3, luminous LAEs are known to have counterparts in multiwavelength data (e.g., X-ray and radio) and/or extended Lyα 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α LFs at log L(Lyα) [erg s−1] ≳ 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 ∼ 6.6 (e.g., Himiko by Ouchi et al. 2009; CR7 and MASOSA by Sobral et al. 2015; 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. 2018b for CR7). In this paper, we present the Lyα LFs at z = 5.7 and 6.6 based on the Subaru/Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP: Aihara et al. 2018a). 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α luminosity more efficiently than Suprime-Cam. In our HSC SSP survey, a total of ∼13.8 deg2 and ∼21.2 deg2 sky areas are covered by NB816 and NB921 observations, respectively (for details, see also subsection 2.1; Ouchi et al. 2018; Shibuya et al. 2018a). These wide-field HSC NB data sets allow us to determine the Lyα LFs at z = 5.7 and 6.6 with unprecedented accuracy. By examining the evolution of these Lyα LFs at z = 5.7–6.6, we can constrain the $$x_\mathrm{H\,{\small I}}$$ value at z = 6.6. Moreover, based on these HSC SSP data, we can push the Lyα 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α LFs at these redshifts, and compare the Lyα LFs with those of previous studies in section 3. We examine the Lyα 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. 2018), the photometric and spectroscopic properties of LAEs at these redshifts (Shibuya et al. 2018a, 2018b), a systematic survey for LAE overdense regions (R. Higuchi et al. in preparation), and our Lyα LF studies. The other program is the high-z dropout galaxy study, Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH: Ono et al. 2018; Harikane et al. 2018; Toshikawa et al. 2018). Throughout this paper, we use magnitudes in the AB system (Oke 1974). We adopt ΛCDM cosmology with a parameter set of (h, Ωm, $$\Omega _\Lambda$$, σ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 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. 2018; Furusawa et al. 2018; Kawanomoto et al. 2017; Komiyama et al. 2018). 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. 2018a). The HSC SSP has been allocated 300 nights over 5 yr, and started in 2014 March. The HSC SSP survey has three layers with different sets of area and depth: the Wide, Deep, and UltraDeep layers. These layers will cover sky areas of ∼1400 deg2, ∼30 deg2, and ∼4 deg2 with the 5 σ limiting magnitudes (in the r band) of ∼26 mag, ∼27 mag, and ∼28 mag, respectively. While the broadband images are taken in all 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 2014 March to 2016 April (S16A), where all additional data taken in 2016 January to April have been merged with the data of Public Data Release 1 (Aihara et al. 2018b). With the NB816 filter, the HSC SSP survey has observed two blank fields in the Deep layer, the D-DEEP2-3 (23h30m00s, +00°00΄00$${^{\prime\prime}_{.}}$$0) and D-ELAIS-N1 (16h10m00s, +54°00΄00$${^{\prime\prime}_{.}}$$0) fields, and two blank fields in the UltraDeep layer, the UD-COSMOS (10h00m29s, +02°12΄21$${^{\prime\prime}_{.}}$$0) and UD-SXDS (02h18m00s, −05°00΄00$${^{\prime\prime}_{.}}$$0) fields. For the NB921 filter, a blank field of the D-COSMOS (10h00m29s, +02°12΄21$${^{\prime\prime}_{.}}$$0) field in the Deep layer has also been observed in addition to 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 are listed in table 1. Fig. 1. View largeDownload slide Filter response curves for the broadband and narrowband filters of Subaru/HSC. The red lines at the wavelengths of ∼8100 Å and ∼9200 Å are the transmission curves of NB816 and NB921, respectively. The black solid lines denote the curves for the broadband filters (g, r, i, z, and y). These response curves take account of the CCD quantum efficiency of HSC (black dotted line), airmass, the transmittance of the dewar window and primary focus unit, and the reflectivity of the primary mirror. For reference, we also plot the filter response curve of the z΄-band filter of Subaru/Suprime-Cam (blue solid curve). The peaks of these curves are normalized to 1.0 for clarity. The upper x-axis shows the redshift of Lyα. (Color online) Fig. 1. View largeDownload slide Filter response curves for the broadband and narrowband filters of Subaru/HSC. The red lines at the wavelengths of ∼8100 Å and ∼9200 Å are the transmission curves of NB816 and NB921, respectively. The black solid lines denote the curves for the broadband filters (g, r, i, z, and y). These response curves take account of the CCD quantum efficiency of HSC (black dotted line), airmass, the transmittance of the dewar window and primary focus unit, and the reflectivity of the primary mirror. For reference, we also plot the filter response curve of the z΄-band filter of Subaru/Suprime-Cam (blue solid curve). The peaks of these curves are normalized to 1.0 for clarity. The upper x-axis shows the redshift of Lyα. (Color online) Table 1. Summary of HSC/NB816 and NB921 data.* Field  Area (NB816)  Area (NB921)  g†  r†  i†  NB816†  z†  NB921†  y†    (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  —  —  —  —  —  —  —  Field  Area (NB816)  Area (NB921)  g†  r†  i†  NB816†  z†  NB921†  y†    (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  —  —  —  —  —  —  —  *The narrowband and broadband data are obtained in the HSC SSP survey. †The 5 σ limiting magnitude in a circular aperture with a diameter of 1$${^{\prime\prime}_{.}}$$5. View Large The HSC data are reduced by the HSC SSP survey team with hscPipe (Bosch et al. 2018), 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. 2018). 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, 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 ∼1.16 × 107 Mpc3 and ∼1.91 × 107 Mpc3 for z = 5.7 and 6.6 LAEs, respectively. The narrowband images reach the 5 σ limiting magnitudes in a 1$${^{\prime\prime}_{.}}$$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 point-spread function (PSF) sizes of the HSC images are typically <0$${^{\prime\prime}_{.}}$$8, which is sufficiently smaller than the aperture diameter of 1$${^{\prime\prime}_{.}}$$5 (see Aihara et al. 2018b for details). We summarize the 5 σ 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. 2018). 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. 2018b). 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α emission, i − NB816 and z − NB921, respectively, and no detection of blue continuum fluxes. We first select objects with magnitudes brighter than the 5 σ limit in NB816 or NB921 from the HSC SSP database. We then apply similar selection criteria to those of Ouchi et al. (2008, 2010):   \begin{equation} \begin{array}{l}i - {\it NB816} \ge 1.2, \\ g > g_{3 \, \sigma }, \ \mathrm{and} \\ {[} ( r \le r_{3 \, \sigma }\ \mathrm{and} \ r - i \ge 1.0) \ \mathrm{or} \ (r > r_{3 \, \sigma }) ] \end{array} \end{equation} (1)for z = 5.7 LAEs, and   \begin{equation} \begin{array}{l}z - {\it NB921} \ge 1.0, \\ g > g_{3 \, \sigma }, \\ r > r_{3 \, \sigma }, \ \mathrm{and} \\ {[} ( z \le z_{3 \, \sigma }\ \mathrm{and} \ i - z \ge 1.0) \ \mathrm{or} \ (z > z_{3 \, \sigma })] \end{array} \end{equation} (2)for z = 6.6 LAEs, where (g3 σ, r3 σ, and z3 σ) are the 3 σ limiting magnitudes of the (g, r, z) bands. Note that the criterion in the first parentheses of the third criterion in equation (1) and the fourth criterion in equation (2) are used to select bright objects whose spectral energy distribution 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 ≥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 visual inspection to reject spurious sources and cosmic rays, we identify 1081 and 1273 LAE candidates at z = 5.7 and 6.6, respectively (Shibuya et al. 2018a). The samples of these LAE candidates are referred to as the “LAE All” samples. The LAE All samples are ∼2–6 times larger than photometric samples in previous studies (e.g., Ouchi et al. 2008, 2010; Matthee et al. 2015; Santos et al. 2016). This sample is used for clustering analyses in our companion paper (Ouchi et al. 2018). Details of the sample construction, including the color–magnitude diagrams of NB−BB vs. NB, are presented in Shibuya et al. (2018a). In this Lyα 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 - \mathit {NB921}$$ color criterion for z = 6.6 LAEs. The color selection criterion for z = 5.7 LAEs [i.e., i − NB816 > 1.2 in equation (1)] corresponds to the rest-frame Lyα equivalent width (EW), EW0, of EW0 ≳ 10 Å in the case of a flat UV continuum (i.e., fν = const.) with IGM attenuation (Madau 1995). This EW limit is similar to those of previous studies (EW0 ≳ 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α LF results. On the other hand, the color criterion of z − NB921 > 1.0 in equation (2) for z = 6.6 LAEs corresponds to an 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 the z΄ (or z) band filter is different between Suprime-Cam and HSC (figure 1). Specifically, the central wavelength of the HSC z-band filter is ∼160 Å shorter than that of the Suprime-Cam z΄-band filter. For consistency of comparison, we adopt a more stringent color criterion of $$z - \mathit {NB921} > 1.8$$. This criterion corresponds to EW0 > 14 Å (fν = 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α LF” samples. We use the LAE Lyα LF samples to derive surface number densities and color distributions (subsection 3.3), and Lyα LFs at z = 5.7 and 6.6 (subsection 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 twice as large as 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 in the numbers of z = 6.6 LAEs between UD-COSMOS and UD-SXDS. Table 2. Photometric sample of z = 5.7 and 6.6 LAEs. Field  LAE All sample*  LAE Lyα LF sample†  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  Field  LAE All sample*  LAE Lyα LF sample†  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  *The numbers of LAE candidates selected based on the color selection criteria [equations (1) and (2)] and the contamination rejection process (Shibuya et al. 2018a). †The numbers of LAE candidates used in our Lyα LF measurements. For the z = 6.6 sample, we adopt a more stringent z − NB921 color criterion (subsection 2.2). View Large 3 Lyα 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. 2018; R. Murata et al. in preparation). Using the SynPipe software, we distribute ∼18000 pseudo-LAEs with various magnitudes in NB816 and NB921 images. These pseudo-LAEs have a Sérsic profile with Sérsic index of n = 1.5, and a half-light radius of re ∼ 0.9 kpc, which corresponds to 0$${^{\prime\prime}_{.}}$$15 and 0$${^{\prime\prime}_{.}}$$17 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 ∼ 6 Lyman-break galaxies with LUV = 0.3–$$1\, L^{*}_{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 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 ≳80% for bright objects with NB ≲24.5 mag, and ∼40% at the 5 σ limiting magnitudes of these narrowband images. We correct for detection completeness to derive the surface number densities and the Lyα LFs of LAEs in subsections 3.3 and 3.4. For the D-DEEP2-3 field, we use the detection completeness shown in figure 2; for the other fields, we shift it along the magnitude considering the limiting magnitudes of the narrowband images. Fig. 2. View largeDownload slide Detection completeness, fdet, of our NB816 and NB921 images taken with Subaru/HSC. The filled circles in the top and bottom figures represent the completeness values in a magnitude bin of Δm = 0.5 mag as a function of narrowband magnitude in the D-DEEP2-3 field. The 5 σ limiting magnitudes of the NB816 and NB921 images are 25.2 mag and 24.9 mag, respectively. Fig. 2. View largeDownload slide Detection completeness, fdet, of our NB816 and NB921 images taken with Subaru/HSC. The filled circles in the top and bottom figures represent the completeness values in a magnitude bin of Δm = 0.5 mag as a function of narrowband magnitude in the D-DEEP2-3 field. The 5 σ limiting magnitudes of the NB816 and NB921 images are 25.2 mag and 24.9 mag, respectively. 3.2 Contamination In our companion paper Shibuya et al. (2018b), 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, fcont, is estimated to be fcont = 4/28 ≃ 14%(4/53 ≃ 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 [O iii] 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 fcont ≃ 33% (=4/12) and ≃ 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 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 σ confidence intervals of the Lyα LFs by 30% (see subsection 3.4). Note that our estimated fcont values are similar to those obtained by Ouchi et al. (2008, 2010) and Kashikawa et al. (2011), fcont = 0%–30%, who 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 (subsection 2.1). These surface number densities are corrected for detection completeness (subsection 3.1). The 1 σ error bars of the surface number densities are calculated based on 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 σ 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) in figure 3 and Ouchi et al. (2010) in figure 4. These previous studies 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, 2010). Fig. 3. View largeDownload slide Surface number densities of z = 5.7 LAEs. The red filled and open circles represent our surface number densities at z = 5.7 in each field with and without detection completeness correction (subsection 3.1), respectively. The black circles denote the z = 5.7 surface number densities of Ouchi et al. (2008). (Color online) Fig. 3. View largeDownload slide Surface number densities of z = 5.7 LAEs. The red filled and open circles represent our surface number densities at z = 5.7 in each field with and without detection completeness correction (subsection 3.1), respectively. The black circles denote the z = 5.7 surface number densities of Ouchi et al. (2008). (Color online) Fig. 4. View largeDownload slide As figure 3, but for z = 6.6. The black circles denote the z = 6.6 surface number densities of Ouchi et al. (2010). (Color online) Fig. 4. View largeDownload slide As figure 3, but for z = 6.6. The black circles denote the z = 6.6 surface number densities of Ouchi et al. (2010). (Color online) Figures 5 and 6 show the color distributions of $$\mathit {i} - \mathit {NB816}$$ and $$z - \mathit {NB921}$$ for z = 5.7 and 6.6 LAEs, respectively. Magnitudes with a detection significance below 2 σ are replaced with the 2 σ limiting magnitudes. Based on figures 3–6, we estimate the best-fit Schechter functions and Lyα EW0 distributions by Monte Carlo simulations in subsection 3.5. Fig. 5. View largeDownload slide $$i - \mathit {NB816}$$ color distribution of z = 5.7 LAEs in each field. Fig. 5. View largeDownload slide $$i - \mathit {NB816}$$ color distribution of z = 5.7 LAEs in each field. Fig. 6. View largeDownload slide $$z -\mathit {NB921}$$ color distribution of z = 6.6 LAEs in each field. Fig. 6. View largeDownload slide $$z -\mathit {NB921}$$ color distribution of z = 6.6 LAEs in each field. 3.4 Lyα luminosity functions at z = 5.7 and 6.6 We present Lyα LFs at z = 5.7 and 6.6 based on our HSC Lyα LF samples constructed in subsection 2.2. We derive the Lyα LFs in the same manner as Ouchi et al. (2008, 2010). We calculate the Lyα EW0 values of z = 5.7 (6.6) LAEs from the magnitudes of NB816 (NB921) and the z band, and estimate the Lyα 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α line and a flat UV continuum (i.e., fν = const.) 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α luminosities are calculated assuming that Lyα emission is placed at the central wavelength of the narrowbands. The uncertainties of the Lyα 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α luminosity bin by dividing the number of observed LAEs in each bin by our survey volume (subsection 2.1). We correct these number densities for the detection completeness estimated in subsection 3.1. The 1 σ uncertainties of the Lyα 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α LFs, because the survey areas for z = 5.7 and 6.6 LAEs are very large (see subsection 2.1). This procedure of Lyα LF derivation is known as the classical method. We first show our derived Lyα 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α 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. Fig. 7. View largeDownload slide Lyα LFs of LAEs at z = 5.7 (top) and z = 6.6 (bottom). The filled circles represent our results averaged over the separate fields. The open circles, squares, pentagons, triangles, and diamonds are the Lyα LFs of the separate fields of D-DEEP2-3, D-ELAIS-N1, UD-COSMOS, UD-SXDS, and D-COSMOS, respectively. The filled squares are the results of Ouchi et al. (2008, 2010). Fig. 7. View largeDownload slide Lyα LFs of LAEs at z = 5.7 (top) and z = 6.6 (bottom). The filled circles represent our results averaged over the separate fields. The open circles, squares, pentagons, triangles, and diamonds are the Lyα LFs of the separate fields of D-DEEP2-3, D-ELAIS-N1, UD-COSMOS, UD-SXDS, and D-COSMOS, respectively. The filled squares are the results of Ouchi et al. (2008, 2010). In figure 8, we show our Lyα LF at z = 5.7 derived with the classical method and previous results. The filled circles represent our z = 5.7 Lyα LF, which is derived from the HSC SSP data. Our Lyα LF covers a Lyα luminosity range of log L(Lyα) [erg s−1] = 42.9–43.8. The wide area of the HSC SSP survey allows us to probe this luminosity range, which is brighter 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 (subsection 3.2) in the calculations of the Lyα LF uncertainties by increasing the lower 1 σ confidence intervals by 30%. Similarly, in figure 9, we show our z = 6.6 Lyα LF from the HSC SSP data derived with the classical method. The uncertainties from the fcont value (subsection 3.2) are considered. Our z = 6.6 Lyα LF covers a bright Lyα luminosity range of log L(Lyα) [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α LFs at z = 5.7 and z = 6.6. Fig. 8. View largeDownload slide Top: Ratio of the luminosity from the classical method (Lc) to that from the simulations (Ls) at the same number density as a function of Lc based on comparisons of the best-fit LFs. Middle: Ratio of the number density from the classical method (nc) to that from the simulations (ns) at a given luminosity. The filled circles compare the LF data points derived by the classical method (filled circles in the bottom panel) with the best-fit LF derived by the simulations (thick curve in the bottom panel). The solid curve is based on a comparison of the two best-fit LFs obtained by the classical method and the simulations. Bottom: Lyα LFs of z = 5.7 LAEs. The filled circles represent our z = 5.7 Lyα LF results based on the HSC SSP data, where we consider the contamination fraction fcont = 0%–30% in the LF uncertainties. The filled squares denote the Lyα LF given by Ouchi et al. (2008). The best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s studies are shown with the thin (dashed) curve, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is considered. We also show the end-to-end Monte Carlo simulation result, as described in subsection 3.4, with the thick curve. The open diamonds, triangles, pentagons, and crosses are the Lyα LFs of Shimasaku et al. (2006), Murayama et al. (2007), Hu et al. (2010), and Santos, Sobral, and Matthee (2016), respectively. Fig. 8. View largeDownload slide Top: Ratio of the luminosity from the classical method (Lc) to that from the simulations (Ls) at the same number density as a function of Lc based on comparisons of the best-fit LFs. Middle: Ratio of the number density from the classical method (nc) to that from the simulations (ns) at a given luminosity. The filled circles compare the LF data points derived by the classical method (filled circles in the bottom panel) with the best-fit LF derived by the simulations (thick curve in the bottom panel). The solid curve is based on a comparison of the two best-fit LFs obtained by the classical method and the simulations. Bottom: Lyα LFs of z = 5.7 LAEs. The filled circles represent our z = 5.7 Lyα LF results based on the HSC SSP data, where we consider the contamination fraction fcont = 0%–30% in the LF uncertainties. The filled squares denote the Lyα LF given by Ouchi et al. (2008). The best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s studies are shown with the thin (dashed) curve, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is considered. We also show the end-to-end Monte Carlo simulation result, as described in subsection 3.4, with the thick curve. The open diamonds, triangles, pentagons, and crosses are the Lyα LFs of Shimasaku et al. (2006), Murayama et al. (2007), Hu et al. (2010), and Santos, Sobral, and Matthee (2016), respectively. Fig. 9. View largeDownload slide As figure 8, but for z = 6.6. In the bottom panel, the filled squares denote the Lyα LF given by Ouchi et al. (2010). The thin (dashed) curve is the best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s results, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is taken into account. The simulation result is also shown with the thick curve. The open pentagons and crosses are the Lyα LFs of Hu et al. (2010). The small crosses are the Lyα LFs of Matthee et al. (2015) derived for the UDS and COSMOS fields. The large crosses are taken from Santos, Sobral, and Matthee (2016), who have revised the Lyα LFs of Matthee et al. (2015). The open triangles are the Lyα LF obtained by Kashikawa et al. (2011). The open pentagon represents the Lyα LF at z = 6.4 by Bagley et al. (2017), who conducted the HST WFC3 infrared spectroscopic parallel survey. Fig. 9. View largeDownload slide As figure 8, but for z = 6.6. In the bottom panel, the filled squares denote the Lyα LF given by Ouchi et al. (2010). The thin (dashed) curve is the best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s results, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is taken into account. The simulation result is also shown with the thick curve. The open pentagons and crosses are the Lyα LFs of Hu et al. (2010). The small crosses are the Lyα LFs of Matthee et al. (2015) derived for the UDS and COSMOS fields. The large crosses are taken from Santos, Sobral, and Matthee (2016), who have revised the Lyα LFs of Matthee et al. (2015). The open triangles are the Lyα LF obtained by Kashikawa et al. (2011). The open pentagon represents the Lyα LF at z = 6.4 by Bagley et al. (2017), who conducted the HST WFC3 infrared spectroscopic parallel survey. Table 3. Lyα LFs at z = 5.7 and 6.6 from this work. log L(Lyα)*  log ϕ†  (erg s−1)  ([Δlog L(Lyα)]−1 Mpc−3)  z = 5.7  42.95  $${-3.478}^{+0.038}_{-0.193}$$  43.05  $${-3.735}^{+0.044}_{-0.199}$$  43.15  $${-3.953}^{+0.043}_{-0.198}$$  43.25  $${-4.163}^{+0.055}_{-0.210}$$  43.35  $${-4.427}^{+0.076}_{-0.231}$$  43.45  $${-4.970}^{+0.147}_{-0.308}$$  43.55  $${-5.170}^{+0.187}_{-0.355}$$  43.65  $${-5.318}^{+0.224}_{-0.401}$$  43.75  $${-5.717}^{+0.365}_{-0.606}$$  z = 6.6  43.15  $${-4.194}^{+0.154}_{-0.317}$$  43.25  $${-4.407}^{+0.101}_{-0.258}$$  43.35  $${-4.748}^{+0.087}_{-0.243}$$  43.45  $${-5.132}^{+0.140}_{-0.300}$$  43.55  $${-5.433}^{+0.203}_{-0.374}$$  43.65  $${-5.609}^{+0.253}_{-0.438}$$  43.75  $${-6.212}^{+0.519}_{-0.917}$$  43.85  $${-6.226}^{+0.519}_{-0.917}$$  log L(Lyα)*  log ϕ†  (erg s−1)  ([Δlog L(Lyα)]−1 Mpc−3)  z = 5.7  42.95  $${-3.478}^{+0.038}_{-0.193}$$  43.05  $${-3.735}^{+0.044}_{-0.199}$$  43.15  $${-3.953}^{+0.043}_{-0.198}$$  43.25  $${-4.163}^{+0.055}_{-0.210}$$  43.35  $${-4.427}^{+0.076}_{-0.231}$$  43.45  $${-4.970}^{+0.147}_{-0.308}$$  43.55  $${-5.170}^{+0.187}_{-0.355}$$  43.65  $${-5.318}^{+0.224}_{-0.401}$$  43.75  $${-5.717}^{+0.365}_{-0.606}$$  z = 6.6  43.15  $${-4.194}^{+0.154}_{-0.317}$$  43.25  $${-4.407}^{+0.101}_{-0.258}$$  43.35  $${-4.748}^{+0.087}_{-0.243}$$  43.45  $${-5.132}^{+0.140}_{-0.300}$$  43.55  $${-5.433}^{+0.203}_{-0.374}$$  43.65  $${-5.609}^{+0.253}_{-0.438}$$  43.75  $${-6.212}^{+0.519}_{-0.917}$$  43.85  $${-6.226}^{+0.519}_{-0.917}$$  *The luminosity bin of our Lyα LFs at z = 5.7 and 6.6. The bin size is Δlog L(Lyα) = 0.1. †The number densities corrected for the detection completeness (see subsection 3.1). View Large We fit a Schechter function (Schechter 1976) to our z = 5.7 and z = 6.6 Lyα LFs by minimum χ2 fitting. The Schechter function is defined by   \begin{eqnarray} \phi (L_{\mathrm{Ly}\alpha }) dL_{\mathrm{Ly}\alpha } = \phi ^{*}_{\mathrm{Ly}\alpha } \left( \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right)^{\alpha } \exp \left( - \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right) d \left( \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right), \nonumber\\ \end{eqnarray} (3)where $$L^{*}_{\mathrm{Ly}\alpha }$$ is the characteristic Lyα luminosity, $$\phi ^{*}_{\mathrm{Ly}\alpha }$$ is the normalization, and α is the faint-end slope. We consider two cases. In one case, we use our Lyα LF measurements at log L(Lyα) [erg s−1] < 43.5, where AGN contamination is not significant in lower-z LAE studies (Ouchi et al. 2008; Konno et al. 2016).1 In the other case, we include the bright-end LF results at log L(Lyα) [erg s−1] ≥ 43.5. In both of these cases, we also use the faint-end Lyα LFs of Ouchi et al. (2008) for z = 5.7 and Ouchi et al. (2010) for z = 6.6. This is because the faint-end Lyα LFs of these studies cover faint Lyα luminosity ranges that we do not reach. Specifically, we include the z = 5.7 Lyα LF data points of Ouchi et al. (2008) in the range log L(Lyα) [erg s−1] = 42.4–42.9 and the z = 6.6 Lyα LF data points of Ouchi et al. (2010) in the range log L(Lyα) [erg s−1] = 42.4–43.0, both of which do not overlap 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). Table 4. Best-fit Schechter parameters and Lyα luminosity densities. Redshift  L*  ϕ*  α  ρLyαobs†    (1043 erg s−1)  (10−4 Mpc−3)    (1039 erg s−1 Mpc−3)  Classical method for log L(Lyα) [erg s−1] = 42.4–43.5†  5.7  $$1.07^{+0.77}_{-0.38}$$  $$2.46^{+3.48}_{-1.86}$$  $$-2.26^{+0.76}_{-0.44}$$  3.39  6.6  $$0.82^{+0.86}_{-0.30}$$  $$2.83^{+3.52}_{-2.38}$$  $$-1.86^{+0.79}_{-0.67}$$  1.96  Classical method for log L(Lyα) [erg s−1] = 42.4–44.0§  5.7  $$1.64^{+2.16}_{-0.62}$$  $$0.849^{+1.87}_{-0.771}$$  $$-2.56^{+0.53}_{-0.45}$$  3.49  6.6  $$1.66^{+0.30}_{-0.69}$$  $$0.467^{+1.44}_{-0.442}$$  $$-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  Redshift  L*  ϕ*  α  ρLyαobs†    (1043 erg s−1)  (10−4 Mpc−3)    (1039 erg s−1 Mpc−3)  Classical method for log L(Lyα) [erg s−1] = 42.4–43.5†  5.7  $$1.07^{+0.77}_{-0.38}$$  $$2.46^{+3.48}_{-1.86}$$  $$-2.26^{+0.76}_{-0.44}$$  3.39  6.6  $$0.82^{+0.86}_{-0.30}$$  $$2.83^{+3.52}_{-2.38}$$  $$-1.86^{+0.79}_{-0.67}$$  1.96  Classical method for log L(Lyα) [erg s−1] = 42.4–44.0§  5.7  $$1.64^{+2.16}_{-0.62}$$  $$0.849^{+1.87}_{-0.771}$$  $$-2.56^{+0.53}_{-0.45}$$  3.49  6.6  $$1.66^{+0.30}_{-0.69}$$  $$0.467^{+1.44}_{-0.442}$$  $$-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  †The Lyα luminosity densities obtained by integrating the Lyα LF down to log LLyα [erg s−1] = 42.4, which corresponds to $$\sim 0.3 \times L^{*}_{\mathrm{Ly}\alpha }\ (z=3\mbox{--}6)$$. ‡The best-fit parameters are derived with the classical method for the Lyα LF measurements in the luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5. §The best-fit parameters are derived with the classical method for the Lyα LF measurements in the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. View Large 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 the following two systematic uncertainties in Lyα LF estimates by the classical method: (I) A Lyα flux of an LAE at a given narrowband magnitude depends on the redshift of the LAE. (II) The minimum EW0 value that corresponds to a given BB − 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α 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 }$$, α) and the standard deviation (σ) of a Gaussian Lyα EW0 probability distribution. LAEs in the mock catalog are uniformly distributed in a comoving volume over the redshift range that the 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 σ of Gaussian Lyα EW0 probability distributions. The faint-end slope α is fixed at α = −2.6 for z = 5.7 and α = −2.5 for z = 6.6, which are the same as those obtained with the classical method for the Lyα LF measurements in the range log L(Lyα) [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 χ2. The best-fit Schechter parameters are summarized in table 4, and examples of the fitting results are shown in figure 10. Fig. 10. View largeDownload slide Examples of the results of our end-to-end Monte Carlo simulations. Top: Best-fit surface number densities for the D-DEEP2-3 LAEs at z = 5.7 (left) and z = 6.6 (right) are shown with solid lines. The other symbols are as in figures 3 and 4. Bottom: Best-fit BB − NB color distributions of LAEs. In the left and right panels, the solid lines denote the results for our z = 5.7 and z = 6.6 LAEs in the D-DEEP2-3 field, respectively. The other symbols are as in figures 5 and 6. (Color online) Fig. 10. View largeDownload slide Examples of the results of our end-to-end Monte Carlo simulations. Top: Best-fit surface number densities for the D-DEEP2-3 LAEs at z = 5.7 (left) and z = 6.6 (right) are shown with solid lines. The other symbols are as in figures 3 and 4. Bottom: Best-fit BB − NB color distributions of LAEs. In the left and right panels, the solid lines denote the results for our z = 5.7 and z = 6.6 LAEs in the D-DEEP2-3 field, respectively. The other symbols are as in figures 5 and 6. (Color online) We show the best-fit functions from the Monte Carlo simulations for our Lyα 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α LFs derived by the classical method. Similar conclusions are obtained by Shimasaku et al. (2006) and Ouchi et al. (2008), who derived the Lyα LFs at z ∼ 3–6 with Subaru/Suprime-Cam. We confirm that the classical method for the Lyα LF calculations gives a good approximation to the true Lyα LF even in the case of our HSC SSP data. The top panel of figure 8 compares the luminosities from the classical method (Lc) and from the simulations (Ls) at the same number densities as a function of Lc. We find that the difference between these two luminosities is only ≲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 ∼1 σ uncertainties shown with the error bars. Moreover, we also find that the classical-method data points are not always underestimated (figure 8 vs. figure 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α LF measurements in the wide luminosity range. If this is true, the faint-end slopes of Lyα LFs are very steep. The best-fit faint-end slope values are α = −2.5–−2.6 (table 4), which may indicate that the faint-end slopes of Lyα LFs are steeper than those of the UV LFs at similar redshifts (e.g., Bouwens et al. 2015a). Note that our best-fit faint-end slopes are steeper than obtained in previous work on the z = 5.7 Lyα LF (Dressler et al. 2015). It should be noted that if we compare our Lyα LF measurements with the best-fit Schechter function results obtained from the classical method, where we consider only the fainter Lyα luminosity range of log LLyα [erg s−1] ≳ 42.5–43.5, we find that there is a significant bright-end excess of the z = 5.7 and z = 6.6 Lyα LF measurements at log LLyα [erg s−1] ≳ 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 ∼3 σ (2.6 σ for z = 5.7 and 3.2 σ 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α 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 subsection 4.1. 3.5 Comparison with previous studies In this section, we compare our Lyα LFs at z = 5.7 and 6.6 with those obtained by previous studies. As shown in figures 8 and 9, our Lyα LFs are generally consistent with those of the previous results. However, our Lyα LF results do not agree with the high number densities of LAEs recently claimed by Matthee et al. (2015) and Santos, Sobral, and Matthee (2016). The reason for this discrepancy is unclear. This study and most of the previous studies have derived the Lyα 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α LF estimates (subsection 3.4). Matthee et al. (2015) and Santos, Sobral, and Matthee (2016) also appear to have considered these two uncertainties; they have applied filter profile correction for Lyα 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α LFs and a Gaussian for Lyα EW0 probability distributions, and have determined their best-fit functions simultaneously based on χ2 fitting to the observed surface number densities and the BB − NB color distributions (subsection 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) estimated the effects of the two uncertainties separately in their subsections 4.1 and 4.3 (see also Santos et al. 2016), which might cause overcorrections due to the redundancy. Another possibility is the difference in the Lyα EW0 distributions. In our simulations we adopted a Gaussian Lyα 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α EW0 distribution in their calcula-tions 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α LF measurements. 4 Discussion 4.1 Systematic effects in the Lyα LF measurements As shown in subsection 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α luminosity range of log LLyα [erg s−1] ≳ 42.5–44.0, are fitted well to the Lyα LF measurements both at the bright end and fainter magnitude bins. However, the best-fit values of the faint-end slope α are very steep compared to the shallower slopes of the UV LFs at similar redshifts (e.g., Bouwens et al. 2015a). Although our results may imply that the wide luminosity range of our Lyα LFs allow us to reveal the true shapes of the Lyα 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 ∼ 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 that are blended at ground-based resolution appear as very bright LAEs (e.g., Bowler et al. 2017a). First, we discuss the possibility of AGNs. Although the number densities of AGNs rapidly decrease from z ∼ 3 toward higher redshift (e.g., Haardt & Madau 2012), some previous studies suggest the existence of (faint) AGNs at z ∼ 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α 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 LLyα [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 (subsection 2.1), we obtain their number densities of 8.6 × 10−7 Mpc−3 and 6.8 × 10−7 Mpc−3 at z = 5.7 and 6.6, respectively. Since the UV magnitudes of the bright-end LAEs are MUV ≳ −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 ∼ 6, which indicates that bright-end LAEs with log LLyα [erg s−1] ≳ 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α emission line and strong highly-ionized metal lines, e.g., N v and C iv (Shibuya et al. 2018b). 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α photons can easily escape into the IGM in the case that the galaxy is surrounded by an ionized bubble large enough to allow the Lyα 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, Zaldarriaga, and Hernquist (2006)— see also Furlanetto and 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 subsection 4.3, we estimate the neutral hydrogen fraction at z = 6.6 to be $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ from the evolution of the Lyα LFs at z = 5.7–6.6. Based on the $$x_\mathrm{H\,{\small I}}$$ value and the top panel of figure 1 of Furlanetto, Zaldarriaga, and Hernquist (2006), we obtain the typical size of ionized bubbles at z = 6.6 of ∼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 ∼15 comoving Mpc. To estimate the sizes of ionized bubbles around bright-end LAEs, we use the following formula for the Strömgren radius RS of an ionized bubble around a source at z = 6.6 by Haiman (2002): RS = 0.8 × (SFR/10 M⊙ yr−1)1/3(t*/100 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–120 M⊙ 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., MUV ≳ −21 mag), we calculate the size of the ionized bubbles of RS ≲ 7 comoving Mpc.2 This size is smaller than that estimated from the analytic model of Furlanetto, Zaldarriaga, and Hernquist (2006) (∼15 comoving Mpc). This result implies that if the bright end of the Lyα 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 that 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. 2018). In addition, we check the positions of massive (Mstar > 1010.3 M⊙) red galaxies with photometric redshift of zphoto = 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 μ ∼ 1.2–1.7. Thus, we conclude that the impact of gravitational lensing on the shapes of the Lyα LFs is small. Finally, we discuss the possibility of blended merging galaxies. Recently, Bowler et al. (2017a) found that multi-component systems account for more than 40% of their bright z ∼ 7 galaxies based on 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. 2018a). 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α 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α LF at z = 5.7–6.6 We investigate the evolution of the Lyα LF at z = 5.7–6.6. In figure 11, we show our Lyα 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 LLyα [erg s−1] ≳ 42.4–44.0 derived with the classical method, which are good approximations to the true LFs (subsection 3.4). We also present the Lyα LF at z = 7.3 derived by Konno et al. (2014) in this figure; they conducted the ultradeep z = 7.3 LAE survey with Subaru/Suprime-Cam. Ouchi et al. (2008, 2010) derived the Lyα LFs at z = 5.7 and 6.6 based on their ∼1 deg2 narrowband imaging data taken with Subaru/Suprime-Cam, and found a decrease of the Lyα LF from z = 5.7 to 6.6. The same results have been obtained by other previous studies (e.g., Kashikawa et al. 2006, 2011; Hu et al. 2010; Santos et al. 2016). We find such evolution from our Lyα 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α LFs shown with the blue and red ovals, respectively. We also show the error contours for the Lyα LF at z = 7.3 of Konno et al. (2014). From this figure, the Schechter parameters of the z = 6.6 Lyα LF are different from those of the z = 5.7 Lyα LF, and the Lyα LF decreases from z = 5.7 to 6.6 at the >90% confidence level. Note that the evolution of the Lyα LFs that we derive is similar to that reported by Santos, Sobral, and Matthee (2016), although our best-fit $$L^{*}_{\mathrm{Ly}\alpha }$$ values are smaller than theirs. The decreasing trend of the Lyα LFs with increasing redshift obtained in this study is also consistent with those of Ouchi et al. (2010), who investigated the evolution of LFs in the faint Lyα range, log L(Lyα) [erg s−1] ≲ 43, as shown in figure 11. It should be noted that the best-fit Lyα 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 α values. In our Schechter function fitting with the classical method, the slope α 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 α value has been fixed at −1.5. Fig. 11. View largeDownload slide Evolution of the Lyα LFs. The blue and red filled circles are our z = 5.7 and 6.6 Lyα LF measurements, respectively, which are derived from the HSC SSP data. The blue and red filled squares denote the z = 5.7 and 6.6 Lyα LF measurements with the Subaru/Suprime-Cam data given by Ouchi et al. (2008, 2010), respectively. The best-fit Schechter function for the Lyα LF at z = 5.7 (6.6) is shown with the blue (red) solid curve, which are derived for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange circles represent the Lyα LF at z = 7.3 derived by Konno et al. (2014). (Color online) Fig. 11. View largeDownload slide Evolution of the Lyα LFs. The blue and red filled circles are our z = 5.7 and 6.6 Lyα LF measurements, respectively, which are derived from the HSC SSP data. The blue and red filled squares denote the z = 5.7 and 6.6 Lyα LF measurements with the Subaru/Suprime-Cam data given by Ouchi et al. (2008, 2010), respectively. The best-fit Schechter function for the Lyα LF at z = 5.7 (6.6) is shown with the blue (red) solid curve, which are derived for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange circles represent the Lyα LF at z = 7.3 derived by Konno et al. (2014). (Color online) Fig. 12. View largeDownload slide Schechter parameter 1 σ and 2 σ confidence intervals, $$L^{*}_{\mathrm{Ly}\alpha }$$ and $$\phi ^{*}_{\mathrm{Ly}\alpha }$$. The blue (red) contours correspond to z = 5.7 (6.6). The blue and red crosses are the best-fit Schechter parameters for the Lyα LFs at z = 5.7 and 6.6, respectively. These results are obtained with the classical method for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange contours show the results for the z = 7.3 Lyα LF (Konno et al. 2014). (Color online) Fig. 12. View largeDownload slide Schechter parameter 1 σ and 2 σ confidence intervals, $$L^{*}_{\mathrm{Ly}\alpha }$$ and $$\phi ^{*}_{\mathrm{Ly}\alpha }$$. The blue (red) contours correspond to z = 5.7 (6.6). The blue and red crosses are the best-fit Schechter parameters for the Lyα LFs at z = 5.7 and 6.6, respectively. These results are obtained with the classical method for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange contours show the results for the z = 7.3 Lyα LF (Konno et al. 2014). (Color online) 4.3 Estimation of $$x_\mathrm{H\,{\small I}}$$ at z = 6.6 We estimate the neutral hydrogen fraction, $$x_\mathrm{H\,{\small I}}$$, at z = 6.6 based on our Lyα 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 TIGMLyα, z = 6.6/TIGMLyα, z = 5.7, where TIGMLyα, z is the Lyα transmission through the IGM at a redshift z. The observed Lyα LD, ρLyα, can be obtained from   \begin{equation} \rho ^{\mathrm{Ly}\alpha } = \kappa \ T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z} \ f^{\mathrm{esc}}_{\mathrm{Ly}\alpha } \ \rho ^{\mathrm{UV}}, \end{equation} (4)where κ is the conversion factor from UV to Lyα fluxes, fescLyα is the Lyα escape fraction through the ISM of a galaxy, and ρUV is the intrinsic UV LD. Based on the equation, we can estimate the Lyα transmission fraction TIGMLyα, z = 6.6/TIGMLyα, z = 5.7 by   \begin{equation} \frac{T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z = 6.6}}{T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z = 5.7}} = \frac{\kappa _{z = 5.7}}{\kappa _{z = 6.6}} \frac{f^{\mathrm{esc}}_{\mathrm{Ly}\alpha , z = 5.7}}{f^{\mathrm{esc}}_{\mathrm{Ly}\alpha , z = 6.6}} \frac{\rho ^{\mathrm{Ly}\alpha }_{z = 6.6} / \rho ^{\mathrm{Ly}\alpha }_{z = 5.7}}{\rho ^{\mathrm{UV}}_{z = 6.6} / \rho ^{\mathrm{UV}}_{z = 5.7}}. \end{equation} (5)To calculate $$\rho ^{\mathrm{Ly}\alpha , \mathrm{tot}}_{z = 6.6} / \rho ^{\mathrm{Ly}\alpha , \mathrm{tot}}_{z = 5.7}$$, we use the Lyα LD results in subsection 3.4. We adopt the Lyα LDs derived for the Lyα LF measurements in the luminosity range of log L(Lyα) [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. (2015a), ρUVz = 6.6/ρUVz = 5.7 = 0.74 ± 0.10 is obtained. Under the assumption of κz=5.7/κz = 6.6 = 1 and fescLyα, z = 5.7/fescLyα, z = 6.6 = 1, we obtain TIGMLyα, z = 6.6/TIGMLyα, z = 5.7 = 0.70 ± 0.15 from equation (5). We obtain constraints on $$x_\mathrm{H\,{\small I}}$$ based on comparisons of our results with theoretical models. Santos (2004) calculated the IGM Lyα transmission fraction as a function of $$x_\mathrm{H\,{\small I}}$$ in two cases of galactic outflow: Lyα 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α emission is ∼200 km s−1 for LAEs at z ∼ 2 (e.g., Hashimoto et al. 2013; Shibuya et al. 2014). Based on figure 25 of Santos (2004), our Lyα transmission fraction result is consistent with $$x_\mathrm{H\,{\small I}}$$ ∼0.0–0.2, considering the two cases. Next, we compare our Lyα LF result with the theoretical results of McQuinn et al. (2007), who derived z = 6.6 Lyα LFs for various $$x_\mathrm{H\,{\small I}}$$ values based on their radiative transfer simulations. From figure 4 of McQuinn et al. (2007), we obtain constraints of $$x_\mathrm{H\,{\small I}} \sim 0.3$$–0.5. Finally, we compare our result with a combination of two theoretical models. Dijkstra, Wyithe, and Haiman (2007b) derived expected Lyα 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{H\,{\small I}}$$ has been calculated by Furlanetto, Zaldarriaga, and Hernquist (2006) based on their analytic model. A comparison of our Lyα 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{H\,{\small I}} \sim 0.1$$–0.3. Based on the results described above, we conclude that the neutral hydrogen fraction is estimated to be $$x_\mathrm{H\,{\small I}} = 0.1$$–0.5, i.e., $$x_\mathrm{H\,{\small I}} = 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α transmission fraction estimates are considered. Figure 13 shows our $$x_\mathrm{H\,{\small I}}$$ estimate at z = 6.6 and those taken from the previous studies. The previous results of the z ≳ 7 Lyα LFs imply $$x_\mathrm{H\,{\small I}} = 0.3$$–0.8 at z = 7.3 (Konno et al. 2014) and $$x_\mathrm{H\,{\small I}} < 0.63$$ at z = 7.0 (Ota et al. 2010). The studies of Lyα emitting fractions indicate $$x_\mathrm{H\,{\small I}} \gtrsim 0.5$$ at z ∼ 7 (e.g., Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The Lyα damping wing absorption measurements of QSOs suggest $$x_\mathrm{H\,{\small I}} \gtrsim 0.1$$ at z = 7.1 (Mortlock et al. 2011; Bolton et al. 2011). Fig. 13. View largeDownload slide Evolution of the IGM neutral hydrogen fraction. The top and bottom panels are the same, except for the scales of the ordinate axes (top: linear; bottom: logarithmic). The red filled circles show the $$x_\mathrm{H\,{\small I}}$$ estimates from the Lyα LFs at z = 6.6 (this study) and 7.3 (Konno et al. 2014). The blue filled triangle, square, diamond, and pentagon are the $$x_\mathrm{H\,{\small I}}$$ estimates based on the evolution of the Lyα LF obtained by Malhotra and Rhoads (2004), Kashikawa et al. (2011), Ouchi et al. (2010), and Ota et al. (2010), respectively. The blue open diamond and circle are the constraints on $$x_\mathrm{H\,{\small I}}$$ from the clustering analyses of LAEs (Ouchi et al. 2010) and the Lyα emitting galaxy fraction (Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The previous results from the GRB optical afterglow spectrum analyses are shown with magenta filled triangles (Totani et al. 2006, 2014). The green filled squares and open triangle are the results from the GP test of QSOs (Fan et al. 2006) and the size of QSO near zones (Mortlock et al. 2011; Bolton et al. 2011), respectively. The gray and hatched regions are the 1 σ confidence intervals for the instantaneous reionization redshifts obtained by Planck (Planck Collaboration 2016b) and nine-year WMAP (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The models A, B, and C of Choudhury, Ferrara, and Gallerani (2008) are shown with dotted, dashed, and solid lines, respectively. (Color online) Fig. 13. View largeDownload slide Evolution of the IGM neutral hydrogen fraction. The top and bottom panels are the same, except for the scales of the ordinate axes (top: linear; bottom: logarithmic). The red filled circles show the $$x_\mathrm{H\,{\small I}}$$ estimates from the Lyα LFs at z = 6.6 (this study) and 7.3 (Konno et al. 2014). The blue filled triangle, square, diamond, and pentagon are the $$x_\mathrm{H\,{\small I}}$$ estimates based on the evolution of the Lyα LF obtained by Malhotra and Rhoads (2004), Kashikawa et al. (2011), Ouchi et al. (2010), and Ota et al. (2010), respectively. The blue open diamond and circle are the constraints on $$x_\mathrm{H\,{\small I}}$$ from the clustering analyses of LAEs (Ouchi et al. 2010) and the Lyα emitting galaxy fraction (Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The previous results from the GRB optical afterglow spectrum analyses are shown with magenta filled triangles (Totani et al. 2006, 2014). The green filled squares and open triangle are the results from the GP test of QSOs (Fan et al. 2006) and the size of QSO near zones (Mortlock et al. 2011; Bolton et al. 2011), respectively. The gray and hatched regions are the 1 σ confidence intervals for the instantaneous reionization redshifts obtained by Planck (Planck Collaboration 2016b) and nine-year WMAP (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The models A, B, and C of Choudhury, Ferrara, and Gallerani (2008) are shown with dotted, dashed, and solid lines, respectively. (Color online) As already pointed out in our previous work (Konno et al. 2014), the decrease of the Lyα LF from z = 6.6 to 7.3 is larger than that from z = 5.7 to 6.6. This accelerated evolution could be also found in figure 13, although the uncertainties are large. The Lyα LF evolves from z = 6.6 to 7.3 at the >90% confidence level, while the difference of $$x_\mathrm{H\,{\small I}}$$ between z = 6.6 and 7.3 is only within 1 σ. This is because, in our $$x_\mathrm{H\,{\small I}}$$ 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α LFs (see Konno et al. 2014 for details). Here, we investigate whether the $$x_\mathrm{H\,{\small I}}$$ evolution obtained by this and previous studies can explain the Thomson scattering optical depth, τel, value obtained from the latest Planck 2016 data. Because one needs to know τel from a given $$x_\mathrm{H\,{\small I}}$$ evolution, we use the semi-analytic models of Choudhury, Ferrara, and Gallerani (2008). They derived $$x_\mathrm{H\,{\small I}}$$ and τel evolutions by considering three models which differ in the minimum halo masses for reionization sources to cover typical scenarios for cosmic reionization history. These three models are referred to as models A, B, and C, corresponding to the minimum halo masses of ∼109, ∼108, and ∼5 × 105 M⊙, respectively, at z = 6. We present the $$x_\mathrm{H\,{\small I}}$$ evolutions of the three models in figure 13, and their τel evolutions in figure 14. The gray (hatched) region in figure 14 shows the 1 σ range of τel obtained by Planck (WMAP). The latest results from the Planck observations indicate that the Thomson scattering optical depth is τel = 0.058 ± 0.009 (Planck Collaboration 2016b), which is significantly lower than that obtained from the WMAP data. In figure 13, models A and B are consistent with our $$x_\mathrm{H\,{\small I}}$$ 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. Model C can barely explain our $$x_\mathrm{H\,{\small I}}$$ value at z = 7.3, but is placed above the τel of Planck beyond the 1 σ error (figure 14). Thus, these results show that cosmic reionization histories such as models A and B can explain both the $$x_\mathrm{H\,{\small I}}$$ estimates and the Planck 2016 τel value simultaneously. Similar conclusions are reached by Robertson et al. (2015) and Bouwens et al. (2015b), who discussed the UV LF evolution of reionization sources independent of our Lyα LF study. Fig. 14. View largeDownload slide Thomson scattering optical depth, τel, as a function of redshift. The gray and hatched regions correspond to the 1 σ confidence intervals of the τel measurements obtained by the Planck 2016 data (Planck Collaboration 2016b) and the nine-year WMAP data (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The dotted, dashed, and solid curves are models A, B, and C of Choudhury, Ferrara, and Gallerani (2008), respectively. Fig. 14. View largeDownload slide Thomson scattering optical depth, τel, as a function of redshift. The gray and hatched regions correspond to the 1 σ confidence intervals of the τel measurements obtained by the Planck 2016 data (Planck Collaboration 2016b) and the nine-year WMAP data (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The dotted, dashed, and solid curves are models A, B, and C of Choudhury, Ferrara, and Gallerani (2008), respectively. 5 Summary We have derived the Lyα 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 ∼13.8 deg2 and ∼21.2 deg2, respectively. The 5 σ limiting magnitudes of the narrowband images are ∼25.0 mag and ∼25.5 mag in the Deep and UltraDeep layers, respectively. Using these narrowband images, we have identified, in total, ∼2000 LAEs at z = 5.7 and 6.6 with a bright Lyα luminosity range of log L(Lyα) [erg s−1] ≃ 42.9–43.8. Our HSC LAE sample is ∼2–6 times larger than those of previous studies of z ∼ 6–7 LAEs. Based on the LAE samples, we have derived the Lyα 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α 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α LF, if we consider the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. Our Lyα 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. We have confirmed the decrease of the Lyα LF from z = 5.7 to 6.6. This evolution is caused by the Lyα damping wing absorption of neutral hydrogen in the IGM. Based on the decrease of the Lyα LF at z = 5.7–6.6, we have estimated the IGM neutral hydrogen fraction of $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ at z = 6.6. The $$x_\mathrm{H\,{\small I}}$$ evolution obtained from this and previous studies can explain the Thomson scattering optical depth measurement of the latest Planck 2016. Acknowledgements We thank Mamoru Doi, Kentaro Motohara, Toshitaka Kajino, and Masafumi Ishigaki for useful discussion and comments. We acknowledge 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 the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), the 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, 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 the Japan Society for the Promotion of Science. N.K. acknowledges support from JSPS grant 15H03645. Footnotes 1 As mentioned in subsection 4.1, Shibuya et al. (2018b) found no clear signature of AGNs for several bright LAEs with log L(Lyα) [erg s−1] > 43.5. Their bright LAEs show narrow Lyα line widths of <400 km s−1 and no clear detection of UV lines such as N v and C iv. However, their investigation was based on rest-frame UV spectroscopic observations and they could not 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α LF fit results for the two cases where we include and exclude the bright-end bins of log L(Lyα) [erg s−1] > 43.5 for a conservative discussion. 2 The SFRs can be estimated from UV luminosities with the following equation (Madau et al. 1998): SFR(M⊙ yr−1) = LUV (erg s−1 Hz−1)/(8 × 1027). From this equation, the SFR corresponding to MUV = −21 is 13.6 M⊙ yr−1. We estimate the ionized bubble size under the assumption that that these bright LAEs have a constant SFR of 13.6 M⊙ yr−1, and emit ionizing photons during their age of 100 Myr. References Adams J. J. et al.   2011, ApJS , 192, 5 CrossRef Search ADS   Aihara H. et al.   2018a, PASJ , 70, S4 Aihara H. et al.   2018b, PASJ , 70, S8 Ajiki M. et al.   2002, ApJ , 576, L25 CrossRef Search ADS   Axelrod T., Kantor J., Lupton R. H., Pierfederici F. 2010, in Proc. SPIE, 7740, Software and Cyberinfrastructure for Astronomy , ed. Radziwill N. M., Bridger A. ( Bellingham, WA: SPIE), 774015 Bagley M. B. et al.   2017, ApJ , 837, 11 CrossRef Search ADS   Barone-Nugent R. L., Wyithe J. S. B., Trenti M., Treu T., Oesch P., Bouwens R., Illingworth G. D., Schmidt K. B. 2015, MNRAS , 450, 1224 CrossRef Search ADS   Bennett C. L. et al.   2013, ApJS , 208, 20 CrossRef Search ADS   Bolton J. S., Haehnelt M. G., Warren S. J., Hewett P. C., Mortlock D. J., Venemans B. P., McMahon R. G., Simpson C. 2011, MNRAS , 416, L70 CrossRef Search ADS   Bosch J. et al.   2018, PASJ , 70, S5 Bouwens R. J. et al.   2015a, ApJ , 803, 34 CrossRef Search ADS   Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B., Smit R., Wilkins S. 2015b, ApJ , 811, 140 CrossRef Search ADS   Bowler R. A. A., Dunlop J. S., McLure R. J., McLeod D. J. 2017a, MNRAS , 466, 3612 CrossRef Search ADS   Bowler R. A. A., McLure R. J., Dunlop J. S., McLeod D. J., Stanway E. R., Eldridge J. J., Jarvis M. J. 2017b, MNRAS , 469, 448 CrossRef Search ADS   Cai Z. et al.   2017, ApJ , 837, 71 CrossRef Search ADS   Cantalupo S., Arrigoni-Battaia F., Prochaska J. X., Hennawi J. F., Madau P. 2014, Nature , 506, 63 CrossRef Search ADS PubMed  Caruana J., Bunker A. J., Wilkins S. M., Stanway E. R., Lacy M., Jarvis M. J., Lorenzoni S., Hickey S. 2012, MNRAS , 427, 3055 CrossRef Search ADS   Caruana J., Bunker A. J., Wilkins S. M., Stanway E. R., Lorenzoni S., Jarvis M. J., Ebert H. 2014, MNRAS , 443, 2831 CrossRef Search ADS   Cassata P. et al.   2015, A&A , 573, A24 CrossRef Search ADS   Castellano M. et al.   2016, ApJ , 818, L3 CrossRef Search ADS   Choudhury T. R., Ferrara A., Gallerani S. 2008, MNRAS , 385, L58 CrossRef Search ADS   Coupon J., Czakon N., Bosch J., Komiyama Y., Medezinski E., Miyazaki S., Oguri M. 2018, PASJ , 70, S7 Cowie L. L., Hu E. M. 1998, AJ , 115, 1319 CrossRef Search ADS   Dayal P., Maselli A., Ferrara A. 2011, MNRAS , 410, 830 CrossRef Search ADS   Deharveng J.-M. et al.   2008, ApJ , 680, 1072 CrossRef Search ADS   Dijkstra M., Lidz A., Wyithe J. S. B. 2007a, MNRAS , 377, 1175 CrossRef Search ADS   Dijkstra M., Wyithe J. S. B., Haiman Z. 2007b, MNRAS , 379, 253 CrossRef Search ADS   Dressler A. et al.   2011, PASP , 123, 288 CrossRef Search ADS   Dressler A., Henry A., Martin C. L., Sawicki M., McCarthy P., Villaneuva E. 2015, ApJ , 806, 19 CrossRef Search ADS   Fan X. et al.   2006, AJ , 132, 117 CrossRef Search ADS   Finkelstein S. L. et al.   2013, Nature , 502, 524 CrossRef Search ADS PubMed  Furlanetto S. R., Oh S. P. 2005, MNRAS , 363, 1031 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L. 2006, MNRAS , 365, 1012 CrossRef Search ADS   Furusawa H. et al.   2018, PASJ , 70, S3 Gehrels N. 1986, ApJ , 303, 336 CrossRef Search ADS   Giallongo E. et al.   2015, A&A , 578, A83 CrossRef Search ADS   Gronwall C. et al.   2007, ApJ , 667, 79 CrossRef Search ADS   Guaita L. et al.   2010, ApJ , 714, 255 CrossRef Search ADS   Haardt F., Madau P. 2012, ApJ , 746, 125 CrossRef Search ADS   Haiman Z. 2002, ApJ , 576, L1 CrossRef Search ADS   Haiman Z., Spaans M. 1999, ApJ , 518, 138 CrossRef Search ADS   Hall P. B. et al.   2004, AJ , 127, 3146 CrossRef Search ADS   Harikane Y. et al.   2018, PASJ , 70, S11 Hashimoto T., Ouchi M., Shimasaku K., Ono Y., Nakajima K., Rauch M., Lee J., Okamura S. 2013, ApJ , 765, 70 CrossRef Search ADS   Hayashino T. et al.   2004, AJ , 128, 2073 CrossRef Search ADS   Hayes M., Schaerer D., Östlin G., Mas-Hesse J. M., Atek H., Kunth D. 2011, ApJ , 730, 8 CrossRef Search ADS   Hinshaw G. et al.   2013, ApJS , 208, 19 CrossRef Search ADS   Hu E. M., Cowie L. L., Barger A. J., Capak P., Kakazu Y., Trouille L. 2010, ApJ , 725, 394 CrossRef Search ADS   Hu E. M., Cowie L. L., McMahon R. G. 1998, ApJ , 502, L99 CrossRef Search ADS   Hu E. M., Cowie L. L., Songaila A., Barger A. J., Rosenwasser B., Wold I. G. B. 2016, ApJ , 825, L7 CrossRef Search ADS   Huang S. et al.   2018, PASJ , 70, S6 Inoue A. K. et al.   2016, Science , 352, 1559 CrossRef Search ADS PubMed  Ivezic Z. et al.   2008, arXiv:0805.2366 Iye M. et al.   2006, Nature , 443, 186 CrossRef Search ADS PubMed  Jiang L. et al.   2016, ApJ , 833, 222 CrossRef Search ADS   Jurić M. et al.   2015, arXiv:1512.07914 Kashikawa N. et al.   2002, PASJ , 54, 819 CrossRef Search ADS   Kashikawa N. et al.   2006, ApJ , 648, 7 CrossRef Search ADS   Kashikawa N. et al.   2011, ApJ , 734, 119 CrossRef Search ADS   Kashikawa N. et al.   2015, ApJ , 798, 28 CrossRef Search ADS   Kawanomoto S. et al.   2017, PASJ , submitted Kobayashi M. A. R., Totani T., Nagashima M. 2007, ApJ , 670, 919 CrossRef Search ADS   Komiyama Y. et al.   2018, PASJ , 70, S2 Konno A. et al.   2014, ApJ , 797, 16 CrossRef Search ADS   Konno A. et al.   2016, ApJ , 823, 20 CrossRef Search ADS   Kusakabe H., Shimasaku K., Nakajima K., Ouchi M. 2015, ApJ , 800, L29 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S. 2007, MNRAS , 381, 75 CrossRef Search ADS   Madau P. 1995, ApJ , 441, 18 CrossRef Search ADS   Madau P., Pozzetti L., Dickinson M. 1998, ApJ , 498, 106 CrossRef Search ADS   Magnier E. A. et al.   2013, ApJS , 205, 20 CrossRef Search ADS   Malhotra S., Rhoads J. E. 2002, ApJ , 565, L71 CrossRef Search ADS   Malhotra S., Rhoads J. E. 2004, ApJ , 617, L5 CrossRef Search ADS   Mao J., Lapi A., Granato G. L., de Zotti G., Danese L. 2007, ApJ , 667, 655 CrossRef Search ADS   Martínez-Sansigre A., Rawlings S., Lacy M., Fadda D., Jarvis M. J., Marleau F. R., Simpson C., Willott C. J. 2006, MNRAS , 370, 1479 CrossRef Search ADS   Mason C. A. et al.   2015, ApJ , 805, 79 CrossRef Search ADS   Matsuda Y. et al.   2004, AJ , 128, 569 CrossRef Search ADS   Matthee J., Sobral D., Santos S., Röttgering H., Darvish B., Mobasher B. 2015, MNRAS , 451, 400 CrossRef Search ADS   Mesinger A., Furlanetto S. R. 2008, MNRAS , 386, 1990 CrossRef Search ADS   Miyazaki S. et al.   2012, in Proc. SPIE, 8446, Ground-based and Airborne Instrumentation for Astronomy IV , ed. McLean I. S. et al.   ( Bellingham, WA: SPIE), 84460Z Miyazaki S. et al.   2018, PASJ , 70, S1 Mortlock D. J. et al.   2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murayama T. et al.   2007, ApJS , 172, 523 CrossRef Search ADS   Nakajima K., Ouchi M. 2014, MNRAS , 442, 900 CrossRef Search ADS   Oesch P. A. et al.   2015, ApJ , 804, L30 CrossRef Search ADS   Oguri M. 2014, MNRAS , 444, 147 CrossRef Search ADS   Oguri M. et al.   2018, PASJ , 70, S20 Oke J. B. 1974, ApJS , 27, 21 CrossRef Search ADS   Ono Y. et al.   2010a, MNRAS , 402, 1580 CrossRef Search ADS   Ono Y. et al.   2012, ApJ , 744, 83 CrossRef Search ADS   Ono Y. et al.   2018, PASJ , 70, S10 Ono Y., Ouchi M., Shimasaku K., Dunlop J., Farrah D., McLure R., Okamura S. 2010b, ApJ , 724, 1524 CrossRef Search ADS   Ota K. et al.   2010, ApJ , 722, 803 CrossRef Search ADS   Ota K. et al.   2017, ApJ , 844, 85 CrossRef Search ADS   Ouchi M. et al.   2003, ApJ , 582, 60 CrossRef Search ADS   Ouchi M. et al.   2008, ApJS , 176, 301 CrossRef Search ADS   Ouchi M. et al.   2009, ApJ , 696, 1164 CrossRef Search ADS   Ouchi M. et al.   2010, ApJ , 723, 869 CrossRef Search ADS   Ouchi M. et al.   2013, ApJ , 778, 102 CrossRef Search ADS   Ouchi M. et al.   2018, PASJ , 70, S13 Pacucci F., Pallottini A., Ferrara A., Gallerani S. 2017, MNRAS , 468, L77 CrossRef Search ADS   Parsa S., Dunlop J. S., McLure R. J. 2017, arXiv:1704.07750 Pentericci L. et al.   2011, ApJ , 743, 132 CrossRef Search ADS   Pentericci L. et al.   2014, ApJ , 793, 113 CrossRef Search ADS   Planck Collaboration 2016a, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration 2016b, A&A , 596, A107 CrossRef Search ADS   Reddy N. A., Steidel C. C. 2009, ApJ , 692, 778 CrossRef Search ADS   Rhoads J. E., Malhotra S., Dey A., Stern D., Spinrad H., Jannuzi B. T. 2000, ApJ , 545, L85 CrossRef Search ADS   Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S. 2015, ApJ , 802, L19 CrossRef Search ADS   Santos M. R. 2004, MNRAS , 349, 1137 CrossRef Search ADS   Santos S., Sobral D., Matthee J. 2016, MNRAS , 463, 1678 CrossRef Search ADS   Schaerer D., Hayes M., Verhamme A., Teyssier R. 2011, A&A , 531, A12 CrossRef Search ADS   Schechter P. 1976, ApJ , 203, 297 CrossRef Search ADS   Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P. 2014, ApJ , 795, 20 CrossRef Search ADS   Schenker M. A., Stark D. P., Ellis R. S., Robertson B. E., Dunlop J. S., McLure R. J., Kneib J.-P., Richard J. 2012, ApJ , 744, 179 CrossRef Search ADS   Schiminovich D. et al.   2005, ApJ , 619, L47 CrossRef Search ADS   Schlafly E. F. et al.   2012, ApJ , 756, 158 CrossRef Search ADS   Schlegel D. J., Finkbeiner D. P., Davis M. 1998, ApJ , 500, 525 CrossRef Search ADS   Shibuya T. et al.   2014, ApJ , 788, 74 CrossRef Search ADS   Shibuya T. et al.   2018a, PASJ , 70, S14 Shibuya T. et al.   2018b, PASJ , 70, S15 Shibuya T., Kashikawa N., Ota K., Iye M., Ouchi M., Furusawa H., Shimasaku K., Hattori T. 2012, ApJ , 752, 114 CrossRef Search ADS   Shibuya T., Ouchi M., Harikane Y. 2015, ApJS , 219, 15 CrossRef Search ADS   Shimasaku K. et al.   2006, PASJ , 58, 313 CrossRef Search ADS   Sobral D., Matthee J., Darvish B., Schaerer D., Mobasher B., Röttgering H. J. A., Santos S., Hemmati S. 2015, ApJ , 808, 139 CrossRef Search ADS   Song M., Finkelstein S. L., Livermore R. C., Capak P. L., Dickinson M., Fontana A. 2016, ApJ , 826, 113 CrossRef Search ADS   Stark D. P. et al.   2017, MNRAS , 464, 469 CrossRef Search ADS   Steidel C. C., Adelberger K. L., Shapley A. E., Pettini M., Dickinson M., Giavalisco M. 2000, ApJ , 532, 170 CrossRef Search ADS   Takahashi R., Oguri M., Sato M., Hamana T. 2011, ApJ , 742, 15 CrossRef Search ADS   Taniguchi Y. et al.   2005, PASJ , 57, 165 CrossRef Search ADS   Tonry J. L. et al.   2012, ApJ , 750, 99 CrossRef Search ADS   Toshikawa J. et al.   2018, PASJ , 70, S12 Totani T. et al.   2014, PASJ , 66, 63 CrossRef Search ADS   Totani T., Kawai N., Kosugi G., Aoki K., Yamada T., Iye M., Ohta K., Hattori T. 2006, PASJ , 58, 485 CrossRef Search ADS   Treu T., Trenti M., Stiavelli M., Auger M. W., Bradley L. D. 2012, ApJ , 747, 27 CrossRef Search ADS   Verhamme A., Schaerer D., Maselli A. 2006, A&A , 460, 397 CrossRef Search ADS   Willott C. J. et al.   2010, AJ , 139, 906 CrossRef Search ADS   Wyithe J. S. B., Yan H., Windhorst R. A., Mao S. 2011, Nature , 469, 181 CrossRef Search ADS PubMed  Yamada T., Nakamura Y., Matsuda Y., Hayashino T., Yamauchi R., Morimoto N., Kousai K., Umemura M. 2012, AJ , 143, 79 CrossRef Search ADS   Zabl J., Nørgaard-Nielsen H. U., Fynbo J. P. U., Laursen P., Ouchi M., Kjærgaard P. 2015, MNRAS , 451, 2050 CrossRef Search ADS   Zheng Z.-Y. et al.   2017, ApJ , 842, L22 CrossRef Search ADS   Zitrin A. et al.   2015, ApJ , 810, L12 CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of the Astronomical Society of Japan. All rights reserved. For Permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Publications of the Astronomical Society of Japan Oxford University Press

Loading next page...
 
/lp/ou_press/silverrush-iv-ly-luminosity-functions-at-z-5-7-and-6-6-studied-with-vTNvct7w1x
Publisher
Oxford University Press
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of the Astronomical Society of Japan. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
0004-6264
eISSN
2053-051X
D.O.I.
10.1093/pasj/psx131
Publisher site
See Article on Publisher Site

Abstract

Abstract We present the Lyα luminosity functions (LFs) at z = 5.7 and 6.6 derived from a new large sample of 1266 Lyα emitters (LAEs) identified in total areas of 14 and 21 deg2, respectively, based on the early narrowband data of the Subaru/Hyper Suprime-Cam 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α LFs with unprecedentedly small statistical and systematic uncertainties in a wide Lyα luminosity range of 1042.8–43.8 erg s−1. We obtain best-fit Schechter parameters of $$L^{*}_{\mathrm{Ly}\alpha } = 1.6^{+2.2}_{-0.6} \ (1.7^{+0.3}_{-0.7}) \times 10^{43}\:\mathrm{erg}\:\mathrm{s}^{-1}$$, $$\phi ^{*}_{\mathrm{Ly}\alpha } = 0.85^{+1.87}_{-0.77} \ (0.47^{+1.44}_{-0.44}) \times 10^{-4}\:\mathrm{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α LFs are consistent with the majority of the previous studies, but find that our Lyα 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α LFs at z = 5.7 and 6.6 show an indication that the faint-end slope is very steep (α ≃ −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α LF measurements with four independent reionization models, we estimate the neutral hydrogen fraction of the intergalactic medium to be $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ at z = 6.6, which is consistent with the small Thomson scattering optical depth obtained by Planck 2016. 1 Introduction Lyα emission lines are one of the key properties of galaxies for exploring a high-z universe. Lyα emitters (LAEs), which generally have a spectrum of a luminous Lyα 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 a highly ionized state (e.g., Ono et al. 2010a, 2010b; Nakajima & Ouchi 2014; Kusakabe et al. 2015; Inoue et al. 2016). Lyα 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α LFs increase from z ∼ 0 to z ∼ 3, show a moderate plateau between z ∼ 3 and z ∼ 6, and decrease toward z ≳ 6 (e.g., Deharveng et al. 2008; Ouchi et al. 2008; Kashikawa et al. 2011). The evolution of Lyα LFs is different from that of UV LFs, which increase from z ∼ 0 to z ∼ 2, and turn to decrease beyond z ≳ 3 (e.g., Schiminovich et al. 2005; Reddy & Steidel 2009; Bouwens et al. 2015a; see also figure 7 of Konno et al. 2016). The difference in the evolutionary trend between Lyα and UV LFs could be related to the escaping process of Lyα photons not only from the H i interstellar medium (ISM) of a galaxy, but also from the H i intergalactic medium (IGM). The Lyα 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α luminosity densities (LDs) to those estimated from intrinsic UV LDs, largely increases from z ∼ 0 to z ∼ 6 by two orders of magnitude, and turns to decrease beyond z ≳ 6 (e.g., Hayes et al. 2011). The rapid evolution of the Lyα escape fraction from z ∼ 6 to z ∼ 0 can be explained by the combination of Lyα attenuation by dust and the Lyα resonance scattering effect by H i in the ISM. In the case that the ISM H i density of a galaxy is large, the path lengths of Lyα photons become longer due to the resonant scattering, and these Lyα photons are subject to attenuation by dust. Konno et al. (2016) used simple expanding shell models, which compute the Lyα radiative transfer by Monte Carlo simulations (MCLya: Verhamme et al. 2006; Schaerer et al. 2011), and suggested that the large increase of the Lyα escape fraction at z = 0–6 can be reproduced by the combination of an H i column density decrease (by two orders of magnitude) and the average dust extinction values. The decrease of the Lyα LFs at z ≳ 6 is related to the cosmic reionization, because the Lyα damping wing of H i in the IGM attenuates Lyα photons from a galaxy. Previous studies have found that Lyα LFs at z ∼ 7 significantly decrease from those at z ∼ 6 (e.g., Kashikawa et al. 2006; Ouchi et al. 2010; Hu et al. 2010; Santos et al. 2016), and especially at z ≳ 7, Lyα LFs decrease rapidly (e.g., Konno et al. 2014). The neutral hydrogen fraction of IGM, $$x_\mathrm{H\,{\small I}}$$, can be estimated by the Lyα LD evolution subtracting the galaxy evolution effect. Ouchi et al. (2010) have constrained $$x_\mathrm{H\,{\small I}} = 0.2 \pm 0.2$$ at z = 6.6 from the Lyα LF evolution at z = 5.7–6.6 (see also Malhotra & Rhoads 2004; Kashikawa et al. 2006). Similarly, the neutral hydrogen fractions at z ≳ 7 have also been estimated from the Lyα LF evolution (Ota et al. 2010, 2017; Konno et al. 2014). These $$x_\mathrm{H\,{\small I}}$$ estimates could constrain the history of cosmic reionization by comparison with the Thomson scattering optical depth of the cosmic microwave background (CMB). Recently, a large number of wide-field narrowband imaging surveys have been conducted not only to spread the Lyα luminosity ranges of Lyα LFs, but also to reveal physical properties for luminous LAEs. At z ∼ 2–3, luminous LAEs are known to have counterparts in multiwavelength data (e.g., X-ray and radio) and/or extended Lyα 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α LFs at log L(Lyα) [erg s−1] ≳ 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 ∼ 6.6 (e.g., Himiko by Ouchi et al. 2009; CR7 and MASOSA by Sobral et al. 2015; 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. 2018b for CR7). In this paper, we present the Lyα LFs at z = 5.7 and 6.6 based on the Subaru/Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP: Aihara et al. 2018a). 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α luminosity more efficiently than Suprime-Cam. In our HSC SSP survey, a total of ∼13.8 deg2 and ∼21.2 deg2 sky areas are covered by NB816 and NB921 observations, respectively (for details, see also subsection 2.1; Ouchi et al. 2018; Shibuya et al. 2018a). These wide-field HSC NB data sets allow us to determine the Lyα LFs at z = 5.7 and 6.6 with unprecedented accuracy. By examining the evolution of these Lyα LFs at z = 5.7–6.6, we can constrain the $$x_\mathrm{H\,{\small I}}$$ value at z = 6.6. Moreover, based on these HSC SSP data, we can push the Lyα 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α LFs at these redshifts, and compare the Lyα LFs with those of previous studies in section 3. We examine the Lyα 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. 2018), the photometric and spectroscopic properties of LAEs at these redshifts (Shibuya et al. 2018a, 2018b), a systematic survey for LAE overdense regions (R. Higuchi et al. in preparation), and our Lyα LF studies. The other program is the high-z dropout galaxy study, Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH: Ono et al. 2018; Harikane et al. 2018; Toshikawa et al. 2018). Throughout this paper, we use magnitudes in the AB system (Oke 1974). We adopt ΛCDM cosmology with a parameter set of (h, Ωm, $$\Omega _\Lambda$$, σ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 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. 2018; Furusawa et al. 2018; Kawanomoto et al. 2017; Komiyama et al. 2018). 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. 2018a). The HSC SSP has been allocated 300 nights over 5 yr, and started in 2014 March. The HSC SSP survey has three layers with different sets of area and depth: the Wide, Deep, and UltraDeep layers. These layers will cover sky areas of ∼1400 deg2, ∼30 deg2, and ∼4 deg2 with the 5 σ limiting magnitudes (in the r band) of ∼26 mag, ∼27 mag, and ∼28 mag, respectively. While the broadband images are taken in all 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 2014 March to 2016 April (S16A), where all additional data taken in 2016 January to April have been merged with the data of Public Data Release 1 (Aihara et al. 2018b). With the NB816 filter, the HSC SSP survey has observed two blank fields in the Deep layer, the D-DEEP2-3 (23h30m00s, +00°00΄00$${^{\prime\prime}_{.}}$$0) and D-ELAIS-N1 (16h10m00s, +54°00΄00$${^{\prime\prime}_{.}}$$0) fields, and two blank fields in the UltraDeep layer, the UD-COSMOS (10h00m29s, +02°12΄21$${^{\prime\prime}_{.}}$$0) and UD-SXDS (02h18m00s, −05°00΄00$${^{\prime\prime}_{.}}$$0) fields. For the NB921 filter, a blank field of the D-COSMOS (10h00m29s, +02°12΄21$${^{\prime\prime}_{.}}$$0) field in the Deep layer has also been observed in addition to 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 are listed in table 1. Fig. 1. View largeDownload slide Filter response curves for the broadband and narrowband filters of Subaru/HSC. The red lines at the wavelengths of ∼8100 Å and ∼9200 Å are the transmission curves of NB816 and NB921, respectively. The black solid lines denote the curves for the broadband filters (g, r, i, z, and y). These response curves take account of the CCD quantum efficiency of HSC (black dotted line), airmass, the transmittance of the dewar window and primary focus unit, and the reflectivity of the primary mirror. For reference, we also plot the filter response curve of the z΄-band filter of Subaru/Suprime-Cam (blue solid curve). The peaks of these curves are normalized to 1.0 for clarity. The upper x-axis shows the redshift of Lyα. (Color online) Fig. 1. View largeDownload slide Filter response curves for the broadband and narrowband filters of Subaru/HSC. The red lines at the wavelengths of ∼8100 Å and ∼9200 Å are the transmission curves of NB816 and NB921, respectively. The black solid lines denote the curves for the broadband filters (g, r, i, z, and y). These response curves take account of the CCD quantum efficiency of HSC (black dotted line), airmass, the transmittance of the dewar window and primary focus unit, and the reflectivity of the primary mirror. For reference, we also plot the filter response curve of the z΄-band filter of Subaru/Suprime-Cam (blue solid curve). The peaks of these curves are normalized to 1.0 for clarity. The upper x-axis shows the redshift of Lyα. (Color online) Table 1. Summary of HSC/NB816 and NB921 data.* Field  Area (NB816)  Area (NB921)  g†  r†  i†  NB816†  z†  NB921†  y†    (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  —  —  —  —  —  —  —  Field  Area (NB816)  Area (NB921)  g†  r†  i†  NB816†  z†  NB921†  y†    (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  —  —  —  —  —  —  —  *The narrowband and broadband data are obtained in the HSC SSP survey. †The 5 σ limiting magnitude in a circular aperture with a diameter of 1$${^{\prime\prime}_{.}}$$5. View Large The HSC data are reduced by the HSC SSP survey team with hscPipe (Bosch et al. 2018), 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. 2018). 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, 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 ∼1.16 × 107 Mpc3 and ∼1.91 × 107 Mpc3 for z = 5.7 and 6.6 LAEs, respectively. The narrowband images reach the 5 σ limiting magnitudes in a 1$${^{\prime\prime}_{.}}$$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 point-spread function (PSF) sizes of the HSC images are typically <0$${^{\prime\prime}_{.}}$$8, which is sufficiently smaller than the aperture diameter of 1$${^{\prime\prime}_{.}}$$5 (see Aihara et al. 2018b for details). We summarize the 5 σ 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. 2018). 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. 2018b). 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α emission, i − NB816 and z − NB921, respectively, and no detection of blue continuum fluxes. We first select objects with magnitudes brighter than the 5 σ limit in NB816 or NB921 from the HSC SSP database. We then apply similar selection criteria to those of Ouchi et al. (2008, 2010):   \begin{equation} \begin{array}{l}i - {\it NB816} \ge 1.2, \\ g > g_{3 \, \sigma }, \ \mathrm{and} \\ {[} ( r \le r_{3 \, \sigma }\ \mathrm{and} \ r - i \ge 1.0) \ \mathrm{or} \ (r > r_{3 \, \sigma }) ] \end{array} \end{equation} (1)for z = 5.7 LAEs, and   \begin{equation} \begin{array}{l}z - {\it NB921} \ge 1.0, \\ g > g_{3 \, \sigma }, \\ r > r_{3 \, \sigma }, \ \mathrm{and} \\ {[} ( z \le z_{3 \, \sigma }\ \mathrm{and} \ i - z \ge 1.0) \ \mathrm{or} \ (z > z_{3 \, \sigma })] \end{array} \end{equation} (2)for z = 6.6 LAEs, where (g3 σ, r3 σ, and z3 σ) are the 3 σ limiting magnitudes of the (g, r, z) bands. Note that the criterion in the first parentheses of the third criterion in equation (1) and the fourth criterion in equation (2) are used to select bright objects whose spectral energy distribution 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 ≥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 visual inspection to reject spurious sources and cosmic rays, we identify 1081 and 1273 LAE candidates at z = 5.7 and 6.6, respectively (Shibuya et al. 2018a). The samples of these LAE candidates are referred to as the “LAE All” samples. The LAE All samples are ∼2–6 times larger than photometric samples in previous studies (e.g., Ouchi et al. 2008, 2010; Matthee et al. 2015; Santos et al. 2016). This sample is used for clustering analyses in our companion paper (Ouchi et al. 2018). Details of the sample construction, including the color–magnitude diagrams of NB−BB vs. NB, are presented in Shibuya et al. (2018a). In this Lyα 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 - \mathit {NB921}$$ color criterion for z = 6.6 LAEs. The color selection criterion for z = 5.7 LAEs [i.e., i − NB816 > 1.2 in equation (1)] corresponds to the rest-frame Lyα equivalent width (EW), EW0, of EW0 ≳ 10 Å in the case of a flat UV continuum (i.e., fν = const.) with IGM attenuation (Madau 1995). This EW limit is similar to those of previous studies (EW0 ≳ 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α LF results. On the other hand, the color criterion of z − NB921 > 1.0 in equation (2) for z = 6.6 LAEs corresponds to an 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 the z΄ (or z) band filter is different between Suprime-Cam and HSC (figure 1). Specifically, the central wavelength of the HSC z-band filter is ∼160 Å shorter than that of the Suprime-Cam z΄-band filter. For consistency of comparison, we adopt a more stringent color criterion of $$z - \mathit {NB921} > 1.8$$. This criterion corresponds to EW0 > 14 Å (fν = 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α LF” samples. We use the LAE Lyα LF samples to derive surface number densities and color distributions (subsection 3.3), and Lyα LFs at z = 5.7 and 6.6 (subsection 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 twice as large as 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 in the numbers of z = 6.6 LAEs between UD-COSMOS and UD-SXDS. Table 2. Photometric sample of z = 5.7 and 6.6 LAEs. Field  LAE All sample*  LAE Lyα LF sample†  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  Field  LAE All sample*  LAE Lyα LF sample†  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  *The numbers of LAE candidates selected based on the color selection criteria [equations (1) and (2)] and the contamination rejection process (Shibuya et al. 2018a). †The numbers of LAE candidates used in our Lyα LF measurements. For the z = 6.6 sample, we adopt a more stringent z − NB921 color criterion (subsection 2.2). View Large 3 Lyα 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. 2018; R. Murata et al. in preparation). Using the SynPipe software, we distribute ∼18000 pseudo-LAEs with various magnitudes in NB816 and NB921 images. These pseudo-LAEs have a Sérsic profile with Sérsic index of n = 1.5, and a half-light radius of re ∼ 0.9 kpc, which corresponds to 0$${^{\prime\prime}_{.}}$$15 and 0$${^{\prime\prime}_{.}}$$17 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 ∼ 6 Lyman-break galaxies with LUV = 0.3–$$1\, L^{*}_{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 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 ≳80% for bright objects with NB ≲24.5 mag, and ∼40% at the 5 σ limiting magnitudes of these narrowband images. We correct for detection completeness to derive the surface number densities and the Lyα LFs of LAEs in subsections 3.3 and 3.4. For the D-DEEP2-3 field, we use the detection completeness shown in figure 2; for the other fields, we shift it along the magnitude considering the limiting magnitudes of the narrowband images. Fig. 2. View largeDownload slide Detection completeness, fdet, of our NB816 and NB921 images taken with Subaru/HSC. The filled circles in the top and bottom figures represent the completeness values in a magnitude bin of Δm = 0.5 mag as a function of narrowband magnitude in the D-DEEP2-3 field. The 5 σ limiting magnitudes of the NB816 and NB921 images are 25.2 mag and 24.9 mag, respectively. Fig. 2. View largeDownload slide Detection completeness, fdet, of our NB816 and NB921 images taken with Subaru/HSC. The filled circles in the top and bottom figures represent the completeness values in a magnitude bin of Δm = 0.5 mag as a function of narrowband magnitude in the D-DEEP2-3 field. The 5 σ limiting magnitudes of the NB816 and NB921 images are 25.2 mag and 24.9 mag, respectively. 3.2 Contamination In our companion paper Shibuya et al. (2018b), 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, fcont, is estimated to be fcont = 4/28 ≃ 14%(4/53 ≃ 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 [O iii] 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 fcont ≃ 33% (=4/12) and ≃ 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 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 σ confidence intervals of the Lyα LFs by 30% (see subsection 3.4). Note that our estimated fcont values are similar to those obtained by Ouchi et al. (2008, 2010) and Kashikawa et al. (2011), fcont = 0%–30%, who 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 (subsection 2.1). These surface number densities are corrected for detection completeness (subsection 3.1). The 1 σ error bars of the surface number densities are calculated based on 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 σ 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) in figure 3 and Ouchi et al. (2010) in figure 4. These previous studies 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, 2010). Fig. 3. View largeDownload slide Surface number densities of z = 5.7 LAEs. The red filled and open circles represent our surface number densities at z = 5.7 in each field with and without detection completeness correction (subsection 3.1), respectively. The black circles denote the z = 5.7 surface number densities of Ouchi et al. (2008). (Color online) Fig. 3. View largeDownload slide Surface number densities of z = 5.7 LAEs. The red filled and open circles represent our surface number densities at z = 5.7 in each field with and without detection completeness correction (subsection 3.1), respectively. The black circles denote the z = 5.7 surface number densities of Ouchi et al. (2008). (Color online) Fig. 4. View largeDownload slide As figure 3, but for z = 6.6. The black circles denote the z = 6.6 surface number densities of Ouchi et al. (2010). (Color online) Fig. 4. View largeDownload slide As figure 3, but for z = 6.6. The black circles denote the z = 6.6 surface number densities of Ouchi et al. (2010). (Color online) Figures 5 and 6 show the color distributions of $$\mathit {i} - \mathit {NB816}$$ and $$z - \mathit {NB921}$$ for z = 5.7 and 6.6 LAEs, respectively. Magnitudes with a detection significance below 2 σ are replaced with the 2 σ limiting magnitudes. Based on figures 3–6, we estimate the best-fit Schechter functions and Lyα EW0 distributions by Monte Carlo simulations in subsection 3.5. Fig. 5. View largeDownload slide $$i - \mathit {NB816}$$ color distribution of z = 5.7 LAEs in each field. Fig. 5. View largeDownload slide $$i - \mathit {NB816}$$ color distribution of z = 5.7 LAEs in each field. Fig. 6. View largeDownload slide $$z -\mathit {NB921}$$ color distribution of z = 6.6 LAEs in each field. Fig. 6. View largeDownload slide $$z -\mathit {NB921}$$ color distribution of z = 6.6 LAEs in each field. 3.4 Lyα luminosity functions at z = 5.7 and 6.6 We present Lyα LFs at z = 5.7 and 6.6 based on our HSC Lyα LF samples constructed in subsection 2.2. We derive the Lyα LFs in the same manner as Ouchi et al. (2008, 2010). We calculate the Lyα EW0 values of z = 5.7 (6.6) LAEs from the magnitudes of NB816 (NB921) and the z band, and estimate the Lyα 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α line and a flat UV continuum (i.e., fν = const.) 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α luminosities are calculated assuming that Lyα emission is placed at the central wavelength of the narrowbands. The uncertainties of the Lyα 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α luminosity bin by dividing the number of observed LAEs in each bin by our survey volume (subsection 2.1). We correct these number densities for the detection completeness estimated in subsection 3.1. The 1 σ uncertainties of the Lyα 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α LFs, because the survey areas for z = 5.7 and 6.6 LAEs are very large (see subsection 2.1). This procedure of Lyα LF derivation is known as the classical method. We first show our derived Lyα 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α 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. Fig. 7. View largeDownload slide Lyα LFs of LAEs at z = 5.7 (top) and z = 6.6 (bottom). The filled circles represent our results averaged over the separate fields. The open circles, squares, pentagons, triangles, and diamonds are the Lyα LFs of the separate fields of D-DEEP2-3, D-ELAIS-N1, UD-COSMOS, UD-SXDS, and D-COSMOS, respectively. The filled squares are the results of Ouchi et al. (2008, 2010). Fig. 7. View largeDownload slide Lyα LFs of LAEs at z = 5.7 (top) and z = 6.6 (bottom). The filled circles represent our results averaged over the separate fields. The open circles, squares, pentagons, triangles, and diamonds are the Lyα LFs of the separate fields of D-DEEP2-3, D-ELAIS-N1, UD-COSMOS, UD-SXDS, and D-COSMOS, respectively. The filled squares are the results of Ouchi et al. (2008, 2010). In figure 8, we show our Lyα LF at z = 5.7 derived with the classical method and previous results. The filled circles represent our z = 5.7 Lyα LF, which is derived from the HSC SSP data. Our Lyα LF covers a Lyα luminosity range of log L(Lyα) [erg s−1] = 42.9–43.8. The wide area of the HSC SSP survey allows us to probe this luminosity range, which is brighter 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 (subsection 3.2) in the calculations of the Lyα LF uncertainties by increasing the lower 1 σ confidence intervals by 30%. Similarly, in figure 9, we show our z = 6.6 Lyα LF from the HSC SSP data derived with the classical method. The uncertainties from the fcont value (subsection 3.2) are considered. Our z = 6.6 Lyα LF covers a bright Lyα luminosity range of log L(Lyα) [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α LFs at z = 5.7 and z = 6.6. Fig. 8. View largeDownload slide Top: Ratio of the luminosity from the classical method (Lc) to that from the simulations (Ls) at the same number density as a function of Lc based on comparisons of the best-fit LFs. Middle: Ratio of the number density from the classical method (nc) to that from the simulations (ns) at a given luminosity. The filled circles compare the LF data points derived by the classical method (filled circles in the bottom panel) with the best-fit LF derived by the simulations (thick curve in the bottom panel). The solid curve is based on a comparison of the two best-fit LFs obtained by the classical method and the simulations. Bottom: Lyα LFs of z = 5.7 LAEs. The filled circles represent our z = 5.7 Lyα LF results based on the HSC SSP data, where we consider the contamination fraction fcont = 0%–30% in the LF uncertainties. The filled squares denote the Lyα LF given by Ouchi et al. (2008). The best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s studies are shown with the thin (dashed) curve, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is considered. We also show the end-to-end Monte Carlo simulation result, as described in subsection 3.4, with the thick curve. The open diamonds, triangles, pentagons, and crosses are the Lyα LFs of Shimasaku et al. (2006), Murayama et al. (2007), Hu et al. (2010), and Santos, Sobral, and Matthee (2016), respectively. Fig. 8. View largeDownload slide Top: Ratio of the luminosity from the classical method (Lc) to that from the simulations (Ls) at the same number density as a function of Lc based on comparisons of the best-fit LFs. Middle: Ratio of the number density from the classical method (nc) to that from the simulations (ns) at a given luminosity. The filled circles compare the LF data points derived by the classical method (filled circles in the bottom panel) with the best-fit LF derived by the simulations (thick curve in the bottom panel). The solid curve is based on a comparison of the two best-fit LFs obtained by the classical method and the simulations. Bottom: Lyα LFs of z = 5.7 LAEs. The filled circles represent our z = 5.7 Lyα LF results based on the HSC SSP data, where we consider the contamination fraction fcont = 0%–30% in the LF uncertainties. The filled squares denote the Lyα LF given by Ouchi et al. (2008). The best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s studies are shown with the thin (dashed) curve, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is considered. We also show the end-to-end Monte Carlo simulation result, as described in subsection 3.4, with the thick curve. The open diamonds, triangles, pentagons, and crosses are the Lyα LFs of Shimasaku et al. (2006), Murayama et al. (2007), Hu et al. (2010), and Santos, Sobral, and Matthee (2016), respectively. Fig. 9. View largeDownload slide As figure 8, but for z = 6.6. In the bottom panel, the filled squares denote the Lyα LF given by Ouchi et al. (2010). The thin (dashed) curve is the best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s results, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is taken into account. The simulation result is also shown with the thick curve. The open pentagons and crosses are the Lyα LFs of Hu et al. (2010). The small crosses are the Lyα LFs of Matthee et al. (2015) derived for the UDS and COSMOS fields. The large crosses are taken from Santos, Sobral, and Matthee (2016), who have revised the Lyα LFs of Matthee et al. (2015). The open triangles are the Lyα LF obtained by Kashikawa et al. (2011). The open pentagon represents the Lyα LF at z = 6.4 by Bagley et al. (2017), who conducted the HST WFC3 infrared spectroscopic parallel survey. Fig. 9. View largeDownload slide As figure 8, but for z = 6.6. In the bottom panel, the filled squares denote the Lyα LF given by Ouchi et al. (2010). The thin (dashed) curve is the best-fit Schechter function for the Lyα LFs of our and Ouchi et al.’s results, where the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5 (44.0) is taken into account. The simulation result is also shown with the thick curve. The open pentagons and crosses are the Lyα LFs of Hu et al. (2010). The small crosses are the Lyα LFs of Matthee et al. (2015) derived for the UDS and COSMOS fields. The large crosses are taken from Santos, Sobral, and Matthee (2016), who have revised the Lyα LFs of Matthee et al. (2015). The open triangles are the Lyα LF obtained by Kashikawa et al. (2011). The open pentagon represents the Lyα LF at z = 6.4 by Bagley et al. (2017), who conducted the HST WFC3 infrared spectroscopic parallel survey. Table 3. Lyα LFs at z = 5.7 and 6.6 from this work. log L(Lyα)*  log ϕ†  (erg s−1)  ([Δlog L(Lyα)]−1 Mpc−3)  z = 5.7  42.95  $${-3.478}^{+0.038}_{-0.193}$$  43.05  $${-3.735}^{+0.044}_{-0.199}$$  43.15  $${-3.953}^{+0.043}_{-0.198}$$  43.25  $${-4.163}^{+0.055}_{-0.210}$$  43.35  $${-4.427}^{+0.076}_{-0.231}$$  43.45  $${-4.970}^{+0.147}_{-0.308}$$  43.55  $${-5.170}^{+0.187}_{-0.355}$$  43.65  $${-5.318}^{+0.224}_{-0.401}$$  43.75  $${-5.717}^{+0.365}_{-0.606}$$  z = 6.6  43.15  $${-4.194}^{+0.154}_{-0.317}$$  43.25  $${-4.407}^{+0.101}_{-0.258}$$  43.35  $${-4.748}^{+0.087}_{-0.243}$$  43.45  $${-5.132}^{+0.140}_{-0.300}$$  43.55  $${-5.433}^{+0.203}_{-0.374}$$  43.65  $${-5.609}^{+0.253}_{-0.438}$$  43.75  $${-6.212}^{+0.519}_{-0.917}$$  43.85  $${-6.226}^{+0.519}_{-0.917}$$  log L(Lyα)*  log ϕ†  (erg s−1)  ([Δlog L(Lyα)]−1 Mpc−3)  z = 5.7  42.95  $${-3.478}^{+0.038}_{-0.193}$$  43.05  $${-3.735}^{+0.044}_{-0.199}$$  43.15  $${-3.953}^{+0.043}_{-0.198}$$  43.25  $${-4.163}^{+0.055}_{-0.210}$$  43.35  $${-4.427}^{+0.076}_{-0.231}$$  43.45  $${-4.970}^{+0.147}_{-0.308}$$  43.55  $${-5.170}^{+0.187}_{-0.355}$$  43.65  $${-5.318}^{+0.224}_{-0.401}$$  43.75  $${-5.717}^{+0.365}_{-0.606}$$  z = 6.6  43.15  $${-4.194}^{+0.154}_{-0.317}$$  43.25  $${-4.407}^{+0.101}_{-0.258}$$  43.35  $${-4.748}^{+0.087}_{-0.243}$$  43.45  $${-5.132}^{+0.140}_{-0.300}$$  43.55  $${-5.433}^{+0.203}_{-0.374}$$  43.65  $${-5.609}^{+0.253}_{-0.438}$$  43.75  $${-6.212}^{+0.519}_{-0.917}$$  43.85  $${-6.226}^{+0.519}_{-0.917}$$  *The luminosity bin of our Lyα LFs at z = 5.7 and 6.6. The bin size is Δlog L(Lyα) = 0.1. †The number densities corrected for the detection completeness (see subsection 3.1). View Large We fit a Schechter function (Schechter 1976) to our z = 5.7 and z = 6.6 Lyα LFs by minimum χ2 fitting. The Schechter function is defined by   \begin{eqnarray} \phi (L_{\mathrm{Ly}\alpha }) dL_{\mathrm{Ly}\alpha } = \phi ^{*}_{\mathrm{Ly}\alpha } \left( \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right)^{\alpha } \exp \left( - \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right) d \left( \frac{L_{\mathrm{Ly}\alpha }}{L^{*}_{\mathrm{Ly}\alpha }} \right), \nonumber\\ \end{eqnarray} (3)where $$L^{*}_{\mathrm{Ly}\alpha }$$ is the characteristic Lyα luminosity, $$\phi ^{*}_{\mathrm{Ly}\alpha }$$ is the normalization, and α is the faint-end slope. We consider two cases. In one case, we use our Lyα LF measurements at log L(Lyα) [erg s−1] < 43.5, where AGN contamination is not significant in lower-z LAE studies (Ouchi et al. 2008; Konno et al. 2016).1 In the other case, we include the bright-end LF results at log L(Lyα) [erg s−1] ≥ 43.5. In both of these cases, we also use the faint-end Lyα LFs of Ouchi et al. (2008) for z = 5.7 and Ouchi et al. (2010) for z = 6.6. This is because the faint-end Lyα LFs of these studies cover faint Lyα luminosity ranges that we do not reach. Specifically, we include the z = 5.7 Lyα LF data points of Ouchi et al. (2008) in the range log L(Lyα) [erg s−1] = 42.4–42.9 and the z = 6.6 Lyα LF data points of Ouchi et al. (2010) in the range log L(Lyα) [erg s−1] = 42.4–43.0, both of which do not overlap 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). Table 4. Best-fit Schechter parameters and Lyα luminosity densities. Redshift  L*  ϕ*  α  ρLyαobs†    (1043 erg s−1)  (10−4 Mpc−3)    (1039 erg s−1 Mpc−3)  Classical method for log L(Lyα) [erg s−1] = 42.4–43.5†  5.7  $$1.07^{+0.77}_{-0.38}$$  $$2.46^{+3.48}_{-1.86}$$  $$-2.26^{+0.76}_{-0.44}$$  3.39  6.6  $$0.82^{+0.86}_{-0.30}$$  $$2.83^{+3.52}_{-2.38}$$  $$-1.86^{+0.79}_{-0.67}$$  1.96  Classical method for log L(Lyα) [erg s−1] = 42.4–44.0§  5.7  $$1.64^{+2.16}_{-0.62}$$  $$0.849^{+1.87}_{-0.771}$$  $$-2.56^{+0.53}_{-0.45}$$  3.49  6.6  $$1.66^{+0.30}_{-0.69}$$  $$0.467^{+1.44}_{-0.442}$$  $$-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  Redshift  L*  ϕ*  α  ρLyαobs†    (1043 erg s−1)  (10−4 Mpc−3)    (1039 erg s−1 Mpc−3)  Classical method for log L(Lyα) [erg s−1] = 42.4–43.5†  5.7  $$1.07^{+0.77}_{-0.38}$$  $$2.46^{+3.48}_{-1.86}$$  $$-2.26^{+0.76}_{-0.44}$$  3.39  6.6  $$0.82^{+0.86}_{-0.30}$$  $$2.83^{+3.52}_{-2.38}$$  $$-1.86^{+0.79}_{-0.67}$$  1.96  Classical method for log L(Lyα) [erg s−1] = 42.4–44.0§  5.7  $$1.64^{+2.16}_{-0.62}$$  $$0.849^{+1.87}_{-0.771}$$  $$-2.56^{+0.53}_{-0.45}$$  3.49  6.6  $$1.66^{+0.30}_{-0.69}$$  $$0.467^{+1.44}_{-0.442}$$  $$-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  †The Lyα luminosity densities obtained by integrating the Lyα LF down to log LLyα [erg s−1] = 42.4, which corresponds to $$\sim 0.3 \times L^{*}_{\mathrm{Ly}\alpha }\ (z=3\mbox{--}6)$$. ‡The best-fit parameters are derived with the classical method for the Lyα LF measurements in the luminosity range of log L(Lyα) [erg s−1] = 42.4–43.5. §The best-fit parameters are derived with the classical method for the Lyα LF measurements in the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. View Large 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 the following two systematic uncertainties in Lyα LF estimates by the classical method: (I) A Lyα flux of an LAE at a given narrowband magnitude depends on the redshift of the LAE. (II) The minimum EW0 value that corresponds to a given BB − 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α 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 }$$, α) and the standard deviation (σ) of a Gaussian Lyα EW0 probability distribution. LAEs in the mock catalog are uniformly distributed in a comoving volume over the redshift range that the 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 σ of Gaussian Lyα EW0 probability distributions. The faint-end slope α is fixed at α = −2.6 for z = 5.7 and α = −2.5 for z = 6.6, which are the same as those obtained with the classical method for the Lyα LF measurements in the range log L(Lyα) [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 χ2. The best-fit Schechter parameters are summarized in table 4, and examples of the fitting results are shown in figure 10. Fig. 10. View largeDownload slide Examples of the results of our end-to-end Monte Carlo simulations. Top: Best-fit surface number densities for the D-DEEP2-3 LAEs at z = 5.7 (left) and z = 6.6 (right) are shown with solid lines. The other symbols are as in figures 3 and 4. Bottom: Best-fit BB − NB color distributions of LAEs. In the left and right panels, the solid lines denote the results for our z = 5.7 and z = 6.6 LAEs in the D-DEEP2-3 field, respectively. The other symbols are as in figures 5 and 6. (Color online) Fig. 10. View largeDownload slide Examples of the results of our end-to-end Monte Carlo simulations. Top: Best-fit surface number densities for the D-DEEP2-3 LAEs at z = 5.7 (left) and z = 6.6 (right) are shown with solid lines. The other symbols are as in figures 3 and 4. Bottom: Best-fit BB − NB color distributions of LAEs. In the left and right panels, the solid lines denote the results for our z = 5.7 and z = 6.6 LAEs in the D-DEEP2-3 field, respectively. The other symbols are as in figures 5 and 6. (Color online) We show the best-fit functions from the Monte Carlo simulations for our Lyα 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α LFs derived by the classical method. Similar conclusions are obtained by Shimasaku et al. (2006) and Ouchi et al. (2008), who derived the Lyα LFs at z ∼ 3–6 with Subaru/Suprime-Cam. We confirm that the classical method for the Lyα LF calculations gives a good approximation to the true Lyα LF even in the case of our HSC SSP data. The top panel of figure 8 compares the luminosities from the classical method (Lc) and from the simulations (Ls) at the same number densities as a function of Lc. We find that the difference between these two luminosities is only ≲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 ∼1 σ uncertainties shown with the error bars. Moreover, we also find that the classical-method data points are not always underestimated (figure 8 vs. figure 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α LF measurements in the wide luminosity range. If this is true, the faint-end slopes of Lyα LFs are very steep. The best-fit faint-end slope values are α = −2.5–−2.6 (table 4), which may indicate that the faint-end slopes of Lyα LFs are steeper than those of the UV LFs at similar redshifts (e.g., Bouwens et al. 2015a). Note that our best-fit faint-end slopes are steeper than obtained in previous work on the z = 5.7 Lyα LF (Dressler et al. 2015). It should be noted that if we compare our Lyα LF measurements with the best-fit Schechter function results obtained from the classical method, where we consider only the fainter Lyα luminosity range of log LLyα [erg s−1] ≳ 42.5–43.5, we find that there is a significant bright-end excess of the z = 5.7 and z = 6.6 Lyα LF measurements at log LLyα [erg s−1] ≳ 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 ∼3 σ (2.6 σ for z = 5.7 and 3.2 σ 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α 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 subsection 4.1. 3.5 Comparison with previous studies In this section, we compare our Lyα LFs at z = 5.7 and 6.6 with those obtained by previous studies. As shown in figures 8 and 9, our Lyα LFs are generally consistent with those of the previous results. However, our Lyα LF results do not agree with the high number densities of LAEs recently claimed by Matthee et al. (2015) and Santos, Sobral, and Matthee (2016). The reason for this discrepancy is unclear. This study and most of the previous studies have derived the Lyα 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α LF estimates (subsection 3.4). Matthee et al. (2015) and Santos, Sobral, and Matthee (2016) also appear to have considered these two uncertainties; they have applied filter profile correction for Lyα 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α LFs and a Gaussian for Lyα EW0 probability distributions, and have determined their best-fit functions simultaneously based on χ2 fitting to the observed surface number densities and the BB − NB color distributions (subsection 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) estimated the effects of the two uncertainties separately in their subsections 4.1 and 4.3 (see also Santos et al. 2016), which might cause overcorrections due to the redundancy. Another possibility is the difference in the Lyα EW0 distributions. In our simulations we adopted a Gaussian Lyα 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α EW0 distribution in their calcula-tions 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α LF measurements. 4 Discussion 4.1 Systematic effects in the Lyα LF measurements As shown in subsection 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α luminosity range of log LLyα [erg s−1] ≳ 42.5–44.0, are fitted well to the Lyα LF measurements both at the bright end and fainter magnitude bins. However, the best-fit values of the faint-end slope α are very steep compared to the shallower slopes of the UV LFs at similar redshifts (e.g., Bouwens et al. 2015a). Although our results may imply that the wide luminosity range of our Lyα LFs allow us to reveal the true shapes of the Lyα 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 ∼ 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 that are blended at ground-based resolution appear as very bright LAEs (e.g., Bowler et al. 2017a). First, we discuss the possibility of AGNs. Although the number densities of AGNs rapidly decrease from z ∼ 3 toward higher redshift (e.g., Haardt & Madau 2012), some previous studies suggest the existence of (faint) AGNs at z ∼ 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α 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 LLyα [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 (subsection 2.1), we obtain their number densities of 8.6 × 10−7 Mpc−3 and 6.8 × 10−7 Mpc−3 at z = 5.7 and 6.6, respectively. Since the UV magnitudes of the bright-end LAEs are MUV ≳ −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 ∼ 6, which indicates that bright-end LAEs with log LLyα [erg s−1] ≳ 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α emission line and strong highly-ionized metal lines, e.g., N v and C iv (Shibuya et al. 2018b). 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α photons can easily escape into the IGM in the case that the galaxy is surrounded by an ionized bubble large enough to allow the Lyα 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, Zaldarriaga, and Hernquist (2006)— see also Furlanetto and 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 subsection 4.3, we estimate the neutral hydrogen fraction at z = 6.6 to be $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ from the evolution of the Lyα LFs at z = 5.7–6.6. Based on the $$x_\mathrm{H\,{\small I}}$$ value and the top panel of figure 1 of Furlanetto, Zaldarriaga, and Hernquist (2006), we obtain the typical size of ionized bubbles at z = 6.6 of ∼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 ∼15 comoving Mpc. To estimate the sizes of ionized bubbles around bright-end LAEs, we use the following formula for the Strömgren radius RS of an ionized bubble around a source at z = 6.6 by Haiman (2002): RS = 0.8 × (SFR/10 M⊙ yr−1)1/3(t*/100 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–120 M⊙ 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., MUV ≳ −21 mag), we calculate the size of the ionized bubbles of RS ≲ 7 comoving Mpc.2 This size is smaller than that estimated from the analytic model of Furlanetto, Zaldarriaga, and Hernquist (2006) (∼15 comoving Mpc). This result implies that if the bright end of the Lyα 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 that 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. 2018). In addition, we check the positions of massive (Mstar > 1010.3 M⊙) red galaxies with photometric redshift of zphoto = 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 μ ∼ 1.2–1.7. Thus, we conclude that the impact of gravitational lensing on the shapes of the Lyα LFs is small. Finally, we discuss the possibility of blended merging galaxies. Recently, Bowler et al. (2017a) found that multi-component systems account for more than 40% of their bright z ∼ 7 galaxies based on 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. 2018a). 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α 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α LF at z = 5.7–6.6 We investigate the evolution of the Lyα LF at z = 5.7–6.6. In figure 11, we show our Lyα 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 LLyα [erg s−1] ≳ 42.4–44.0 derived with the classical method, which are good approximations to the true LFs (subsection 3.4). We also present the Lyα LF at z = 7.3 derived by Konno et al. (2014) in this figure; they conducted the ultradeep z = 7.3 LAE survey with Subaru/Suprime-Cam. Ouchi et al. (2008, 2010) derived the Lyα LFs at z = 5.7 and 6.6 based on their ∼1 deg2 narrowband imaging data taken with Subaru/Suprime-Cam, and found a decrease of the Lyα LF from z = 5.7 to 6.6. The same results have been obtained by other previous studies (e.g., Kashikawa et al. 2006, 2011; Hu et al. 2010; Santos et al. 2016). We find such evolution from our Lyα 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α LFs shown with the blue and red ovals, respectively. We also show the error contours for the Lyα LF at z = 7.3 of Konno et al. (2014). From this figure, the Schechter parameters of the z = 6.6 Lyα LF are different from those of the z = 5.7 Lyα LF, and the Lyα LF decreases from z = 5.7 to 6.6 at the >90% confidence level. Note that the evolution of the Lyα LFs that we derive is similar to that reported by Santos, Sobral, and Matthee (2016), although our best-fit $$L^{*}_{\mathrm{Ly}\alpha }$$ values are smaller than theirs. The decreasing trend of the Lyα LFs with increasing redshift obtained in this study is also consistent with those of Ouchi et al. (2010), who investigated the evolution of LFs in the faint Lyα range, log L(Lyα) [erg s−1] ≲ 43, as shown in figure 11. It should be noted that the best-fit Lyα 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 α values. In our Schechter function fitting with the classical method, the slope α 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 α value has been fixed at −1.5. Fig. 11. View largeDownload slide Evolution of the Lyα LFs. The blue and red filled circles are our z = 5.7 and 6.6 Lyα LF measurements, respectively, which are derived from the HSC SSP data. The blue and red filled squares denote the z = 5.7 and 6.6 Lyα LF measurements with the Subaru/Suprime-Cam data given by Ouchi et al. (2008, 2010), respectively. The best-fit Schechter function for the Lyα LF at z = 5.7 (6.6) is shown with the blue (red) solid curve, which are derived for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange circles represent the Lyα LF at z = 7.3 derived by Konno et al. (2014). (Color online) Fig. 11. View largeDownload slide Evolution of the Lyα LFs. The blue and red filled circles are our z = 5.7 and 6.6 Lyα LF measurements, respectively, which are derived from the HSC SSP data. The blue and red filled squares denote the z = 5.7 and 6.6 Lyα LF measurements with the Subaru/Suprime-Cam data given by Ouchi et al. (2008, 2010), respectively. The best-fit Schechter function for the Lyα LF at z = 5.7 (6.6) is shown with the blue (red) solid curve, which are derived for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange circles represent the Lyα LF at z = 7.3 derived by Konno et al. (2014). (Color online) Fig. 12. View largeDownload slide Schechter parameter 1 σ and 2 σ confidence intervals, $$L^{*}_{\mathrm{Ly}\alpha }$$ and $$\phi ^{*}_{\mathrm{Ly}\alpha }$$. The blue (red) contours correspond to z = 5.7 (6.6). The blue and red crosses are the best-fit Schechter parameters for the Lyα LFs at z = 5.7 and 6.6, respectively. These results are obtained with the classical method for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange contours show the results for the z = 7.3 Lyα LF (Konno et al. 2014). (Color online) Fig. 12. View largeDownload slide Schechter parameter 1 σ and 2 σ confidence intervals, $$L^{*}_{\mathrm{Ly}\alpha }$$ and $$\phi ^{*}_{\mathrm{Ly}\alpha }$$. The blue (red) contours correspond to z = 5.7 (6.6). The blue and red crosses are the best-fit Schechter parameters for the Lyα LFs at z = 5.7 and 6.6, respectively. These results are obtained with the classical method for the luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. The orange contours show the results for the z = 7.3 Lyα LF (Konno et al. 2014). (Color online) 4.3 Estimation of $$x_\mathrm{H\,{\small I}}$$ at z = 6.6 We estimate the neutral hydrogen fraction, $$x_\mathrm{H\,{\small I}}$$, at z = 6.6 based on our Lyα 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 TIGMLyα, z = 6.6/TIGMLyα, z = 5.7, where TIGMLyα, z is the Lyα transmission through the IGM at a redshift z. The observed Lyα LD, ρLyα, can be obtained from   \begin{equation} \rho ^{\mathrm{Ly}\alpha } = \kappa \ T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z} \ f^{\mathrm{esc}}_{\mathrm{Ly}\alpha } \ \rho ^{\mathrm{UV}}, \end{equation} (4)where κ is the conversion factor from UV to Lyα fluxes, fescLyα is the Lyα escape fraction through the ISM of a galaxy, and ρUV is the intrinsic UV LD. Based on the equation, we can estimate the Lyα transmission fraction TIGMLyα, z = 6.6/TIGMLyα, z = 5.7 by   \begin{equation} \frac{T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z = 6.6}}{T^{\mathrm{IGM}}_{\mathrm{Ly}\alpha , z = 5.7}} = \frac{\kappa _{z = 5.7}}{\kappa _{z = 6.6}} \frac{f^{\mathrm{esc}}_{\mathrm{Ly}\alpha , z = 5.7}}{f^{\mathrm{esc}}_{\mathrm{Ly}\alpha , z = 6.6}} \frac{\rho ^{\mathrm{Ly}\alpha }_{z = 6.6} / \rho ^{\mathrm{Ly}\alpha }_{z = 5.7}}{\rho ^{\mathrm{UV}}_{z = 6.6} / \rho ^{\mathrm{UV}}_{z = 5.7}}. \end{equation} (5)To calculate $$\rho ^{\mathrm{Ly}\alpha , \mathrm{tot}}_{z = 6.6} / \rho ^{\mathrm{Ly}\alpha , \mathrm{tot}}_{z = 5.7}$$, we use the Lyα LD results in subsection 3.4. We adopt the Lyα LDs derived for the Lyα LF measurements in the luminosity range of log L(Lyα) [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. (2015a), ρUVz = 6.6/ρUVz = 5.7 = 0.74 ± 0.10 is obtained. Under the assumption of κz=5.7/κz = 6.6 = 1 and fescLyα, z = 5.7/fescLyα, z = 6.6 = 1, we obtain TIGMLyα, z = 6.6/TIGMLyα, z = 5.7 = 0.70 ± 0.15 from equation (5). We obtain constraints on $$x_\mathrm{H\,{\small I}}$$ based on comparisons of our results with theoretical models. Santos (2004) calculated the IGM Lyα transmission fraction as a function of $$x_\mathrm{H\,{\small I}}$$ in two cases of galactic outflow: Lyα 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α emission is ∼200 km s−1 for LAEs at z ∼ 2 (e.g., Hashimoto et al. 2013; Shibuya et al. 2014). Based on figure 25 of Santos (2004), our Lyα transmission fraction result is consistent with $$x_\mathrm{H\,{\small I}}$$ ∼0.0–0.2, considering the two cases. Next, we compare our Lyα LF result with the theoretical results of McQuinn et al. (2007), who derived z = 6.6 Lyα LFs for various $$x_\mathrm{H\,{\small I}}$$ values based on their radiative transfer simulations. From figure 4 of McQuinn et al. (2007), we obtain constraints of $$x_\mathrm{H\,{\small I}} \sim 0.3$$–0.5. Finally, we compare our result with a combination of two theoretical models. Dijkstra, Wyithe, and Haiman (2007b) derived expected Lyα 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{H\,{\small I}}$$ has been calculated by Furlanetto, Zaldarriaga, and Hernquist (2006) based on their analytic model. A comparison of our Lyα 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{H\,{\small I}} \sim 0.1$$–0.3. Based on the results described above, we conclude that the neutral hydrogen fraction is estimated to be $$x_\mathrm{H\,{\small I}} = 0.1$$–0.5, i.e., $$x_\mathrm{H\,{\small I}} = 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α transmission fraction estimates are considered. Figure 13 shows our $$x_\mathrm{H\,{\small I}}$$ estimate at z = 6.6 and those taken from the previous studies. The previous results of the z ≳ 7 Lyα LFs imply $$x_\mathrm{H\,{\small I}} = 0.3$$–0.8 at z = 7.3 (Konno et al. 2014) and $$x_\mathrm{H\,{\small I}} < 0.63$$ at z = 7.0 (Ota et al. 2010). The studies of Lyα emitting fractions indicate $$x_\mathrm{H\,{\small I}} \gtrsim 0.5$$ at z ∼ 7 (e.g., Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The Lyα damping wing absorption measurements of QSOs suggest $$x_\mathrm{H\,{\small I}} \gtrsim 0.1$$ at z = 7.1 (Mortlock et al. 2011; Bolton et al. 2011). Fig. 13. View largeDownload slide Evolution of the IGM neutral hydrogen fraction. The top and bottom panels are the same, except for the scales of the ordinate axes (top: linear; bottom: logarithmic). The red filled circles show the $$x_\mathrm{H\,{\small I}}$$ estimates from the Lyα LFs at z = 6.6 (this study) and 7.3 (Konno et al. 2014). The blue filled triangle, square, diamond, and pentagon are the $$x_\mathrm{H\,{\small I}}$$ estimates based on the evolution of the Lyα LF obtained by Malhotra and Rhoads (2004), Kashikawa et al. (2011), Ouchi et al. (2010), and Ota et al. (2010), respectively. The blue open diamond and circle are the constraints on $$x_\mathrm{H\,{\small I}}$$ from the clustering analyses of LAEs (Ouchi et al. 2010) and the Lyα emitting galaxy fraction (Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The previous results from the GRB optical afterglow spectrum analyses are shown with magenta filled triangles (Totani et al. 2006, 2014). The green filled squares and open triangle are the results from the GP test of QSOs (Fan et al. 2006) and the size of QSO near zones (Mortlock et al. 2011; Bolton et al. 2011), respectively. The gray and hatched regions are the 1 σ confidence intervals for the instantaneous reionization redshifts obtained by Planck (Planck Collaboration 2016b) and nine-year WMAP (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The models A, B, and C of Choudhury, Ferrara, and Gallerani (2008) are shown with dotted, dashed, and solid lines, respectively. (Color online) Fig. 13. View largeDownload slide Evolution of the IGM neutral hydrogen fraction. The top and bottom panels are the same, except for the scales of the ordinate axes (top: linear; bottom: logarithmic). The red filled circles show the $$x_\mathrm{H\,{\small I}}$$ estimates from the Lyα LFs at z = 6.6 (this study) and 7.3 (Konno et al. 2014). The blue filled triangle, square, diamond, and pentagon are the $$x_\mathrm{H\,{\small I}}$$ estimates based on the evolution of the Lyα LF obtained by Malhotra and Rhoads (2004), Kashikawa et al. (2011), Ouchi et al. (2010), and Ota et al. (2010), respectively. The blue open diamond and circle are the constraints on $$x_\mathrm{H\,{\small I}}$$ from the clustering analyses of LAEs (Ouchi et al. 2010) and the Lyα emitting galaxy fraction (Pentericci et al. 2011, 2014; Schenker et al. 2012, 2014; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014). The previous results from the GRB optical afterglow spectrum analyses are shown with magenta filled triangles (Totani et al. 2006, 2014). The green filled squares and open triangle are the results from the GP test of QSOs (Fan et al. 2006) and the size of QSO near zones (Mortlock et al. 2011; Bolton et al. 2011), respectively. The gray and hatched regions are the 1 σ confidence intervals for the instantaneous reionization redshifts obtained by Planck (Planck Collaboration 2016b) and nine-year WMAP (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The models A, B, and C of Choudhury, Ferrara, and Gallerani (2008) are shown with dotted, dashed, and solid lines, respectively. (Color online) As already pointed out in our previous work (Konno et al. 2014), the decrease of the Lyα LF from z = 6.6 to 7.3 is larger than that from z = 5.7 to 6.6. This accelerated evolution could be also found in figure 13, although the uncertainties are large. The Lyα LF evolves from z = 6.6 to 7.3 at the >90% confidence level, while the difference of $$x_\mathrm{H\,{\small I}}$$ between z = 6.6 and 7.3 is only within 1 σ. This is because, in our $$x_\mathrm{H\,{\small I}}$$ 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α LFs (see Konno et al. 2014 for details). Here, we investigate whether the $$x_\mathrm{H\,{\small I}}$$ evolution obtained by this and previous studies can explain the Thomson scattering optical depth, τel, value obtained from the latest Planck 2016 data. Because one needs to know τel from a given $$x_\mathrm{H\,{\small I}}$$ evolution, we use the semi-analytic models of Choudhury, Ferrara, and Gallerani (2008). They derived $$x_\mathrm{H\,{\small I}}$$ and τel evolutions by considering three models which differ in the minimum halo masses for reionization sources to cover typical scenarios for cosmic reionization history. These three models are referred to as models A, B, and C, corresponding to the minimum halo masses of ∼109, ∼108, and ∼5 × 105 M⊙, respectively, at z = 6. We present the $$x_\mathrm{H\,{\small I}}$$ evolutions of the three models in figure 13, and their τel evolutions in figure 14. The gray (hatched) region in figure 14 shows the 1 σ range of τel obtained by Planck (WMAP). The latest results from the Planck observations indicate that the Thomson scattering optical depth is τel = 0.058 ± 0.009 (Planck Collaboration 2016b), which is significantly lower than that obtained from the WMAP data. In figure 13, models A and B are consistent with our $$x_\mathrm{H\,{\small I}}$$ 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. Model C can barely explain our $$x_\mathrm{H\,{\small I}}$$ value at z = 7.3, but is placed above the τel of Planck beyond the 1 σ error (figure 14). Thus, these results show that cosmic reionization histories such as models A and B can explain both the $$x_\mathrm{H\,{\small I}}$$ estimates and the Planck 2016 τel value simultaneously. Similar conclusions are reached by Robertson et al. (2015) and Bouwens et al. (2015b), who discussed the UV LF evolution of reionization sources independent of our Lyα LF study. Fig. 14. View largeDownload slide Thomson scattering optical depth, τel, as a function of redshift. The gray and hatched regions correspond to the 1 σ confidence intervals of the τel measurements obtained by the Planck 2016 data (Planck Collaboration 2016b) and the nine-year WMAP data (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The dotted, dashed, and solid curves are models A, B, and C of Choudhury, Ferrara, and Gallerani (2008), respectively. Fig. 14. View largeDownload slide Thomson scattering optical depth, τel, as a function of redshift. The gray and hatched regions correspond to the 1 σ confidence intervals of the τel measurements obtained by the Planck 2016 data (Planck Collaboration 2016b) and the nine-year WMAP data (Hinshaw et al. 2013; Bennett et al. 2013), respectively. The dotted, dashed, and solid curves are models A, B, and C of Choudhury, Ferrara, and Gallerani (2008), respectively. 5 Summary We have derived the Lyα 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 ∼13.8 deg2 and ∼21.2 deg2, respectively. The 5 σ limiting magnitudes of the narrowband images are ∼25.0 mag and ∼25.5 mag in the Deep and UltraDeep layers, respectively. Using these narrowband images, we have identified, in total, ∼2000 LAEs at z = 5.7 and 6.6 with a bright Lyα luminosity range of log L(Lyα) [erg s−1] ≃ 42.9–43.8. Our HSC LAE sample is ∼2–6 times larger than those of previous studies of z ∼ 6–7 LAEs. Based on the LAE samples, we have derived the Lyα 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α 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α LF, if we consider the Lyα luminosity range of log L(Lyα) [erg s−1] = 42.4–44.0. Our Lyα 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. We have confirmed the decrease of the Lyα LF from z = 5.7 to 6.6. This evolution is caused by the Lyα damping wing absorption of neutral hydrogen in the IGM. Based on the decrease of the Lyα LF at z = 5.7–6.6, we have estimated the IGM neutral hydrogen fraction of $$x_\mathrm{H\,{\small I}} = 0.3 \pm 0.2$$ at z = 6.6. The $$x_\mathrm{H\,{\small I}}$$ evolution obtained from this and previous studies can explain the Thomson scattering optical depth measurement of the latest Planck 2016. Acknowledgements We thank Mamoru Doi, Kentaro Motohara, Toshitaka Kajino, and Masafumi Ishigaki for useful discussion and comments. We acknowledge 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 the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), the 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, 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 the Japan Society for the Promotion of Science. N.K. acknowledges support from JSPS grant 15H03645. Footnotes 1 As mentioned in subsection 4.1, Shibuya et al. (2018b) found no clear signature of AGNs for several bright LAEs with log L(Lyα) [erg s−1] > 43.5. Their bright LAEs show narrow Lyα line widths of <400 km s−1 and no clear detection of UV lines such as N v and C iv. However, their investigation was based on rest-frame UV spectroscopic observations and they could not 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α LF fit results for the two cases where we include and exclude the bright-end bins of log L(Lyα) [erg s−1] > 43.5 for a conservative discussion. 2 The SFRs can be estimated from UV luminosities with the following equation (Madau et al. 1998): SFR(M⊙ yr−1) = LUV (erg s−1 Hz−1)/(8 × 1027). From this equation, the SFR corresponding to MUV = −21 is 13.6 M⊙ yr−1. We estimate the ionized bubble size under the assumption that that these bright LAEs have a constant SFR of 13.6 M⊙ yr−1, and emit ionizing photons during their age of 100 Myr. References Adams J. J. et al.   2011, ApJS , 192, 5 CrossRef Search ADS   Aihara H. et al.   2018a, PASJ , 70, S4 Aihara H. et al.   2018b, PASJ , 70, S8 Ajiki M. et al.   2002, ApJ , 576, L25 CrossRef Search ADS   Axelrod T., Kantor J., Lupton R. H., Pierfederici F. 2010, in Proc. SPIE, 7740, Software and Cyberinfrastructure for Astronomy , ed. Radziwill N. M., Bridger A. ( Bellingham, WA: SPIE), 774015 Bagley M. B. et al.   2017, ApJ , 837, 11 CrossRef Search ADS   Barone-Nugent R. L., Wyithe J. S. B., Trenti M., Treu T., Oesch P., Bouwens R., Illingworth G. D., Schmidt K. B. 2015, MNRAS , 450, 1224 CrossRef Search ADS   Bennett C. L. et al.   2013, ApJS , 208, 20 CrossRef Search ADS   Bolton J. S., Haehnelt M. G., Warren S. J., Hewett P. C., Mortlock D. J., Venemans B. P., McMahon R. G., Simpson C. 2011, MNRAS , 416, L70 CrossRef Search ADS   Bosch J. et al.   2018, PASJ , 70, S5 Bouwens R. J. et al.   2015a, ApJ , 803, 34 CrossRef Search ADS   Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B., Smit R., Wilkins S. 2015b, ApJ , 811, 140 CrossRef Search ADS   Bowler R. A. A., Dunlop J. S., McLure R. J., McLeod D. J. 2017a, MNRAS , 466, 3612 CrossRef Search ADS   Bowler R. A. A., McLure R. J., Dunlop J. S., McLeod D. J., Stanway E. R., Eldridge J. J., Jarvis M. J. 2017b, MNRAS , 469, 448 CrossRef Search ADS   Cai Z. et al.   2017, ApJ , 837, 71 CrossRef Search ADS   Cantalupo S., Arrigoni-Battaia F., Prochaska J. X., Hennawi J. F., Madau P. 2014, Nature , 506, 63 CrossRef Search ADS PubMed  Caruana J., Bunker A. J., Wilkins S. M., Stanway E. R., Lacy M., Jarvis M. J., Lorenzoni S., Hickey S. 2012, MNRAS , 427, 3055 CrossRef Search ADS   Caruana J., Bunker A. J., Wilkins S. M., Stanway E. R., Lorenzoni S., Jarvis M. J., Ebert H. 2014, MNRAS , 443, 2831 CrossRef Search ADS   Cassata P. et al.   2015, A&A , 573, A24 CrossRef Search ADS   Castellano M. et al.   2016, ApJ , 818, L3 CrossRef Search ADS   Choudhury T. R., Ferrara A., Gallerani S. 2008, MNRAS , 385, L58 CrossRef Search ADS   Coupon J., Czakon N., Bosch J., Komiyama Y., Medezinski E., Miyazaki S., Oguri M. 2018, PASJ , 70, S7 Cowie L. L., Hu E. M. 1998, AJ , 115, 1319 CrossRef Search ADS   Dayal P., Maselli A., Ferrara A. 2011, MNRAS , 410, 830 CrossRef Search ADS   Deharveng J.-M. et al.   2008, ApJ , 680, 1072 CrossRef Search ADS   Dijkstra M., Lidz A., Wyithe J. S. B. 2007a, MNRAS , 377, 1175 CrossRef Search ADS   Dijkstra M., Wyithe J. S. B., Haiman Z. 2007b, MNRAS , 379, 253 CrossRef Search ADS   Dressler A. et al.   2011, PASP , 123, 288 CrossRef Search ADS   Dressler A., Henry A., Martin C. L., Sawicki M., McCarthy P., Villaneuva E. 2015, ApJ , 806, 19 CrossRef Search ADS   Fan X. et al.   2006, AJ , 132, 117 CrossRef Search ADS   Finkelstein S. L. et al.   2013, Nature , 502, 524 CrossRef Search ADS PubMed  Furlanetto S. R., Oh S. P. 2005, MNRAS , 363, 1031 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L. 2006, MNRAS , 365, 1012 CrossRef Search ADS   Furusawa H. et al.   2018, PASJ , 70, S3 Gehrels N. 1986, ApJ , 303, 336 CrossRef Search ADS   Giallongo E. et al.   2015, A&A , 578, A83 CrossRef Search ADS   Gronwall C. et al.   2007, ApJ , 667, 79 CrossRef Search ADS   Guaita L. et al.   2010, ApJ , 714, 255 CrossRef Search ADS   Haardt F., Madau P. 2012, ApJ , 746, 125 CrossRef Search ADS   Haiman Z. 2002, ApJ , 576, L1 CrossRef Search ADS   Haiman Z., Spaans M. 1999, ApJ , 518, 138 CrossRef Search ADS   Hall P. B. et al.   2004, AJ , 127, 3146 CrossRef Search ADS   Harikane Y. et al.   2018, PASJ , 70, S11 Hashimoto T., Ouchi M., Shimasaku K., Ono Y., Nakajima K., Rauch M., Lee J., Okamura S. 2013, ApJ , 765, 70 CrossRef Search ADS   Hayashino T. et al.   2004, AJ , 128, 2073 CrossRef Search ADS   Hayes M., Schaerer D., Östlin G., Mas-Hesse J. M., Atek H., Kunth D. 2011, ApJ , 730, 8 CrossRef Search ADS   Hinshaw G. et al.   2013, ApJS , 208, 19 CrossRef Search ADS   Hu E. M., Cowie L. L., Barger A. J., Capak P., Kakazu Y., Trouille L. 2010, ApJ , 725, 394 CrossRef Search ADS   Hu E. M., Cowie L. L., McMahon R. G. 1998, ApJ , 502, L99 CrossRef Search ADS   Hu E. M., Cowie L. L., Songaila A., Barger A. J., Rosenwasser B., Wold I. G. B. 2016, ApJ , 825, L7 CrossRef Search ADS   Huang S. et al.   2018, PASJ , 70, S6 Inoue A. K. et al.   2016, Science , 352, 1559 CrossRef Search ADS PubMed  Ivezic Z. et al.   2008, arXiv:0805.2366 Iye M. et al.   2006, Nature , 443, 186 CrossRef Search ADS PubMed  Jiang L. et al.   2016, ApJ , 833, 222 CrossRef Search ADS   Jurić M. et al.   2015, arXiv:1512.07914 Kashikawa N. et al.   2002, PASJ , 54, 819 CrossRef Search ADS   Kashikawa N. et al.   2006, ApJ , 648, 7 CrossRef Search ADS   Kashikawa N. et al.   2011, ApJ , 734, 119 CrossRef Search ADS   Kashikawa N. et al.   2015, ApJ , 798, 28 CrossRef Search ADS   Kawanomoto S. et al.   2017, PASJ , submitted Kobayashi M. A. R., Totani T., Nagashima M. 2007, ApJ , 670, 919 CrossRef Search ADS   Komiyama Y. et al.   2018, PASJ , 70, S2 Konno A. et al.   2014, ApJ , 797, 16 CrossRef Search ADS   Konno A. et al.   2016, ApJ , 823, 20 CrossRef Search ADS   Kusakabe H., Shimasaku K., Nakajima K., Ouchi M. 2015, ApJ , 800, L29 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S. 2007, MNRAS , 381, 75 CrossRef Search ADS   Madau P. 1995, ApJ , 441, 18 CrossRef Search ADS   Madau P., Pozzetti L., Dickinson M. 1998, ApJ , 498, 106 CrossRef Search ADS   Magnier E. A. et al.   2013, ApJS , 205, 20 CrossRef Search ADS   Malhotra S., Rhoads J. E. 2002, ApJ , 565, L71 CrossRef Search ADS   Malhotra S., Rhoads J. E. 2004, ApJ , 617, L5 CrossRef Search ADS   Mao J., Lapi A., Granato G. L., de Zotti G., Danese L. 2007, ApJ , 667, 655 CrossRef Search ADS   Martínez-Sansigre A., Rawlings S., Lacy M., Fadda D., Jarvis M. J., Marleau F. R., Simpson C., Willott C. J. 2006, MNRAS , 370, 1479 CrossRef Search ADS   Mason C. A. et al.   2015, ApJ , 805, 79 CrossRef Search ADS   Matsuda Y. et al.   2004, AJ , 128, 569 CrossRef Search ADS   Matthee J., Sobral D., Santos S., Röttgering H., Darvish B., Mobasher B. 2015, MNRAS , 451, 400 CrossRef Search ADS   Mesinger A., Furlanetto S. R. 2008, MNRAS , 386, 1990 CrossRef Search ADS   Miyazaki S. et al.   2012, in Proc. SPIE, 8446, Ground-based and Airborne Instrumentation for Astronomy IV , ed. McLean I. S. et al.   ( Bellingham, WA: SPIE), 84460Z Miyazaki S. et al.   2018, PASJ , 70, S1 Mortlock D. J. et al.   2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murayama T. et al.   2007, ApJS , 172, 523 CrossRef Search ADS   Nakajima K., Ouchi M. 2014, MNRAS , 442, 900 CrossRef Search ADS   Oesch P. A. et al.   2015, ApJ , 804, L30 CrossRef Search ADS   Oguri M. 2014, MNRAS , 444, 147 CrossRef Search ADS   Oguri M. et al.   2018, PASJ , 70, S20 Oke J. B. 1974, ApJS , 27, 21 CrossRef Search ADS   Ono Y. et al.   2010a, MNRAS , 402, 1580 CrossRef Search ADS   Ono Y. et al.   2012, ApJ , 744, 83 CrossRef Search ADS   Ono Y. et al.   2018, PASJ , 70, S10 Ono Y., Ouchi M., Shimasaku K., Dunlop J., Farrah D., McLure R., Okamura S. 2010b, ApJ , 724, 1524 CrossRef Search ADS   Ota K. et al.   2010, ApJ , 722, 803 CrossRef Search ADS   Ota K. et al.   2017, ApJ , 844, 85 CrossRef Search ADS   Ouchi M. et al.   2003, ApJ , 582, 60 CrossRef Search ADS   Ouchi M. et al.   2008, ApJS , 176, 301 CrossRef Search ADS   Ouchi M. et al.   2009, ApJ , 696, 1164 CrossRef Search ADS   Ouchi M. et al.   2010, ApJ , 723, 869 CrossRef Search ADS   Ouchi M. et al.   2013, ApJ , 778, 102 CrossRef Search ADS   Ouchi M. et al.   2018, PASJ , 70, S13 Pacucci F., Pallottini A., Ferrara A., Gallerani S. 2017, MNRAS , 468, L77 CrossRef Search ADS   Parsa S., Dunlop J. S., McLure R. J. 2017, arXiv:1704.07750 Pentericci L. et al.   2011, ApJ , 743, 132 CrossRef Search ADS   Pentericci L. et al.   2014, ApJ , 793, 113 CrossRef Search ADS   Planck Collaboration 2016a, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration 2016b, A&A , 596, A107 CrossRef Search ADS   Reddy N. A., Steidel C. C. 2009, ApJ , 692, 778 CrossRef Search ADS   Rhoads J. E., Malhotra S., Dey A., Stern D., Spinrad H., Jannuzi B. T. 2000, ApJ , 545, L85 CrossRef Search ADS   Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S. 2015, ApJ , 802, L19 CrossRef Search ADS   Santos M. R. 2004, MNRAS , 349, 1137 CrossRef Search ADS   Santos S., Sobral D., Matthee J. 2016, MNRAS , 463, 1678 CrossRef Search ADS   Schaerer D., Hayes M., Verhamme A., Teyssier R. 2011, A&A , 531, A12 CrossRef Search ADS   Schechter P. 1976, ApJ , 203, 297 CrossRef Search ADS   Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P. 2014, ApJ , 795, 20 CrossRef Search ADS   Schenker M. A., Stark D. P., Ellis R. S., Robertson B. E., Dunlop J. S., McLure R. J., Kneib J.-P., Richard J. 2012, ApJ , 744, 179 CrossRef Search ADS   Schiminovich D. et al.   2005, ApJ , 619, L47 CrossRef Search ADS   Schlafly E. F. et al.   2012, ApJ , 756, 158 CrossRef Search ADS   Schlegel D. J., Finkbeiner D. P., Davis M. 1998, ApJ , 500, 525 CrossRef Search ADS   Shibuya T. et al.   2014, ApJ , 788, 74 CrossRef Search ADS   Shibuya T. et al.   2018a, PASJ , 70, S14 Shibuya T. et al.   2018b, PASJ , 70, S15 Shibuya T., Kashikawa N., Ota K., Iye M., Ouchi M., Furusawa H., Shimasaku K., Hattori T. 2012, ApJ , 752, 114 CrossRef Search ADS   Shibuya T., Ouchi M., Harikane Y. 2015, ApJS , 219, 15 CrossRef Search ADS   Shimasaku K. et al.   2006, PASJ , 58, 313 CrossRef Search ADS   Sobral D., Matthee J., Darvish B., Schaerer D., Mobasher B., Röttgering H. J. A., Santos S., Hemmati S. 2015, ApJ , 808, 139 CrossRef Search ADS   Song M., Finkelstein S. L., Livermore R. C., Capak P. L., Dickinson M., Fontana A. 2016, ApJ , 826, 113 CrossRef Search ADS   Stark D. P. et al.   2017, MNRAS , 464, 469 CrossRef Search ADS   Steidel C. C., Adelberger K. L., Shapley A. E., Pettini M., Dickinson M., Giavalisco M. 2000, ApJ , 532, 170 CrossRef Search ADS   Takahashi R., Oguri M., Sato M., Hamana T. 2011, ApJ , 742, 15 CrossRef Search ADS   Taniguchi Y. et al.   2005, PASJ , 57, 165 CrossRef Search ADS   Tonry J. L. et al.   2012, ApJ , 750, 99 CrossRef Search ADS   Toshikawa J. et al.   2018, PASJ , 70, S12 Totani T. et al.   2014, PASJ , 66, 63 CrossRef Search ADS   Totani T., Kawai N., Kosugi G., Aoki K., Yamada T., Iye M., Ohta K., Hattori T. 2006, PASJ , 58, 485 CrossRef Search ADS   Treu T., Trenti M., Stiavelli M., Auger M. W., Bradley L. D. 2012, ApJ , 747, 27 CrossRef Search ADS   Verhamme A., Schaerer D., Maselli A. 2006, A&A , 460, 397 CrossRef Search ADS   Willott C. J. et al.   2010, AJ , 139, 906 CrossRef Search ADS   Wyithe J. S. B., Yan H., Windhorst R. A., Mao S. 2011, Nature , 469, 181 CrossRef Search ADS PubMed  Yamada T., Nakamura Y., Matsuda Y., Hayashino T., Yamauchi R., Morimoto N., Kousai K., Umemura M. 2012, AJ , 143, 79 CrossRef Search ADS   Zabl J., Nørgaard-Nielsen H. U., Fynbo J. P. U., Laursen P., Ouchi M., Kjærgaard P. 2015, MNRAS , 451, 2050 CrossRef Search ADS   Zheng Z.-Y. et al.   2017, ApJ , 842, L22 CrossRef Search ADS   Zitrin A. et al.   2015, ApJ , 810, L12 CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of the Astronomical Society of Japan. All rights reserved. For Permissions, please email: journals.permissions@oup.com

Journal

Publications of the Astronomical Society of JapanOxford University Press

Published: Jan 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off