# Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH). I. Program strategy and clustering properties of ∼2000 Lyα emitters at z = 6–7 over the 0.3–0.5 Gpc2 survey area

Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru... Abstract We present the SILVERRUSH program strategy and clustering properties investigated with ∼2000 Lyα emitters (LAEs) at z = 5.7 and 6.6 found in the early data of the Hyper Suprime-Cam (HSC) Subaru Strategic Program survey exploiting the carefully designed narrow-band filters. We derive angular correlation functions with the unprecedentedly large samples of LAEs at z = 6–7 over the large total area of 14–21 deg2 corresponding to 0.3–0.5 comoving Gpc2. We obtain the average large-scale bias values of bavg = 4.1 ± 0.2 (4.5 ± 0.6) at z = 5.7 (z = 6.6) for ≳ L* LAEs, indicating a weak evolution of LAE clustering from z = 5.7 to 6.6. We compare the LAE clustering results with two independent theoretical models that suggest an increase of an LAE clustering signal by the patchy ionized bubbles at the epoch of reionization (EoR), and estimate the neutral hydrogen fraction to be $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ at z = 6.6. Based on the halo occupation distribution models, we find that the ≳ L* LAEs are hosted by dark-matter halos with an average mass of $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) =11.1^{+0.2}_{-0.4}$$ ($$10.8^{+0.3}_{-0.5}$$) at z = 5.7 (6.6) with a Lyα duty cycle of 1% or less, where the results of z = 6.6 LAEs may be slightly biased, due to the increase of the clustering signal at the EoR. Our clustering analysis reveals the low-mass nature of ≳ L* LAEs at z = 6–7, and that these LAEs probably evolve into massive super-L* galaxies in the present-day universe. 1 Introduction Lyα emitters (LAEs) are star-forming galaxies (and active galactic nuclei, AGNs) with a strong Lyα emission line. In the 1960s, it was theoretically predicted that such galaxies are candidates of very young galaxies residing at z ∼ 10–30 (Partridge & Peebles 1967). About 30 years after these predictions, deep observations identified a few LAEs around AGN regions at z = 4.6 (Hu & McMahon 1996) and z = 2.4 (Pascarelle et al. 1996). Subsequently, LAEs in blank fields have routinely been found by deep and wide-field narrow-band observations conducted under the Hawaii Survey (Cowie & Hu 1998), the Large Area Lyman Alpha Survey (Rhoads et al. 2000), and Subaru Surveys (e.g., Ouchi et al. 2003). These blank-field LAE surveys reveal various statistical properties of LAEs up to z ∼ 7 including the evolution of Lyα luminosity functions (LFs; Malhotra & Rhoads 2004; Kashikawa et al. 2006, 2011; Deharveng et al. 2008; Ouchi et al. 2008, 2010; Cowie et al. 2010; Hu et al. 2010; Konno et al. 2014; Matthee et al. 2015; Santos et al. 2016; Ota et al. 2017; Zheng et al. 2017). Moreover, multi-wavelength follow-up imaging and spectroscopy of LAEs reveal that typical (∼L*) LAEs (LLyα = 1042–1043 erg s−1) at z ∼ 2–3 have a stellar mass of 108–109 M⊙ (Gawiser et al. 2007; Hagen et al. 2014), a star-formation rate (SFR) of ∼10 M⊙ (Nakajima et al. 2012), and a gas-phase metallicity of ∼0.1 Z⊙ (Finkelstein et al. 2011; Guaita et al. 2013; Kojima et al. 2017). LAEs have strong high-ionization lines such as [O iii]5007, indicative of the high ionization state of the inter-stellar medium (ISM) given by intense ionizing radiation from young massive stars (Nakajima & Ouchi 2014). Outside the highly-ionized ISM, LAEs show diffuse extended Lyα halos in the circum-galactic medium (CGM) up to a few 10 kpc (Hayashino et al. 2004; Steidel et al. 2011). The bright Lyα halos extending over a few hundred kpc are known as Lyα blobs and their total Lyα luminosities are ∼1044 erg s−1 (Steidel et al. 2000; Matsuda et al. 2004; cf. Saito et al. 2006), more than an order of magnitude brighter than the diffuse Lyα halos (Momose et al. 2014, 2016). Deep spectroscopic observations identified large spatial overdensities of LAEs up to z ∼ 6 (Shimasaku et al. 2003; Ouchi et al. 2005a) that are possibly progenitors of massive galaxy clusters found today. Recent deep spectroscopy also detected strong Lyα emission of ∼10 LAEs at z = 7.0–8.7, pushing the redshift frontier of LAEs to-date (e.g., Pentericci et al. 2011; Ono et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Schenker et al. 2014; Oesch et al. 2015; Zitrin et al. 2015). There are three major scientific goals in recent LAE studies. The first goal is characterizing the nature of low-mass young galaxies at high-z. Because LAEs are low stellar mass galaxies with strong Lyα indicative of a starburst with young massive stars producing a large amount of ionizing photons (Ouchi et al. 2013; Sobral et al. 2015; Stark et al. 2015), LAEs are used as probes of young galaxies. This motivation is the same as the one behind the original predictions of LAEs given by Partridge and Peebles (1967). The second goal is to understand the density and dynamics of H i clouds in the ISM and the CGM of star-forming galaxies. Due to the resonance nature of Lyα emission, Lyα photons are scattered by H i gas in the ISM and the CGM, and the dynamical properties of Lyα sources are not easily understood with the observed Lyα velocity fields. Instead, both the density and kinematics of H i gas in the ISM and the CGM are encoded in the observed Lyα line velocity and spatial profiles. The combination of the Lyα line profile observations and models provides important information about the H i gas of the ISM and the CGM (Verhamme et al. 2008; Zheng et al. 2011; Shibuya et al. 2014; Hashimoto et al. 2015; Momose et al. 2016). The third goal is to reveal the cosmic reionization history and its properties. Due to the strong damping absorption of Lyα by the IGM H i gas at the epoch of reionization (EoR: z > 6), Lyα emission of LAEs is absorbed in the partially neutral IGM. The Lyα absorption depress the observed Lyα luminosities of LAEs, which causes a decrease in the Lyα LF towards the early stage of the EoR (Malhotra & Rhoads 2004; Kashikawa et al. 2006, 2011; Ouchi et al. 2010; Hu et al. 2010; Santos et al. 2016; Ota et al. 2017). Moreover, the clustering signal of observed LAEs is boosted by the existence of the ionized bubbles (Furlanetto et al. 2006; McQuinn et al. 2007; Ouchi et al. 2010), because the Lyα absorption is selectively weak for LAEs residing in the ionized bubbles of the IGM at the EoR. The Lyα LFs and clustering of LAEs are important quantities in characterizing the cosmic reionization history and ionized bubble topologies at the EoR. There are three new large instruments that are used for LAE observations; the Very Large Telescope/Multi-Unit Spectroscopic Explorer (MUSE: Bacon et al. 2010), the Hobby–Eberly Telescope Dark Energy Experiment/Visible Integral-Field Replicable Unit Spectrograph (HETDEX/VIRUS: Hill et al. 2012), and the Subaru/Hyper Suprime-Cam (HSC: Miyazaki et al. 2012). These three instruments can cover the complementary parameter space of LAEs in redshift, depth, and survey volume. Specifically, Subaru/HSC has the capability to take large-area (>10 deg2) deep images with custom-made narrow-bands (section 2) targeting LAEs at z ∼ 2–7. The combination of Subaru/HSC imaging and deep follow-up spectroscopy allows us to address the key issues of LAEs for accomplishing the three major goals described above. One of the most unique studies realized with Subaru/HSC is clustering of LAEs at the EoR (i.e., z ≳ 6), which requires high-sensitivity observations over a large area of the sky, as we show in this study. LAE clustering signals depend on the hosting dark-matter halo properties of LAEs, and the distribution of the patchy ionized IGM (i.e., ionized bubbles) that allows Lyα photons to escape from the partially neutral IGM at the EoR (Furlanetto et al. 2006; McQuinn et al. 2007; Ouchi et al. 2010). Theoretical models suggest that a Lyα LF measurement, a popular statistical quantity of LAEs, has a degeneracy between the Lyα escape fraction $$f_{\rm esc}^{\rm Ly \alpha }$$ (a flux ratio of observed to intrinsic Lyα emission: Le Delliou et al. 2006; Mao et al. 2007; Kobayashi et al. 2010) and duty cycle $$f_{\rm duty}^{\rm LAE}$$ (an abundance ratio of Lyα emitting galaxies to dark-matter haloes: Nagamine et al. 2010). Moreover, at the EoR, there is another degeneracy between the ionized bubble topology and neutral hydrogen fraction xH i of the IGM (Kakiichi et al. 2016). Theoretical studies show that these degeneracies can be resolved with a combination of Lyα LF and clustering measurements (e.g., Nagamine et al. 2010; Sobacchi & Mesinger 2015; Hutter et al. 2015; Kakiichi et al. 2016). The LAE samples at the post-EoR, e.g., z = 5.7, (the EoR, e.g., z = 6.6) are useful to determine $$f_{\rm esc}^{\rm Ly \alpha }$$ and $$f_{\rm duty}^{\rm LAE}$$. In this paper, we describe our narrow-band filters and the program strategy of the HSC LAE studies named “Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC” (SILVERRUSH: section 2), and present our LAE samples (section 3). We obtain measurements of LAE clustering (section 4) that are calculated with the early samples of HSC LAEs at z = 5.7 and 6.6 (Shibuya et al. 2018a) based on the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP: HSC et al. in preparation) data whose first data release is presented in Aihara et al. (2018). Based on the LAE clustering measurements, we discuss the cosmic reionization history and dark-matter halo properties of LAEs by a comparison with the ΛCDM structure formation models in section 5. Throughout this paper, magnitudes are in the AB system. We adopt a cosmology parameter set of (h, Ωm, ΩΛ, ns, σ8) = (0.7, 0.3, 0.7, 1.0, 0.8), consistent with the latest Planck results (Planck Collaboration 2016). 2 Narrow-band filters and the program strategy for the LAE studies In our LAE studies, we use four custom-made narrow-band filters that are fabricated for the HSC-SSP survey. We design the HSC narrow-band filters targeting the OH sky windows at 816, 921, and 1010 nm; these correspond to NB816 (PI: Y. Taniguchi), NB921 (PI: M. Ouchi), and NB101 (PI: K. Shimasaku) filters, which identify the redshifted Lyα of LAEs at z = 5.7, 6.6, and 7.3, respectively. These three narrow-band filters densely cover LAEs from the heart of the EoR to the post-EoR. Another HSC narrow-band filter is NB387 (PI: K. Shimasaku), the central wavelength of which is 387 nm. The NB387 band is the bluest narrow-band among the existing HSC filters, targeting the lowest-redshift LAEs identified by the HSC observations. The central wavelength of 387 nm is specifically chosen, because this band identifies LAEs at z = 2.2 that have major strong lines of [O ii]3727 and Hα6563 which fall in the OH sky windows of 1.19 and 2.10 μm, respectively. The other major strong lines, Hβ4861 and [O iii]5007, of z = 2.2 LAEs are placed at the H band, where the effects of the OH sky emission are unavoidable, due to no clear OH windows. We carefully determine the full width at half maximum (FWHM) value for each narrow-band filter transmission curve. First, we determine the FWHM of NB921 that properly covers a clear OH window at the reddest band, allowing the typical transmission curve errors of the filter fabrication.1 We aim to accomplish the same detection limit of the rest-frame Lyα equivalent width EW0 in our samples of LAEs at z = 2.2, 5.7, and 6.6. An observed-frame Lyα equivalent width, EWobs, has the relation of EWobs = (1 + z)EW0. Because a narrow-band flux excess is proportional (inversely proportional) to $$\mathit {EW}_{\rm obs}$$ (a narrow-band FWHM value), we scale the FWHM value of the NB921 band (FWHMNB921) by (1 + z) to design the FWHM values of NB387 (FWHMNB387) and NB816 (FWHMNB816). More specifically, the determined FWHM values keep the relations of FWHMNB387 = FWHMNB921(1 + 2.18)/(1 + 6.58) and FWHMNB816 = FWHMNB921(1 + 5.71)/(1 + 6.58), where the numbers of 2.18, 5.71, and 6.58 correspond to the redshifts of the Lyα lines (121.567 nm) that fall in the centers of the NB387, NB816, and NB921 bandpasses, respectively. Here, the FWHM of the NB101 band does not follow this scaling relation. Instead, we choose an FWHM that is the narrowest limit that can be reasonably accomplished by the present technology of the HSC narrow-band filter fabrication. This is because Lyα lines of the NB101 LAEs fall in a wavelength range where the HSC imaging cannot reach a very deep magnitude limit, due to the relatively low CCD quantum efficiency and the unclean OH window. In this way, we obtain the designed specifications of (central wavelength, FWHM wavelength) = (387.0 nm, 5.5 nm) for NB387, (816.0 nm, 11.6 nm) for NB816, (921.0 nm, 13.1 nm) for NB921, and (1009.5 nm, 9.0 nm) for NB101. The HSC narrow-band filters are made of a B270 or quartz glass with multi-layer interferometric coatings that make the sharp cut-on and cut-off bandpasses for the narrow-band filters. The diameter of the HSC narrow-band filter is large; 600 mm. Although uniform coatings are key to accomplishing the spatially homogeneous transmission of the narrow-band filters, small non-uniformities are included in the fabrication processes which give a filter characterization that is slightly different from the design. Nevertheless, we have tried to make the coating distribution as uniform as possible to achieve spatially homogeneous transmission curves for the narrow-band filters. The picture of NB921 is shown in figure 1. Fig. 1. View largeDownload slide HSC NB921 filter in a filter holder. (Color online) Fig. 1. View largeDownload slide HSC NB921 filter in a filter holder. (Color online) Figures 2 and 3 present the transmission curves of the HSC NB816 and NB921 filters that are used in the HSC first-data release (Aihara et al. 2018). We measure the transmission curves at 21 points over the 600-mm-diameter circular filter. The detectors of HSC cover the circular area from the center r = 0 to the radius of r = 280 mm. On the filter, we define the axis of the 0° angle for the arbitrary direction from the filter center to the edge. At an angle of 0°, we measure transmission curves at five positions; r = 50, 100, 150, 200, and 270 mm. We change the angle to 90°, 180°, and 270° in a counter-clockwise direction, and thereby obtain transmission curves at 20 positions (= 5 × 4). Including the measurement at the center (r = 0 mm), we have transmission curves at a total of 21 positions. Note that the measurement position of NB921 is slightly different from that of NB816. The r = 270 mm position is replaced with r = 265 mm for the NB921 filter. We find that the area-weighted average central wavelength and FWHM of NB816 (NB921) are 817.7 nm and 11.3 nm (921.5 nm and 13.5 nm), respectively, and that the central-position central wavelength and FWHM of NB816 (NB921) are 816.8 nm and 11.0 nm (920.5 nm and 13.3 nm), respectively. At any positions within r ≤ 265 mm, the peak transmissions of NB816 and NB921 are high; >90%. The deviations of the central and FWHM wavelengths are typically within ≃ 0.3% and ≃ 10%, respectively. Fig. 2. View largeDownload slide Filter transmission curves of NB816. The gray lines represent the transmission curves at the 21 positions (see text). The black line is the same as the gray lines, but for the transmission curve at the central position. The red line represents the area-averaged mean transmission curve that is shown in the Subaru website ⟨http://www.naoj.org/Observing/Instruments/HSC/sensitivity.html⟩. (Color online) Fig. 2. View largeDownload slide Filter transmission curves of NB816. The gray lines represent the transmission curves at the 21 positions (see text). The black line is the same as the gray lines, but for the transmission curve at the central position. The red line represents the area-averaged mean transmission curve that is shown in the Subaru website ⟨http://www.naoj.org/Observing/Instruments/HSC/sensitivity.html⟩. (Color online) Fig. 3. View largeDownload slide Same as figure 2, but for NB921. (Color online) Fig. 3. View largeDownload slide Same as figure 2, but for NB921. (Color online) Using the HSC narrow-band filters, we study LAEs at z = 2.2, 5.7, and 6.6 (7.3) over the large areas of 26 (3.5) deg2, which are about an order of magnitude larger than those of the previous studies for a given redshift. Exploiting samples of z = 2.2 and 5.7 LAEs, we study the evolution of LAEs at z = 2.2–5.7, from the lowest redshift accessible by ground-based observations (z ∼ 2) to the high-z edge of the post-EoR (z ∼ 6). The LAEs at the three redshifts of z = 5.7, 6.6, and 7.3 allow us to study the evolution of galaxy formation as well as cosmic reionization with the Lyα damping wing absorption of the neutral IGM. The LAEs at the post-EoR of z = 5.7 play the role of a baseline for the properties of LAEs with no cosmic reionization effects given by the IGM Lyα damping wing absorption. The comparisons of z = 5.7 and 6.6 LAEs over the 26 deg2 sky provide large statistical results, while the z = 7.3 LAEs in the moderately large 3.5 deg2 area allows us to investigate galaxy formation and cosmic reionization at the heart of the EoR. Based on the data sets of the HSC-SSP survey (see Aihara et al. 2018 for the details of the first-data release) and our extensive spectroscopic follow-up observations with Subaru, Keck, and Magellan telescopes, we start the program named SILVERRUSH (section 1). This program is one of a set of twin programs; the other is is a study for dropouts, called the Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH), and is detailed in Ono et al. (2018), Harikane et al. (2018), and Toshikawa et al. (2018). In a series of SILVERRUSH papers, we present the SILVERRUSH program strategy (this paper), LAE sample selections (Shibuya et al. 2018a) and spectroscopic follow-up observations (Shibuya et al. 2018b), Lyα LF evolution (Konno et al. 2018), LAE clustering evolution (this paper), and comparisons with numerical simulation results (A. K. Inoue et al. in preparation). These are early papers of the SILVERRUSH program, the observations of which are still underway. In these early papers, we present the results of the LAE studies based on the images taken until 2016 April that include neither NB387 (for z = 2.2 LAEs) nor NB101 (for z = 7.3 LAEs) data, but only NB816 and NB921 data for z = 5.7 and 6.6 LAEs, respectively (section 3). Here, our aim is that SILVERRUSH results would serve as the baselines of LAE properties useful for the ongoing LAE studies including MUSE, HETDEX, and the forthcoming z > 7 galaxy studies of the James Webb Space Telescope (JWST) and extremely large telescopes (ELTs). 3 Observations and sample The HSC-SSP survey started observations in 2014 March, taking deep broad- and narrow-band images on large areas of the sky (PI: S. Miyazaki). In our study, we use the HSC-SSP survey (HSC et al. in preparation) data taken until 2016 April with two narrow-bands (NB816 and NB921) and five broad-bands (grizy) useful for our LAE studies. Because the data of NB387 and NB101 (for detections of LAEs at z = 2.2 and 7.3) were not taken until 2016 April, we present the results of imaging with the two narrow-band NB816 and NB921 (for LAEs at z = 5.7 and 6.6), the observations of which have been partly conducted. The data are reduced with the HSC pipeline software, hscPipe version 4.0.2 (Bosch et al. 2018). The 5σ depths of the imaging data are typically ≃ 25–25.5 and ≃ 26–27 mag in narrow- and broad-bands, respectively (see Shibuya et al. 2018a for more details). Note that the HSC-SSP survey data set used in this study is notably larger than that of the first-data release (Aihara et al. 2018), which is composed of images taken only in 2014 March and 2015 November. Our LAEs are selected with the combinations of broad- and narrow-band colors down to the 5σ detection limits. The LAEs should have a narrow-band excess, the existence of a Gunn–Peterson trough, and no detection of blue continuum fluxes. The color selection criteria are similar to those of Ouchi et al. (2008, 2010), which are defined as   \begin{eqnarray} i-{\textit NB}{\textit 8}{\textit 1}{\textit 6}\ge 1.2\ {\rm and}\ g>g_{3\sigma }\ {\rm and}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \nonumber \\ \ \ \ \ \ \left[(r\le r_{3\sigma }\ {\rm and}\ r-i\ge 1.0)\ {\rm or}\ (r>r_{3\sigma })\right] \end{eqnarray} (1)and   \begin{eqnarray} z-{\textit NB}{\textit 9}{\textit 2}{\textit 1}\ge 1.0\ {\rm and}\ g>g_{3\sigma }\ {\rm and}\ r>r_{3\sigma }\ {\rm and}\ \nonumber \\ \ \ \ \ \ \left[(z\le z_{3\sigma }\ {\rm and}\ i-z\ge 1.3)\ {\rm or}\ (z>z_{3\sigma })\right] \end{eqnarray} (2)for z = 5.7 and 6.6 LAEs, respectively, where g3σ (r3σ) is the 3σ limiting magnitude of the g (r) band that ensures no detection of a continuum bluer than the Lyman break. Similarly, z3σ is the 3σ detection limit of the z band. After applying the candidate screening on the basis of hscPipe parameters+flags and visual inspection, we obtain a total of 2354 LAEs (1081 and 1273 LAEs at z = 5.7 and 6.6, respectively) from the LAE ALL catalogs (Shibuya et al. 2018a). We investigate the 2354 LAEs, and make homogeneous samples over the survey areas that consist of 1832 LAEs (959 and 873 LAEs at z = 5.7 and 6.6, respectively) with a common narrow-band limiting magnitude of NB816 < 25.0 or NB921 < 25.0 that is surely brighter than the 5σ detection levels. The homogenous samples of z = 5.7 and 6.6 LAEs have very similar threshold luminosities (the rest-frame equivalent width) of $$L_{\rm Ly\alpha }^{\rm th}>6.3 \times 10^{42}$$ and 7.9 × 1042 erg s−1 ($$\mathit {EW}_0 \gtrsim 20\,\,$$Å) at z = 5.7 and 6.6, respectively. The threshold luminosities correspond to ≃ L* luminosities (e.g., Ouchi et al. 2008, 2010; Konno et al. 2018). The LAE sample selection is detailed in Shibuya et al. (2018a). These LAEs are found in an area totalling 13.8 (21.2) deg2 consisting of the fields named D-ELAIS-N1, D-DEEP2-3, UD-COSMOS, and UD-SXDS for the z = 5.7 LAE sample (D-ELAIS-N1, D-DEEP2-3, D-COSMOS, UD-COSMOS, and UD-SXDS for the z = 6.6 LAE sample). Because, for the z = 5.7 and 6.6 LAE samples, the redshift ranges are z = 5.726 ± 0.046 and z = 6.580 ± 0.056 in the case that top-hat selection functions of LAEs with the FWHMs of the narrow-band filters are assumed, the total survey volumes covered with the narrow-band transmission are 1.2 × 107 and 1.9 × 107 comoving Mpc3 at z = 5.7 and 6.6, respectively. Note that the total areas of 13.8 and 21.2 deg2 for the z = 5.7 and 6.6 samples correspond to 0.3 and 0.5 comoving Gpc2 areas, respectively. These large survey areas allow us to study average properties of LAEs at z = 6–7, observabilities of which are influenced by patchy (10–100 Mpc) ionized bubbles probably existing at the end of the EoR (Furlanetto et al. 2006). In our LAE samples, a total of 97 LAEs at z = 5.7 and 6.6 are confirmed by our spectroscopic observations and cross-matching the existing spectra (Shibuya et al. 2018b). Because there are 81 LAE candidates that have spectroscopic identifications that are obtained via our past and present programs, we estimate the contamination rates of our LAE samples using these 81 LAE candidates, which consist of 53 and 28 LAEs at z = 5.7 and 6.6, respectively. We find that four out of 53 (four out of 28) candidates are foreground contamination objects, and estimate the contamination rates to be ∼8% and ∼14% for the samples of LAEs at z = 5.7 and 6.6, respectively. We also investigate contamination rates of bright LAEs that are brighter than 24 magn in the narrow-band. There are six and 12 bright LAE candidates with spectroscopic identifications. The spectroscopic results indicate that one out of six (four out of 12) LAE candidates are foreground interlopers, which corresponds to a contamination rate of ∼17% (∼33%) for the sample of LAEs at z = 5.7 (6.6). Although the contamination rates depend on magnitude, the contamination rates indicated by the spectroscopic confirmation range around 0%–30% in our z = 5.7 and 6.6 LAE samples. Although the contamination rates include large uncertainties due to the small number of our spectroscopically confirmed sources at this early stage of our LAE studies, we assume that the contamination rates are 0%–30% in the LAE samples that are used in our analysis below. 4 Clustering analysis and results Figures 4 and 5 present the sky distribution of the LAEs at z = 5.7 and 6.6, respectively. In figure 4, we find that the density of z = 5.7 LAEs in the D-ELAIS-N1 field is largely biased towards the southern region, indicating the possibilities of the detection incompleteness and the real large overdensity. We have carefully investigated our D-ELAIS-N1 images, but found no significant spatial fluctuations in the detection completeness produced by the inhomogeneity of the data handling and the data qualities, including depths and seeing sizes. However, there remains a possibility that unknown effects of the data reduction may produce a large overdensity. Because the origin of this large overdensity is still unknown, we do not use the data of D-ELAIS-N1 in this study. Instead, we discuss this overdensity in R. Higuchi et al. (in preparation). We thus use the remaining 734 LAEs at z = 5.7 in our analysis. Fig. 4. View largeDownload slide Sky distribution of the LAEs at z = 5.7. The red squares, magenta diamonds, and black circles represent positions of narrow-band bright (<23.5 mag), medium-bright (23.5–24.0 mag), and faint (24.0–25.0 mag) LAEs, respectively. The large red open square indicates the LAEs with spatially extended Lyα emission (Shibuya et al. 2018b). The gray shades denote either the areas with no HSC data or the masked regions with a bad data quality. The scale on the map is marked in angles (degrees) and the projected distances (comoving megaparsecs). (Color online) Fig. 4. View largeDownload slide Sky distribution of the LAEs at z = 5.7. The red squares, magenta diamonds, and black circles represent positions of narrow-band bright (<23.5 mag), medium-bright (23.5–24.0 mag), and faint (24.0–25.0 mag) LAEs, respectively. The large red open square indicates the LAEs with spatially extended Lyα emission (Shibuya et al. 2018b). The gray shades denote either the areas with no HSC data or the masked regions with a bad data quality. The scale on the map is marked in angles (degrees) and the projected distances (comoving megaparsecs). (Color online) Fig. 5. View largeDownload slide Same as figure 4, but for the LAEs z = 6.6. The large red open squares indicate the LAEs with spatially extended Lyα emission including Himiko (Ouchi et al. 2009) and CR7 (Sobral et al. 2015). See Shibuya et al. (2018b) for more details. (Color online) Fig. 5. View largeDownload slide Same as figure 4, but for the LAEs z = 6.6. The large red open squares indicate the LAEs with spatially extended Lyα emission including Himiko (Ouchi et al. 2009) and CR7 (Sobral et al. 2015). See Shibuya et al. (2018b) for more details. (Color online) We quantify clustering properties based on the measurements of the angular correlation functions in the following sections. 4.1 Angular correlation function We derive angular two-point correlation functions (ACFs) of our LAEs in the same manner as Ouchi et al. (2003, 2005b, 2010) and Harikane et al. (2016). We use the estimator of Landy and Szalay (1993),   $$\omega _{\rm obs}(\theta ) = \frac{DD(\theta )-2DR(\theta )+RR(\theta )}{RR(\theta )},$$ (3)where DD(θ), DR(θ), and RR(θ) are the numbers of galaxy–galaxy, galaxy–random, and random–random pairs normalized by the total number of pairs in each of the samples of the pairs. We use a random catalog (Coupon et al. 2018) that has a surface number density of 100 arcmin−2. The random catalog has the same geometrical constraint as the distribution of our LAEs, representing our survey areas. Statistical errors are estimated using Jackknife resampling with subsamples that each have an area of ∼10002 arcsec2. Figure 6 shows the observed ACFs ωobs(θ) of our LAEs at z = 5.7 and 6.6. The ACF measurements cover the scale of ∼1–100 comoving Mpc that is indicated in the upper axes of figure 6. We find that these ACFs are consistent with those obtained with the previous Subaru/Suprime-Cam data (Ouchi et al. 2010), and that our present HSC data provide the large-scale ACFs with uncertainties smaller than those of the previous data. Fig. 6. View largeDownload slide Angular correlation function (ACF) and bias of the LAEs at z = 5.7 (left) and 6.6 (right). The top and bottom panels show the ACFs and the bias. The black-filled and gray-open squares represent the ACFs of our LAEs with and without the IC correction, respectively, while the filled circles denote the ACFs with the IC correction derived by Ouchi et al. (2010). For presentation purposes, we slightly shift the gray open squares along the abscissa. In the top panels, the solid and dotted lines present the best-fitting power-law functions of our ACFs and the ACFs of the underlying dark matter predicted by the linear theory (e.g., Peacock & Dodds 1994). Because the power-law spatial correlation function is projected on sky with the method in Simon (2007), the best-fitting power-law functions are curves. In the bottom panels, the solid and dashed lines indicate the average bias and the 1σ error values, respectively. The crosses in the top left-hand panel show the ACFs obtained by Murayama et al. (2007). The top axis denotes the projected distance in comoving megaparsecs. Fig. 6. View largeDownload slide Angular correlation function (ACF) and bias of the LAEs at z = 5.7 (left) and 6.6 (right). The top and bottom panels show the ACFs and the bias. The black-filled and gray-open squares represent the ACFs of our LAEs with and without the IC correction, respectively, while the filled circles denote the ACFs with the IC correction derived by Ouchi et al. (2010). For presentation purposes, we slightly shift the gray open squares along the abscissa. In the top panels, the solid and dotted lines present the best-fitting power-law functions of our ACFs and the ACFs of the underlying dark matter predicted by the linear theory (e.g., Peacock & Dodds 1994). Because the power-law spatial correlation function is projected on sky with the method in Simon (2007), the best-fitting power-law functions are curves. In the bottom panels, the solid and dashed lines indicate the average bias and the 1σ error values, respectively. The crosses in the top left-hand panel show the ACFs obtained by Murayama et al. (2007). The top axis denotes the projected distance in comoving megaparsecs. Although our survey area is large, we evaluate the integral constraint (IC) (Groth & Peebles 1977). The IC value corresponds to the observational offset in ωobs(θ) originated by a limited survey area. Including the correction for the number of objects in the sample, N, the true ACF is given by   $$\omega (\theta )=\omega _\mathrm{obs}(\theta )+\mathrm{IC}+\frac{1}{N}.$$ (4)We evaluate the integral constraint with   $$\mathrm{IC}=\frac{\Sigma _i RR(\theta _i)\omega (\theta _i)}{\Sigma _i RR(\theta _i)}.$$ (5)We adopt the model ACF of Simon (2007) for the function of ω(θ) that is detailed in subsection 4.2. We show IC values in table 1, and plot the ACFs corrected for IC in figure 6. Note that the IC values are very small, because our survey areas are very large. Table 1. LAE samples and the clustering measurements. z  $$L_{\rm Ly\alpha }^{\rm th}$$  ng  N  IC  r0  bavg  $$r_0^{\rm max}$$*  $$b_{\rm avg}^{\rm max}$$  χ2/dof    (1042 erg s−1)  (10−5 Mpc−3)    (10−3)  ($$h_{100}^{-1}\:$$Mpc)    ($$h_{100}^{-1}\:$$Mpc)      (1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10)  5.7  6.3  7.3  959†  1.83  $$3.01^{+0.35}_{-0.35}$$  4.13 ± 0.17  $$4.47^{+0.52}_{-0.52}$$  5.90 ± 0.24  5.69/5  6.6  7.9  2.4  873  1.35  $$2.66^{+0.49}_{-0.70}$$  4.54 ± 0.63  $$3.95^{+0.73}_{-1.04}$$  6.49 ± 0.90  2.70/2  z  $$L_{\rm Ly\alpha }^{\rm th}$$  ng  N  IC  r0  bavg  $$r_0^{\rm max}$$*  $$b_{\rm avg}^{\rm max}$$  χ2/dof    (1042 erg s−1)  (10−5 Mpc−3)    (10−3)  ($$h_{100}^{-1}\:$$Mpc)    ($$h_{100}^{-1}\:$$Mpc)      (1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10)  5.7  6.3  7.3  959†  1.83  $$3.01^{+0.35}_{-0.35}$$  4.13 ± 0.17  $$4.47^{+0.52}_{-0.52}$$  5.90 ± 0.24  5.69/5  6.6  7.9  2.4  873  1.35  $$2.66^{+0.49}_{-0.70}$$  4.54 ± 0.63  $$3.95^{+0.73}_{-1.04}$$  6.49 ± 0.90  2.70/2  * For comparison, r0 is expressed with h100. † This number includes sources in D-ELAIS-N1. View Large A clustering signal of galaxies is diluted by contamination objects in a galaxy sample. If the galaxy sample includes randomly distributed contamination objects with a fraction fc, the value of the ACF is reduced by a factor of (1 − fc)2. The ACF corrected for the randomly-distributed contamination ωmax is written as   $$\omega _{\rm max} = \frac{\omega }{(1-f_{\rm c})^2}.$$ (6)This is the maximum reduction of the ACF, because the real contamination objects are not randomly distributed but spatially correlated. For reference, we evaluate the possible maximum values of ωmax, using a contamination fraction with an upper limit of fc = 0.3 (section 3). 4.2 Correlation length and bias We fit the ACFs with a simple power-law model of the spatial correlation function,   $$\xi (r)=\left(\frac{r}{r_0}\right)^{-\gamma },$$ (7)where γ, r0, and r are the slope of the power law, the correlation length, and the spatial separation between two galaxies, respectively. The spatial correlation function is related to the ACF by the Limber equation (Peebles 1980; Efstathiou et al. 1991), which is an integral equation of the (three-dimensional) spatial correlation function connecting with the (two-dimensional) ACF. However, Simon (2007) claims that the Limber equation does not provide accurate values in a very large separation of galaxies that have redshift-distribution distances narrower than the transverse distance in the case of narrow-band-selected LAEs. We thus adopt the method that Simon (2007) suggests, and derive r0 of our LAEs, fitting the power-law functions to our ACFs. Because no meaningful constraints on γ are obtained with the present samples, we adopt the fiducial γ value of γ = 1.8 that is adopted in the previous clustering analyses for LAEs (Ouchi et al. 2003; Gawiser et al. 2007; Kovač et al. 2007; Ouchi et al. 2010). To investigate the dependencies of our results on the value of γ, we also use the other fixed γ values of γ = 1.6 and 2.0 that bracket the possible range of the power-law index of high-z galaxies at z ∼ 4–6 (Lee et al. 2006; McLure et al. 2009). We have found that neither r0 nor bias values change by over 10%, and that r0 and bias values fall well within the errors. Here, we use the redshift distribution of our LAEs estimated with the narrow-band transmission curve. The best-fitting power-law functions are presented in figure 6. The correlation lengths thus obtained are $$r_0=3.01^{+0.35}_{-0.35}$$ and $$2.66^{+0.49}_{-0.70}$$$$h^{-1}_{100}\:$$Mpc for LAEs at z = 5.7 and 6.6, respectively, as summarized in table 1. To make comparisons with the previous results, we express the correlation lengths using h100, the Hubble constant in units of 100 km s−1, instead of 70 km s−1. Table 1 also shows the upper limit of the correlation length, $$r_0^{\rm max}$$, which is given by   $$r_0^{\rm max}=r_0\left(\frac{1}{1-f_\mathrm{C}}\right)^{2/\gamma }.$$ (8)We obtain $$r_0^{\rm max} = 4.47^{+0.52}_{-0.52}$$ and $$3.95^{+0.73}_{-1.04}\,\,h_{100}^{-1}\:$$Mpc for our z = 5.7 and 6.6 LAEs, respectively. These correlation length values are consistent with the previous measurements within the errors. Ouchi et al. (2010) obtain (r0, $$r_0^{\rm max}$$) of ($$3.12^{+0.33}_{-0.36}$$, $$4.30^{+0.45}_{-0.50}$$) and ($$2.31^{+0.65}_{-0.85}$$, $$3.60^{+1.02}_{-1.32}$$) for LAEs at z = 5.7 and 6.6, respectively. In the framework of the ΛCDM model, the galaxy dark-matter bias is defined by   $$b_{\rm g}(\theta )^2 \equiv \omega (\theta )/\omega _{\rm dm}(\theta ),$$ (9)where ωdm(θ) is the ACF of dark matter predicted by the linear theory model (e.g., Peacock & Dodds 1994). The model calculations are made with the same survey volumes as those of our LAE samples. Note that bg is the bias value equivalent to the one given by the three-dimensional spatial correlation functions, $$b_{\rm g}^2=\xi (r)/\xi _{\rm dm}(r)$$, where ξdm(r) is the spatial correlation function of dark matter. In figure 6, the top and bottom panels show ωdm(θ) and bg(θ), respectively. We estimate the average bias bavg by averaging bg(θ) with error weighting over the angular range presented in figure 6. We also calculate the upper limits of bias values, $$b_{\rm avg}^{\rm max}$$, with the maximal contamination correction that corresponds to $$A_{\omega }^{\rm max}$$ [equation (6); see subsection 4.1]. We obtain bavg = 4.13 ± 0.17 and 4.54 ± 0.63 ($$b_{\rm avg}^{\rm max} = 5.90\,\,\pm \,\, 0.24$$ and 6.49 ± 0.90) for our z = 5.7 and 6.6 LAEs, respectively. We summarize bavg and $$b_{\rm avg}^{\rm max}$$ values thus obtained in table 1. We compare these bias values with those obtained in the previous studies. We find that the bias value of our LAEs are consistent with those derived by Ouchi et al. (2005a, 2010) within the ∼1σ errors. Figure 7 presents bias of LAEs as a function of redshift. In figure 7, we plot bias measurements of LAEs at z = 2–7 derived in the previous studies. We correct these previous study measurements of bg for the difference of σ8 (Ouchi et al. 2010), and show the corrected bg values in figure 7. Fig. 7. View largeDownload slide Bias of ≳ L* LAEs as a function of redshift. The red squares represent bias of our LAEs at z = 5.7 and 6.6, while the filled circles denote those of LAEs obtained in Ouchi et al. (2010). The open circles, pentagon, and hexagon indicate LAE bias values obtained by Ouchi et al. (2010, 2005a, 2003), respectively. The filled hexagon, pentagon, and diamond are LAE bias values measured by Guaita et al. (2010), Gawiser et al. (2007), and Kovač et al. (2007), respectively. For the previous study results, we correct for the difference of σ8. To clarify the bias measurements, we give small offsets along the abscissa axis to the data points of the previous studies. The solid lines indicate bias of dark-matter halos with a halo mass of 108, 109, 1010, 1011, and 1012 M⊙ in the case of one-to-one correspondence between galaxies and dark-matter halos (Ouchi et al. 2010). The gray region shows the dark-matter halo mass range of 1010–1012 M⊙. The dotted lines are evolutionary tracks of bias in the case of the galaxy-conserving model [equation (13)]. (Color online) Fig. 7. View largeDownload slide Bias of ≳ L* LAEs as a function of redshift. The red squares represent bias of our LAEs at z = 5.7 and 6.6, while the filled circles denote those of LAEs obtained in Ouchi et al. (2010). The open circles, pentagon, and hexagon indicate LAE bias values obtained by Ouchi et al. (2010, 2005a, 2003), respectively. The filled hexagon, pentagon, and diamond are LAE bias values measured by Guaita et al. (2010), Gawiser et al. (2007), and Kovač et al. (2007), respectively. For the previous study results, we correct for the difference of σ8. To clarify the bias measurements, we give small offsets along the abscissa axis to the data points of the previous studies. The solid lines indicate bias of dark-matter halos with a halo mass of 108, 109, 1010, 1011, and 1012 M⊙ in the case of one-to-one correspondence between galaxies and dark-matter halos (Ouchi et al. 2010). The gray region shows the dark-matter halo mass range of 1010–1012 M⊙. The dotted lines are evolutionary tracks of bias in the case of the galaxy-conserving model [equation (13)]. (Color online) This work and the previous studies measure bias of LAEs that have similar Lyα luminosity detection limits; a few times 1042 erg s−1 corresponding to ∼L* luminosities over z ∼ 2–7 (e.g., Ouchi et al. 2008, 2010; Konno et al. 2016, 2018). Thus, we can omit the difference of bias depending on Lyα luminosity. Figure 7 indicates that the bias of ≳ L* LAEs significantly increases from z ∼ 2–3 to z ∼ 6–7. The physical origins of this increase are discussed in section 5. 5 Discussion 5.1 Cosmic reionization 5.1.1 Constraints on xH i from the clustering results Clustering signals of observed LAEs are strengthened by the patchy distribution of the neutral hydrogen at the EoR, because Lyα photons of LAEs in the ionized bubbles can selectively escape from the patchy neutral IGM with a small amount of Lyα damping wing absorption (Furlanetto et al. 2006; McQuinn et al. 2007; Lidz et al. 2009; Iliev et al. 2008; Sobacchi & Mesinger 2015). Based on the fact that the Gunn–Peterson optical depth largely increases at z ∼ 6 (Fan et al. 2006), one would expect that the clustering amplitude of the observed LAEs increases from z = 5.7 (the ionized universe) to z = 6.6 (the partly-neutral universe). Figure 7 indicates no significant rise of bias from z = 5.7 to 6.6 beyond the moderately small errors (i.e., by a factor of ∼20%). This result suggests that the clustering of z = 6.6 LAEs is not largely affected by cosmic reionization effects, where the bias evolution of the hosting dark-matter halos towards high z may also be involved. Based on this bias evolution result, we place constraints on cosmic reionization parameters with the help of theoretical models. We compare our observational results with multiple theoretical models, because there is a possibility that our conclusions may depend on the model chosen for the comparisons. We thus try to avoid model-dependent conclusions as much as possible, and to obtain more objective interpretations for our observational results. Note that the arguments below follow those of Ouchi et al. (2010), with our HSC clustering measurements. In the top panel of figure 8, we compare our z = 6.6 LAE clustering measurements with those of theoretical predictions (McQuinn et al. 2007; Furlanetto et al. 2006). The model of McQuinn et al. (2007) is presented in the top panel of figure 8. McQuinn et al. (2007) conducted radiative transfer simulations and predicted clustering of LAEs at z = 6.6. Their models assign a Lyα flux to a dark-matter halo whose mass is greater than a minimum dark-matter halo mass. We use the minimum dark-matter halo mass of our LAEs at the post-EoR of z = 5.7 to estimate the intrinsic LAE clustering at z = 6.6, assuming no redshift evolution (z = 5.7–6.6) of the minimum dark-matter halo mass (section 3). Although the minimum dark-matter halo mass is 3 × 109 M⊙ for the z = 5.7 LAEs (subsection 5.2), McQuinn et al. (2007) do not calculate the models of the minimum dark-matter halo mass as low as 3 × 109 M⊙. We thus correct the angular correlation functions of McQuinn et al.’s lowest mass (3 × 1010 M⊙) model for the difference between 3 × 1010 and 3 × 109 M⊙, using equation (9). Here we calculate the bias values of the 3 × 1010 and 3 × 109 M⊙ cases using the best-fitting halo model (subsection 5.2), the parameters of which are exactly the same as those determined at z = 5.7. In the top panel of figure 8, we show McQuinn et al.’s models for 3 × 109 M⊙ with the correction, together with the original 3 × 1010 M⊙ model. Because McQuinn et al.’s models do not correct for the integral constraint, we do not compare the model predictions below ω ∼ 0.1, where the contribution of the integral-constraint correction term is large. The top panel of figure 8 indicates that our z = 6.6 LAE data points fall in the range of xH i ≃ 0–0.3, indicating xH i = 0.15 ± 0.15. Fig. 8. View largeDownload slide Comparisons of the theoretical models with the ACF and bias measurements at z = 6.6. In the top (bottom) panel, the filled squares and circles denote the ACF (bias) measurements with the IC correction obtained by this work and Ouchi et al. (2010), respectively, which are the same as those shown in the top (bottom) right-hand panel of figure 6. Top: The solid and dotted lines represent McQuinn et al.’s (2007) models of z = 6.6 LAEs with a dark-matter halo masses of 3 × 109 M⊙ and 3 × 1010 M⊙, respectively. The model of 3 × 109 M⊙ is the modification of the model of 3 × 1010 M⊙ (see text). From the bottom to top, the solid (dotted) lines, showing the neutral hydrogen fractions of the IGM, are xH i = 0.0, 0.3, 0.5, and 0.8. Bottom: The ticks at the right-hand side indicate bias values predicted by Furlanetto, Zaldarriaga, and Hernquist’s (2006) models (see text). The solid and dashed lines indicate the average bias and the 1σ uncertainty values, respectively, that are the same as those shown in the bottom right-hand panel of figure 6. Fig. 8. View largeDownload slide Comparisons of the theoretical models with the ACF and bias measurements at z = 6.6. In the top (bottom) panel, the filled squares and circles denote the ACF (bias) measurements with the IC correction obtained by this work and Ouchi et al. (2010), respectively, which are the same as those shown in the top (bottom) right-hand panel of figure 6. Top: The solid and dotted lines represent McQuinn et al.’s (2007) models of z = 6.6 LAEs with a dark-matter halo masses of 3 × 109 M⊙ and 3 × 1010 M⊙, respectively. The model of 3 × 109 M⊙ is the modification of the model of 3 × 1010 M⊙ (see text). From the bottom to top, the solid (dotted) lines, showing the neutral hydrogen fractions of the IGM, are xH i = 0.0, 0.3, 0.5, and 0.8. Bottom: The ticks at the right-hand side indicate bias values predicted by Furlanetto, Zaldarriaga, and Hernquist’s (2006) models (see text). The solid and dashed lines indicate the average bias and the 1σ uncertainty values, respectively, that are the same as those shown in the bottom right-hand panel of figure 6. In the bottom panel of figure 8, we compare our bias results with those of the analytical models of Furlanetto, Zaldarriaga, and Hernquist (2006) at z = 6.6. In this comparison, we adopt the bias values of the small scale (≃1−10 Mpc) of Furlanetto, Zaldarriaga, and Hernquist (2006) that is similar to the major angular scale of our bias measurements. For the intrinsic bias value of the z = 6.6 LAEs, we use b = 4.6, which is obtained under the assumption of no evolution from z = 5.7 to 6.6 for the minimum dark-matter halo mass (3 × 109 M⊙) and the halo-model parameters (subsection 5.2). The bottom panel of figure 8 shows that our average-bias measurement agrees with the low-xH i models of Furlanetto, Zaldarriaga, and Hernquist (2006), and suggests xH i ≲ 0.13 including the uncertainty of the average-bias measurement of the z = 6.6 LAEs. The comparisons of our observational results with the models of McQuinn et al. (2007) and Furlanetto, Zaldarriaga, and Hernquist (2006) indicate that the neutral hydrogen fraction at z = 6.6 falls in the range of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$. 5.1.2 Cosmic reionization history Figure 9 shows our xH i value of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ that is estimated with our LAE clustering measurements. This xH i estimate suggests a moderate neutral hydrogen fraction at z = 6.6. In figure 9, we also present various results of xH i constrained by a similar LAE clustering analysis (Ouchi et al. 2010) and the independent xH i estimates from QSO Gunn–Peterson absorptions (Fan et al. 2006), as well as Lyα damping wing (DW) absorption measurements of QSOs (Mortlock et al. 2011; Bolton et al. 2011), gamma-ray bursts (GRBs; Totani et al. 2006, 2014, 2016; see also Greiner et al. 2009), and galaxies including LAEs (Malhotra & Rhoads 2004; Ouchi et al. 2010; Ota et al. 2010; Kashikawa et al. 2011; Pentericci et al. 2011; Ono et al. 2012; Treu et al. 2013; Caruana et al. 2014; Konno et al. 2014, 2018; Schenker et al. 2014). Our results agree with the previous LAE clustering results of Ouchi et al. (2010). Our constraint on xH i is stronger than that of Ouchi et al. (2010), due to the fact that the statistical uncertainties of the ACF measurements of this study are smaller than those of Ouchi et al. (2010). This is because the number of LAEs in this study is larger than that of Ouchi et al. (2010). Moreover, our xH i estimate agrees with many of the estimates and constraints obtained by the other independent objects (QSOs, GRBs, and galaxies) and methods (Gunn–Peterson and Lyα DW absorptions) in the previous studies. It should be noted that the results for GRBs have a large scatter. This scatter is probably made by the systematic uncertainty raised by sightline variance effects. The numerical simulations of McQuinn et al. (2008) suggest that the patchiness of reionization gives the systematic error of δxH i ∼ 0.3 in an xH i measurement of a single GRB sightline. Fig. 9. View largeDownload slide Neutral hydrogen fraction xH i of the IGM as a function of redshift. Note that the top and bottom panels are the same, but with an ordinate axis of linear and log scales, respectively. The red filled square is the xH i estimate obtained by our HSC LAE clustering analysis. The black filled square and circle are the xH i estimates from the LAE LF evolution of Konno et al. (2018, 2014), respectively. The open circles are the constraints at z = 6.6 obtained by Ouchi et al. (2010) from the evolution of Lyα LF (left circle) and clustering (right circle), while the open diamond and the open pentagon represent the upper limits from the Lyα LF evolution to z = 6.5 given by Malhotra and Rhoads (2004) and Kashikawa et al. (2006). Here, we add small offsets along redshift to the positions of the filled square, the open circles, and the open diamond, avoiding overlapping symbols. The filled hexagon and the filled pentagons show the constraints from a spectrum of a GRB (Gallerani et al. 2008b) and statistics of QSO dark-gaps (Gallerani et al. 2008a), respectively. The open hexagons are the constraints calculated from the Lyα damping wing absorption of GRBs at z = 6.3 (Totani et al. 2006) and z = 5.9 (Totani et al. 2016). The filled diamonds indicate the QSO Gunn–Peterson optical depth measurement results (Fan et al. 2006). The triangle denotes the 1σ lower-limit of redshift obtained by Planck 2015 (Planck Collaboration 2016) in the case of instantaneous reionization. The solid line and the gray shade indicate the best estimate and the uncertainty of the xH i evolution (Ishigaki et al. 2017) that agrees the evolutions of τe and ρUV with free parameters including the ionizing photon escape fraction. The dotted, dashed, and dot–dashed lines are the evolution of xH i for the reionizing sources down to the massive halos, the moderately massive halos, and the mini-halos, respectively, in the model of Choudhury, Ferrara, and Gallerani (2008). The dashed–double-dotted line indicates the prediction of the double reionization model (Cen 2003). (Color online) Fig. 9. View largeDownload slide Neutral hydrogen fraction xH i of the IGM as a function of redshift. Note that the top and bottom panels are the same, but with an ordinate axis of linear and log scales, respectively. The red filled square is the xH i estimate obtained by our HSC LAE clustering analysis. The black filled square and circle are the xH i estimates from the LAE LF evolution of Konno et al. (2018, 2014), respectively. The open circles are the constraints at z = 6.6 obtained by Ouchi et al. (2010) from the evolution of Lyα LF (left circle) and clustering (right circle), while the open diamond and the open pentagon represent the upper limits from the Lyα LF evolution to z = 6.5 given by Malhotra and Rhoads (2004) and Kashikawa et al. (2006). Here, we add small offsets along redshift to the positions of the filled square, the open circles, and the open diamond, avoiding overlapping symbols. The filled hexagon and the filled pentagons show the constraints from a spectrum of a GRB (Gallerani et al. 2008b) and statistics of QSO dark-gaps (Gallerani et al. 2008a), respectively. The open hexagons are the constraints calculated from the Lyα damping wing absorption of GRBs at z = 6.3 (Totani et al. 2006) and z = 5.9 (Totani et al. 2016). The filled diamonds indicate the QSO Gunn–Peterson optical depth measurement results (Fan et al. 2006). The triangle denotes the 1σ lower-limit of redshift obtained by Planck 2015 (Planck Collaboration 2016) in the case of instantaneous reionization. The solid line and the gray shade indicate the best estimate and the uncertainty of the xH i evolution (Ishigaki et al. 2017) that agrees the evolutions of τe and ρUV with free parameters including the ionizing photon escape fraction. The dotted, dashed, and dot–dashed lines are the evolution of xH i for the reionizing sources down to the massive halos, the moderately massive halos, and the mini-halos, respectively, in the model of Choudhury, Ferrara, and Gallerani (2008). The dashed–double-dotted line indicates the prediction of the double reionization model (Cen 2003). (Color online) Figure 9 presents the evolution of xH i as determined by Ishigaki et al. (2017), who use all of the observational results related to reionization, CMB Thomson scattering optical depth, τe, UV luminosity density, ρUV, and the ionized fraction, QH ii, obtained to-date, based on the standard analytical model (Madau et al. 1999; Robertson et al. 2015) fitting that allows a wide parameter range of unknown parameters such as the ionizing-photon escape fraction. It should be noted that Ishigaki et al.’s (2017) result of xH i evolution agrees with the τe measurement of Planck2016 (Planck Collaboration 2016). We find in figure 9 that the evolution of xH i (the solid curve) is consistent with our xH i estimate from the LAE clustering, supportng the moderately late reionization scenario suggested by Ishigaki et al. (2017). 5.2 Hosting dark-matter halos 5.2.1 Halo occupation distribution modeling We carry out halo occupation distribution (HOD) modeling for our ACFs where the hosting dark-matter halo properties are encoded. Note that the ACF of the LAEs at z = 5.7 is not affected by the effects of patchy ionized bubbles of cosmic reionization (section 1), because z = 5.7 corresponds to the post-EoR. However, the LAE ACF at z = 6.6 may be strengthened by the selective escape of Lyα photons from patchy ionized bubbles that probably exist at this redshift, although the z = 6.6 LAE ACF appears to be weakly affected by cosmic reionization (sub-subsection 5.1.1). In this section, we perform HOD modeling for the ACF of the LAEs at z = 5.7 for the secure results. Moreover, we conduct the same analysis for the LAEs at z = 6.6 under the assumption that the reionization effect is negligibly small, following our results that the influence of cosmic reionization is not large at z = 6.6 (sub-subsection 5.1.1). We adopt the HOD models of Harikane et al. (2016) based on the halo mass function of Behroozi, Wechsler, and Conroy (2013) that is a modification of the Tinker et al. (2008) halo mass function. In the HOD models of Harikane et al. (2016), there are a total of five parameters: Mmin, $$M_1^{\prime }$$, α, M0, and $$\sigma _{\rm \log M}$$. Here, Mmin is the mass of the dark-matter halo that hosts a galaxy at a possibility of 50%. The terms $$M_1^{\prime }$$ and α represent the normalization and the slope of the power law, respectively, for the occupation number of the satellite galaxy. We define M0 and σlog  M as the halo mass cut-off and the halo mass transition width, respectively, of the central and satellite galaxies. We do not include the parameter of the duty cycle in our HOD modeling, because the number density is not used in our fitting. Instead, we compare the number densities given by our HOD modeling and our observations to constrain the duty cycle of LAEs, which are detailed in sub-subsection 5.2.2. Following Harikane et al. (2016), we adopt the two relations $$\log M_0 = 0.76 \log M_1^{\prime } + 2.3$$ and $$\log M_1^{\prime }=1.18 \log M_{\rm min} - 1.28$$. We assume $$\sigma _{\rm \log M}=0.2$$ and α = 1, which are suggested by the fitting results of Zheng, Coil, and Zehavie (2007), although we confirm that the choice of these two parameters does not change our conclusions. We thus obtain the best-fitting HOD models (i.e., Mmin values) by χ2 minimization. Figure 10 presents the best-fitting HOD models, and table 2 summarizes the best-fitting Mmin values. We define the average dark-matter halo mass ⟨Mh⟩ and the average galaxy bias $$b_{\rm g}^{\rm HOD}$$ from the HOD modeling results;   \begin{eqnarray} \left\langle M_{\rm h} \right\rangle & = & \frac{\int _{M_{\rm h}^{\rm min}}^{\infty } M N_{\rm g} n dM}{\int _{M_{\rm h}^{\rm min}}^{\infty }N_{\rm g} n dM},\ \ \ {\rm and} \end{eqnarray} (10)  \begin{eqnarray} b_{\rm g}^{\rm HOD} & = & \frac{\int _{M_{\rm h}^{\rm min}}^{\infty } b N_{\rm g} n dM}{\int _{M_{\rm h}^{\rm min}}^{\infty } N_{\rm g} n dM}, \end{eqnarray} (11)where n and b are the number density and the bias of dark-matter halos, respectively, with dark-matter halo mass Mh (Behroozi et al. 2013; Tinker et al. 2010). Ng is the galaxy occupation function at Mh that is determined by the HOD modeling. The values of ⟨Mh⟩ and $$b_{\rm g}^{\rm HOD}$$ for our LAEs are summarized in table 2. We find that the average dark-matter halo masses of the ≳ L* LAEs are moderately small, $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) = 11.1^{+0.2}_{-0.4}$$ at z = 5.7 and $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) = 10.8^{+0.3}_{-0.5}$$ at z = 6.6, and are consistent with the previous estimates [log (⟨Mh⟩/M⊙) ∼ 10–11] found via clustering analysis (Ouchi et al. 2010). The values of $$b_{\rm g}^{\rm HOD}$$ are consistent with those of bavg (table 1) within the errors. Fig. 10. View largeDownload slide Best-fitting HOD models of the LAEs at z = 5.7 (left) and 6.6 (right). The solid lines represent the best-fitting HOD models, while the dashed lines denote the breakdowns of the best-fitting HOD models, 1 and 2 halo terms. The squares and circles are the ACF measurements that are the same as those in figure 6. Fig. 10. View largeDownload slide Best-fitting HOD models of the LAEs at z = 5.7 (left) and 6.6 (right). The solid lines represent the best-fitting HOD models, while the dashed lines denote the breakdowns of the best-fitting HOD models, 1 and 2 halo terms. The squares and circles are the ACF measurements that are the same as those in figure 6. Table 2. Properties of dark-matter halos of the LAEs estimated with the HOD model. z  $$L_{\rm Ly\alpha }^{\rm th}$$  log (Mmin/M⊙)  log (⟨M⟩/M⊙)  $$b_{\rm g}^{\rm HOD}$$  $$f_{\rm duty}^{\rm LAE}$$  χ2/dof    (1042 ergz >s−1)        (10−4)    (1)  (2)  (3)  (4)  (5)  (6)  (7)  5.7  6.3  $$9.5^{+0.5}_{-1.2}$$  $$11.1^{+0.2}_{-0.4}$$  $$3.9^{+0.7}_{-1.0}$$  $$7.7^{+53.0}_{-7.6}$$  1.1/6  6.6*  7.9  $$9.1^{+0.7}_{-1.9}$$  $$10.8^{+0.3}_{-0.5}$$  $$4.1^{+1.0}_{-1.4}$$  $$0.9^{+13.5}_{-0.9}$$  1.1/4  z  $$L_{\rm Ly\alpha }^{\rm th}$$  log (Mmin/M⊙)  log (⟨M⟩/M⊙)  $$b_{\rm g}^{\rm HOD}$$  $$f_{\rm duty}^{\rm LAE}$$  χ2/dof    (1042 ergz >s−1)        (10−4)    (1)  (2)  (3)  (4)  (5)  (6)  (7)  5.7  6.3  $$9.5^{+0.5}_{-1.2}$$  $$11.1^{+0.2}_{-0.4}$$  $$3.9^{+0.7}_{-1.0}$$  $$7.7^{+53.0}_{-7.6}$$  1.1/6  6.6*  7.9  $$9.1^{+0.7}_{-1.9}$$  $$10.8^{+0.3}_{-0.5}$$  $$4.1^{+1.0}_{-1.4}$$  $$0.9^{+13.5}_{-0.9}$$  1.1/4  * The results of z = 6.6 LAEs may be weakly affected by the reionization effects (see text). View Large 5.2.2 Duty cycle It is suggested that LAEs do not exist in any dark-matter halos that are more massive than $$M_{\rm h}^{\rm min}$$ (e.g., Ouchi et al. 2010) because LAEs should meet two physical conditions: (i) LAEs should be active star-forming galaxies producing Lyα photons, and (ii) the Lyα photons of LAEs should largely escape from the ISM. In other words, not all galaxies in a dark-matter halo can be LAEs. We evaluate the duty cycle of LAEs (i.e., $$f_{\rm duty}^{\rm LAE}$$) that is defined by the probability that a dark-matter halo with a mass greater than the minimum-halo mass hosts an LAE(s). With fduty, one can calculate the galaxy number density,   $$n_{\rm g} = \int _{M_{\rm h}^{\rm min}}^{\infty } f_{\rm duty}^{\rm LAE}\ N_{\rm g}\ n\ dM.$$ (12)We derive $$f_{\rm duty}^{\rm LAE}$$ with the best-fitting HOD parameters, assuming that $$f_{\rm duty}^{\rm LAE}$$ does not depend on the dark-matter halo mass. The observational measurements of ng are taken from the HSC SILVERRUSH paper of Konno et al. (2018). We thus obtain $$f_{\rm duty}^{\rm LAE}=7.7^{+53.0}_{-0.75} \times 10^{-4}$$ and $$0.9^{+13.5}_{-0.9} \times 10^{-4}$$ for our LAEs at z = 5.7 and 6.6, respectively, which are listed in table 2. Although $$f_{\rm duty}^{\rm LAE}$$ is very poorly constrained with large uncertainties, we find that the duty cycle of LAEs is ∼1% or less at z = 5.7 and 6.6. These small duty cycle values for z = 6–7 LAEs are consistent with the previous estimates found via a similar technique (Ouchi et al. 2010). The small duty cycle values indicate that dark-matter halos hosting LAEs are rare, and that the majority of LAEs are located in abundant low-mass dark-matter halos with a mass close to the minimum-halo mass. 5.2.3 Galaxy formation and LAEs Figure 7 presents the evolutionary tracks of dark-matter halos for the galaxy-conserving model and the constant-mass model. For the galaxy-conserving models, we assume that gravity drives the motion of galaxies under the condition of no mergers. The bias evolution in the galaxy-conserving model (Fry 1996) is described as   $$b_{\rm g}=1+(b_{\rm g}^0-1)/D(z),$$ (13)where D(z) and $$b_{\rm g}^0$$ are the growth factor and the bias at z = 0, respectively. If one assumes the galaxy-conserving evolution, our LAEs at z = 5.7 and 6.6 evolve into the present-day galaxies with $$b_{\rm g}^0\sim 1.6$$–1.7. Galaxies with $$b_{\rm g}^0\sim 1.6$$–1.7 are massive, bright ∼6L* galaxies today (Zehavi et al. 2005). The galaxy-conserving evolution suggests that our LAEs at z = 5.7 and 6.6 are the progenitors of very massive bright galaxies in the present-day universe. 6 Summary We investigate clustering of Lyα emitters (LAEs) at z = 5.7 and 6.6 with the early data from the Hyper Suprime-Cam (HSC) Subaru Strategic Program (first-data release shown in Aihara et al. 2018), and introduce the program of Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH). From the early data, we obtain an unprecedentedly large sample of 1832 ≳ L* LAEs at z = 6–7 over the total area of 14–21 deg2, which is about an order of magnitude larger than the previous z = 6–7 LAE clustering studies. Based on the LAE clustering measurements, we study cosmic reionization and galaxy formation, comparing theoretical models. In this study, there are two major results, which are summarized below. 1. We calculate angular correlation functions of the z = 6–7 LAEs. The correlation lengths are estimated to be $$r_0=3.01^{+0.35}_{-0.35}$$ and $$2.66^{+0.49}_{-0.70}\,\,h^{-1}_{100}\:$$Mpc for the ≳ L* LAEs at z = 5.7 and 6.6, respectively. The average of the large-scale bias value is bavg = 4.13 ± 0.17 (4.54 ± 0.63) at z = 5.7 (z = 6.6) for the LAEs. Because Lyα photons emitted from LAEs in ionized bubbles can selectively escape from the partly neutral IGM at the EoR, observed LAEs at z = 6.6 should have clustering signals stronger than the intrinsic clustering. Based on this physical picture, we obtain the constraint of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ at z = 6.6 from the comparisons between our clustering measurements and two independent theoretical models. 2. We study the ≳ L* LAE clustering via HOD modeling. (Here, for the LAE clustering at the EoR of z = 6.6, we assume that the LAE clustering is not largely impacted by the cosmic reionization.) The best-fitting models indicate that the ≳ L* LAEs are hosted by the dark-matter halos with average masses of $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) =11.1^{+0.2}_{-0.4}$$ and $$10.8^{+0.3}_{-0.5}$$ at z = 5.7 and 6.6, respectively. Comparing the number densities of observed LAEs and those suggested from the HOD modeling, we find that dark-matter halos of only 1% (or less) down to the minimum-halo mass can host the ≳ L* LAEs. With the standard structure formation models, the bias evolution of the dark-matter halos indicates that the ≳ L* LAEs at z = 6–7 are progenitors of massive ∼6L* galaxies in the present-day universe. Acknowledgements We are grateful to useful discussion with Mark Dijikstra, Richard Ellis, Andrea Ferrara, Martin Haehnelt, Alex Hagen, Koki Kakiichi, Andrei Mesinger, Naveen Reddy, and Zheng Zheng. We acknowledge Jirong Mao and Anne Hutter for their comments. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. The NB816 filter was supported by Ehime University (PI: Y. Taniguchi). The NB921 filter was supported by KAKENHI (23244025) Grant-in-Aid for Scientific Research (A) through the Japan Society for the Promotion of Science (PI: M. Ouchi). This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory. This work is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (15H02064) Grant-in-Aid for Scientific Research (A) through Japan Society for the Promotion of Science. 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 at National Astronomical Observatory of Japan. Footnotes † Based on data obtained with the Subaru Telescope. The Subaru Telescope is operated by the National Astronomical Observatory of Japan. 1 There are some more OH windows redder than 921 nm, which include the 1010-nm window for the NB101 filter. However, these OH windows in the red band are contaminated by moderately strong OH lines, unlike the one at 921 nm. References Aihara H. et al.   2018, PASJ , 70, S8 Bacon R. et al.   2010, SPIE Proc. , 7735, 773508 Behroozi P. S., Wechsler R. H., Conroy C. 2013, ApJ , 770, 57 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 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   Cen R. 2003, ApJ , 591, 12 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., Barger A. J., Hu E. M. 2010, ApJ , 711, 928 CrossRef Search ADS   Cowie L. L., Hu E. M. 1998, AJ , 115, 1319 CrossRef Search ADS   Deharveng J.-M. et al.   2008, ApJ , 680, 1072 CrossRef Search ADS   Efstathiou G., Bernstein G., Tyson J. A., Katz N., Guhathakurta P. 1991, ApJ , 380, L47 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  Finkelstein S. L., Cohen S. H., Moustakas J., Malhotra S., Rhoads J. E., Papovich C. 2011, ApJ , 733, 117 CrossRef Search ADS   Fry J. N. 1996, ApJ , 461, L65 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L. 2006, MNRAS , 365, 1012 CrossRef Search ADS   Gallerani S., Ferrara A., Fan X., Choudhury T. R. 2008, MNRAS , 386, 359 CrossRef Search ADS   Gallerani S., Salvaterra R., Ferrara A., Choudhury T. R. 2008, MNRAS , 388, L84 CrossRef Search ADS   Gawiser E. et al.   2007, ApJ , 671, 278 CrossRef Search ADS   Greiner J. et al.   2009, ApJ , 693, 1610 CrossRef Search ADS   Groth E. J., Peebles P. J. E. 1977, ApJ , 217, 385 CrossRef Search ADS   Guaita L. et al.   2010, ApJ , 714, 255 CrossRef Search ADS   Guaita L., Francke H., Gawiser E., Bauer F. E., Hayes M., Östlin G., Padilla N. 2013, A&A , 551, A93 CrossRef Search ADS   Hagen A. et al.   2014, ApJ , 786, 59 CrossRef Search ADS   Harikane Y. et al.   2016, ApJ , 821, 123 CrossRef Search ADS   Harikane Y. et al.   2018, PASJ , 70, S11 Hashimoto T. et al.   2015, ApJ , 812, 157 CrossRef Search ADS   Hayashino T. et al.   2004, AJ , 128, 2073 CrossRef Search ADS   Hill G. J. et al.   2012, SPIE Proc. , 8446, 84460N 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., McMahon R. G. 1996, Nature , 382, 231 CrossRef Search ADS   Hutter A., Dayal P., Müller V. 2015, MNRAS , 450, 4025 CrossRef Search ADS   Iliev I. T., Shapiro P. R., McDonald P., Mellema G., Pen U.-L. 2008, MNRAS , 391, 63 CrossRef Search ADS   Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K. 2017, arXiv:1702.04867 Kakiichi K., Dijkstra M., Ciardi B., Graziani L. 2016, MNRAS , 463, 4019 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   Kobayashi M. A. R., Totani T., Nagashima M. 2010, ApJ , 708, 1119 CrossRef Search ADS   Kojima T., Ouchi M., Nakajima K., Shibuya T., Harikane Y., Ono Y. 2017, PASJ , 69, 44 CrossRef Search ADS   Konno A. et al.   2014, ApJ , 797, 16 CrossRef Search ADS   Konno A. et al.   2018, PASJ , 70, S16 Konno A., Ouchi M., Nakajima K., Duval F., Kusakabe H., Ono Y., Shimasaku K. 2016, ApJ , 823, 20 CrossRef Search ADS   Kovač K., Somerville R. S., Rhoads J. E., Malhotra S., Wang J. 2007, ApJ , 668, 15 CrossRef Search ADS   Landy S. D., Szalay A. S. 1993, ApJ , 412, 64 CrossRef Search ADS   Le Delliou M., Lacey C. G., Baugh C. M., Morris S. L. 2006, MNRAS , 365, 712 CrossRef Search ADS   Lee K.-S., Giavalisco M., Gnedin O. Y., Somerville R. S., Ferguson H. C., Dickinson M., Ouchi M. 2006, ApJ , 642, 63 CrossRef Search ADS   Lidz A., Zahn O., Furlanetto S. R., McQuinn M., Hernquist L., Zaldarriaga M. 2009, ApJ , 690, 252 CrossRef Search ADS   McLure R. J., Cirasuolo M., Dunlop J. S., Foucaud S., Almaini O. 2009, MNRAS , 395, 2196 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S. 2007, MNRAS , 381, 75 CrossRef Search ADS   McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Dutta S. 2008, MNRAS , 388, 1101 Madau P., Haardt F., Rees M. J. 1999, ApJ , 514, 648 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   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   Miyazaki S. et al.   2012, Proc. SPIE , 8446, 84460Z Momose R. et al.   2014, MNRAS , 442, 110 CrossRef Search ADS   Momose R. et al.   2016, MNRAS , 457, 2318 CrossRef Search ADS   Mortlock D. J. et al.   2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murayama T. et al.   2007, ApJS , 172, 523 CrossRef Search ADS   Nagamine K., Ouchi M., Springel V., Hernquist L. 2010, PASJ , 62, 1455 CrossRef Search ADS   Nakajima K. et al.   2012, ApJ , 745, 12 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   Ono Y. et al.   2012, ApJ , 744, 83 CrossRef Search ADS   Ono Y. et al.   2018, PASJ , 70, S10 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.   2005a, ApJ , 620, L1 CrossRef Search ADS   Ouchi M. et al.   2005b, ApJ , 635, L117 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   Partridge R. B., Peebles P. J. E. 1967, ApJ , 147, 868 CrossRef Search ADS   Pascarelle S. M., Windhorst R. A., Keel W. C., Odewahn S. C. 1996, Nature , 383, 45 CrossRef Search ADS   Peacock J. A., Dodds S. J. 1994, MNRAS , 267, 1020 CrossRef Search ADS   Peebles P. J. E. 1980, The Large-Scale Structure of the Universe  ( Princeton: Princeton University Press) Pentericci L. et al.   2011, ApJ , 743, 132 CrossRef Search ADS   Planck Collaboration 2016, A&A , 594, A13 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   Saito T., Shimasaku K., Okamura S., Ouchi M., Akiyama M., Yoshida M. 2006, ApJ , 648, 54 CrossRef Search ADS   Santos S., Sobral D., Matthee J. 2016, MNRAS , 463, 1678 CrossRef Search ADS   Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P. 2014, ApJ , 795, 20 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   Shimasaku K. et al.   2003, ApJ , 586, L111 CrossRef Search ADS   Simon P. 2007, A&A , 473, 711 CrossRef Search ADS   Sobacchi E., Mesinger A. 2015, MNRAS , 453, 1843 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   Stark D. P. et al.   2015, MNRAS , 454, 1393 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   Steidel C. C., Bogosavljević M., Shapley A. E., Kollmeier J. A., Reddy N. A., Erb D. K., Pettini M. 2011, ApJ , 736, 160 CrossRef Search ADS   Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E. 2008, ApJ , 688, 709 CrossRef Search ADS   Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S. 2010, ApJ , 724, 878 CrossRef Search ADS   Toshikawa J. et al.   2018, PASJ , 70, S12 Totani T. et al.   2014, PASJ , 66, 63 CrossRef Search ADS   Totani T., Aoki K., Hattori T., Kawai N. 2016, PASJ , 68, 15 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., Schmidt K. B., Trenti M., Bradley L. D., Stiavelli M. 2013, ApJ , 775, L29 CrossRef Search ADS   Verhamme A., Schaerer D., Atek H., Tapken C. 2008, A&A , 491, 89 CrossRef Search ADS   Zehavi I. et al.   2005, ApJ , 630, 1 CrossRef Search ADS   Zheng Z., Cen R., Weinberg D., Trac H., Miralda-Escudé J. 2011, ApJ , 739, 62 CrossRef Search ADS   Zheng Z., Coil A. L., Zehavi I. 2007, ApJ , 667, 760 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 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

# Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH). I. Program strategy and clustering properties of ∼2000 Lyα emitters at z = 6–7 over the 0.3–0.5 Gpc2 survey area

16 pages

/lp/ou_press/systematic-identification-of-laes-for-visible-exploration-and-k9oIr8a63e
Publisher
Oxford University Press
ISSN
0004-6264
eISSN
2053-051X
D.O.I.
10.1093/pasj/psx074
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present the SILVERRUSH program strategy and clustering properties investigated with ∼2000 Lyα emitters (LAEs) at z = 5.7 and 6.6 found in the early data of the Hyper Suprime-Cam (HSC) Subaru Strategic Program survey exploiting the carefully designed narrow-band filters. We derive angular correlation functions with the unprecedentedly large samples of LAEs at z = 6–7 over the large total area of 14–21 deg2 corresponding to 0.3–0.5 comoving Gpc2. We obtain the average large-scale bias values of bavg = 4.1 ± 0.2 (4.5 ± 0.6) at z = 5.7 (z = 6.6) for ≳ L* LAEs, indicating a weak evolution of LAE clustering from z = 5.7 to 6.6. We compare the LAE clustering results with two independent theoretical models that suggest an increase of an LAE clustering signal by the patchy ionized bubbles at the epoch of reionization (EoR), and estimate the neutral hydrogen fraction to be $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ at z = 6.6. Based on the halo occupation distribution models, we find that the ≳ L* LAEs are hosted by dark-matter halos with an average mass of $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) =11.1^{+0.2}_{-0.4}$$ ($$10.8^{+0.3}_{-0.5}$$) at z = 5.7 (6.6) with a Lyα duty cycle of 1% or less, where the results of z = 6.6 LAEs may be slightly biased, due to the increase of the clustering signal at the EoR. Our clustering analysis reveals the low-mass nature of ≳ L* LAEs at z = 6–7, and that these LAEs probably evolve into massive super-L* galaxies in the present-day universe. 1 Introduction Lyα emitters (LAEs) are star-forming galaxies (and active galactic nuclei, AGNs) with a strong Lyα emission line. In the 1960s, it was theoretically predicted that such galaxies are candidates of very young galaxies residing at z ∼ 10–30 (Partridge & Peebles 1967). About 30 years after these predictions, deep observations identified a few LAEs around AGN regions at z = 4.6 (Hu & McMahon 1996) and z = 2.4 (Pascarelle et al. 1996). Subsequently, LAEs in blank fields have routinely been found by deep and wide-field narrow-band observations conducted under the Hawaii Survey (Cowie & Hu 1998), the Large Area Lyman Alpha Survey (Rhoads et al. 2000), and Subaru Surveys (e.g., Ouchi et al. 2003). These blank-field LAE surveys reveal various statistical properties of LAEs up to z ∼ 7 including the evolution of Lyα luminosity functions (LFs; Malhotra & Rhoads 2004; Kashikawa et al. 2006, 2011; Deharveng et al. 2008; Ouchi et al. 2008, 2010; Cowie et al. 2010; Hu et al. 2010; Konno et al. 2014; Matthee et al. 2015; Santos et al. 2016; Ota et al. 2017; Zheng et al. 2017). Moreover, multi-wavelength follow-up imaging and spectroscopy of LAEs reveal that typical (∼L*) LAEs (LLyα = 1042–1043 erg s−1) at z ∼ 2–3 have a stellar mass of 108–109 M⊙ (Gawiser et al. 2007; Hagen et al. 2014), a star-formation rate (SFR) of ∼10 M⊙ (Nakajima et al. 2012), and a gas-phase metallicity of ∼0.1 Z⊙ (Finkelstein et al. 2011; Guaita et al. 2013; Kojima et al. 2017). LAEs have strong high-ionization lines such as [O iii]5007, indicative of the high ionization state of the inter-stellar medium (ISM) given by intense ionizing radiation from young massive stars (Nakajima & Ouchi 2014). Outside the highly-ionized ISM, LAEs show diffuse extended Lyα halos in the circum-galactic medium (CGM) up to a few 10 kpc (Hayashino et al. 2004; Steidel et al. 2011). The bright Lyα halos extending over a few hundred kpc are known as Lyα blobs and their total Lyα luminosities are ∼1044 erg s−1 (Steidel et al. 2000; Matsuda et al. 2004; cf. Saito et al. 2006), more than an order of magnitude brighter than the diffuse Lyα halos (Momose et al. 2014, 2016). Deep spectroscopic observations identified large spatial overdensities of LAEs up to z ∼ 6 (Shimasaku et al. 2003; Ouchi et al. 2005a) that are possibly progenitors of massive galaxy clusters found today. Recent deep spectroscopy also detected strong Lyα emission of ∼10 LAEs at z = 7.0–8.7, pushing the redshift frontier of LAEs to-date (e.g., Pentericci et al. 2011; Ono et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Schenker et al. 2014; Oesch et al. 2015; Zitrin et al. 2015). There are three major scientific goals in recent LAE studies. The first goal is characterizing the nature of low-mass young galaxies at high-z. Because LAEs are low stellar mass galaxies with strong Lyα indicative of a starburst with young massive stars producing a large amount of ionizing photons (Ouchi et al. 2013; Sobral et al. 2015; Stark et al. 2015), LAEs are used as probes of young galaxies. This motivation is the same as the one behind the original predictions of LAEs given by Partridge and Peebles (1967). The second goal is to understand the density and dynamics of H i clouds in the ISM and the CGM of star-forming galaxies. Due to the resonance nature of Lyα emission, Lyα photons are scattered by H i gas in the ISM and the CGM, and the dynamical properties of Lyα sources are not easily understood with the observed Lyα velocity fields. Instead, both the density and kinematics of H i gas in the ISM and the CGM are encoded in the observed Lyα line velocity and spatial profiles. The combination of the Lyα line profile observations and models provides important information about the H i gas of the ISM and the CGM (Verhamme et al. 2008; Zheng et al. 2011; Shibuya et al. 2014; Hashimoto et al. 2015; Momose et al. 2016). The third goal is to reveal the cosmic reionization history and its properties. Due to the strong damping absorption of Lyα by the IGM H i gas at the epoch of reionization (EoR: z > 6), Lyα emission of LAEs is absorbed in the partially neutral IGM. The Lyα absorption depress the observed Lyα luminosities of LAEs, which causes a decrease in the Lyα LF towards the early stage of the EoR (Malhotra & Rhoads 2004; Kashikawa et al. 2006, 2011; Ouchi et al. 2010; Hu et al. 2010; Santos et al. 2016; Ota et al. 2017). Moreover, the clustering signal of observed LAEs is boosted by the existence of the ionized bubbles (Furlanetto et al. 2006; McQuinn et al. 2007; Ouchi et al. 2010), because the Lyα absorption is selectively weak for LAEs residing in the ionized bubbles of the IGM at the EoR. The Lyα LFs and clustering of LAEs are important quantities in characterizing the cosmic reionization history and ionized bubble topologies at the EoR. There are three new large instruments that are used for LAE observations; the Very Large Telescope/Multi-Unit Spectroscopic Explorer (MUSE: Bacon et al. 2010), the Hobby–Eberly Telescope Dark Energy Experiment/Visible Integral-Field Replicable Unit Spectrograph (HETDEX/VIRUS: Hill et al. 2012), and the Subaru/Hyper Suprime-Cam (HSC: Miyazaki et al. 2012). These three instruments can cover the complementary parameter space of LAEs in redshift, depth, and survey volume. Specifically, Subaru/HSC has the capability to take large-area (>10 deg2) deep images with custom-made narrow-bands (section 2) targeting LAEs at z ∼ 2–7. The combination of Subaru/HSC imaging and deep follow-up spectroscopy allows us to address the key issues of LAEs for accomplishing the three major goals described above. One of the most unique studies realized with Subaru/HSC is clustering of LAEs at the EoR (i.e., z ≳ 6), which requires high-sensitivity observations over a large area of the sky, as we show in this study. LAE clustering signals depend on the hosting dark-matter halo properties of LAEs, and the distribution of the patchy ionized IGM (i.e., ionized bubbles) that allows Lyα photons to escape from the partially neutral IGM at the EoR (Furlanetto et al. 2006; McQuinn et al. 2007; Ouchi et al. 2010). Theoretical models suggest that a Lyα LF measurement, a popular statistical quantity of LAEs, has a degeneracy between the Lyα escape fraction $$f_{\rm esc}^{\rm Ly \alpha }$$ (a flux ratio of observed to intrinsic Lyα emission: Le Delliou et al. 2006; Mao et al. 2007; Kobayashi et al. 2010) and duty cycle $$f_{\rm duty}^{\rm LAE}$$ (an abundance ratio of Lyα emitting galaxies to dark-matter haloes: Nagamine et al. 2010). Moreover, at the EoR, there is another degeneracy between the ionized bubble topology and neutral hydrogen fraction xH i of the IGM (Kakiichi et al. 2016). Theoretical studies show that these degeneracies can be resolved with a combination of Lyα LF and clustering measurements (e.g., Nagamine et al. 2010; Sobacchi & Mesinger 2015; Hutter et al. 2015; Kakiichi et al. 2016). The LAE samples at the post-EoR, e.g., z = 5.7, (the EoR, e.g., z = 6.6) are useful to determine $$f_{\rm esc}^{\rm Ly \alpha }$$ and $$f_{\rm duty}^{\rm LAE}$$. In this paper, we describe our narrow-band filters and the program strategy of the HSC LAE studies named “Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC” (SILVERRUSH: section 2), and present our LAE samples (section 3). We obtain measurements of LAE clustering (section 4) that are calculated with the early samples of HSC LAEs at z = 5.7 and 6.6 (Shibuya et al. 2018a) based on the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP: HSC et al. in preparation) data whose first data release is presented in Aihara et al. (2018). Based on the LAE clustering measurements, we discuss the cosmic reionization history and dark-matter halo properties of LAEs by a comparison with the ΛCDM structure formation models in section 5. Throughout this paper, magnitudes are in the AB system. We adopt a cosmology parameter set of (h, Ωm, ΩΛ, ns, σ8) = (0.7, 0.3, 0.7, 1.0, 0.8), consistent with the latest Planck results (Planck Collaboration 2016). 2 Narrow-band filters and the program strategy for the LAE studies In our LAE studies, we use four custom-made narrow-band filters that are fabricated for the HSC-SSP survey. We design the HSC narrow-band filters targeting the OH sky windows at 816, 921, and 1010 nm; these correspond to NB816 (PI: Y. Taniguchi), NB921 (PI: M. Ouchi), and NB101 (PI: K. Shimasaku) filters, which identify the redshifted Lyα of LAEs at z = 5.7, 6.6, and 7.3, respectively. These three narrow-band filters densely cover LAEs from the heart of the EoR to the post-EoR. Another HSC narrow-band filter is NB387 (PI: K. Shimasaku), the central wavelength of which is 387 nm. The NB387 band is the bluest narrow-band among the existing HSC filters, targeting the lowest-redshift LAEs identified by the HSC observations. The central wavelength of 387 nm is specifically chosen, because this band identifies LAEs at z = 2.2 that have major strong lines of [O ii]3727 and Hα6563 which fall in the OH sky windows of 1.19 and 2.10 μm, respectively. The other major strong lines, Hβ4861 and [O iii]5007, of z = 2.2 LAEs are placed at the H band, where the effects of the OH sky emission are unavoidable, due to no clear OH windows. We carefully determine the full width at half maximum (FWHM) value for each narrow-band filter transmission curve. First, we determine the FWHM of NB921 that properly covers a clear OH window at the reddest band, allowing the typical transmission curve errors of the filter fabrication.1 We aim to accomplish the same detection limit of the rest-frame Lyα equivalent width EW0 in our samples of LAEs at z = 2.2, 5.7, and 6.6. An observed-frame Lyα equivalent width, EWobs, has the relation of EWobs = (1 + z)EW0. Because a narrow-band flux excess is proportional (inversely proportional) to $$\mathit {EW}_{\rm obs}$$ (a narrow-band FWHM value), we scale the FWHM value of the NB921 band (FWHMNB921) by (1 + z) to design the FWHM values of NB387 (FWHMNB387) and NB816 (FWHMNB816). More specifically, the determined FWHM values keep the relations of FWHMNB387 = FWHMNB921(1 + 2.18)/(1 + 6.58) and FWHMNB816 = FWHMNB921(1 + 5.71)/(1 + 6.58), where the numbers of 2.18, 5.71, and 6.58 correspond to the redshifts of the Lyα lines (121.567 nm) that fall in the centers of the NB387, NB816, and NB921 bandpasses, respectively. Here, the FWHM of the NB101 band does not follow this scaling relation. Instead, we choose an FWHM that is the narrowest limit that can be reasonably accomplished by the present technology of the HSC narrow-band filter fabrication. This is because Lyα lines of the NB101 LAEs fall in a wavelength range where the HSC imaging cannot reach a very deep magnitude limit, due to the relatively low CCD quantum efficiency and the unclean OH window. In this way, we obtain the designed specifications of (central wavelength, FWHM wavelength) = (387.0 nm, 5.5 nm) for NB387, (816.0 nm, 11.6 nm) for NB816, (921.0 nm, 13.1 nm) for NB921, and (1009.5 nm, 9.0 nm) for NB101. The HSC narrow-band filters are made of a B270 or quartz glass with multi-layer interferometric coatings that make the sharp cut-on and cut-off bandpasses for the narrow-band filters. The diameter of the HSC narrow-band filter is large; 600 mm. Although uniform coatings are key to accomplishing the spatially homogeneous transmission of the narrow-band filters, small non-uniformities are included in the fabrication processes which give a filter characterization that is slightly different from the design. Nevertheless, we have tried to make the coating distribution as uniform as possible to achieve spatially homogeneous transmission curves for the narrow-band filters. The picture of NB921 is shown in figure 1. Fig. 1. View largeDownload slide HSC NB921 filter in a filter holder. (Color online) Fig. 1. View largeDownload slide HSC NB921 filter in a filter holder. (Color online) Figures 2 and 3 present the transmission curves of the HSC NB816 and NB921 filters that are used in the HSC first-data release (Aihara et al. 2018). We measure the transmission curves at 21 points over the 600-mm-diameter circular filter. The detectors of HSC cover the circular area from the center r = 0 to the radius of r = 280 mm. On the filter, we define the axis of the 0° angle for the arbitrary direction from the filter center to the edge. At an angle of 0°, we measure transmission curves at five positions; r = 50, 100, 150, 200, and 270 mm. We change the angle to 90°, 180°, and 270° in a counter-clockwise direction, and thereby obtain transmission curves at 20 positions (= 5 × 4). Including the measurement at the center (r = 0 mm), we have transmission curves at a total of 21 positions. Note that the measurement position of NB921 is slightly different from that of NB816. The r = 270 mm position is replaced with r = 265 mm for the NB921 filter. We find that the area-weighted average central wavelength and FWHM of NB816 (NB921) are 817.7 nm and 11.3 nm (921.5 nm and 13.5 nm), respectively, and that the central-position central wavelength and FWHM of NB816 (NB921) are 816.8 nm and 11.0 nm (920.5 nm and 13.3 nm), respectively. At any positions within r ≤ 265 mm, the peak transmissions of NB816 and NB921 are high; >90%. The deviations of the central and FWHM wavelengths are typically within ≃ 0.3% and ≃ 10%, respectively. Fig. 2. View largeDownload slide Filter transmission curves of NB816. The gray lines represent the transmission curves at the 21 positions (see text). The black line is the same as the gray lines, but for the transmission curve at the central position. The red line represents the area-averaged mean transmission curve that is shown in the Subaru website ⟨http://www.naoj.org/Observing/Instruments/HSC/sensitivity.html⟩. (Color online) Fig. 2. View largeDownload slide Filter transmission curves of NB816. The gray lines represent the transmission curves at the 21 positions (see text). The black line is the same as the gray lines, but for the transmission curve at the central position. The red line represents the area-averaged mean transmission curve that is shown in the Subaru website ⟨http://www.naoj.org/Observing/Instruments/HSC/sensitivity.html⟩. (Color online) Fig. 3. View largeDownload slide Same as figure 2, but for NB921. (Color online) Fig. 3. View largeDownload slide Same as figure 2, but for NB921. (Color online) Using the HSC narrow-band filters, we study LAEs at z = 2.2, 5.7, and 6.6 (7.3) over the large areas of 26 (3.5) deg2, which are about an order of magnitude larger than those of the previous studies for a given redshift. Exploiting samples of z = 2.2 and 5.7 LAEs, we study the evolution of LAEs at z = 2.2–5.7, from the lowest redshift accessible by ground-based observations (z ∼ 2) to the high-z edge of the post-EoR (z ∼ 6). The LAEs at the three redshifts of z = 5.7, 6.6, and 7.3 allow us to study the evolution of galaxy formation as well as cosmic reionization with the Lyα damping wing absorption of the neutral IGM. The LAEs at the post-EoR of z = 5.7 play the role of a baseline for the properties of LAEs with no cosmic reionization effects given by the IGM Lyα damping wing absorption. The comparisons of z = 5.7 and 6.6 LAEs over the 26 deg2 sky provide large statistical results, while the z = 7.3 LAEs in the moderately large 3.5 deg2 area allows us to investigate galaxy formation and cosmic reionization at the heart of the EoR. Based on the data sets of the HSC-SSP survey (see Aihara et al. 2018 for the details of the first-data release) and our extensive spectroscopic follow-up observations with Subaru, Keck, and Magellan telescopes, we start the program named SILVERRUSH (section 1). This program is one of a set of twin programs; the other is is a study for dropouts, called the Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH), and is detailed in Ono et al. (2018), Harikane et al. (2018), and Toshikawa et al. (2018). In a series of SILVERRUSH papers, we present the SILVERRUSH program strategy (this paper), LAE sample selections (Shibuya et al. 2018a) and spectroscopic follow-up observations (Shibuya et al. 2018b), Lyα LF evolution (Konno et al. 2018), LAE clustering evolution (this paper), and comparisons with numerical simulation results (A. K. Inoue et al. in preparation). These are early papers of the SILVERRUSH program, the observations of which are still underway. In these early papers, we present the results of the LAE studies based on the images taken until 2016 April that include neither NB387 (for z = 2.2 LAEs) nor NB101 (for z = 7.3 LAEs) data, but only NB816 and NB921 data for z = 5.7 and 6.6 LAEs, respectively (section 3). Here, our aim is that SILVERRUSH results would serve as the baselines of LAE properties useful for the ongoing LAE studies including MUSE, HETDEX, and the forthcoming z > 7 galaxy studies of the James Webb Space Telescope (JWST) and extremely large telescopes (ELTs). 3 Observations and sample The HSC-SSP survey started observations in 2014 March, taking deep broad- and narrow-band images on large areas of the sky (PI: S. Miyazaki). In our study, we use the HSC-SSP survey (HSC et al. in preparation) data taken until 2016 April with two narrow-bands (NB816 and NB921) and five broad-bands (grizy) useful for our LAE studies. Because the data of NB387 and NB101 (for detections of LAEs at z = 2.2 and 7.3) were not taken until 2016 April, we present the results of imaging with the two narrow-band NB816 and NB921 (for LAEs at z = 5.7 and 6.6), the observations of which have been partly conducted. The data are reduced with the HSC pipeline software, hscPipe version 4.0.2 (Bosch et al. 2018). The 5σ depths of the imaging data are typically ≃ 25–25.5 and ≃ 26–27 mag in narrow- and broad-bands, respectively (see Shibuya et al. 2018a for more details). Note that the HSC-SSP survey data set used in this study is notably larger than that of the first-data release (Aihara et al. 2018), which is composed of images taken only in 2014 March and 2015 November. Our LAEs are selected with the combinations of broad- and narrow-band colors down to the 5σ detection limits. The LAEs should have a narrow-band excess, the existence of a Gunn–Peterson trough, and no detection of blue continuum fluxes. The color selection criteria are similar to those of Ouchi et al. (2008, 2010), which are defined as   \begin{eqnarray} i-{\textit NB}{\textit 8}{\textit 1}{\textit 6}\ge 1.2\ {\rm and}\ g>g_{3\sigma }\ {\rm and}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \nonumber \\ \ \ \ \ \ \left[(r\le r_{3\sigma }\ {\rm and}\ r-i\ge 1.0)\ {\rm or}\ (r>r_{3\sigma })\right] \end{eqnarray} (1)and   \begin{eqnarray} z-{\textit NB}{\textit 9}{\textit 2}{\textit 1}\ge 1.0\ {\rm and}\ g>g_{3\sigma }\ {\rm and}\ r>r_{3\sigma }\ {\rm and}\ \nonumber \\ \ \ \ \ \ \left[(z\le z_{3\sigma }\ {\rm and}\ i-z\ge 1.3)\ {\rm or}\ (z>z_{3\sigma })\right] \end{eqnarray} (2)for z = 5.7 and 6.6 LAEs, respectively, where g3σ (r3σ) is the 3σ limiting magnitude of the g (r) band that ensures no detection of a continuum bluer than the Lyman break. Similarly, z3σ is the 3σ detection limit of the z band. After applying the candidate screening on the basis of hscPipe parameters+flags and visual inspection, we obtain a total of 2354 LAEs (1081 and 1273 LAEs at z = 5.7 and 6.6, respectively) from the LAE ALL catalogs (Shibuya et al. 2018a). We investigate the 2354 LAEs, and make homogeneous samples over the survey areas that consist of 1832 LAEs (959 and 873 LAEs at z = 5.7 and 6.6, respectively) with a common narrow-band limiting magnitude of NB816 < 25.0 or NB921 < 25.0 that is surely brighter than the 5σ detection levels. The homogenous samples of z = 5.7 and 6.6 LAEs have very similar threshold luminosities (the rest-frame equivalent width) of $$L_{\rm Ly\alpha }^{\rm th}>6.3 \times 10^{42}$$ and 7.9 × 1042 erg s−1 ($$\mathit {EW}_0 \gtrsim 20\,\,$$Å) at z = 5.7 and 6.6, respectively. The threshold luminosities correspond to ≃ L* luminosities (e.g., Ouchi et al. 2008, 2010; Konno et al. 2018). The LAE sample selection is detailed in Shibuya et al. (2018a). These LAEs are found in an area totalling 13.8 (21.2) deg2 consisting of the fields named D-ELAIS-N1, D-DEEP2-3, UD-COSMOS, and UD-SXDS for the z = 5.7 LAE sample (D-ELAIS-N1, D-DEEP2-3, D-COSMOS, UD-COSMOS, and UD-SXDS for the z = 6.6 LAE sample). Because, for the z = 5.7 and 6.6 LAE samples, the redshift ranges are z = 5.726 ± 0.046 and z = 6.580 ± 0.056 in the case that top-hat selection functions of LAEs with the FWHMs of the narrow-band filters are assumed, the total survey volumes covered with the narrow-band transmission are 1.2 × 107 and 1.9 × 107 comoving Mpc3 at z = 5.7 and 6.6, respectively. Note that the total areas of 13.8 and 21.2 deg2 for the z = 5.7 and 6.6 samples correspond to 0.3 and 0.5 comoving Gpc2 areas, respectively. These large survey areas allow us to study average properties of LAEs at z = 6–7, observabilities of which are influenced by patchy (10–100 Mpc) ionized bubbles probably existing at the end of the EoR (Furlanetto et al. 2006). In our LAE samples, a total of 97 LAEs at z = 5.7 and 6.6 are confirmed by our spectroscopic observations and cross-matching the existing spectra (Shibuya et al. 2018b). Because there are 81 LAE candidates that have spectroscopic identifications that are obtained via our past and present programs, we estimate the contamination rates of our LAE samples using these 81 LAE candidates, which consist of 53 and 28 LAEs at z = 5.7 and 6.6, respectively. We find that four out of 53 (four out of 28) candidates are foreground contamination objects, and estimate the contamination rates to be ∼8% and ∼14% for the samples of LAEs at z = 5.7 and 6.6, respectively. We also investigate contamination rates of bright LAEs that are brighter than 24 magn in the narrow-band. There are six and 12 bright LAE candidates with spectroscopic identifications. The spectroscopic results indicate that one out of six (four out of 12) LAE candidates are foreground interlopers, which corresponds to a contamination rate of ∼17% (∼33%) for the sample of LAEs at z = 5.7 (6.6). Although the contamination rates depend on magnitude, the contamination rates indicated by the spectroscopic confirmation range around 0%–30% in our z = 5.7 and 6.6 LAE samples. Although the contamination rates include large uncertainties due to the small number of our spectroscopically confirmed sources at this early stage of our LAE studies, we assume that the contamination rates are 0%–30% in the LAE samples that are used in our analysis below. 4 Clustering analysis and results Figures 4 and 5 present the sky distribution of the LAEs at z = 5.7 and 6.6, respectively. In figure 4, we find that the density of z = 5.7 LAEs in the D-ELAIS-N1 field is largely biased towards the southern region, indicating the possibilities of the detection incompleteness and the real large overdensity. We have carefully investigated our D-ELAIS-N1 images, but found no significant spatial fluctuations in the detection completeness produced by the inhomogeneity of the data handling and the data qualities, including depths and seeing sizes. However, there remains a possibility that unknown effects of the data reduction may produce a large overdensity. Because the origin of this large overdensity is still unknown, we do not use the data of D-ELAIS-N1 in this study. Instead, we discuss this overdensity in R. Higuchi et al. (in preparation). We thus use the remaining 734 LAEs at z = 5.7 in our analysis. Fig. 4. View largeDownload slide Sky distribution of the LAEs at z = 5.7. The red squares, magenta diamonds, and black circles represent positions of narrow-band bright (<23.5 mag), medium-bright (23.5–24.0 mag), and faint (24.0–25.0 mag) LAEs, respectively. The large red open square indicates the LAEs with spatially extended Lyα emission (Shibuya et al. 2018b). The gray shades denote either the areas with no HSC data or the masked regions with a bad data quality. The scale on the map is marked in angles (degrees) and the projected distances (comoving megaparsecs). (Color online) Fig. 4. View largeDownload slide Sky distribution of the LAEs at z = 5.7. The red squares, magenta diamonds, and black circles represent positions of narrow-band bright (<23.5 mag), medium-bright (23.5–24.0 mag), and faint (24.0–25.0 mag) LAEs, respectively. The large red open square indicates the LAEs with spatially extended Lyα emission (Shibuya et al. 2018b). The gray shades denote either the areas with no HSC data or the masked regions with a bad data quality. The scale on the map is marked in angles (degrees) and the projected distances (comoving megaparsecs). (Color online) Fig. 5. View largeDownload slide Same as figure 4, but for the LAEs z = 6.6. The large red open squares indicate the LAEs with spatially extended Lyα emission including Himiko (Ouchi et al. 2009) and CR7 (Sobral et al. 2015). See Shibuya et al. (2018b) for more details. (Color online) Fig. 5. View largeDownload slide Same as figure 4, but for the LAEs z = 6.6. The large red open squares indicate the LAEs with spatially extended Lyα emission including Himiko (Ouchi et al. 2009) and CR7 (Sobral et al. 2015). See Shibuya et al. (2018b) for more details. (Color online) We quantify clustering properties based on the measurements of the angular correlation functions in the following sections. 4.1 Angular correlation function We derive angular two-point correlation functions (ACFs) of our LAEs in the same manner as Ouchi et al. (2003, 2005b, 2010) and Harikane et al. (2016). We use the estimator of Landy and Szalay (1993),   $$\omega _{\rm obs}(\theta ) = \frac{DD(\theta )-2DR(\theta )+RR(\theta )}{RR(\theta )},$$ (3)where DD(θ), DR(θ), and RR(θ) are the numbers of galaxy–galaxy, galaxy–random, and random–random pairs normalized by the total number of pairs in each of the samples of the pairs. We use a random catalog (Coupon et al. 2018) that has a surface number density of 100 arcmin−2. The random catalog has the same geometrical constraint as the distribution of our LAEs, representing our survey areas. Statistical errors are estimated using Jackknife resampling with subsamples that each have an area of ∼10002 arcsec2. Figure 6 shows the observed ACFs ωobs(θ) of our LAEs at z = 5.7 and 6.6. The ACF measurements cover the scale of ∼1–100 comoving Mpc that is indicated in the upper axes of figure 6. We find that these ACFs are consistent with those obtained with the previous Subaru/Suprime-Cam data (Ouchi et al. 2010), and that our present HSC data provide the large-scale ACFs with uncertainties smaller than those of the previous data. Fig. 6. View largeDownload slide Angular correlation function (ACF) and bias of the LAEs at z = 5.7 (left) and 6.6 (right). The top and bottom panels show the ACFs and the bias. The black-filled and gray-open squares represent the ACFs of our LAEs with and without the IC correction, respectively, while the filled circles denote the ACFs with the IC correction derived by Ouchi et al. (2010). For presentation purposes, we slightly shift the gray open squares along the abscissa. In the top panels, the solid and dotted lines present the best-fitting power-law functions of our ACFs and the ACFs of the underlying dark matter predicted by the linear theory (e.g., Peacock & Dodds 1994). Because the power-law spatial correlation function is projected on sky with the method in Simon (2007), the best-fitting power-law functions are curves. In the bottom panels, the solid and dashed lines indicate the average bias and the 1σ error values, respectively. The crosses in the top left-hand panel show the ACFs obtained by Murayama et al. (2007). The top axis denotes the projected distance in comoving megaparsecs. Fig. 6. View largeDownload slide Angular correlation function (ACF) and bias of the LAEs at z = 5.7 (left) and 6.6 (right). The top and bottom panels show the ACFs and the bias. The black-filled and gray-open squares represent the ACFs of our LAEs with and without the IC correction, respectively, while the filled circles denote the ACFs with the IC correction derived by Ouchi et al. (2010). For presentation purposes, we slightly shift the gray open squares along the abscissa. In the top panels, the solid and dotted lines present the best-fitting power-law functions of our ACFs and the ACFs of the underlying dark matter predicted by the linear theory (e.g., Peacock & Dodds 1994). Because the power-law spatial correlation function is projected on sky with the method in Simon (2007), the best-fitting power-law functions are curves. In the bottom panels, the solid and dashed lines indicate the average bias and the 1σ error values, respectively. The crosses in the top left-hand panel show the ACFs obtained by Murayama et al. (2007). The top axis denotes the projected distance in comoving megaparsecs. Although our survey area is large, we evaluate the integral constraint (IC) (Groth & Peebles 1977). The IC value corresponds to the observational offset in ωobs(θ) originated by a limited survey area. Including the correction for the number of objects in the sample, N, the true ACF is given by   $$\omega (\theta )=\omega _\mathrm{obs}(\theta )+\mathrm{IC}+\frac{1}{N}.$$ (4)We evaluate the integral constraint with   $$\mathrm{IC}=\frac{\Sigma _i RR(\theta _i)\omega (\theta _i)}{\Sigma _i RR(\theta _i)}.$$ (5)We adopt the model ACF of Simon (2007) for the function of ω(θ) that is detailed in subsection 4.2. We show IC values in table 1, and plot the ACFs corrected for IC in figure 6. Note that the IC values are very small, because our survey areas are very large. Table 1. LAE samples and the clustering measurements. z  $$L_{\rm Ly\alpha }^{\rm th}$$  ng  N  IC  r0  bavg  $$r_0^{\rm max}$$*  $$b_{\rm avg}^{\rm max}$$  χ2/dof    (1042 erg s−1)  (10−5 Mpc−3)    (10−3)  ($$h_{100}^{-1}\:$$Mpc)    ($$h_{100}^{-1}\:$$Mpc)      (1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10)  5.7  6.3  7.3  959†  1.83  $$3.01^{+0.35}_{-0.35}$$  4.13 ± 0.17  $$4.47^{+0.52}_{-0.52}$$  5.90 ± 0.24  5.69/5  6.6  7.9  2.4  873  1.35  $$2.66^{+0.49}_{-0.70}$$  4.54 ± 0.63  $$3.95^{+0.73}_{-1.04}$$  6.49 ± 0.90  2.70/2  z  $$L_{\rm Ly\alpha }^{\rm th}$$  ng  N  IC  r0  bavg  $$r_0^{\rm max}$$*  $$b_{\rm avg}^{\rm max}$$  χ2/dof    (1042 erg s−1)  (10−5 Mpc−3)    (10−3)  ($$h_{100}^{-1}\:$$Mpc)    ($$h_{100}^{-1}\:$$Mpc)      (1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10)  5.7  6.3  7.3  959†  1.83  $$3.01^{+0.35}_{-0.35}$$  4.13 ± 0.17  $$4.47^{+0.52}_{-0.52}$$  5.90 ± 0.24  5.69/5  6.6  7.9  2.4  873  1.35  $$2.66^{+0.49}_{-0.70}$$  4.54 ± 0.63  $$3.95^{+0.73}_{-1.04}$$  6.49 ± 0.90  2.70/2  * For comparison, r0 is expressed with h100. † This number includes sources in D-ELAIS-N1. View Large A clustering signal of galaxies is diluted by contamination objects in a galaxy sample. If the galaxy sample includes randomly distributed contamination objects with a fraction fc, the value of the ACF is reduced by a factor of (1 − fc)2. The ACF corrected for the randomly-distributed contamination ωmax is written as   $$\omega _{\rm max} = \frac{\omega }{(1-f_{\rm c})^2}.$$ (6)This is the maximum reduction of the ACF, because the real contamination objects are not randomly distributed but spatially correlated. For reference, we evaluate the possible maximum values of ωmax, using a contamination fraction with an upper limit of fc = 0.3 (section 3). 4.2 Correlation length and bias We fit the ACFs with a simple power-law model of the spatial correlation function,   $$\xi (r)=\left(\frac{r}{r_0}\right)^{-\gamma },$$ (7)where γ, r0, and r are the slope of the power law, the correlation length, and the spatial separation between two galaxies, respectively. The spatial correlation function is related to the ACF by the Limber equation (Peebles 1980; Efstathiou et al. 1991), which is an integral equation of the (three-dimensional) spatial correlation function connecting with the (two-dimensional) ACF. However, Simon (2007) claims that the Limber equation does not provide accurate values in a very large separation of galaxies that have redshift-distribution distances narrower than the transverse distance in the case of narrow-band-selected LAEs. We thus adopt the method that Simon (2007) suggests, and derive r0 of our LAEs, fitting the power-law functions to our ACFs. Because no meaningful constraints on γ are obtained with the present samples, we adopt the fiducial γ value of γ = 1.8 that is adopted in the previous clustering analyses for LAEs (Ouchi et al. 2003; Gawiser et al. 2007; Kovač et al. 2007; Ouchi et al. 2010). To investigate the dependencies of our results on the value of γ, we also use the other fixed γ values of γ = 1.6 and 2.0 that bracket the possible range of the power-law index of high-z galaxies at z ∼ 4–6 (Lee et al. 2006; McLure et al. 2009). We have found that neither r0 nor bias values change by over 10%, and that r0 and bias values fall well within the errors. Here, we use the redshift distribution of our LAEs estimated with the narrow-band transmission curve. The best-fitting power-law functions are presented in figure 6. The correlation lengths thus obtained are $$r_0=3.01^{+0.35}_{-0.35}$$ and $$2.66^{+0.49}_{-0.70}$$$$h^{-1}_{100}\:$$Mpc for LAEs at z = 5.7 and 6.6, respectively, as summarized in table 1. To make comparisons with the previous results, we express the correlation lengths using h100, the Hubble constant in units of 100 km s−1, instead of 70 km s−1. Table 1 also shows the upper limit of the correlation length, $$r_0^{\rm max}$$, which is given by   $$r_0^{\rm max}=r_0\left(\frac{1}{1-f_\mathrm{C}}\right)^{2/\gamma }.$$ (8)We obtain $$r_0^{\rm max} = 4.47^{+0.52}_{-0.52}$$ and $$3.95^{+0.73}_{-1.04}\,\,h_{100}^{-1}\:$$Mpc for our z = 5.7 and 6.6 LAEs, respectively. These correlation length values are consistent with the previous measurements within the errors. Ouchi et al. (2010) obtain (r0, $$r_0^{\rm max}$$) of ($$3.12^{+0.33}_{-0.36}$$, $$4.30^{+0.45}_{-0.50}$$) and ($$2.31^{+0.65}_{-0.85}$$, $$3.60^{+1.02}_{-1.32}$$) for LAEs at z = 5.7 and 6.6, respectively. In the framework of the ΛCDM model, the galaxy dark-matter bias is defined by   $$b_{\rm g}(\theta )^2 \equiv \omega (\theta )/\omega _{\rm dm}(\theta ),$$ (9)where ωdm(θ) is the ACF of dark matter predicted by the linear theory model (e.g., Peacock & Dodds 1994). The model calculations are made with the same survey volumes as those of our LAE samples. Note that bg is the bias value equivalent to the one given by the three-dimensional spatial correlation functions, $$b_{\rm g}^2=\xi (r)/\xi _{\rm dm}(r)$$, where ξdm(r) is the spatial correlation function of dark matter. In figure 6, the top and bottom panels show ωdm(θ) and bg(θ), respectively. We estimate the average bias bavg by averaging bg(θ) with error weighting over the angular range presented in figure 6. We also calculate the upper limits of bias values, $$b_{\rm avg}^{\rm max}$$, with the maximal contamination correction that corresponds to $$A_{\omega }^{\rm max}$$ [equation (6); see subsection 4.1]. We obtain bavg = 4.13 ± 0.17 and 4.54 ± 0.63 ($$b_{\rm avg}^{\rm max} = 5.90\,\,\pm \,\, 0.24$$ and 6.49 ± 0.90) for our z = 5.7 and 6.6 LAEs, respectively. We summarize bavg and $$b_{\rm avg}^{\rm max}$$ values thus obtained in table 1. We compare these bias values with those obtained in the previous studies. We find that the bias value of our LAEs are consistent with those derived by Ouchi et al. (2005a, 2010) within the ∼1σ errors. Figure 7 presents bias of LAEs as a function of redshift. In figure 7, we plot bias measurements of LAEs at z = 2–7 derived in the previous studies. We correct these previous study measurements of bg for the difference of σ8 (Ouchi et al. 2010), and show the corrected bg values in figure 7. Fig. 7. View largeDownload slide Bias of ≳ L* LAEs as a function of redshift. The red squares represent bias of our LAEs at z = 5.7 and 6.6, while the filled circles denote those of LAEs obtained in Ouchi et al. (2010). The open circles, pentagon, and hexagon indicate LAE bias values obtained by Ouchi et al. (2010, 2005a, 2003), respectively. The filled hexagon, pentagon, and diamond are LAE bias values measured by Guaita et al. (2010), Gawiser et al. (2007), and Kovač et al. (2007), respectively. For the previous study results, we correct for the difference of σ8. To clarify the bias measurements, we give small offsets along the abscissa axis to the data points of the previous studies. The solid lines indicate bias of dark-matter halos with a halo mass of 108, 109, 1010, 1011, and 1012 M⊙ in the case of one-to-one correspondence between galaxies and dark-matter halos (Ouchi et al. 2010). The gray region shows the dark-matter halo mass range of 1010–1012 M⊙. The dotted lines are evolutionary tracks of bias in the case of the galaxy-conserving model [equation (13)]. (Color online) Fig. 7. View largeDownload slide Bias of ≳ L* LAEs as a function of redshift. The red squares represent bias of our LAEs at z = 5.7 and 6.6, while the filled circles denote those of LAEs obtained in Ouchi et al. (2010). The open circles, pentagon, and hexagon indicate LAE bias values obtained by Ouchi et al. (2010, 2005a, 2003), respectively. The filled hexagon, pentagon, and diamond are LAE bias values measured by Guaita et al. (2010), Gawiser et al. (2007), and Kovač et al. (2007), respectively. For the previous study results, we correct for the difference of σ8. To clarify the bias measurements, we give small offsets along the abscissa axis to the data points of the previous studies. The solid lines indicate bias of dark-matter halos with a halo mass of 108, 109, 1010, 1011, and 1012 M⊙ in the case of one-to-one correspondence between galaxies and dark-matter halos (Ouchi et al. 2010). The gray region shows the dark-matter halo mass range of 1010–1012 M⊙. The dotted lines are evolutionary tracks of bias in the case of the galaxy-conserving model [equation (13)]. (Color online) This work and the previous studies measure bias of LAEs that have similar Lyα luminosity detection limits; a few times 1042 erg s−1 corresponding to ∼L* luminosities over z ∼ 2–7 (e.g., Ouchi et al. 2008, 2010; Konno et al. 2016, 2018). Thus, we can omit the difference of bias depending on Lyα luminosity. Figure 7 indicates that the bias of ≳ L* LAEs significantly increases from z ∼ 2–3 to z ∼ 6–7. The physical origins of this increase are discussed in section 5. 5 Discussion 5.1 Cosmic reionization 5.1.1 Constraints on xH i from the clustering results Clustering signals of observed LAEs are strengthened by the patchy distribution of the neutral hydrogen at the EoR, because Lyα photons of LAEs in the ionized bubbles can selectively escape from the patchy neutral IGM with a small amount of Lyα damping wing absorption (Furlanetto et al. 2006; McQuinn et al. 2007; Lidz et al. 2009; Iliev et al. 2008; Sobacchi & Mesinger 2015). Based on the fact that the Gunn–Peterson optical depth largely increases at z ∼ 6 (Fan et al. 2006), one would expect that the clustering amplitude of the observed LAEs increases from z = 5.7 (the ionized universe) to z = 6.6 (the partly-neutral universe). Figure 7 indicates no significant rise of bias from z = 5.7 to 6.6 beyond the moderately small errors (i.e., by a factor of ∼20%). This result suggests that the clustering of z = 6.6 LAEs is not largely affected by cosmic reionization effects, where the bias evolution of the hosting dark-matter halos towards high z may also be involved. Based on this bias evolution result, we place constraints on cosmic reionization parameters with the help of theoretical models. We compare our observational results with multiple theoretical models, because there is a possibility that our conclusions may depend on the model chosen for the comparisons. We thus try to avoid model-dependent conclusions as much as possible, and to obtain more objective interpretations for our observational results. Note that the arguments below follow those of Ouchi et al. (2010), with our HSC clustering measurements. In the top panel of figure 8, we compare our z = 6.6 LAE clustering measurements with those of theoretical predictions (McQuinn et al. 2007; Furlanetto et al. 2006). The model of McQuinn et al. (2007) is presented in the top panel of figure 8. McQuinn et al. (2007) conducted radiative transfer simulations and predicted clustering of LAEs at z = 6.6. Their models assign a Lyα flux to a dark-matter halo whose mass is greater than a minimum dark-matter halo mass. We use the minimum dark-matter halo mass of our LAEs at the post-EoR of z = 5.7 to estimate the intrinsic LAE clustering at z = 6.6, assuming no redshift evolution (z = 5.7–6.6) of the minimum dark-matter halo mass (section 3). Although the minimum dark-matter halo mass is 3 × 109 M⊙ for the z = 5.7 LAEs (subsection 5.2), McQuinn et al. (2007) do not calculate the models of the minimum dark-matter halo mass as low as 3 × 109 M⊙. We thus correct the angular correlation functions of McQuinn et al.’s lowest mass (3 × 1010 M⊙) model for the difference between 3 × 1010 and 3 × 109 M⊙, using equation (9). Here we calculate the bias values of the 3 × 1010 and 3 × 109 M⊙ cases using the best-fitting halo model (subsection 5.2), the parameters of which are exactly the same as those determined at z = 5.7. In the top panel of figure 8, we show McQuinn et al.’s models for 3 × 109 M⊙ with the correction, together with the original 3 × 1010 M⊙ model. Because McQuinn et al.’s models do not correct for the integral constraint, we do not compare the model predictions below ω ∼ 0.1, where the contribution of the integral-constraint correction term is large. The top panel of figure 8 indicates that our z = 6.6 LAE data points fall in the range of xH i ≃ 0–0.3, indicating xH i = 0.15 ± 0.15. Fig. 8. View largeDownload slide Comparisons of the theoretical models with the ACF and bias measurements at z = 6.6. In the top (bottom) panel, the filled squares and circles denote the ACF (bias) measurements with the IC correction obtained by this work and Ouchi et al. (2010), respectively, which are the same as those shown in the top (bottom) right-hand panel of figure 6. Top: The solid and dotted lines represent McQuinn et al.’s (2007) models of z = 6.6 LAEs with a dark-matter halo masses of 3 × 109 M⊙ and 3 × 1010 M⊙, respectively. The model of 3 × 109 M⊙ is the modification of the model of 3 × 1010 M⊙ (see text). From the bottom to top, the solid (dotted) lines, showing the neutral hydrogen fractions of the IGM, are xH i = 0.0, 0.3, 0.5, and 0.8. Bottom: The ticks at the right-hand side indicate bias values predicted by Furlanetto, Zaldarriaga, and Hernquist’s (2006) models (see text). The solid and dashed lines indicate the average bias and the 1σ uncertainty values, respectively, that are the same as those shown in the bottom right-hand panel of figure 6. Fig. 8. View largeDownload slide Comparisons of the theoretical models with the ACF and bias measurements at z = 6.6. In the top (bottom) panel, the filled squares and circles denote the ACF (bias) measurements with the IC correction obtained by this work and Ouchi et al. (2010), respectively, which are the same as those shown in the top (bottom) right-hand panel of figure 6. Top: The solid and dotted lines represent McQuinn et al.’s (2007) models of z = 6.6 LAEs with a dark-matter halo masses of 3 × 109 M⊙ and 3 × 1010 M⊙, respectively. The model of 3 × 109 M⊙ is the modification of the model of 3 × 1010 M⊙ (see text). From the bottom to top, the solid (dotted) lines, showing the neutral hydrogen fractions of the IGM, are xH i = 0.0, 0.3, 0.5, and 0.8. Bottom: The ticks at the right-hand side indicate bias values predicted by Furlanetto, Zaldarriaga, and Hernquist’s (2006) models (see text). The solid and dashed lines indicate the average bias and the 1σ uncertainty values, respectively, that are the same as those shown in the bottom right-hand panel of figure 6. In the bottom panel of figure 8, we compare our bias results with those of the analytical models of Furlanetto, Zaldarriaga, and Hernquist (2006) at z = 6.6. In this comparison, we adopt the bias values of the small scale (≃1−10 Mpc) of Furlanetto, Zaldarriaga, and Hernquist (2006) that is similar to the major angular scale of our bias measurements. For the intrinsic bias value of the z = 6.6 LAEs, we use b = 4.6, which is obtained under the assumption of no evolution from z = 5.7 to 6.6 for the minimum dark-matter halo mass (3 × 109 M⊙) and the halo-model parameters (subsection 5.2). The bottom panel of figure 8 shows that our average-bias measurement agrees with the low-xH i models of Furlanetto, Zaldarriaga, and Hernquist (2006), and suggests xH i ≲ 0.13 including the uncertainty of the average-bias measurement of the z = 6.6 LAEs. The comparisons of our observational results with the models of McQuinn et al. (2007) and Furlanetto, Zaldarriaga, and Hernquist (2006) indicate that the neutral hydrogen fraction at z = 6.6 falls in the range of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$. 5.1.2 Cosmic reionization history Figure 9 shows our xH i value of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ that is estimated with our LAE clustering measurements. This xH i estimate suggests a moderate neutral hydrogen fraction at z = 6.6. In figure 9, we also present various results of xH i constrained by a similar LAE clustering analysis (Ouchi et al. 2010) and the independent xH i estimates from QSO Gunn–Peterson absorptions (Fan et al. 2006), as well as Lyα damping wing (DW) absorption measurements of QSOs (Mortlock et al. 2011; Bolton et al. 2011), gamma-ray bursts (GRBs; Totani et al. 2006, 2014, 2016; see also Greiner et al. 2009), and galaxies including LAEs (Malhotra & Rhoads 2004; Ouchi et al. 2010; Ota et al. 2010; Kashikawa et al. 2011; Pentericci et al. 2011; Ono et al. 2012; Treu et al. 2013; Caruana et al. 2014; Konno et al. 2014, 2018; Schenker et al. 2014). Our results agree with the previous LAE clustering results of Ouchi et al. (2010). Our constraint on xH i is stronger than that of Ouchi et al. (2010), due to the fact that the statistical uncertainties of the ACF measurements of this study are smaller than those of Ouchi et al. (2010). This is because the number of LAEs in this study is larger than that of Ouchi et al. (2010). Moreover, our xH i estimate agrees with many of the estimates and constraints obtained by the other independent objects (QSOs, GRBs, and galaxies) and methods (Gunn–Peterson and Lyα DW absorptions) in the previous studies. It should be noted that the results for GRBs have a large scatter. This scatter is probably made by the systematic uncertainty raised by sightline variance effects. The numerical simulations of McQuinn et al. (2008) suggest that the patchiness of reionization gives the systematic error of δxH i ∼ 0.3 in an xH i measurement of a single GRB sightline. Fig. 9. View largeDownload slide Neutral hydrogen fraction xH i of the IGM as a function of redshift. Note that the top and bottom panels are the same, but with an ordinate axis of linear and log scales, respectively. The red filled square is the xH i estimate obtained by our HSC LAE clustering analysis. The black filled square and circle are the xH i estimates from the LAE LF evolution of Konno et al. (2018, 2014), respectively. The open circles are the constraints at z = 6.6 obtained by Ouchi et al. (2010) from the evolution of Lyα LF (left circle) and clustering (right circle), while the open diamond and the open pentagon represent the upper limits from the Lyα LF evolution to z = 6.5 given by Malhotra and Rhoads (2004) and Kashikawa et al. (2006). Here, we add small offsets along redshift to the positions of the filled square, the open circles, and the open diamond, avoiding overlapping symbols. The filled hexagon and the filled pentagons show the constraints from a spectrum of a GRB (Gallerani et al. 2008b) and statistics of QSO dark-gaps (Gallerani et al. 2008a), respectively. The open hexagons are the constraints calculated from the Lyα damping wing absorption of GRBs at z = 6.3 (Totani et al. 2006) and z = 5.9 (Totani et al. 2016). The filled diamonds indicate the QSO Gunn–Peterson optical depth measurement results (Fan et al. 2006). The triangle denotes the 1σ lower-limit of redshift obtained by Planck 2015 (Planck Collaboration 2016) in the case of instantaneous reionization. The solid line and the gray shade indicate the best estimate and the uncertainty of the xH i evolution (Ishigaki et al. 2017) that agrees the evolutions of τe and ρUV with free parameters including the ionizing photon escape fraction. The dotted, dashed, and dot–dashed lines are the evolution of xH i for the reionizing sources down to the massive halos, the moderately massive halos, and the mini-halos, respectively, in the model of Choudhury, Ferrara, and Gallerani (2008). The dashed–double-dotted line indicates the prediction of the double reionization model (Cen 2003). (Color online) Fig. 9. View largeDownload slide Neutral hydrogen fraction xH i of the IGM as a function of redshift. Note that the top and bottom panels are the same, but with an ordinate axis of linear and log scales, respectively. The red filled square is the xH i estimate obtained by our HSC LAE clustering analysis. The black filled square and circle are the xH i estimates from the LAE LF evolution of Konno et al. (2018, 2014), respectively. The open circles are the constraints at z = 6.6 obtained by Ouchi et al. (2010) from the evolution of Lyα LF (left circle) and clustering (right circle), while the open diamond and the open pentagon represent the upper limits from the Lyα LF evolution to z = 6.5 given by Malhotra and Rhoads (2004) and Kashikawa et al. (2006). Here, we add small offsets along redshift to the positions of the filled square, the open circles, and the open diamond, avoiding overlapping symbols. The filled hexagon and the filled pentagons show the constraints from a spectrum of a GRB (Gallerani et al. 2008b) and statistics of QSO dark-gaps (Gallerani et al. 2008a), respectively. The open hexagons are the constraints calculated from the Lyα damping wing absorption of GRBs at z = 6.3 (Totani et al. 2006) and z = 5.9 (Totani et al. 2016). The filled diamonds indicate the QSO Gunn–Peterson optical depth measurement results (Fan et al. 2006). The triangle denotes the 1σ lower-limit of redshift obtained by Planck 2015 (Planck Collaboration 2016) in the case of instantaneous reionization. The solid line and the gray shade indicate the best estimate and the uncertainty of the xH i evolution (Ishigaki et al. 2017) that agrees the evolutions of τe and ρUV with free parameters including the ionizing photon escape fraction. The dotted, dashed, and dot–dashed lines are the evolution of xH i for the reionizing sources down to the massive halos, the moderately massive halos, and the mini-halos, respectively, in the model of Choudhury, Ferrara, and Gallerani (2008). The dashed–double-dotted line indicates the prediction of the double reionization model (Cen 2003). (Color online) Figure 9 presents the evolution of xH i as determined by Ishigaki et al. (2017), who use all of the observational results related to reionization, CMB Thomson scattering optical depth, τe, UV luminosity density, ρUV, and the ionized fraction, QH ii, obtained to-date, based on the standard analytical model (Madau et al. 1999; Robertson et al. 2015) fitting that allows a wide parameter range of unknown parameters such as the ionizing-photon escape fraction. It should be noted that Ishigaki et al.’s (2017) result of xH i evolution agrees with the τe measurement of Planck2016 (Planck Collaboration 2016). We find in figure 9 that the evolution of xH i (the solid curve) is consistent with our xH i estimate from the LAE clustering, supportng the moderately late reionization scenario suggested by Ishigaki et al. (2017). 5.2 Hosting dark-matter halos 5.2.1 Halo occupation distribution modeling We carry out halo occupation distribution (HOD) modeling for our ACFs where the hosting dark-matter halo properties are encoded. Note that the ACF of the LAEs at z = 5.7 is not affected by the effects of patchy ionized bubbles of cosmic reionization (section 1), because z = 5.7 corresponds to the post-EoR. However, the LAE ACF at z = 6.6 may be strengthened by the selective escape of Lyα photons from patchy ionized bubbles that probably exist at this redshift, although the z = 6.6 LAE ACF appears to be weakly affected by cosmic reionization (sub-subsection 5.1.1). In this section, we perform HOD modeling for the ACF of the LAEs at z = 5.7 for the secure results. Moreover, we conduct the same analysis for the LAEs at z = 6.6 under the assumption that the reionization effect is negligibly small, following our results that the influence of cosmic reionization is not large at z = 6.6 (sub-subsection 5.1.1). We adopt the HOD models of Harikane et al. (2016) based on the halo mass function of Behroozi, Wechsler, and Conroy (2013) that is a modification of the Tinker et al. (2008) halo mass function. In the HOD models of Harikane et al. (2016), there are a total of five parameters: Mmin, $$M_1^{\prime }$$, α, M0, and $$\sigma _{\rm \log M}$$. Here, Mmin is the mass of the dark-matter halo that hosts a galaxy at a possibility of 50%. The terms $$M_1^{\prime }$$ and α represent the normalization and the slope of the power law, respectively, for the occupation number of the satellite galaxy. We define M0 and σlog  M as the halo mass cut-off and the halo mass transition width, respectively, of the central and satellite galaxies. We do not include the parameter of the duty cycle in our HOD modeling, because the number density is not used in our fitting. Instead, we compare the number densities given by our HOD modeling and our observations to constrain the duty cycle of LAEs, which are detailed in sub-subsection 5.2.2. Following Harikane et al. (2016), we adopt the two relations $$\log M_0 = 0.76 \log M_1^{\prime } + 2.3$$ and $$\log M_1^{\prime }=1.18 \log M_{\rm min} - 1.28$$. We assume $$\sigma _{\rm \log M}=0.2$$ and α = 1, which are suggested by the fitting results of Zheng, Coil, and Zehavie (2007), although we confirm that the choice of these two parameters does not change our conclusions. We thus obtain the best-fitting HOD models (i.e., Mmin values) by χ2 minimization. Figure 10 presents the best-fitting HOD models, and table 2 summarizes the best-fitting Mmin values. We define the average dark-matter halo mass ⟨Mh⟩ and the average galaxy bias $$b_{\rm g}^{\rm HOD}$$ from the HOD modeling results;   \begin{eqnarray} \left\langle M_{\rm h} \right\rangle & = & \frac{\int _{M_{\rm h}^{\rm min}}^{\infty } M N_{\rm g} n dM}{\int _{M_{\rm h}^{\rm min}}^{\infty }N_{\rm g} n dM},\ \ \ {\rm and} \end{eqnarray} (10)  \begin{eqnarray} b_{\rm g}^{\rm HOD} & = & \frac{\int _{M_{\rm h}^{\rm min}}^{\infty } b N_{\rm g} n dM}{\int _{M_{\rm h}^{\rm min}}^{\infty } N_{\rm g} n dM}, \end{eqnarray} (11)where n and b are the number density and the bias of dark-matter halos, respectively, with dark-matter halo mass Mh (Behroozi et al. 2013; Tinker et al. 2010). Ng is the galaxy occupation function at Mh that is determined by the HOD modeling. The values of ⟨Mh⟩ and $$b_{\rm g}^{\rm HOD}$$ for our LAEs are summarized in table 2. We find that the average dark-matter halo masses of the ≳ L* LAEs are moderately small, $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) = 11.1^{+0.2}_{-0.4}$$ at z = 5.7 and $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) = 10.8^{+0.3}_{-0.5}$$ at z = 6.6, and are consistent with the previous estimates [log (⟨Mh⟩/M⊙) ∼ 10–11] found via clustering analysis (Ouchi et al. 2010). The values of $$b_{\rm g}^{\rm HOD}$$ are consistent with those of bavg (table 1) within the errors. Fig. 10. View largeDownload slide Best-fitting HOD models of the LAEs at z = 5.7 (left) and 6.6 (right). The solid lines represent the best-fitting HOD models, while the dashed lines denote the breakdowns of the best-fitting HOD models, 1 and 2 halo terms. The squares and circles are the ACF measurements that are the same as those in figure 6. Fig. 10. View largeDownload slide Best-fitting HOD models of the LAEs at z = 5.7 (left) and 6.6 (right). The solid lines represent the best-fitting HOD models, while the dashed lines denote the breakdowns of the best-fitting HOD models, 1 and 2 halo terms. The squares and circles are the ACF measurements that are the same as those in figure 6. Table 2. Properties of dark-matter halos of the LAEs estimated with the HOD model. z  $$L_{\rm Ly\alpha }^{\rm th}$$  log (Mmin/M⊙)  log (⟨M⟩/M⊙)  $$b_{\rm g}^{\rm HOD}$$  $$f_{\rm duty}^{\rm LAE}$$  χ2/dof    (1042 ergz >s−1)        (10−4)    (1)  (2)  (3)  (4)  (5)  (6)  (7)  5.7  6.3  $$9.5^{+0.5}_{-1.2}$$  $$11.1^{+0.2}_{-0.4}$$  $$3.9^{+0.7}_{-1.0}$$  $$7.7^{+53.0}_{-7.6}$$  1.1/6  6.6*  7.9  $$9.1^{+0.7}_{-1.9}$$  $$10.8^{+0.3}_{-0.5}$$  $$4.1^{+1.0}_{-1.4}$$  $$0.9^{+13.5}_{-0.9}$$  1.1/4  z  $$L_{\rm Ly\alpha }^{\rm th}$$  log (Mmin/M⊙)  log (⟨M⟩/M⊙)  $$b_{\rm g}^{\rm HOD}$$  $$f_{\rm duty}^{\rm LAE}$$  χ2/dof    (1042 ergz >s−1)        (10−4)    (1)  (2)  (3)  (4)  (5)  (6)  (7)  5.7  6.3  $$9.5^{+0.5}_{-1.2}$$  $$11.1^{+0.2}_{-0.4}$$  $$3.9^{+0.7}_{-1.0}$$  $$7.7^{+53.0}_{-7.6}$$  1.1/6  6.6*  7.9  $$9.1^{+0.7}_{-1.9}$$  $$10.8^{+0.3}_{-0.5}$$  $$4.1^{+1.0}_{-1.4}$$  $$0.9^{+13.5}_{-0.9}$$  1.1/4  * The results of z = 6.6 LAEs may be weakly affected by the reionization effects (see text). View Large 5.2.2 Duty cycle It is suggested that LAEs do not exist in any dark-matter halos that are more massive than $$M_{\rm h}^{\rm min}$$ (e.g., Ouchi et al. 2010) because LAEs should meet two physical conditions: (i) LAEs should be active star-forming galaxies producing Lyα photons, and (ii) the Lyα photons of LAEs should largely escape from the ISM. In other words, not all galaxies in a dark-matter halo can be LAEs. We evaluate the duty cycle of LAEs (i.e., $$f_{\rm duty}^{\rm LAE}$$) that is defined by the probability that a dark-matter halo with a mass greater than the minimum-halo mass hosts an LAE(s). With fduty, one can calculate the galaxy number density,   $$n_{\rm g} = \int _{M_{\rm h}^{\rm min}}^{\infty } f_{\rm duty}^{\rm LAE}\ N_{\rm g}\ n\ dM.$$ (12)We derive $$f_{\rm duty}^{\rm LAE}$$ with the best-fitting HOD parameters, assuming that $$f_{\rm duty}^{\rm LAE}$$ does not depend on the dark-matter halo mass. The observational measurements of ng are taken from the HSC SILVERRUSH paper of Konno et al. (2018). We thus obtain $$f_{\rm duty}^{\rm LAE}=7.7^{+53.0}_{-0.75} \times 10^{-4}$$ and $$0.9^{+13.5}_{-0.9} \times 10^{-4}$$ for our LAEs at z = 5.7 and 6.6, respectively, which are listed in table 2. Although $$f_{\rm duty}^{\rm LAE}$$ is very poorly constrained with large uncertainties, we find that the duty cycle of LAEs is ∼1% or less at z = 5.7 and 6.6. These small duty cycle values for z = 6–7 LAEs are consistent with the previous estimates found via a similar technique (Ouchi et al. 2010). The small duty cycle values indicate that dark-matter halos hosting LAEs are rare, and that the majority of LAEs are located in abundant low-mass dark-matter halos with a mass close to the minimum-halo mass. 5.2.3 Galaxy formation and LAEs Figure 7 presents the evolutionary tracks of dark-matter halos for the galaxy-conserving model and the constant-mass model. For the galaxy-conserving models, we assume that gravity drives the motion of galaxies under the condition of no mergers. The bias evolution in the galaxy-conserving model (Fry 1996) is described as   $$b_{\rm g}=1+(b_{\rm g}^0-1)/D(z),$$ (13)where D(z) and $$b_{\rm g}^0$$ are the growth factor and the bias at z = 0, respectively. If one assumes the galaxy-conserving evolution, our LAEs at z = 5.7 and 6.6 evolve into the present-day galaxies with $$b_{\rm g}^0\sim 1.6$$–1.7. Galaxies with $$b_{\rm g}^0\sim 1.6$$–1.7 are massive, bright ∼6L* galaxies today (Zehavi et al. 2005). The galaxy-conserving evolution suggests that our LAEs at z = 5.7 and 6.6 are the progenitors of very massive bright galaxies in the present-day universe. 6 Summary We investigate clustering of Lyα emitters (LAEs) at z = 5.7 and 6.6 with the early data from the Hyper Suprime-Cam (HSC) Subaru Strategic Program (first-data release shown in Aihara et al. 2018), and introduce the program of Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH). From the early data, we obtain an unprecedentedly large sample of 1832 ≳ L* LAEs at z = 6–7 over the total area of 14–21 deg2, which is about an order of magnitude larger than the previous z = 6–7 LAE clustering studies. Based on the LAE clustering measurements, we study cosmic reionization and galaxy formation, comparing theoretical models. In this study, there are two major results, which are summarized below. 1. We calculate angular correlation functions of the z = 6–7 LAEs. The correlation lengths are estimated to be $$r_0=3.01^{+0.35}_{-0.35}$$ and $$2.66^{+0.49}_{-0.70}\,\,h^{-1}_{100}\:$$Mpc for the ≳ L* LAEs at z = 5.7 and 6.6, respectively. The average of the large-scale bias value is bavg = 4.13 ± 0.17 (4.54 ± 0.63) at z = 5.7 (z = 6.6) for the LAEs. Because Lyα photons emitted from LAEs in ionized bubbles can selectively escape from the partly neutral IGM at the EoR, observed LAEs at z = 6.6 should have clustering signals stronger than the intrinsic clustering. Based on this physical picture, we obtain the constraint of $$x_{\rm H\,{\small I}}=0.15^{+0.15}_{-0.15}$$ at z = 6.6 from the comparisons between our clustering measurements and two independent theoretical models. 2. We study the ≳ L* LAE clustering via HOD modeling. (Here, for the LAE clustering at the EoR of z = 6.6, we assume that the LAE clustering is not largely impacted by the cosmic reionization.) The best-fitting models indicate that the ≳ L* LAEs are hosted by the dark-matter halos with average masses of $$\log (\left\langle M_{\rm h} \right\rangle /M_\odot ) =11.1^{+0.2}_{-0.4}$$ and $$10.8^{+0.3}_{-0.5}$$ at z = 5.7 and 6.6, respectively. Comparing the number densities of observed LAEs and those suggested from the HOD modeling, we find that dark-matter halos of only 1% (or less) down to the minimum-halo mass can host the ≳ L* LAEs. With the standard structure formation models, the bias evolution of the dark-matter halos indicates that the ≳ L* LAEs at z = 6–7 are progenitors of massive ∼6L* galaxies in the present-day universe. Acknowledgements We are grateful to useful discussion with Mark Dijikstra, Richard Ellis, Andrea Ferrara, Martin Haehnelt, Alex Hagen, Koki Kakiichi, Andrei Mesinger, Naveen Reddy, and Zheng Zheng. We acknowledge Jirong Mao and Anne Hutter for their comments. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. The NB816 filter was supported by Ehime University (PI: Y. Taniguchi). The NB921 filter was supported by KAKENHI (23244025) Grant-in-Aid for Scientific Research (A) through the Japan Society for the Promotion of Science (PI: M. Ouchi). This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory. This work is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (15H02064) Grant-in-Aid for Scientific Research (A) through Japan Society for the Promotion of Science. 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 at National Astronomical Observatory of Japan. Footnotes † Based on data obtained with the Subaru Telescope. The Subaru Telescope is operated by the National Astronomical Observatory of Japan. 1 There are some more OH windows redder than 921 nm, which include the 1010-nm window for the NB101 filter. However, these OH windows in the red band are contaminated by moderately strong OH lines, unlike the one at 921 nm. References Aihara H. et al.   2018, PASJ , 70, S8 Bacon R. et al.   2010, SPIE Proc. , 7735, 773508 Behroozi P. S., Wechsler R. H., Conroy C. 2013, ApJ , 770, 57 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 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   Cen R. 2003, ApJ , 591, 12 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., Barger A. J., Hu E. M. 2010, ApJ , 711, 928 CrossRef Search ADS   Cowie L. L., Hu E. M. 1998, AJ , 115, 1319 CrossRef Search ADS   Deharveng J.-M. et al.   2008, ApJ , 680, 1072 CrossRef Search ADS   Efstathiou G., Bernstein G., Tyson J. A., Katz N., Guhathakurta P. 1991, ApJ , 380, L47 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  Finkelstein S. L., Cohen S. H., Moustakas J., Malhotra S., Rhoads J. E., Papovich C. 2011, ApJ , 733, 117 CrossRef Search ADS   Fry J. N. 1996, ApJ , 461, L65 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L. 2006, MNRAS , 365, 1012 CrossRef Search ADS   Gallerani S., Ferrara A., Fan X., Choudhury T. R. 2008, MNRAS , 386, 359 CrossRef Search ADS   Gallerani S., Salvaterra R., Ferrara A., Choudhury T. R. 2008, MNRAS , 388, L84 CrossRef Search ADS   Gawiser E. et al.   2007, ApJ , 671, 278 CrossRef Search ADS   Greiner J. et al.   2009, ApJ , 693, 1610 CrossRef Search ADS   Groth E. J., Peebles P. J. E. 1977, ApJ , 217, 385 CrossRef Search ADS   Guaita L. et al.   2010, ApJ , 714, 255 CrossRef Search ADS   Guaita L., Francke H., Gawiser E., Bauer F. E., Hayes M., Östlin G., Padilla N. 2013, A&A , 551, A93 CrossRef Search ADS   Hagen A. et al.   2014, ApJ , 786, 59 CrossRef Search ADS   Harikane Y. et al.   2016, ApJ , 821, 123 CrossRef Search ADS   Harikane Y. et al.   2018, PASJ , 70, S11 Hashimoto T. et al.   2015, ApJ , 812, 157 CrossRef Search ADS   Hayashino T. et al.   2004, AJ , 128, 2073 CrossRef Search ADS   Hill G. J. et al.   2012, SPIE Proc. , 8446, 84460N 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., McMahon R. G. 1996, Nature , 382, 231 CrossRef Search ADS   Hutter A., Dayal P., Müller V. 2015, MNRAS , 450, 4025 CrossRef Search ADS   Iliev I. T., Shapiro P. R., McDonald P., Mellema G., Pen U.-L. 2008, MNRAS , 391, 63 CrossRef Search ADS   Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K. 2017, arXiv:1702.04867 Kakiichi K., Dijkstra M., Ciardi B., Graziani L. 2016, MNRAS , 463, 4019 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   Kobayashi M. A. R., Totani T., Nagashima M. 2010, ApJ , 708, 1119 CrossRef Search ADS   Kojima T., Ouchi M., Nakajima K., Shibuya T., Harikane Y., Ono Y. 2017, PASJ , 69, 44 CrossRef Search ADS   Konno A. et al.   2014, ApJ , 797, 16 CrossRef Search ADS   Konno A. et al.   2018, PASJ , 70, S16 Konno A., Ouchi M., Nakajima K., Duval F., Kusakabe H., Ono Y., Shimasaku K. 2016, ApJ , 823, 20 CrossRef Search ADS   Kovač K., Somerville R. S., Rhoads J. E., Malhotra S., Wang J. 2007, ApJ , 668, 15 CrossRef Search ADS   Landy S. D., Szalay A. S. 1993, ApJ , 412, 64 CrossRef Search ADS   Le Delliou M., Lacey C. G., Baugh C. M., Morris S. L. 2006, MNRAS , 365, 712 CrossRef Search ADS   Lee K.-S., Giavalisco M., Gnedin O. Y., Somerville R. S., Ferguson H. C., Dickinson M., Ouchi M. 2006, ApJ , 642, 63 CrossRef Search ADS   Lidz A., Zahn O., Furlanetto S. R., McQuinn M., Hernquist L., Zaldarriaga M. 2009, ApJ , 690, 252 CrossRef Search ADS   McLure R. J., Cirasuolo M., Dunlop J. S., Foucaud S., Almaini O. 2009, MNRAS , 395, 2196 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S. 2007, MNRAS , 381, 75 CrossRef Search ADS   McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Dutta S. 2008, MNRAS , 388, 1101 Madau P., Haardt F., Rees M. J. 1999, ApJ , 514, 648 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   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   Miyazaki S. et al.   2012, Proc. SPIE , 8446, 84460Z Momose R. et al.   2014, MNRAS , 442, 110 CrossRef Search ADS   Momose R. et al.   2016, MNRAS , 457, 2318 CrossRef Search ADS   Mortlock D. J. et al.   2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murayama T. et al.   2007, ApJS , 172, 523 CrossRef Search ADS   Nagamine K., Ouchi M., Springel V., Hernquist L. 2010, PASJ , 62, 1455 CrossRef Search ADS   Nakajima K. et al.   2012, ApJ , 745, 12 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   Ono Y. et al.   2012, ApJ , 744, 83 CrossRef Search ADS   Ono Y. et al.   2018, PASJ , 70, S10 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.   2005a, ApJ , 620, L1 CrossRef Search ADS   Ouchi M. et al.   2005b, ApJ , 635, L117 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   Partridge R. B., Peebles P. J. E. 1967, ApJ , 147, 868 CrossRef Search ADS   Pascarelle S. M., Windhorst R. A., Keel W. C., Odewahn S. C. 1996, Nature , 383, 45 CrossRef Search ADS   Peacock J. A., Dodds S. J. 1994, MNRAS , 267, 1020 CrossRef Search ADS   Peebles P. J. E. 1980, The Large-Scale Structure of the Universe  ( Princeton: Princeton University Press) Pentericci L. et al.   2011, ApJ , 743, 132 CrossRef Search ADS   Planck Collaboration 2016, A&A , 594, A13 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   Saito T., Shimasaku K., Okamura S., Ouchi M., Akiyama M., Yoshida M. 2006, ApJ , 648, 54 CrossRef Search ADS   Santos S., Sobral D., Matthee J. 2016, MNRAS , 463, 1678 CrossRef Search ADS   Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P. 2014, ApJ , 795, 20 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   Shimasaku K. et al.   2003, ApJ , 586, L111 CrossRef Search ADS   Simon P. 2007, A&A , 473, 711 CrossRef Search ADS   Sobacchi E., Mesinger A. 2015, MNRAS , 453, 1843 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   Stark D. P. et al.   2015, MNRAS , 454, 1393 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   Steidel C. C., Bogosavljević M., Shapley A. E., Kollmeier J. A., Reddy N. A., Erb D. K., Pettini M. 2011, ApJ , 736, 160 CrossRef Search ADS   Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E. 2008, ApJ , 688, 709 CrossRef Search ADS   Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S. 2010, ApJ , 724, 878 CrossRef Search ADS   Toshikawa J. et al.   2018, PASJ , 70, S12 Totani T. et al.   2014, PASJ , 66, 63 CrossRef Search ADS   Totani T., Aoki K., Hattori T., Kawai N. 2016, PASJ , 68, 15 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., Schmidt K. B., Trenti M., Bradley L. D., Stiavelli M. 2013, ApJ , 775, L29 CrossRef Search ADS   Verhamme A., Schaerer D., Atek H., Tapken C. 2008, A&A , 491, 89 CrossRef Search ADS   Zehavi I. et al.   2005, ApJ , 630, 1 CrossRef Search ADS   Zheng Z., Cen R., Weinberg D., Trac H., Miralda-Escudé J. 2011, ApJ , 739, 62 CrossRef Search ADS   Zheng Z., Coil A. L., Zehavi I. 2007, ApJ , 667, 760 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 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

## 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
that matters to you.

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. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations