MMT spectroscopy of Lyman-alpha at z ≃ 7: evidence for accelerated reionization around massive galaxies

ABSTRACT

1 INTRODUCTION

Within the past decade, a number of concerted efforts have aimed at better understanding how and when the process of hydrogen reionization occurred. Such information provides insight into not only the growth of structure in the very early Universe, but also the nature of ionizing sources responsible for reionization (e.g. Loeb & Barkana 2001; Fan, Carilli & Keating 2006; Bouwens et al. 2015b; Robertson et al. 2015; Finkelstein 2016; Stark 2016; Dayal & Ferrara 2018; Naidu et al. 2020). Observationally, there has been great progress in identifying star-forming galaxies at z ∼ 6–10 (e.g. Ellis et al. 2013; McLure et al. 2013; Bowler et al. 2014, 2020; Atek et al. 2015; Bouwens et al. 2015a, 2019; Finkelstein et al. 2015; Livermore, Finkelstein & Lotz 2017; Oesch et al. 2018; Ono et al. 2018; Stefanon et al. 2017, 2019), revealing a population capable of generating copious amounts of ionizing photons in the early Universe.

The impact of these early ionizing agents on the IGM can be tracked independently using measurements of the timeline of the reionization process. Prominent detections of the Ly α and Ly β forests in quasar spectra at z ≲ 6 provide a model-independent constraint suggesting that reionization was largely complete by z = 5.9 (neutral fraction x_{H i} ≲ 10 per cent; McGreer, Mesinger & D’Odorico 2015). However, the presence of a long neutral patch identified at z = 5.5 (Becker et al. 2015) suggests that reionization may have ended as late as z ≃ 5.2–5.3 (Kulkarni et al. 2019; Keating et al. 2020; Nasir & D’Aloisio 2020). At yet higher redshifts, strong damping wing signatures in the earliest known quasars (z ≃ 7.0–7.5) provide evidence for significant neutral hydrogen fractions (x_H i ∼ 40–60 per cent; Greig et al. 2017; Bañados et al. 2018; Davies et al. 2018; Wang et al. 2020). Measurements of the Thomson scattering optical depth from the CMB are consistent with this reionization timeline, suggesting a reionization midpoint of z ∼ 7.7 (Planck Collaboration XLVII 2016; Planck Collaboration VI 2020).

High-redshift (z ≳ 6) galaxies have also long enabled complementary constraints on the timeline of reionization owing to the resonant nature of the Ly α emission line (e.g. Miralda-Escudé 1998; Malhotra & Rhoads 2004; Santos 2004; Mesinger & Furlanetto 2008; Dijkstra 2014). A number of studies have uncovered a marked decline in the fraction of star-forming galaxies showing strong Ly α (rest-frame equivalent width EW > 25 Å) at z > 6, consistent with expectations of an increasingly neutral IGM (e.g. Fontana et al. 2010; Stark et al. 2010; Vanzella et al. 2011; Ono et al. 2012; Treu et al. 2013; Caruana et al. 2014; Pentericci et al. 2014, 2018; Schenker et al. 2014; De Barros et al. 2017; Jung et al. 2017; Hoag et al. 2019; Mason et al. 2019). In a similar manner, the volume density of narrow-band selected Ly α emitters has been found to decrease rapidly from z ∼ 6 to z ∼ 7 (e.g. Malhotra & Rhoads 2004; Hu et al. 2010; Ouchi et al. 2010; Kashikawa et al. 2011; Konno et al. 2014; Ota et al. 2017; Zheng et al. 2017). Both these downturns suggest a highly neutral Universe (x_H i > 0.4) at z ≃ 7 (e.g. Ota et al. 2017; Zheng et al. 2017; Mason et al. 2018a; Whitler et al. 2020), consistent with inferences from quasar spectra.

In the last 5 yr, attention has begun to focus on the Ly α properties of the most massive galaxies known in the reionization era. Recent spectroscopic observations have revealed a 100 per cent Ly α detection rate among four of the brightest¹ (3–4 L|$_{\mathrm{UV}}^{\ast }$|⁠) and most massive z = 7–9 galaxies selected over the CANDELS fields (Oesch et al. 2015; Zitrin et al. 2015; Roberts-Borsani et al. 2016, hereafter the RB16 sample; Stark et al. 2017), in marked contrast to the ≲10–20 per cent Ly α detection rate among typical (<L|$_{\mathrm{UV}}^{\ast }$|⁠) galaxies at z ≳ 7 (e.g. Ono et al. 2012; Treu et al. 2013; Pentericci et al. 2014, 2018; Schenker et al. 2014). It is not yet clear why the Ly α photons emerging from the RB16 galaxies are so readily detectable given the highly neutral state of the IGM at z > 7. One likely possibility is that these massive galaxies trace strong overdensities (Zitrin et al. 2015) and are hence situated in large ionized regions in the IGM (e.g. Furlanetto, Zaldarriaga & Hernquist 2004; Wyithe & Loeb 2005; Lee et al. 2007; McQuinn et al. 2007; Weinberger et al. 2018). These bubbles enable Ly α photons to redshift further into the damping wing before encountering H i, boosting the transmission (e.g. Mesinger, Haiman & Cen 2004; Mason & Gronke 2020). On the other hand, the high-EW [O iii]+H β emission (>800 Å) of the RB16 sample may suggest that their Ly α detections were driven more by physics internal to each of the four galaxies (Stark et al. 2017). Strong [O iii]+H β typically signals a recent rapid upturn or burst in star formation activity, i.e. high specific star formation rate (sSFR; Tang et al. 2019; Endsley et al. 2021, hereafter E21), giving rise to very young stellar populations. The spectra of such galaxies are therefore dominated by extremely hot stars that produce ionizing photons very efficiently with respect to the far-UV continuum. Assuming these ionizing photons are reprocessed into recombination lines, this should also boost the production rate of Ly α relative to the rest-UV continuum.

To better understand the origin of strong Ly α from massive reionization-era galaxies, we have initiated an MMT/Binospec campaign targeting Ly α in a much larger sample of bright (1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) z ≃ 7 systems. Here, we use the first results from our spectroscopic campaign to investigate to what extent strong [O iii]+H β (and hence large sSFRs) enhance Ly α detectability among bright reionization-era galaxies. Because our sources were selected over the very wide-area COSMOS and XMM3 fields (≈1.5 deg² each), we are able to assemble a much larger sample of these rare systems relative to CANDELS (≈0.2 deg² total; Grogin et al. 2011; Koekemoer et al. 2011).

We furthermore test whether our results are consistent with a picture wherein bright (and therefore likely massive; e.g. Barone-Nugent et al. 2014; Song et al. 2016a) z ≃ 7 galaxies often reside in large highly ionized bubbles, as would be expected if they preferentially trace overdensities. If this is the case, we would expect the transmission of Ly α photons from bright galaxies to show a substantially weaker decline at z ≳ 6 relative to much fainter sources. Recent observations suggest a sharp (factor of ∼10) decline in Ly α transmission for faint (∼0.1 L|$_{\mathrm{UV}}^{\ast }$|⁠) lensed galaxies between z ≃ 6–7 (Hoag et al. 2019; Fuller et al. 2020). Our goal is to test whether our bright (1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) z ≃ 7 galaxy sample exhibit a less rapid decline in Ly α transmission, building on earlier work targeting massive galaxies (e.g. Ono et al. 2012; Furusawa et al. 2016). For this investigation, we complement our z ≃ 7 targets with a collection of similarly bright z ≃ 6 sources that were observed simultaneously with the multiplex Binospec instrument.

This paper is organized as follows. In Section 2, we describe our observations including our sample selection, spectroscopic results, and the inferred [O iii]+H β EWs of each source. We present our analysis in Section 3, including an investigation of whether Ly α detectability is enhanced in bright z ≃ 7 systems with strong [O iii]+H β. We then discuss what our results imply for the presence of ionized regions around massive reionization-era galaxies in Section 4. Our main conclusions are summarized in Section 5.

All magnitudes are quoted in the AB system (Oke & Gunn 1983) and we adopt a flat ΛCDM cosmology with h = 0.7, Ω_M = 0.3, and Ω_Λ = 0.7, consistent with Planck results (Planck Collaboration VI 2020). All distances are quoted in physical units unless otherwise stated.

2 OBSERVATIONS

We have initiated a large spectroscopic Ly α survey of bright (≃1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) galaxies at z ≃ 7 with MMT/Binospec. Our selection criteria are described in Section 2.1. We then detail our spectroscopic results and the photometric properties of galaxies in our spectroscopic sample in Section 2.2. Finally, we infer the physical properties (e.g. stellar mass and [O iii]+H β EW) of each spectroscopic target in Section 2.3.

2.1 Source selection and photometry

In this work, we aim to accomplish two goals focused on targeting Ly α in bright (≳L|$_{\mathrm{UV}}^{\ast }$|⁠) galaxies at z ≃ 6–7. First, we seek to investigate whether strong [O iii]+H β (and hence large sSFRs) lead to enhanced Ly α emission among bright reionization-era galaxies. Secondly, we aim to quantify the evolution in the Ly α EW distribution of bright galaxies between z ≃ 6 and 7 to test whether our results are consistent with a picture wherein massive z ≃ 7 systems often reside in large ionized bubbles. Given the rarity of such luminous high-redshift galaxies, we select these sources over very wide-area fields, specifically the COSMOS and XMM3 fields that span 1.5 and 1.8 deg², respectively. We describe our z ≃ 7 and z ≃ 6 source selection criteria in turn below.

2.1.1 Selection of z ≃ 7 sources

To ensure we can reliably identify z ≃ 7 galaxies and infer their [O iii]+H β EWs, we photometrically select our targets over wide-area (>1 deg²) fields with deep optical through mid-infrared imaging. Because of the exceptionally deep imaging over COSMOS, the majority of our sources are selected within this field. The Subaru/Hyper Surpime-Cam (HSC) Subaru Strategic Program (HSCSSP; Aihara et al. 2018, 2019) provides optical (0.4–1μm) imaging across COSMOS in the g, r, i, z, y broad-bands as well as the nb816 and nb921 narrow-bands. The near-infrared data come from the UltraVISTA survey (McCracken et al. 2012) that delivers imaging in VISTA/VIRCam Y, J, H, and K_s broad-bands. We use PDR2 and DR4 data from HSCSSP and UltraVISTA, respectively, both of which are already astrometrically calibrated to the Gaia frame that we adopt throughout this work.

Our z ≃ 7 source selection over COSMOS largely follows that described in E21 that utilizes, in part, the nb921 photometry to limit the selection to z ≳ 6.6 with high confidence. The specific adopted Lyman-break colour cuts are z−y > 1.5, z−Y > 1.5, nb921−Y > 1.0, and y−Y < 0.4 where we set fluxes in the z and nb921 dropout filters to their 1σ values in cases of non-detections, consistent with previous literature (e.g. Bouwens et al. 2015a; Stefanon et al. 2019). These z ≃ 7 selection criteria were designed by simulating the colours of z = 6–8 galaxies with flat rest-UV slopes (β = −2 where F_ν ∝ λ^{β + 2}) using the Inoue et al. (2014) IGM transmission function. As with any Lyman-break selection utilizing broad-band photometry, the exact redshift selection window depends slightly on the assumed Ly α EW. Galaxies with weak Ly α emission (≃0 Å EW) are only likely to satisfy the above criteria from z ≃ 6.6–6.9 because at higher redshifts the y−Y colour becomes too red. On the other hand, sources with strong Ly α satisfy our selection criteria to larger maximum redshifts – up to z ≃ 7.1 and z ≃ 7.2 for Ly α EW = 25 Å and 100 Å, respectively. This is because strong Ly α will boost the y-band flux resulting in bluer y−Y colours at a given redshift. We refer the interested reader to fig. 2 of E21 for a visualization of the impact of Ly α on our redshift selection window. Because only ∼10 per cent of bright z ≃ 7 galaxies show strong Ly α emission (EW > 25 Å; Ono et al. 2012; Schenker et al. 2014; Pentericci et al. 2018), the large majority of our sample will likely lie at z ≃ 6.6–6.9. None the less, we directly account for how our redshift selection window depends on Ly α strength when inferring the Ly α EW distribution as described in Section 2.2.4.

Along with the above colour cuts, we ensure that each z ≃ 7 source is real by requiring a > 3σ detection in y, Y, and J as well as a >5σ detection in at least one of those three bands. We also enforce non-detections (<2σ) in g and r since both these bands probe flux blueward of the Lyman-continuum limit at z ≥ 6.6. Finally, we clean our sample of T-type brown dwarfs (which can exhibit very strong z and nb921-drops) by preserving only those with Y−J < 0.45 or (J−H > 0 and J−K_s >0). These cuts are guided by brown dwarf spectra in the SPEX library (Burgasser 2014) that demonstrate that T-type brown drafts have red Y−J colours but blue J−H and J−K_s colours.

To build statistics on the brightest end, we also select galaxies over the XMM3 field. While the optical through mid-infrared coverage over XMM3 is significantly shallower relative to COSMOS, we are still able to identify the most luminous (J ≲ 24.5) z ≃ 7 galaxies within this field. Near-infrared imaging over XMM3 comes from DR3 of the VIDEO survey² (Jarvis et al. 2013) that provides data from the Z, Y, J, H, and K_s VISTA/VIRCam broad-bands. For the optical imaging, we again use PDR2 of the HSCSSP that provides coverage in all the same bands as in COSMOS except for nb921. To compensate for the lack of nb921 imaging over XMM3, we enforce a much stronger dropout in the Z band, specifically requiring Z−y > 2.5 and Z−Y > 2.5 while maintaining the y−Y < 0.4 cut we applied in COSMOS. We found that such cuts select z ≃ 7 sources in a very similar redshift window as COSMOS after performing colour simulations of z = 6–8 galaxies similar to those described in E21. Here, we use VIRCam Z as the dropout band in XMM3 because the imaging in Z is often slightly deeper relative to z and both bands have very similar (normalized) transmission curves. All other selection criteria are equivalent to that over COSMOS.

To identify sources within each field, we run SExtractor (Bertin & Arnouts 1996) on a yYJHK_s χ² detection image (Szalay, Connolly & Szokoly 1999). We then calculate the optical and near-infrared photometry of each source in 1.2 arcsec diameter apertures which is ≈1.5 × the seeing in all bands. Aperture corrections are calculated from the median curve-of-growth of nearby isolated, unsaturated stars and photometric errors are determined using the standard deviation of flux within apertures randomly placed in nearby empty locations. As reported in E21, the typical 5σ depths in the COSMOS field are m = 27.7, 27.4, 27.3, 26.1, 26.9, 26.2, and 26.1 in g, r, i, nb816, z, nb921, and y, respectively. For the UltraVISTA data, the typical depths are m = 24.7, 24.5, 24.2, and 24.5 in Y, J, H, and K_s, respectively, for the deep stripes and m = 25.9, 25.8, 25.6, and 25.2, respectively, for the ultra-deep stripes. We are therefore able to identify z ≃ 7 galaxies as faint as J ∼ 26 in the ultra-deep stripes and J ∼ 25 in the deep stripes. In XMM3, the typical 5σ depths are m = 26.7, 26.3, 25.7, 25.3, 25.5, 24.3, 25.6, 25.0, 24.7, 24.3, and 23.9 in g, r, i, nb816, z, y, Z, Y, J, H, and K_s, respectively. We can therefore identify z ≃ 7 galaxies as faint as J ∼ 24.5 across XMM3.

2.1.2 Selection of z ≃ 6 sources

To quantify the evolution in the Ly α transmission for massive galaxies between z ≃ 6 and 7, we also observed bright z ≃ 6 sources that were photometrically selected as follows. Similar to our z ≃ 7 procedure, we design optimal z ≃ 6 selection criteria by simulating the HSC and VIRCam colours of z = 5–7 galaxies with rest-UV slope β = −2 and adopting the IGM transmission function from Inoue et al. (2014).

Our resulting z ≃ 6 selection criteria over COSMOS includes the following colour cuts:

• nb816−z > 1.5

• r− z > 1.5

• z− y < 0.5

• −0.5 <nb921−y < 0.5

Similar to our z ≃ 7 selection, fluxes in nb816 and r are set to their 1σ value in cases of non-detections. With these criteria, most selected bright z ≃ 6 sources (i.e. those with weak Ly α; Stark, Ellis & Ouchi 2011; De Barros et al. 2017) will lie at z ≃ 5.75–6.25 while strong Ly α emitters can be selected up to z ≃ 6.5. We have enforced the cut −0.5 <nb921−y < 0.5 to minimize the likelihood that strong z > 6.5 Ly α emitters are scattering into our z ≃ 6 selection. We additionally require a >5σ detection in z as well as a >3σ detection in nb921 and y to ensure each source is real. Finally, we enforce non-detections (<2σ) in g given that this band lies blueward of the Lyman-continuum break at z > 5.

In XMM3, our z ≃ 6 selection is largely similar. Specifically, we enforce nb816−Z > 1.5, r−Z > 1.5, and Z−Y < 0.5 (as detailed in Section 2.1.1, the Z and Y imaging are generally deeper than z and y, respectively, across XMM3). We also require a >5σ detection in Z as well as a >3σ detection in z and Y to ensure each source is real, in addition to non-detections (<2σ) in g. As in COSMOS, these cuts result in a redshift selection window of z ≃ 5.75–6.25 for sources with weak Ly α emission (EW ≃ 0 Å). Due to the lack of nb921 imaging across this field, sources with very strong Ly α emission (EW∼100 Å) can be selected up to z ≃ 6.6. However, our spectroscopic data confirm that all three of our z ≃ 6 XMM3 targets lie at z < 6.5 (Section 2.2) so this is not a concern for this work.

These z ≃ 6 sources are identified after running SExtractor on a znb921yYJHK_s and zZyYJHK|$_s\, \chi ^2$| detection image over COSMOS and XMM3, respectively. Optical and near-infrared photometry are calculated in the same way as for the z ≃ 7 sources (1.2 arcsec diameter apertures). Given the depths in each field, we are able to identify z ≃ 6 sources down to y ∼ 26.5 in COSMOS³ and Y ∼ 25.5 in XMM3. Note that throughout this work, for the z ≃ 6 galaxies, we quote y-band magnitudes for those in COSMOS and Y-band magnitudes for those in XMM3 given the depth differences noted above.

2.1.3 IRAC photometry

To infer the [O iii]+H β EWs of our z ≃ 7 sample, we take advantage of the Spitzer/IRAC imaging over both COSMOS and XMM3. A full description of our procedure for generating IRAC mosaics in both the 3.6 μm and 4.5 μm filters is provided in E21. Briefly, we use the mopex software (Makovoz & Khan 2005) to coadd background-subtracted images (using SExtractor) and astrometrically match the output mosaics to the Gaia reference frame using the iraf package ccmap. In COSMOS, the IRAC data comes from a multitude of surveys: the Spitzer Extended Deep Survey (SEDS; Ashby et al. 2013), S-CANDELS (Ashby et al. 2015), Star Formation at 4<z<6 from the Spitzer Large Area Survey with Hyper Suprime-Cam (SPLASH; Steinhardt et al. 2014), Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS; Ashby et al. 2018), and Completing the Legacy of Spitzer/IRAC over COSMOS (P.I. I. Labbé). The IRAC data over XMM3 largely comes from the Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012) with deeper imaging on our z ≃ 7 targets from late-2019 observations led by P.I. M. Stefanon.

To compensate for the considerably broader IRAC PSF relative to the optical/near-infrared seeing, we measure IRAC photometry in 2.8 arcsec diameter apertures and utilize a deconfusion algorithm to remove contaminating flux from neighbouring sources. In XMM3, our deconfusion approach is equivalent to that detailed in E21. To summarize, we convolve the flux profile of every nearby source detected in the yYJHK|$_s\, \chi ^2$| image with a 2D Gaussian having full width at half-maximum (FWHM) equal to the quadrature difference of the IRAC FWHM and the median seeing from each band in our χ² detection images. Flux profiles are calculated as the square root of the χ² image using the SExtractor segmentation map to determine source footprints. The convolved flux profiles are fit to the IRAC image with total fluxes of each source as free parameters. Once this is done, the best-fitting flux profile of each neighbouring source is subtracted before measuring the IRAC photometry. We note that our XMM3 targets are not strongly confused and that residuals from the deconfusion algorithm are acceptably smooth. That is the S/N of the pixels inside the aperture do not appear to be systematically offset due to poor flux-profile fitting of the neighbouring sources upon visual inspection. We show the deconfused IRAC postage stamps of both sources in Fig. A1 (see Appendix) to illustrate this point.

In COSMOS, we take advantage of the very high-resolution HST F814W imaging across this field (Scoville et al. 2007) to calculate the flux profile of each neighbouring source. By convolving these flux profiles with IRAC PSFs calculated using unsaturated stars near each source (<3 arcmin separation), we are able to obtain much smoother residual images from the deconfusion algorithm for sources lying in crowded regions. For the purposes of this work, we remove sources from our sample that have poor IRAC residuals after employing our deconfusion algorithm, as is common in the literature (e.g. Labbé et al. 2013; Bouwens et al. 2015a; Smit et al. 2015; De Barros et al. 2019). One of our COSMOS targets (COS-862541) lies outside the F814W imaging and we therefore employ the same deconfusion approach as for XMM3, noting that this source is also not strongly confused. The deconfused IRAC postage stamp images of all our targets are shown in Fig. A1 (see Appendix).

2.2 MMT/Binospec spectroscopy

We have followed up a subset of our z ≃ 6 and z ≃ 7 galaxy samples described in Section 2.1 using the Binospec spectrograph (Fabricant et al. 2019) installed at the MMT. Binospec is a wide-field (240 arcmin²) and multi-object (up to ∼150 sources) spectrograph enabling wavelength coverage up to ≈1 μm with moderately high resolution (R ≈ 4400). As such, Binospec is an ideal instrument to target Ly α in large numbers of z ≃ 6–7 galaxies selected over wide-area fields. So far, we have observed seven masks with five in COSMOS and two in XMM3, totaling just over 31 h of integration time across all masks. In Table 1, we report the central coordinate, position angle, central wavelength, total exposure time, and average seeing for each mask. When choosing which sources would be assigned slits for each mask, we assigned higher priority to sources with brighter UV magnitudes to obtain a more homogeneous sample across luminosity. Because galaxy rest-UV slopes typically become bluer at higher redshifts (e.g. Bouwens et al. 2014), we also gave higher priority to z ≃ 6 sources with bluer rest-UV slopes to obtain z ≃ 6 and z ≃ 7 samples with more similar properties.

2.2.1 Sample description

Across these seven masks, we have targeted 22 z ≃ 7 sources satisfying the selection criteria described in Section 2.1. In Table 2, we report the coordinates, J-band magnitudes, 3.6 and 4.5 μm photometry, as well as the [3.6]−[4.5] colour for each targeted z ≃ 7 galaxy. We also note in Table 2 the mask containing each source, where four sources were placed on both the COSc and COSd masks for a total integration time of 26 100 s. The J-band magnitudes of our z ≃ 7 targets range from J = 24.3 to 26.3. The vast majority of these sources show blue [3.6]−[4.5] colours (<0), as expected given that z ≳ 7 galaxies typically possess strong [O iii]+H β emission (Labbé et al. 2013; Smit et al. 2014; De Barros et al. 2019; E21). Furthermore, 8 of our 22 z ≃ 7 targets have very blue [3.6]−[4.5] colours (⁠|${\lt}-$|0.8) that translate to strong [O iii]+H β emission (EW ≳ 800 Å) assuming a flat rest-optical continuum slope and z = 6.75. This includes two sources (COS-221419 and COS-940214) that have 3.6 μm fluxes substantially above the 4.5 μm 2σ upper limit. Thus, even though these sources are not detected in 4.5 μm, the data still suggest that they likely possess strong [O iii]+H β emission.

Within these same seven masks, we have also targeted 30 z ≃ 6 galaxies satisfying the selection criteria from Section 2.1. The rest-UV apparent magnitudes⁴ of these z ≃ 6 targets range from |$m_{\mathrm{UV}}^{ }$| = 24.3–26.4, very similar to that of the z ≃ 7 targets (Fig. 1 top). The range of rest-UV slopes spanned by our z ≃ 6 spectroscopic targets (−1.2 ≤ β ≤ −3.5) is also very similar to that of the z ≃ 7 targets (−1.1 ≤ β ≤ −3.4). Furthermore, the typical rest-UV slopes of each sample are almost exactly equal – the median β of the z ≃ 7 sample is −2.07 and that of the z ≃ 6 sample is −2.09 (Fig. 1 bottom). For the z ≃ 7 sources, we calculate rest-UV slopes using YJHK_s, while for the z ≃ 6 sources we use nb921yYJHK_s and yYJHK_s in COSMOS and XMM3, respectively.

Figure 1.

Figures comparing the distribution of apparent UV magnitudes (top) and rest-UV slopes (bottom) for our spectroscopic z ≃ 6 and z ≃ 7 samples. We show the z ≃ 6 distributions with a shaded orange histogram while the z ≃ 7 distributions are shown with a thick blue line. We show median values for each sample with vertical dashed lines. The range of UV magnitudes and rest-UV slopes, as well as their typical values, are very similar for these two spectroscopic samples.

Open in new tabDownload slide

2.2.2 Data reduction

To design our Binospec multi-object slit masks, we used the binomask software. We adopted a slit width of 1.0 arcsec and the 600 l/mm grating, yielding a resolving power of R ≈ 4360. Central wavelengths for each mask range from 8500 to 8720 Å (see Table 1) with the exact value chosen to optimize the red-end wavelength coverage of our z ≃ 7 targets. For all z ≃ 7 targets, the maximum wavelength coverage was ≥9635 Å meaning that our observations covered Ly α up to at least |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 6.93. Furthermore, the typical wavelength coverage of our z ≃ 7 targets was 0.75–1.00 μm (⁠|${z}_{_{\mathrm{Ly\,\alpha }}}$| = 5.17–7.23). For the z ≃ 6 targets, our observations covered at least 7890–9510 Å and therefore (assuming they all lie at z > 5) fully encompassed the Ly α redshift range |${z}_{_{\mathrm{Ly\,\alpha }}}$| ≃ 5.75–6.50 set by our selection criteria (Section 2.1.2).

We adopted a slit length of at least 20 arcsec (7 arcsec) for the z ≃ 7 (z ≃ 6) targets that we found leads to sufficient modelling of the sky background from ≈0.9–1 μm (≈0.8–0.9 μm). Each mask also included at least five stars off the primary field of views⁵ (FoVs) for guiding and monitoring seeing throughout the observations. In addition, we placed three stars within the primary FoVs to determine exposure weighting (see below) and absolute flux calibration.

Each individual exposure was reduced separately using the publicly available Binospec data reduction pipeline (Kansky et al. 2019). Exposures for each mask were then co-added using the weighting scheme from Kriek et al. (2015) that takes the height of the telluric-corrected 1D flux profile of bright stars on the mask as the relative weight. This accounts for both the average sky transmission and relative seeing in each exposure. The average seeing of each mask (calculated from the co-added spectra) is listed in Table 1 and ranges from 0.9 to 1.2 arcsec.

We extract 1D spectra using optimal extraction (Horne 1986) after fitting a Gaussian⁶ to the observed emission line profile along the spatial axis. As in E21, absolute flux calibration is determined by calculating the average scaling factor that matches the 1D spectra of the three stars placed on the mask to their mean PSF z-band magnitudes from the Pan-STARRS survey (Chambers et al. 2016). Given the relatively narrow wavelength range covered by these observations (≈0.75-1 μm), we assume that this factor does not evolve with wavelength.

To estimate slit loss factors, we adopt the size–luminosity relation of bright z = 6–7 galaxies found by Curtis-Lake et al. (2016) that assumes a Sérsic profile with n = 1.0 (consistent with the approach of Ono et al. 2013). The modelled Sérsic profile of each source is convolved with a 2D Gaussian with FWHM equal to the seeing of the respective mask, and the fraction of flux within the 1 arcsec Binospec slit is compared to that for a point source. The resulting relative slit loss correction factors range from 4 to 10 per cent for all z ≃ 6 and z ≃ 7 targets, with brighter sources having larger corrections due to their larger assumed half-light radius. We note that our results are not significantly altered if we instead calculate slit loss correction factors assuming the size–luminosity relation from Bowler et al. (2017) that mainly yields larger sizes (and hence larger correction factors) for the brightest (⁠|$m_{\mathrm{UV}}^{ }$| < 25) sources.

We identify emission lines by first visually inspecting the 2D spectra of all sources. For emission features relatively clear of strong skylines, we calculate the flux in an aperture where the width along the wavelength axis is set by visual inspection of the 2D spectra with flux errors computed as described in E21. For these relatively clean emission features, we estimate that the fraction⁷ of flux potentially obscured by strong skylines is small (15–20 per cent) on average in both the z ≃ 6 and z ≃ 7 galaxies. We therefore conclude that our results on the evolution of the Ly α EW distribution between these two redshifts (Section 3.2) are not significantly impacted by this slight obscuration. Such small obscuration would furthermore only serve to slightly strengthen our conclusions⁸ on differences in the Ly α EW distribution for moderate versus strong [O iii]+H β emitters at z ≃ 7 (Section 3.1). Because the exact skyline obscuration fraction for a given source depends on the assumed line profile and accounting for it does not significantly alter our conclusions, we do not fold an obscuration correction for these relatively clean sources into our analysis.

For the few features that overlap more significantly with moderate-strength skylines, we adopt a slightly different approach to minimize the impact of skylines on the recovered flux (these sources are marked with asterisks in Tables 3 and 4). We fit the 1D spectra with a half-Gaussian (red-side only) convolved with the spectral resolution of our instrument (e.g. Hu et al. 2010) after masking portions of the spectrum contaminated by skylines. In this fitting procedure, we adopt a grid of three parameters describing the half-Gaussian (amplitude, standard deviation, and central wavelength) and compute the χ² value for each point in the grid. The likelihood for a given set of parameters is then calculated as |$P(A,\sigma ,\lambda _0) \propto \mathrm{ e}^{-\chi ^2/2}$| which we convert to a probability distribution on the flux of the emission line. The flux and its uncertainty are then computed as the median value and standard deviation from this probability distribution. For all emission features, we subtract the continuum flux estimated from the photometry⁹ but compute the significance of the feature prior to this subtraction.

2.2.3 Detected emission lines

We detect confident (>7σ) emission features in 9 of the 22 z ≃ 7 galaxies in our sample (see Fig. 2). In all cases, we interpret these features as Ly α because all are fully consistent with such a solution given the expected redshift range of our sample and none are consistent with an [O ii]λ3727,3729 solution (two narrow peaks of roughly similar strength separated by ≈6.8 Å in the observed frame). We consider [O ii]λ3727,3729 as the primary alternative solution because it would be very difficult for anything but a Balmer break to mimic the strong nb921 and Z drops required by our selection criteria. We also find no other convincing features in the spectra of these sources (aside from the tentative features in COS-469110 consistent with NV emission; see below). The measured line fluxes and rest-frame EWs for sources with a Ly α detection are reported in Table 3. We also calculate the corresponding Ly α redshifts, |${z}_{_{\mathrm{Ly\,\alpha }}}$|⁠, using the observed wavelength of peak flux in the 1D extraction and a rest-frame Ly α wavelength of 1215.67 Å. The detection for COS-862541 was previously reported in E21.

Figure 2.

2D signal-to-noise ratio maps (black is positive) of our nine confident (>7σ) Ly α detections at z ≃ 7. These detections span a redshift range of |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 6.650–7.093. Each sub-figure spans ±100 Å along the x-axis and ±7.3 arcsec along the y-axis. We mask the bright continuum from nearby sources in the spectra of COS-955126 and XMM3-227436 to improve clarity.

Open in new tabDownload slide

Our z ≃ 7 Ly α detections span redshifts of |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 6.650–7.093. Measured fluxes and EWs range from (5.1–18.6) × 10⁻¹⁸ erg s⁻¹ cm⁻² and 3.7–43.1 Å, respectively. None of the detections lie at |${z}_{_{\mathrm{Ly\,\alpha }}}$| < 6.6, consistent with the strong nb921 and Z drops used in our selection. Furthermore, the large majority of detected targets lie at |${z}_{_{\mathrm{Ly\,\alpha }}}$| ≃ 6.6–6.9. While two of our sources (COS-1205190 and XMM3-227436) lie at z = 7.05–7.09, these redshifts are still consistent with our selection given the moderate-strength Ly α emission for these sources (EW = 15–29 Å; see Section 2.1.1).

We also detect confident (>7σ) emission features in 13 of our 30 z ≃ 6 targets. In Table 4, we report the coordinates, rest-UV slopes, Ly α redshifts, fluxes, EWs, and y/Y-band magnitudes of each detected z ≃ 6 source in COSMOS/XMM3. The 2D spectra of all z ≃ 6 detections are shown in Fig. 3 and span redshifts of |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 5.824–6.303, consistent with expectations given their selection criteria (Section 2.1.2). We measure Ly α fluxes ranging from (6.8–62.3) × 10⁻¹⁸ erg s⁻¹ cm⁻² and EWs ranging from 5.5–113 Å.

Figure 3.

2D signal-to-noise ratio maps of our 13 Ly α detections at z ≃ 6 where black is positive. These detections span a redshift range of |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 5.824–6.303. Each sub-Fig. spans ±100 Å along the x-axis and ±2.5 arcsec along the y-axis.

Open in new tabDownload slide

Motivated by recent detections in similarly luminous z ≳ 7 galaxies (Tilvi et al. 2016; Laporte et al. 2017; Mainali et al. 2018; see also Hu et al. 2017), we search for NVλ1238.8,1242.8 emission in the spectra of our z ≃ 7 Ly α emitters. Given the very high ionization potential of this line (77 eV), any such detection would likely signal significant active galactic nucleus (AGN) activity. As shown in Fig. 4, we tentatively identify the NV doublet in COS-469110 (⁠|$\mathrm{\mathit{ M}}_{\mathrm{UV}}^{ }$| = −21.7) where the emission is located at the exact expected spatial position of the source. We regard both these features as tentative because each is detected at <7σ significance, specifically 6.4σ and 5.2σ for the 1238.8 and 1242.8 Å components, respectively. However, we note that if we split our data into two distinct stacks,¹⁰ the NVλ1238.8 feature is significantly detected in each stack (5.2σ and 3.8σ), adding evidence that this feature is likely real.

Figure 4.

MMT/Binospec spectra of COS-469110 where we identify Ly α as well as tentative detections of the NVλ1238.8,1242.8 doublet, a signpost of significant AGN activity. The top and middle panels show the 2D and smoothed 2D signal-to-noise ratio maps, respectively, where black is positive. In the bottom panel, we show the 1D extraction centred on the expected spatial position of the source with the 1σ noise level in grey. The two NV doublet components are measured with a significance of 6.4σ and 5.2σ, respectively, and have peak wavelengths corresponding to the exact same redshift of z = 6.645. This translates to a Ly α velocity offset of ≈200 km s⁻¹ that may help explain the exceptionally high NV to Ly α flux ratio of this source (see text).

Open in new tabDownload slide

The NVλ1238.8 component is clearly distinct from strong skylines and has a measured flux of (5.4 ± 1.0) × 10⁻¹⁸ erg s⁻¹ cm⁻² which, using the Y-band photometry for the continuum, corresponds to an EW of 8.4 ± 3.0 Å. The flux of the NVλ1242.8 component is measured to be (3.8 ± 0.8) × 10⁻¹⁸ erg s⁻¹ cm⁻² that corresponds to an EW of 5.9 ± 2.2 Å. While the 1242.8 Å component does rest up against a skyline, the peak wavelengths of both components correspond to the exact same redshift of z = 6.645 that is well consistent with the Ly α redshift of |${z}_{_{\mathrm{Ly\,\alpha }}}$| = 6.650 (we discuss the implied Ly α velocity offset below). Furthermore, the measured EWs of both components are consistent with that recently reported in three other similarly luminous z = 7–9 galaxies (Tilvi et al. 2016; Laporte et al. 2017; Mainali et al. 2018) as well as a 1:1 to 2:1 flux ratio for the doublet (Bickel 1969; Torres-Peimbert & Pena 1984). However, we do note that calculating this flux ratio is complicated by the partial skyline masking of the 1242.8 Å component.

While we do not detect any significant NV features in any of our other z ≃ 7 Ly α emitters, we are able to place 5σ NVλ1238.8 EW limits¹¹ of ≤10 Å in six of our other 12 z ≃ 7 Ly α emitters, suggesting that the (tentative) NV emission from COS-469110 is likely exceptional among the bright z ≃ 7 population. This is perhaps further supported by the fact that the total NV flux we measure implies a line flux ratio of |$f_{_{\mathrm{NV}}}/f_{_{\mathrm{Ly}\,\alpha {}}}$| = 0.95 ± 0.19, much larger than typical upper limits recently placed on several other bright z > 6.5 Ly α emitters (⁠|$f_{_{\mathrm{NV}}}/f_{_{\mathrm{Ly}\,\alpha {}}}\, \lesssim$| 0.2; Mainali et al. 2018; Shibuya et al. 2018). Because of the partial skyline obscuration of the 1242.8 Å component, this flux ratio for COS-469110 is likely a lower limit.

The Ly α velocity offset implied for COS-469110 is ≈200 km s⁻¹, assuming that NV well traces the systemic redshift (Laporte et al. 2017). This velocity offset falls well within the range previously reported for similarly luminous (⁠|$\mathrm{\mathit{ M}}_{\mathrm{UV}}^{ }$|  |$\lt \, -$|21.5) galaxies at z > 6 (110–500 km s⁻¹; Willott et al. 2015; Inoue et al. 2016; Pentericci et al. 2016; Laporte et al. 2017; Stark et al. 2017; Mainali et al. 2018; Matthee et al. 2020), though it does sit on the lower end of that range. This may help explain the exceptionally large NV to Ly α line flux ratio mentioned above as Ly α emission is more susceptible to strong scattering by the partially neutral IGM at lower velocity offsets (e.g. Miralda-Escudé 1998; Mason et al. 2018b).

2.2.4 Completeness simulations

When inferring the Ly α EW distribution in Section 3.2, we fold in constraints from non-detections. We do so by calculating the Ly α detection completeness of each non-detected source as a function of Ly α EW from EW = 0–100 Å in 1 Å steps. We insert 10 000 simulated 1D flux profiles for each value of EW into the fully reduced 1D spectrum of each source at its expected spatial position and determine what fraction of simulated features would be detected at >7σ. The observed wavelengths of these simulated features are randomly sampled using the redshift completeness distributions as a function of Ly α EW and rest-UV magnitude described in Section 2.1. In doing so, we account for the dependence of our redshift selection window on Ly α EW and for the impact of skyline obscuration. The profile of each simulated Ly α feature is a Gaussian with FWHM set to 220 km s⁻¹, equal to the mean value measured from the two z ≃ 7 Ly α features detected with very high significance (>15σ) in our sample, namely COS-788571 (FWHM = 180 km s⁻¹) and COS-1009842 (FWHM = 260 km s⁻¹). This FWHM is also in excellent agreement with that found by Pentericci et al. (2018) using stacked spectra of fainter z ≃ 7 galaxies.

We plot the simulated spectroscopic completeness for a representative subset of the non-detected z ≃ 7 (z ≃ 6) galaxies in the top (bottom) panel of Fig. 5. The spectral constraints on most of these z ≃ 6–7 targets yield ≳50 per cent completeness at EW = 20 Å, increasing to ≳80 per cent at EW = 60 Å. For a small subset of the non-detected z ≃ 7 sources, the spectral completeness declines at high Ly α EW due to their lack of very red (≳9750 Å) spectral coverage that prohibits detecting Ly α at z ≳ 7.0. Such a high-redshift solution is only likely for strong Ly α emitters given our selection criteria (Section 2.1.1). All of this information is folded into our analysis when inferring the Ly α EW distributions (Section 3).

Figure 5.

Figures showing our simulated spectroscopic completeness as a function of Ly α EW for non-detected sources at z ≃ 7 (top) and z ≃ 6 (bottom). To improve clarity, we only show the completeness curves for a representative subset of non-detected sources. For a small subset of the non-detected z ≃ 7 sources, the spectral completeness declines at high Ly α EW due to their lack of very red (≳9750 Å) spectral coverage that prohibits detecting Ly α at z ≳ 7.0 and such a high-redshift solution is only likely for strong Ly α emitters (Section 2.1.1). All of this information is folded into our analysis when inferring the Ly α EW distributions.

Open in new tabDownload slide

2.3 Photoionization modelling

We now infer the [O iii]+H β EWs of each of our z ≃ 7 targets to test whether we see a correlation with Ly α visibility. To infer the [O iii]+H β emission strength (as well as other physical properties such as stellar mass), we use the BayEsian Analysis of GaLaxy sEds (beagle; Chevallard & Charlot 2016) SED fitting code. beagle computes both the stellar and nebular emission of star-forming galaxies by adopting the photoionization models from Gutkin, Charlot & Bruzual (2016), which are derived by incorporating the latest version of Bruzual & Charlot (2003) stellar population synthesis models into cloudy (Ferland et al. 2013). Our SED fitting procedure with beagle matches that described in E21. Briefly, we assume a delayed star formation history (SFR ∝ t e^−t/τ) with an allowed recent (<10 Myr) burst, a minimum age of 1 Myr, an sSFR ranging from 0.1 Gyr⁻¹ to 1000 Gyr⁻¹, and an SMC dust prescription.

We note that some of our spectroscopically confirmed sources lie at |${z}_{_{\mathrm{Ly\,\alpha }}}$| ≈ 6.85–6.9 (COS-788571, COS-862541, and XMM3-504799) where the transmission of [O iii]λ5007 through [3.6] (and hence the inferred [O iii]+H β EW) is very sensitive to the exact systemic redshift, z_sys. At z ≳ 6, visible Ly α emission is often redshifted relative to systemic due to complex radiative transfer effects in the ISM, CGM, and IGM (see e.g. Dijkstra 2014 for a review). The current maximum observed velocity offset at z > 6 is 500 km s⁻¹ (Willott et al. 2015) that translates to z_sys = |${z}_{_{\mathrm{Ly\,\alpha }}}$|−0.013 at z = 6.6–7.1 where we have z ≃ 7 Ly α detections. Therefore, during the beagle fitting process, we allow z_sys to range between |${z}_{_{\mathrm{Ly\,\alpha }}}$|−0.013 and |${z}_{_{\mathrm{Ly\,\alpha }}}$| for z ≃ 7 sources with a Ly α detection.

For sources with a Ly α detection, we furthermore exclude bands blueward of the Ly α break during the fitting process. We also do not fit to bands significantly impacted by Ly α emission for these sources as it is currently not possible to leave the effective Ly α transmission fraction through the IGM as a free parameter in beagle. This means we fit to all bands redder than and including Y for sources with |${z}_{_{\mathrm{Ly\,\alpha }}}$| < 6.9 and all bands redder than including J for sources with |${z}_{_{\mathrm{Ly\,\alpha }}}$| > 6.9. Sources without a Ly α detection are fit using all optical through mid-infrared photometry with a uniform redshift prior of z = 6–8 and we remove Ly α emission from the nebular templates given the low (∼10 per cent) fraction of strong Ly α emitters (EW > 25 Å) found among the bright z ≃ 7 population in previous works (e.g. Ono et al. 2012; Schenker et al. 2014; Pentericci et al. 2018). We report the inferred [O iii]+H β EWs of each z ≃ 7 source in Table 5, along with their inferred absolute UV magnitudes (at 1600 Å rest-frame), stellar masses, V-band optical depths, and sSFRs. In the final column, we report the inferred ionizing photon production efficiencies, |$\xi _{\mathrm{ion}}^{\ast }$|⁠, computed using the intrinsic UV luminosity (at 1500 Å rest-frame) of the stellar population before processing through dust and gas (see Chevallard et al. 2018 for details). We also quote photometric redshifts for sources without a Ly α detection and |${z}_{_{\mathrm{Ly\,\alpha }}}$| for sources with a Ly α detection in Table 5.

Overall, we find a diverse range of galaxy properties within our sample as illustrated in Fig. 6. In particular, the inferred [O iii]+H β EWs range from 200 to 4000 Å and the sSFRs span ∼1–150 Gyr⁻¹. Galaxies with inferred [O iii]+H β EW > 800 Å (and therefore similar to those in the RB16 sample) tend to possess very large sSFRs (≳30 Gyr⁻¹), consistent with expectations of galaxies experiencing a recent strong upturn or burst in star formation activity (Atek et al. 2011; Amorín et al. 2015; Maseda et al. 2018; Tang et al. 2019; E21). As expected (Chevallard et al. 2018; Tang et al. 2019), we also find that |$\xi _{\mathrm{ion}}^{\ast }$| tends to increase with [O iii]+H β EW among our sample where sources with [O iii]+H β EW > 800 Å are inferred to produce ionizing photons 2.2× as efficiently as those with weaker [O iii]+H β. The median inferred [O iii]+H β EW of all sources is 710 Å, suggesting that our sample is reasonably well representative of the global z ≃ 7 galaxy population (typical EW = 692 Å; E21).

Figure 6.

Spectral energy distributions of a subset of z ≃ 7 sources in our spectroscopic sample. We select sources to display the variety of absolute UV magnitudes, |$\mathrm{\mathit{ M}}_{\mathrm{UV}}^{ }$|⁠, and inferred [O iii]+H β EWs (generally increasing to the right) in this sample. In each panel, the median fit beagle SED model is overlaid in black with grey shading showing the inner 68 per cent confidence interval from the posterior output by beagle. Blue diamonds show the fitted photometry (with 2σ upper limits in cases of non-detections) while the red empty circles show the median model photometry from beagle with errors enclosing the 68 per cent confidence interval. We only show model photometry for bands used in the fitting process where we ignore bands impacted by Ly α emission or bluewards of the Ly α break for sources with a >7σ Ly α detection.

Open in new tabDownload slide

The absolute UV magnitudes of our z ≃ 7 targets also span −22.5 ≤ |$\mathrm{M}_{\mathrm{UV}}^{ }$|  |$\le \, -$|20.4 that corresponds to a luminosity range of 0.8–5.8 L|$_{\mathrm{UV}}^{\ast }$| when adopting the z ≃ 7 |$\mathrm{M_{UV}^{\ast }}$| value¹² of −20.6 from Bowler et al. (2017). Our survey is therefore capturing a subset of the most luminous reionization-era galaxies. Furthermore, over half the galaxies in our z ≃ 7 sample are >2 L|$_{\mathrm{UV}}^{\ast }$|⁠, a population that has been largely missed by previous spectroscopy targeting this epoch. Considering the published literature of z ≃ 7 Ly α observations from Fontana et al. (2010), Pentericci et al. (2011, 2014), Vanzella et al. (2011), and Schenker et al. (2012, 2014), only 4 of 64 targeted galaxies were >2 L|$_{\mathrm{UV}}^{\ast }$|⁠. The only previous spectroscopic z ≃ 7 studies with a comparable luminosity range to our sample are Ono et al. (2012) and Furusawa et al. (2016). Even so, our sample is at least twice as large as either of these previous studies.

The bright UV luminosities of our z ≃ 7 targets suggest that they are among the most massive galaxies present in the reionization era, given both clustering measurements (Barone-Nugent et al. 2014) as well as inferences of the |$\mathrm{\mathit{ M}}_{\mathrm{UV}}^{ }$|–M_* relation (Song et al. 2016a). Stellar masses derived from our SED fitting procedure suggest a similar picture. Most of our z ≃ 7 galaxies are inferred to have log₁₀(M_*/M_⊙) ∼ 9–10, which is significantly more massive than fainter z ≃ 7 galaxies identified over the GOODS fields (E21). Therefore, in what follows, we interchangeably use the terms ‘bright’ and ‘massive’ when describing our sample.

3 ANALYSIS

In this section, we first investigate how Ly α detectability correlates with [O iii]+H β emission (and hence sSFR) within our sample of 22 bright (L|$_{\mathrm{UV}}^{}\, \simeq$| 1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) z ≃ 7 galaxies (Section 3.1). We then quantify evolution in the Ly α EW distribution of bright galaxies between z ≃ 6–7 to test whether our results are consistent with accelerated reionization around massive z ≃ 7 systems (Section 3.2). Finally, we investigate the spatial separations of our z ≃ 7 Ly α emitters to identify any potentially large ionized bubbles within our observed fields (Section 3.3).

3.1 The Connection between [O iii]+H β and Ly α at z ≃ 7

The number of Ly α detections among UV-continuum selected galaxies in the reionization era (z ≳ 7) has grown rapidly over the past decade (e.g. Vanzella et al. 2011; Ono et al. 2012; Schenker et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Pentericci et al. 2014, 2018; Song et al. 2016b; Jung et al. 2017, 2018; Castellano et al. 2018; Larson et al. 2018; Hoag et al. 2019; Fuller et al. 2020; Tilvi et al. 2020). One of the most surprising results from these campaigns was the 100 per cent Ly α detection rate among the four brightest (L|$_{\mathrm{UV}}^{}$| = 3–4 L|$_{\mathrm{UV}}^{\ast }$|⁠) z > 7 galaxies selected over the CANDELS fields (i.e. the RB16 sample). While this led some to suggest that massive reionization-era systems may trace accelerated sites of reionization (Zitrin et al. 2015; Stark et al. 2017), this interpretation was complicated by the fact that the RB16 galaxies were selected to not only be very massive, but also to possess strong [O iii]+H β emission (EW > 800 Å). Such strong nebular emission signals the presence of intense radiation fields likely powered by a recent rapid upturn in star formation activity (i.e. high sSFR; Tang et al. 2019; E21). This therefore raises the possibility that the Ly α detections within the RB16 sample were not necessarily due to large ionized regions, but perhaps instead driven by physics internal to each of the four galaxies (Tang et al. 2019).

To better understand the origin of strong Ly α from the RB16 sample, we here analyse the Ly α constraints from our larger sample (N = 22) of similarly massive (L|$_{\mathrm{UV}}^{}\, \simeq$| 1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) z ≃ 7 galaxies possessing a wide range of inferred [O iii]+H β EWs (200–4000 Å). This diversity of nebular emission strengths enables us to explore how the sSFRs of massive reionization-era galaxies impact their Ly α detectability. To this end, we divide our sample into galaxies with strong (EW > 800 Å) versus more moderate (EW = 200–800 Å) inferred [O iii]+H β emission. By this classification, galaxies in the strong [O iii]+H β emitter sub-sample will possess large sSFRs comparable to those in the RB16 sample while the moderate [O iii]+H β emitters are more representative of the typical massive z ≃ 7 population (E21). Here, we do not include COS-469110 in either sub-sample given its tentative NV detection (Fig. 4) that is a likely signpost of significant AGN activity. We do, however, discuss the potential implications of this unique source below.

An essential question for interpreting the RB16 result is whether our z ≃ 7 targets with strong [O iii]+H β emission show an increased Ly α detection rate. From our observations, we confidently (>7σ) detect Ly α in seven of nine (78 per cent) galaxies with strong [O iii]+H β emission, yet in only one of 12 (8 per cent) galaxies with more moderate [O iii]+H β (see Fig. 7). Our results therefore suggest that Ly α is indeed more easily detectable from massive z ≃ 7 galaxies possessing very large sSFRs (≳30 Gyr⁻¹). We note that the two strong [O iii]+H β emitters in our sample that went undetected may exhibit prominent Ly α emission that is masked by skylines. This possibility is accounted for in our analysis below.

Figure 7.

Plot showing the inferred [O iii]+H β EWs and measured rest-UV slopes for each of our 22 massive (L|$_{\mathrm{UV}}^{}\, \simeq$| 1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) z ≃ 7 systems. We mark those with and without confident (>7σ) Ly α detections as shown in the legend. As indicated by the prominence of blue markers to the right of the plot, we confidently detect Ly α at a much higher rate from strong [O iii]+H β emitters (>800 Å EW; 78 per cent detection rate) than more moderate [O iii]+H β emitters (200–800 Å EW; 8 per cent detection rate). Our results therefore suggest that Ly α is more readily detectable from massive z ≃ 7 galaxies experiencing a rapid upturn in star formation activity (i.e. high sSFR). We also use a red star to show COS-469110 which likely harbors an AGN given its tentative NV detection (see Fig. 4).

Open in new tabDownload slide

We find that massive z ≃ 7 galaxies with strong [O iii]+H β emission typically exhibit much stronger Ly α (4.7× higher EW) relative to more moderate [O iii]+H β emitters when comparing the best-fitting μ_EW values for each sub-sample (see top panel of Fig. 8). Therefore, while uncertainties remain significant,¹⁵ our results suggest that galaxy sSFRs play a significant role in regulating Ly α emission from massive reionization-era galaxies. The recent results of Castellano et al. (2017) provide additional empirical support of this picture. In their study, they find that the stacked IRAC colours of z ≃ 6.8 Ly α emitters selected over the CANDELS fields indicate much stronger [O iii]+H β emission relative to galaxies that went undetected in Ly α. This is consistent with our results using [O iii]+H β EWs inferred for individual (and generally more massive) galaxies.

Figure 8.

Comparison of the inferred Ly α EW distributions for our z ≃ 7 sample as a function of physical properties. In the top panel, we compare the Ly α EW distributions for moderate (orange dashed curve) versus strong [O iii]+H β emitters (solid blue). In the bottom panel, we compare the inferred Ly α EW distributions for our z ≃ 7 galaxies with rest-UV slopes of β > −2.1 (dashed orange) versus the bluest sources with β < −2.1 (solid blue curve). These plots demonstrate that massive z ≃ 7 galaxies with strong [O iii]+H β emission and/or very blue rest-UV slopes show evidence for systematically higher EW Ly α emission.

Open in new tabDownload slide

Interestingly, two sources in our sample, COS-469100 and COS-1205190, both show fairly strong Ly α (EW = 15–29 Å), yet possess only moderate [O iii]+H β emission (EW = 200–800 Å). It is therefore of interest to explore whether there are other reasons to believe that these two sources may still possess unusually powerful radiation fields as expected in very large sSFR systems. Indeed, the tentative NV detection within COS-469110 (see Fig. 4) signals that this source likely harbors an AGN. We also find reason to believe that COS-1205190 may be powering an intense radiation field given its extremely blue rest-UV slope of β = −3.44 ± 1.03. Such a blue slope (the bluest in our sample; Table 6) is consistent with not only extremely low dust content, but also very low metallicity. This possible dearth of metals in COS-1205190 would naturally lead to relatively weak [O iii] emission¹⁶ even if it recently experienced a burst of star formation activity (sSFR ≳ 30 Gyr⁻¹) as typically inferred for the strong [O iii]+H β emitters in our sample. If our above suspicions of COS-469100 and COS-1205190 are correct, all nine of our z ≃ 7 Ly α detections are from sources with intense radiation fields and we would infer that such systems typically exhibit substantially stronger (⁠|$5.9_{-3.1}^{+4.3} \times $| higher EW) Ly α relative to the more typical massive z ≃ 7 population that has moderate [O iii]+H β (EW = 200–800 Å), moderately low metallicity (∼0.2 Z_⊙; E21), and emission dominated by star formation.

The strong Ly α emission (EW = 29 Å) from our bluest source (COS-1205190) also motivates an investigation into how the observed Ly α EW is related to the rest-UV slope among massive z ≃ 7 galaxies. While it has been shown that Ly α becomes stronger in z ≃ 4 galaxies from β = −1.4 to β = −1.8 (Stark et al. 2010), it is not necessarily clear that such a prominent trend would continue to exist when the bulk of the galaxy population is very blue (β ≲ −2) as is the case at z ≃ 7 (Bouwens et al. 2014). To test for any such association, we split our sample (again ignoring COS-469110) into sources with β < −2.1 and β > −2.1, where we adopt a dividing point equal to the approximate median rest-UV slope of the entire sample (β = −2.07). Following the Bayesian approach used for the [O iii]+H β versus Ly α analysis above, we find that our bluest galaxies (−3.4 < β < −2.1) typically exhibit much stronger Ly α (⁠|$4.0^{+4.1}_{-1.9} \times$| higher EW) relative to those more representative of the massive z ≃ 7 population (−1.1 < β < −2.1; E21) as illustrated in the bottom panel of Fig. 8. This trend likely reflects lower dust content (and therefore relatively little Ly α attenuation) within bluer sources as is the case at lower redshifts (e.g. Shapley et al. 2003; Pentericci et al. 2009; Kornei et al. 2010; Stark et al. 2010; Hathi et al. 2016; Trainor et al. 2016). Notably, our z ≃ 7 galaxies also follow a similar behaviour to the z ≃ 6 sources presented in De Barros et al. (2017) where our strongest Ly α emitters (EW > 25 Å) all possess very blue rest-UV slopes (β < −2.1; see Table 6) likely signaling particularly low dust content along the line of sight.

The overall picture emerging from these results is that the observed Ly α EW is substantially enhanced by a recent strong burst of star formation (or presence of an AGN) among massive reionization-era galaxies. Low dust content along the line of sight will further enhance visible Ly α similar to trends at lower redshifts. While we defer a more detailed physical interpretation until Section 4, we emphasize that our results demonstrate that Ly α can be detected with high success rate in massive z ≃ 7 galaxies possessing strong [O iii]+H β emission (>800 Å EW) and very blue rest-UV slopes (β < −2.1).

3.2 Evolution in the Ly α EW distribution of massive galaxies from z ≃ 6 to z ≃ 7

Over the past decade, a variety of observational campaigns have presented evidence that the IGM rapidly transitions from highly ionized at z ≃ 6 (x_H i ∼ 10 per cent; e.g. McGreer et al. 2015) to substantially neutral at z ≃ 7 (x_H i ∼ 50 per cent; e.g. Greig & Mesinger 2017; Zheng et al. 2017; Davies et al. 2018; Wang et al. 2020; Whitler et al. 2020). Such a rapid transition is naturally expected to reduce Ly α transmission between z ≃ 6–7 for the typical galaxy population. This is perhaps evidenced by the factor of ∼10 decline in the strong Ly α (EW > 25 Å) emitter fraction among the low-mass lensed population (∼0.1 L|$_{\mathrm{UV}}^{\ast }$|⁠; Hoag et al. 2019; Fuller et al. 2020). One of the primary goals of our spectroscopic campaign is to build a sufficiently large sample of very bright z ≃ 6–7 galaxies to explicitly test whether Ly α transmission from massive systems is evolving at a slower pace, as would be expected if they commonly reside in large ionized bubbles (e.g. Wyithe & Loeb 2005; McQuinn et al. 2007; Weinberger et al. 2018). While this campaign is still ongoing, we here report our current constraints on the Ly α EW distribution evolution between z ≃ 6–7 using our sample of 30 and 22 massive (L|$_{\mathrm{UV}}^{}\, \simeq$| 1–6 L|$_{\mathrm{UV}}^{\ast }$|⁠) galaxies at z ≃ 6 and z ≃ 7, respectively.

At each redshift, we infer the Ly α EW distribution following the Bayesian approach presented in Section 3.1 where we assume a lognormal function parametrized by a median EW, μ_EW, and standard deviation, σ_EW. Here, we also add a third parameter, A, quantifying the fraction of sources with Ly α EW > 0 Å because some galaxies may not show any Ly α in emission (e.g. Schenker et al. 2014; Mason et al. 2018a). Using the same grid for μ_EW and σ_EW as in Section 3.1 and allowing A to vary between 0 and 1 (with a spacing of 0.01), we calculate the probability for each set of parameters P(μ_EW, σ_EW, A) as in equation (1).

The Ly α EW distribution of our 22 massive z ≃ 7 galaxies is well characterized by a lognormal distribution with a median EW μ_EW = |$11.0^{+3.8}_{-3.2}$| Å and standard deviation σ_EW = |$0.37^{+0.15}_{-0.10}$| dex. The fraction of EW > 0 Å sources is found to be very high (A = |$0.88^{+0.09}_{-0.15}$|⁠) with a best-fitting value of A = 1, suggesting that a large majority of massive z ≃ 7 galaxies show Ly α in emission. This EW distribution further implies that approximately 15 per cent of our z ≃ 7 galaxies exhibit strong (EW > 25 Å) Ly α, in agreement with previous studies quantifying this Ly α emitter fraction among the bright (⁠|$\mathrm{\mathit{ M}}_{\mathrm{UV}}^{ }$|  |$\lt \, -$|20.25) population (e.g. Ono et al. 2012; Schenker et al. 2014; Pentericci et al. 2018). Interestingly, the Ly α EW distribution inferred from our 30 massive z ≃ 6 galaxies is remarkably similar to that at z ≃ 7 with all three parameters agreeing well within 1σ uncertainties (see Table 7).

Figure 9.

Plot showing the inferred Ly α EW distributions of our massive galaxy samples at z ≃ 7 (solid blue) and z ≃ 6 (dashed orange) using the Bayesian approach of equation (2). These two distributions are equivalent within uncertainties, suggesting that Ly α transmission is evolving less rapidly between z ≃ 6–7 for the massive population relative to low-mass lensed systems (Hoag et al. 2019; Fuller et al. 2020). This is consistent with a scenario wherein massive z ≃ 7 galaxies often reside in large highly ionized bubbles.

Open in new tabDownload slide

Our goal is to quantify evolution in the transmission of Ly α for massive galaxies between z ≃ 6–7. While the value of T inferred above does account for evolution in transmission, it also includes any evolution in the physical conditions that may be impacting Ly α production and escape within galaxies. Fortunately, we can take the first steps towards decoupling these internal factors using the insight gained from our analysis in Section 3.1. Therein, we found that the observed Ly α EW in massive z ≃ 7 galaxies correlates with rest-UV slope and [O iii]+H β EW. As detailed in Section 2.2, the typical rest-UV slopes of our z ≃ 6 (β = −2.09) and z ≃ 7 (β = −2.07) samples are essentially equal. Given this, we assume that the dust content of our two samples are similar enough to not cause a significant evolution in the Ly α EW distribution. As discussed further in Section 4, the trend between [O iii]+H β EW and Ly α EW seen in our sample is likely driven by the relationship between [O iii]+H β EW and the stellar ionizing photon production efficiency, |$\xi _{\mathrm{ion}}^{\ast }$| (Chevallard et al. 2018; Tang et al. 2019). Our SED fits with beagle suggests that the typical ionizing photon production efficiencies of our z ≃ 6 (⁠|$\log _{10}[\xi _{\mathrm{ion}}^{\ast } / (\mathrm{erg}^{-1}\ \mathrm{Hz})]$| = 25.62) and z ≃ 7 (⁠|$\log _{10}[\xi _{\mathrm{ion}}^{\ast } / (\mathrm{erg}^{-1}\ \mathrm{Hz})]$| = 25.59) samples are also similar enough to not significantly impact the evolution of the Ly α EW distribution. We therefore assume that our inferred value of T is dominated by the evolution in Ly α transmission.

Equipped with this knowledge, we find that our spectroscopic data rule out a strong decline in Ly α transmission (T_{z ≃ 7} / T_{z ≃ 6} < 0.5) with 98.5 per cent (≈2σ) confidence. This is in contrast to the ten-fold decrease in strong Ly α emitters among low-mass lensed galaxies over the same redshift interval (Hoag et al. 2019; Fuller et al. 2020). Therefore, while uncertainties are significant given our modest sample size, our results currently suggest that Ly α transmission is evolving less rapidly between z ≃ 6–7 for the massive population.

There are a number of possible explanations for this difference in Ly α transmission evolution as we enter the epoch of reionization. One is that massive z ≃ 7 galaxies may tend to possess larger Ly α velocity offsets, helping push the photons into the damping wing before escaping the host galaxy (Mason et al. 2018b). This may be expected given the positive correlation between Ly α velocity offset and rest-UV luminosity at z ∼ 2–3 (Erb et al. 2014), likely reflecting (in part) the ability of more massive galaxies to drive stronger outflows. Another explanation is that massive z ≃ 7 galaxies are more likely to reside in large ionized bubbles (Furlanetto et al. 2004; Wyithe & Loeb 2005; Lee et al. 2007; McQuinn et al. 2007; Weinberger et al. 2018) which enable Ly α photons to cosmologically redshift far into the damping wing and thus transmit more easily through a partially neutral IGM.

3.3 A possible large ionized region at z ≃ 7

As discussed in Section 3.2, our results are consistent with a scenario wherein massive z ≃ 7 systems often reside in large ionized regions. In this same picture, we might expect to find some of our z ≃ 7 Ly α emitters nearby one another if they reside in the same bubble. Interestingly, we do indeed find an instance where three of our massive z ≃ 7 Ly α emitters possess very similar redshifts (Δz ∼ 0.05) and are also closely separated in angular space (<5 arcmin).

Our spectroscopic results indicate that COS-940214 (⁠|${z}_{_{\mathrm{Ly\alpha }}}$| = 6.748), COS-1009842 (⁠|${z}_{_{\mathrm{Ly\alpha }}}$| = 6.761), and COS-995126 (⁠|${z}_{_{\mathrm{Ly\,\alpha }}}$| = 6.813) all lie within a spherical volume with radius R = 1.7 physical Mpc (Fig. 10), consistent with the expected sizes of H ii regions at z = 7 (e.g. Furlanetto et al. 2004; Mesinger & Furlanetto 2007; Lin et al. 2016). These three sources have luminosities ranging from 0.8–2.3 L|$^{\ast }_{\mathrm{UV}}$| and therefore are likely among the most massive galaxies present at z ≃ 7. All three are also inferred to be actively forming stars with SFRs of 7–28 M_⊙ yr⁻¹ using the SED fitting procedure described in Section 2.3. Their very blue IRAC colours ([3.6]−[4.5] |$\lt \, -$|0.6) further suggest large sSFRs ranging from approximately 10–130 Gyr⁻¹ (see Table 5). Because the typical sSFR at z ≃ 7 is 5 Gyr⁻¹ (E21), these three galaxies are likely undergoing a burst of star formation.

Figure 10.

Illustration of three closely separated Ly α emitters in our z ≃ 7 sample. All three sources fall within a spherical region with radius R = 1.7 physical Mpc, consistent with the expected sizes of H ii regions at z ≃ 7 (e.g. Lin et al. 2016). The proximity of these Ly α emitters thus may plausibly reflect the presence of a large ionized bubble.

Open in new tabDownload slide

We estimate that the three aforementioned sources could individually power bubbles of sizes R = 0.29–0.41 physical Mpc, consistent with other estimates recently reported at z ≃ 7–8 (Castellano et al. 2018; Tilvi et al. 2020). The largest region (R = 0.41 Mpc) is estimated to come from COS-955126 given that it is inferred to produce more than twice as many ionizing photons as the other two sources, mainly due to its much higher luminosity. COS-940214 and COS-1009842 have very similar estimated bubble sizes (R = 0.30 and 0.29 Mpc, respectively) due to the similarity in their inferred |$\dot{N}_{\mathrm{ion}}$|⁠. Notably, these estimated bubble sizes (⁠|$R\, \sim$| 0.3 Mpc) are much smaller than the observed separation between the three galaxies (R = 1.7 Mpc).

One way in which a larger bubble may have grown around these sources is through the impact of a local overdensity (e.g. Furlanetto et al. 2004; Wyithe & Loeb 2005; Lee et al. 2007; McQuinn et al. 2007). Such a scenario may already be hinted at by the proximity of these ≳L|$^{\ast }_{\mathrm{UV}}$| galaxies. Using the z ≃ 7 luminosity function from Bowler et al. (2017), we would expect only N = 0.7 bright (≥0.8 L|$^{\ast }_{\mathrm{UV}}$|⁠) galaxies on average within a 1.7 Mpc radius sphere. Because we know there are at least three such galaxies within this volume, our data suggest this region is likely a factor of ≳4 × overdense. Further evidence of an overdensity comes from investigating the surface density of photometrically selected z ≃ 7 galaxies over this region. In an ≈60 arcmin² rectangular area encompassing the three Ly α emitters, we identify a total of nine galaxies satisfying our z ≃ 7 selection criteria (Section 2.1), a factor of 6× above that expected on average¹⁹ in the same area (N = 1.5). Future Ly α observations will enable us to better quantify the spectroscopic overdensity within this region.

4 DISCUSSION

Numerous observational studies over the past decade have demonstrated a low (≲10–20 per cent) Ly α detection rate among typical (<L|$^{\ast }_{\mathrm{UV}}$|⁠) z ≳ 7 galaxies (e.g. Ono et al. 2012; Treu et al. 2013; Pentericci et al. 2014, 2018; Schenker et al. 2014; Tilvi et al. 2014). It was therefore a surprise when Ly α was detected in all four of the brightest (3–4 L|$^{\ast }_{\mathrm{UV}}$|⁠) z > 7 galaxies selected over the CANDELS fields (i.e. the RB16 sample), particularly given that the IGM is thought to be highly neutral at these epochs (x_H i > 40 per cent; e.g. Davies et al. 2018; Planck Collaboration VI 2020; Wang et al. 2020; Whitler et al. 2020). The unusual Ly α detectability of these sources hence suggested that their Ly α photons may be less sensitive to strong scattering by the IGM, as would be expected if these objects reside in large, highly ionized bubbles (e.g. Zitrin et al. 2015; Weinberger et al. 2018). However, it has been proposed that intense radiation fields may also be driving their enhanced Ly α detectability given that all four RB16 galaxies are inferred to possess strong [O iii]+H β emission (EW > 800 Å; Stark et al. 2017). In this section, we use our observations of a larger population (N = 22) of similarly bright (L|$^{}_{\mathrm{UV}}\, \simeq$| 1–6 L|$^{\ast }_{\mathrm{UV}}$|⁠) z ≃ 7 galaxies with Spitzer/IRAC constraints on [O iii]+H β EWs to better understand the origin of strong Ly α within the RB16 sample.

We find that the detectability of Ly α depends strongly on the [O iii]+H β EW at z ≃ 7. We detect Ly α in 78 per cent (7/9) of bright z ≃ 7 galaxies possessing strong [O iii]+H β (EW > 800 Å) as opposed to only 8 per cent (1/12) of galaxies with more moderate [O iii]+H β (EW = 200–800 Å; Section 3.1). The much (4.7×) higher Ly α EWs of the strong [O iii]+H β population likely reflects, in part, a larger ionizing photon production efficiency associated with young stellar populations (Chevallard et al. 2018; Tang et al. 2019). The high sSFRs (≳30 Gyr⁻¹) of these galaxies may also help drive low H i column density channels through the ISM/CGM (e.g. Clarke & Oey 2002; Ma et al. 2020) through which Ly α photons can efficiently escape (Jaskot et al. 2019; Gazagnes et al. 2020).

While we have shown that efficient ionizing photon production significantly boosted the Ly α detectability of the RB16 sample, our results also suggest that the IGM ionization may have played a significant role in enhancing their visibility as well. In particular, we find no evidence of strong evolution in the Ly α EW distribution of bright (L|$^{}_{\mathrm{UV}}\, \simeq$| 1–6 L|$^{\ast }_{\mathrm{UV}}$|⁠) galaxies between z ≃ 6–7 (Section 3.2). This result stands in contrast to the rapid downturn in the fraction of strong Ly α emitters (>25 Å EW) among low-mass (∼0.1 L|$^{\ast }_{\mathrm{UV}}$|⁠) lensed galaxies between z ≃ 6–7 (Hoag et al. 2019; Fuller et al. 2020), suggesting that Ly α in the bright population may evolve more slowly. Such findings are consistent with theoretical expectation that the brightest (and hence the most massive; e.g. Barone-Nugent et al. 2014) z ≃ 7 galaxies trace overdensities that are the first to create large, highly ionized bubbles (e.g. Furlanetto et al. 2004; Wyithe & Loeb 2005; Lee et al. 2007; McQuinn et al. 2007; Weinberger et al. 2018).

If massive z ≃ 7 galaxies do trace overdense regions, we may expect to find Ly α emitting galaxies in our sample nearby one another within the same ionized structure. Indeed, our results revealed three z ≃ 7 Ly α emitters separated by relatively small angular distances (<5 arcmin) and with similar redshifts (Section 3.3). The corresponding spatial separation of these Ly α emitters (R = 1.7 physical Mpc) is consistent with the expected sizes of H ii regions at z ≃ 7 (Furlanetto et al. 2004; Mesinger & Furlanetto 2007; Lin et al. 2016), suggesting that they may lie within the same large, highly ionized bubble. None the less, with our present data, it is not yet possible to characterize the extent of the ionized regions surrounding these z ≃ 7 Ly α emitters. There are two distinct possibilities consistent with our current data. On one hand, these massive z ≃ 7 galaxies may be tracing a large ionized region, with a size comparable to their spatial separation (R ≳ 2 Mpc). This would be expected if the galaxies trace a strongly overdense structure with larger-than-average ionizing photon output. In this case, we would expect to see enhanced Ly α from all systems within the ionized region. However, it is also possible that these systems are situated within distinct, moderate-sized (R ∼ 0.3 Mpc) bubbles in which they are the dominant contributors of ionizing photons. The smaller bubbles would lead to larger damping wing attenuation of Ly α from the IGM. But the corresponding reduction in the Ly α flux is countered by efficient production of nebular emission, as signaled by the large EW [O iii]+H β emission (and hence high sSFRs). The transmission could be further boosted if these massive systems have large velocity offsets (e.g. Erb et al. 2014; Stark et al. 2017; Mason et al. 2018b), redshifting the line further into the damping wing before it encounters hydrogen.

Which of these two pictures is true depends largely on whether the nearby Ly α emitting galaxies we have identified trace a strong overdensity in galaxies. In Section 3.3, we presented evidence that this may indeed be the case. The number density of spectroscopically confirmed >0.8 L|$^{\ast }_{\mathrm{UV}}$| systems in this region (N = 3) is ≈4 times the average (N = 0.7) expected from z ≃ 7 luminosity functions (e.g. Bowler et al. 2017). Furthermore, the surface density of photometrically selected z ≃ 7 galaxies surrounding the nearby Ly α emitters is 6× the average. This suggests that an overdense population of neighbouring galaxies may potentially help power a large ionized region.

These results add to two similar instances of grouped Ly α emitting galaxies that have been reported at z = 7.0–7.7 (Vanzella et al. 2011; Castellano et al. 2018; Tilvi et al. 2020). Both these previously identified regions contain at least one bright (>L|$^{\ast }_{\mathrm{UV}}$|⁠) galaxy. Furthermore, one of these regions (the BDF) shows evidence of a host overdensity similar to that in our sample (Castellano et al. 2016). Interestingly, deep spectroscopic follow-up of relatively faint (<0.7 L|$^{\ast }_{\mathrm{UV}}$|⁠) z ≃ 7 systems identified photometrically in the BDF revealed no Ly α detections (Castellano et al. 2018). As stated therein, this may suggest that many of the fainter systems reside outside the ionized region, though it is also possible that Ly α is preferentially seen in the brightest systems because of their larger velocity offsets. Ultimately, future Ly α spectroscopy of fainter sources (with known spatial positions via e.g. JWST rest-optical line detections) will be able to discern whether the Ly α EW distribution in the vicinity of these massive Ly α emitters is enhanced towards larger values as would be expected if they reside in large bubbles.

5 SUMMARY

We present the first results from a new spectroscopic MMT/Binospec program aimed at targeting Ly α in a large sample of bright (L|$^{}_{\mathrm{UV}}\, \simeq$| 1–6 L|$^{\ast }_{\mathrm{UV}}$|⁠) z ≃ 7 galaxies selected over very wide-area fields (∼3 deg² total). We use these results to investigate to what extent strong [O iii]+H β emission (and hence large sSFRs) boost Ly α detectability from bright reionization-era galaxies. We secondly test whether the decline in the visibility of Ly α emission is less rapid for massive galaxies, as might be expected if they trace strong overdensities that are situated in large ionized bubbles. Our conclusions are as follows:

• From our 22 targeted bright z ≃ 7 galaxies, we confidently (>7σ) detect Ly α from nine sources. The redshifts of our z ≃ 7 Ly α emitters range from |${z}_{_{\mathrm{Ly\alpha }}}$| = 6.650–7.093, consistent with expectations given our photometric selection criteria. We measure Ly α fluxes ranging from (5.1–18.6)× 10⁻¹⁸ erg s⁻¹ cm⁻² and rest-frame EWs spanning 3.7–43.1 Å. For sources lacking a confident detection, the data typically place a (7σ) Ly α EW limit of ≲10 Å in skyline-free regions of the spectra.

• We find that the detectability of Ly α depends strongly on [O iii]+H β EW among luminous z ≃ 7 galaxies. We confidently detect Ly α in 78 per cent (7/9) of sources with strong [O iii]+H β emission (EW > 800 Å) as opposed to only 8 per cent (1/12) of sources with more moderate [O iii]+H β (EW = 200–800 Å). The much (4.7×) higher Ly α EWs of the strong [O iii]+H β population likely reflect a larger ionizing photon production efficiency (and hence likely also a larger L|$_{\mathrm{Ly}\alpha}$|/L|$_{\text{UV}}$|⁠) owing to extremely young, hot stars formed in a recent rapid upturn or burst in star formation activity (Chevallard et al. 2018; Tang et al. 2019). The high sSFRs (≳30 Gyr⁻¹) of such galaxies may also help create low H i column density channels through the ISM/CGM (e.g. Clarke & Oey 2002; Ma et al. 2020) through which Ly α photons can efficiently escape (e.g. Jaskot et al. 2019; Gazagnes et al. 2020). Given the large variations in Ly α detectability within our data, it will be necessary to ensure that future samples be well matched in sSFR across cosmic time to well infer the evolution of the IGM neutral state.

• We tentatively detect the nebular NVλ1238.8,1242.8 doublet, a signpost of significant AGN activity, in one of our z ≃ 7 Ly α emitters. This is the fifth tentative detection of NV so far reported in a z ≳ 7 Ly α emitter (Tilvi et al. 2016; Hu et al. 2017; Laporte et al. 2017; Mainali et al. 2018) suggesting that low-luminosity AGN are present in a subset of the most massive reionization-era galaxies. Further investigation of this population will ultimately help clarify the contribution of such AGN to cosmic reionization (Madau & Haardt 2015; Giallongo et al. 2019; Grazian et al. 2020).

• We find no evidence for strong evolution in the transmission of Ly α emission for the massive (≃1–6 L|$^{\ast }_{\mathrm{UV}}$|⁠) population between z ≃ 6–7 (T_{z ≃ 7}/T_{z ≃ 6} = 1.04|$^{+0.39}_{-0.29}$|⁠). This is in contrast to observations of low-mass (∼0.1 L|$^{\ast }_{\mathrm{UV}}$|⁠) lensed galaxies which suggest a factor of ∼10 decline in transmission (Hoag et al. 2019; Fuller et al. 2020). With our current sample size, we can rule out a factor of ≥2 decline in Ly α transmission with 98.5 per cent (≈2σ) confidence among massive galaxies. We discuss a number of possible explanations for these findings, including the expectation that massive z ≃ 7 galaxies often reside within large, highly ionized bubbles (e.g. Furlanetto et al. 2004; McQuinn et al. 2007; Weinberger et al. 2018).

• We find three Ly α emitters in our z ≃ 7 sample separated by relatively small angular distances (<5 arcmin) and with nearly identical redshifts (Δz ≈ 0.05). These small angular separations may be a signpost of an ionized bubble enhancing Ly α transmission. Indeed, the spatial separations of these Ly α emitters (R = 1.7 physical Mpc) are consistent with the expected sizes of H ii regions at z ≃ 7 (e.g. Furlanetto et al. 2004; Mesinger & Furlanetto 2007; Lin et al. 2016). With our present data, we cannot yet characterize the full extent of the ionized regions surrounding these closely separated Ly α emitters. However, we estimate that these Ly α emitters are individually capable of powering ionized bubbles with radii of R∼0.3 Mpc, consistent with other estimates recently reported at z ≃ 7–8 (Castellano et al. 2018; Tilvi et al. 2020). Future work targeting fainter galaxies in the region should be able to determine if the ionized region extends beyond these radii.

• We find tentative evidence of an overdensity surrounding these closely separated Ly α emitters. The number density of spectroscopically confirmed >0.8 L|$^{\ast }_{\mathrm{UV}}$| systems in these regions (N = 3) is ≃4 times the average (N = 0.7) expected from z ≃ 7 luminosity functions. Furthermore, the surface density of photometrically selected z ≃ 7 galaxies surrounding the nearby Ly α emitters is 6× the average. Such an overdensity may help facilitate the growth of a large (R ≳ 2 physical Mpc) bubble around these sources.

As a next step towards understanding Ly α emission from massive reionziation-era galaxies, in future work we will combine this Binospec data set with recent results from a Cycle 7 ALMA Large Program, the Reionization Era Bright Emission Line Survey (REBELS; 2019.1.01634.L). REBELS has begun providing systemic redshifts via the [C ii]158 μm emission line as well as constraints on far-infrared dust continuum emission for a number of our z ≃ 7 targets. With this more complete rest-UV through far-infrared view, we will begin to better characterize the Ly α velocity offsets of massive z ≃ 7 galaxies as well as understand how their observed [C ii] and dust continuum emission relates to Ly α.

ACKNOWLEDGEMENTS

RE and DPS acknowledge funding from JWST/NIRCam contract to the University of Arizona, NAS5-02015. JC and SC acknowledge financial support from the European Research Council (ERC) via an Advanced Grant under grant agreement no. 321323 – NEOGAL. RJB and MS acknowledge support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek via TOP grant TOP1.16.057. BER was supported in part by NASA program HST-GO-14747, contract NNG16PJ25C, and grant 80NSSC18K0563, and NSF award 1828315. Observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution.

This research has benefited from the SpeX Prism Library (and/or SpeX Prism Library Analysis Toolkit), maintained by Adam Burgasser at http://www.browndwarfs.org/spexprism. This research also made use of astropy, a community-developed core python package for Astronomy (Astropy Collaboration 2013; Price-Whelan et al. 2018); matplotlib (Hunter 2007); numpy (Van Der Walt, Colbert & Varoquaux 2011); and scipy (Jones et al. 2001).

DATA AVAILABILITY

The optical through mid-infrared imaging data underlying this article are available through their respective data repositories. See https://hsc-release.mtk.nao.ac.jp/doc/ for HSC data, http://www.eso.org/rm/publicAccess#/dataReleases for UltraVISTA and VIDEO data, and https://sha.ipac.caltech.edu/applications/Spitzer/SHA/ for IRAC data. The MMT/Binospec data will be shared upon reasonable request to the corresponding author.

Footnotes

REFERENCES