SILVERRUSH. XI. Constraints on the Lyα Luminosity Function and Cosmic Reionization at z = 7.3 with Subaru/Hyper Suprime-Cam

Cosmic reionization is a key process in the early universe where massive stars and/or active galactic nuclei ionized the intergalactic medium (IGM) hydrogen that had been neutral after recombination at redshift (z) ∼ 1100. Understanding how and when this process occurred is one of the major goals of modern cosmology and astronomy.

The Thomson scattering optical depth of the cosmic microwave background (CMB) suggests a midpoint reionization redshift of ${z}_{\mathrm{re}}=7.7\pm 0.7$ (Planck Collaboration et al. 2020). Furthermore, recent observations of various kinds of distant objects beyond z ∼ 6 have constrained the period of reionization by estimating the neutral hydrogen fraction in the IGM, x_{H I}, as a function of redshift. Gunn–Peterson troughs in quasi-stellar object (QSO) spectra suggest that cosmic reionization has been completed by z ∼ 6 (e.g., Fan et al. 2006). Damping-wing signatures in QSO spectra (e.g., Schroeder et al. 2013; Greig et al. 2017; Bañados et al. 2018; Davies et al. 2018; Greig et al. 2019; Wang et al. 2020) and gamma-ray burst (GRB) spectra (e.g., Totani et al. 2006, 2014) have also placed constraints on x_{H I} at z ≃ 6–7.5, although the total number of observed sources is very limited.

Lyα emission from galaxies is also a powerful probe of x_{H I} because galaxies are much more numerous than QSOs and GRBs. Methods using galaxies' Lyα emission include the fraction of Lyman-break galaxies (LBGs) that emit Lyα (e.g., Stark et al. 2011; Ono et al. 2012; Mesinger et al. 2015), the Lyα equivalent-width (EW) distribution of LBGs (e.g., Mason et al. 2018a; Hoag et al. 2019; Mason et al. 2019; Jung et al. 2020; Whitler et al. 2020), and the Lyα luminosity function (LF) of Lyα emitters (LAEs; i.e., galaxies with strong Lyα emission; e.g., Kashikawa et al. 2006; Ouchi et al. 2010; Kashikawa et al. 2011; Konno et al. 2014; Ota et al. 2017; Zheng et al. 2017; Itoh et al. 2018; Konno et al. 2018; Hu et al. 2019; Wold et al. 2021).

The observations mentioned above collectively suggest that the universe is significantly neutral at z ≳ 7. However, the x_{H I} estimates at z > 7 still have a large scatter, spanning x_{H I} ∼0.2–0.9, perhaps suggesting field-to-field variation or the presence of a systematic uncertainty in each method. To further constrain the time evolution of x_{H I}, we need to accumulate estimates by individual methods.

In this study, we focus on the Lyα LF method. This method estimates x_{H I} by comparing an observed Lyα LF of LAEs at a target redshift with that after completion of reionization (e.g., z = 5.7). This method's advantage is that x_{H I} can be estimated with a negligibly small redshift uncertainty for a large cosmic volume if an LAE sample from a large-area narrowband (NB) survey is used. A drawback is that the effect of galaxy evolution on observed LFs has to be eliminated using the UV LF of LBGs or a theoretical model.

The Lyα LF in the reionization era has been obtained at z = 5.7, 6.6, 7.0, and 7.3. At z = 7.3, the highest redshift that can be probed with optical CCD detectors, Konno et al. (2014) have obtained x_{H I} = 0.3–0.8 from an accelerated decline of the Lyα LF from z = 5.7. However, because of a relatively small survey area, they have derived only the faint (L_Lyα <10^42.9 erg s⁻¹) part of the Lyα LF. This is contrasted to the studies at lower redshifts that cover up to L_Lyα ≳ 10^43.5 erg s⁻¹ thanks to a large survey volume of >1 × 10⁶ Mpc³. To obtain a robust x_{H I} estimate at z = 7.3, we need to compare the entire LF including the bright part with that at z = 5.7.

In this paper, we present the results of a new survey of bright z = 7.3 LAEs with Subaru/Hyper Suprime-Cam (HSC; Miyazaki et al. 2012; Furusawa et al. 2018; Kawanomoto et al. 2018; Komiyama et al. 2018; Miyazaki et al. 2018), conducted as part of the SILVERRUSH project (Harikane et al.2018; Inoue et al. 2018; Konno et al. 2018; Ouchi et al. 2018; Shibuya et al. 2018a, 2018b; Harikane et al. 2019; Higuchi et al. 2019; Kakuma et al. 2021; Y. Ono et al. 2021, in preparation) that uses four HSC NB filters to study LAEs. Because our survey volume is as large as ∼ 2 × 10⁶ Mpc³, which is seven times larger than that of Konno et al. (2014), our constraint on x_{H I} should also be robust against the uncertainty due to spatially inhomogeneous reionization (see Section 2.2). To infer x_{H I}, we first need to calculate the Lyα transmission of the IGM, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$. We use a new method to calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ and apply it also to the previous studies' Lyα LFs at z = 6.6 (Konno et al. 2018), 7.0 (Itoh et al. 2018; Hu et al. 2019), and 7.3 (Konno et al. 2014).

This paper is structured as follows. Section 2 describes the HSC imaging data used in this study. Section 3 describes our LAE selection. Section 4 presents new constraints on the Lyα LF, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$, and x_{H I}, including a comparison with x_{H I} estimates by other methods. Section 5 is devoted to conclusions.

Throughout this paper, we assume a flat ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. Magnitudes are given in the AB system (Oke & Gunn 1983). Distances are in comoving units unless otherwise noted.

In this work, we use the images in the HSC Subaru Strategic Program ¹¹ (SSP; Aihara et al. 2018) S19A release data, which are reduced with hscPipe 7 (Bosch et al. 2018). ¹² To search for LAEs at z = 7.3, we use the custom-made NB filter NB1010.

2.1. NB1010 Filter

The NB1010 filter has a central wavelength of λ_c = 10092 Å and an FWHM of 91 Å to identify the Lyα emission at z = 7.30 ± 0.04. At a given observing time, a filter with a narrower FWHM can detect fainter Lyα emission. With the interference coating technique when this filter was manufactured, the above FWHM value was practically the narrowest that could be achieved around 10000 Å in the very fast (F/2.25) incoming beam of the HSC. Figure 1 shows the transmission curves of the NB and broadband (BB) filters used in this study (NB1010, y, z, i, r, g, NB921, and NB816).

Zoom In Zoom Out Reset image size

2.2. Images

The details of the NB and BB imaging data used in this study are summarized in Table 1. The NB1010 observations were carried out between 2018 February and 2019 January in two fields, COSMOS and SXDS. The total exposure time is 13.6 hr in COSMOS and 14.0 hr in SXDS, respectively.

Table 1. Summary of Imaging Data

Field	Band	Area a	Exposure Time	PSF Size b	${m}_{\mathrm{lim}}$ c	Dates of Observation
		(deg2)	(s)	(arcsec)	(mag)
COSMOS	NB1010	1.55	49,122	0.69	24.7	2018 Feb 12, 13, 21; 2019 Jan 5
	y			0.70	26.4
	z			0.59	27.1
	i			0.64	27.5
	r			0.67	27.7
	g			0.80	28.2
	NB921			0.64	26.4
	NB816			0.63	26.5
SXDS	NB1010	1.47	50,400	0.72	24.5	2018 Feb 12, 13, 21; 2019 Jan 4, 5, 9
	y			0.58	25.5
	z			0.56	26.3
	i			0.62	27.0
	r			0.62	27.1
	g			0.66	27.6
	NB921			0.69	26.2
	NB816			0.58	26.3

Notes.

^a The effective survey area after removing masked regions (Section 2.2 and Figure 2). ^b Values calculated at patches in the central region of the field of view, 9813–5, 4 (tract-patch) in COSMOS and 8523–2,6 in SXDS (see Figure 3). ^c The 5σ limiting magnitude in an aperture with a diameter of two times the PSF FWHM, 141 (COSMOS) and 144 (SXDS).

Download table as:  ASCIITypeset image

To mask out regions around bright stars, we use the mask images provided by the pipeline. In our analysis, we do not use pixels that have either the SAT, BRIGHT_OBJECT, or NO_DATA flag. ¹³ We also remove low signal-to-noise ratio (S/N) regions near the edges of the images. After removal of these masked regions, the effective survey areas are 1.55 deg² and 1.47 deg² in the COSMOS and SXDS fields, respectively. Figure 2 shows the effective survey areas of the two fields. Assuming a top-hat NB filter with the FWHM of NB1010, the survey volumes are 8.79 × 10⁵ Mpc³ and 8.35 × 10⁵ Mpc³ in COSMOS and SXDS, respectively. At z = 7.3 and with x_{H I} ∼ 0.1, the volume of typical ionized bubbles is estimated to be ∼ 3 × 10⁵ Mpc³ using an analytic model by Furlanetto & Oh (2005). Because our total survey volume is ∼5 times larger than this, our constraint on x_{H I} should be robust against the uncertainty due to spatially inhomogeneous reionization.

Zoom In Zoom Out Reset image size

In the HSC-SSP data processing, the sky is divided into grids called tracts, and each tract is further divided into subareas called patches. Each patch covers approximately $12^{\prime} \times 12^{\prime}$ of the sky (Aihara et al. 2018). We conduct LAE selection at each patch, using the local limiting magnitude estimated at that patch. To do so, we estimate limiting magnitudes at all patches, as shown in Figure 3, by placing in the unmasked region random apertures whose diameter is two times the point-spread function (PSF) FWHM averaged over the image, 141 (COSMOS) and 144 (SXDS). For each field, the limiting magnitude gradually becomes brighter toward the edge of the image. At patches in the central region, the NB1010 images have seeing sizes of 069 (COSMOS) and 072 (SXDS) and reach 5σ limiting magnitudes of 24.7 mag (COSMOS) and 24.5 mag (SXDS).

Zoom In Zoom Out Reset image size

3.1. Source Detection and Photometry

3.2. LAE Selection

We select z = 7.3 LAE candidates based on (1) significant detection in the NB1010 image, (2) NB color excess due to the Lyα emission, y − NB1010, and (3) no detection in the bluer bands to exclude foreground galaxies. The exact selection criteria are as follows:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{NB}1010 & \lt & \mathrm{NB}{1010}_{5\sigma },\\ y-\mathrm{NB}1010 & \gt & 1.9,\\ z\gt {z}_{3\sigma },i\gt {i}_{3\sigma },r & \gt & {r}_{3\sigma },g\gt {g}_{3\sigma },\\ \mathrm{NB}921 & \gt & \mathrm{NB}{921}_{3\sigma },\mathrm{NB}816\gt \mathrm{NB}{816}_{3\sigma },\end{array}\end{eqnarray} \tag{ 1 }$$

where NB1010_5σ is the 5σ limiting magnitude of NB1010 and [z_3σ , i_3σ , r_3σ , g_3σ , NB921_3σ , and NB816_3σ ] are the 3σ limiting magnitudes of the [z, i, r, g, NB921, and NB816] bands. Note that we use the limiting magnitude estimated at the patch in which the object exists (see Section 2.2). We use aperture magnitudes, MAG_APER, (Section 3.1), to measure S/Ns and colors. To measure colors accurately, we convolve the y image of the SXDS field to have the same PSF size as the NB1010 image.

To determine the y − NB1010 color criterion above, we calculate the expected colors of z = 7.3 LAEs. We assume a simple model spectrum that has a flat continuum (f_ν = const., i.e., the UV continuum slope β = − 2) ¹⁴ and δ-function Lyα emission with rest-frame EWs of EW₀ = 0, 10, 20, 30, 50, 150, and 300 Å. Then we redshift the spectra and apply IGM absorption (Madau 1995). ¹⁵ The colors of the spectra are calculated with the transmission curves of the HSC filters shown in Figure 1. Figure 4 shows the results of the expected colors as a function of redshift. Based on Figure 4, we adopt y − NB1010 > 1.9 as our color criterion for z = 7.3 LAEs, corresponding to EW₀ ≳ 10 Å. Because we want to see the evolution from z = 5.7 to estimate x_{H I}, we adopt the same EW limit as Konno et al.'s (2018) z = 5.7 LAEs. This EW limit is also the same as those adopted in Konno et al. (2018) for z = 6.6 LAEs and Itoh et al. (2018) and Hu et al. (2019) for z = 7.0 LAEs.

Zoom In Zoom Out Reset image size

We apply the selection criteria to the objects detected in Section 3.1. Then we perform visual inspection of the objects that pass the selection criteria. Spurious sources such as cosmic rays, CCD artifacts, and artificial diffuse objects outside the masked regions are removed. Example images of the spurious sources are shown in Appendix A. After the visual inspection, there are no LAE candidates left in either the COSMOS or SXDS field.

3.3. Sample Incompleteness

To estimate what fraction of real LAEs pass our selection, we insert pseudo-LAEs into the NB1010 image of each field, and then calculate "detection completeness" and "selection completeness".

3.3.1. Pseudo-LAEs

We use GALSIM (Rowe et al. 2015) to simulate pseudo-LAEs. The pseudo-LAEs have a Sérsic index of n = 1.0 and a half-light radius of r_e ∼ 0.8 kpc (physical units), which corresponds to 016 at z = 7.3. These values are consistent with those of z ∼ 7 LBGs (e.g., Shibuya et al. 2015; Kawamata et al. 2018) at M_UV ≲ −21, corresponding to the luminosity limit of this study, $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]\gtrsim 43.2$, and typical rest-frame Lyα EWs at this redshift, EW₀ ≲ 100 Å (e.g., Hashimoto et al. 2019). Previous studies have also adopted similar values (Itoh et al. 2018; Konno et al. 2018; Hu et al. 2019).

Most LAEs have an extended Lyα halo component (e.g., Momose et al. 2016; Leclercq et al. 2017) that cannot be detected in NB images. Hu et al. (2019) have simulated pseudo-LAEs with larger half-light radii of 0.9, 1.2, and 1.5 kpc (physical units), taking account of the total Lyα emission (main body plus halo) based on MUSE observations of z = 3–6 LAEs by Leclercq et al. (2017) and found negligible differences in the completeness measurements with these radii. Because the sizes and luminosities of Lyα halo components at z > 6 are yet to be examined, we adopt r_e ∼ 0.8 kpc (physical units), which is consistent with the values assumed in the previous studies. We apply PSF convolution to the pseudo-LAEs and randomly insert them into the NB1010 images avoiding the masked regions.

Drake et al. (2017) have shown with MUSE data of LAEs over 3 ≲ z ≲ 6 that extended Lyα emission affects the detection completeness of LAEs. However, if the ratio of the extended component to the total luminosity does not evolve with redshift, all NB surveys will be underestimating the total luminosity and incompleteness in the same way, which does not affect x_{H I} estimates. Indeed, Figure 4 of Drake et al. (2017) shows that the contribution of the extended component is almost the same regardless of redshift and brightness.

3.3.2. Detection Completeness and Selection Completeness

We perform source detection and photometry for the pseudo-LAEs with SExtractor in exactly the same manner as in Section 3.1 to calculate detection completeness. We define detection completeness as the fraction of detected pseudo-LAEs to the input pseudo-LAEs.

We also calculate selection completeness, which has been introduced in Hu et al. (2019), to account for the effects of not meeting the LAE selection criteria because of foreground sources. Some LAEs may be blended with foreground sources that are not bright enough in NB1010 to hinder the LAEs' detection but bright enough in y or bluer bands to prevent them from passing the selection defined as lines 2–4 of Equation (1). To calculate this selection completeness, we assume underlying broadband fluxes to be zero and insert pseudo-LAEs only into the NB1010 images, following Hu et al. (2019; see Section 4.1 of Hu et al. 2019 for more details). Selection completeness is defined as the fraction of pseudo-LAEs that meet the selection criteria (lines 2–4 of Equation (1)) to the detected pseudo-LAEs.

As an example, Figure 5 shows the results for the 9813–4, 4 (COSMOS) and 8523–2,6 (SXDS) patches in the central regions. Their detection completeness and selection completeness are ≳90% and ∼70%–80%, respectively, at magnitudes brighter than the 5σ limiting magnitude. Note that selection completeness, which has not been considered in previous studies except in Hu et al. (2019), is more dominant (i.e., lower) than detection completeness at input magnitudes ≲25 mag. The results of our selection completeness, ∼70%–80% at magnitudes brighter than the 5σ magnitude, are similar to those of Hu et al. (2019) and are mainly due to foreground contamination in bluer bands. Hu et al. (2019) have estimated the effect of blending with foreground sources in bluer bands by random aperture photometry. For example, they have found that only a 74.9% area of the COSMOS HSC g-band image has S/N < 3σ with an aperture size of 135.

Zoom In Zoom Out Reset image size

As found in Figure 5, selection completeness has a mild peak near the 5σ limiting magnitude. As explained by Hu et al. (2019), fainter LAEs would not be detected in the detection image (NB1010 in our case) if blended with a foreground source, which results in lower detection completeness. Consequently, selection completeness gradually increases toward fainter magnitudes because nondetected pseudo-LAEs, blended with a foreground source, are preexcluded from the calculation, i.e., most of the detected pseudo-LAEs are located in sparse regions. On the other hand, selection completeness drops at magnitudes fainter than 5σ magnitude, as most of these LAEs are detected simply because they happen to be injected on top of a foreground source.

Note that we only consider detection completeness when calculating upper limits of the cumulative Lyα LF (Section 4.1) to directly compare them with the Lyα LF measurements in previous studies, which have not considered selection completeness. We use the detection completeness averaged over the effective area, 0.95 (COSMOS) and 0.96 (SXDS).

4.1. Cumulative Lyα Luminosity Function

From the result of no detection of z = 7.3 LAEs (Section 3.2), we calculate the upper limits of the cumulative Lyα LF. The upper limit of the cumulative number density at a given L_Lyα is calculated as

$$\begin{eqnarray}&&n(\gt {L}_{\mathrm{Ly}\alpha })\lt \displaystyle \frac{1.15}{{V}_{\mathrm{eff},1}\,{f}_{\det ,1}+{V}_{\mathrm{eff},2}\,{f}_{\det ,2}},\end{eqnarray} \tag{ 2 }$$

where 1.15 corresponds to the 68% upper limit for no detection assuming Poisson statistics, V_eff is the total survey volume of patches whose limiting luminosity is fainter than this L_Lyα (i.e., patches that allow an LAE search at > L_Lyα ), and f_det is the detection completeness derived in Section 3.3.2, 0.95 (COSMOS) and 0.96 (SXDS). The subscripts 1 and 2 in Equation (2) represent the COSMOS and SXDS fields, respectively.

In Figure 6, we show upper limits for three L_Lyα values: $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]=43.19$ (23 patches), 43.23 (56), and 43.27 (75). The number in the parentheses is the number of patches used, which is dependent on L_Lyα because different patches have different limiting magnitudes (Figure 3). We have searched for L_Lyα (and the corresponding effective survey area) that gives the most stringent upper limit of ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$, finding that $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]=43.23$ (the data point in the middle) is the one. The results of the other two luminosities are plotted to show how much the upper limit changes with a slight change of L_Lyα .

Zoom In Zoom Out Reset image size

To estimate Lyα-limiting luminosities from the limiting magnitudes, we assume spectral energy distributions that have a flat (f_ν = const.) continuum, δ-function Lyα emission with EW₀ = 100 Å, and zero flux at the wavelength bluer than Lyα due to the IGM absorption. This EW₀ value results in a conservative estimate of L_Lyα (see, e.g., Hashimoto et al. 2019 for the EW distribution of z ∼ 6–8 LAEs).

Figure 6 shows the upper limits from this study together with the cumulative Lyα LFs of previous studies. Our results are the first constraints on the bright ($\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]\gtrsim 43$) part of the z = 7.3 LF, making it possible to evaluate the IGM transmission using bright LAEs. Our upper limits show a decrease from the Lyα LFs at z = 7.0 derived by Itoh et al. (2018) (pink solid line and filled circles) and Hu et al. (2019) (pink dashed line and open triangles).

4.2. IGM Transmission to Lyα Photons

In this section, we derive the transmission of Lyα through the IGM, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)$, from the luminosity decrease of the Lyα LF.

The evolution of the Lyα LF is a combination of two effects: galaxy evolution (i.e., the intrinsic evolution of LAEs) and the change in ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ due to cosmic reionization. To obtain implications for cosmic reionization, we need to resolve the degeneracy of these two effects. Ouchi et al. (2010) evaluated the effect of galaxy evolution using the UV LF evolution of LBGs. The UV LF of LBGs also decreases from z ∼ 6 to z ∼ 8 (e.g., Bouwens et al. 2015; Finkelstein et al. 2015), suggesting that the cosmic star formation rate of galaxies declines over this redshift range. In this study, we also estimate the effect of galaxy evolution using the same idea.

We assume that the observed L_Lyα of galaxies can be written using their L_UV as

$$\begin{eqnarray}&&{L}_{\mathrm{Ly}\alpha }={T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z){f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }\kappa {L}_{\mathrm{UV}},\end{eqnarray} \tag{ 3 }$$

where ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ is the Lyα escape fraction through the interstellar medium of galaxies and κ is the Lyα production rate per UV luminosity. The assumption that ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ is independent of intrinsic Lyα luminosity means that the observed Lyα luminosities of galaxies are uniformly decreased by IGM absorption, i.e., IGM absorption does not change the shape of the Lyα LF and only changes the characteristic luminosity, L^*. We also assume that ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ and κ do not change with redshift or UV luminosity, implying that the intrinsic Lyα LF of LAEs evolves in the same manner as the UV LF of LBGs.

To estimate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$, we first predict the Lyα LF with the fully ionized IGM from the evolution of the UV LF. A Schechter function (Schechter 1976) is defined by

$$\begin{eqnarray}&&\phi (L){dL}={\phi }^{* }{\left(L/{L}^{* }\right)}^{\alpha }\exp (-L/{L}^{* })d(L/{L}^{* }),\end{eqnarray} \tag{ 4 }$$

where L^* is the characteristic luminosity, ϕ^* is the characteristic number density, and α is the faint-end slope. We calculate the Schechter parameters of the predicted Lyα LF with the fully ionized IGM as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{Ly}\alpha }^{* \mathrm{pred}}(z) & = & {L}_{\mathrm{Ly}\alpha }^{* \mathrm{obs}}(5.7)\times {L}_{\mathrm{UV}}^{* \mathrm{obs}}(z)/{L}_{\mathrm{UV}}^{* \mathrm{obs}}(5.7),\\ {\phi }_{\mathrm{Ly}\alpha }^{* \mathrm{pred}}(z) & = & {\phi }_{\mathrm{Ly}\alpha }^{* \mathrm{obs}}(5.7)\times {\phi }_{\mathrm{UV}}^{* \mathrm{obs}}(z)/{\phi }_{\mathrm{UV}}^{* \mathrm{obs}}(5.7),\\ {\alpha }_{\mathrm{Ly}\alpha }^{\mathrm{pred}}(z) & = & {\alpha }_{\mathrm{Ly}\alpha }^{\mathrm{obs}}(5.7)+{\alpha }_{\mathrm{UV}}^{\mathrm{obs}}(z)-{\alpha }_{\mathrm{UV}}^{\mathrm{obs}}(5.7),\end{array}\end{eqnarray} \tag{ 5 }$$

where the superscripts "pred" and "obs" mean predicted and observed values, respectively. This equation assumes that the UV LF evolves as described by the empirical model of Bouwens et al. (2015) (the first equation in their Section 5.1). Specifically, we assume that ${L}_{\mathrm{Ly}\alpha }^{* }$ and ${\phi }_{\mathrm{Ly}\alpha }^{* }$ increase or decrease in the same ratio as those of the UV LF and that α_Lyα increases or decreases additively in the same way as the UV LF. For these calculations, we use the Schechter parameters of Konno et al. (2018) for the observed Lyα LF at z = 5.7 and Bouwens et al. (2015) for the observed UV LFs. Figure 7 shows a comparison between the predicted and observed Lyα LFs at z = 6.6, 7.0, and 7.3. We find that the observed Lyα LF follows the predicted one (and hence the UV LF) up to z = 7.0 and then moves down at z = 7.3.

Zoom In Zoom Out Reset image size

We then calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ by measuring the luminosity decrease between the predicted and observed Lyα LFs. Previous studies have used the luminosity density to evaluate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ (e.g., Itoh et al. 2018; Konno et al. 2018; Hu et al. 2019), but this method systematically underestimates the luminosity decrease because of a fixed integration range of the LF, as explained in Appendix B. Therefore, we directly measure the luminosity decrease using another method (Figure 8). First, we set a reference cumulative number density, n_ref. Then, we look for the Lyα luminosities (${L}_{\mathrm{Ly}\alpha }^{\mathrm{obs}}(z)$ for the observed Lyα LF and ${L}_{\mathrm{Ly}\alpha }^{\mathrm{pred}}(z)$ for the predicted one) that satisfy ${\int }_{{L}_{\mathrm{Ly}\alpha }^{\mathrm{obs}}}^{\infty }{\phi }^{\mathrm{obs}}(L){dL}={\int }_{{L}_{\mathrm{Ly}\alpha }^{\mathrm{pred}}}^{\infty }{\phi }^{\mathrm{pred}}(L){dL}={n}_{\mathrm{ref}}$. We calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ as

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)}{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)}=\displaystyle \frac{{L}_{\mathrm{Ly}\alpha }^{\mathrm{obs}}(z)}{{L}_{\mathrm{Ly}\alpha }^{\mathrm{pred}}(z)}.\end{eqnarray} \tag{ 6 }$$

In this way, we calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ from the most stringent upper limit of this study at z = 7.3. We also apply the same calculation to the Lyα LFs derived by previous studies at z = 6.6 (Konno et al. 2018), 7.0 (Itoh et al. 2018; Hu et al. 2019), and 7.3 (Konno et al. 2014). ¹⁶ At each redshift, we calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ from a bright part and a faint part of the Lyα LF to examine if different parts of the LF give consistent ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ values. If not, it implies either that the actual Lyα LF does not obey Equation (5) or that ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ is not independent of Lyα luminosity. We adopt n_ref = 1 × 10⁻⁶ Mpc⁻³ for a bright part because, at this value, our z = 7.3 data can place the most stringent upper limit on the IGM transmission. For the faint part, we adopt n_ref = 1 ×10⁻⁴ Mpc⁻³, the highest value where LF measurements are available for all four redshifts. These n_ref values are common to z = 6.6, 7.0, and 7.3.

Zoom In Zoom Out Reset image size

Figure 9 and Table 2 show the values of ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ thus obtained. Error bars include uncertainties from the Lyα LFs at that redshift and z = 5.7 and the UV LF evolution. To estimate the uncertainties from the UV LF evolution, we use the first equation in Section 5.1 of Bouwens et al. (2015). From the upper limit of this study, we obtain ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)\lt 0.77$, and from the Lyα LFs of previous studies, we obtain ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(6.6,7.0)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)\simeq 1$ and ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)={0.53}_{-0.22}^{+0.18}$. The bright and faint parts give almost the same results at z = 6.6 and 7.0, ¹⁷ which is consistent with our assumption that the intrinsic Lyα LF of LAEs evolves in the same way as the UV LF of LBGs and that the effect of IGM absorption, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$, does not depend on Lyα luminosity (Equation (3)). We also plot ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ calculated in previous studies using luminosity densities. The measurements of Konno et al. (2018), Itoh et al. (2018), and Hu et al. (2019), who adopted an integration range of $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]=42.4\mbox{--}44$, are lower than our results as expected (see Appendix B).

Zoom In Zoom Out Reset image size

Table 2. Summary of ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ and x_{H I} Estimates

z	Lyα LF	L*	ϕ*	α	${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$	xH I
		(1043 erg s−1)	(10−4 Mpc−3)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
7.3	This study	⋯	⋯	⋯	<0.77	>0.28
	Konno et al. (2014) b	${0.55}_{-0.33}^{+9.45}$	${0.94}_{-0.93}^{+12.03}$	−2.5 (fixed)	${0.53}_{-0.22}^{+0.18}$ a	${0.39}_{-0.08}^{+0.12}$ a
7.0	Itoh et al. (2018)	${1.50}_{-0.31}^{+0.42}$	${0.45}_{-0.18}^{+0.26}$	−2.5 (fixed)	${0.94}_{-0.17}^{+0.12}$ a	${0.16}_{-0.16}^{+0.14}$ a
	Hu et al. (2019)	${1.20}_{-0.27}^{+0.46}$	${0.65}_{-0.33}^{+0.52}$	−2.5 (fixed)	${0.90}_{-0.14}^{+0.11}$ a	${0.20}_{-0.20}^{+0.11}$ a
6.6	Konno et al. (2018)	${1.66}_{-0.69}^{+0.30}$	${0.467}_{-0.442}^{+1.44}$	$-{2.49}_{-0.50}^{+0.50}$	${0.95}_{-0.19}^{+0.14}$ a	${0.15}_{-0.15}^{+0.16}$ a
5.7	Konno et al. (2018)	${1.64}_{-0.62}^{+2.16}$	${0.849}_{-0.771}^{+1.87}$	$-{2.56}_{-0.45}^{+0.53}$	⋯	⋯

Notes. (1) Redshift. (2) Lyα LF used to calculate ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$. (3) Characteristic luminosity of the Lyα LF. (4) Characteristic number density of the Lyα LF. (5) Faint-end slope of the Lyα LF. (6) Transmission of Lyα through the IGM obtained in Section 4.2. (7) Volume-averaged neutral hydrogen fraction in the IGM obtained in Section 4.3.1.

^a Values calculated from the faint part of the Lyα LF (n_ref = 1 × 10⁻⁴ Mpc⁻³; see Section 4.2). ^b We use the Lyα LF calculated by Itoh et al. (2018) using the data given by Konno et al. (2014).

Download table as:  ASCIITypeset image

4.3. IGM Neutral Hydrogen Fraction

4.3.1. Estimation of x_{H I}

From the ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ obtained in Section 4.2, we estimate the volume-averaged neutral hydrogen fraction in the IGM, x_{H I}, ¹⁸ in the same manner as Jung et al. (2020) assuming inhomogeneous reionization. Theoretically, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)$ is described as

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z) & = & {e}^{-{\tau }_{\mathrm{IGM}}},\\ {\tau }_{\mathrm{IGM}} & = & {\tau }_{{\rm{D}}}+{\tau }_{{\rm{H}}\,{\rm\small{II}}},\end{array}\end{eqnarray} \tag{ 7 }$$

where τ_IGM is the total optical depth of the IGM, τ_D is the optical depth of neutral patches, and τ_{H II} is the optical depth of ionized bubbles (Dijkstra 2014). We assume that τ_{H II} does not change with redshift, which leads to

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)}{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)}={e}^{-{\tau }_{{\rm{D}}}(z)},\end{eqnarray} \tag{ 8 }$$

assuming τ_D(5.7) = 0.

To obtain x_{H I}, we use the analytical approach of Dijkstra (2014) that considers inhomogeneous reionization (his Equation (30)):

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{D}({z}_{g},{\rm{\Delta }}v) & \approx & 2.3{x}_{{\rm{H}}\,{\rm\small{I}}}{\left(\displaystyle \frac{{\rm{\Delta }}{v}_{b}}{600\ \mathrm{km}\ {{\rm{s}}}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{1+{z}_{g}}{10}\right)}^{3/2},\\ {\rm{\Delta }}{v}_{b} & = & {\rm{\Delta }}v+H({z}_{g}){R}_{b}/(1+{z}_{g}),\end{array}\end{eqnarray} \tag{ 9 }$$

where z_g is the systemic redshift of a galaxy, Δv is the velocity offset of the galaxy's Lyα emission from the systemic redshift, Δv_b is the velocity offset from line resonance when the Lyα photons from the galaxy first enter a neutral patch, H(z_g ) is the Hubble constant at z_g , and R_b is the comoving distance to the surface of the neutral patch. If we adopt for Δv the typical range for z ∼ 6–8 LAEs obtained by Hashimoto et al. (2019), ${\rm{\Delta }}v={100}_{-100}^{+100}\ \mathrm{km}\ {{\rm{s}}}^{-1}$, then the unknown quantities are x_{H I} and R_b .

We then use the characteristic size of ionized bubbles, R_b , as a function of z predicted by Furlanetto & Oh (2005) with an analytic model of patchy reionization; we calculate the x_{H I}–R_b relation at z_g = 6.6, 7.0, and 7.3 by interpolating the relations at z = 6 and 9 in Figure 1 of Furlanetto & Oh (2005).

Figure 10 shows the x_{H I}–R_b relation from Dijkstra (2014) for ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)=0.77$ as an example. The blue, black, and red lines represent the calculation from Equations (8) and (9) (i.e., Dijkstra 2014) with Δv = 0, 100, and 200 km s⁻¹, respectively, which indicates that larger Δv and R_b result in higher x_{H I} because of an easier escape of Lyα photons. On the other hand, the green line is the prediction by Furlanetto & Oh (2005), which shows a larger R_b in a more ionized (lower x_{H I}) universe. We obtain x_{H I} as the intersection of these two lines and conservatively evaluate its uncertainty following Jung et al. (2020), allowing a range of Δv = 0–200 km s⁻¹ (see Figure 12 of Jung et al. 2020). In this way, we calculate x_{H I}(z) for all the ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ measurements obtained in Section 4.2.

Zoom In Zoom Out Reset image size

From our new data, we obtain x_{H I} > 0.28 at z = 7.3. From the literature Lyα LFs, we obtain ${x}_{{\rm{H}}\,{\rm\small{I}}}={0.39}_{-0.08}^{+0.12}$ at z = 7.3 and x_{H I} consistent with zero within the errors at z = 6.6 and 7.0 (Figure 11 and Table 2). Because these x_{H I} estimates are based on specific models of Lyα transmission in the IGM and the evolution of ionized bubbles, we also estimate x_{H I} using two other models and obtain consistent results (Appendix C).

Zoom In Zoom Out Reset image size

4.3.2. Possible Uncertainties in the x_{H I} Estimates

In this section, we discuss possible uncertainties in our x_{H I} estimates. First, if κ and/or ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ in Equation (3), which we assume to be constant in Section 4.2, increases at z > 5.7, the actual x_{H I} would be larger than our results. If κ or ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ increases, the emitted Lyα luminosity also increases; thus, the luminosity decrease due to neutral hydrogen needs to be greater to reproduce the observed Lyα LF, which results in higher x_{H I}. Indeed, Hayes et al. (2011) have found that ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ increases with redshift over 0 < z < 6. It is, however, not clear whether ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ and κ significantly increase from z = 5.7 to z = 7.0 (or to z = 7.3), a period shorter than 300 Myr. We also note that the Lyα EW method also adopts essentially the same assumption.

Second, bright LAEs targeted in this study may be in larger ionized bubbles than the average ones adopted in Section 4.3.1, implying that we may be underestimating x_{H I}. Figure 2 of Furlanetto & Oh (2005) shows bubble size distributions for different total masses, with the most massive regions having three times larger sizes than the average. Using the three times larger size in Section 4.3.1 will give x_{H I}(7.3) ≳ 0.4.

Finally, previous simulations (Mesinger & Furlanetto 2008; Dijkstra et al. 2011; Mason et al. 2018b; Weinberger et al. 2019) of LAEs during reionization demonstrate that ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ has a broad distribution at a given x_{H I} due to a broad range of ionized bubble sizes and a sight-line-to-sight-line scatter. This effect gives an additional uncertainty to x_{H I} estimates using ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$. However, applying this effect to Lyα LF-based x_{H I} estimates is complicated and beyond the scope of this paper.

We also note that because our analysis uses the integrated number density, any information on the shape of the Lyα LF is lost. An accurate determination of the LF shape from a deeper and larger LAE survey may place some constraints on the topology of reionization through, e.g., the dependence of bubble sizes on Lyα luminosity.

4.3.3. Comparison with Previous Studies

Figure 11 and Table 2 show the estimates of x_{H I} from our new data (red symbol in Figure 11) and the previous studies' Lyα LFs (blue symbols). Also plotted in Figure 11 are other estimates in the literature (black symbols).

First, we focus on the Lyα LF-based x_{H I} estimates. At z = 7.3, we obtain x_{H I} > 0.28 from our new data. This lower limit is consistent with our estimate from Konno et al.'s (2014) LF (the faint part of the z = 7.3 Lyα LF) and indicates that cosmic reionization is ongoing at z ∼ 7.3. On the other hand, the x_{H I}(6.6) and x_{H I}(7.0) values obtained from both the bright and faint parts of the corresponding LFs are consistent with full ionization within the errors, indicating that the universe is completing reionization around these redshifts. The estimates from this study are also consistent with those by Inoue et al. (2018), ${x}_{{\rm{H}}\,{\rm\small{I}}}(7.3)={0.5}_{-0.3}^{+0.1}$ and x_{H I}(5.7, 6.6, 7.0) < 0.4. They have predicted Lyα LFs in the fully ionized IGM not from observed UV LFs but by a physically motivated analytic model of LAEs that calculates Lyα luminosity as a function of dark halo mass. Their model reproduces observed Lyα LFs, LAE angular autocorrelation functions, and LAE fractions in LBGs at z ∼ 6–7. Very recently, Morales et al. (2021) have predicted Lyα LFs for various x_{H I} in a partially ionized universe with an analytic model of the UV LF and infer the IGM neutral fraction at z = 6.6, 7.0, and 7.3 from a comparison with observed Lyα LFs. Their x_{H I}(6.6) and x_{H I}(7.0) are consistent with our results within the errors, but their x_{H I}(7.3) is higher than ours from Konno et al.'s (2014) LF. The cause of the difference at z = 7.3 has not been fully identified, but it is partly because of the conversion from the decrease in the Lyα LF to x_{H I}.

Next, we compare these Lyα LF-based results with the other methods' results. At z ∼ 7.0, our constraints, ${x}_{{\rm{H}}\,{\rm\small{I}}}(7.0)={0.16}_{-0.16}^{+0.14}$ and ${0.20}_{-0.20}^{+0.11}$ from the Lyα LFs of Itoh et al. (2018) and Hu et al. (2019), respectively, are lower than the other results, ${x}_{{\rm{H}}\,{\rm\small{I}}}(7.0)={0.70}_{-0.23}^{+0.20}$ (Wang et al. 2020; QSO damping wing), ${x}_{{\rm{H}}\,{\rm\small{I}}}(\sim 7)\,={0.55}_{-0.13}^{+0.11}$ (Mason et al. 2018a; Whitler et al. 2020; LBG EW distribution), and x_{H I}(∼7) > 0.4 (Mesinger et al. 2015; LBG Lyα fraction) despite a large uncertainty in each estimate. If the Lyα LF-based estimates are underestimating x_{H I}, the cause could be the assumption of constant κ and ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ and/or the use of the average size of ionized bubbles as mentioned in Section 4.3.2.

Our constraint of x_{H I}(7.3) > 0.28 is broadly consistent with the other estimates at z ∼ 7.0–7.6 that span 0.2 ≲ x_{H I} ≲ 0.9, thus adding further evidence of reionization being still underway around z = 7.3. The estimate by Greig et al. (2019), ${x}_{{\rm{H}}\,{\rm\small{I}}}(7.54)={0.21}_{-0.19}^{+0.17}$, is the lowest among the all estimates including ours (although within the errors). Their source, QSO J1342, is the same as of Bañados et al. (2018) and Davies et al. (2018), who have obtained x_{H I} ∼ 0.6, but Greig et al. (2019) have analyzed only the red side of the observed Lyα spectrum to avoid complicated modeling of the near-zone transmission. Therefore, different analyses can lead to largely different results even for the same source. If the result of Greig et al. (2019) is correct, it is possible that this QSO resides in a large H ii region. Indeed, it has been suggested that QSOs inhabit highly biased overdense regions that were reionized early (e.g., Mesinger 2010; Dijkstra 2014). A similarly large difference among the four Lyα EW-based estimates over 7 ≲ z ≲ 8 may also be partly attributed to the presence or not of a highly ionized region as suggested by Jung et al. (2020), although these estimates might be detecting a real change in x_{H I} with a coarse resolution of Δz ∼ 1.

In Figure 11, we also plot semiempirical models of reionization by Finkelstein et al. (2019) and Naidu et al. (2020). Finkelstein et al. (2019) have predicted early and smooth reionization driven by faint galaxies, with a steep faint-end slope (α_UV < −2) of the UV LF and higher escape fractions of ionizing photons (${f}_{\mathrm{esc}}^{\mathrm{ion}}$) in fainter galaxies. On the other hand, Naidu et al. (2020) have predicted late and rapid reionization driven by bright galaxies, with a shallow faint-end slope (α_UV > −2). For ${f}_{\mathrm{esc}}^{\mathrm{ion}}$, Naidu et al. (2020) have examined two cases: one assuming a constant ${f}_{\mathrm{esc}}^{\mathrm{ion}}$ across all galaxies (Model I) and the other assuming ${f}_{\mathrm{esc}}^{\mathrm{ion}}$ to be dependent on the SFR surface density, Σ_SFR, of galaxies (Model II). The difference in x_{H I} between these two cases is relatively small. Our x_{H I} estimates seem to prefer Finkelstein et al.'s (2019) model to Naidu et al.'s (2020).

In summary, we provide a new constraint of x_{H I}(7.3) > 0.28 from NB-selected LAEs' Lyα LF. This is a constraint from a large (∼2 × 10⁶ Mpc³) cosmic volume with a negligibly small redshift uncertainty. If the possible underestimation of the Lyα LF method due, for example, to an increase in κ or ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ is true, then the actual x_{H I} will be even higher.

We have derived a new constraint on the Lyα LF of z = 7.3 LAEs based on a large-area NB imaging survey with Subaru/Hyper Suprime-Cam whose effective survey volume is ∼2 ×10⁶ Mpc³. Using this constraint, we have calculated the Lyα transmission in the IGM, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$, and then the volume-averaged neutral hydrogen fraction in the IGM, x_{H I}, at z = 7.3. In the calculation of ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$, we have applied a new method that directly measures the luminosity decrease between an observed LF and a predicted LF (i.e., LF for the fully ionized IGM predicted by the evolution of the UV LF). We have also applied this method to previous studies' Lyα LFs at z = 6.6, 7.0, and 7.3. Our main results are summarized below.

• 1.  
  We have detected no z = 7.3 LAEs in either the COSMOS or SXDS field, which results in a decrease in the bright part of the Lyα LF from z = 7.0 (Itoh et al. 2018; Hu et al. 2019) to z = 7.3 (Figure 6).
• 2.  
  To estimate x_{H I}, we have predicted the Lyα LF in the case of the fully ionized IGM on the assumption that the intrinsic Lyα LF evolves in the same way as the observed UV LF. We have found that the observed Lyα LF follows the predicted one (and hence the UV LF) up to z = 7.0 and then moves down at z = 7.3 (Figure 7).
• 3.  
  We have estimated ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$ in a new method that directly measures the luminosity decrease between an observed Lyα LF and a predicted one. We have obtained ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)\lt 0.77$ from our new data, and ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(6.6,7.0)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)\simeq 1$ and ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(7.3)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)={0.53}_{-0.22}^{+0.18}$ from the previous studies' Lyα LFs (Konno et al. 2014; Itoh et al. 2018; Konno et al. 2018; Hu et al. 2019; Figure 9 and Table 2). Bright and faint parts of the Lyα LF give almost the same ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$.
• 4.  
  Using the obtained ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}$, we have estimated x_{H I} in the same manner as Jung et al. (2020). The constraint of x_{H I}(7.3) > 0.28 estimated from our new data is broadly consistent with the other estimates in the literature and indicates that cosmic reionization is still ongoing at z ∼ 7.3. On the other hand, the x_{H I}(6.6) and x_{H I}(7.0) estimated from the previous studies' Lyα LFs are consistent with full ionization but are lower than the other estimates at the same redshift (Mesinger et al. 2015; Mason et al. 2018a; Wang et al. 2020; Whitler et al. 2020; Figure 11 and Table 2). If this implies underestimation of our calculation due, for example, to an increase with redshift in the Lyα escape fraction of galaxies or the Lyα production rate per UV luminosity, then the lower limit to x_{H I}(7.3) will also become higher than 0.28.

We thank the anonymous referee for the constructive comments that greatly improved the manuscript. We thank Rikako Ishimoto, Kei Ito, Ryohei Itoh, Ryota Kakuma, Shotaro Kikuchihara, Haruka Kusakabe, and Yongming Liang for useful comments and discussions. K.S. is supported by the Toray Science Foundation. T.H. is supported by Leading Initiative for Excellent Young Researchers, MEXT, Japan (HJH02007) and JSPS KAKENHI grant number 20K22358. R.M. acknowledges a Japan Society for the Promotion of Science (JSPS) Fellowship at Japan and JSPS KAKENHI grant No. JP18J40088. This work is supported by JSPS KAKENHI grant No. 17H01114 (A.K.I., S.Y., K.S., and M.O.).

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

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

This paper is 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 (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by 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, the 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 grant No. AST-1238877, the University of Maryland, Eötvös Loránd University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

In Figure 12, we present example images of spurious sources removed in our visual inspection (Section 3.2). For comparison, we also show example images of pseudo-LAEs used in the calculation of completeness (Section 3.3). The two spurious sources in the top row of this figure are either a cosmic ray or a CCD artifact. Both have a well-outlined and highly concentrated light distribution, with their total luminosity being contributed by only a small number (≲10) of pixels. The two spurious sources in the bottom row are diffuse, elongated structures. Each of them is due to a different, very bright star and is located outside the mask for the star.

Zoom In Zoom Out Reset image size

To obtain the relative transmission at a certain redshift, ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$, we assume that observed Lyα luminosities are uniformly decreased in proportion to the IGM transmission (Equation (3)). We then estimate this luminosity decrease by directly comparing the observed and predicted Lyα luminosities at a fixed cumulative number density (Equation (6)). Previous studies have, however, estimated ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ from a decrease in the Lyα luminosity density as

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)}{{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)}=\displaystyle \frac{{\rho }_{\mathrm{Ly}\alpha }(z)/{\rho }_{\mathrm{Ly}\alpha }(5.7)}{{\rho }_{\mathrm{UV}}(z)/{\rho }_{\mathrm{UV}}(5.7)},\end{eqnarray} \tag{ B1 }$$

where ρ_Lyα and ρ_UV are the Lyα and UV luminosity densities, respectively, calculated by integrating the corresponding LFs. Comparison with Equation (6) finds that ρ_Lyα (z) is used in place of ${L}_{\mathrm{Ly}\alpha }^{\mathrm{obs}}(z)$ in Equation (6), and ρ_Lyα (5.7)ρ_UV(z)/ρ_UV(5.7) in ${L}_{\mathrm{Ly}\alpha }^{\mathrm{pred}}(z)$. However, Equation (B1) is correct only when all four corresponding LFs are integrated down to zero luminosity. Specifically, Equation (B1) overestimates the Lyα luminosity decrease and hence underestimates ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ if the integration range of the two Lyα LFs is limited, as has been done by most previous studies: e.g., $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]$ =42.4–4 has been adopted by Konno et al. (2018), Itoh et al. (2018), and Hu et al. (2019).

To see why Equation (B1) overestimates the luminosity decrease with limited integration ranges of the Lyα LFs, let us assume a simple case that Lyα luminosities are uniformly decreased due to IGM absorption (the same assumption as in this study), and the UV LF does not evolve with redshift. In this case, the right-hand side of Equation (B1) is reduced to ρ_Lyα (z)/ρ_Lyα (5.7). If the Lyα LF at z = 5.7 has the characteristic luminosity L*(5.7), then that at a redshift before completion of reionization will have ${L}^{* }(z)={L}^{* }(5.7)\times {T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$, with the remaining two Schechter parameters (Schechter 1976), ϕ* and α, being the same as those of the z = 5.7 LF because we have assumed a uniform luminosity decrease due to IGM absorption. The question now is whether the equation ρ_Lyα (z)/ρ_Lyα (5.7) = L*(z)/L*(5.7) is correct. Using the Schechter parameters above, ρ_Lyα (z)/ρ_Lyα (5.7) is written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rho }_{\mathrm{Ly}\alpha }(z)}{{\rho }_{\mathrm{Ly}\alpha }(5.7)} & = & \displaystyle \frac{{\displaystyle \int }_{{L}_{\mathrm{lim}}}^{\infty }L\phi ({L}^{* }(z),{\phi }^{* },\alpha ;L){dL}}{{\displaystyle \int }_{{L}_{\mathrm{lim}}}^{\infty }L\phi ({L}^{* }(5.7),{\phi }^{* },\alpha ;L){dL}}\\ & = & \displaystyle \frac{{\displaystyle \int }_{{L}_{\mathrm{lim}}}^{\infty }L{\phi }^{* }{\left(\tfrac{L}{{L}^{* }(z)}\right)}^{\alpha }\exp \left(-\tfrac{L}{{L}^{* }(z)}\right)\tfrac{{dL}}{{L}^{* }(z)}}{{\displaystyle \int }_{{L}_{\mathrm{lim}}}^{\infty }L{\phi }^{* }{\left(\tfrac{L}{{L}^{* }(5.7)}\right)}^{\alpha }\exp \left(-\tfrac{L}{{L}^{* }(5.7)}\right)\tfrac{{dL}}{{L}^{* }(5.7)}}\\ & = & \displaystyle \frac{{L}^{* }(z)}{{L}^{* }(5.7)}\times \displaystyle \frac{{\displaystyle \int }_{{L}_{\mathrm{lim}}/{L}^{* }(z)}^{\infty }{x}^{\alpha +1}\exp (-x){dx}}{{\displaystyle \int }_{{L}_{\mathrm{lim}}/{L}^{* }(5.7)}^{\infty }{x}^{\alpha +1}\exp (-x){dx}},\end{array}\end{eqnarray} \tag{ B2 }$$

where L_lim is the integration limit and x ≡ L/L^*(z). ¹⁹ The right-hand side of the last line of this equation is lower than ${L}^{* }(z)/{L}^{* }(5.7)(={T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7))$, because the integral at the numerator is smaller than that at the denominator owing to the narrower integration range of x (because of ${L}_{\mathrm{lim}}/{L}^{* }(z)\,\gt {L}_{\mathrm{lim}}/{L}^{* }(5.7)$).

As an example, let us take $\mathrm{log}{L}^{* }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]=43.2$ and α = −2.56 at z = 5.7 (the values obtained by Konno et al. 2018) and assume a 50% luminosity decrease due to IGM absorption, i.e., ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)={L}^{* }(z)/{L}^{* }(5.7)=0.5$. In this case, Equation (B2) with an integration range of $\mathrm{log}{L}_{\mathrm{Ly}\alpha }\ [\mathrm{erg}\ {{\rm{s}}}^{-1}]=$ 42.4–44 gives ρ_Lyα (z)/ρ_Lyα (5.7) = 0.23, i.e., 77% decrease.

On the other hand, Ouchi et al.'s (2010) and Konno et al.'s (2014) results include no systematic bias because they integrated the LF down to zero luminosity. However, their strategy instead leads to extremely large uncertainty in the luminosity density due to a large extrapolation of the LF from the observed luminosity range.

Estimating x_{H I} using the Lyα luminosity of galaxies requires a theoretical model that relates observed Lyα luminosities, or Lyα LFs, with x_{H I}. To mitigate model dependence, we also estimate x_{H I} using methods other than Jung et al.'s (2020), as in previous studies. First, using Lyα LFs for several x_{H I} values simulated by Inoue et al. (2018) (their Figure 18; see also Section 4.3.3), we obtain x_{H I}(7.3) > 0.2 from our new data. Second, we use the analytic model of Santos (2004) that calculates observed Lyα emission from an isolated galaxy at z = 6.5 by varying many parameters associated with the galaxy and the IGM around it. By comparing our ${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$ estimates with their Figure 25, we obtain x_{H I} for a galactic wind with an Lyα velocity offset of 0 and 360 km s⁻¹, as shown in Table 3. These values are consistent with our results from the Jung et al. (2020) method within the errors.

Table 3.  x_{H I} Estimates from Santos (2004)

z	Lyα LF	${T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(z)/{T}_{\mathrm{Ly}\alpha }^{\mathrm{IGM}}(5.7)$	xH I
			v = 0 km s−1	v = 360 km s−1
7.3	this study	<0.77	>0.2	>0.2
	Konno et al. (2014)	${0.53}_{-0.22}^{+0.18}$	${0.4}_{-0.2}^{+0.2}$	${0.5}_{-0.2}^{+0.3}$
7.0	Itoh et al. (2018)	${0.94}_{-0.17}^{+0.12}$	0+0.2	${0.1}_{-0.1}^{+0.1}$
	Hu et al. (2019)	${0.90}_{-0.14}^{+0.11}$	${0.1}_{-0.1}^{+0.1}$	${0.1}_{-0.1}^{+0.1}$
6.6	Konno et al. (2018)	${0.95}_{-0.19}^{+0.14}$	0+0.2	${0.1}_{-0.1}^{+0.1}$

Download table as:  ASCIITypeset image