SILVERRUSH. IV. Lyα Luminosity Functions at z = 5.7 and 6.6 Studied with sim 1,300 LAEs on the 14-21 deg² Sky

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