Line Emitting Galaxies Beyond a Redshift of 7: An Improved Method for Estimating the Evolving Neutrality of the Intergalactic Medium
Matthew A. Schenker, Richard S. Ellis, Nick P. Konidaris, Daniel P. Stark
aa r X i v : . [ a s t r o - ph . C O ] A p r Draft version September 24, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
LINE EMITTING GALAXIES BEYOND A REDSHIFT OF 7: AN IMPROVED METHOD FOR ESTIMATINGTHE EVOLVING NEUTRALITY OF THE INTERGALACTIC MEDIUM
Matthew A Schenker , Richard S Ellis Nick P Konidaris , Daniel P Stark , Draft version September 24, 2018
ABSTRACTThe redshift-dependent fraction of color-selected galaxies revealing Lyman alpha emission, x Ly α hasbecome the most valuable constraint on the evolving neutrality of the early intergalactic medium.However, in addition to resonant scattering by neutral gas, the visibility of Lyman alpha is alsodependent on the intrinsic properties of the host galaxy, including its stellar population, dust contentand the nature of outflowing gas. Taking advantage of significant progress we have made in determiningthe line emitting properties of z ≃ − z > < z < z ≃
8. As a byproduct, we also present a new line emitting galaxy at a redshift z = 7 . < z < Subject headings: galaxies: evolution INTRODUCTIONThe transition from a neutral intergalactic medium(IGM) to one that is ionized, and therefore transpar-ent to ultraviolet photons, represents the latest frontierin our overall understanding of cosmic history. In addi-tion to determining when this ‘cosmic reionization’ oc-curred, a key question is the role of star-forming galax-ies in governing the process. Structure in the polariza-tion of the cosmic microwave background suggests thereionization process occurred within the redshift inter-val 6 < z <
20 (Hinshaw et al. 2013) and deep in-frared imaging with Hubble Space Telescope has pro-vided the first opportunity to conduct a census of galax-ies during the latter half of this period (Ellis et al. 2013;Oesch et al. 2013). Recent progress in this area has beenreviewed by Robertson et al. (2013) and Bromm (2013).In the absence of significant numbers of high redshiftQSOs or gamma ray bursts, the most immediately avail-able probe of the evolving neutrality of the IGM beyond z ≃ − α ) emis-sion line in controlled samples of color-selected galaxies.Although a prominent line in star-forming galaxies at z ≤
6, as Ly α is a resonant transition, it is readily sup-pressed by neutral gas, both in the host galaxy and, ifpresent, in the surrounding IGM. First proposed as apractical experiment using Lyman break galaxies (LBGs)by Stark et al. (2010), the idea followed earlier theoreti-cal work by Miralda-Escud´e et al. (2000), Santos (2004)and others.Ground-based near-infrared spectroscopic surveys have Department of Astrophysics, California Instituteof Technology, MC 249-17, Pasadena, CA 91125;[email protected] Department of Astronomy and Steward Observatory, Uni-versity of Arizona, Tucson AZ 85721 Hubble Fellow now targeted various samples of color-selected Lymanbreak galaxies over 6 < z < α fraction, x Ly α ,which falls sharply from a value of ≃
50% at z ≃ z ≃ HI , is uncertain(Bolton & Haehnelt 2013), the prospects for improvingthe statistics of this test are promising given the arrivalof multi-object instruments such as MOSFIRE on theKeck 1 telescope (McLean et al. 2012).So far, this important measure of late reionization hasbeen applied by adopting an empirical description of thedemographics of Ly α emission in LBGs, parameterizedaccording to the equivalent width (EW) distribution forvarious UV luminosities over the redshift range 3 < z < < z < α in the vicinity of the hostgalaxy. Our new approach aims to predict its visibilityin a high redshift galaxy on the basis of its measured UVcontinuum slope that, in turn, contains information onthe dust content, and stellar population, which both di-rectly influence the strength of any Ly α emission. Thisapproach has the distinct advantage that, for the new z > α fraction test and applying it to thefirst comprehensive set of spectroscopic data emergingfrom MOSFIRE. In addition to incorporating the ear-lier surveys conducted with Keck (Schenker et al. 2012;Ono et al. 2012; Treu et al. 2012), and FORS2 on theVLT (Pentericci et al. 2011), we present the first resultsfrom a survey of high quality Ultra Deep Field (UDF)targets which provides a valuable extension of the afore-mentioned studies. As part of this survey, we demon-strate a new Ly α -emitting galaxy at a redshift z =7.62extending once again the frontier of spectroscopically-confirmed HST sources.A plan of the paper follows. In Section 2, we introduceour new method for the Lyman alpha fraction test. Sec-tion 3 introduces the new compilation of 3 < z < LYMAN ALPHA FRACTION TEST - A NEWAPPROACHAlthough the traditional Ly α fraction test as first pro-posed by Stark et al. (2010) has already provided mean-ingful constraints on the evolution of the IGM beyond z ∼ .
5, there are two limitations in the current method-ology. Firstly, as inferring the presence of neutral gas inthe intergalactic medium represents a differential mea-surement, it is necessary to assume a form of the distri-bution of the equivalent widths of Ly α emission unpro-cessed by the IGM for the galaxies at z ≥
7. Comparingthis to the observed distribution allows the extinctionimposed by the IGM, and through the application of the-oretical models, the IGM neutrality to be derived. Thecurrent methodology splits the sample into UV-luminousand UV-faint bins, and tracks the Ly α fraction in eachbin as a function of redshift. The IGM unprocessed dis-tribution at z ∼ z ∼ < z < z ≥ α distribution has been characterized in manydifferent ways, including an exponential (Dijkstra et al.2011), a direct histogram (Schenker et al. 2012), a halfgaussian (Treu et al. 2012) and a half gaussian with aconstant probability tail (Pentericci et al. 2011). Thoughthe distributions are largely similar, no detailed compar-ison has been performed to determine which one opti-mally represents the 3 < z < rest-frame UV luminosity is the optimum parameterto predict the visibility of Ly α in the absence of any sup-pression by a partly neutral IGM. The approach, basedon correlations first noted by Shapley et al. (2003), wasadopted by Stark et al. (2010) as M UV can be readily de- termined from the available photometry of distant galax-ies together with a photometric redshift. However, M UV is likely to be a coarse predictor of the Ly α EW as itignores second order parameters such as metallicity, thestellar initial mass function and dust content.The
UV continuum slope is a more natural choice as abasic variable as it encodes each of these physical prop-erties (Meurer et al. 1999). Lower metallicity and hotterstars produce more ionizing photons per unit UV con-tinuum, thus driving EW Ly α upwards. Dust very effi-ciently absorbs Ly α photons given their large effectivepath lengths from the many scatterings required to es-cape an HII region. These changes also result in a blueror redder UV continuum slope, respectively. Thus, asthe UV slope reflects more of the parameters that likelygovern EW Ly α compared to M UV , we should expect itto be a more robust predictor of the visibility of the linein high redshift samples.Until recently, determining the UV continuum slopewas only possible for a restricted subset of z < B -dropouts. Stark et al. (2010) showed that within thissubset, strong Ly α emitting galaxies have bluer UV con-tinuum slopes than their non-emitting counterparts, butas there existed no high-quality infrared photometry inthe GOODS fields at this time, it was necessary toparametrize distributions at higher redshift by their ab-solute magnitude. However, in addition to the now com-pleted Keck spectroscopic survey of LBGs over 3 < z < Ly α . The ad-dition of Y , J , and H photometric data enablesthe derivation of accurate UV continuum slopes for cata-logued galaxies, given for each source there is a minimumof 3 broad-band filters longward of the Lyman break.As such UV continuum slopes are now available for thegrowing body of z > α fraction test that overcomes several of the issuesassociated with the earlier approach.In the following, we discuss the new data for the Keck3 < z < z ≃ IMPROVED POST-REIONIZATION DATA3.1.
DEIMOS/FORS2 Spectroscopy
As discussed in Stark et al. (2010, 2011), the 3 < z < B -, v -, and i -drops, andthe spectroscopically-confirmed sample spans a redshiftrange of 3 . < z < .
99. Typical 10 σ limiting Ly α fluxes for these targets ranged between 1.0-1.5 × − erg cm − s − .The complementary FORS2 campaign (Vanzella et al.2009) targeted 214 LBG candidates in GOODS-S be-tween 2002 and 2006. These targets were, on average,brighter than those studied at Keck (see Stark et al.2010, Figure 2), and the confirmed galaxies span a red-shift range 3 . < z < .
28. In total, the sample com-prises 607 galaxies, 269 of which are spectroscopicallyconfirmed. We direct the reader to Stark et al. (2014)for further details. 3.2.
Photometry
The primary advance in our analysis of the equiva-lent width distribution of Ly α in the above spectroscopicsample relates to the combination of the earlier HSTACS optical imaging data with new, deep WFC3/IRnear-infrared data critical to assessing how Ly α emis-sion correlates with the measured UV continuum slope.To reliably bring together the various imaging datasets,it is necessary to account for the significantly differentpoint-spread functions (PSFs) between the ACS (FWHM ∼ ∼ ColorPro program (Coe et al. 2006). APSF was constructed for each filter by shifting and stack-ing ∼
20 bright unsaturated stars. All objects were de-tected using the H image, and colors were determinedusing matched isophotal apertures after degrading theresolutions of all other images to that of the H image. ANALYSIS4.1.
Lyman alpha and the UV continuum
The basis of our analysis relies on accurate determina-tions of both the Ly α equivalent widths (EW Ly α ), andultraviolet slopes of our sample. Thus, we detail here themethodology used in determining both these quantitiesfor use in our analysis.In order to measure the UV slope for each object in oursample, we first used a custom photometric redshift codeto determine the approximate redshifts of those galax-ies without spectroscopic confirmation. The code fits asuite of synthetic fluxes from the Bruzual-Charlot (BC03,Bruzual & Charlot 2003) spectra to the available pho-tometry. To determine the best-fit redshift, we marginal-ized across all other parameters (mass, dust extinction, and age), and used the maximum likelihood value.We measure the UV continuum slope using the β for-malism first introduced by Calzetti (1994), where theflux is parametrized as f λ ∝ λ β . In our fitting pro-cess, we include photometric filters with central restwavelengths within the range defined by Calzetti (1994),1350 < λ/ ˚A < β = − − . < β < − . − . < β < . β = 0 .
01 was fit to the observed photometry, andthe relative likelihood of each computed using p ( β i ) ∝ exp( − χ / p ( β ) were manually inspected foreach galaxy, flagging and removing objects with clearlydeviant photometry or incorrect solutions. After thisprocess, 297 of our sample of 393 galaxies in GOODS-N, and 154 of our 214 galaxies in GOODS-S remained.Galaxies removed were typically those adjacent to brightobjects, for which accurate photometry could not be as-sured, and faint, distant z ≥ α EW were taken fromStark et al. (2014) with errors computed using both the1 σ flux errors in the spectrum, as well as errors in M UV ,added in quadrature to produce a likelihood curve for p (EW Ly α ). For cases where Ly α was not spectroscopi-cally detected, we assume one of three cases: (1) the lineflux falls below the 10 σ limit (typically 1.0-1.5 × . − erg cm − s − ), (2) the line emission, though brighter thanthe 10 σ limit is missed, due to poor sky subtraction orobscuration by skylines, or (3) the object is a contami-nant outside the expected redshift range. We discuss theimplementation of this approach in Section 4.3.1.With these results in hand, we now provide the basicevidence that the UV continuum slope of a galaxy, β , isa reliable predictor of its Ly α EW. In Figure 1, we plotthe best fit β for each galaxy in our final sample againsteither its EW Ly α , or its 10 σ upper-limit, if a line is notdetected. The red crosses denote the mean EW Ly α in binsspanning ∆ β = 0 .
25, where undetected objects are set tohave a value of 0. Error bars for each bin are calculatedby bootstrap sampling.As shown from earlier work by Shapley et al. (2003)and Stark et al. (2010), a clear trend of increasing meanEW Ly α with bluer (more negative) UV slopes is visible.With our large spectroscopic sample, this can be seen di-rectly through the measures of individual objects, ratherthan via stacked spectra or consideration of average β values. Encouraged by this trend, we now develop amodel that can predict the probability distribution ofEW Ly α given a measured value of β , for example for a z > Ly α than the absolute UV magnitude. 4.2. The UV slope-dependent EW distribution
We now seek to establish a formalism to predict theprobability distribution for EW Ly α given a particularmeasured UV slope. As our goal is to improve the accu-racy of the Lyman alpha fraction test, we first considerthose galaxies with blue UV slopes, likely to have strongLy α emission before processing through the IGM.4.2.1. Equivalent Width Distributions for a Fixed UVContinuum Slope
As an illustration of our new method we begin by ex-amining the distribution of the Ly α EW for a fixed UVcontinuum slope, β . A natural choice is β = − . UV ≥ − z ≥
6) galaxies have slopes that asymp-tote to this value. Near-IR spectrographs such as MOS-FIRE have begun to target these galaxies in earnest (thiswork, Finkelstein et al. 2013, Treu et al. 2013), and it isbecoming increasingly important to characterize their ex-pected IGM unprocessed Ly α emission. To construct asample of galaxies for this task, we limit our overall sam-ple of 468 galaxies with 3 . < z < .
28 to those galaxieswith a best fit value within ∆ β = 0 .
25 of − .
3, resultingin a total of 131 objects.We now require a model to represent the IGM unpro-cessed EW distribution. In Appendix A, we review foursuch options using the methodology outlined here, andfind that a lognormal distribution provides a significantlybetter fit than any of the others. In this case the nat-ural logarithm of EW Ly α obeys a normal distribution.The two relevant parameters of the distribution, µ and σ , are typically referred to as the location parameter and scale parameter , respectively. They denote the mean ofthe natural log of EW Ly α and its variance. However,while the median of the distribution is given, as mightbe expected, by exp( µ ), the mean is slightly larger atexp( µ + σ / A em , determines thefraction of galaxies that have EW Ly α >
0, as there is noreason a priori to expect all galaxies to display Ly α inemission. The resulting distribution can be written as: p ( EW ) = A em × √ πσ EW exp( − (ln(EW) − µ ) σ )+ (1 . − A em ) × δ ( EW ) (1)In Figure 2, we illustrate how the EW Ly α probabil-ity distribution function, p (EW), and the Ly α fraction,x Ly α , change as these parameters are varied.We now describe the Bayesian formalism we developedto evaluate the likelihood of the underlying parametersfor our lognormal distribution, and determine which pro-vides the best fit to the data. The entire set of spectro-scopic Ly α observations is denoted as Obs ; this con-tains the information for observations of each individualgalaxy, Obs i . We can then denote the parameters for themodel being fit as θ ≡ [ µ, σ, A em ]. Our overall goal is todetermine the probability distributions for the underly-ing parameters of each model, given our observations, i.e. p ( θ | Obs ). Using Bayes’ theorem, we can rewrite this as p ( θ | Obs ) ∝ p ( θ ) × p ( Obs | θ ) (2)Here, p ( θ ) represents our uniform priors for the under- lying parameters, while the term on the right hand siderepresents the probability of our observations given themodel parameters. For any single object in which wemeasure a definite EW Ly α , this posterior probability canbe expressed as: p ( Obs i | θ ) = Z ∞ p ( EW Obs,i ) p ( EW | θ ) dEW (3)In the case of an object for which Ly α remains unde-tected above our (10 σ ) limit, we compute the posteriorprobability as: p ( Obs i | θ ) = p ( EW < EW σ | θ )+ p ( EW > EW σ | θ ) × C + C (4)Here, the first term represents the probability that theobject intrinsically posses an EW Ly α below our detectionlimits, while the C term takes into account incomplete-ness in the sample (caused by skylines or, in some cases,poor background subtraction). Contamination by lowredshift sources is taken into account through the finalterm, C . We assume modest values for our contamina-tion terms of C = C = 0.05, motivated by the com-pleteness simulations of Stark et al. (2010), and the lackof low-redshift interlopers found in other spectroscopicfollow-up surveys (Pentericci et al. 2011). The full pos-terior distribution for the parameters can then be com-puted by multiplying the individual posterior probabil-ities for each object. This allows us to infer the mostlikely values on parameters, as well as their marginalizedand un-marginalized errors.We display the best-fit distribution for our sample of131 galaxies with 3 . < z < .
28, overplotted on a his-togram of their EW Ly α detections and upper limits inFigure 3. The best fit parameters are µ = 2 . ± . σ = 1 . +0 . − . , and A em = 1 . +0 . − . . As we show in Ap-pendix A, this formalism is the best fit to our post-reionization data.4.2.2. A Generalized Approach
Now that we have determined the distribution whichbest fits the data at the key UV slope value of β ∼ − . β . Although faintgalaxies at z ≥ β ∼ − .
3, the UV-bright galaxies almost certainly do not (Bouwens et al.2013). In order to fully leverage the x Ly α test for themore luminous objects, we must use a model that deter-mines the EW distribution across a wide range of β .To achieve this goal, we extend the Bayesian formalismintroduced above. The differences are twofold. Firstly,we now include our entire sample of spectroscopically ob-served galaxies when fitting, rather than just those nar-rowly clustered around a particular value of β . Secondly,we must reconsider the nature of Eqn. 1 since it is clearthat the EW distribution varies as a function of β fromFigure 1.It is most reasonable to consider that the location pa-rameter , µ , varies with β since it is this parameter thatgoverns the redistribution of EWs (see Figure 2). Forconvenience we assume that µ varies linearly with β ,whence µ ( β ) = µ a + µ s × ( β + 2), where µ a represents the Reanalysis of the Lyman Alpha Fraction 5location parameter at β = − .
0. Prior to selecting thismodel, we performed similar fits to a narrow range in β as in the above section, with different ranges in β , andfound that while the best fit for µ varied strongly withthe UV slope, both σ and A em did not.Figure 4 shows a representation of the product of theresulting formalism, while we detail the full results inAppendix B. Because our model treats the location pa-rameter as a linear function of β , we can generate an EWprobability distribution function for any UV slope, al-though our sample only provides meaningful constraintsin the observed range − . < β < − .
25. In this fig-ure, we plot 3 examples and their associated errors. Ofcourse, a measured UV slope for a high redshift galaxyhas some error uncertainty, and thus its own probabilitydistribution, p ( β ). To obtain p ( EW ) given our model,and an observation of β , we can marginalize over β foreach galaxy in the sample: p ( EW i ) = Z p ( EW | β ) p ( β i ) dβ (5)4.3. UV Slope versus UV Luminosity
We now return to one of the assumptions motivatingthis paper. Does parametrizing the likelihood of Ly α emission via the UV slope represent a statistically betteroption than the combination of absolute UV magnitudeand redshift used in previous high-redshift Ly α studies?Since we can measure β , M UV , EW Ly α and a photo-metric or spectroscopic redshift for 451 objects in ourspectroscopic sample, we can directly address this ques-tion.In the widely-used UV luminosity method, the depen-dence on M UV is handled in a discrete manner, assign-ing galaxies into one of two bins, depending on whetherthey are greater or less than M UV = − .
25. x Ly α is also calculated in discrete bins at z = 4 , , and 6.For purposes of comparison, we use the Bayesian for-malism developed in this paper to generalize this to acontinuous model. To do so, we alter the definition of µ from µ ( β ) = µ a + µ s × ( β + 2 .
0) to µ ( M uv , z ) = µ a + µ s,M uv × ( M uv + 19 .
5) + µ s,z × ( z − . µ , on UV slope, with a bivariate linear dependenceon UV magnitude and redshift. We then use the samefitting process to determine the optimal values for allparameters in the model.To compare how well each of these two models fits theavailable data, we use the Bayesian evidence ratio, orBayes factor. The evidence is a measure of how likelythe data are, given the model, and can be expressed asan integration of the likelihood function over all possiblevalues for each parameter in the model E = Z p ( θ | Obs ) dθ (6)Evaluating this for both models yields a significant gainvia a ratio of E β /E M uv ,z = 29, showing convincing pref-erence for the model based on β compared to the earliermethod. FIRST APPLICATION TO DATA WITHIN THEREIONIZATION ERA Although the body of spectroscopic data targetinggalaxies beyond z ≃ A New MOSFIRE Survey
As part of a long term survey targeting z >
The GOODS-South / Ultra Deep Field
On the night of Nov 5, 2013, we secured a total of 3.5hours exposure in the Hubble Ultra Deep Field regionof GOODS-South. Observations were taken using 0.7”slits through intermittent high cirrus cloud and ∼ . z > z > . Y -drops outside the UDF properfrom Oesch et al. (2012). We used the photometric red-shift code, described in Section 4, to compute a redshiftprobability distribution for each object.The UDF 2012 dataset(Koekemoer et al. 2013) (GO12498, PI; Ellis) offers many distinct advantages forthis first application of our method. Foremost, byvirtue of the extraordinarily deep optical and F105Wdata available, contamination by foreground objects, asdetermined by the photometric scatter simulations inSchenker et al. (2013a) is ∼ J = 28 . ≥
10% con-tamination rate affecting galaxies in the CANDELS fields(Oesch et al. 2012). Secondly, as a result of a strategicdeployment of near-infrared filters, our UDF 2012 can-didates have better defined redshift probability distribu-tions, allowing us to more confidently exclude the pos-sibility of Ly α emission in the event of a non-detection.The median 68% confidence interval in photometric red-shift for the UDF 2012 objects on our MOSFIRE mask issmaller by ∆ z = 0 . § Ly α . As a first attempt, we used the z ∼ Ly α for each object, as a function of its UV magni-tude. The fractional number of expected detections wascalculated for each object, taking into account the pho-tometric redshift distribution (as our spectral coverage isincomplete), UV magnitude, and expected limiting fluxfor a likely MOSFIRE exposure. This exercise resultedin the final list of 16 candidates presented in Table 1. 5.1.2. CLASH Lensing Sample
Over the course of our November and March observa-tions, we also targeted 3 candidates with a photomet-ric redshift z ≥ . λ = 7750 ˚A theyhave sharp redshift probability distributions, and well-determined UV slopes making them ideally suited forour new method.We first targeted the z ∼ . Data reduction
The data was reduced using the publicly availableMOSFIRE data reduction pipeline . This pipeline firstcreates a median, cosmic-ray subtracted flat field imagefor each mask. Wavelength solutions for each slit arefit interactively for the central pixel in each slit, thenpropagated outwards to the slit edges to derive a fullwavelength grid for each slit. The sky background is es-timated as a function of wavelength and time using aseries of B-splines and subtracted from each frame. Thenodded A-B frames are differenced, stacked, rectified andoutput for use along with inverse variance-weighted im-ages used for error estimation.Examination of our GOODS-S data revealed a grad-ual 0 .
6” ( ∼ J ∼
19 star con-veniently placed on one slit. The intensity of the starallowed us to follow the extinction for each frame andeliminating those frames affected by significant extinc-tion or drift, we secured 2.35 hours of useful exposure.Due to this drift, the star on our original reduction dis-played a somewhat greater FWHM of ∼ .
2” than any ofour individual exposures, which typically had a FWHM ∼ . ∼ ∼ . σ limiting sensitivity between skylinesof ∼ . × − erg cm − s − . We note that approxi-mately ∼
33% of the Y-band spectral range is obscured https://code.google.com/p/mosfire/ by skylines at the MOSFIRE resolution of R ∼ A New z=7.62 Lyman Alpha Emitting Galaxy
We inspected the reduced, two-dimensional spectra ofall 16 objects by eye to search for Ly α emission. Fromthis, we were able to locate only a single candidate emis-sion line. Surprisingly, this emission line is located inone of our faintest targets, UDF12-3313-6545 (first iden-tified by McLure et al. 2010, Bouwens et al. 2011), witha measured flux of 5 . × − erg cm − s − . If the lineis indeed Ly α , the galaxy lies at a spectroscopic redshiftof z = 7 .
62, making it a promising candidate for themost distant spectroscopically-confirmed galaxy to date.We present the relevant details and HST cutouts of thisgalaxy in Figure 6.Given the faintness (the emission line is detected at4 . σ ), line asymmetry, commonly used to distinguishLy α , is not detectable. However, the fact that the linedisplays a clear positive signal flanked by two negativepeaks indicates that the signal was present in both the Aand B exposure positions. Although the line has a sur-prisingly large rest-frame equivalent width of 160 ± z ∼ α emitting galaxies (Ouchi et al. 2010). Notably, the spec-troscopic redshift lies well within to the 1 σ confidenceinterval of our derived photometric redshift distributionwhen line emission is accounted for, instilling further con-fidence in the redshift.5.3. Additional Data from Published Surveys
In order to achieve the most up-to-date and precisemeasurement of the Ly α opacity at z ≥ .
5, we havecompiled a comprehensive sample of other near-infraredsurveys for Ly α at high redshift, which we will uti-lize in our analysis. This includes our own prior workwith Keck’s NIRSPEC (Schenker et al. 2012), as well asa number of other surveys, using red-sensitive opticaldetectors on the VLT and Keck (Fontana et al. 2010;Pentericci et al. 2011; Ono et al. 2012; Pentericci et al.2014), as well as independent work by Treu and collabo-rators using NIRSPEC (Treu et al. 2012) and MOSFIRE(Treu et al. 2013).In total, this literature sample comprises 83 z ≥ . z ∼ z ∼
8. Themanner in which targets were assigned to each bin re-quired careful consideration given the limited wavelengthresponse of each instrument with respect to the photo-metric redshift likelihood distribution P ( z ). Rather thanbinning on photometric redshift alone, we carefully con-sidered the redshift range within which a null detectioncould be determined. If the median redshift for which anull detection could be determined was less than (greaterthan) z = 7 .
5, we place the object in the z ∼ Monte Carlo simulation
To predict the number of detections expected in anIGM with no additional opacity to Ly α , we use a similarMonte Carlo method to that developed in Schenker et al.(2012). This simulation has three key inputs for each Reanalysis of the Lyman Alpha Fraction 7object: flux limits as a function of wavelength from thespectroscopic observations, which also take into accountthe night sky emission, a photometric redshift probabilitydistribution, and a prediction for the IGM unprocessedEW Ly α (and thus f Lyα ) distribution.For objects observed in either this paper, orSchenker et al. (2012), flux limits as a function of wave-length were calculated directly from the reduced spec-tra by computing the variance in an aperture of 10 ˚Aspectral extent. For the data in Treu et al. (2012) andTreu et al. (2013) we rescaled our flux limits from NIR-SPEC and MOSFIRE, respectively, to match the quotedlimits in the paper for each object. For Pentericci et al.(2011); Ono et al. (2012); Pentericci et al. (2014), wedid the same with our LRIS limits, as presented inSchenker et al. (2012).We used the published photometry from each paper(and our own here) in conjunction with our photometricredshift code described previously to determine a photo-metric redshift distribution for each object. The only ex-ceptions are for the samples from Pentericci et al. (2011)and Ono et al. (2012), for which either photometry or co-ordinates were not available. For these objects, we usedthe photometric redshift distribution for ground-based z -drops from Ouchi et al. (2009).Finally, for the objects in our new MOSFIRE survey,we generated the prediction for the IGM unprocessedEW Ly α distribution using the observed UV slope, as de-scribed in Section 4. Ideally, we would prefer to usethis new method for all objects in the combined sam-ple, in order to eliminate the potential bias of simplyusing M UV as a predictor. However, with the exceptionof galaxies in the UDF 2012 field and CLASH lensedsample, the requisite 3 infrared photometric data pointslongward of the Lyman break essential for achieving anaccurate measure of β are not available. Thus, for allother objects, we must predict EW Ly α as a function ofM UV from Treu et al. (2012), using the data presentedin Stark et al. (2011). As an illustrative exercise, we alsogenerated a prediction for the IGM unprocessed EW Ly α distribution of these objects using their M UV to calcu-late a β derived from the M UV - β relation at z ∼ α equivalent width from the EW Ly α distribution. Fromthe observed photometry, this EW Ly α is converted to aflux. We then sample the spectroscopic flux limit at theredshift drawn to determine if the emission line wouldbe observed at ≥ σ . This process is then repeated with N = 10000 trials for each object.5.3.2. Comparison between UV slope and UV luminositypredictions
Before considering the total sample (i.e. including pre-vious data from the literature), we compare the differencein the expected Ly α statistics for the high redshift sam-ple using either M UV or β as the basis for predicting theIGM unprocessed EW Ly α distribution. Since only theUDF 2012 and CLASH surveys currently have accurateindividual measurements of β , this comparison can onlybe done for the 13 targets from our recent MOSFIRE survey.The number of expected detections is compared in Fig-ure 7. Although hampered by limited statistics, the dif-ference is still significant. Using β as a basis, we predictan average of 1.4 more detections than using M UV . Thisdifference represents an important correction of a system-atic error in the prior x Ly α tests. Our new results showthat the M UV method, for this specific sample of 13 tar-gets, significantly underpredicts the incidence of IGM un-processed Ly α emission, which results in an overestimateof the IGM transmission. The difference in predicted de-tections is dependent upon the properties of the sampleconsidered but, as the objects probed from the both theUDF and CLASH are intrinsically faint, with blue UVslopes, it is not surprising that the difference is so great.This change in predicted Ly α emission has impor-tant consequences for the transmission fraction of theIGM implying a lower limit on the neutral fraction thatis a factor 0.16 larger. Clearly, for a given survey,our new β method for predicting the IGM unprocessedEW Ly α distribution has significant advantages, and re-duces a key systematic error. This will be especiallyrelevant for spectroscopic follow-up of the HST FrontierFields (GO:13496, PI: Lotz), and their parallels, giventhese fields will have full coverage with the same fourWFC3/IR filters used for the UDF.5.3.3. Analysing the Entire Sample
We can now combine the various subsamples, our ownMOSFIRE survey and that from diverse sources in theliterature. The net result is a histogram of the numberof 5 σ detections overall. These histograms are displayed,both for the z ∼ z ∼ z ∼ z = 7 . σ flux limit.As an illustration, we can convert these observational-based results to an IGM extinction of Ly α by adoptingthe model used in Schenker et al. (2012) appropriate forpatchy reionization. In this model, the IGM is partiallyopaque, such that Ly α escapes to the observer unattenu-ated from a fraction, f , of galaxies, while it is completelyextinguished by the IGM in a fraction 1 − f . Given thehistogram of expected detections and the number actu-ally observed, we infer a probability distribution for thistransmission fraction, f . At z ∼
7, we find f = 0 . ± . z ∼
8, a 1 σ upper limit of f < .
19. The full re-sults can be found in Table 2 and are plotted along withthe lower redshift data on x Ly α in Figure 9.Discussing the uncertainties in the transformation fromtransmission fraction to x HI is beyond the scope of thepresent paper. However, clearly this conversion is de-pendent upon a number of physical parameters, someinternal to the galaxy, and others from the IGM stateitself. These include the velocity offset of Ly α fromthe galaxy’s systemic velocity (e.g. Hashimoto et al.2013; Schenker, et al. 2013b), the ionizing photon es-cape fraction (Dijkstra et al. 2014), as well as the pos-sible presence of optically thick absorption systems(Bolton & Haehnelt 2013). Until the theoretical mod-els converge, and/or observations of these key quantitiesare available, absolute measures of the neutral fractionwill still be subject to systematic errors. Nonetheless,we have demonstrated substantial observational progresswith our new survey and improved methodology, reduc-ing one of the key systematic errors. Using the modelsof McQuinn et al. (2007) here to provide an estimate ofx HI , we find x HI = 0 . +0 . − . at z ∼
7, and x HI > .
65 at z ∼ CONCLUSIONSUsing our sample of 451 3 < z < β , and its Ly α emission strength. Giventhe availability of deep WFC3 photometry for both theGOODS-N and S fields, this progress follows measure-ments for many individual galaxies in this redshift range,rather than via stacked or averaged UV slopes, as in ear-lier work (Shapley et al. 2003; Stark et al. 2010).We demonstrate that this correlation with the presenceof Ly α is stronger and more physically-motivated thanthat based on the UV luminosity and thus provides anatural basis for an improved model for the Ly α fraction test, now widely used to measure the evolving neutral-ity of the z > . < z < α emission at z ≃ z ≃ . . σ confirmationof Ly α in a galaxy at z = 7 .
62, likely the most distantspectroscopically-confirmed galaxy.We thank Chuck Steidel and Ian McLain for their ster-ling efforts in developing the highly successful MOSFIREinstrument. We also wish to recognize and acknowledgethe very signicant cultural role and reverence that thesummit of Mauna Kea has always had within the indige-nous Hawaiian community. We are most fortunate tohave the opportunity to conduct observations from thismountain.
REFERENCESBertin, E., Mellier, Y., Radovich, M., Missonnier, G., Didelon, P.,& Morin, B. 2002, Astronomical Data Analysis Software andSystems XI, 281, 228Bolton, J. S., & Haehnelt, M. G. 2013, MNRAS, 429, 1695Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2011, ApJ,737, 90Bouwens, R. J., et al. 2013, arXiv:1306.2950Bradley, L. D., et al. 2013, arXiv:1308.1692Bromm, V. 2013, Reports on Progress in Physics, 76, 112901Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ,429, 582Coe, D., Ben´ıtez, N., S´anchez, S. F., Jee, M., Bouwens, R., &Ford, H. 2006, AJ, 132, 926Dijkstra, M., Mesinger, A., & Wyithe, J. S. B. 2011, MNRAS,414, 2139Dijkstra, M., Wyithe, S., Haiman, Z., Mesinger, A., & Pentericci,L. 2014, arXiv:1401.7676Dunlop, J. S., et al. 2013, MNRAS, 432, 3520Ellis, R. S., et al. 2013, ApJ, 763, L7Finkelstein, S. L., et al. 2013, Nature, 502, 524Fontana, A., et al. 2010, ApJ, 725, L205Guo, Y., et al. 2013, ApJS, 207, 24Hashimoto, T., Ouchi, M., Shimasaku, K., Ono, Y., Nakajima,K., Rauch, M., Lee, J., & Okamura, S. 2013, ApJ, 765, 70Hinshaw, G., et al. 2013, ApJS, 208, 19Koekemoer, A. M., et al. 2013, ApJS, 209, 3McLean, I. S., et al. 2012, Proc. SPIE, 8446McLure, R. J., et al. 2010, MNRAS, 403, 960McLure, R. J., et al. 2013, MNRAS, 432, 2696McQuinn, M., Hernquist, L., Zaldarriaga, M., & Dutta, S. 2007,MNRAS, 381, 75 Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64Miralda-Escud´e, J., Haehnelt, M., & Rees, M. J. 2000, ApJ, 530, 1Oesch, P. A., et al. 2012, ApJ, 759, 135Oesch, P. A., et al. 2013, arXiv:1309.2280Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713Ono, Y., et al. 2012, ApJ, 744, 83Ouchi, M., et al. 2009, ApJ, 706, 1136Ouchi, M., et al. 2010, ApJ, 723, 869Pentericci, L., et al. 2011, ApJ, 743, 132Pentericci, L., et al. 2014, arXiv:1403.5466Robertson, B. E., et al. 2013, ApJ, 768, 71Rogers, A. B., et al. 2013, arXiv:1312.4975Santos, M. R. 2004, MNRAS, 349, 1137Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L.2003, ApJ, 588, 65Schenker, M. A., Stark, D. P., Ellis, R. S., Robertson, B. E.,Dunlop, J. S., McLure, R. J., Kneib, J.-P., & Richard, J. 2012,ApJ, 744, 179Schenker, M. A., et al. 2013, ApJ, 768, 196Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P.2013, ApJ, 777, 67Stark, D. P., Ellis, R. S., Bunker, A., Bundy, K., Targett, T.,Benson, A., & Lacy, M. 2009, ApJ, 697, 1493Stark, D. P., Ellis, R. S., Chiu, K., Ouchi, M., & Bunker, A. 2010,MNRAS, 408, 1628Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, ApJ, 728, L2Stark, D. P., Ellis, R., et al. 2014, in prep.Treu, T., Trenti, M., Stiavelli, M., Auger, M. W., & Bradley,L. D. 2012, ApJ, 747, 27Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli,M. 2013, ApJ, 775, L29Vanzella, E., et al. 2009, ApJ, 695, 1163
APPENDIX
A. MODELS FOR P(EW— β ) The maximum likelihoods inferred from each of the four distributions are noted in Table 3. These results demonstratethat the lognormal distribution provides the best fit to the available data - its likelihood surpasses that of any othermodel by two orders of magnitude. Thus, we use this distribution as the basis for the more general form of p ( EW | β )we consider next. B. RESULTS OF FULL MODELING PROCEDURE
For reference, and such that they are available for use in future work, we list here the final values for our generalizedlognormal fit to the EW Ly α distribution at 3 < z <
6. They are µ a = 3 . +0 . − . , µ s = − . ± . σ = 1 . ± .
1, Reanalysis of the Lyman Alpha Fraction 9 -2.5 -2.0 -1.5 -1.0UV slope ( β )050100150200250300 E W L y α Ly α detectionLy α upper limitBinned mean EW Ly α Fig. 1.—
Compilation of our entire GOODS catalog of Ly α equivalent widths as a function of UV slope, β . Red triangles show theaverage EW, binned in steps of 0.25 in β , displaying a strong increase toward bluer slopes. This dataset forms the basis of our predictivemodel for Ly α emission incidence as a function of UV slope. Ly α p ( E W L y α ) µ = 3.0 µ = 2.0 µ = 4.0 0 20 40 60 80 100EW Ly α p ( E W L y α > E W ) Ly α p ( E W L y α ) σ = 1.0 σ = 0.5 σ = 1.5 0 20 40 60 80 100EW Ly α p ( E W L y α > E W ) Ly α p ( E W L y α ) A em = 1.0A em = 0.8A em = 0.6 0 20 40 60 80 100EW Ly α p ( E W L y α > E W ) Fig. 2.—
Example curves for how our lognormal model of EW Ly α distribution varies with each parameter. Left: Probability distributionsfor EW Ly α . In each panel, the black curve has the same parameter values: µ = 3 . σ = 1 . A em = 1 .
0. From top to bottom, the twocolored curves each display the affect of a change in a single parameter on the distribution. Right: Complementary cumulative distributionfunctions for the same parameters used in each left panel. This method of display is especially useful, as the Lyman alpha fraction, x Ly α ,for any EW Ly α can simply be read off the plot by finding the value of the curve at the desired EW Ly α along the x-axis. TABLE 1
Summary of MOSFIRE survey for Ly α Name RA Dec J z phot t exp [hr] 5 σ EW limitA611-0193 8:01:00.32 36:04:24.3 26.0 7.9 2.6 12.9MACS0647-1411 6:47:40.91 70:14:41.6 26.1 7.6 1.0 20.7RXJ1347-0943 13:47:33.90 -11:45:09.4 26.4 7.5 1.0 27.7CANDY-2272447364 3:32:27.24 -27:47:36.4 27.1 9.1 2.35 60CANDY-2243349150 3:32:24.33 -27:49:15.0 27.2 8.5 2.35 57GS-2098-8535 3:32:20.98 -27:48:53.5 27.2 9.7 2.35 64GS-2533-8541 3:32:25.33 -27:48:54.1 26.4 7.6 2.35 30GS-2779-5141 3:32:27.79 -27:45:14.1 27.2 9.6 2.35 64GS-4402-7273 3:32:44.02 -27:47:27.3 26.9 7.4 2.35 50UDF12-3313-6545 3:32:33.13 -27:46:54.5 28.5 7.4 2.35 200UDF12-3722-8061 3:32:37.22 -27:48:06.1 27.8 7.7 2.35 110UDF12-3846-7326 3:32:38.46 -27:47:32.6 29.0 7.0 2.35 330UDF12-3880-7072 3:32:38.80 -27:47:07.2 26.9 7.5 2.35 220UDF12-3939-7040 3:32:39.39 -27:47:04.0 28.6 7.7 2.35 180UDF12-3762-6011 3:32:37.62 -27:46:01.1 28.4 8.1 2.35 150UDF12-3813-5540 3:32:38.13 -27:45:54.0 28.2 8.2 2.35 70UDF12-4256-7314 3:32:42.56 -27:47:31.4 27.3 7.1 2.35 50UDF12-4308-6277 3:32:43.08 -27:46:27.7 28.5 8.0 2.35 200UDF12-4470-6443 3:32:44.70 -27:46:44.3 27.3 7.7 2.35 70
Note . — MOSFIRE survey targets. 5 σ limiting sensitivities are calculated using themedian limiting flux limit between skylines, and assuming a spectroscopic redshift of z =7 . TABLE 2
Summary of MOSFIRE survey for Ly α Survey Observed 5 σ Detections Transmission fraction ( f ) x HI This work MOSFIRE M UV
19 0 < . > . β
19 0 < . > . z ∼ . +0 . − . . +0 . − . Composite z ∼ < . > . z ∼ a
72 11 0 . +0 . − . . +0 . − . Composite z ∼ a
27 0 < . > . Note . — Monte Carlo results for transmission fraction, f , and x HI . a For these results, where individual UV slopes are not available, we instead use individual valuesof M UV to predict a value of β , which in turn is used to generate the IGM unprocessed EW Ly α distribution. TABLE 3 Ly α functional forms Distribution name Equation Free parameters Reference Log max likelihoodLognormal A em √ πσ EW exp ( − (ln(EW) − µ ) / σ ) A em , µ , σ This work -75.4Half-gaussian √ πσ exp ( − (EW − µ ) / σ ) A em , σ Treu et al. (2012) -80.7’ w/ high-EW tail A em,g √ πσ exp ( − (EW − µ ) / σ ) + A em,c A em,g , A em,c , µ , σ Pentericci et al. (2011) - 80.8Declining exponential A em exp ( − EW/EW ) A em , EW Dijkstra et al. (2011) -77.9
Note . — List of the mathematical distributions used to fit the EW Ly α distribution at β ∼ − .
3, and the calculated maximum likelihoods.Having a maximum likelihood two orders of magnitude greater than any other distribution considered demonstrates the lognormal distributionprovides the best fit. a This is a note and A em = 1 . +0 . − . . We also provide a plot of the posterior probability distribution below so the reader is able toappreciate the sometimes non-negligible covariances between parameters. Reanalysis of the Lyman Alpha Fraction 11 β = -2.30510152025 ga l β = -2.30510152025 ga l Fig. 3.—
Histogram of our observed EW Ly α detections (solid blue) and upper limits (unfilled black) for galaxies with best fit β = − . ± .
25. Overplotted in solid black is our best fit lognormal model as described in Section 4.2.1. Ly α p ( E W L y α > E W ) β = -2.5 β = -2.0 β = -1.5 Ly α p ( E W L y α > E W ) Fig. 4.—
Inverse cumulative distribution functions of EW Ly α for our best fit model, plotted as a function of β . For a desired β , theLyman alpha fraction, x Ly α , for an arbitrary equivalent width is defined by the y-value of the curve at the given EW Ly α . λ [A]-1•10 -18 -18 -18 -18 f λ [ e r g / c m / s / A ] λ [A]-1•10 -18 -18 -18 -18 f λ [ e r g / c m / s / A ] Fig. 5.—
MOSFIRE observations of our lone target with visible line emission, UDF12 40242. The full two dimensional spectrum isshown at top, with the one dimensional spectrum, plotted at the bottom, along with the error spectrum in grey. Given our A − B reductionstrategy (described in Section 2), our 2D spectrum shows the expected positive signal (red line) flanked by two negative signals (blue lines),each separated by the amplitude of the dither pattern.
ACS
ACS p ( z ) λ [ µ m]27.028.029.030.027.028.029.030.0 m AB Fig. 6.—
Right: 1.5 arcsecond diameter cutouts and total magnitudes (see data section) of our MOSFIRE target, UDF12 40242. Asexpected for a z > z = 7 .
62 denoted in red. The solid black curve displaysthe pdf from the raw photometry, while the dashed grey curve shows the pdf after the observed MOSFIRE line flux has been subtractedfrom the Y data point. Left, Bottom: Best fitting SED for the galaxy, along with HST WFC3 Photometry. Reanalysis of the Lyman Alpha Fraction 13 det p ( n de t ) det p ( n de t ) det p ( n de t ) det p ( n de t ) Fig. 7.—
Left: Predicted number of detections for our MOSFIRE survey using the EW Ly α probability distribution from Treu et al.(2012), which uses M UV as the predictor. The observed number of 0 detections is indicated by the red crosshatch. Right: Predictednumber of detections for same survey, but using β as the predictor for EW Ly α , as outlined in Section 4. In this case, the average numberof expected detections is increased by a factor of 0.4, highlighting the importance of using a model that accurately predicts the IGMunprocessed equivalent width distribution. The equivalent upper limit on the transmission fraction is also decreased by a factor of 0.23. det p ( n de t ) det p ( n de t ) p ( f ) p ( f ) det p ( n de t ) det p ( n de t ) p ( f ) p ( f ) Fig. 8.—
Results from our new MOSFIRE campaign, combined with data from the literature. Top left: Histogram of expected 5 σ detections of Ly α , computed using the Monte Carlo method described in the text for our z ∼ α , assuming a patchy opacity. Dark blue and light blue shading encompass the 1 σ and 2 σ confidenceintervals, respectively. Bottom left and right: Same as above, but for our z ∼ x L y α , M UV > -20.25M UV < -20.25 Stark et al. (2011)This work Fig. 9.—
The fraction of Lyman break galaxies that display Ly α in emission at an EW ≥
25 ˚A, plotted as a function of redshift. Thevalues at z = 7 and 8 reflect differential measurements with the data at z = 6, as described in the text. Thus, these data points and errorsare simply the convolution of the x Ly α PDF at z = 6 and the transmission fraction PDF at z = 7 and 8. µ a A e m µ a A e m µ a σ µ a σ µ a -2.0-1.5-1.0-0.5 µ s µ a -2.0-1.5-1.0-0.5 µ s µ a µ s A e m -2.0 -1.5 -1.0 -0.5 µ s A e m -2.0 -1.5 -1.0 -0.5 µ s σ -2.0 -1.5 -1.0 -0.5 µ s σ -2.0 -1.5 -1.0 -0.5 µ s σ A e m σ A e m σ em Fig. 10.—
Posterior probability distribution for our full model, p (EW Ly α — ββ