KLASS z\sim 8

Inferences on the Timeline of Reionization at
From the KMOS Lens-Amplified Spectroscopic Survey

Charlotte A. Mason, Adriano Fontana, Tommaso Treu, Kasper B. Schmidt,Austin Hoag, Louis Abramson, Ricardo Amorin, Maruša Bradač,Lucia Guaita, Tucker Jones, Alaina Henry, Matthew A. Malkan,Laura Pentericci, Michele Trenti, and Eros Vanzella
Department of Physics and Astronomy, UCLA, Los Angeles, CA, 90095-1547, USA
Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA, 02138, USA
INAF Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monteporzio (RM), Italy
Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
Department of Physics, University of California, Davis, CA, 95616, USA
Carnegie Observatories, 813 Santa Barbara St., Pasadena, California 91101-1292, USA
Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile
Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile
Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218
School of Physics, University of Melbourne, Parkville, Victoria, Australia
ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)
INAF Osservatorio Astronomico di Bologna, via Gobetti 93/3, 40129 Bologna, Italy
Hubble FellowE-mail: charlotte.mason@cfa.harvard.edu
Accepted XXX. Received YYY; in original form ZZZ
Abstract

Detections and non-detections of Lyman alpha (Ly) emission from galaxies ( Gyr after the Big Bang) can be used to measure the timeline of cosmic reionization. Of key interest to measuring reionization’s mid-stages, but also increasing observational challenge, are observations at , where Ly redshifts to near infra-red wavelengths. Here we present a search for Ly emission in 53 intrinsically faint Lyman Break Galaxy candidates, gravitationally lensed by massive galaxy clusters, in the KMOS Lens-Amplified Spectroscopic Survey (KLASS). With integration times of hours, we detect no Ly emission with S/N in our sample. We determine our observations to be 80% complete for spatially and spectrally unresolved emission lines with integrated line flux  erg s cm. We define a photometrically selected sub-sample of 29 targets at , with a median Ly EW limit of 58 Å. We perform a Bayesian inference of the average intergalactic medium (IGM) neutral hydrogen fraction using their spectra. Our inference accounts for the wavelength sensitivity and incomplete redshift coverage of our observations, and the photometric redshift probability distribution of each target. These observations, combined with samples from the literature, enable us to place a lower limit on the average IGM neutral hydrogen fraction of at , providing further evidence of rapid reionization at . We show that this is consistent with reionization history models extending the galaxy luminosity function to , with low ionizing photon escape fractions, .

keywords:
dark ages, reionization, first stars – galaxies: high-redshift – galaxies: evolution – intergalactic medium
pubyear: 2019pagerange: Inferences on the Timeline of Reionization at From the KMOS Lens-Amplified Spectroscopic SurveyC.4

1 Introduction

The reionization of intergalactic hydrogen in the universe’s first billion years is likely linked to the formation of the first stars and galaxies: considered to be the primary producers of hydrogen-ionizing photons (e.g., Lehnert & Bremer, 2003; Bouwens et al., 2003; Yan & Windhorst, 2004; Bunker et al., 2004; Shull et al., 2012; Bouwens et al., 2015). Accurately measuring the timeline of reionization enables us to constrain properties of these first sources (e.g., Robertson et al., 2013; Robertson et al., 2015; Greig & Mesinger, 2017).

Measurements of the reionization timeline are challenging, however, due to the rarity of bright quasars at (Fan et al., 2001; Manti et al., 2016; Parsa et al., 2018), which have historically provided strong constraints on the end stages of reionization (e.g., Fan et al., 2006; McGreer et al., 2014; Greig et al., 2017; Bañados et al., 2017). In the coming decade 21 cm observations are expected to provide information about the IGM and the nature of the first galaxies (e.g., Liu & Parsons, 2016; Mirocha et al., 2016), but current progress has been driven by observations of Ly (rest-frame 1216 Å) emission in galaxies, using near infra-red (NIR) spectroscopy.

Ly is a highly resonant line, and strongly scattered by intervening neutral hydrogen as it travels to our telescopes. Whilst young star-forming galaxies, selected with a Lyman Break (Lyman Break Galaxies – LBGs) show Ly emission in abundance up to (e.g., Stark et al., 2011; Hayes et al., 2011; Curtis-Lake et al., 2012; Cassata et al., 2015; De Barros et al., 2017), at higher redshifts the fraction of galaxies detected with Ly emission, and the scale length of the Ly rest-frame equivalent width (EW) distribution, decreases rapidly (e.g., Fontana et al., 2010; Pentericci et al., 2011; Caruana et al., 2012; Treu et al., 2012, 2013; Ono et al., 2012; Caruana et al., 2014; Pentericci et al., 2014; Schenker et al., 2014; Tilvi et al., 2014; Faisst et al., 2014; Jung et al., 2018). This rapid decline of detected Ly emission is most plausibly due to absorption in an increasingly neutral IGM (Dijkstra et al., 2011; Dijkstra, 2014; Mesinger et al., 2015).

Large spectroscopic surveys of LBG candidates are being assembled out to (Pentericci et al., 2011, 2014, 2018; Hoag et al., 2019) but exploring the earliest stages of reionization requires us to observe Ly at even higher redshifts. Only a handful of Ly emitters have been confirmed at (Zitrin et al., 2015; Oesch et al., 2015; Roberts-Borsani et al., 2016; Stark et al., 2017; Hoag et al., 2017), where the dominance of sky emission in the NIR makes observations of faint sources even more challenging. Additionally, because Ly emission can be spatially extended and/or offset from the UV continuum emission (Wisotzki et al., 2016; Leclercq et al., 2017), it is likely that slit-based spectroscopy is not capturing the full Ly flux. Hence, the observed decline in Ly emission at could be partially due to redshift-dependent slit-losses as well as reionization.

In this paper we present a search for Ly emission in NIR spectroscopy of 53 intrinsically faint LBG candidates (), gravitationally lensed behind 6 massive galaxy clusters, including 4 of the Frontier Fields (Lotz et al., 2017), selected from the Grism Lens-Amplified Survey from Space (hereafter GLASS, Schmidt et al., 2014a; Treu et al., 2015). We also present observations of Civ emission in 3 images of a previously confirmed multiply-imaged galaxy (Boone et al., 2013; Balestra et al., 2013; Monna et al., 2014).

The observations presented in this work were carried out with the ESO Very Large Telescope (VLT) K-band Multi Object Spectrometer (hereafter KMOS, Sharples et al., 2013). This work presents the first results of observations with KMOS. KMOS is an integral field unit (IFU) instrument, and we demonstrate here that our observations are more complete to spatially extended and/or offset Ly emission than traditional slit spectrographs.

We use our new deep spectroscopic observations to infer the average IGM neutral hydrogen fraction () at Mason et al. (2018a, hereafter M18a) presented a flexible Bayesian framework to directly infer from detections and non-detections of Ly emission from LBGs. The framework combines realistic inhomogeneous reionization simulations and models of galaxy properties. That work measured ( confidence intervals) at . Building on Treu et al. (2012) and M18a we extend this framework to use the full spectra obtained in our observations for the Bayesian inference, accounting for the incomplete wavelength coverage and spectral variation of the noise, and marginalising over emission linewidth. Our framework uses the photometric redshift probability distribution, obtained from deep photometry including new Spitzer/IRAC data, of each object to robustly account for uncertainties in redshift determination.

The paper is structured as follows: Section 2 describes our KMOS observations and the target selection from the GLASS parent sample; Section 3 describes the search for Ly emission in our KMOS data cubes, and the purity and completeness of our survey; and Section 4 describes the Bayesian inference of the neutral fraction and presents our limit on at . We discuss our findings in Section 5, including an assessment of the performance of KMOS for background-limited observations using our deep observations, and summarise in Section 6.

We use the Planck Collaboration et al. (2015) cosmology where (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s Mpc). All magnitudes are given in the AB system.

2 Observations

2.1 The KMOS Lens-Amplified Spectroscopic Survey

KLASS is an ESO VLT KMOS Large Program (196.A-0778, PI: A. Fontana) which targeted the fields of six massive galaxy clusters: Abell 2744 (hereafter A2744); MACS J0416.1-2403 (M0416); MACS J1149.6+2223 (M1149); MACS J2129.4-0741 (M2129); RXC J1347.5-1145 (RXJ1347); and RXC J2248.7-4431 (RXJ2248, aka Abell S1063). A2744, M0416, M1149 and RXJ2248 are all Frontier Fields (hereafter HFF, Lotz et al., 2017). Observations were carried out in Service Mode during Periods (October 2015 - October 2017).

KMOS is a multi-object IFU spectrograph, with 24 movable IFUs, split between 3 different spectrographs (Sharples et al., 2013). Each IFU is field of view, with pixel size , and 2048 pixels along the wavelength axis111We use the following definitions for describing 3D spectra in this paper. Pixel: 2D spatial pixel (size ). Spaxel: the 1D spectrum in a single spatial pixel (spanning the spectral range m, in 2048 spectral pixels). Voxel: 3D pixel in the data cube with both spatial and spectral indices..

The key science drivers of KLASS are:

  1. To probe the internal kinematics of galaxies at , with superior spatial resolution compared to surveys in blank fields (Mason et al., 2017).

  2. To investigate Ly emission from the GLASS sample, independently of the HST spectroscopic observations, providing validation and cross-calibration of the results and enabling us to constrain the timeline and topology of reionization (Treu et al., 2012, 2013; Schmidt et al., 2016; Mason et al., 2018a).

Mason et al. (2017) addressed the first science driver by presenting spatially resolved kinematics in 4 of the 6 KLASS clusters from our early data, including five of the lowest mass galaxies with IFU kinematics at , and provided evidence of mass-dependent disk settling at high redshift (Simons et al., 2017). The KLASS kinematic data were combined with metallicity gradients from the HST GLASS data to enable the study of metallicity gradients as a diagnostic of gas inflows and outflows (Wang et al., 2016).

This paper addresses the second science driver by presenting our candidate targets with complete exposures. We use the YJ observing band, giving us access to Ly emission at .

The choice of an IFU instrument for high-redshift Ly observations was motivated by indications that ground-based slit-spectroscopy measures lower Ly flux than HST slit-less grism spectroscopy (Tilvi et al., 2016; Huang et al., 2016b; Hoag et al., 2017), which, as well as reionization, could contribute to the observed decline in Ly emission at . Ly emission can be spatially extended and/or offset from the UV continuum emission (Feldmeier et al., 2013; Momose et al., 2014; Wisotzki et al., 2016; Leclercq et al., 2017), so it is likely that slit-based spectrographs do not capture the full Ly flux.

By using IFUs our observations should be more complete to spatially extended and/or offset Ly than traditional slit spectrographs. Mason et al. (2017) showed that only of emission line flux was contained in simulated slits (a typical slit-width used for Ly observations, e.g., Hoag et al., 2017) on KMOS spectra, whereas the full flux is captured within the KMOS field of view. Thus we expect most Ly flux to be captured within the KMOS IFUs. The wide IFUs cover proper kpc at , while the UV effective radii of galaxies at these redshifts is only proper kpc (Shibuya et al., 2015). We demonstrate in Section 3.3 that our KMOS observations have good completeness for spatially extended and/or offset Ly emission.

2.2 Target selection

Cluster Run ID DIT [s] NDITs Exposure [hrs] Number of targets
Category 1 Category 2
A2744 A 900 25 6.25 3 7
M0416 B 900 43 10.75 2 5
M1149 C 900 40 10.00 2 6
M2129 E 450 85 10.625 3 7
RXJ1347 D 450 88 11.00 3 8
RXJ2248 F 300 93 7.75 1 6
Note. The number of Detector Integration Times (DITs) used in this analysis: we discarded DITs if the seeing was as measured by stars observed in each DIT. The total exposure time = DIT NDITs. We had to discard our initial 4 hours of observations of A2744 due to irreparable flexure issues due to rotating the instrument between science and sky DITs. All subsequent observations were performed with no rotation of the instrument between science and sky DITs. Target selection in these clusters was primarily done from preliminary versions of the ASTRODEEP catalogues (Castellano et al., 2016; Merlin et al., 2016; Di Criscienzo et al., 2017), which did not include Spitzer/IRAC photometry. Due to a high proper motion reference star, some of the observations of RXJ2248 were taken at a slight offset from the required target centre, reducing the total exposure at that position. RXJ2248 also included 3 targets (Appendix A).
Table 1: KLASS cluster targets

KLASS targets were selected from the GLASS survey222http://glass.astro.ucla.edu (Schmidt et al., 2014a; Treu et al., 2015), a large Hubble Space Telescope (HST) slit-less grism spectroscopy program. GLASS obtained spectroscopy of the fields of 10 massive galaxy clusters, including the HFF and 8 CLASH clusters (Postman et al., 2012). The Wide Field Camera 3 (WFC3) grisms G102 and G141 were used to cover the wavelength range m with spectral resolution . We refer the reader to Schmidt et al. (2014a) and Treu et al. (2015) for full details of GLASS.

KLASS observations aimed to provide the high spectral resolution necessary to measure the purity and completeness of the grism spectra, to measure lines that were unresolved in HST, and to obtain velocity information for targets which the low resolution grisms cannot provide. In combination with additional GLASS follow-up observations at Keck (Huang et al., 2016b; Hoag et al., 2017; Hoag et al., 2019) we will address the purity and completeness of the HST grisms in a future work. In this work we present our high-redshift candidate targets and our inferences about reionization obtained from the KLASS data.

Two categories of high-redshift candidate KLASS targets were selected from the GLASS data:

  1. Category 1: 14 objects with marginal (S/N ) candidate Ly emission in the HST grisms, identified by visual inspection of the GLASS data, which fall within the KMOS YJ spectral coverage (m). 4 candidates were selected from a list of candidates in a preliminary census of GLASS data by Schmidt et al. (2016). The remaining candidates were selected in a similar method to the procedure followed by Schmidt et al. (2016).

  2. Category 2: 39 LBG candidates selected with , from an ensemble of photometric catalogues described by Schmidt et al. (2016). This includes three LBGs which were spectroscopically confirmed via sub-mm emission lines after our survey began: A2744_YD4 (A2744_2248 in this paper), at (Laporte et al., 2017a, discussed in Section 4.2 and 4.3), M0416_Y1 (M0416_99 in this paper), at (Tamura et al., 2018, discussed in Section 4.2), and M1149_JD1, at (Hashimoto et al., 2018, discussed in Section 5.2).

An additional three targets were multiple images of the system in RXJ2248 (Boone et al., 2013; Balestra et al., 2013; Monna et al., 2014) where we targeted Civ1548,1551 emission. This object is discussed in Appendix A.

We ranked objects in order of the number of inspectors who reported a candidate emission line for our Category 1 targets, and then by the number of independent photometric catalogues the target appeared in (for both categories). Our observations were planned prior to the release of the full HFF datasets, so the photometric catalogues we used to select candidates did not contain the full photometry now available. In particular, deep Spitzer/IRAC data did not exist, which can be useful for distinguishing between high-redshift star forming galaxies and passive galaxies. Nor were sophisticated intra-cluster light (ICL) removal techniques developed at that point (e.g., Merlin et al., 2016; Morishita et al., 2017; Livermore et al., 2017).

Thus our LBG selection was heterogeneous, but in this paper we now add in the new deep and extended photometry to define a homogeneous photometric selection. We expect some faint candidates may have been spurious in the initial photometry and may not appear in the final deep catalogues. Additionally we expect that with the inclusion of Spitzer/IRAC photometry some of the objects originally selected to be may be low redshift contaminants. In our reionization analysis we use catalogues built using the final HFF datasets to define a selection function for a photometrically-selected sample for our inference (described in Section 4.2). We demonstrate that this KLASS sub-sample is not a biased sample of the final parent catalogues in Appendix B.

The GLASS median flux limit is  erg s cm (Schmidt et al., 2016), and we tried to be as inclusive as possible when assigning candidates to the KMOS IFUs from the GLASS parent catalogue. Most of the candidates were only significance in GLASS data and we were aiming to provide confirmation of those tentative targets our deep KMOS observations - though for our ground-based observations, at least 50% of the wavelength range is dominated by sky emission and the low spectral resolution () of the HST grisms means that the line position is uncertain by  Å so the lines could be in bad sky regions. Additionally, in planning our observations we likely overestimated the sensitivity of KMOS YJ using the online exposure time calculator, especially at the blue end of the detectors. We discuss this in more detail in Section 5.4.

As we describe below in Section 3.4, our KLASS observations are complete for lines with flux  erg s cm, which suggests we should have confirmed the majority of the 14 GLASS candidate Ly emission targets we observed (Category 1, described in Section 2.2). However, we did not detect any emission in the cubes containing these candidates, suggesting at least some of the GLASS candidates were spurious noise peaks in the HST grisms. A more thorough comparison of the GLASS HST grism and ground-based follow-up observations (including KLASS and Keck observations, Huang et al., 2016b; Hoag et al., 2017; Hoag et al., 2019) to recover the grism purity and completeness will be left to a future work.

53 candidate targets across the 6 clusters were assigned to 51 KMOS IFUs (two IFUs contained two nearby candidates). The cluster list and number of high-redshift candidate targets per cluster is shown in Table 1.

2.3 KLASS observing strategy and reduction

KLASS observations were carried out with KMOS YJ (m). The spectral resolution is sufficient to distinguish Ly from potential low redshift contaminants with the [Oii] emission doublet at .

Observations were carried out in service mode and executed in one hour observing blocks with repeating ABA science-sky-science integration units (detector integration times – DITs). Each observing block comprised 1800 s of science integration, and 900 s on sky. Pixel dither shifts were included between science frames. A star was observed in 1 IFU in every observing block to monitor the point spread function (PSF) and the accuracy of dither offsets. The PSF was well-described by a circular Gaussian and the median seeing of our observations was FWHM .

In each cluster, the 3 top priority targets were observed for the average exposure time by assigning 2 IFUs per target and nodding between them during A and B modes.

2.4 Reduction

Data were reduced using the ESO KMOS pipeline v.1.4.3 (Davies et al., 2013). We apply a correction for known readout channel level offsets before running the pipeline. We run the pipeline including optimised sky subtraction routines sky_tweak (Davies, 2007) and sky_scale.

To improve the sky subtraction in the pipeline-reduced ‘A-B’ cubes we produced master sky residual spectra by median combining all IFUs on each spectrograph into a 1D master sky residual spectrum for each DIT, excluding cubes containing targets with bright emission lines and/or continua. We then subtract these 1D sky residual spectra from the ‘A-B’ cubes on the same spectrograph for each DIT, rescaling the 1D spectra in each spaxel so that the resulting 1D spectrum in each spaxel is minimised.

Similar techniques to improve sky subtraction are described by Stott et al. (2016). This method worked best for our survey design. We note that this method performed better than in-IFU sky residual subtraction (i.e. subtracting a median sky residual spectrum produced from ‘empty’ spaxels in each IFU) as it preserved emission line flux in the modestly sized KMOS IFUs.

Cube frames from each DIT are combined via sigma clipping, using spatial shifts determined by the position of the star observed in the same IFU in each DIT, to produce flux and noise cubes. For this work we used only frames with seeing (as measured by the star observed in our science frames). The median seeing was . DIT length, observing pattern and total integration times used for this paper are listed in Table 1. We note that due to the failure of one of the KMOS arms, no star was observed in the A2744 observations. We used a bright target to estimate the dither offsets for this cluster.

For pure Gaussian noise, the pixel distribution of S/N should be normally distributed. We tested this by selecting non-central regions of cubes containing high-redshift candidate targets (i.e. where we expect very little source flux) and found the pixel distribution of S/N to have standard deviation , suggesting the noise is underestimated by the pipeline.

We therefore apply a rescaling to the noise of the combined cubes. We create an average 1D noise spectrum in a single spaxel for each cluster by taking the root-mean-square (RMS) at every wavelength of every spaxel from the cubes containing high-redshift candidate targets. Since the cubes are predominantly noise, taking the RMS of the flux at each wavelength across multiple cubes should give the appropriate noise. We find this RMS spectrum is higher than the pipeline average 1D noise spectrum (taking the average of the noise cubes across the same set of high-redshift targets). We rescale the pipeline noise in every cube by this ratio of the cluster RMS noise spectrum to the cluster average noise spectrum.

Finally, we rescale the noise in each cube by a constant value so that the S/N distribution of all pixels has standard deviation 1 (clipping pixels within 99.9% to remove spurious peaks). We find the S/N distribution is well-described by a Gaussian distribution, with non-Gaussian tails only beyond the confidence regions, due to bad sky subtraction residuals.

3 Emission Line Search, Purity and Completeness

In this section we describe our search for Ly emission in our KMOS observations. We give our algorithm for line detection in Section 3.1, and calculate the purity and completeness of our observations in Sections 3.2 and 3.3. Given that we detect no convincing Ly emission lines in our sample we present our flux and EW upper limits in Section 3.4.

3.1 Emission line detection technique

To search for emission lines in the KMOS cubes, to robustly determine the completeness and purity of our survey, and determine the flux limits of our observations, we used the following algorithm to flag potential lines:

  1. Create a circular aperture with pixels (using our median seeing ), which will capture of the total flux for spatially unresolved emission line at the centre of the aperture.

  2. Sum the flux, and take the RMS of the noise of all spaxels in the aperture to create 1D data and noise spectra.

  3. Rescale the 1D noise spectrum so the S/N in all pixels (excluding the 0.1% most extreme S/N values) is Normal.

  4. Scan through in wavelength and flag a detection if 3 adjacent wavelength pixels have . This corresponds to a detection of the integrated line flux.

  5. Iterate over 25 apertures centred within 3 pixels () of the IFU centre, i.e. , where is the IFU centre.

Our search covers potential emission line positions in each cube. As our detection threshold is we would expect a false positive rate of , i.e. false detections per cube for Gaussian noise. As discussed in Section 2.4 the S/N has small non-Gaussian tails due to sky subtraction residuals so we expect a slightly higher false detection rate than this.

3.2 Candidate emission lines and sample purity

We ran the detection algorithm described in Section 3.1 on the 54 cubes containing our high-redshift candidate targets (including the 3 cubes containing the images). 9 unique candidate lines were flagged (combining candidates at the same wavelength identified in different apertures). Each of these candidate lines was then visually inspected to determine whether it was a true emission line or a spurious noise peak. For our inspections we use both 1D spectra extracted in the detection apertures as well as 2D collapsed images of the candidate line obtained by summing cube voxels in the wavelength direction. The 2D images are helpful for determining plausible spatially compact emission from the uniform emission produced by sky residuals.

Our algorithm correctly identifies the Civ emission at 11023.7 Å in the brightest image of the multiply-imaged system, demonstrating the depth of our KMOS observations and the fidelity of our algorithm. Another detection is flagged in this object at 13358.6 Å but the emission appears diffuse and the wavelength is not consistent with other expected UV emission lines so we deem this spurious. We describe this object in more detail in Appendix A.

Of the remaining 7 lines flagged, 6 are deemed to be spurious detections as they are at the spectral edges of the detector, or immediately adjacent to strong skylines and appear to have P-Cygni profiles, indicating extreme sky subtraction failures. Whilst it could be possible to add a cut to e.g. downweight flagged lines adjacent to skylines, given the relatively low spectral resolution of our observation () we were wary that many true emission lines could be overlapping with skylines, thus visual inspection was necessary. This is clearly demonstrated in our detection of Civ emission where both doublet components overlap with sky lines (see Figure 9).

The remaining candidate emission line at 12683.7 Å is spatially offset from the LBG candidate in the cube. We determine the detected emission to be associated with a nearby () galaxy with , which has bright continuum emission in the GLASS data. The candidate line appears in a particularly bad spectral region of telluric absorption, and we determine the detection to be due to inadequate continuum subtraction of the source.

In our reductions we subtract a sky residual spectrum to minimise the flux in each spaxel of the high-redshift candidate cubes (Section 2.4). During that process most of the continuum emission from the object was poorly subtracted by scaling the sky residual spectrum to high values. Some residual flux is left, which correlates with the positions of sky residuals. We note that the LBG candidate targeted in this IFU is not present in the final deep photometric catalogues and is excluded from our reionization (it was likely a spurious detection in the original shallow photometry, Section 4.2). We remove this cube from further analysis.

Thus we determine our algorithm has detected 1 real emission line, and 7 spurious detections (excluding the continuum object described above), allowing us to define the purity of our spectral sample:

(1)

where is the total number of spurious flags (8 unique false detections which were sometimes flagged in multiple apertures) and is the number of possible emission line positions in the 53 useful cubes, removing wavelength pixels not covered by certain detectors, in 25 apertures. We measure . Our spurious detection rate is higher than that expected for fluctuations in the noise, which was expected due to the non-Gaussian tail in our S/N distribution due to sky subtraction residuals. To verify that the S/N distribution is symmetrical we also ran the detection algorithm to look for negative peaks (S/N) which should occur at the same rate. We found 12 flagged negative S/N detections, comparable to our 7 flagged spurious detections with positive S/N.

We ran the algorithm on our Category 1 sources with a lower S/N threshold: S/N per wavelength pixel, corresponding to in the integrated line. We found no convincing detections with this lower threshold and are thus unable to confirm any of the candidate GLASS emission lines. Given that most of the GLASS Ly candidates were of low significance in the GLASS HST data these candidates may have been spurious noise peaks in the grism data.

In Section 4.2 below we list the Ly flux and EW limits for our most likely LBGs candidates. We discuss our limits on other UV lines in Section 5.3.

3.3 Completeness

Figure 1: Completeness as a function of line flux (top left), spatial offset from IFU centre (top right), spectral linewidth (lower left) and spatial extent (lower right). Each colour corresponds to a separate cluster target. Dashed lines show the completeness across the entire wavelength range, solid lines show the completeness in wavelength regions where the noise level is below the median. FWHM velocities were calculated assuming . Spatial extent, FWHM is the extent of the source (excluding the PSF). We create simulated lines with total spatial extent . In each plot the parameter of interest and wavelength are varied, while the other parameters are held constant. The fiducial parameters are: line flux =  erg s cm, observed FWHM Å (i.e. unresolved), line centred at the IFU centre, and FWHM (i.e. unresolved point source).

To evaluate the completeness of our emission line search we carry out comprehensive Monte Carlo simulations: inserting simulated lines into cubes with varied total flux, spectral FWHM, spatial position, spatial extent FWHM, and wavelength, and testing whether they are detected by our detection algorithm (Section 3.1). Traditionally, these types of simulations are carried out by inserting simulated lines into real raw data and then running through the full reduction pipeline (Fontana et al., 2010; Pentericci et al., 2014; De Barros et al., 2017), however, due to the complexity of the KMOS pipeline which constructs 3D cubes from 2D frames we instead create simulated cubes and add Gaussian noise drawn from an average noise cube for each cluster, mimicking completeness simulations traditionally done in imaging.

We create simulated flux cubes with a 3D Gaussian emission line with varied properties and add noise to each voxel drawn from a Gaussian distribution with mean zero and standard deviation for each cluster. The cubes are constructed by taking the RMS at every voxel of all the final sky-subtracted cubes which do not contain bright sources ( cubes per cluster). As each ‘empty’ cube is expected to be pure noise, taking the RMS at each voxel across the cubes should give an estimate of the noise per voxel, .

We calculate completeness as a function of flux, spatial offset from the IFU centre, spectral linewidth and spatial extent. For each simulation we vary the parameter of interest and wavelength, and fix the other three parameters. Our fiducial values for the parameters are: line flux =  erg s cm, observed line FWHM Å (the spectral resolution, i.e. unresolved lines), line centred at the IFU centre, with source spatial extent FWHM (i.e. unresolved point source, the emission will have observed spatial extent with ). We draw 1000 realizations of an emission line with noise at every tested value of a parameter. The resulting completeness is the fraction of these simulated lines detected by our detection algorithm.

Our fiducial simulations assume Ly emission will be spatially unresolved. These assumptions are reasonable for the intrinsically UV faint LBGs we are observing (Schmidt et al., 2016; Marchi et al., 2018). Typical slit spectrograph observations of Ly emission centre slits on the UV continuum and use slit-widths , thus in KLASS we are more complete to Ly emission that may be spatially extended and/or offset from the UV continuum.

Figure 1 shows the results of our completeness simulations for all clusters. We reach completeness over the full wavelength range for lines  erg s cm, centred within of the IFU centre and with intrinsic line FWHM km s, assuming to calculate FWHM (median over all clusters). For wavelength ranges where the noise level is below the median across the whole spectrum, we reach completeness for lines  erg s cm, centred with of the IFU centre and with intrinsic line FWHM km s. The completeness is fairly flat for Ly spatial extent (total extent ) demonstrating our good completeness for spatially extended Ly emission, with the normalisation of the completeness as a function of FWHM scaling with the completeness at a given total line flux.

3.4 Flux and equivalent width limits

Figure 2: Average flux limits for each cluster as a function of wavelength, assuming emission lines are spatially unresolved. We use the 1D RMS noise spectrum for each cluster as described in Section 3.4 to obtain the flux limits. Each plot corresponds to the median across all IFUs containing high-redshift candidates, for the different cluster targets. The dashed horizontal lines mark the median flux limit for each cluster.
Figure 3: rest-frame Ly EW limits in RXJ1347 as a function of wavelength, for 3 values of UV apparent magnitude , assuming emission lines are spatially unresolved. We use the 5 flux limit for RXJ1347 shown in Figure 2 and divide by the continuum flux and at each wavelength to obtain the EW limit.

To calculate average flux limits for each cluster we take the average 3D noise spectrum for each cluster, (created by taking the RMS at every voxel across the IFUs observing high redshift candidates in each cluster). We then create a 1D noise spectrum, , by summing the average noise at each wavelength pixel in a circular aperture with radius (where we use our median seeing ). At each wavelength pixel, the flux limit in erg s cm is given by:

(2)

Here we obtaining an estimate of the integrated noise for an emission line with observed  Å or  km s (the instrumental resolution), and use a threshold integrated . The term in the denominator accounts for the fact that the apertures only capture a fraction of the flux. For this results in a rescaling of 1.16. The spectral pixel width of KMOS YJ is  Å. The above calculation assumes the emission is spatially and spectrally unresolved by KMOS, which is reasonable given the expectation that Ly emission from UV faint galaxies is likely to be more spatially compact and have lower linewidth than Ly from UV bright galaxies (e.g., Schmidt et al., 2016; Marchi et al., 2018). We note that the flux limit for wider lines can be estimated as .

The flux limits for all clusters are shown in Figure 2. The median flux limit is  erg s cm, and the range of medians for each cluster is  erg s cm.

Rest-frame Ly equivalent widths are , where (with  Å), and we define the continuum flux:

(3)

where erg s Hz cm, is the apparent magnitude of the UV continuum, is the speed of light, is the rest-frame wavelength of the UV continuum (usually 1500 Å), and is the UV slope. We assume , consistent with observations (e.g., Stanway et al., 2005; Blanc et al., 2011; Wilkins et al., 2011; Castellano et al., 2012; Bouwens et al., 2012, 2014). We use the magnitude measured in HST WFC3/F160W for the apparent magnitude (automag). Example EW limits for objects with a given apparent magnitude, using the RXJ1347 average flux limit, are plotted in Figure 3.

4 Reionization inference

In this section we describe the extension to the M18a Bayesian inference framework to include the full spectra, robustly including the uncertainties in redshift via the photometric redshift distribution (Section 4.1), and marginalising over the linewidth of potential emission lines. Using the observations described above we now define a clear selection function for a photometrically-selected sample of LBGs within our survey (Section 4.2), and perform the inference of the IGM neutral fraction using these data (Section 4.3).

4.1 Bayesian inference framework

To use our observations to make inferences about the neutral hydrogen fraction at we use the method described by M18a. This forward-models the observed rest-frame Ly EW distribution as a function of the neutral fraction and galaxy UV magnitude, , using a combination of reionization simulations with realistic inhomogeneous IGM structure (Mesinger et al., 2016), and empirical and semi-analytic models of galaxy properties.

The models assume the observed Ly EW distribution is the ‘emitted’ distribution (i.e. the distribution without IGM attenuation due to reionization) and use that to forward-model the observed distribution, including the impact of Ly velocity offsets. Here, as in M18a, we use the recent comprehensive Ly EW observations from De Barros et al. (2017). We use the public Evolution of 21cm Structure (EoS) suite of reionization simulations described by Mesinger et al. (2015); Mesinger et al. (2016)333http://homepage.sns.it/mesinger/EOS.html to generate Ly optical depths along millions of sightlines in simulated IGM cubes for a grid of volume-averaged values. As the size of ionised regions during reionization is expected to be nearly independent of redshift at fixed (as there is little difference in the matter power spectrum from , McQuinn et al., 2007), we use the same cubes as used by M18a rather than generating new cubes.

We refer the reader to M18a for more details of the forward-modelling approach. Here we describe the modifications we have made to our Bayesian inference to make use of the spectral coverage and sensitivity of our observations. We account for the incomplete redshift coverage and for the gravitational lensing magnification of the objects by the foreground clusters. We marginalise over a range of potential linewidths for the Ly emission lines. We also marginalise over the photometric redshift distribution for each galaxy, which we obtain from comprehensive photometry (Section 4.2), to robustly account for uncertainties and degeneracies in redshift determination.

We want to obtain the posterior distribution for the neutral fraction: for each galaxy, where is an observed flux density spectrum as a function of wavelength, is the observed apparent UV magnitude, and is the magnification. A full derivation of the posterior is shown in Appendix C, and we summarise it here.

Our inference framework calculates the likelihood of an emission line emitted at redshift with observed rest-frame EW, , being present in an observed flux density spectrum. To calculate this likelihood we must assume a lineshape for the observed emission line. In previous inferences by Treu et al. (2013); Pentericci et al. (2014) and Tilvi et al. (2014) treated emission lines as unresolved: lines were modelled as Dirac Delta functions, with all the flux contained in a single spectral pixel. However, motivated by recent observations of Ly emission with linewidth  km s, several times greater than the instrumental resolution (Ono et al., 2012; Finkelstein et al., 2013; Vanzella et al., 2011; Oesch et al., 2015; Zitrin et al., 2015; Song et al., 2016b), here we improve the method by including the effect of linewidth.

The inference is quite sensitive to linewidth as at fixed EW a broader line will have lower S/N in our observations. By assuming unresolved emission lines the lower limits on the reionization ‘patchiness’ parameter inferred by Treu et al. (2013); Pentericci et al. (2014) and Tilvi et al. (2014) will be slightly overestimated compared to a more realistic treatment of linewidth. We note that the neutral fraction inference by (M18a) used EW limits calculated assuming a range of realistic linewidths so their result does not need revision. We discuss the impact of linewidth in more detail in Appendix C.4 but note that our results are robust for FWHM in a realistic range  km s.

To modify our inference to account for linewidth, we assume Gaussian emission lines for simplicity so can write the model emission line flux density as a function of EW:

(4)

where , with  Å, is the redshift of an emission line, is the rest-frame equivalent width of the emission line, is the flux density of the continuum calculated using Equation 3 using the observed continuum apparent magnitude , and is the spectral linewidth.

The likelihood of observing a 1D flux density spectrum for an individual galaxy (where is the wavelength pixel index), given our model where the true EW is drawn from the conditional probability distribution is:

(5)

where is the uncertainty in flux density at wavelength pixel and there are a total of wavelength pixels in the spectrum. is the probability distribution for the observed rest-frame EW as a function of the neutral fraction, and galaxy properties – UV apparent magnitude, magnification, and redshift. This PDF is obtained by convolving the model outputs from M18a with the probability distribution for each galaxy’s absolute UV magnitude, including errors on and (Equation 13).

We note that the range of neutral fraction in the EoS simulations is . In order to correctly calculate posteriors and confidence intervals we set the likelihood at such that we expect to observe no Ly flux at all (in a fully neutral universe). I.e. .

Given our relatively small sample size, we choose to restrict our inference to , thus for ease of computation we evaluate at ; this has a negligible impact on the final likelihood. We keep free in the rest of the inference. This product of likelihoods over the wavelength range of the spectrum accounts for the wavelength sensitivity of our observations, i.e. high noise regions are weighted lower than low noise regions.

We also note that EW is independent of magnification. Therefore, our inferences should be quite robust to magnification, which enters only through the dependency on of the assumed intrinsic EW distribution.

Using Bayes’ theorem, the posterior distribution for , and FWHM is:

(6)

We use a uniform prior on between 0 and 1, , and use the photometric redshift distribution for the prior . As we are only interested in the posterior probability of we can marginalise over FWHM and for each galaxy. We use a log-normal prior on FWHM with mean depending on derived through empirical relations and 0.3 dex width; we discuss our choice of FWHM priors in more detail in Appendix C.4 but find our results to be negligibly changed if we had used a uniform prior spanning the range of observed Ly FWHM at ( km s). To account for the incomplete wavelength coverage, we use the fact that if the object has Ly outside of the KMOS wavelength range (covering , ) we would measure a non-detection in our data. Thus the posterior for from one galaxy is:

(7)

We assume all galaxies observed are independent, so that the final posterior is the product of the normalised posteriors (Equation 7) for each object.

Using the photometric redshift distributions as a prior on the redshift allows us to incorporate the probability of each galaxy truly being at high redshift (rather than a low redshift contaminant) in a statistically rigorous way. In combining the posteriors in Equation 7 for each galaxy, the photometric redshift distribution weights the individual posteriors based on the probability of the source being within our redshift range. LBGs usually have degeneracies in their photometry which make it difficult to determine whether they are high redshift star-forming galaxies or mature galaxies. Thus with our method we are able to obtain reionization inferences from sources even when the photometric redshift distribution has multiple and/or broad peaks.

Whilst here we have carried out the inference at only, with larger samples, it will be possible to measure directly, for example by parametrising its evolution with redshift and inferring the values of its redshift-dependent parameters, or in a Markov Chain Monte Carlo exploration of IGM simulations to also infer relevant astrophysical parameters (Greig & Mesinger, 2015, 2017).

4.2 Defining a selection function for a photometric sample

Figure 4: Photometric redshift distributions centred on the KMOS observable range. We show the KMOS YJ range for Ly with the solid blue horizontal line. Black lines show the for the 29 sources which have probability of (marked by blue dashed vertical lines) which we use for the inference. 14 sources have , including the galaxy M1149_JD1, recently spectroscopically confirmed at with ALMA by Hashimoto et al. (2018). In our photometric catalogue this galaxy is correctly found to be outside of our redshift range of interest (shown here as the red curve with ), so we do not use it for our reionization analysis but discuss it in Section 5.2. Note – the remaining 13 objects have photometric redshift distributions outside of the range plotted here.
Object ID R.A. Dec.
[deg] [deg] [erg s cm] [Å]
A2744_2036 0.97
A2744_2346 1.00
A2744_2345 0.99
A2744_2261 0.79
A2744_2503 0.36
A2744_2257 0.54
A2744_20236 0.42
A2744_1040 0.04
A2744_2248 0.96
M0416_99 0.78
M0416_286 0.66
M0416_743 0.07
M0416_1956 0.91
M0416_1997 0.90
M0416_22746 0.62
M1149_23695 0.77
M1149_3343 0.04
M1149_1428 0.25
M1149_945 0.16
M2129_2633 0.20
M2129_2661 0.07
M2129_1556 0.01
RXJ1347_1831 0.06
RXJ1347_656 0.72
RXJ1347_101 0.20
RXJ1347_1368 0.34
RXJ1347_1280 0.03
RXJ2248_1006 0.92
RXJ2248_2086 0.48
Note. IDs for A2744, M0416 and M1149 match the ASTRODEEP catalogue IDs (Merlin et al., 2016; Di Criscienzo et al., 2017). These listed intrinsic magnitudes are calculated using and the listed magnifications and errors. This is the photometric redshift from EAzY integrated between and , i.e. the total probability of the object to have. Flux and EW limits are . All EW are rest-frame. We stress that the EW limits only hold if the Ly is actually in the KMOS range, which has probability given by . This object was spectroscopically confirmed by Laporte et al. (2017a) at . This object was spectroscopically confirmed by Tamura et al. (2018) at .
Table 2: KLASS targets with solutions

To make accurate inferences for reionization it is important to have uniform and well-understood target selection functions for the sources we use. At the time of target selection for KLASS not all deep HFF data were available, nor were sophisticated ICL removal techniques developed (e.g., Merlin et al., 2016; Morishita et al., 2017; Livermore et al., 2017). This led to heterogeneous target selections. However, for this analysis we now use the most up-to-date photometry available to create a sub-sample for analysis with a homogeneous selection function. We demonstrate in Appendix B that this sub-sample is not a biased selection from the final parent catalogues.

Deep, multi-band HST, Spitzer-IRAC and HAWK-I photometry is now available for all our targets through the CLASH, SURFSUP, and HFF programs (Postman et al., 2012; Bradač et al., 2014; Huang et al., 2016a; Lotz et al., 2017). For A2744, M0416 and M1149 we used the ASTRODEEP photometric catalogue which removed foreground intra-cluster light (Castellano et al., 2016; Merlin et al., 2016; Di Criscienzo et al., 2017). For M2129, RXJ1347 and RXJ2248 we created our own catalogues based on the ASTRODEEP methodology (M. Bradač et al., in prep). Of the 56 high-redshift candidate targets we assigned to KMOS IFUs, 46 have matches in these final deep catalogues (including the 3 images of a multiply-imaged system in RXJ2248).

To determine why 10 targets had no match in the final photometric catalogues we examined our target selection catalogues. We used preliminary versions of the ASTRODEEP catalogues for A2744, M0416 and M1149 in our initial selection, so all the objects targeted in A2744 and M0416 have matches in the final catalogues. 3 targets do not appear in the final M1149 catalogue, these objects were never in the preliminary ASTRODEEP catalogue but were selected from alternative preliminary HFF catalogues. 3 targets from M2129, 3 targets from RXJ1347 and 1 target from RXJ2248 have no matches in the final catalogues, which was expected as they were selected from an ensemble of preliminary photometric catalogues with shallower photometry, and narrower wavelength coverage compared to our final catalogues. 3 of the unmatched objects were Category 1 targets. These missing targets were likely faint in the initial photometry and so turn out to be spurious in deep photometry.

Photometric redshift distributions were obtained from the final catalogues with the EAzY code (Brammer et al., 2008). We perform the EAzY fit to the entire photometric dataset, and obtain photometric redshift posteriors without the magnitude prior (which weights bright objects to lower redshifts based on observations of field galaxies and may be inappropriate for our lensed sources). As described in Section 4.1 our inference framework uses the full photometric redshift distribution thus we can robustly use all objects with non-zero probability of being in our redshift range of interest for our inferences.

Taking the 43 high redshift KMOS targets matched in the catalogues (excluding the three images of the galaxy described in Appendix A) we then use the photometric redshift distributions to select objects which could be in the KMOS YJ range (). We calculate using the normalised EAzY photometric redshift distribution for each object to find the total probability of the object being within that redshift range. We select 30 objects with (though the majority have a much higher probability of being in that redshift range). The photometric redshift distributions of these objects within the KMOS YJ range are plotted in Figure 4.

We examined the final deep photometry of the 13 objects which dropped out of the KMOS YJ range in this selection, which include 6 Category 1 targets. As expected, the selection of these objects shifts to lower redshifts now the full photometry is available. The majority of them have detections in the bluest bands which would negate a Lyman Break, and several are clearly passive galaxies when the IRAC bands are included.

Due to the relatively small sample size, we choose to perform our inference at , so we select only objects with some probability to have . We calculate . We select 29 objects with probability of being within this redshift range (21 have probability, and 13 probability). One object has and is excluded from our inference. This is M1149_JD1, recently spectroscopically confirmed at by Hashimoto et al. (2018) in ALMA, who also show a tentative Ly detection from X-shooter. As this galaxy’s photometric redshift distribution clearly puts it at we do not include it in our reionization inferences. Its can be seen in Figure 4 (red line) and we discuss our observations of it, including a tentative detection of flux at the same wavelength as the X-shooter flux, in Section 5.2.

Our inference uses the full distribution, to robustly account for any probability of an object being a lower redshift contaminant. The median and standard deviation of best-fit photometric redshifts over this range for the sub-sample of 29 objects is . These objects and their observed properties, including are listed in Table 2. We demonstrate that this sub-sample is not a biased sample of the final photometric catalogues in Appendix B.

We also cross-checked our Table 2 with publicly available spectroscopic catalogues from ground-based follow-up at optical wavelengths for clusters A2744 (Mahler et al., 2018), M0416 (Balestra et al., 2016; Caminha et al., 2017), M1149 (Grillo et al., 2016), M2129 (Monna et al., 2017) and RXJ2248 (Karman et al., 2015, 2016). We found no matches in those catalogues for any of the objects in Table 2. Non-detections of these objects in optical spectroscopy lends credence to their selection as candidates.

Two objects have been spectroscopically confirmed at by groups. A2744_2248 (a.k.a. A2744_YD4) was confirmed at via [Oiii]88m emission in ALMA, and a tentative Ly emission line was also reported with line flux  erg s cm and (Laporte et al., 2017a), which is well below our limit for that object. As discussed in Section 4.3 we find that treating the object as a detection in our inference has a negligible impact on our inferred limits on the neutral fraction. M0416_99 (a.k.a M0416_Y1) was also confirmed via [Oiii]88m emission in ALMA observations at by Tamura et al. (2018). They also observed the object with X-shooter and found no rest-frame UV emission lines, with a Ly flux limit of  erg s cm (if the line if offset by up to 250 km s). Our KMOS flux median limits are of a comparable depth ( erg s cm).

We obtain magnification estimates for each object using the publicly available HFF lens models444https://archive.stsci.edu/prepds/frontier/lensmodels/. We take the best-fit magnifications from the most recent versions of all available lens models for each object, drop the highest and lowest magnifications to produce an approximate range of estimated magnifications, . We list the median magnification from this sub-sample, and the upper and lower bounds in Table 2. For the inference, we assume magnifications are log-normally distributed with mean given by the median and standard deviation given by half the range of , which is a reasonable fit to the distribution of magnifications from the models. For M2129 and RXJ1347, the only non-HFF clusters, we use the magnification distribution from the Bradač group lens models (Huang et al., 2016b; Hoag et al., 2019) and obtain mean and standard deviation log magnifications. As discussed in Section 4.1 by using the EW in our inference, which is independent of magnification (as opposed to flux), our results are quite robust to magnification uncertainties.

We calculate flux and Ly EW limits for individual objects as in Section 3.4, using Equations 2 and 3. The median intrinsic UV absolute magnitude (i.e., corrected for magnification) of the sample is . The median observed flux upper limit in this sub-sample is  erg s cm, and the median rest-frame Ly EW upper limit is  Å.

4.3 Inference on the IGM neutral fraction

We use 1D spectra and uncertainties as a function of wavelength for the 29 objects described above to infer the IGM neutral fraction at using Equation 7 to calculate the posterior distribution of .

We obtain the flux density spectra using the cubes for each object, extracting flux and noise in a circular aperture with , and apply a rescaling to both to account for the incomplete recovery of flux in the aperture, and a constant rescaling to the noise spectrum to ensure the S/N distribution of pixels in each spectrum is a Normal distribution.

In Figure 5 we plot the posterior distribution for obtained using our observations of the 29 KLASS targets, as well as Keck/MOSFIRE observations of 8 LBGs from the Brightest of Reionizing Galaxies survey (BoRG, Trenti et al., 2011; Bradley et al., 2012; Schmidt et al., 2014b) described by Treu et al. (2013). Using the BoRG sample allows us to cover a broader range in intrinsic magnitudes spanning opposite ends of the galaxy UV luminosity function: the IGM attenuation of Ly from UV bright and UV faint galaxies is expected to be different due to differing Ly escape paths through their interstellar and circumgalactic media (e.g., Stark et al., 2010; Stark et al., 2017; Mason et al., 2018b).

These two sets of independent observations, both indicate a predominantly neutral IGM at . The BoRG data alone produce a lower limit of (68%) and for the KLASS data alone (68%). Lower limits from the combined dataset are (68%) and (95%).

By exploiting gravitational lensing, the KLASS sample sets much lower limits on the Ly EW for intrinsically UV faint galaxies (which produce the strongest constraints on reionization’s mid-stages, M18a, ) than is possible in blank fields. Our new KLASS sample also demonstrates how increasing the number of sources for the inference produces much tighter constraints on the IGM neutral fraction compared to the 8 BoRG sources.

To test whether the inclusion of objects with candidate Ly emission in GLASS data biased our sample, we tested the inference with and without including the Category 1 targets (which were specifically targeted in KLASS because they had candidate Ly emission in the HST data). We found no significant difference in the posteriors. We also tested the inference with and without including the marginal detection of Ly in object A2744_2248 by Laporte et al. (2017a) with spectroscopic confirmation from Oiii emission in ALMA observations. We use the EW reported by Laporte et al. (2017a)  Å, which is well below our limit for that object ( Å). Despite the potential detection, the posterior distribution for this single object strongly favours a mostly neutral IGM due to its very low EW and low significance. We did our inference using both our KMOS spectra and the Laporte et al. (2017a) measurement for this object and found it to have a negligible impact on our final posterior (changing the inferred limit by only ), demonstrating that deep limits on non-detections have a lot of power in our inferences. Our quoted posterior limits include the object as a non-detection.

Figure 5: Posterior probability distribution for the IGM neutral fraction at obtained using Equation 7 and the EW spectra from the KLASS sample described in Section 4.2 and the BoRG sample described by Treu et al. (2013). The blue line and shaded regions show the posterior from the combined datasets, and its 68% and 95% confidence regions (the darkest region is the 68% confidence range).

5 Discussion

We discuss our new lower limit on the neutral fraction and the implications for the timeline of reionization in Section 5.1, and show it favours reionization driven by UV faint galaxies with a low ionizing photon escape fraction. In Section 5.2 we discuss the recent tentative detection of Ly at by (Hashimoto et al., 2018) and show it is not inconsistent with our results. In Section 5.3 we discuss our EW limits on NV and Civ emission. Finally, in Section 5.4 we present a comparison of the KMOS ETC and our achieved S/N for background-limited observations.

5.1 The timeline of reionization

Figure 6: The redshift evolution of the volume average neutral hydrogen fraction of the IGM. Our new lower limit is shown in orange, with the horizontal errorbar at the 68% confidence level. We also plot measurements derived from observations of: the evolving Ly EW distribution at (orange filled star M18a, ) previous estimates from the fraction of LBGs emitting Ly (open black star, Mesinger et al., 2015); the clustering of Ly emitting galaxies (square, Ouchi et al., 2010; Sobacchi & Mesinger, 2015); Ly and Ly forest dark fraction (circle - 68% limits, McGreer et al., 2014); and QSO damping wings (diamond, Davies et al., 2018; Greig et al., 2018). We offset the constraints at by for clarity. We also plot the Planck Collaboration et al. (2016) redshift range of instantaneous reionization (black pentagon). We show median model reionization histories derived from the Mason et al. (2015) UV luminosity function models as coloured lines. We plot models obtained from integrating the luminosity function down to two magnitude limits – (purple dashed line) and (darkest blue solid line) and drawing from uniform distributions for the ionizing photon escape fraction % () and clumping factor , and log-normal distribution for the ionizing efficiency with mean and standard deviation  dex. Comparing reionization histories with ionizing escape fraction drawn from a uniform distribution (light green, ) and (medium teal, ), integrating LFs down to in both cases and using the same distribution for the clumping factor and as above.
Figure 7: The redshift evolution of the ‘Ly fraction’ for UV faint galaxies, the fraction of LBGs observed with Ly EW  Å. We plot literature measurements from Stark et al. (2011); Pentericci et al. (2014); Treu et al. (2013); Tilvi et al. (2014); Schenker et al. (2014); Caruana et al. (2014) and De Barros et al. (2017). We add small offsets in redshift for measurements at the same redshifts to ease the display of the data. We also plot the predicted Ly fraction from M18a calculating using galaxies and the neutral fraction constraint ( confidence intervals) as the orange star. We plot the upper limits recovered in this paper as orange lines, with the solid line showing our 68% confidence limit, and the dotted line extending to the 95% confidence limit. We calculate again using . Our constraint is consistent with literature values at the same redshift.

We plot our new limit on the reionization timeline in Figure 6. We also plot other statistically robust constraints from Ouchi et al. (2010); McGreer et al. (2014); Sobacchi & Mesinger (2015); Mesinger et al. (2015); Davies et al. (2018); Greig et al. (2018); Mason et al. (2018a) and the Planck Collaboration et al. (2016). While no increase in compared to the Mason et al. (2018a) constraint is statistically possible within the 95% confidence interval, our new limit, combined with the other recent statistical measurements at , and other estimates (e.g., Caruana et al., 2014; Zheng et al., 2017), provides increasing evidence for the bulk of hydrogen reionization occurring (Greig et al., 2017; Bañados et al., 2017; Davies et al., 2018; Mason et al., 2018a), late in the Planck Collaboration et al. (2016) confidence range.

Accurate measurements of the reionization timeline can help constrain properties of early galaxies. In Figure 6 we show model reionization histories obtained from integrating the Mason et al. (2015) UV luminosity functions, varying the typical reionization parameters: the minimum UV luminosity of galaxies, and the average ionizing photon escape fraction. We see that late reionization is most consistent with either a high minimum UV luminosity of galaxies () and moderate escape fraction (), or with including ultra-faint galaxies ) with low escape fractions ().

There are many degeneracies between these reionization parameters, and certainly the escape fraction is unlikely to be constant for all galaxies at all times (Trebitsch et al., 2017), but non-detections of high-redshift GRB host galaxies, and observations of lensed high-redshift galaxies, and local dwarfs, indicate galaxies fainter than likely exist at (e.g., Kistler et al., 2009; Tanvir et al., 2012; Trenti et al., 2012; Alavi et al., 2014; Weisz & Boylan-Kolchin, 2017; Livermore et al., 2017; Bouwens et al., 2017; Ishigaki et al., 2018). If ultra-faint galaxies do contribute significantly to reionization our result suggests reionization can be completed with low escape fractions, consistent with low redshift estimates of the average escape fraction (Marchi et al., 2018; Rutkowski et al., 2017; Naidu et al., 2018; Steidel et al., 2018).

For comparison with previous high redshift Ly spectroscopic surveys we plot the so-called ‘Ly fraction’, the fraction of LBGs emitting Ly with EW  Å in Figure 7. We compare our new upper limits on the Ly fraction with literature measurements from Stark et al. (2011); Pentericci et al. (2011); Treu et al. (2013); Tilvi et al. (2014); Schenker et al. (2014); Caruana et al. (2014) and De Barros et al. (2017). We also plot the predicted Ly fraction from M18a. Using the M18a model EW distributions we can calculate the Ly fraction as the probability of EW  Å given our constraint on the neutral fraction.

As noted by M18a and Mason et al. (2018b) the Ly EW distribution is likely a function of at least UV magnitude as well as the neutral fraction (see Oyarzún et al., 2017, for a thorough analysis of Ly EW dependencies on galaxy properties), so it can be difficult to compare Ly fraction from samples with different . Hence, when converting from the neutral fraction measurement in this work and M18a we use the model Ly EW distribution for galaxies to compare more easily with the literature values for which that is the typical median UV magnitude. For our Ly fraction limits are (68%), (95%). Using our sample median magnitude, , the limits are not significantly different: (68%), (95%). Our measurements are consistent with the literature values.

We note that our inference assumes no evolution in the emitted Ly EW distribution at fixed UV magnitude from , i.e. the only evolution in the observed EW distribution is due to reionization. Whilst there may be evolution in the amount of Ly escaping the ISM of galaxies with increasing redshift, it is probably increasing as dust masses and HI covering fractions may decrease at higher redshifts and facilitate Ly escape at fixed galaxy mass (Hayes et al., 2011; Oyarzún et al., 2016). In this case we expect our model to underestimate the observed EW distribution, which would suggest an even higher neutral hydrogen fraction given our non-detections. Our model also assumes no significant evolution in the dust spatial distribution and/or CGM opacity between , which could both reduce the Ly EW before the photons reach the IGM. If these effects do significantly decrease Ly EW between , this could lower our constraint on the neutral fraction. In modelling the emitted Ly EW distribution we assume a Gaussian plus Dirac Delta function parameterisation, which has been shown to describe the Ly EW distribution well (Oyarzún et al., 2017). However, choosing another functional form for the distribution will not significantly change the results (Treu et al., 2012; Schenker et al., 2014).

More accurate models of Ly emerging from the ISM are required to improve our inferences. Whilst it is increasingly difficult to directly observe all of the emitted Ly from galaxies, because of the intervening neutral gas, other emission lines could be used as a diagnostic of emerging Ly. For example, Henry et al. (2018) showed that Mgii emission line profiles and escape fractions closely trace those of Ly in Green Peas, low-redshift analogues of high redshift galaxies (Jaskot & Oey, 2014; Yang et al., 2016). As the IGM optical depth to Mgii is much lower than for Ly, observations of Mgii at (which will be possible with JWST) could be used infer the nature of Ly emission at these redshifts.

Additionally, better knowledge of Ly line profiles at are necessary to provide more informative priors on the observed FWHM for our inferences. In particular, high resolution spectroscopy () is needed to resolve the narrow lines expected for UV faint galaxies (Verhamme et al., 2015), and could provide additional constraints via the evolving Ly profile (Pentericci et al., 2018) and the prevalence of double-peaked Ly in the late stages of reionization (Matthee et al., 2018).

We also assume the fraction of low redshift contaminants in our photometric sample is the same as our reference sample from De Barros et al. (2017). Whilst the selection techniques for the two samples are different (ours is based on photometric redshifts, De Barros et al. (2017) uses a colour selection) our targets have extensive multi-wavelength photometry which help rule out low redshift contaminants (e.g., Vulcani et al., 2017; Livermore et al., 2018). Additionally, we use the full photometric redshift distribution from EAzY in our inference which will weight the most convincing high redshift candidates most strongly in our inference, and robustly account for contamination. With the final GLASS Ly candidate sample it will be possible to use the same selections for both the reference EW distribution and the samples for reionization inferences (K. B. Schmidt et al., in prep).

As our inference weights sources by their photometric redshift distribution, the tightest constraints on will be obtained from samples with robust redshift estimates or, ideally, spectroscopic redshifts obtained from other emission lines, and deep Ly EW limits. We note that the objects which constribute the most to our posterior are the objects with the highest probability () of having photometric redshift at due to their SEDs, and we expect these to have consistent high redshift solutions even if the photometric redshift fitting priors are changed. The prospects for large spectroscopic samples at these redshifts is increasing: ALMA is enabling spectroscopic confirmation of galaxies in the sub-mm (e.g., Bradač et al., 2017; Laporte et al., 2017a; Smit et al., 2018; Hashimoto et al., 2018; Tamura et al., 2018), and other UV emission lines have also been confirmed (Stark et al., 2015; Schmidt et al., 2016; Stark et al., 2017; Mainali et al., 2017; Mainali et al., 2018). Future observations with JWST slitless and slit spectroscopy will be able to build large and deep spectroscopic samples of galaxies, ideal for this type of analysis.

Understanding the differing evolution of Ly emission as a function of galaxy properties and environment will be key to understanding how reionization progresses. Here we have shown that a sample of intrinsically UV faint systems at (more likely to live in low density environments) show no significant Ly emission, and favour a mostly neutral IGM. However, Ly has been observed in a handful of UV bright galaxies at (Zitrin et al., 2015; Oesch et al., 2015; Roberts-Borsani et al., 2016; Stark et al., 2017). Mason et al. (2018b) showed that the observed Ly fraction for UV bright galaxies at could not be reproduced with standard reionization models (using the EoS simulations, Mesinger et al., 2016), even when placing them in overdense regions (which reionize early) and giving them high Ly velocity offsets to facilitate Ly IGM transmission. Mason et al. (2018b) proposed those objects have detectable Ly because they have unusually high emitted Ly EW (they were certainly selected to have high nebular line EW, Roberts-Borsani et al., 2016).

Fluctuations in the UV background during reionization, for example, due to the inhomogeneous distribution of ionizing sources, could also contribute to the differing evolution of Ly emission from UV bright and UV faint galaxies by boosting the IGM opacity (transparency) in underdense (overdense) regions (Davies & Furlanetto, 2016; Becker et al., 2018). One important missing piece in our inference is the halo environment of the LBGs. This work assumes a simple mapping between UV luminosity and halo mass. This works well in an average sense (Mason et al., 2015), but deep imaging with JWST could measure the clustering strength and scatter of galaxies in the reionization epoch (Ren et al., 2018), and be used to inform more realistic IGM simulations.

5.2 M1149_JD1 – Ly emission at ?

One target in our observations (known as M1149_JD1, Zheng et al., 2012; Hoag et al., 2018) was recently spectroscopically confirmed at via [Oiii]m emission with ALMA observations (Hashimoto et al., 2018). Our EAzY photometric redshift distribution for this galaxy put it outside of our inference redshift range (all of the is at , see Figure 4), so it was not used in our reionization inference. However, Hashimoto et al. (2018) also report a tentative 4 detection of Ly emission from this galaxy in X-shooter observations at 12271.5 Å with total line flux  erg s cm