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

TL;DR

This study uses spectroscopic observations of faint, lens-amplified galaxies at z~8 to constrain the timeline of cosmic reionization, finding a high neutral hydrogen fraction indicating rapid reionization during this epoch.

Contribution

It provides new constraints on the IGM neutral hydrogen fraction at z~8 using deep spectroscopic data and Bayesian inference, extending understanding of reionization's progression.

Findings

No Lyα emission detected in the sample.

Lower limit on neutral hydrogen fraction > 0.76 at z~8.

Supports models with rapid reionization and low ionizing photon escape fractions.

Abstract

Detections and non-detections of Lyman alpha (Ly$\alpha$) emission from $z>6$ galaxies ($<1$ 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 z > 7, where Ly$\alpha$ redshifts to near infra-red wavelengths. Here we present a search for z > 7.2 Ly$\alpha$ 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 ~7-10 hours, we detect no Ly$\alpha$ emission with S/N>5 in our sample. We determine our observations to be 80% complete for 5$\sigma$ spatially and spectrally unresolved emission lines with integrated line flux $>5.7\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$. We define a photometrically selected sub-sample of 29 targets at $z=7.9\pm0.6$, with a median 5$\sigma$ Ly$\alpha$ EW limit of 58A. 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 $> 0.76 \; (68\%), \; > 0.46 \; (95\%)$ at z ~ 8, providing further evidence of rapid reionization at z~6-8. We show that this is consistent with reionization history models extending the galaxy luminosity function to $M_\textrm{UV} \lesssim -12$, with low ionizing photon escape fractions, $f_\textrm{esc} \lesssim 15\%$.

Tables

Table 1. Table 1: KLASS cluster targets

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 $>0\stackrel{\prime \prime }{.}8$ as measured by stars observed in each DIT. The total exposure time = DIT $\times$ 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 $z=6.11$ targets (Appendix A).

Table 2. Table 2: KLASS targets with P ​ ( 7.2 ≤ z phot ≤ 8.8 ) 𝑃 7.2 subscript 𝑧 phot 8.8 P(7.2\leq z_{\textrm{phot}}\leq 8.8) solutions

Object ID∗	R.A.	Dec.	$m_{\mathrm{F160W}}$	$\mu$	$M_{\text{uv}}^{†}$	$P{\left(z_{\text{phot}}\right)}^{\diamond }$	$f_{\text{lim}}^{‡}\times {10}^{-18}$	$E{W}_{\text{Ly}\alpha }^{‡,\star }$
	[deg]	[deg]					[erg s-1 cm-2]	[Å]
A2744_2036	$3.596087$	$-30.385836$	$26.95\pm 0.07$	${2.4}_{-0.5}^{+7.4}$	$-19.27\pm 0.89$	0.97
A2744_2346	$3.606460$	$-30.380995$	$26.78\pm 0.06$	${1.6}_{-0.5}^{+0.8}$	$-19.89\pm 0.41$	1.00
A2744_2345	$3.606572$	$-30.380932$	$26.49\pm 0.06$	${1.6}_{-0.5}^{+0.8}$	$-20.19\pm 0.41$	0.99
A2744_2261	$3.603996$	$-30.382309$	$27.29\pm 0.10$	${1.7}_{-0.5}^{+1.1}$	$-19.34\pm 0.47$	0.79
A2744_2503	$3.588979$	$-30.378668$	$27.27\pm 0.12$	${2.2}_{-0.7}^{+0.9}$	$-19.04\pm 0.39$	0.36
A2744_2257	$3.598123$	$-30.382393$	$28.62\pm 0.18$	${1.9}_{-0.4}^{+0.8}$	$-17.87\pm 0.36$	0.54
A2744_20236	$3.572523$	$-30.413267$	$28.61\pm 0.24$	${1.8}_{-0.5}^{+1.0}$	$-17.94\pm 0.48$	0.42
A2744_1040	$3.592505$	$-30.401482$	$27.52\pm 0.15$	${14.2}_{-6.3}^{+11.2}$	$-16.79\pm 0.65$	0.04
A2744_2248∗∗	$3.603863$	$-30.382261$	$26.57\pm 0.07$	${1.7}_{-0.5}^{+1.1}$	$-20.06\pm 0.47$	0.96
M0416_99∗∗∗	$64.039162$	$-24.093182$	$26.28\pm 0.05$	${1.5}_{-0.3}^{+0.5}$	$-20.49\pm 0.30$	0.78
M0416_286	$64.037567$	$-24.088116$	$28.20\pm 0.17$	${1.9}_{-0.5}^{+0.3}$	$-18.29\pm 0.31$	0.66
M0416_743	$64.048058$	$-24.081427$	$26.56\pm 0.06$	${1.7}_{-0.2}^{+0.3}$	$-20.07\pm 0.19$	0.07
M0416_1956	$64.060333$	$-24.064962$	$28.16\pm 0.16$	${1.9}_{-0.6}^{+0.2}$	$-18.33\pm 0.28$	0.91
M0416_1997	$64.049583$	$-24.064596$	$27.56\pm 0.17$	${6.3}_{-1.5}^{+39.3}$	$-17.64\pm 1.23$	0.90
M0416_22746	$64.046509$	$-24.061630$	$27.77\pm 0.23$	${8.1}_{-3.0}^{+4.3}$	$-17.15\pm 0.53$	0.62
M1149_23695	$177.382996$	$22.412041$	$28.11\pm 0.14$	${3.6}_{-2.1}^{+0.7}$	$-17.69\pm 0.58$	0.77
M1149_3343	$177.392715$	$22.384718$	$28.64\pm 0.28$	${1.7}_{-0.5}^{+0.4}$	$-17.96\pm 0.42$	0.04
M1149_1428	$177.412216$	$22.394894$	$28.34\pm 0.17$	${7.5}_{-2.8}^{+0.9}$	$-16.67\pm 0.36$	0.25
M1149_945	$177.412079$	$22.389055$	$27.92\pm 0.13$	${9.2}_{-3.2}^{+14.4}$	$-16.87\pm 0.76$	0.16
M2129_2633	$322.345232$	$-7.671373$	$25.65\pm 0.12$	${1.6}_{-0.1}^{+0.1}$	$-21.06\pm 0.13$	0.20
M2129_2661	$322.350848$	$-7.675239$	$26.38\pm 0.17$	${1.7}_{-0.0}^{+0.0}$	$-20.25\pm 0.17$	0.07
M2129_1556	$322.344535$	$-7.688473$	$27.53\pm 0.26$	${4.2}_{-0.2}^{+0.2}$	$-18.11\pm 0.27$	0.01
RXJ1347_1831	$206.896270$	$-11.742338$	$26.30\pm 0.26$	${9.2}_{-0.4}^{+0.4}$	$-18.49\pm 0.26$	0.06
RXJ1347_656	$206.891246$	$-11.752607$	$26.43\pm 0.24$	${20.4}_{-1.2}^{+1.6}$	$-17.49\pm 0.25$	0.72
RXJ1347_101	$206.880973$	$-11.769816$	$25.16\pm 0.15$	${43.9}_{-5.4}^{+10.2}$	$-17.92\pm 0.26$	0.20
RXJ1347_1368	$206.893076$	$-11.760230$	$27.92\pm 0.43$	${16.6}_{-1.1}^{+1.1}$	$-16.22\pm 0.43$	0.34
RXJ1347_1280	$206.896921$	$-11.763833$	$27.28\pm 0.28$	${4.8}_{-0.5}^{+0.7}$	$-18.22\pm 0.31$	0.03
RXJ2248_1006	$342.208379$	$-44.537520$	$25.83\pm 0.17$	${1.6}_{-0.4}^{+0.4}$	$-20.88\pm 0.32$	0.92
RXJ2248_2086	$342.179829$	$-44.525664$	$26.88\pm 0.13$	${41.0}_{-25.5}^{+72.3}$	$-16.28\pm 1.09$	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 $z=8$ and the listed magnifications and errors. ⋄ This is the photometric redshift from EAzY integrated between $z=7.2$ and $z=8.8$, i.e. the total probability of the object to have. ‡ Flux and EW limits are $5\sigma$. ⋆ All EW are rest-frame. We stress that the EW limits only hold if the Ly$\alpha$ is actually in the KMOS range, which has probability given by $P\left(7.2\le z_{\text{phot}}\le 8.8\right)$. ∗∗ This object was spectroscopically confirmed by Laporte et al. (2017a) at $z=8.38$. ∗∗∗ This object was spectroscopically confirmed by Tamura et al. (2018) at $z=8.31$.

Equations

$$P=1-\frac{N_{\text{spurious}}}{N_{\text{pos}}}$$

$$f_{\text{lim},i}=5\times \frac{1}{1-e^{-\frac{r^2}{2{\sigma }_{\text{textsc}psf}^2}}}\sqrt{\frac{2FWH{M}_{\mathrm{res}}}{\Delta \lambda }}{\sigma }_i\times \Delta \lambda$$

$$f_{\text{cont}}(m,z)=f_{0}10^{-0.4m}\frac{c}{{\lambda }_{\alpha }^2(1+z{)}^2}{\left(\frac{{\lambda }_{\text{textsc}uv}}{{\lambda }_{\alpha }}\right)}^{-\beta -2}$$

$$\begin{array}{rl}{f}_{\mathrm{mod}}& (\lambda ,W,m,z_d,\mathrm{FWHM})=\\ & \frac{W{f}_{cont}(m,z_d)(1+z_d)}{\sqrt{2\pi }{\sigma }_{\lambda }}{e}^{-\frac{1}{2}{\left(\frac{\lambda -{\lambda }_d}{{\sigma }_{\lambda }}\right)}^2}\end{array}$$

$$\begin{split}p&(\{f\}\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d},\mathrm{FWHM})=\\ &\prod_{i}^{N}\int_{0}^{\infty}dW\,\biggl{[}\frac{1}{\sqrt{2\pi}\sigma_{i}}e^{-\frac{1}{2}\left(\frac{f_{i}-f_{\mathrm{mod}}(\lambda_{i},W,m,z_{d},\mathrm{FWHM})}{\sigma_{i}}\right)^{2}}\\ &\quad\quad\quad\quad\quad\quad\times p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d})\biggr{]}\end{split}$$

$$\begin{array}{rl}p({\overline{x}}_{\text{hi}},z_d,\mathrm{FWHM}|\{f\},m,\mu )\propto & p(\{f\}|{\overline{x}}_{\text{hi}},m,\mu ,z_d,\mathrm{FWHM})\\ & \times p({\overline{x}}_{\text{hi}})p(z_d)p(\mathrm{FWHM})\end{array}$$

$$\begin{array}{rl}p({\overline{x}}_{\text{hi}}|\{f\},m,\mu )& \propto {\int }_{z_{\text{min}}}^{z_{\text{max}}}d{z}_{d}p(\{f\}|{\overline{x}}_{\text{hi}},m,\mu ,z_d)p(z_d)\\ & +\prod _i^N\frac{1}{\sqrt{2\pi }{\sigma }_i}{e}^{-\frac{f_i^2}{2{\sigma }_i^2}}\left(1-{\int }_{z_{\text{min}}}^{z_{\text{max}}}d{z}_{d}p(z_d)\right)\end{array}$$

$$\frac{S}{N}=\frac{\sqrt{NDIT}\times S_{\text{source}}}{\sqrt{S_{\text{source}}+S_{\text{bkg}}+n_{\text{spat}}(\text{DC}\times \text{DIT}+{\text{RON}}^2)}}$$

$$\begin{array}{rl}p& (f_i|{\overline{x}}_{\text{hi}},m,\mu ,z_d,\mathrm{FWHM})=\\ & {\int }_0^{\infty }dWp(f_i|W,m,z_d,\mathrm{FWHM})p(W|{\overline{x}}_{\text{hi}},m,\mu ,z_d)\end{array}$$

$$\begin{array}{rl}p& (f_i|W,m,z_d,\mathrm{FWHM})=\\ & \frac{1}{\sqrt{2\pi }{\sigma }_i}\exp \left[-\frac{{\left(f_i-f_{\mathrm{mod}}({\lambda }_i,W,m,z_d,\mathrm{FWHM})\right)}^2}{2{\sigma }_i^2}\right]\end{array}$$

$$\begin{array}{rl}{f}_{\mathrm{mod}}& ({\lambda }_i,W,m,z_d,\mathrm{FWHM})=\\ & \frac{W{f}_{cont}(m,z_d)(1+z_d)}{\sqrt{2\pi }{\sigma }_{\lambda }}\exp \left[-\frac{{\left({\lambda }_i-{\lambda }_d\right)}^2}{2{\sigma }_{\lambda }^2}\right]\end{array}$$

$$p(W|{\overline{x}}_{\text{hi}},m,\mu ,z_d)={\int }_{-\infty }^{\infty }d{M}_{\text{textsc}uv}p(W|{\overline{x}}_{\text{hi}},M_{\text{textsc}uv})p(M_{\text{textsc}uv}|m,\mu ,z_d)$$

$$p(M_{\text{textsc}uv}|m,\mu ,z_d)=\frac{1}{\sqrt{2\pi {\sigma }_M^2}}\exp \left[-\frac{{\left(M_{\text{textsc}uv}-M_{\text{textsc}uv,\text{mod}}(\mu ,z_d)\right)}^2}{2{\sigma }_M^2}\right]$$

$$\begin{array}{rl}& p(\{f\}|{\overline{x}}_{\text{hi}},m,\mu ,z_d,\mathrm{FWHM})=\\ & \prod _i^N{\int }_0^{\infty }dW\frac{1}{\sqrt{2\pi }{\sigma }_i}{e}^{-\frac{1}{2}{\left(\frac{f_i-f_{\mathrm{mod},\mathrm{i}}}{{\sigma }_i}\right)}^2}p(W|{\overline{x}}_{\text{hi}},m,\mu ,z_d)\end{array}$$

$$\begin{array}{rl}p({\overline{x}}_{\text{hi}},z_d,\mathrm{FWHM}|\{f\},m,\mu )\propto & p(\{f\}|{\overline{x}}_{\text{hi}},m,\mu ,z_d,\mathrm{FWHM})\\ & \times p(z_d)p({\overline{x}}_{\text{hi}})p(\mathrm{FWHM})\end{array}$$

$$p({\overline{x}}_{\text{hi}}|\{f\},m,\mu )=\int d\mathrm{FWHM}\int d{z}_{d}p({\overline{x}}_{\text{hi}},z_d,\mathrm{FWHM}|\{f\},m,\mu )$$

$$\begin{array}{rl}p({\overline{x}}_{\text{hi}}|\{f\},m,\mu )\propto & {\int }_{z_{\text{min}}}^{z_{\text{max}}}d{z}_{d}p(\{f\}|{\overline{x}}_{\text{hi}},m,\mu ,z_d)p(z_d)\\ & +\prod _{i}p(\{f\}=0)\left(1-{\int }_{z_{\text{min}}}^{z_{\text{max}}}d{z}_{d}p(z_d)\right)\end{array}$$

Full text

††pubyear: 2019††pagerange: Inferences on the Timeline of Reionization at $z\sim 8$ From the KMOS Lens-Amplified Spectroscopic Survey–C.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 $z>6$ (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 $z>6$ 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$\alpha$ (rest-frame 1216 Å) emission in galaxies, using near infra-red (NIR) spectroscopy.

Ly$\alpha$ 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$\alpha$ emission in abundance up to $z\sim 6$ (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$\alpha$ emission, and the scale length of the Ly$\alpha$ 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$\alpha$ 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 $z\sim 7$ (Pentericci et al., 2011, 2014, 2018; Hoag et al., 2019) but exploring the earliest stages of reionization requires us to observe Ly$\alpha$ at even higher redshifts. Only a handful of Ly$\alpha$ emitters have been confirmed at z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7.5 (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$\alpha$ 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$\alpha$ flux. Hence, the observed decline in Ly$\alpha$ emission at $z>6$ could be partially due to redshift-dependent slit-losses as well as reionization.

In this paper we present a search for z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7.2 Ly$\alpha$ emission in NIR spectroscopy of 53 intrinsically faint LBG candidates (M_{\textsc{uv}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,-20), 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 $z=6.11$ 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 $z>3.8$ 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$\alpha$ emission than traditional slit spectrographs.

We use our new deep spectroscopic observations to infer the average IGM neutral hydrogen fraction (${\overline{x}_{\textrm{{{hi}}}}}$) at $z\sim 8.$ Mason et al. (2018a, hereafter M18a) presented a flexible Bayesian framework to directly infer ${\overline{x}_{\textrm{{{hi}}}}}$ from detections and non-detections of Ly$\alpha$ emission from LBGs. The framework combines realistic inhomogeneous reionization simulations and models of galaxy properties. That work measured ${\overline{x}_{\textrm{{{hi}}}}}=0.59_{-0.15}^{+0.11}$ ($16-84\%$ confidence intervals) at $z\sim 7$. 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$\alpha$ 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 ${\overline{x}_{\textrm{{{hi}}}}}$ at $z\sim 8$. 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 $(\Omega_{\Lambda},\Omega_{\textrm{m}},\Omega_{\textrm{b}},n,\sigma_{8},H_{0})=$ (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s*-1* Mpc*-1*). 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 $96-99$ (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 $2\aas@@fstack{\prime\prime}8\times 2\aas@@fstack{\prime\prime}8$ field of view, with pixel size $0\aas@@fstack{\prime\prime}2\times 0\aas@@fstack{\prime\prime}2$, and 2048 pixels along the wavelength axis111We use the following definitions for describing 3D spectra in this paper. Pixel: 2D spatial pixel (size $0\aas@@fstack{\prime\prime}2\times 0\aas@@fstack{\prime\prime}2$). Spaxel: the 1D spectrum in a single spatial pixel (spanning the spectral range $\sim 1-1.35\,\mu$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:

To probe the internal kinematics of galaxies at $z\sim 1-3$, with superior spatial resolution compared to surveys in blank fields (Mason et al., 2017). 2. 2.

To investigate z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7 Ly$\alpha$ 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 $z>1$, 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 $z>7$ candidate targets with complete exposures. We use the YJ observing band, giving us access to Ly$\alpha$ emission at $z\sim 7.2-10.1$.

The choice of an IFU instrument for high-redshift Ly$\alpha$ observations was motivated by indications that ground-based slit-spectroscopy measures lower Ly$\alpha$ 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$\alpha$ emission at $z>6$. Ly$\alpha$ 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$\alpha$ flux.

By using IFUs our observations should be more complete to spatially extended and/or offset Ly$\alpha$ than traditional slit spectrographs. Mason et al. (2017) showed that only $\sim 60\%$ of emission line flux was contained in $\sim 0\aas@@fstack{\prime\prime}7$ simulated slits (a typical slit-width used for Ly$\alpha$ observations, e.g., Hoag et al., 2017) on KMOS spectra, whereas the full flux is captured within the $2\aas@@fstack{\prime\prime}8\times 2\aas@@fstack{\prime\prime}8$ KMOS field of view. Thus we expect most Ly$\alpha$ flux to be captured within the KMOS IFUs. The $2\aas@@fstack{\prime\prime}8$ wide IFUs cover $\sim 14$ proper kpc at $z\sim 8$, while the UV effective radii of galaxies at these redshifts is only \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,1 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$\alpha$ emission.

2.2 Target selection

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 $0.8-1.6\,\mu$m with spectral resolution $R\sim 150$. 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 $z\sim 1$ 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:

Category 1: 14 objects with marginal (S/N $\sim 3$) candidate Ly$\alpha$ emission in the HST grisms, identified by visual inspection of the GLASS data, which fall within the KMOS YJ spectral coverage ($\sim 1-1.35\,\mu$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. 2.

Category 2: 39 LBG candidates selected with $z_{\textrm{phot}}>7.2$, 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 $z=8.38$ (Laporte et al., 2017a, discussed in Section 4.2 and 4.3), M0416_Y1 (M0416_99 in this paper), at $z=8.12$ (Tamura et al., 2018, discussed in Section 4.2), and M1149_JD1, at $z=9.11$ (Hashimoto et al., 2018, discussed in Section 5.2).

An additional three targets were multiple images of the $z=6.11$ system in RXJ2248 (Boone et al., 2013; Balestra et al., 2013; Monna et al., 2014) where we targeted Civ$\lambda$1548,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 $z\sim 1-2$ 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 $z>7$ 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 $1\sigma$ flux limit is $5\times 10^{-18}$ erg s*-1* cm*-2* (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 $3\sigma$ 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 ($R\sim 100$) of the HST grisms means that the line position is uncertain by $\pm 25$ Å 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 $\sim 80\%$ complete for lines with flux $>5.7\times 10^{-18}$ erg s*-1* cm*-2*, which suggests we should have confirmed the majority of the 14 GLASS candidate Ly$\alpha$ 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 $z_{\mathrm{phot}}>7.2$ 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 ($\sim 1-1.35\,\mu$m). The spectral resolution $R\sim 3400$ is sufficient to distinguish Ly$\alpha$ from potential low redshift contaminants with the [Oii]$\lambda 3726,3729$ emission doublet at $z\sim 2$.

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 $\sim 0\aas@@fstack{\prime\prime}6$.

In each cluster, the 3 top priority targets were observed for $1.5\times$ 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 z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,2 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 $\leq 0\aas@@fstack{\prime\prime}8$ (as measured by the star observed in our science frames). The median seeing was $\sim 0\aas@@fstack{\prime\prime}6$. 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 $z\sim 1$ 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 $>1$, 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 $\sim 1.2\times$ 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 \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7\sigma confidence regions, due to bad sky subtraction residuals.

3 Emission Line Search, Purity and Completeness

In this section we describe our search for Ly$\alpha$ 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$\alpha$ 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:

Create a circular aperture with $r=2\sigma_{\textsc{psf}}\sim 0\aas@@fstack{\prime\prime}5\sim 2.5$ pixels (using our median seeing $FWHM_{\mathrm{psf}}=0\aas@@fstack{\prime\prime}6$), which will capture $86\%$ of the total flux for spatially unresolved emission line at the centre of the aperture. 2. 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. 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. 4.

Scan through in wavelength and flag a detection if 3 adjacent wavelength pixels have $S/N>3$. This corresponds to a S/N\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,5 detection of the integrated line flux. 5. 5.

Iterate over 25 apertures centred within 3 pixels ($0\aas@@fstack{\prime\prime}6$) of the IFU centre, i.e. $x=[-3,-1.5,0,1.5,3]$, $y=[-3,-1.5,0,1.5,3]$ where $(x,y)=(0,0)$ is the IFU centre.

Our search covers $\sim 25\times 2000=50,000$ potential emission line positions in each cube. As our detection threshold is $5\sigma$ we would expect a false positive rate of $6\times 10^{-7}$, i.e. $\sim 0.03$ 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 $z=6.11$ 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$\lambda 1551$ emission at 11023.7 Å in the brightest image of the multiply-imaged $z=6.11$ 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 ($R\sim 3400$) 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 $z>7$ LBG candidate in the cube. We determine the detected emission to be associated with a nearby ($\sim 1.1\arcsec$) galaxy with $z_{\textrm{phot}}=4.2$, 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 $z\sim 4$ 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 $z\sim 4$ 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 $z\sim 4$ continuum object described above), allowing us to define the purity of our spectral sample:

[TABLE]

where $N_{\textrm{spurious}}=17$ is the total number of spurious flags (8 unique false detections which were sometimes flagged in multiple apertures) and $N_{\textrm{pos}}=101763\times 25$ 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 $P=1-7\times 10^{-6}$. Our spurious detection rate is $\sim 10\times$ higher than that expected for $5\sigma$ 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\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,-5) 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 $>2.5$ per wavelength pixel, corresponding to S/N\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,4 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$\alpha$ 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$\alpha$ flux and EW limits for our most likely $z_{\mathrm{phot}}$ LBGs candidates. We discuss our limits on other UV lines in Section 5.3.

3.3 Completeness

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${}_{\textrm{spec}}$, spatial position, spatial extent FWHM${}_{\textrm{spat}}$, 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 $\sigma_{x,y,\lambda}$ for each cluster. The $\sigma_{x,y,\lambda}$ cubes are constructed by taking the RMS at every voxel of all the final sky-subtracted cubes which do not contain bright z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,2 sources ($\sim 10$ 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, $\sigma_{x,y,\lambda}$.

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 = $1\times 10^{-17}$ erg s*-1* cm*-2*, observed line FWHM${}_{\textrm{spec}}=4$ Å (the spectral resolution, i.e. unresolved lines), line centred at the IFU centre, with source spatial extent FWHM${}_{\textrm{spat}}=0\arcsec$ (i.e. unresolved point source, the emission will have observed spatial extent with $\textrm{FWHM}_{\textrm{spat,tot}}=\sqrt{\textrm{FWHM}_{\textrm{PSF}}^{2}+\textrm{FWHM}_{\textrm{spat}}^{2}}$). 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$\alpha$ 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$\alpha$ emission centre slits on the UV continuum and use slit-widths $\sim 0\aas@@fstack{\prime\prime}7$, thus in KLASS we are more complete to Ly$\alpha$ 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 $80\%$ completeness over the full wavelength range for lines \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,5.7\times 10^{-18} erg s*-1* cm*-2*, centred within $<0\aas@@fstack{\prime\prime}8$ of the IFU centre and with intrinsic line FWHM{}_{\textrm{spec}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,250 km s*-1*, assuming $z=8$ to calculate FWHM${}_{\textrm{spec}}$ (median over all clusters). For wavelength ranges where the noise level is below the median across the whole spectrum, we reach $80\%$ completeness for $5\sigma$ lines \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,3.2\times 10^{-18} erg s*-1* cm*-2*, centred with $<0\aas@@fstack{\prime\prime}9$ of the IFU centre and with intrinsic line FWHM{}_{\textrm{spec}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,550 km s*-1*. The completeness is fairly flat for Ly$\alpha$ spatial extent \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,0\aas@@fstack{\prime\prime}6 (total extent \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,0\aas@@fstack{\prime\prime}8) demonstrating our good completeness for spatially extended Ly$\alpha$ emission, with the normalisation of the completeness as a function of FWHM${}_{\textrm{spat}}$ scaling with the completeness at a given total line flux.

3.4 Flux and equivalent width limits

To calculate average flux limits for each cluster we take the average 3D noise spectrum for each cluster, $\sigma_{x,y,\lambda}$ (created by taking the RMS at every voxel across the $\sim 10$ IFUs observing high redshift candidates in each cluster). We then create a 1D noise spectrum, $\sigma_{\lambda}$, by summing the average noise at each wavelength pixel in a circular aperture with radius $r=2\sigma_{\textsc{psf}}$ (where we use our median seeing $\mathrm{FWHM}_{\textsc{psf}}=0\aas@@fstack{\prime\prime}6$). At each wavelength pixel, $i$ the flux limit in erg s*-1* cm*-2* is given by:

[TABLE]

Here we obtaining an estimate of the integrated noise for an emission line with observed $FWHM_{\mathrm{res}}=4$ Å or $\approx 110$ km s*-1* (the instrumental resolution), and use a threshold integrated $S/N=5$. The term in the denominator accounts for the fact that the apertures only capture a fraction of the flux. For $r=2\sigma_{\textsc{psf}}$ this results in a rescaling of 1.16. The spectral pixel width of KMOS YJ is $\Delta\lambda=1.75$ Å. The above calculation assumes the emission is spatially and spectrally unresolved by KMOS, which is reasonable given the expectation that Ly$\alpha$ emission from UV faint galaxies is likely to be more spatially compact and have lower linewidth than Ly$\alpha$ 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 $f_{\textrm{lim}}\propto\sqrt{FWHM/4\,\textrm{\AA}}$.

The $5\sigma$ flux limits for all clusters are shown in Figure 2. The median flux limit is $4.5\times 10^{-18}$ erg s*-1* cm*-2*, and the range of medians for each cluster is $3.9-6.4\times 10^{-18}$ erg s*-1* cm*-2*.

Rest-frame Ly$\alpha$ equivalent widths are $W=(1+z)^{-1}f(\lambda)/f_{\textrm{cont}}$, where $z=\lambda/\lambda_{\alpha}-1$ (with $\lambda_{\alpha}=1216$ Å), and we define the continuum flux:

[TABLE]

where $f_{0}=3.631\times 10^{-20}$ erg s*-1* Hz*-1* cm*-2*, $m$ is the apparent magnitude of the UV continuum, $c$ is the speed of light, $\lambda_{\textsc{uv}}$ is the rest-frame wavelength of the UV continuum (usually 1500 Å), and $\beta$ is the UV slope. We assume $\beta=-2$, consistent with $z\sim 7$ 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 $z\sim 8$ we use the method described by M18a. This forward-models the observed rest-frame Ly$\alpha$ EW distribution as a function of the neutral fraction and galaxy UV magnitude, $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$, 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 $z\sim 6$ Ly$\alpha$ 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$\alpha$ velocity offsets. Here, as in M18a, we use the recent comprehensive $z\sim 6$ Ly$\alpha$ 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$\alpha$ optical depths along millions of sightlines in simulated IGM cubes for a grid of volume-averaged ${\overline{x}_{\textrm{{{hi}}}}}$ values. As the size of ionised regions during reionization is expected to be nearly independent of redshift at fixed ${\overline{x}_{\textrm{{{hi}}}}}$ (as there is little difference in the matter power spectrum from $z\sim 7-11$, McQuinn et al., 2007), we use the same $z\sim 7$ cubes as used by M18a rather than generating new $z\sim 8$ 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$\alpha$ 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: $p({\overline{x}_{\textrm{{{hi}}}}}\,|\,\{f\},m,\mu)$ for each galaxy, where $\{f\}$ is an observed flux density spectrum as a function of wavelength, $m$ is the observed apparent UV magnitude, and $\mu$ 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 $z_{d}$ with observed rest-frame EW, $=W$, 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 $z>6$ Ly$\alpha$ emission with linewidth $\mathrm{FWHM}_{\mathrm{spec}}\sim 200-450$ km s*-1*, 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 $z\sim 7$ 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 $\leavevmode\nobreak\ \sim 100-400$ km s*-1*.

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:

[TABLE]

where $z_{d}=\lambda_{d}/\lambda_{\alpha}-1$, with $\lambda_{\alpha}=1216$ Å, is the redshift of an emission line, $W$ is the rest-frame equivalent width of the emission line, $f_{\mathrm{cont}}$ is the flux density of the continuum calculated using Equation 3 using the observed continuum apparent magnitude $m$, and $\sigma_{\lambda}=FWHM/2.355$ is the spectral linewidth.

The likelihood of observing a 1D flux density spectrum $\{f\}=f(\lambda_{i})$ for an individual galaxy (where $i$ is the wavelength pixel index), given our model where the true EW is drawn from the conditional probability distribution $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d})$ is:

[TABLE]

where $\sigma_{i}$ is the uncertainty in flux density at wavelength pixel $i$ and there are a total of $N$ wavelength pixels in the spectrum. $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d})$ 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 $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$ model outputs from M18a with the probability distribution for each galaxy’s absolute UV magnitude, including errors on $m$ and $\mu$ (Equation 13).

We note that the range of neutral fraction in the EoS simulations is ${\overline{x}_{\textrm{{{hi}}}}}=0.01-0.95$. In order to correctly calculate posteriors and confidence intervals we set the likelihood at ${\overline{x}_{\textrm{{{hi}}}}}$ such that we expect to observe no Ly$\alpha$ flux at all (in a fully neutral universe). I.e. $p(\{f\}\,|\,{\overline{x}_{\textrm{{{hi}}}}}=1,m,\mu,z_{d},\mathrm{FWHM})=\prod_{i}^{N}\frac{1}{\sqrt{2\pi}\sigma_{i}}\exp{(-f_{i}^{2}/2\sigma_{i}^{2})}$.

Given our relatively small sample size, we choose to restrict our inference to $z\sim 8$, thus for ease of computation we evaluate $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d})$ at $z_{d}=8$; this has a negligible impact on the final likelihood. We keep $z_{d}$ 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 $M_{\textsc{uv}}$ of the assumed intrinsic EW distribution.

Using Bayes’ theorem, the posterior distribution for ${\overline{x}_{\textrm{{{hi}}}}}$, $z_{d}$ and FWHM is:

[TABLE]

We use a uniform prior on ${\overline{x}_{\textrm{{{hi}}}}}$ between 0 and 1, $p({\overline{x}_{\textrm{{{hi}}}}})$, and use the photometric redshift distribution for the prior $p(z_{d})$. As we are only interested in the posterior probability of ${\overline{x}_{\textrm{{{hi}}}}}$ we can marginalise over FWHM and $z_{d}$ for each galaxy. We use a log-normal prior on FWHM with mean depending on $M_{\textsc{uv}}$ 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$\alpha$ FWHM at $z>7$ ($\sim 100-400$ km s*-1*). To account for the incomplete wavelength coverage, we use the fact that if the object has Ly$\alpha$ outside of the KMOS wavelength range (covering $[z_{\textrm{min}}=7.2$, $z_{\textrm{max}}=10.1]$) we would measure a non-detection in our data. Thus the posterior for ${\overline{x}_{\textrm{{{hi}}}}}$ from one galaxy is:

[TABLE]

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 $z\sim 1-2$ 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 $z\sim 8$ only, with larger samples, it will be possible to measure ${\overline{x}_{\textrm{{{hi}}}}}(z)$ 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

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 $z=6.11$ 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 $z=6.11$ galaxy described in Appendix A) we then use the photometric redshift distributions to select objects which could be in the KMOS YJ range ($7.2<z_{\textrm{phot}}<10.1$). We calculate $P(7.2<z_{\textrm{phot}}<10.1)=\int_{7.2}^{10.1}p(z_{\textrm{phot}})dz_{\textrm{phot}}$ 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 $P(7.2<z_{\textrm{phot}}<10.1)>0.01$ (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 $z>7$ Lyman Break, and several are clearly $z\sim 1$ passive galaxies when the IRAC bands are included.

Due to the relatively small sample size, we choose to perform our inference at $z\sim 8$, so we select only objects with some probability to have $7.2<z_{\textrm{phot}}<8.8$. We calculate $P(7.2<z_{\textrm{phot}}<8.8)=\int_{7.2}^{8.8}p(z_{\textrm{phot}})dz_{\textrm{phot}}$. We select 29 objects with $>1\%$ probability of being within this redshift range (21 have $>10\%$ probability, and 13 $>60\%$ probability). One object has $z_{\textrm{phot}}>9$ and is excluded from our inference. This is M1149_JD1, recently spectroscopically confirmed at $z=9.11$ by Hashimoto et al. (2018) in ALMA, who also show a tentative Ly$\alpha$ detection from X-shooter. As this galaxy’s photometric redshift distribution clearly puts it at $z>9$ we do not include it in our $z\sim 8$ reionization inferences. Its $p(z_{\textrm{phot}})$ can be seen in Figure 4 (red line) and we discuss our observations of it in Section 5.2.

Our inference uses the full $p(z_{\textrm{phot}})$ 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 $z_{\textrm{phot}}=7.9\pm 0.6$. These objects and their observed properties, including $P(7.2<z_{\textrm{phot}}<8.8)$ 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 $z>7$ candidates.

Two objects have been spectroscopically confirmed at $7.2<z<8.8$ by groups. A2744_2248 (a.k.a. A2744_YD4) was confirmed at $z=8.38$ via [Oiii]88$\mu$m emission in ALMA, and a tentative Ly$\alpha$ emission line was also reported with line flux $(1.82\pm 0.64)\times 10^{-18}$ erg s*-1* cm*-2* and $EW=10.7\pm 2.7$ (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]88$\mu$m emission in ALMA observations at $z=8.31$ by Tamura et al. (2018). They also observed the object with X-shooter and found no rest-frame UV emission lines, with a $5\sigma$ Ly$\alpha$ flux limit of $<8.0\times 10^{-}18$ erg s*-1* cm*-2* (if the line if offset by up to 250 km s*-1*). Our KMOS flux median limits are of a comparable depth ($<3.4\times 10^{-}18$ erg s*-1* cm*-2*).

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 $1\sigma$ range of estimated magnifications, $\{\mu\}$. 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 $\log_{10}\{\mu\}$ and standard deviation given by half the range of $\log_{10}\{\mu\}$, 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$\alpha$ 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 $M_{\textsc{uv}}=-18.2$. The median observed flux $5\sigma$ upper limit in this sub-sample is $<3.6\times 10^{-18}$ erg s*-1* cm*-2*, and the median rest-frame Ly$\alpha$ EW $5\sigma$ upper limit is $<58$ Å.

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 $z\sim 8$ using Equation 7 to calculate the posterior distribution of ${\overline{x}_{\textrm{{{hi}}}}}$.

We obtain the flux density spectra using the cubes for each object, extracting flux and noise in a circular aperture with $r=2\sigma_{psf}$, 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 ${\overline{x}_{\textrm{{{hi}}}}}$ obtained using our observations of the 29 $z\sim 8$ KLASS targets, as well as Keck/MOSFIRE observations of 8 $z\sim 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$\alpha$ from UV bright and UV faint galaxies is expected to be different due to differing Ly$\alpha$ 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 $z\sim 8$. The BoRG data alone produce a lower limit of ${\overline{x}_{\textrm{{{hi}}}}}>0.34$ (68%) and for the KLASS data alone ${\overline{x}_{\textrm{{{hi}}}}}>0.76$ (68%). Lower limits from the combined dataset are ${\overline{x}_{\textrm{{{hi}}}}}>0.76$ (68%) and ${\overline{x}_{\textrm{{{hi}}}}}>0.46$ (95%).

By exploiting gravitational lensing, the KLASS sample sets much lower limits on the Ly$\alpha$ 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$\alpha$ 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$\alpha$ emission in the HST data). We found no significant difference in the posteriors. We also tested the inference with and without including the $z=8.38$ marginal detection of Ly$\alpha$ 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) $W=10.7\pm 2.7$ Å, which is well below our $5\sigma$ limit for that object ($<53$ Å). 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 $\Delta{\overline{x}_{\textrm{{{hi}}}}}\sim 0.01$), 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.

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$\alpha$ at $z=9.11$ 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

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 ${\overline{x}_{\textrm{{{hi}}}}}$ compared to the Mason et al. (2018a) constraint is statistically possible within the 95% confidence interval, our new limit, combined with the other recent ${\overline{x}_{\textrm{{{hi}}}}}$ statistical measurements at $z\sim 7$, and other estimates (e.g., Caruana et al., 2014; Zheng et al., 2017), provides increasing evidence for the bulk of hydrogen reionization occurring $z\sim 6-8$ (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 ($M_{\textsc{uv}}<-17$) and moderate escape fraction ($\langle f_{\textrm{esc}}\rangle=20\%$), or with including ultra-faint galaxies $M_{\textsc{uv}}<-12$) with low escape fractions (\langle f_{\textrm{esc}}\rangle\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,15\%).

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 $M_{\textsc{uv}}=-17$ likely exist at $z\sim 8$ (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$\alpha$ spectroscopic surveys we plot the so-called ‘Ly$\alpha$ fraction’, the fraction of LBGs emitting Ly$\alpha$ with EW $\geq 25$ Å in Figure 7. We compare our new upper limits on the Ly$\alpha$ 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$\alpha$ fraction from M18a. Using the M18a model EW distributions $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$ we can calculate the Ly$\alpha$ fraction as the probability of EW $\geq 25$ Å given our constraint on the neutral fraction.

As noted by M18a and Mason et al. (2018b) the Ly$\alpha$ 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$\alpha$ EW dependencies on galaxy properties), so it can be difficult to compare Ly$\alpha$ fraction from samples with different $M_{\textsc{uv}}$. Hence, when converting from the neutral fraction measurement in this work and M18a we use the model Ly$\alpha$ EW distribution for $M_{\textsc{uv}}=-20$ galaxies to compare more easily with the literature values for which that is the typical median UV magnitude. For $M_{\textsc{uv}}=-20$ our Ly$\alpha$ fraction limits are $f_{\textrm{Ly}\alpha}<0.11$ (68%), $<0.27$ (95%). Using our sample median magnitude, $M_{\textsc{uv}}=-18.2$, the limits are not significantly different: $f_{\textrm{Ly}\alpha}<0.08$ (68%), $<0.24$ (95%). Our measurements are consistent with the literature values.

We note that our inference assumes no evolution in the emitted Ly$\alpha$ EW distribution at fixed UV magnitude from $z\sim 6-8$, i.e. the only evolution in the observed EW distribution is due to reionization. Whilst there may be evolution in the amount of Ly$\alpha$ 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$\alpha$ 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 $z\sim 6-8$, which could both reduce the Ly$\alpha$ EW before the photons reach the IGM. If these effects do significantly decrease Ly$\alpha$ EW between $z\sim 6-8$, this could lower our constraint on the neutral fraction. In modelling the emitted Ly$\alpha$ EW distribution we assume a Gaussian plus Dirac Delta function parameterisation, which has been shown to describe the Ly$\alpha$ 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$\alpha$ emerging from the $z>6$ ISM are required to improve our inferences. Whilst it is increasingly difficult to directly observe all of the emitted Ly$\alpha$ from $z>6$ galaxies, because of the intervening neutral gas, other emission lines could be used as a diagnostic of emerging Ly$\alpha$. For example, Henry et al. (2018) showed that Mgii emission line profiles and escape fractions closely trace those of Ly$\alpha$ 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$\alpha$, observations of Mgii at $z>6$ (which will be possible with JWST) could be used infer the nature of Ly$\alpha$ emission at these redshifts.

Additionally, better knowledge of Ly$\alpha$ line profiles at z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,5 are necessary to provide more informative priors on the observed FWHM for our inferences. In particular, high resolution spectroscopy ($R>4000$) 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$\alpha$ profile (Pentericci et al., 2018) and the prevalence of double-peaked Ly$\alpha$ 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 $z\sim 6$ 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$\alpha$ candidate sample it will be possible to use the same selections for both the $z\sim 6$ reference EW distribution and the $z>6$ 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 ${\overline{x}_{\textrm{{{hi}}}}}$ will be obtained from samples with robust redshift estimates or, ideally, spectroscopic redshifts obtained from other emission lines, and deep Ly$\alpha$ EW limits. We note that the objects which constribute the most to our posterior are the objects with the highest probability ($>60\%$) of having photometric redshift at $z\sim 8$ 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 z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7 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 z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7 galaxies, ideal for this type of analysis.

Understanding the differing evolution of Ly$\alpha$ 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 $z\sim 8$ (more likely to live in low density environments) show no significant Ly$\alpha$ emission, and favour a mostly neutral IGM. However, Ly$\alpha$ has been observed in a handful of UV bright galaxies at z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7.5 (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$\alpha$ fraction for UV bright galaxies at $z\sim 8$ 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$\alpha$ velocity offsets to facilitate Ly$\alpha$ IGM transmission. Mason et al. (2018b) proposed those objects have detectable Ly$\alpha$ because they have unusually high emitted Ly$\alpha$ 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$\alpha$ 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$\alpha$ emission at $z=9.11$?

One target in our observations (known as M1149_JD1, Zheng et al., 2012; Hoag et al., 2018) was recently spectroscopically confirmed at $z=9.11$ via [Oiii]$88\,\mu$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 $p(z)$ is at $z>9$, see Figure 4), so it was not used in our reionization inference. However, Hashimoto et al. (2018) also report a tentative 4$\sigma$ detection of Ly$\alpha$ emission from this galaxy in X-shooter observations at 12271.5 Å with total line flux $(4.3\pm 1.1)\times 10^{-18}$ erg s*-1* cm*-2*. Hoag et al. (2018) also targeted this galaxy with low resolution HST grism spectroscopy, including GLASS data, which covered the Lya wavelength at $z=9.11$. While they did not claim a detection, their spectra show a $\sim 2.5\sigma$ feature at approximately the same wavelength and flux as Hashimoto et al. (2018). We examined our KMOS cube and find no evidence of a feature at this wavelength. Our median $1\sigma$ flux limit for $z>9$ Ly$\alpha$ in the cube is $>1.1\times 10^{-18}$ erg s*-1* cm*-2*, and $>0.8\times 10^{-18}$ erg s*-1* cm*-2* at 12271.5 Å.

As noted by Hashimoto et al. (2018), if their candidate line is Ly$\alpha$, it is blueshifted by $\sim 450$ km/s with respect to the [Oiii] emission. Hashimoto et al. (2018) suggest that Ly$\alpha$ photons scattered off inflowing gas, causing it to emerge blueshifted from the galaxy’s systemic velocity. Whilst blueshifts due to inflows are expected and observed for Ly$\alpha$ (e.g., Verhamme et al., 2006; Dijkstra et al., 2006; Trainor et al., 2015), at $z>6$ the IGM is opaque to emission $<1216\,$Å, thus no Ly$\alpha$ emitted bluer than its source galaxy’s systemic redshift should be transmitted through the IGM (Dayal et al., 2011; Dijkstra et al., 2011).

Observing blueshifted Ly$\alpha$ requires the galaxy to sit in a large ionized bubble (\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,500 km/s or \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,400 kpc in radius, Haiman, 2002). Alternatively, the Ly$\alpha$ emission could arise in a different component or merging companion of the [Oiii] emitting galaxy, similar to a $z=7.1$ galaxy observed by Carniani et al. (2017). The tentative emission we observe in our KMOS cube does appear spatially offset from the predicted position of the UV continuum and [Oiii] by $\sim 0\aas@@fstack{\prime\prime}4$, which could provide evidence for the multi-component/merger scenario. This may also account for our slightly lower flux measurement as the emission extends to the edge of the IFU. However, the weakness of the detection and some general astrometric uncertainty in KMOS make a thorough analysis difficult. Deeper near-IR IFU observations of this galaxy would be extremely interesting to confirm and determine the nature of the Ly$\alpha$ emission, and will be possible in the future with JWST NIRSpec.

We calculate the probability of observing Ly$\alpha$ emission from such an object in a mostly neutral IGM using the framework of M18a, which modelled $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$. Using $m_{\textrm{F160W}}=25.7$ (Zheng et al., 2012) we obtain $M_{\textsc{uv}}=19.2-2.5\log_{10}(10/\mu)$, $EW=4\pm 2$ Å for our measured flux and $EW=11\pm 3$ Å from the measurement by Hashimoto et al. (2018). Using these measurements we calculate $p(W=4\pm 2\,\textrm{\AA}\,|\,{\overline{x}_{\textrm{{{hi}}}}}>0.76,M_{\textsc{uv}}=-19.2)<0.05$, while $p(W=11\pm 3\,\textrm{\AA}\,|\,{\overline{x}_{\textrm{{{hi}}}}}>0.76,M_{\textsc{uv}}=-19.2)<0.03$. In fact, the total probability of observing Ly$\alpha$ from this galaxy with EW $>4\pm 2$ Å if ${\overline{x}_{\textrm{{{hi}}}}}>0.76$ is \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,0.5: low Ly$\alpha$ EW are expected and consistent with a mostly neutral IGM.

We note that our calculations assume the Ly$\alpha$ is emitted close to systemic velocity (i.e., assuming that the Ly$\alpha$ comes from another component). Obviously if the galaxy does sit in an ionized bubble the probability of seeing emission would be higher. But we note that assuming emission is emitted at systemic velocity the probability of detecting the emission is not negligible, and thus this detection is still consistent with a mostly neutral IGM at $z>8$.

5.3 Other UV emission lines at $z\sim 8$

With Ly$\alpha$ increasingly suppressed at $z>6$, rest-frame UV emission lines can be used to spectroscopically confirm high-redshift LBGs. These lines can also be used as diagnostics for the stellar populations and physical conditions present in these high-redshift galaxies. Our KMOS observations cover the wavelength range where NV$\lambda 1238,1242$ and Civ$\lambda 1558,1551$ can be observed, and we briefly discuss our upper limits on the EW of these lines.

NV$\lambda 1238,1242$ can arise due to stellar winds, particularly from very young stars (Shapley et al., 2003; Jones et al., 2012), or from Hii regions if powered by an AGN or radiative shocks. Of the three $z>7$ galaxies detected to-date with tentative NV emission ($S/N\sim 4$) all have been UV bright galaxies, where AGN activity could plausibly be powering NV emission (Tilvi et al., 2016; Laporte et al., 2017b; Mainali et al., 2018). In our KLASS $7.2<z_{\mathrm{phot}}<8.8$ sub-sample (Section 4.2), the median NV EW upper limit is $<60$ Å. As our sample comprises intrinsically faint galaxies, which are less likely to have strong AGN activity, it is not surprising we do not detect strong NV emission.

Nebular Civ$\lambda 1558,1551$ emission has been observed in two galaxies at $z>6$ (Stark et al., 2015; Schmidt et al., 2017; Mainali et al., 2017). We observed Civ in the $z=6.11$ galaxy with KMOS and describe our observations in more detail in Appendix A. The Civ emission can be powered by either AGN activity or extremely metal poor stars. Limits on other UV lines in these objects find low metallicity stars are a more likely source of the hard photons needed to produce Civ emission, rather than AGN. The two galaxies are also both UV faint galaxies (M_{\textsc{uv}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,-20.2) and Mainali et al. (2018) has suggested that there is anti-correlation between UV luminosity and Civ EW, which could arise if the lowest luminosity (mass) systems are more metal-poor.

Our KLASS observation provide a large additional sample of UV faint galaxies which can place new limits on Civ emission. In our KLASS $7.2<z_{\mathrm{phot}}<8.8$ sub-sample (Section 4.2), the median Civ EW upper limit is $<74$ Å. In the three most UV faint systems with $P(7.2<z_{\mathrm{phot}}<8.8)>0.6$, M0416_22746, RXJ1347_656, M0416_1997 (all with $M_{\textsc{uv}}\sim-17.5$), the Civ upper limits are $<62$ Å, $<22$ Å, and $<51$ Å respectively. These upper limits are comparable to, and in one case below, the Civ detection presented by Stark et al. (2015) in a $M_{\textsc{uv}}\sim-19$ galaxy (with $EW_{\textsc{civ}}\approx 38$ Å), and so suggest that the proposed anti-correlation between UV luminosity and Civ EW may not be so simple.

5.4 Background limited observations with KMOS

Optical and near-IR IFU observations have provided revolutionary 3D information about the structure and kinematics of galaxies out to $z\sim 2$ (Förster Schreiber et al., 2009; Epinat et al., 2009; Wisnioski et al., 2015; Stott et al., 2016; Genzel et al., 2017) and revealed diffuse Ly$\alpha$ halos around z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,6 galaxies (Bacon et al., 2014; Karman et al., 2016; Wisotzki et al., 2016; Leclercq et al., 2017).

In KLASS we have provided the first deep NIR IFU observations of z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7 galaxy candidates. Whilst we did not make any $5\sigma$ detections of Ly$\alpha$ it is important to understand how this depended not only on the selection of our targets and the opacity of the IGM to Ly$\alpha$ at $z\sim 8$, but on the sensitivity of KMOS. In our long integrations we have pushed KMOS to the limits of its sensitivity to search for faint emission lines in near-IR IFU cubes, in wavelength regions dominated by OH sky emission lines. Using our deep observations we provide an assessment of the performance of KMOS for background-limited observations.

As described in Section 2.4 we performed additional sky subtraction routines after running the ESO pipeline to reduce residuals around bright OH lines. We also found the pipeline underestimated the noise in cubes by a factor $\sim 1.2$ and performed additional rescaling of the noise as a function of wavelength using the RMS noise obtained from the flux cubes.

One key question is how well the instrument performs relative to the predictions based on its instrumental capabilities. We can compare S/N estimated by the KMOS ETC555https://www.eso.org/observing/etc/bin/gen/form?INS.NAME=KMOS+INS.MODE=lwspectr to our achieved S/N to assess its performance. We take the $5\sigma$ flux density limits as a function of wavelength for our deepest exposure, 11 hours in RXJ1347, (shown in Figure 2) and calculate the S/N as estimated by the ETC. We use our flux calibration based on observations of standard stars to convert flux to e-/s and rescaled by the wavelength-dependent sky transmission and KMOS throughput curve (both obtained through the KMOS ETC webpage). We use the following ETC settings which are comparable to those of our observations: line FWHM${}_{\textrm{spec}}=4$ Å (unresolved); point source; seeing $0\aas@@fstack{\prime\prime}6$; airmass: 1.50, Moon illumination FLI: 0.50, Moon-target separation: 45 degrees, PWV: $<2.5$ mm. We calculate the S/N in an aperture with radius equal to the seeing FWHM $\sim 0\aas@@fstack{\prime\prime}6$.

At every wavelength, the estimated S/N is:

[TABLE]

where for RXJ1347 NDIT $=88$ is the number of DITs, of length DIT $=450$ seconds. The KMOS dark current (DC) is $0.01$ e*-/pixel/s and the read-out noise is $3.5$ e-/pixel/DIT. The aperture corresponds to $n_{\textrm{spat}}=25$ spatial pixels and the calculation is done at the peak wavelength pixel. We use the online ETC to generate the background flux $S_{\textrm{bkg}}$ in e-*/DIT as a function of wavelength, convolved with the instrumental resolution, given our input settings described above. We then calculate the estimated S/N using Equation 8 at every wavelength using our 5$\sigma$ flux density limits as the source flux.

In Figure 8 we show a comparison of the ETC estimated S/N as a function of wavelength for the line fluxes corresponding to our $5\sigma$ limits. We plot the S/N estimated by the pipeline divided by 5 to show how the achieved S/N compares to the predicted S/N from the ETC. The public ETC does not account for noise due to sky subtraction routines. Assuming all DITs have equal noise $\sigma$, for ‘A-B’ frames the noise should be $\sqrt{2}\sigma$. Thus in Figure 8 we also divide the ETC estimate by a factor $\sqrt{2}$ for a fairer comparison with our data. We find that the ETC S/N is a median $\sim 1.4\times$ higher than our achieved values, and this overestimate is highest for wavelengths \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,11500 Å, where the ETC estimate can be $\sim 1.6-1.8\times$ higher.

As shown in Figure 8, the KMOS YJ throughput is known to decrease at \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,11500 Å but our results suggests that the YJ grating is less sensitive in the blue for background-limited observations than expected.

Unfortunately this corresponds to Ly$\alpha$ redshifts z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,8.5, where we expect to find the majority of our targets. Using the S/N estimated from the ETC in planning our observations likely led us to overestimate the line sensitivity of KMOS for our targets. Most of the GLASS Ly$\alpha$ candidates we assigned to KMOS IFUs had tentative detections in the HST grisms. Thus a key aim of the deeper KMOS observations was to confirm these emission lines. While our deepest $1\sigma$ flux limit in our KMOS sample is $0.8\times 10^{-18}$ erg s*-1* cm*-2*, deeper than the $1\sigma$ flux limit in GLASS ($5\times 10^{-18}$ erg s*-1* cm*-2*), we did not detect any emission from the tentative GLASS Ly$\alpha$ candidates with KMOS, suggesting that some of the HST grism lines were spurious noise fluctuations. A thorough comparison of the GLASS and KLASS observations, in combination with other follow-up at Keck, to determine the HST grism purity and completeness will be discussed in a future paper.

We advise any future KMOS users planning observations of faint targets to take into consideration both the additional noise from sky subtraction when using the KMOS ETC, and the lower than expected performance at the blue end of YJ. However, we find better agreement with the ETC estimates at redder wavelengths, demonstrating that KMOS YJ is performing well at \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,11500 Å.

6 Summary and Conclusions

We have presented an analysis of reionization epoch targets from KLASS, a large ground-based ESO VLT/KMOS program following up sources studied in the HST grism survey GLASS. Our main conclusions are as follows:

The median $5\sigma$ flux limit of our survey is $4.5\times 10^{-18}$ erg s*-1* cm*-2*. We determine our spectroscopic survey to be 80% complete over the full wavelength range for 7.2\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,10.1 spatially unresolved Ly$\alpha$ emission lines with flux \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,5.7\times 10^{-18} erg s*-1* cm*-2*, centred within $0\aas@@fstack{\prime\prime}8$ of the IFU centre and with intrinsic line FWHM{}_{\textrm{spec}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,250 km s*-1*. Our observations are more complete to Ly$\alpha$ emission that may be spatially offset and/or extended compared to the UV continuum than typical slit-spectroscopy surveys. 2. 2.

Of the 52 z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,7 candidate targets observed, none have confirmed Ly$\alpha$ emission, including those with candidate lines detected in the HST grisms. No other UV emission lines are detected at $z>7$. We detect Civ emission in one image of a previously known Civ emitter at $z=6.11$. 3. 3.

We define a sub-sample of 29 targets with a homogeneous photometric selection of $7.2<z_{\textrm{phot}}<8.8$ for a Bayesian inference of the IGM neutral hydrogen fraction. The median Ly$\alpha$ flux limit for our sample is $3.6\times 10^{-18}$ erg s*-1* cm*-2* and the median Ly$\alpha$ EW upper limit is $58$ Å. Combining our sub-sample with 8 previously observed $z\sim 8$ LBGs from the BoRG survey (Trenti et al., 2011; Treu et al., 2013; Schmidt et al., 2014b) we obtain a lower limit on the IGM neutral hydrogen fraction at $z=7.9\pm 0.6$, ${\overline{x}_{\textrm{{{hi}}}}}>0.76$ (68%) and ${\overline{x}_{\textrm{{{hi}}}}}>0.46$ (95%). 4. 4.

Our constraint favours a late reionization consistent with models where ultra-faint galaxies contribute significantly to reionization, with an ionizing photon escape fraction \langle f_{\textrm{esc}}\rangle\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,15\%.

Our KMOS observations provide more evidence of a predominantly neutral IGM at $z\sim 8$. To make more precise constraints on the timeline of reionization will require larger samples of LBGs with precise photometric (or even better, spectroscopic) redshift estimates, more informative priors on Ly$\alpha$ FWHM, and deep spectroscopic limits on Ly$\alpha$. Forthcoming deep spectroscopic observations with JWST (e.g., Treu et al., 2017) will provide ideal samples for future inferences on reionization.

Acknowledgements

The authors thank the referee for their constructive comments. We thank Trevor Mendel and Owen Turner for useful discussions related to KMOS reductions for faint sources, and T. Mendel for sharing the readout correction code. We thank Andrei Mesinger for providing Ly$\alpha$ optical depths from the EoS simulations.

C.M. acknowledges support by NASA Headquarters through the NASA Earth and Space Science Fellowship Program Grant NNX16AO85H, and through the NASA Hubble Fellowship grant HST-HF2-51413.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This work was supported by the HST GLASS grant GO-13459, the HST BoRG grants GO-12572, 12905, 13767 and 15212, and HST-AR-13235 and HST-AR-14280.

This work was based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program 196.A-0778; and on observations made with the NASA/ESA Hubble Space Telescope, obtained at STScI. We are very grateful to all ESO and STScI staff who have assisted in these observations.

This work utilises gravitational lensing models produced by PIs Bradač, Natarajan & Kneib (CATS), Merten & Zitrin, Sharon, Williams, Keeton, Bernstein and Diego, and the GLAFIC group. This lens modelling was partially funded by the HST Frontier Fields program conducted by STScI. The lens models were obtained from the Mikulski Archive for Space Telescopes (MAST).

Software: IPython (Pérez & Granger, 2007), matplotlib (Hunter, 2007), NumPy (Van Der Walt et al., 2011), SciPy (Oliphant, 2007), Astropy (Robitaille et al., 2013), QFitsView (http://www.mpe.mpg.de/~ott/QFitsView/).

Appendix A Independent confirmation of CIV emission at $z=6.11$

As well as Ly$\alpha$ candidates we also targeted 3 images of a multiply-imaged $z=6.11$ galaxy in RXJ2248 to observe Civ$\lambda$1548,1551 emission. This multiple-image system has been detected with Ly$\alpha$ emission by Boone et al. (2013); Balestra et al. (2013); Monna et al. (2014); Karman et al. (2015); Schmidt et al. (2016) and Mainali et al. (2017). Detections of Civ and Oiii]$\lambda$1666 emission in one of the images were also presented by Mainali et al. (2017), and Civ detections and limits in all of the images by Schmidt et al. (2017). The presence of these highly ionised lines and lack of observed Heii$\lambda$1640 indicate this object has a hard ionizing spectrum, but unlikely to be dominated by an AGN.

In Figure 9 we show our GLASS and KMOS spectra for the brightest image E. An emission line is seen in the KMOS data at the same wavelength as Mainali et al. (2017), and appears to be spatially compact with size approximately equal to our seeing ($\sim 0\aas@@fstack{\prime\prime}6$). However, due to the lower spectral resolution of our data ($R\sim 3400$ compared to $R\sim 6000$ in Mainali et al., 2017) the emission line is overlapping with an adjacent skyline, adding some uncertainty to the extracted parameters.

We fit a model with 2 Gaussian emission lines to an extracted 1D spectrum, weighted by the inverse variance, fixing FWHM${}_{\textrm{spec}}=5$ Å (close to the instrumental resolution, 4 Å) and allowing flux and the wavelength position of the doublet centre to vary. We note that the spatial aperture we use to extract the spectrum (1.2 sq arcsec) is slightly larger than the one used to plot Figure 9 (0.8 sq arcsec), to ensure we recover the full flux. We use the smaller aperture in the plot to maximise the plotted S/N and for comparison with the line identification technique described in Section 3.1 which also uses 0.8 sq arcsec apertures.

The strongest peak is fit at 11023.7 Å which Mainali et al. (2017) assigned to the 1551 Å emitted peak. There is a weaker peak at 11002.5 Å which would corresponds to the 1548 Å emitted peak. We note that in the 1D spectrum there does appear to be a peak redward of the strong peak, around 11040 Å which could alternatively be the 1551 Å emitted peak. However, the peak separation to the 11040 Å line is too low for this to be part of the doublet and by inspection of the 2D emission postage stamps in the KMOS cubes we see the 11002.5 Å flux is spatially coincident with the 11023.7 Å emission, whilst the 11040 Å flux is more spatially uniform and thus likely to be spurious noise peak. Thus our assignment of the doublet wavelengths is consistent with the observations of Mainali et al. (2017).

Using our model we find a total line flux of $(1.6\pm 0.3)\times 10^{-17}$ erg s*-1* cm*-2*. We find a flux ratio of Civ$\lambda$1548:Civ$\lambda$1551 = 0.7:1 which is much lower than the theoretically motivated expected value of 2:1 (e.g., Flower et al., 1979). Low mass metal-poor galaxies at $z\sim 2-3$ have been observed with flux ratios both comparable to the theoretical value (Vanzella et al., 2016; Caminha et al., 2016) and closer to a 1:1 ratio (Christensen et al., 2012; Stark et al., 2014).

However, the confusion with the adjacent skylines makes it difficult to measure an accurate flux. Our measured total flux is lower than those measured by Mainali et al. (2017) and Schmidt et al. (2017), but consistent with the value measured by Schmidt et al. (2017) within $2\sigma$. We also find a tentative (S/N $\sim 3$) detection at the wavelength where Oiii]$\lambda$1666 was identified by Mainali et al. (2017) in our data (11837.1 Å).

Appendix B Testing the selection function of our final sub-sample

As described in Section 4.2 we define a selection function based on the photometric redshift distributions to select a sub-sample of KLASS targets to use in the reionization inference. In this appendix we demonstrate the the sub-sample was not a biased selection from the final parent catalogues.

Objects in the sub-sample must have a match in our final deep photometric catalogues: for A2744, M0416 and M1149 we used the ASTRODEEP catalogues (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. And the objects must have $P(7.2\geq z_{\mathrm{phot}}\geq 8.8)>0.01$ based on photometric redshift probability distributions calculated using EAzY (Brammer et al., 2008). We use the full photometric redshift distribution in our inference to robustly weight objects based on their probability of being in our redshift range of interest (Section 4.1 and Appendix C).

However, given that we construct this sub-sample after the observations were taken, we must check that the objects we observed were not a biased sample from the final catalogues. There are many more objects in the final catalogues which were not observed in KMOS so it could be that the KLASS targets are a biased sample of the final catalogue.

In Figure 10 we show the distribution of $P(7.2\geq z_{\mathrm{phot}}\geq 8.8)$ values of our KLASS sub-sample and the parent catalogues. The parent catalogues have many more objects with low $P(7.2\geq z_{\mathrm{phot}}\geq 8.8)$, showing that our sub-sample is skewed towards objects which are most likely to be within that redshift range, demonstrating that our initial $z_{\mathrm{phot}}>7.2$ target selection was good. In the right panel we show the probability distribution of objects in the parent catalogues and our sub-sample, using only objects with $P(7.2\geq z_{\mathrm{phot}}\geq 8.8)>0.6$. As noted previously, we do not expect all the LBGs to have $P(7.2\geq z_{\mathrm{phot}}\geq 8.8)=1$ as there are often degeneracies in the photometry that make $z_{\mathrm{phot}}\sim 1-2$ solutions possible. Here we see that the distributions are consistent and demonstrate that the best targets in KLASS sub-sample (which weight most in the reionization inference) are drawn randomly from the best targets in the parent catalogues.

We also want to be sure our sub-sample is drawn randomly from the parent sample in terms of photometry. In particular, if we had selected only the reddest objects for our KMOS targets, our sub-sample could be biased towards dustier galaxies. As dust can significantly attenuate Ly$\alpha$ emission (Hayes et al., 2011) this would mean our non-detections of Ly$\alpha$ could be due to stronger dust attenuation as well as reionization. However, we demonstrate in Figure 11 that the $\mathrm{F125W}-\mathrm{F160W}$ colour distribution and the UV slope ($\beta$, where flux $f_{\lambda}\propto\lambda^{\beta}$) distribution of our sub-sample is consistent with that of the parent catalogues. We measure UV slopes for the GLASS and KLASS objects by fitting to the $\mathrm{F125W}$ and $\mathrm{F160W}$ magnitudes (Castellano et al., 2012). The parent catalogues contain a small fraction ($<10\%$) of redder objects that are not present in our KLASS sub-sample, but given the small number of objects in our sub-sample missing this small fraction of redder objects is expected.

In addition, in the right panel of Figure 11 we compare the UV slope distribution of the $z\sim 6$ EW calibration sources used for the inference EW models (De Barros et al., 2017) to our GLASS and KLASS catalogues. Again, we want to be sure that the (De Barros et al., 2017) sample and KLASS sample are similar as we assume that these $z\sim 6$ sources are a good proxy for the $z\sim 8$ galaxies in our inference. We plot the UV slope distributions and the median errors on the $\beta$ measurements. For the KLASS $z\sim 8$ objects we can only fit the slopes using $\mathrm{F125W}$ and $\mathrm{F160W}$ magnitudes, so these measurements have large uncertainties. Using a Kolmogorov–Smirnov test we find the De Barros et al. (2017) sample, KLASS high-$z$ sample, and GLASS parent catalogue $\beta$ distributions are all consistent with being drawn from the same population. Therefore, differences in the Ly$\alpha$ transmission from $z\sim 6$ to $z\sim 8$ due to dust absorption are likely to be negligible. Spectral coverage at $>2\,\mu$m with JWST will hugely improve constraints on UV slopes for z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,6 galaxies and enable a better understanding of how dust mediates Ly$\alpha$ escape at high redshifts.

Appendix C Reionization Inference

This appendix derives the posterior distribution for the IGM neutral fraction, ${\overline{x}_{\textrm{{{hi}}}}}$, extending the framework of M18a from single EW measurements to an input flux density spectrum as a function of wavelength, and galaxy UV apparent magnitude and gravitational lensing magnification.

C.1 Likelihood at one spectral pixel

We want to obtain the likelihood of observing a flux density spectrum $\{f\}$ with flux $f_{i}$ at wavelength pixel $i$ given our the neutral fraction ${\overline{x}_{\textrm{{{hi}}}}}$ and the properties of that observed galaxy.

The likelihood of measuring flux density $f_{i}$ at wavelength $\lambda_{i}$ given that the photons originate at $z_{d}$ from a galaxy with apparent magnitude $m$ and travel through an IGM with neutral fraction ${\overline{x}_{\textrm{{{hi}}}}}$ is

[TABLE]

Including Gaussian errors in the spectra (with error $\sigma_{i}$ at spectral pixel $i$), the probability of measuring flux density $f_{i}$ at spectral pixel $i$ is given by a Gaussian distribution at each spectral pixel, with mean given by the model flux density for a given equivalent width and standard deviation $\sigma_{i}$:

[TABLE]

The total flux of the model emission line is given by $F_{tot}=Wf_{cont}(1+z)$. For simplicity, we model emission lines as Gaussians, so that the flux density produced at a single spectral pixel $i$ by an emission line at pixel $d$ (at wavelength $\lambda_{d}=\lambda_{\alpha}(1+z_{d})$) is:

[TABLE]

where $m$ is the observed apparent magnitude of the source, and $\sigma_{\lambda}=\mathrm{FWHM}/2.355$ is the spectral linewidth.

The strength of our inferred limit on the neutral fraction will depend on the choice of linewidth, as that determines the EW sensitivity. We discuss our choice of these values below in Appendix C.4.

The second term on the right-hand-side of Equation 9 can be expanded as:

[TABLE]

where the integral convolves the simulated $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$ from M18a with the probability distribution of the absolute UV magnitude, $M_{\textsc{uv}}$, given our observed data:

[TABLE]

where $M_{\textsc{uv},\textrm{mod}}(\mu,z_{d})=m-5\log_{10}({D_{L}/10\,\textrm{pc}})+2.5\log_{10}(1+z_{d})+2.5\log_{10}{\mu}$ converts observed magnitudes to rest-frame UV magnitudes, assuming the UV spectral slope $\beta=-2$ to calculate the K-correction (e.g., Blanton & Roweis, 2007). We assume the magnification distribution is log-normally distributed such that we can easily add the uncertainties in magnification and apparent magnitude: $\sigma_{M}^{2}=\sigma_{m}^{2}+(2.5\sigma_{\log{\mu}})^{2}$, where $\sigma_{m}$ is the error on the observed apparent magnitude and $\sigma_{\log{\mu}}$ is the uncertainty in the logarithmic magnifications.

We note that the dependence on $M_{\textsc{uv}}$ in $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},M_{\textsc{uv}})$ is weak compared to the dependence on ${\overline{x}_{\textrm{{{hi}}}}}$, and was parametrised with a smooth transition between two EW distributions for $M_{\textsc{uv}}>-20$ and $M_{\textsc{uv}}<-21$. Thus only the parameters which dominate changes in $M_{\textsc{uv}}$ are important. We note that the distance modulus term changes the magnitude in our redshift range of interest, $z=7.2-8.8$ ($\Delta M_{\textsc{uv}}\sim 0.5)$ less than magnification ($\Delta M_{\textsc{uv}}\sim 0.75$ for $\mu=2$), so for ease of computation we compute $p(W\,|\,{\overline{x}_{\textrm{{{hi}}}}},m,\mu,z_{d})$ ahead of time for each galaxy, setting $z_{d}=8$, rather than having to redo this integral at every spectral pixel.

C.2 Likelihood for a full spectrum

For a full spectrum $\{f\}=f(\lambda_{i})$ the likelihood is just the product of the likelihoods at each wavelength pixel:

[TABLE]

where $d$ is the index of the emission line, and $z_{d}=\lambda_{d}/\lambda_{\alpha}-1$, and $f_{\mathrm{mod,i}}$ is given by Equation 11.

C.3 Posteriors

Using Bayes’ Theorem the posterior distribution for ${\overline{x}_{\textrm{{{hi}}}}}$, FWHM and $z_{d}$ is

[TABLE]

We use a uniform prior on ${\overline{x}_{\textrm{{{hi}}}}}$ between 0 and 1, we use the photometric redshift for the prior $p(z)$, and we use a log-normal prior on FWHM (discussed in Section C.4). As we are only interested in the posterior probability of ${\overline{x}_{\textrm{{{hi}}}}}$ we can marginalise over $z_{d}$ and FWHM:

[TABLE]

To account for the incomplete wavelength coverage, we make use of the fact if the object has Ly$\alpha$ outside of the wavelength range (covering $[z_{\textrm{min}}$, $z_{\textrm{max}}]$) we would measure a non-detection in our data. Thus the integral over $z_{d}$ becomes:

[TABLE]

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

C.4 The impact of linewidth on the inference

The observed linewidth of Ly$\alpha$ emission will impact our inferences: broader lines will have lower significance and so reduce the strength of our inferences (see Figure 12. As all of our observations are non-detections we must be careful to account for this effect, which we do by marginalising our posterior over FWHM using a empirically-motivated FWHM prior.

In Figure 13 we plot $M_{\textsc{uv}}$ versus Ly$\alpha$ FWHM (deconvolved with instrumental resolution) for current detections of Ly$\alpha$ in $z>7$ LBGs (Oesch et al., 2015; Roberts-Borsani et al., 2016; Stark et al., 2017; Zitrin et al., 2015; Song et al., 2016a; Finkelstein et al., 2013; Shibuya et al., 2012; Ono et al., 2012; Schenker et al., 2012; Vanzella et al., 2011; Laporte et al., 2017a). The observed FWHM span $20-450$ km s*-1*. However linewidth measurements only exist for objects with $M_{\textsc{uv}}<-20$, significantly brighter than our KLASS sample, and there is a lot of scatter, motivating future observations to further explore the relationship.

Rather than fitting a relation to this limited current sample we construct FWHM priors using other simple empirical relations. To extrapolate to lower UV luminosities, motivated by observed correlations between Ly$\alpha$ velocity offsets and the linewidths (Verhamme et al., 2018), we also plot the ${\Delta v}(M_{\textsc{uv}})$ model derived empirically by Mason et al. (2018a) (their Equation 3). We then rescale this model using the relation between Ly$\alpha$ FWHM and ${\Delta v}$ derived by (Verhamme et al., 2018): $\mathrm{FWHM}=(34+{\Delta v})/0.9$. The Verhamme et al. (2018) relation is derived from a sample of $z\sim 0-8$ Ly$\alpha$ emitters, and suggests that the Ly$\alpha$ photon scattering in ISM which produces this correlation between line offsets and linewidth is not entirely erased by the IGM.

We use two test FWHM priors: a uniform distribution between 100 and 400 km s*-1*, spanning the currently observed range; and a log-normal distribution with mean given by the Mason et al. (2018a) ${\Delta v}(M_{\textsc{uv}})$ model scaled by the Verhamme et al. (2018) relation, with a 0.3 dex scatter. These priors are plotted in the left panel of Figure 14.

In the right panel of Figure 14 we show the posteriors for the IGM neutral fraction obtained from the full KLASS and BoRG samples by marginalising over FWHM using these two priors, as well as the case of FWHM$=300$ km s*-1* and the case of unresolved lines. We see that expect for the case of unresolved lines the difference between the posteriors is negligible. Thus our results are robust for lines with FWHM$<400$ km s*-1*. We adopt the log-normal prior for our final inference as it is the most physically motivated.

If the lines are truly very narrow ($<100$ km s*-1*) our constraint would be stronger. However, without high-resolution observations of $z>5$ Ly$\alpha$ in UV faint galaxies to provide evidence of narrow linewidths we have decided to be more conservative and use the log-normal prior in our final inferences.