SILVERRUSH. IX. Lyα Intensity Mapping with Star-forming Galaxies at z = 5.7 and 6.6: A Possible Detection of Extended Lyα Emission at ≳100 Comoving Kiloparsecs around and beyond the Virial-radius Scale of Galaxy Dark Matter Halos

For source catalogs of our intensity mapping analysis, we use the LAE catalogs obtained in the program of "Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC" (SILVERRUSH; Ouchi et al. 2018; Shibuya et al. 2018). The SILVERRUSH LAE catalogs are derived from the HSC-SSP survey data whose first data release is presented in Aihara et al. (2018b). This is the largest z ≳ 5 LAE catalog to date. The details of the SILVERRUSH catalogs are listed in Table 1.

Table 1. Summary of Our LAE Samples

Field	Area	${m}_{\mathrm{lim}}$	NLAE	$\mathrm{log}\left({L}_{\mathrm{Ly}\alpha }/\,[\mathrm{erg}\,{{\rm{s}}}^{-1}]\right)$
	(deg2)	(mag)
(1)	(2)	(3)	(4)	(5)
NB816 (z ∼ 5.7)
UD-COSMOS	1.97	25.7	201	⋯
UD-SXDS	1.93	25.5	224	⋯
Total	3.9	⋯	425	42.0–43.8
NB921 (z ∼ 6.6)
UD-COSMOS	2.05	25.6	338	⋯
UD-SXDS	2.02	25.5	58	⋯
Total	4.07	⋯	396	42.3–44.0

Note. (1) Field. (2) Survey Area. (3) Limiting magnitude of the NB image defined by the 5σ detection level in a 15 diameter circular aperture. (4) Number of the LAEs. (5) The range of LAE Lyα luminosity.

Download table as:  ASCIITypeset image

The number of z = 6.6 LAEs in the UD-SXDS field appears to be smaller than that of the UD-COSMOS field. Shibuya et al. (2018) attribute this difference to the seeing size of the NB921 images of the UD-SXDS field (∼07), which is worse than that of the UD-COSMOS field (∼06). Figure 5 of Shibuya et al. (2018) shows that the surface number density of LAEs in the UD-COSMOS field is comparable to those identified in Subaru/Suprime-Cam NB surveys. Though the SILVERRUSH program exploits the HSC-SSP data taken in the D and UD layers, we use only those taken in the UD layer that provides the highest quality data in the HSC SSP survey. Figure 2 displays the sky distribution of the LAEs in the UD layer.

Zoom In Zoom Out Reset image size

For an intensity map of the Lyα emission at a redshift that is the same as those of LAEs at z = 5.7 (6.6), we use NB816 (NB921) images in the HSC-SSP survey S18A data release (Aihara et al. 2019). The NB images in the S18A data release are ∼1 mag deeper than those in the S16A data release. The images in the S18A data release are reduced with the HSC pipeline v6.7 (Bosch et al. 2018) that implements an improved sky-background subtraction approach. The HSC pipeline v6.7 jointly models and subtracts scattered lights from bright objects and instrumental features crossing all CCDs. In addition, the sky-background emission are estimated and subtracted in the mosaic with a grid size of 6000 pixels, corresponding to about 17'. This new sky-subtraction method reduces overfittings and oversubtractions of the sky background made by small-scale fluctuations. Although the latest S18A data was released very recently (2018 June), we have not selected LAEs from this data.

The NB images contain not only Lyα emission at z = 5.7 or 6.6 but also continuum and emission lines that are emitted from low-z sources. The effects from these contaminants can nevertheless be removed by taking a spatial cross-correlation with LAEs because these sources are randomly located on the sky with respect to the LAEs. The low-z sources only add noise on the cross-correlation between the LAEs and Lyα emission.

Note that we need to mask bright sources on the images to improve the signal-to-noise ratio (S/N) of the diffuse extended emission. The HSC pipeline sets some flags to each pixel. We do not use pixels with flags of BRIGHT_OBJECT or DETECTED. The BRIGHT_OBJECT flag is given to pixels where nearby very bright objects would affect the background subtraction or detection. The DETECTED flag is given to the pixels in which sources are detected at the S/N > 5 level. With these flags, the pixels in which LAEs are detected are also masked out. We focus on extended Lyα emission at ≥15 separated from the LAEs. This is about two times larger than the size of LAEs on the images that are marginally resolved.

We derive the angular cross-correlation function between the LAEs and the Lyα intensity ${\xi }_{\mathrm{IM}}(r)$, taking a mean over all LAE-pixel pairs that are separated by r within a certain bin:

$$\begin{eqnarray}&&{\xi }_{\mathrm{IM}}(r)=\displaystyle \frac{1}{N}\sum _{i}{\mu }_{r;i},\end{eqnarray} \tag{ 1 }$$

where μ_r;i is the Lyα intensity in the images at pixel i for the bin r. N is the number of all LAE-pixel pairs separated by r within a certain bin. The obtained ${\xi }_{\mathrm{IM}}(r)$ is not the dimensionless cross-correlation function, but the cross-correlation function with the unit of intensity. Note that this 2D intensity mapping analysis method is basically the same as image stacking. We set six bins logarithmically spaced between 15 and 40''. In our calculation, image pixels can contribute more than once to the cross-correlation functions. However, the contribution of such pixels is unlikely to be significant, since the minimum distance between our LAEs is about 20'' and most of them are separated by more than 100''. Statistical uncertainties on the measurements are computed by a bootstrap resampling method. We randomly select LAEs with numbers equal to that of our LAEs, allowing a duplication. We create 1000 LAE samples by random selection, and derive a 1σ standard deviation that is referred to as a 1σ error.

There are a number of systematic uncertainties that may produce spurious extended sources in the measurement of diffuse emissions. For example, a largely extended point-spread function (PSF) and flat-fielding/sky-subtraction systematics may affect the extended profile. We need to carefully evaluate total uncertainties from these systematics. To test the systematics in the cross-correlation between the Lyα intensity map and the LAEs, we derive a cross-correlation function between the Lyα intensity map and real foreground objects residing at neither z = 5.7 nor 6.6. We refer to these objects as non-LAEs. The non-LAEs are not correlated with Lyα emission at neither z = 5.7 nor 6.6.

For a catalog of non-LAEs, we use the g-band dropout galaxy catalog taken from the "Great Optically Luminous Dropout Research Using Subaru HSC" (GOLDRUSH; Harikane et al. 2018a; Ono et al. 2018) program. These non-LAEs (g-dropout galaxies) are selected mainly by probing the redshifted Lyman break feature and relatively blue UV continuum slope based on the broadband color selection criteria of g − r > 1.0, r − i < 1.0, and g − r > 1.5(r − i) + 0.8. Based on their selection completeness simulations, the expected redshifts of the non-LAEs are z ≃ 3.8 ± 0.5, which is significantly lower than those of our LAEs, indicating that their Lyα emission is not captured in our NB images. Thanks to the wide redshift coverage, the total number of the non-LAEs is about 10,000 in each field. Because the systematics like the largely extended PSF should depend on the magnitude of the sources, we randomly select the non-LAEs with the same NB magnitude distribution as that of our LAEs. Since we should reduce the statistical errors of this analysis originated from the non-LAE samples as much as possible, we construct a sample of non-LAEs whose number is 3–13 times larger than that of the LAEs in each field. We measure the cross-correlation with non-LAEs in the same manner as the one with LAEs, and evaluate the systematic uncertainties.

In Figure 3, we show the cross-correlation between the LAEs (non-LAEs) and the Lyα intensity map. The cross-correlations with the non-LAEs have positive values even at a large scale. This is probably because the images in the S18A data release have small residuals that are left in the sky-background subtraction, while the various techniques are applied to remove the residuals. Comparing the cross-correlations of the LAEs and the non-LAEs, we estimate the amount of the small residuals. Figure 3 indicates that the amplitude of the cross-correlation function of the LAEs is higher than those of non-LAEs at each redshift in each field, although the z = 5.7 LAEs in UD-SXDS show no clear excess especially at a large scale of >200 ckpc. Because the cross-correlation functions of the LAEs mostly exceed those of the non-LAEs, we confirm that the systematic uncertainties alone do not explain the cross-correlation functions of the LAEs, but real signals of spatially extended Lyα emission.

Zoom In Zoom Out Reset image size

To evaluate the spatially extended Lyα emission quantitatively, we subtract the cross-correlation function of the non-LAEs from that of the LAEs as shown in the bottom panels of Figure 3. The calculated values are presented in Table 2. We obtain weighted mean values of the results in the UD-COSMOS and UD-SXDS fields that are calculated with the weights defined by the inverse values of the 1σ errors. The uncertainties of the measurements are estimated on the basis of the error propagation. Our results show spatially extended emission at the scale of ≳100 ckpc.

Table 2. Lyα Flux Cross-correlation Functions for the z = 5.7 and z = 6.6 LAEs

θ	ξ(r)COSMOS	ξ(r)SXDS	ξ(r)	Ξ(r)
(arcsec)	(10−21 erg s−1 cm−2 Å−1 arcsec−2)	(10−21 erg s−1 cm−2 Å−1 arcsec−2)	(10−21 erg s−1 cm−2 Å−1 arcsec−2)	(10−20 erg s−1 cm−2 arcsec−2)
(1)	(2)	(3)	(4)	(5)
z = 5.7
2.05	${1.225}_{-0.529}^{+0.535}$	${0.382}_{-0.643}^{+0.648}$	0.883 ± 0.409	9.98 ± 4.62
3.54	${0.685}_{-0.515}^{+0.509}$	${0.550}_{-0.593}^{+0.591}$	0.627 ± 0.389	7.09 ± 4.39
6.13	${0.662}_{-0.483}^{+0.474}$	${0.140}_{-0.551}^{+0.556}$	0.438 ± 0.363	4.95 ± 4.10
10.6	${0.764}_{-0.462}^{+0.461}$	$-{0.133}_{-0.515}^{+0.515}$	0.364 ± 0.344	4.11 ± 3.88
18.3	${0.368}_{-0.403}^{+0.396}$	$-{0.247}_{-0.483}^{+0.489}$	0.120 ± 0.309	1.36 ± 3.50
31.7	${0.207}_{-0.346}^{+0.335}$	$-{0.221}_{-0.459}^{+0.465}$	0.056 ± 0.276	0.64 ± 3.12
z = 6.6
2.05	${0.600}_{-0.389}^{+0.378}$	${2.326}_{-1.309}^{+1.342}$	0.734 ± 0.373	9.90 ± 5.04
3.54	${0.512}_{-0.369}^{+0.356}$	${1.249}_{-1.382}^{+1.371}$	0.559 ± 0.357	7.55 ± 4.81
6.13	${0.197}_{-0.348}^{+0.338}$	${1.448}_{-1.374}^{+1.369}$	0.270 ± 0.337	3.65 ± 4.55
10.6	${0.150}_{-0.335}^{+0.334}$	${1.746}_{-1.379}^{+1.375}$	0.239 ± 0.326	3.22 ± 4.40
18.3	$-{0.003}_{-0.323}^{+0.320}$	${1.776}_{-1.351}^{+1.348}$	0.092 ± 0.314	1.24 ± 4.24

Note. (1) Projected distances from the LAEs. (2) Cross-correlation between the LAEs and the Lyα intensity for the COSMOS field. (3) Same as (2) but for the SXDS field. (4) Same as (2) but for the average values of the two fields. (5) Lyα flux cross-correlation functions obtained from Equation (2).

Download table as:  ASCIITypeset image

To examine the possibility that this spatially extended emission in the weighted mean of the cross-correlation functions is caused by other systematics that are not corrected in the analysis above, we also investigate the cross-correlation between the LAEs and the intensity map taken with another filter. Here we use g-band images taken from the HSC-SSP survey S18A data release (Section 2.2). The 5σ limiting magnitude of the g-band images is ∼27 mag. Because the g-band filter covers a wavelength range shorter than the observed-frame wavelength of Lyα emission at z = 5.7 and 6.6, the LAEs in our samples should not be detected in the g-band images. We thus expect that there are no signals in the cross-correlation between the LAEs and the g-band images. One can test the reliability of our results with the systematics subtraction, investigating departures from zero in the cross-correlation between the LAEs and the g-band images. Performing this test, we derive the cross-correlation function between the LAEs and the g-band images in the UD-COSMOS and UD-SXDS fields, respectively, that are corrected for the systematics in the same manner as our results. In this test, the non-LAE sample consists of g-dropouts whose g-band fluxes are mostly negligibly small. We confirm that the cross-correlation function is consistent with zero within the errors, and that there is no significant systematics mimicking the spatially extended Lyα emission in our results of Figure 3.

Note that the cross-correlations of the z = 5.7 (z = 6.6) LAEs after subtracting the systematics in COSMOS (SXDS) appear to be larger than that in the other field with large uncertainties. Although we carefully control the systematics by using the cross-correlations of the non-LAEs, this difference might be induced by some unknown uncertainty that is not corrected in our careful analysis. Conceivably this difference might be caused by an environmental effect, although in that case the dependency of the cross-correlation on the Galaxy overdensity appears to be opposite to that reported in previous work (Matsuda et al. 2012). A more detailed examination of this issue is beyond the scope of this paper; we would like to derive cross-correlations of LAEs that are detected with other NBs whose central wavelengths are shorter than ours, to check whether or not similar differences are seen, which may be helpful to address this issue, in a forthcoming paper.

In Figure 4, we show the weighted mean values of the cross-correlation functions at z = 5.7 and 6.6. Here, we convert ξ(r) [erg s⁻¹ cm⁻² Å⁻¹ arcsec⁻²] to the Lyα flux cross-correlation function Ξ(r) [erg s⁻¹ cm⁻² arcsec⁻²] by multiplying the FWHM of NB filter (FWHM_NB).

$$\begin{eqnarray}&&{\rm{\Xi }}(r)={\xi }_{\mathrm{IM}}(r)\times {\mathrm{FWHM}}_{\mathrm{NB}}.\end{eqnarray} \tag{ 2 }$$

The obtained results are presented in Table 2. For comparison, in Figure 4, we also present the Lyα radial profile of relatively luminous LAEs at similar redshifts of z = 5–6 with Lyα luminosities of >10⁴² erg s⁻¹ obtained by MUSE observations (Extended Data Figure 5 in Wisotzki et al. 2018). As shown in Figure 4, the MUSE observations identify Lyα emission at the small scale of <150ckpc beyond the errors. In this scale, the Lyα intensity radial profiles of our results are consistent with those of the MUSE observations both at z = 5.7 and 6.6. We also check previous MUSE results for Lyα radial profiles around individual high-z LAEs reported by Leclercq et al. (2017). In their sample, the #547 LAE has a comparable Lyα luminosity (10^42.77erg s⁻¹) and is located at a similar redshift (z = 5.98) to our samples. We find that our results are also consistent with the Lyα radial profile of this individual LAE at the small scale of <150 ckpc.

Zoom In Zoom Out Reset image size

At the larger scale of ≳100 ckpc, our results indicate the possible existence of very faint spatially extended Lyα emission more largely beyond the errors, which is not probed with the MUSE results. We estimate the detection confidence levels of the spatially extended Lyα emission at the scale of ≳100 ckpc (≃80–400 ckpc) to be the 2.6σ and 2.2σ levels at z = 5.7 and 6.6, respectively, based on Fisher's method. ¹³ Figure 4 also indicates a hint of the extended Lyα emission up to ∼1000 ckpc. Based on the halo occupation distribution models, LAEs in our sample are hosted by the DMHs with average masses of $\mathrm{log}(\langle {M}_{{\rm{h}}}\rangle /{M}_{\odot })={10.8}_{-0.5}^{+0.3}$ (${11.1}_{-0.4}^{+0.2}$) at z = 5.7 (6.6) with a Lyα duty cycle of 1% or less (Ouchi et al. 2018). The DMH virial radius of a galaxy whose halo mass is $\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })=11$ is r_vir ∼ 150 ckpc (Mo & White 2002). By the comparison with the DMH virial radius of r_vir ∼ 150 ckpc, our study has shown the marginal detection of the spatially extended Lyα emission around and beyond the DMH virial radius at the ≃3σ (≃2σ) level at z = 5.7 (6.6). ¹⁴

By the cross-correlation analysis, we tentatively identify the Lyα emission spatially extending beyond the scale of r_vir. Theoretical studies predict that H i gas in the CGM and IGM resonantly scatters Lyα photons that escape from the interstellar medium (ISM) of star-forming galaxies. Zheng et al. (2010) perform a radiation-hydrodynamic reionization simulation in a cosmological volume, and present a physical model of Lyα emission observed around LAEs. The radiative transfer simulation is conducted to explain the observed properties of LAEs at z ∼ 5.7 in the Subaru/XMM-Newton Deep Survey (SXDS; Ouchi et al. 2008). This simulation includes physical processes of Lyα photons scattered by the CGM and IGM. On the basis of the simulation, Zheng et al. (2011) present Lyα radial profiles of stacked images of model LAEs in the simulation. The model LAEs are residing in halos with masses of M_h ∼ 10¹¹ M_⊙. Figure 5 compares Lyα radial profiles of the simulation results and our observational results for the LAEs at z = 5.7. For comparison, we use the simulation results of "the bright half of the sources" presented in Figure 6 of Zheng et al. (2011), because these sources have Lyα luminosities comparable to those of our LAEs identified in the HSC observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zheng et al. (2011) present that the Lyα radial profile of the simulation consists of two components whose Lyα sources are star-forming regions of the LAE (called "central LAE") and clustered sources around the LAE (called "clustering"), i.e., neighboring galaxies including satellites. Figure 5 also shows the total Lyα radial profile, a sum of the two-component Lyα radial profiles. Figure 5 indicates that our observational results agree neither with the total profile nor the clustered-source profile. This is probably because we mask out the pixels in which LAEs are detected at S/N > 5 (Section 2.2). In other words, in our measurements, the contribution from satellite LAEs down to our Lyα luminosity limit of about 10⁴³ erg s⁻¹ is subtracted. Note that, based on the previous results of very deep LAE searches, LAEs whose Lyα luminosities are well below our detection limits do exist (Drake et al. 2017), and such faint LAEs would contribute to some extent to our obtained Lyα cross-correlation function. However, Figure 5 suggests that the Lyα radial profile shapes are different between the clustered-source "clustering" and our observational results. In Figure 5, the "clustering" profile decreases only by a factor less than 2 at 150–1000 ckpc, while our observational result profile changes by an order of magnitude beyond the errors. This profile shape difference suggests that such faint LAEs are unlikely to be dominant contributors to the Lyα cross-correlation. In Figure 5, the Lyα radial profile of our observational results with the error bars prefers the model of "central LAE" in the simulation. However, the Lyα radial profile of our observational results is higher than that of the "central LAE" component, albeit with the large observational uncertainties. The Lyα radial profile of the observations exceeding the "central LAE" component may be produced by other physical processes that are not included in the simulation of Zheng et al. (2010), such as a cold-gas stream and galactic outflow that generate Lyα photons by collisional excitation processes.

To understand physical origins of the spatially extended Lyα emission at the large scale, numerical simulations with various possible physical processes are needed (Sadoun et al. 2019). We also need Lyα intensity mapping measurements with small uncertainties to distinguish the various possible physical models. More Lyα intensity mapping studies should be conducted with the existing instruments including Subaru/HSC. Moreover, we expect to obtain conclusive Lyα intensity mapping results by the on-going and forthcoming projects such as Hobby–Eberly Telescope Dark Energy Experiment (HETDEX; Hill et al. 2008), Wide-Field InfrarRed Survey Telescope (WFIRST; Spergel et al. 2015), and Spectro-Photometer for the History of the Universe Epoch of Reionization and Ices Explorer (SPHEREx; Doré et al. 2014).

Many observational studies have claimed the increase of the IGM neutral hydrogen fraction ${x}_{{\rm\small{}}{{\rm{H}}}_{I}}$ with the decrease of the Lyα luminosity function from z ∼ 6 to 7 and beyond that can be explained by the increase of the Lyα damping wing absorption given by the neutral IGM (e.g., Ouchi et al. 2010; Goto et al. 2011; Kashikawa et al. 2011). Jeeson-Daniel et al. (2012) predict that the Lyα radial profile of LAEs flattens toward the early epoch of cosmic reionization with a high ${x}_{{\rm\small{}}{{\rm{H}}}_{I}}$, due to the IGM neutral hydrogen scattering Lyα photons. Because we have the observational data of the Lyα radial profiles at the epoch of reionization (z = 6.6) and the post reionization epoch (z = 5.7), we examine whether the flattening is found in our observational data. It should be noted that the Lyα luminosity ranges of our LAEs at z = 5.7 and 6.6 are comparable, $\mathrm{log}({L}_{\mathrm{Ly}\alpha }/[\mathrm{erg}\ {{\rm{s}}}^{-1}])$ = 42.0–43.8 and 42.3–44.0, respectively (Konno et al. 2018). The difference between the median values of the Lyα luminosity of the LAEs at z = 5.7 and 6.6 is small (within a factor of 3). Figure 6 compares our observational results at z = 5.7 and 6.6. Here we apply corrections for the surface brightness dimming effect in our z = 6.6 results, and carry out the comparison at z = 5.7. Figure 6 indicates that our observational results at z = 5.7 and 6.6 are similar, and that there is no signature of the flattening of the Lyα radial profile toward high-z beyond the errors. Because the error bars, especially at a large scale, of our observational results are probably too large to identify the Lyα radial profile flattening caused by the cosmic reionization, one should investigate a Lyα radial profile slope with the large observational data that will be taken by on-going and forthcoming programs such as HETDEX, WFIRST, and SPHEREx discussed in Section 4.1.