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