Constraining the Neutral Fraction of Hydrogen in the IGM at Redshift 7.5

The spectroscopic data used in this work were obtained entirely using the MOSFIRE instrument on the Keck I telescope. The majority of the data (14 of 15 nights) come from the program "Dawn of the Galaxies: Spectroscopy of Sources at z ≳ 7" (PI Bradac; 2013 project codes: U032M, U004M; 2014: U004M; 2015: U005M, U031M; 2016: U004M; 2017: U027, U026; 2018: U005). However, we did not use five of 14 nights in our analysis because there was poor weather. In addition, we use data from one night (2016 March 20) from project code Z054M (PI M. Trenti).

We designed MOSFIRE slit masks using the publicly available MAGMA software.⁸ The names and exposure times of the targeted slit masks are listed in Table 1. In designing the masks, we assigned z ≳ 7 galaxies the highest priority, followed by z ≳ 6 galaxies (Section 2.3). We filled the remaining slits with, in decreasing order of priority, gravitationally lensed arcs, line-emitting galaxies, and cluster members, depending on the availability of preexisting spectroscopy for each cluster. There was an average of about six of the highest-priority (z ≳ 7) galaxies on each mask, with each mask containing a total of ∼30 objects.

Table 1.  MOSFIRE Slit Masks Targeted

Cluster	Mask Name	Date of Observation	texp	$\langle \mathrm{seeing}\rangle$	$\langle \mathrm{attenuation}\rangle$	$\langle \mathrm{airmass}\rangle$
		(UTC)	(sec.)	(")	(mag.)
A2744	A2744_MR	2015 Nov 7	4320	0.73	0.05	1.13
A2744	A2744_MR	2015 Nov 8	2520	1.18	0.05	1.17
A370	A370_M	2013 Dec 16a	3060	1.13	⋯	1.10
A370	A370_M	2013 Dec 17a	1800	1.03	0.05	1.12
A370	A370_M	2013 Dec 18a	3600	0.89	0.1	1.09
A370	A370_M	2017 Oct 1	14580	0.67	0.05	1.18
MACS0416	MACS0416_M	2015 Nov 7	9000	0.67	0.04	1.12
MACS0416	MACS0416_M	2015 Nov 8	7200	0.84	0.05	1.12
MACS0454	macs0454_M	2013 Dec 15a	3240	0.64	0.05	1.19
MACS0744	MACS0744_M_backup	2015 Nov 8	7200	0.60	0.08	1.08
MACS0744	MACS0744_M	2016 Feb 22	14220	1.06	0.03	1.14
MACS0744	MACS0744_M	2016 Feb 23a	13680	0.66	0.05	1.12
MACS0744	MACS0744_ATH4Yband_mt_M	2016 Mar 20	5040	0.52	0.08	1.40
MACS1149	MACS1149	2014 Feb 14	2520	0.60	0.9	1.15
MACS1149	MACS1149_rot180	2016 Feb 22	10800	1.40	⋯	1.15
MACS1423	miki14M	2013 Jun 13a	5040	0.90	0.07	1.44
MACS1423	MACS1423_M	2015 Apr 27a	7560	0.68	0.06	1.31
MACS1423	MACS1423_052715_M	2015 May 28	17280	0.69	0.06	1.20
MACS1423	MACS1423_20160319_M	2016 Mar 20	6660	0.85	0.05	1.16
MACS2129	MACS2129_M	2015 Apr 27a	1440	0.80	0.04	1.55
MACS2129	MACS2129_052715_M	2015 may 28	7560	0.55	0.1	1.27
MACS2129	MACS2129_M	2015 Nov 7	8460	0.53	0.05	1.17
MACS2129	MACS2129_M	2015 Nov 8	9540	0.84	0.04	1.20
MACS2214	MACS2214_M	2017 Oct 1	14400	0.84	0.04	1.16
RCS2327	rcs22327_M	2013 Dec 15a	3420	0.51	0.06	1.13
RCS2327	rcs22327_M	2013 Dec 17a	3600	0.71	0.07	1.13
RCS2327	rcs22327_M_1	2013 Dec 18a	3600	1.14	0.08	1.13
RXJ1347	RXJ1347_v4	2018 Jun 1	11700	0.60	0.18	1.13

Note.

^aIndicates a half night of observation. Attenuation values marked by "⋯" mean that attenuation data were not available on these dates.

Download table as:  ASCIITypeset image

All MOSFIRE data were obtained with 07-wide slits in the Y band, which covers Lyα over the redshift range 7 ≲ z ≲ 8.2. We observed with an ABBA dither pattern, using 125 amplitude dithers, which results in a spacing of 25 between the positive signal and each negative shadow. The average FWHM of the atmospheric seeing was generally  ≲1'', and the average attenuation of the data we used (measured in the V band) was generally small (≲0.05 mag). Seeing was either measured directly from stars on the science slit mask or from prescience alignment images taken in the J band.

We reduced the MOSFIRE data using the publicly available data reduction pipeline (DRP).⁹ The DRP produces wavelength-calibrated, rectified, background- and skyline-subtracted 2D signal and noise spectra for each slit on the mask. We extracted 1D signal and noise spectra from the 2D spectra for each slit, using uniform weights on all spatial pixels for the extraction. The spatial aperture we used for extraction was twice the FWHM of the atmospheric seeing. We found that for the majority of our observations, the pipeline-produced noise spectra underestimate the uncertainty derived directly from signal spectra in regions where no objects or sky emission was present. To remedy this, we scaled all 1D noise spectra so that the signal-to-noise ratio (S/N) spectra in regions free from object traces and skyline emission had a distribution with standard deviation equal to 1. We then divided the 1D signal spectra by the scaled 1D noise spectra to obtain properly behaved 1D S/N spectra for each slit on each slit mask. We opted to use the rescaled noise spectrum rather than, for example, the standard deviation within the flux spectrum itself because of potential contamination in the flux spectrum from nearby objects.

Flux calibration of the 1D signal and noise spectra was performed following the steps outlined by Hoag et al. (2015). For all masks observed before 2016 September, we used a spectral type AOV star we observed on 2013 December 15 to flux calibrate our spectra. MOSFIRE was repaired during the interval between 2016 September and 2017 February, potentially altering the response function of the instrument. Therefore, we used an AOV star observed after the repair (2017 October 1) to flux calibrate the masks observed after the repair. We account for differences in the extinction and air mass of the calibration star compared to each mask observation. This is done by dereddening and applying an air-mass correction using the average air-mass value during observations to both the calibration star and each science spectrum before applying the telescope correction to the science spectra. This step effectively applies a slit-loss correction to each science spectrum, assuming our targets are all point sources. Most sources are indeed not extended, so we do not apply an extended-source slit-loss term to our spectra. While Lyα is known to be on average more extended than the rest-frame UV continuum (e.g., Steidel et al. 2011; Wisotzki et al. 2016; Leclercq et al. 2017), an analysis with the KMOS IFU down to flux limits similar to that obtained in this work (Mason et al. 2019) found no evidence that would suggest we are missing Lyα because of these effects. These effects may still, however, decrease our sensitivity, so our flux limits presented in Section 3.6 are effectively lower limits.

Many of the same objects were observed on multiple masks. In these cases, we stacked the calibrated 1D signal and noise spectra from all available individual masks using a simple mean and produced full-depth S/N spectra from these stacks. While masks with somewhat poor conditions were included, this did not result in the loss of candidate emission lines in the final stacks; we inspected all individual spectra for emission lines in addition to our search described in Section 3.3.

The imaging data were primarily used to measure photometric redshifts to select a sample of high-z galaxies for spectroscopic follow-up. The images were also used in constructing the lens models (Section 3.2). The data come from several Hubble Space Telescope (HST) and Spitzer/IRAC programs, summarized in Table 2. HST mosaics of the clusters used in this work are shown in Figures 1 and 2. We targeted several clusters belonging to the Hubble Frontier Fields Initiative (HFF; Lotz et al. 2017). These are Abell 2744 (A2744), M0416.1-2403 (MACS0416), MACSJ0717.5+3745 (MACS0717), MACSJ1149.5+2223 (MACS1149), and Abell 370 (A370). As explained in Section 3.1, MACS0717 had zero targets in our final photometric selection, so it is not included in the rest of this work. For these clusters, we also use Keck/MOSFIRE K_s-band photometry from the K-band Imaging of the Frontier Fields program (Brammer et al. 2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another four clusters, MACS0744, MACS1423, MACS2129, and RXJ1347, are part of the Cluster Lensing And Supernova survey with Hubble (CLASH; Postman et al. 2012). MACS0454 was observed by HST-GO-11591/GO-9836/GO-9722, and the HST images for both MACS2214 and RCS2327 were obtained as part the Spitzer UltRa Faint SUrvey Program (SURFSUP; Bradač et al. 2014). Many of these clusters were observed by the Grism Lens-Amplified Survey from Space (GLASS; Schmidt et al. 2014, 2016; Treu et al. 2015), which obtained HST WFC3/IR imaging and grism spectroscopy of the prime cluster fields. Our target selection was mildly influenced by the Lyα candidates assembled by Schmidt et al. (2016; see Section 3.1).

The depths of the HST images vary among the clusters. For example, the images of the HFF clusters that we targeted all reach limiting magnitudes of ∼29 AB mag (5σ, point source) per band, whereas the remaining clusters are ∼2 mag shallower per band. The images of the HFF and CLASH clusters were obtained from the latest public mosaics available on the Mikulski Archive for Space Telescopes (MAST).^{10 ,11} The HFF data span three Advanced Camera for Surveys (ACS) filters, F435W, F606W, and F814W, and four Wide Field Camera 3 (WFC3) IR filters: F105W, F125W, F140W, and F160W. The CLASH clusters share these filters but have five additional filters that we use: F475W, F625W, F775W, F850LP (ACS), and F110W (WFC3).¹² The three non-CLASH and non-HFF clusters, MACS0454, MACS2214, and RCS2327, were processed as described by Huang et al. (2016a) and have a subset of the CLASH filters.

We also use deep Spitzer/IRAC data in tandem with the HST data. The primary purposes of these data are to more accurately select and characterize the high-z sample (e.g., Huang et al. 2016a; V. Strait et al. 2019, in preparation). The IRAC images were obtained by the SURFSUP (PID 90009; Bradač et al. 2014) and the SFF (PIDs: 83, 137, 10171, 40652, 60034, 80168, 90257, 90258, 90259, 90260). The data have similar depths in all 12 clusters, reaching at least ∼30 hr in the [3.6] μm and [4.5] μm bands (hereafter [3.6] and [4.5]). The processing of the Spitzer/IRAC data is described by Ryan et al. (2014) and Huang et al. (2016a).

While the HST photometry procedure is described in detail by Huang et al. (2016a), we briefly outline the primary steps we take to produce photometric catalogs. First, we produce point-spread function (PSF) matched HST images with a 60 mas pixel size, convolving all images (except F160W) to match the resolution of F160W, the lowest-resolution HST image. We produce a stacked NIR image from the HST WFC3 images to use for object detection. We then run Source Extractor (SExtractor; Bertin & Arnouts 1996) in dual-image mode on each filter using the stacked WFC3 NIR image as the detection image. We correct each filter for galactic extinction using the IR dust maps obtained by Schlegel et al. (1998) and the coefficients in each filter from Postman et al. (2012). We also find from simulations that the flux errors in each filter reported by SExtractor are underestimated. We correct the errors via source simulations using a method similar to Trenti et al. (2011).

For some of the clusters, that is, the HFF clusters, MACS1423, MACS2129, and RXJ1347, we also performed a subtraction of the intracluster light (ICL) and bright cluster members using the method developed by the ASTRODEEP collaboration (Castellano et al. 2016; Merlin et al. 2016). This was done to obtain more accurate photometry and to improve number counts of lensed, high-redshift sources. It is important to perform it on the deeper images, for example, the HFF images, because of the stronger contamination from the ICL. However, the ICL is less of an issue for the CLASH-depth clusters, so the fact that this was not performed on MACS0454, MACS0744, and RCS2327 is therefore of little concern. In particular, the sources that would be unveiled in these clusters from ICL and cluster member subtraction would be far too contaminated in the MOSFIRE slits to obtain constraining Lyα flux limits.

The IRAC photometry is described by Huang et al. (2016a). In brief, we use T-PHOT (Merlin et al. 2015) to measure colors between HST and IRAC. We use priors based on the segmentation map from the HST F160W image, the image closest in wavelength to the IRAC bands, to deblend objects in the lower-resolution IRAC bands. In this way, we can assemble a consistent photometric catalog with IRAC for all of the sources identified in the HST NIR detection image. The end result of the process is a merged HST and IRAC photometric catalog for all detected sources in the field.

The selection criteria for including high-redshift galaxies on our slit masks was very inclusive. Where possible, we incorporated multiple selections from the literature to maximize the number of high-z targets on the mask, using the compilation by Schmidt et al. (2016). Since our observations began in 2013 December, we have obtained significantly improved HST and Spitzer data for many of the clusters in this sample. In particular, the HFF program and much of the SURFSUP and GLASS programs were performed during this time.

In order to use a more consistent photometric selection across all clusters in this work and to make use of the much deeper current data, we performed the same photometric processing of all of the clusters in our sample after all of our spectroscopic observations were taken. We recalculated the photometry for each cluster using the pipeline described in Section 2.3 and refit all of the photometry to obtain photometric redshift probability distributions, P(z)'s, and spectral energy distributions (SEDs) using EAzY (Brammer et al. 2008), adopting the default v1.3 EAzY spectral templates. As in Huang et al. (2016a), we inspected the distribution of best-fit photo-z's for all clusters, ensuring that the distribution peaked at the cluster redshift. For the clusters in the GLASS sample (see Table 2), we were able to anchor the EAzY photo-z distribution by comparing our photo-z's to spectroscopic redshifts of intermediate-z galaxies in the GLASS data set. We found good statistical agreement in all clusters (e.g., Wang et al. 2015; Hoag et al. 2016). This was not possible for the non-GLASS clusters, MACS0454, MACS2214, and RCS2327. These clusters hold less weight in the neutral fraction inference because they generally have fewer objects and shallower observations. Nonetheless, their photo-z's are potentially less reliable than the photo-z's of targets in the GLASS clusters.

We used a flat prior on z, rather than adopting the default magnitude prior when deriving the P(z)'s of all objects with EAzY. We did not adopt the default EAzY magnitude prior because our targets are lensed. From the P(z) of each galaxy, we calculated the probability that Lyα falls in the MOSFIRE Y band, that is, p(7 < z < 8.2). We then selected objects with p(7 < z < 8.2) > 0.01, finding a total of 2828 such objects in all 11 clusters (excluding MACS0717; as explained below). A significant fraction of these were spurious. We cleaned the catalog of objects that were obvious parts of larger galaxies, diffraction spikes, bad pixels, or clearly originated from the edge effects at the boundary of the detector or overlap between epochs. This inspection was performed visually, so it is naturally subjective. However, it is more reliable than our current automatic methods to remove such artifacts. After cleaning the sample, we arrived at a sample of 514 objects that fulfill the p(7 < z < 8.2) > 0.01 condition. We note that the 5σ limiting magnitudes of our parent photometric catalog are ∼28.5 mag (HFF) and ∼26.5 mag (non-HFF), about 0.5 mag brighter than the quoted limiting magnitudes in these bands.

Of these 514 objects, we observed 70 in our Keck/MOSFIRE survey. Two of these objects were spectroscopically confirmed to be at lower redshift (6 < z < 7) by Huang et al. (2016b) and S. Fuller et al. (2019, in preparation) with Keck/DEIMOS during a simultaneous spectroscopic campaign. We note that the spectroscopic redshifts of both galaxies are consistent with the photometric redshifts we obtained for those objects with EAzY. The remaining 68 objects will hereafter be referred to as the MOSFIRE sample. We note that while we observed several candidates from our original photometric catalogs in MACS0717 in 2013 December based on pre-HFF HST and Spitzer/IRAC imaging, the final photometric selection from our full-depth images yielded zero candidates that we had targeted with MOSFIRE. As a result, MACS0717 is not part of the MOSFIRE sample and is not discussed in the rest of this work.

As a sanity check of our LBG selection, we compared the total number of galaxies selected to be in the MOSFIRE range 7 < z < 8.2 with the expected number of galaxies in this volume from two different LF models (Bouwens et al. 2015b; Mason et al. 2015), correcting for the magnification in each cluster field separately. The actual number of sources in our LBG selection is 61, which we estimated by summing the integrated probabilities p(7 < z < 8.2) over all objects in our LBG sample. The predicted numbers from the LF in this same redshift interval are 74 and 73 from the Mason et al. (2015) and Bouwens et al. (2015b) LFs, respectively. The fact that the LF models predict slightly larger numbers most likely reflects some incompleteness in our LBG selection, which is expected based on the crowding in cluster fields. While the photometry for the HFF clusters and MACS1423, MACS2129, and RXJ1347 included some cluster member and ICL subtraction, completeness of 100% is never achieved, especially for faint sources. The fact that we selected a similar number of sources to the LF predictions indicates that the influence of low-z contaminants on our final neutral fraction result is small.

The quantity p(7 < z < 8.2) is used both for the sample selection and as an effective weight in the inference (described in Section 3.7). We explored whether using a different photo-z code would affect our LBG selection and hence the neutral fraction inference. We used the ASTRODEEP team's OAR photo-z's (Fontana et al. 2000) to recalculate the quantity p(7 < z < 8.2) for a subset of our sample (A2744, MACS0416, MACS1149) This quantity is in general not identical for a given object between OAR and EAzY, which is expected because of the differences in the code and templates used, and so on. However, we expect that the distribution of differences across many objects should not be biased. To test this, we compared the difference in this quantity for all objects in the ASTRODEEP photometric catalog for the three clusters for which we have OAR and EAzY redshifts. We found that the distribution of differences in the effective weight factor between codes is strongly peaked at zero and is not biased, indicating that using OAR instead of EAzY would not significantly affect our inference on the neutral fraction.

During some of our observing runs, we preferentially targeted GLASS Lyα emitter candidates from Schmidt et al. (2016). While these represented typically only one or two targets per cluster, we checked that this selection did not introduce a bias into our final result. We ran the neutral fraction inference (described in Section 3.7) with and without these targets and found that the difference was insignificant.

We also compared the distributions of properties of the MOSFIRE sample with the parent photometric sample. Figure 3 shows the comparison for the F160W apparent magnitude, peak photometric redshift (z_peak), intrinsic absolute magnitude corrected for lensing (${M}_{\mathrm{UV}}-2.5{\mathrm{log}}_{10}(\mu /{\mu }_{\mathrm{best}})$), and the rest-frame UV color (measured near-IR color F125W–F160W). We performed a two-sample Kolmogorov–Smirnov test for each of these properties, finding significant differences (>99% probability) between the distributions for the F160W magnitude, z_peak, and the absolute magnitude. These differences simply illustrate our mask design strategy, which was to put brighter candidates that had a higher probability of being high redshift on the masks. Because the magnitudes of our spectroscopic sources are not near the detection limit of the HST images, the contamination fraction is likely lower than in the parent photometry catalog. F160W, z, and M_UV are parameters in the model we use to infer the neutral fraction, so they do not bias our inference (see Section 3.7). The inference is only biased if there are differences in observables that may impact Lyα that are not parameters in our model, for example the rest-UV color. In Figure 3, we show that the rest-UV color is statistically consistent between the two samples.

Zoom In Zoom Out Reset image size

We plot the P(z)'s for all 68 galaxies in the MOSFIRE sample in Figure 4. Most of the probability from the 68 galaxies falls within the MOSFIRE Y-band coverage. We account for the photometric redshift impurity in the neutral fraction inference, as explained in Section 3.7. We list the targets along with some of their properties in the Appendix. Throughout this work, we refer to individual targets via their cluster and ID in the following manner: the first target in Table 5 is ID = A370-000 and the last is ID = MACS0744-069. We note that the median absolute magnitude of our sample is M_UV = −18.25, where we account for magnification for each galaxy before computing the median. This is ∼0.1 L_⋆ at z ∼ 7–8 (Bouwens et al. 2015b), illustrating that the galaxies we are targeting are intrinsically very faint and probably more characteristic of the general galaxy population than ∼L_⋆ galaxies.

Zoom In Zoom Out Reset image size

To link the observed LBGs to halos in reionization simulations for performing the inference on the hydrogen neutral fraction (see Section 3.7), we use the galaxy intrinsic luminosity. Because we observed LBGs in lensed fields, we must correct their observed luminosities by the magnification factor. To obtain the magnification factor for each galaxy, we produce lens models of all of the clusters in our sample. Lens modeling is performed using the free-form code developed by Bradač et al. (2005, 2009). Briefly, the code solves for the gravitational potential on an adaptive pixel grid via χ² minimization. The lensing quantities such as convergence, shear, and magnification are all derived from the best-fit gravitational potential.

The lens models for all clusters were constructed using the HST imaging data described in Section 2.2. HST images are used to identify strongly lensed (multiply imaged) galaxies as well as weakly lensed (singly imaged) galaxies, both of which are used as constraints to the lens model. We also used spectroscopic redshifts from GLASS (Schmidt et al. 2014; Wang et al. 2015; Hoag et al. 2016; Treu et al. 2016), from the CLASH-VLT program (P.I.: P. Rosati; Rosati et al. 2014), and from Limousin et al. (2012), Johnson et al. (2014), Zitrin et al. (2015), Grillo et al. (2016), Monna et al. (2017), and Lagattuta et al. (2017) to improve the precision of the lens models. In general, we used only the set of multiply imaged systems with either secure spectroscopic redshifts or consistent photometric redshifts (as in, e.g., Hoag et al. 2016). To obtain uncertainties on the lensing quantities, including the magnification, we bootstrap resample the weak lensing catalog and rerun the lens modeling code 100 times. The number of weakly lensed galaxies (typically ∼500–1000) is much larger than the number of strongly lensed galaxies, so we opt to use the weak lensing catalog for the bootstraps. In cases where the number of multiply imaged systems is sufficiently large (≳15), we simultaneously resample from the strongly lensed and weakly lensed galaxy catalogs to obtain uncertainties. For more details on this procedure, see Hoag et al. (2016).

Some of the lens models used in this work are described in previous works: the RCS2327 model was obtained by Hoag et al. (2015), A2744 by Wang et al. (2015), MACS0416 by Hoag et al. (2016), MACS2129 by Huang et al. (2016b), MACS1423 by Hoag et al. (2017), MACS1149 by Finney et al. (2018), and A370 by Strait et al. (2018). Lens models for MACS0454, MACS0744, MACS2214, and RXJ1347 are presented here for the first time. We show the critical curves from all of our lens models as well as the location of the high-redshift objects targeted in the MOSFIRE campaign with nonzero probability of being at 7 < z < 8.2 in Figures 1 and 2.

After collecting all of our MOSFIRE observations and selecting the 68 targets that fulfill our high-z selection criteria, we performed an automatic search for emission lines. For each target, we calculated S/N_int, the integrated S/N, at each wavelength in the full-depth 1D S/N spectrum, assembling the integrated signal-to-noise spectrum. We used a spectral bandpass of 6 Å (about six pixels) for the integration because this is twice the FWHM spectral resolution of the instrument. For each object, we flagged any wavelengths where S/N_int ≥ 5, that is, lines with high confidence. This resulted in 22 candidate emission lines.

We visually inspected all 22 candidate lines identified by the automatic detection procedure. Real emission lines have several characteristics that can distinguish them from spurious features:

• 1.  
  Two negative residuals at symmetric positions flanking the central emission line and with approximately half of the intensity of the central emission line, due to our ABBA dither pattern
• 2.  
  At least as spectrally extended as the MOSFIRE resolution, that it, about three pixels in the Y band
• 3.  
  At least as spatially extended as the atmospheric seeing at the time of observation
• 4.  
  Distinguishable from sky emission line features.

Because these tests were performed visually, they are inherently subjective. However, the candidate list was short (22 candidates), and most of the candidates obviously failed multiple tests. There was one ambiguous case, which we describe in Section 3.4.

We identified two emission lines that passed the tests described in Section 3.3. The S/N spectra for these two objects, MACS0744-064 and RXJ1347-018, are shown in Figures 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed MACS0744-064 (Figure 5) on 2016 February 22, 2016 February 23, and 2016 March 20 with slit masks MACS0744_M, MACS0744_M, and MACS0744_ATH4Yband_mt_M, respectively (see Table 1). The emission line is significantly detected on 2016 February 23 and March 20, but hardly detected on 2016 February 22. While the exposure time was comparable for 2016 February 22 and 23, the seeing was much worse on 2016 February 22, explaining the difference in S/N. The seeing and attenuation were both good on 2016 February 22 and March 20, the two dates where the emission line is significantly detected. ID = MACS0744-064 is confidently detected, with an S/N_int = 10.3 in the full-depth stack. The two negative residuals are also clearly detected in the stack at the correct locations and can also be seen on multiple nights of observation. The line is completely separated from sky emission lines. As a result, this line is considered secure.

ID = RXJ1347-018 is detected at S/N_int = 5.1. Unlike MACS0744-064, we only observed this object during a single night of observation, so we do not have multiple nights to corroborate the detection. To test the robustness of the detection, we split the data in half, rereduced each half, and reextracted the spectra to see if the emission line was present in both halves.¹³ We found a peak in the integrated S/N in both halves at the same wavelength as in the full-depth spectrum. Using the same bandpass as in the full-depth spectrum, we found values of S/N_int = 3.82 and S/N_int = 3.11 in the two halves, consistent with a real line of S/N_int ∼ 5 in the combined data set.

Despite this evidence in support of the line, other aspects of the line are concerning. A sky emission line is present just blueward of the line-peak emission wavelength of 9924 Å. Part of the 6 Å spectral bandpass used to calculate the integrated S/N falls on this sky emission line. While the noise spectrum in principle takes into account the noise from the sky emission line, the sky subtraction is imperfect and could introduce spurious features at this level of S/N. The negative residuals are present for this object with ratios of −0.9 ± 0.3 for both top and bottom residuals. These ratios should be consistent with −0.5 as a result of our dither pattern. However, because the central emission line is only detected at S/N_int ∼ 5, the negative residuals themselves are only expected to be detected at $\left|S/{N}_{\mathrm{int}}\right|\sim 3.5$. At such low S/N, the ratios of the negative residuals to the central emission are less robust, and thus this test is less reliable. For example, the proximity to the skyline could influence the deviation of these ratios via imperfect skyline subtraction. Given that these factors reduce our confidence in the emission line, we carry out the neutral fraction inference with and without this detection and compare the results in Section 3.7, finding that its inclusion or exclusion does not change our result in a statistically significant way. Future follow-up of this cluster is needed to confirm or deny the validity of this emission line.

Both of the emission lines presented here are the only lines detected in the spectra of their respective targets. Here we provide additional evidence to properly identify them as Lyα. For both objects, our main evidence in support of Lyα is the photometric redshift distribution, P(z). We show the P(z) plots for MACS0744-064 and RXJ1347-018 in Figures 7 and 8, respectively. In both cases, the spectroscopic redshift is in good agreement with the photometric P(z), falling within the 68% confidence interval. The most common contaminant of Lyα at z ∼ 7 is the [O ii]λλ3726,3729 doublet at z ∼ 1.5. The P(z) plots suggest that [O ii] is very unlikely for both MACS0744-064 and RXJ1347-018. Furthermore, the [O ii] doublet would be resolved given the MOSFIRE resolution, with a peak separation of ∼7 Å. We do not detect any doublets with this separation in our spectra. Within the 95% confidence interval of the P(z)'s of both targets, the C iv and C iii] doublets could fall in the MOSFIRE Y band. However, these lines are typically much weaker in LBGs than the rest-frame equivalent widths that would be required to detect them given the faintness of these targets (e.g., Shapley et al. 2003; Stark et al. 2014; Le Fèvre et al. 2019). Furthermore, in some cases we would expect to see both lines in the doublet. However, we cannot completely rule out that high equivalent width C iv or C iii], such as values observed by Stark et al. (2015, 2017), could explain the emission lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The RXJ1347 emission line looks marginally double peaked, but the separation of these two features is ∼4 Å and cannot be explained by any strong emission line doublets at redshifts consistent with the P(z). The bluer peak falls on a narrow feature in the noise spectrum between a bright skyline and the center of the emission line, but which is not associated with the skyline. The feature originates from two pixels in the 2D noise spectrum, which are identified by the DRP. We tested to see whether this feature was responsible for the bluer bump in the S/N spectrum. We masked out these two pixels and reextracted the S/N spectrum, finding that the bump did not disappear and that the integrated S/N was slightly larger than in the original case. We conclude that the blue bump originates from the flux spectrum, but it could be a result of imperfect sky subtraction of the bright skyline at ∼9917 Å.

We also show the SEDs of both galaxies in Figures 7 and 8. We perform SED fitting using the method described by Huang et al. (2016a). In brief, we use EAzY with the Bruzual & Charlot (2003) stellar population synthesis models, adopting a Chabrier (2003) initial mass function between 0.1 and 100 M_⊙, a Calzetti et al. (2000) dust attenuation law, and an exponentially declining star formation history with a fixed metallicity of Z = 0.2 Z_⊙. Nebular emission lines are added to the templates as described by Huang et al. (2016a).

We find that MACS0744-064 has a considerably bluer rest-frame UV slope than RXJ1347-018. This results in a stark difference in some of their physical properties, which are provided in Tables 3 and 4. In particular, MACS0744-064 is significantly younger (age = ${17.38}_{-6.91}^{+5.53}\,\mathrm{Myr}$) than RXJ1347-018 (age = ${640.47}_{-131.73}^{+78.15}\,\mathrm{Myr}$). Younger stellar ages (∼10–100 Myr) are expected at z ∼ 7 given the relatively brief amount of time at this redshift for star formation to have taken place since the big bang (∼750 Myr). With this in mind, RXJ1347-018 is surprisingly old. However, relatively evolved stellar populations have possibly been observed at similar redshifts and beyond (e.g., Richard et al. 2011; Zheng et al. 2012; Bradač et al. 2014; Hashimoto et al. 2018; V. Strait et al. 2019, in preparation), suggesting the onset of star formation within the first ∼300 Myr after the big bang. Katz et al. (2019) showed that the SEDs of these rare, evolved stellar populations can be reproduced by cosmological simulations, but only in the most massive halos in their simulation.

Table 3.  Properties of MACS0744-064

Quantity	Unit	Value
αJ2000	(degree)	116.24648
δJ2000	(degree)	39.46042
F160W	mag	27.17 ± 0.38
[3.6]	mag	<25.54
[4.5]	mag	<25.37
μbest		3.2 ± 0.1
MUV − 2.5 log10(μ/μbest)	mag	−18.5 ± 0.4
LUV × μ/μbest	L⋆	${0.12}_{-0.04}^{+0.06}$
zphot		${6.9}_{-0.2}^{+0.5}$
zLyα		7.148 ± 0.001
SFR × μ/μbest	M⊙ yr−1	${0.66}_{-0.11}^{+0.18}$
M⋆ × μ/μbest	107M⊙	${1.11}_{-0.28}^{+0.35}$
Age	Myr	${17.38}_{-6.91}^{+5.53}$
E(B − V)	mag	<0.01
${f}_{,\mathrm{Ly}\alpha }$	10−18 erg s−1 cm−2	7.1 ± 0.7
${W}_{0,\mathrm{Ly}\alpha }$	Å	58.3 ± 25.1
FWHMLyα	Å	5.4 ± 1.3
S/Nint, Lyα		10.3

Note. IRAC [3.6] and [4.5] magnitude limits are 3σ. The factor μ/μ_best is introduced to show how the magnification factor enters into some properties, where μ_best is the best-fit magnification value measured in this work. This factor allows one to recalculate the value of a given property provided with a new magnification factor. M_UV and L_UV are calculated using the F160W magnitude to approximate the rest-frame UV luminosity. L_⋆ is calculated from ${M}_{\mathrm{UV}}^{\star }=20.87\pm 0.26$ measured by Bouwens et al. (2015b) at z = 6.8.

Download table as:  ASCIITypeset image

Table 4.  Properties of RXJ1347-018

Quantity	Unit	Value
αJ2000	(degree)	206.89124
δJ2000	(degree)	−11.75261
F160W	mag	26.43 ± 0.14
[3.6]	mag	24.47 ± 0.17
[4.5]	mag	24.07 ± 0.11
μbest		${21.4}_{-1.3}^{+1.7}$
MUV − 2.5 log10(μ/μbest)	mag	−17.2 ± 0.2
LUV × μ/μbest	L⋆	0.03 ± 0.01
zphot		${7.5}_{-0.7}^{+0.3}$
zLyα		7.161 ± 0.001
SFR × μ/μbest	M⊙ yr−1	${3.40}_{-0.59}^{+1.45}$
M⋆ × μ/μbest	107 M⊙	${143.34}_{-35.18}^{+26.38}$
Age	Myr	${640.47}_{-131.73}^{+78.15}$
${(\mathrm{SFR}\times \mu /{\mu }_{\mathrm{best}})}_{{\text{}}\mathrm{HST}}$	M⊙ yr−1	${1.59}_{-1.05}^{+3.63}$
${({M}_{\star }\times \mu /{\mu }_{\mathrm{best}})}_{{\text{}}\mathrm{HST}}$	107 M⊙	${14.53}_{-11.88}^{+41.37}$
(Age)HST	Myr	${321.00}_{-310.52}^{+397.62}$
E(B − V)	mag	<0.35
${f}_{,\mathrm{Ly}\alpha }$	10−18 erg s−1 cm−2	6.6 ± 1.3
${W}_{0,\mathrm{Ly}\alpha }$	Å	27.2 ± 6.5
FWHMLyα	Å	7.6 ± 1.3
S/Nint,Lyα		5.1

Note. Quantities with subscript HST are derived from the best-fit SED to HST photometry only, whereas the other quantities include the IRAC photometry. The uncertainty on FWHM_Lyα does not include additional uncertainty introduced by the sky emission line blueward of the emission peak.

Download table as:  ASCIITypeset image

The red UV slope of RXJ1347-018 is partially driven by the detections in the IRAC [3.6] μm and [4.5] μm bands. For example, in Figure 8 we show our best-fit SED for RXJ1347-018 using only the HST constraints. The UV slope is bluer than when we include the IRAC constraints. The IRAC detections are questionable, due to the presence of a bright neighboring galaxy. In Figure 9, we show a cutout of the HST F160W image as well as the two IRAC bands centered on RXJ1347-018. The galaxy is resolved in HST, but it is heavily blended with the neighboring galaxy in IRAC, due to the much larger PSF. While we attempted to subtract the neighbor using T-PHOT when performing photometry on RXJ1347-018, the proximity to RXJ1347-018 makes the fluxes untrustworthy. As a result, we show the SED with and without the IRAC constraints in Figure 8, and we report the physical properties of the galaxy in both cases in Table 4. Using HST photometry only, the age of the galaxy is very poorly constrained (age = ${321.00}_{-310.52}^{+397.62}\,\mathrm{Myr}$), and the stellar mass is an order of magnitude smaller than when we include the IRAC constraints.

Zoom In Zoom Out Reset image size

In some cases, age and stellar mass inferences can be biased by strong unmodeled rest-frame UV and optical emission lines falling in the near-IR HST and mid-IR IRAC filters (e.g., Schaerer & de Barros 2009; Labbé et al. 2013; Smit et al. 2014). The strongest of such emission lines are the [O ii], [O iii], Hβ, and Hα lines. At the redshift of RXJ1347-018, z = 7.161, none of these strong lines fall in [3.6], so they could not explain the bright flux in this filter even if it were real. While Hβ and [O iii] fall in the [4.5] band at z = 7.161 and could potentially explain this flux, this would not explain the flux in the [3.6] band.

Molino et al. (2017) performed photometry and SED fitting for the CLASH clusters, which include MACS0744 and RXJ1347. Both MACS0744-064 and RXJ1347-018 are present in their publicly available photometric catalogs.¹⁴ We compare to the only two fitted quantities present in their catalog, which are the photo-z and the stellar mass. For RXJ1347-018 (CLASHID = 0911), they obtained a best-fit photometric redshift of $z={7.61}_{-0.81}^{+0.13}$ (95% conf.), in agreement with the Lyα spectroscopic redshift (and the photo-z) that we measured. They report a stellar mass for RXJ1347-018 (uncertainties are not provided) of 3.89 × 10⁹ M_⊙ before accounting for magnification. Adopting our measurement of μ_best = 21.4, we find this results in a stellar mass of M_⋆ = 1.82 × 10⁸ M_⊙. We note that Molino et al. (2017) did not use Spitzer/IRAC data, so, in comparing to our stellar mass measurement of ${M}_{\star }\,={1.45}_{-1.19}^{+4.14}\times {10}^{8}{M}_{\odot }$ derived without IRAC, the two results are in agreement.

For MACS0744-064 (CLASHID = 1176), the catalog photometric redshift is $z={4.14}_{-3.70}^{+1.42}$ (95% conf.). This is significantly different than our photo-z, which does not show any probability near the peak of their distribution, that is, at z ∼ 4 (Figure 7). We investigated different origins for this discrepancy. As for RXJ1347-018, Molino et al. (2017) do not use the IRAC [3.6] and [4.5] measurements as constraints on the photo-z. However, this is not the primary reason for the photo-z discrepancy. We reran EAzY on our photometry without the two IRAC bands, and we found approximately the same photo-z as when we included the IRAC data. We also checked to see if the different photo-z's arose from differences in HST photometry. We reran EAzY using the Molino et al. (2017) photometry and found $z={6.39}_{-5.27}^{+0.51}$, similar to the result using our own photometry. This means that the difference most likely arises from the different codes used to derive the redshifts, rather than the photometry. Molino et al. (2017) use Bayesian photometric redshifts (BPZ; Benítez 2000; Coe et al. 2006), whereas we use EAzY. We note that the brightness of this galaxy is near the detection limit of the HST images; the S/N values in the four bands in which it is detected above S/N = 3 are 3.8 for F160W, 4.0 for F140W, 3.4 for F125W, and 5.3 for F110W. RXJ1347-018 is detected with much higher S/N in HST, and there the photo-z's are in much better agreement. As a result, we expect that the different templates and other slight implementation differences between the two codes are amplified when the S/N of the data is low, giving rise to the different photometric redshift result.

The Lyα properties of the two galaxies are also shown in Tables 3 and 4. We calculated rest-frame Lyα equivalent widths (EW_Lyα) by dividing the line flux by the continuum flux density measured in F160W and then by the scale factor (1 + z). MACS0744-064 has EW_Lyα = 59.4 Å, and RXJ1347-018 has EW_Lyα = 27.8 Å. These are within the range of values for EW_Lyα among the other known Lyα emitters at z = 7–7.5 (see Table 2 of Stark et al. 2017). However, our measurements are for much fainter galaxies than are typically observed at z ∼ 7, due to the lensing magnification. MACS0744-064 and RXJ1347-018 have intrinsic absolute magnitudes of ${M}_{\mathrm{UV}}\,-2.5{\mathrm{log}}_{10}(\mu /{\mu }_{\mathrm{best}})$ = −18.5 ± 0.4 mag ($L={0.12}_{-0.04}^{+0.06}{L}_{\star }$) and M_UV − 2.5 log₁₀(μ/μ_best) = −17.2 ± 0.2 mag ($L=0.03\,\pm 0.01{L}_{\star }$), respectively.¹⁵ Galaxies this faint far outnumber those observed in the field at luminosities of L ∼ L_⋆, which likely live in rare and massive halos (Roberts-Borsani et al. 2016; Mason et al. 2018a). As a result, we can probe the neutral hydrogen fraction from more typical halo masses at z ∼ 7–8, as opposed to fewer halos at the high-mass end. An ideal approach, and one we hope to adopt in future work, is to combine data sets to probe the neutral fraction using the entire available range of halo masses.

The remaining 20 of 22 candidate lines that were flagged by the automatic line detector are much less secure than the two lines presented above. Most fail multiple of the visual inspection tests. However, two of these correspond to real emission lines from a single lower-z galaxy that we did not target, but which happened to fall in the same slit as one of our z ≳ 7 targets, MACS1149-014. The two emission lines are confidently identified as the [O iii]λλ4959,5007 emission line pair at z = 1.225, and they are present at a spatial offset from the central trace of our actual target, consistent with that expected from the position of the slit relative to our target and the lower-z contaminating galaxy in the HST image.

We note that while Hoag et al. (2017) presented a detection of Lyα in ID = MACS1423-040, an object in our MOSFIRE sample, it does not pass the above automatic line detection criteria. This is primarily due to our rescaling (increasing) of the noise as described in Section 2.1. If we perform the same extraction of the line as done by Hoag et al. (2017) but use the rescaled noise, we find S/N_int = 3.4. Because the Lyα EW of this object is so low (9 Å; Hoag et al. 2017), it has an insignificant impact on our neutral fraction results. We demonstrate this in Section 3.7 by inferring the neutral fraction both with this target as a detection and as a nondetection.

While we likely only detect Lyα in two of the 68 galaxies in our sample, our nondetections are useful constraints on the neutral fraction. An individual 1σ flux limit spectrum, f_lim (λ), is obtained from the scaled noise spectrum of each object, σ(λ)_MOS,scaled, via

$$\begin{eqnarray}&&{f}_{\mathrm{lim}}{(\lambda )={\rm{\Delta }}\lambda \times \sqrt{2\ast \mathrm{FWHM}/{\rm{\Delta }}\lambda }\times \sigma (\lambda )}_{\mathrm{MOS},\mathrm{scaled}},\end{eqnarray} \tag{ 1 }$$

where Δλ = 1.086 Å is the spectral pixel size and FWHM is the FWHM of the MOSFIRE Y-band spectral resolution, ∼3 Å. We show the median flux limit spectrum of the sample in Figure 10. An individual 1σ rest-frame Lyα equivalent-width spectrum is obtained from the flux limit spectrum via

$$\begin{eqnarray}&&{\mathrm{EW}}_{\mathrm{Ly}\alpha }={f}_{\mathrm{lim}}(\lambda )/{f}_{\mathrm{cont}}/(1+z),\end{eqnarray} \tag{ 2 }$$

where f_cont is the continuum flux density measured in F160W. We note that F160W samples the rest-UV of galaxies in our redshift range without encompassing the Lyα line. The neutral fraction inference uses the entire flux density spectrum in the calculation of the likelihood for each object, which takes into account the variable sensitivity of the spectrum (Mason et al. 2019).

Zoom In Zoom Out Reset image size

Here we use the formalism introduced by Treu et al. (2012) and extended by Mason et al. (2018b) and Mason et al. (2019) including models by Mesinger et al. (2015) to infer the volume-averaged neutral hydrogen fraction in the IGM at z ∼ 7.5 using our Lyα spectroscopy from MOSFIRE. We adopt z = 7.6 as our fiducial redshift because it is the average Lyα redshift probed by the MOSFIRE Y band. The formalism uses the Evolution of 21 cm Structure cosmological simulations (Mesinger et al. 2015, 2016), that is, cosmological-scale IGM simulations of inhomogeneous reionization, to obtain the global neutral fraction distribution as well as a realistic treatment of the impact of the ISM and CGM on the emerging Lyα line.

The model yields a likelihood that can be used to interpret observables of LBGs in terms of neutral fraction. The input observables are the rest-frame Lyα equivalent-width measurements and limits (for nondetections), as well as galaxy luminosities. For nondetections, we use the entire flux density and noise spectra, taking into account the varying sensitivity as a function of wavelength (see Figure 10).

Our ignorance of the underlying Lyα FWHM distribution of the objects without Lyα detections introduces uncertainty into the inferred neutral fraction. Furthermore, larger assumed FWHMs result in weaker constraints, due to the line flux being spread out over more detector pixels. To parameterize our ignorance of the FWHM, we marginalize over it when calculating the posterior for each object. We explored three different priors for the FWHM distribution: (1) a uniform prior in the range 100–400 km s⁻¹ for each object; (2) a log-normal $p(\mathrm{FWHM}| {M}_{\mathrm{UV}})$ that uses the Mason et al. (2018b) scaling relation between Lyα velocity offset and M_UV and the Verhamme et al. (2018) scaling relation between velocity offset and FWHM, where we adopt the median M_UV of the sample for each object; (3) the same functional form and scaling relations as in (2), but where we use the individual M_UV for each galaxy rather than the median of the sample in the scaling relations. We found that the result was robust against the choice of prior, so we adopted case (3) as it is the most physically motivated. The log-normal tends to peak around 50–100 km s⁻¹ for the M_UV values in the MOSFIRE sample. For more details and the form of the full likelihood, see Mason et al. (2018b).

We also use the photometric redshift probability distributions P(z) as priors in the likelihood. Effectively, we weight each galaxy in the inference, where the effective weight is given by the probability that Lyα falls within the MOSFIRE Y band, that is, p(7 < z < 8.2). Doing so accounts for the photometric redshift impurity of our sample. Galaxy intrinsic luminosities enter into the inference via the halo-mass dependence (and hence luminosity dependence) of the local neutral fraction in the IGM simulations. Furthermore, the model assumes an intrinsic emitted EW distribution at z ∼ 6 before any attenuation. This model is based on real observations at z ∼ 6 compiled by De Barros et al. (2017). For more details, see Mason et al. (2018b) and Mason et al. (2019).

After applying the prior, we obtain the posterior, $p({\overline{x}}_{{\rm{H}}{\rm{I}}}| \{{\rm{f}}(\lambda ),m,\mu \}$), where the data consist of the set of flux density spectra, f(λ), apparent magnitudes, m, and magnifications, μ. For objects where there are detections, we use the measured equivalent width of the line rather than the entire flux spectrum to construct the posterior.

We show our posterior on the volume-averaged neutral fraction ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ in Figure 11. Using our fiducial sample of two Lyα detections (see Section 3.3), we infer a neutral fraction of ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.88}_{-0.10}^{+0.05}$. Our result is robust against the set of Lyα detections we choose, as long as the MACS0744-064 detection is included. Excluding the RXJ1347-018 detection (the one-detection case), we infer a statistically consistent neutral fraction to the fiducial case: ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.89}_{-0.10}^{+0.04}$. If instead we include the additional Lyα detection presented by Hoag et al. (2017; the three-detection case), we also infer a statistically consistent neutral fraction: ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.88}_{-0.09}^{+0.05}$. The difference between the three cases shown is small because in all cases the MACS0744-064 detection is included. This detection is S/N ∼ 10, so it has a much larger statistical weight in the inference than the other detections. While we find the zero-detection case unlikely, we cannot rule it out because there is a small chance that the single emission line in each of the three spectra is not Lyα. If this were the case, we would put a lower limit on the neutral fraction.

Zoom In Zoom Out Reset image size