The Keck Lyman Continuum Spectroscopic Survey (KLCS): The Emergent Ionizing Spectrum of Galaxies at $z\sim3$

We present results of a deep spectroscopic survey designed to quantify the statistics of the escape of ionizing photons from star-forming galaxies at z~3. We measure the ratio of ionizing to non-ionizing UV flux density_obs, where f900 is the mean flux density evaluated over the range [880,910] A. We quantify the emergent ratio of ionizing to non-ionizing UV flux density by analyzing high-S/N composite spectra formed from sub-samples with common observed properties and numbers sufficient to reduce the statistical uncertainty in the modeled IGM+CGM correction to obtain precise values of_out, including a full-sample average_out=$0.057\pm0.006$. We further show that_out increases monotonically with Ly$\alpha$ rest equivalent width, inducing an inverse correlation with UV luminosity as a by-product. We fit the composite spectra using stellar spectral synthesis together with models of the ISM in which a fraction f_c of the stellar continuum is covered by gas with column density N(HI). We show that the composite spectra simultaneously constrain the intrinsic properties of the stars (L900/L1500)_int along with f_c, N(HI), E(B-V), and $f_{esc,abs}$, the absolute escape fraction of ionizing photons. We find a sample-averaged $f_{esc,abs} =0.09\pm0.01$, and that subsamples fall along a linear relation $\langle f_{esc,abs}\rangle \sim 0.75[W(Ly\alpha)/110 A]$. We use the FUV luminosity function, the distribution function $n[W(Ly\alpha)]$, and the relationship between $W(Ly\alpha)$ and_out to estimate the total ionizing emissivity of $z\sim3$ star-forming galaxies with Muv<-19.5: $\epsilon_{LyC}\sim 6\times10^{24}$ ergs/s/Hz/Mpc$^3$, exceeding the contribution of QSOs by a factor of $\sim 3$, and accounting for $\sim50$% of the total $\epsilon_{LyC}$ at $z\sim3$ estimated using indirect methods.


INTRODUCTION
Substantial recent efforts have focused on establishing the demographics of star-forming galaxies in the redshift range 6 < ∼ z < ∼ 10 now believed to be most relevant for cosmic reionization (Planck Collaboration et al. 2016). Nevertheless, a detailed physical understanding of the reionization process remains elusive, due in large part to uncertainties that cannot be reduced simply by identifying a larger number of potential sources of ionizing photons. Crucial missing ingredients include knowledge of the intrinsic ionizing spectra of sources, and, more importantly, the net ionizing spectrum presented to the intergalactic medium (IGM) after passing through layers of gas and dust in the galaxy interstellar medium (ISM.) Unfortunately, these will be impossible to measure from direct photometric or spectroscopic study of reionization-era galaxies, even when the James Webb Space Telescope (JWST) comes on line, due to the rapidly increasing H I opacity with redshift along extended lines of sight to high redshifts, even post-reionization.
However, one can explore the likely behavior of Lyman continuum-producing objects, and perhaps make testable predictions, using observations of analogous sources at lower redshifts, where the opacity of intervening neutral H (H I) is less limiting, and where ancillary multiwavelength observations are more easily obtained. One avenue that has enjoyed recent success is ultraviolet (UV) observations conducted using Hubble Space Telescope (HST) of low-redshift galaxies that have been identified as likely analogs of the high redshift objects (e.g., Borthakur et al. 2014;Izotov et al. 2016a,b;Leitherer et al. 2016), such as 'Lyman break analogs" (LBAs; Overzier et al. 2009), and "Green Peas" (Cardamone et al. 2009.) These objects, although rare in the present-day universe, have many of the same properties typical of high-redshift star-forming galaxies-e.g., high UV lu-minosity, strong nebular emission lines, strong galaxy-scale outflows, compact sizes, and high specific star formation rates (sSFRs). Alternatively, one can obtain very deep observations of larger samples of more distant galaxies, where direct observations of the Lyman continuum (LyC) are possible without necessarily observing from space (e.g., Steidel et al. 2001;Shapley et al. 2006;Iwata et al. 2009;Nestor et al. 2011;Vanzella et al. 2012;Nestor et al. 2013;Mostardi et al. 2013Mostardi et al. , 2015Grazian et al. 2016.) Redshifts z > ∼ 3 bring the restframe Lyman limit of hydrogen at 911.75 Å (13.6 eV) to observed wavelengths well above the atmospheric cutoff near 3100-3200 Å, where the atmosphere is transparent, the terrestrial background is darkest, and the instrumental sensitivity of spectrometers on large ground-based telescopes is high.
In general, ionizing photons produced by massive stars in H II regions are a local phenomenon, with a sphere of influence measured in pc; the "escape fraction" of ionizing photons from an isolated ionization-bounded H II region (i.e., where the ionized region is entirely embedded within a predominantly neutral region, and the extent of the H II region is determined by the production rate of ionizing photons) is zero, by definition. However, when the density of star formation is very high, as is often the case for high redshift star-forming systems, intense episodes of star formation and frequent supernovae can in principle produce ionized bubbles that carve channels through which Lyman continuum photons might escape unimpeded. The net escaping ionizing radiation depends on the geometry of the sites of massive star formation and the surrounding ISM, the lifetimes of the stars that produce the bulk of the ionizing photons (e.g., Ma et al. 2016), and the probability that during that lifetime favorable conditions for LyC photon escape will occur.
Once H-ionizing photons escape the ISM of a parent galaxy, the probability of detection by an observer at z = 0 is governed by the effective opacity of intervening H I along the line of sight. This opacity increases steeply with redshift (e.g., Madau 1995;Steidel et al. 2001;Vanzella et al. 2010Vanzella et al. , 2012Becker et al. 2015), and for the wavelength range most relevant to LyC measurement (i.e., λ 0 < ∼ 912(1 + z s ) for a source with redshift z s ) it is also subject to large fluctuations from sightline to sightline, since it is dominated by the incidence of relatively small numbers of intervening H I systems with 16 < ∼ log (N HI /cm −2 ) < ∼ 18 (e.g., Rudie et al. 2013). As we discuss in more detail below ( §7), these competing factors strongly favor the range 2.75 < ∼ z < ∼ 3.50 for a ground-based survey.
Even within the optimal redshift range for ground-based observations, practical sensitivity limits impose severe restrictions on the dynamic range available for possible detections. The combination of limited sensitivity and detectability dominated by the stochastic behavior of the IGM foreground means that individual detections of LyC signal i.e., those for which significant LyC flux is detected without stacking) are almost guaranteed to be unusual either in their intrinsic properties, in having a fortuitously transparent line of sight through the IGM, or a combination of both. Direct detections of individual sources at high redshift Shapley et al. 2016;Vanzella et al. 2017) are valuable for demonstrating that at least some galaxies produce ionizing radiation that propagates beyond their own ISM, but they do not place strong constraints on more typical galaxies.
In view of these challenges, successful characterization of the propensity for galaxies with particular common proper-ties to "leak" LyC radiation requires observations of an ensemble, in order to marginalize over the fluctuations in the intervening IGM opacity. It should also include 1) the most sensitive possible measurements of individual sources, made as close as possible to the rest-frame Lyman limit of each; 2) a very accurate characterization of the statistics of intervening H I as a function of column density and redshift; 3) control over systematics -those affecting measurement of individual sources (e.g, background subtraction, contamination) and those that would invalidate the statistical IGM correction. The latter suggests observing sources in several independent survey fields and avoiding regions known to harbor unusual large-scale structures, if possible.
Although spectroscopic surveys at z ∼ 3 using 8m-class telescopes first became feasible in the mid-1990s (e.g., Steidel et al. 1996Steidel et al. , 2003, the initially-available instruments were not optimized for high near-UV/blue sensitivity as required for the most effective observations of LyC emission (Steidel et al. 2001). The situation changed substantially for the better with the commissioning on the Keck 1 telescope of the blue channel of the LRIS spectrograph (LRIS-B, Steidel et al. 2004) in 2002, followed by the installation of the Keck 1 Cassegrain Atmospheric Dispersion Corrector (Phillips et al. 2006) in 2007; projects demanding high efficiency in the wavelength range (3100 − 4500 Å) became much more feasible. Shapley et al. (2006) [S06], using pilot data obtained immediately following LRIS-B commissioning, presented what were apparently the first direct detections of LyC emission in the spectra of individual star-forming galaxies at z ∼ 3. S06 observed a sample of 14 LBGs with LRIS in multi-slit mode for a total of ∼ 8 hrs (in the mode sensitive to LyC light), reaching an unprecedented depth for individual galaxy spectra of 3 × 10 −31 ergs s −1 cm −2 Hz −1 at ∼ 3600 Å (3σ detection limit), corresponding to AB 27.6, or ∼ 20 times fainter than the (non-ionizing) continuum flux density of L * galaxies at z ∼ 3. While the spectra presented by S06 were far superior for LyC detection compared to what had been available, it later turned out that 2 of the 3 putative detections of residual LyC flux were due to contamination of the LyC rest-frame spectral region by faint, unrelated foreground galaxies, based on subsequent near-IR spectra and HST imaging . In the years since, it has transpired that most apparent detections of significant LyC emission from z > ∼ 3 galaxies have, on further inspection -particularly using high spatial resolution images obtained with HST -been attributed to similar foreground contamination (Vanzella et al. 2012Siana et al. 2015;Mostardi et al. 2015).
Most of the more recent observational effort toward detecting LyC emission from high redshift galaxies has been invested in imaging surveys, which have an obvious multiplex advantage, particularly when aimed at fields containing known galaxy over-densities. In such fields, narrow or intermediate band filters can be fine-tuned to lie just below the rest-frame Lyman limit at the redshift of interest (Inoue et al. 2005;Iwata et al. 2009;Nestor et al. 2011;Mostardi et al. 2013). Alternatively, one can use extremely deep broad-band UV images to search for LyC emission from galaxies having known spectroscopic redshifts that ensure the band lies entirely shortward of the rest-frame Lyman limit (Grazian et al. , 2017. While imaging surveys obtain LyC measurements for every galaxy known or suspected to lie at high enough redshift in the field of view, putative detections (and the quantification of non-detections) requires both follow-up spectroscopy and/or high-resolution HST imaging (e.g., Vanzella et al. 2010;Nestor et al. 2013;Mostardi et al. 2015).
In the present work, we return to using very deep spectroscopic observations, similar to those presented by S06. The Keck Lyman Continuum Spectroscopic Survey (KLCS) expands and improves on the S06 study in several respects: first, the sample is larger by an order of magnitude, with a total of 136 z 3 galaxies observed. Second, the observations were conducted in 9 independent survey fields, which should drastically reduce the sample variance of the results, particularly if there are large-scale correlations in intergalactic LyC opacity that could have a very strong effect on results based on a single field (see discussion in S06 10 .) Third, KLCS covers a broader range in both redshift (2.72 ≤ z ≤ 3.54), and galaxy luminosity (0.2 < (L uv /L * uv ) < 3) compared to S06. Most importantly, however, KLCS has benefited from the accumulated insight and lessons learned through experience -e.g., the importance of false positive detections due to foreground contamination and the sensitivity required for plausible detections -as well as from advances in the physical interpretation of the far-UV spectra of high redshift galaxies (e.g., Steidel et al. 2016;Eldridge et al. 2017) and in the precision of our statistical knowledge of the foreground IGM+CGM opacity (e.g., Rudie et al. 2013.) In this paper, we show that, through careful control of systematics and concerted efforts to eliminate contamination, ensembles of deep rest-UV spectra can be used to measure the ratio of LyC flux density to non-ionizing UV flux density (hereafter ( f 900 / f 1500 ) obs ) with high precision. The KLCS observations provide individual galaxy spectra of unprecedented quality; composite spectra formed from substantial subsets provide templates that are the most sensitive ever obtained for similar high redshift objects, enabling access to a remarkable range of stellar, interstellar, and nebular spectral features, many of which have not been observed previously beyond the local universe. As well as direct constraints on the leakage of ionizing photons from galaxies, high quality rest-frame far-UV spectra encode the ancillary information on the massive star populations, the geometry and porosity of the ISM, the kinematics, physics, and chemistry of galaxy-scale outflows, and stellar and ionized gas-phase metallicities of the same galaxies -all of which are needed to place LyC leakage within the broader context of galaxies and the diffuse IGM.
The paper is organized as follows: §2 describes the selection of the KLCS sample; §3 details the spectroscopic observations, while §4 describes the data reduction, including the steps taken to minimize residual systematic errors in the sample. §5 defines the final KLCS sample used for subsequent analysis. LyC measurements from individual KLCS spectra are covered in §6. §7 describes modeling of the IGM and CGM transmission used to correct the KLCS LyC observations. In §8, we form a number of KLCS sub-samples and discuss the construction of their stacked (composite) spectra; §9 relates the composite spectra of the sub-samples to the corresponding LyC measurements, while §10 discusses the implications of the measurements. In §11 we evaluate the spectroscopic results in the context of a simple model for LyC escape and its connection to other observed galaxy properties, and propose the most appropriate method for calculating the total ionizing emissivity of the galaxy population at z ∼ 3. Finally, §12 summarizes the principal results and discusses their broader implications and suggestions for future work. Appendix A summarizes data reduction steps taken to minimize residual systematic errors in KLCS as well as tests of their efficacy; Appendix B contains details of the IGM+CGM Monte Carlo transmission model used for the analysis. Where relevant, we assume a CDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and h = 0.7. All spectroscopic flux density measurements used in the paper are expressed as flux per unit frequency ( f ν ), and are generally plotted as a function of wavelength, so that a spectrum with constant f ν , constant m AB , or f λ ∝ λ −β with β = 2, appears "flat".

Target Redshift Range
Ozone in the Earth's atmosphere efficiently blocks ultraviolet (UV) radiation with a sharp transparency cutoff preventing photons with λ ≤ 3100 Å (at elevation of 4200m, as on Mauna Kea) from reaching the surface. A consequence is that ground-based observations of rest-frame LyC photons from celestial objects require observing them at redshifts z > ∼ 2.725, where the rest-frame Lyman limit of H I (λ 0 = 911.75 Å) falls at an observed wavelength of 3400 Å. At higher redshifts, observations of the rest-frame LyC benefit from the generally higher instrumental throughput and atmospheric transmission at longer wavelengths, but the sensitivity for detection of LyC flux escaping from galaxies actually declines precipitously with increasing redshift beyond z 3.5 (e.g., Madau 1995;Steidel et al. 2001;Vanzella et al. 2010;see §7) The decreasing sensitivity with increasing redshift is due to a combination of several effects: first and most obviously, galaxies of a given intrinsic UV luminosity (L uv ) become apparently fainter with redshift; in addition the characteristic L * uv (z) itself dims as redshift increases beyond z ∼ 3 (e.g., Reddy & Steidel 2009;Bouwens et al. 2010;Finkelstein et al. 2015). Although the net instrumental sensitivity at wavelengths 912(1 + z s ) Å increases with redshift from z ∼ 2.7 to z ∼ 3.5, at redshifts z > ∼ 3.5 the throughput gains are more than offset by the increasing intensity of the sky background against which any faint LyC signal must be measured.
Most importantly, line and continuum opacity of neutral hydrogen (H I) in the IGM along the line of sight increases steeply with z s (e.g., Madau 1995). Even if intergalactic H I contributed no net continuum opacity for ionizing photons emitted from a source with z s , the LyC region will be blanketed by the effective opacity caused by Lyman series lines with z < ∼ 911.8/1215.7 × (1 + z s ) −1 , reducing the dynamic range accessible to LyC detection 11 . When one includes the net LyC opacity contributed by gas outside of the galaxy, but at redshifts near enough to z s to impact the net transmission averaged over the LyC detection band, the median transmission in the rest-wavelength interval 880-910 Å decreases by a factor of 7 between z = 3 and z = 5, (see discussion in §7, and Table B2.) Tallying all of the exacerbating factors, the overall difficulty of a (ground-based) detection of LyC signal from a L uv = L * uv (z s ) increases by a factor of ∼ 35 as one moves from z ∼ 3 to z ∼ 5. Thus, there is strong impetus for a ground- 11 According to Becker et al. (2011), continuum blanketing from the Lyman α forest increases ∝ (1 + z) 2.8 over the redshift range 2.1 < ∼ zs < ∼ 5.5.
based LyC survey to focus on sources with 2.75 < ∼ z s < ∼ 3.5, as we have done for KLCS.

Survey Design
Targets for KLCS were selected in 9 separate fields on the sky (Table 1), chosen from among 25 high latitude fields in which we have obtained deep U n GR photometry suitable for selecting LBG candidates at z ∼ 3 (see Reddy et al. 2012 for a nearly-complete list) as well as spectroscopic follow-up observations. The final field selection was based on a combination of visibility time at low airmass during scheduled observing runs and the number of galaxies with previouslyobtained spectroscopic redshifts in the range 2.9 < ∼ z < ∼ 3.2 that could be accommodated on a single slit mask of the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995;Steidel et al. 2004) on the Keck 1 telescope. For six of the nine KLCS fields (see Table 1), selection of star-forming galaxy candidates using rest-UV continuum photometry was performed during the course of a survey for z ∼ 3 Lymanbreak galaxies (LBGs; see Steidel et al. 2003). Three additional fields (Q0100+1300, Q1009+2956 and HS1549+1919) including z ∼ 3 LBGs were observed during the period 2003 to 2009; these comprise part of the Keck Baryonic Structure Survey (KBSS; Steidel et al. 2004Steidel et al. , 2010Rudie et al. 2012;Steidel et al. 2014;Strom et al. 2017.) Selection of z ∼ 3 LBGs in all 9 of the KLCS fields was based on photometric selection using the 3-band (U n GR) photometric system described in detail by Steidel et al. (2003); photometric and spectroscopic catalogs for 6 of these fields (as of 2003) were also presented in that work. Full photometric and spectroscopic survey catalogs for the 3 KBSS fields will be presented elsewhere.
We used a slitmask design strategy that assigned highest priority to comparatively brighter (R ≤ 24.5) LBGs known to have redshifts in the interval 2.9 ≤ z ≤ 3.3. Other star-forming galaxies in a broader redshift range 2.7 < z < 3.6 were assigned somewhat lower priority. If space on a mask was still available, we included additional candidates that were identically selected but had not been previously observed spectroscopically 12 . Objects that had already been classified as AGN based on their existing survey spectra were deliberately given relatively high priority in the KLCS mask design, as little is known about whether ionizing radiation escapes from lowluminosity AGN or QSOs 13 . This small sub-sample will be addressed in future work. Figure 1 presents comparisons between the parent sample from which targets were selected (light histograms), and those that were successfully observed (dark histograms), in terms of redshift, apparent magnitude R, and color G − R. The "parent" sample in this case includes all galaxies with redshifts 2.7 < z < 3.6 located in the same set of fields used in the KLCS. The middle panel of Figure 1 shows that the KLCS sample has a moderate excess of sources with R < 24.5 and a related deficiency of galaxies with 25.0 ≤ R ≤ 25.5, which is an expected consequence of our observing strategy. The redshift distribution of KLCS galaxies also demonstrates our slight preference for targets in 2.9 ≤ z ≤ 3.2 redshift "window". We find no difference in the distribution of G − R color Figure 1. The distribution of redshift (top), R magnitude (center), and G−R color for the KLCS sample (dark shaded histograms) as compared with those of all confirmed LBGs in the 9 survey fields used (light shaded histograms). Apart from the slight preference for brighter objects (R < ∼ 25) and for objects at redshift z ∼ 3.0, the KLCS sample is not significantly biased with respect to the LBG population.
between KLCS and all spectroscopically identified galaxies in these fields (Figure 1, bottom panel). Thus, the sub-sample of galaxies observed for the KLCS is slightly brighter, and has a slightly tighter redshift distribution, than the LBGs in the same fields, but is otherwise representative.

SPECTROSCOPIC OBSERVATIONS
As summarized in Table 1, a total of 136 galaxies was observed, in 9 independent fields using 10 different slitmasks. Of the galaxies observed, 13 were later excluded from the LyC analysis because of uncalibrated slit defects, close companions on the slit, scattered light, or other potential sources of systematic error that would make LyC flux measurements less secure (see the detailed discussion in §5 below).
All 10 slitmasks were designed with 1 . 2 slit widths and individual slit lengths between 10 and 30 ; the median slit length for high priority KLCS targets was 20 . With the exception of the Q0933+2854 field, a single slitmask was observed in each field, containing between 8 and 16 objects known to lie at redshifts 2.7 < z < 3.6. The observations were conducted over the course of 6 separate observing runs using LRIS on the Keck 1 10m telescope between 2006 December and 2008 June. As summarized in Table 1, the total integration time per slitmask was between 8.2 and 12.8 hours. Fields were generally observed within 2 hours of the meridian in order to minimize attenuation by atmospheric extinction and (in the case of the observations made prior to 2007 August) slit losses due to differential atmospheric refraction.
All observations were made using the same configuration of the LRIS double-beamed spectrograph (Oke et al. 1995;  ). g Fields from KBSS (see also Rudie et al. 2012;Steidel et al. 2014). h Two objects are common to masks q0933_L2 and q0933_L3. i Total exposure time, in hours. Steidel et al. 2004), with the incoming beam divided near 5000 Å using the "d500" dichroic beamsplitter. Wavelengths shortward of 5000 Å were recorded by the "blue" spectrograph channel (LRIS-B) and those longward of 5000 Å by the red channel (LRIS-R). LRIS-B was configured with a 400 line/mm grism with first-order blaze at 3400 Å, providing wavelength coverage from 3100 Å to beyond the dichroic split near 5000Å. The LRIS-B detector was binned 1x2 (binned in the dispersion direction) at readout in order to minimize the effect of read noise, which for these devices is 3.8e − pix −1 , resulting in pixels that project to 0 . 135 on the sky in the spatial direction and 2.18 Å per pixel in the dispersion direction. LRIS-R was configured with a 600 line/mm grating with first-order blaze at 5000 Å, with spectra recorded using the (pre-2009) Tektronix 2k x 2k (monolithic) detector with 24 µm pixels. With no on-chip binning, the scale at the detector was 0 . 211 per pixel spatially and 1.28 Å pix −1 in the dispersion direction. The LRIS-R grating was tilted such that all KLCS slits would have wavelength coverage from shortward of the dichroic split to > ∼ 7000 Å depending on the spatial position of the slit within the LRIS field of view.
Individual exposure times were 1800s and all LRIS-B and LRIS-R exposures were obtained simultaneously, resulting in identical observing conditions and integration times for a given mask. Data were collected under mostly photometric observing conditions, and all data used in the KLCS were obtained with seeing < 1 . 0, and typically < ∼ 0 . 7. Prior to the commissioning of the Keck 1 Cassegrain Atmospheric Dispersion Corrector (ADC; Phillips et al. 2006) in the summer of 2007, all KLCS observations were obtained at elevations within 30 • of zenith and position angle close to the parallactic in order to minimize effects of differential refraction. Once commissioned, the ADC was used for all observations, greatly enhancing the efficiency of the survey by allowing position angle to be unconstrained during mask de-sign, thus allowing inclusion of a larger number of high priority targets. By correcting differential refraction before the slitmask, the ADC improved the data quality particularly for LRIS-B, and enables the use of simpler (and more robust) data reduction and extraction techniques (see §4.2).
The spectral resolution achieved varied slightly depending on observing conditions and the angular size of objects within each slit. With a typical seeing-convolved profile size of FWHM 0 . 8 − 0 . 9 for Lyman-break galaxies at z ∼ 3 (Law et al. 2007(Law et al. , 2011, the average resolving power is R 800 (LRIS-B) and R 1400 (LRIS-R) with the dispersers described above (see Steidel et al. 2010).
The dichroic split at 5000 Å typically places the location of Lyα, i.e. 1215.67(1 + z s ) Å, near the transition wavelength between the spectral channels for z ∼ 3. Longward of Lyα our primary objective was to resolve the width of typical interstellar absorption lines (FWHM 500−700 km s −1 ; see, e.g., Pettini et al. 2002;Shapley et al. 2006;Steidel et al. 2010) so that the degeneracy between velocity width and covering fraction could be disentangled. The FWHM 200 km s −1 resolution provided by the 600/5000 grating represented a compromise between resolution and sensitivity. For LRIS-B, high sensitivity (particularly in the wavelength range 3400-4000 Å) was of paramount importance, hence the choice of the 400/3400 grism, which in combination with the LRIS-B optics produces very high UV/blue throughput (see Steidel et al. 2004). The spectral resolution of FWHM 375 km s −1 is modest, but still adequate to resolve typical strong absorption and emission lines.

Calibrations
We obtained spectroscopic flat-field calibration images for LRIS-R and LRIS-B separately. Internal halogen lamp spectra provided adequate flat-fields for all LRIS-R data, but were not suitable for LRIS-B, which for our configuration requires very good flats particularly in the wavelength range 3100-4000 Å, where there are substantial spatial variations in quantum efficiency due to non-uniformities in the thinning of the silicon during manufacture.
We found from experience (for LRIS-B) that slitmask spectra of the twilight sky, obtained at similar elevation and instrument rotator angle to the science observations through the same mask, are the most effective solution; these were obtained at the beginning and end of each observing night. The twilight sky spectral flats produce adequate signal in the UV, but record the scattered solar spectrum rather than a featureless continuum. To remove the G-star spectrum but preserve the pixel-to-pixel sensitivity variations, we divide the raw flats by a spatially median-filtered 1-d spectrum calculated at regular intervals along each slit, producing images normalized to an average value of unity but retaining the desired pixel-topixel sensitivity variations. The issue of scattered light associated with flat fielding is addressed in §A.1 below.
Internal arc lamp spectra (Hg, Ne, Zn, and Cd for LRIS-B, Hg, Ne, Ar, Zn for LRIS-R) were used for wavelength calibration for both LRIS-B and LRIS-R, with 5th order polynomial fits resulting in typical residuals of 0.10 Å and 0.07 Å for LRIS-B and LRIS-R, respectively. The arc-based wavelength solutions were subsequently shifted by small amounts using measurements of night sky emission lines recorded in each science exposure.
Spectrophotometric standard stars from the list of Massey et al. (1988) were observed at the end of each night through 1 . 0 slit oriented at the parallactic angle (for all observations, both pre-and post-installation of the ADC in August 2007), with configuration settings otherwise identical to mask observations. Absolute flux calibration uncertainties are estimated to be of order 20% due to potential variations in seeing conditions (and the associated variation in slit losses) between slitmask and standard star observations. However, red-side and blue-side exposures were always obtained simultaneously, for both science and standard star observations, so that with careful reductions of the standards, the relative spectrophotometry between the blue and red channels is much more precise than the absolute spectrophotometry; the latter is relatively unimportant to our analysis (see §4.2 below).

Data Reduction
The standard LRIS spectroscopic data processing pipeline we have used for previous surveys with LRIS (e.g., Steidel et al. 2003Steidel et al. , 2004Steidel et al. , 2010 was generally followed in processing KLCS LRIS-R data. However, given the challenge of measuring very faint flux densities at levels well below the sky background in the LyC region, we paid particular attention to developing procedures for LRIS-B reductions with a goal of minimizing residual systematic errors wherever possible. Given the potential importance of systematics to the final LyC results, we describe the details of the reduction procedures up to the extraction of 1-D spectra in Appendix A. The 2-D, background-subtracted, stacked spectrograms for each sequence of LRIS-B or LRIS-R observations with a given slitmask were reduced (as described in Appendix A) so that the centroid of the trace for a given object on the slit lies at a constant spatial pixel along the dispersion direction, whether or not the observation was made using the ADC. We adopted, conservatively, a 1 . 35 boxcar extraction window (10 spatial pixels on the LRIS-B detector) to avoid making assumptions about variations of the spatial profile with wavelength -par-ticularly when the trace extends to wavelengths beyond where there is obvious detected flux. Thus, the extraction aperture for every object is a rectangular region of size 1 . 2 by 1 . 35 on the sky.
For each extracted 1-D spectrum, we used an identical extraction region on the 2-D "sky + object" 2-D spectrograms described above ( §A.5.1), i.e., and where k is the spatial position (in j pixels) of the object trace, S[i] is the resulting 1-D spectrum at dispersion point [i], and σ[i] is the corresponding one-dimensional error vector. Section A.6 describes in more detail the extent to which the noise model agrees with the data. Mask data sets obtained on more than one observing run (see Table 1) were reduced to 1-D separately, and then combined using inverse variance weighting to produce final 1-D spectra. The 1-D spectra were wavelength calibrated using the 1-D arc spectra extracted from the same region of the processed 2-D arc frames, zero-pointed using night sky emission lines measured in the extracted 1-D object+sky spectra ( §A.5.1) to correct for small amounts of flexure and slight differences in illumination between the internal arc lamps and the night sky. The final LRIS-B wavelength scales are in the vacuum, heliocentric frame with a linear dispersion of 2.14 Å pix −1 (LRIS-B), covering 3200-5000 Å.
Finally, the extracted LRIS-B and LRIS-R spectra were flux calibrated using standard stars from the list of Massey et al. (1988), obtained on the same night as the data comprising each 1-D extracted science spectrum, and corrected for Galactic extinction assuming the reddening maps of Schlegel et al. (1998). The standard star observations on the blue and red sides were made simultaneously, with the slit oriented at the parallactic angle. Because the LRIS long slit lies at a fixed position in the field of view of the instrument, the sensitivity curves in the wavelength transition region of the dichroic beamsplitter (where the spectral response of the interference coatings change rapidly with increasing wavelength from reflection to transmission) are not perfectly matched for slits located far from the focal plane position of the longslit, due to small differences in angle of incidence of the incoming light. However, through experimentation we found that accurate relative spectrophotometry could be achieved by using only LRIS-B spectra for λ < 5000 Å and only LRIS-R spectra for λ > 5000 Å; i.e., without averaging over the region of overlap.
Continuous LRIS-B+R spectra, covering 3200-7200 Å were produced by remapping the individual calibrated spectra onto a linear wavelength scale chosen so as to preserve the spectral sampling of the higher resolution red-side data, 1.0 Å pix −1 . These were used to produce composite spectra using subsets of the KLCS sample ( §8.) Thus, we made final 1-D flux calibrated spectra for each source.

DEFINING THE KLCS SAMPLE FOR ANALYSIS
In addition to the procedures described above to ensure that the observed sample is as free as possible from background subtraction systematics, we carefully examined all stacked 2-D spectrograms and extracted 1-D spectra for the sample of 136 observed galaxies to check for remaining issues that might compromise measurements of residual LyC flux. The following criteria were considered serious enough to warrant removing galaxies from the analysis sample: 1) Instrumental defects were apparent in the twodimensional spectrograms. As discussed above, the effects of non-uniform scattered light ( §A.1) and/or irreparable slit illumination irregularities ( §A.3) can negatively impact the quality of the flat-fielding and background subtraction stages. There were 5 cases in which such effects were judged to be relevant, 1 in each of the Q0100, Q0256, and B20902 fields, and 2 on mask q0933_L2 (pre-ADC) in the Q0933 field; all 5 were removed from the KLCS analysis sample.
2) The primary target galaxy appeared to have another object near enough on the slit that the light from the two objects could not be separated with confidence in the 1-D extraction; 3 such cases were identified (one object from each of Q0100, B20902, and Q0933 [mask q0933_L2]), and removed from the sample.
3) A target was revealed to be a spectroscopic "blend" of unrelated objects, where the foreground object has the potential to cause a spurious (false-positive) detection of LyC flux if the source redshift is assumed to be the higher of the two redshifts. We identified 7 cases of spectroscopic blends, of which 5 were removed from the sample: Q0100-C1 (z = 3.44/2.21; see Figure 2), Q0256-d4 (z = 3.67/2.63; see Figure 2), Q0256-m11 (z = 3.090/2.09), Q0256-md34 (z = 3.130/1.870), and Q1009-C41 (z = 3.62/3.22/1.9). In the remaining two cases, multiple redshifts were identified in the spectrum but the lower of the two redshifts was sufficiently high to allow a LyC measurement (i.e., z ≥ 2.75): these are Q0933-M23 (z = 3.380/3.289) and Q1549-D25 (z = 2.822/2.755; see Figure 2). In both cases the target was retained in the sample, but the lower of the two redshifts was assigned for subsequent analysis. The 5 discarded objects, had they been analyzed without the spectroscopic identification of the potential contamination, would all have been classified as LyC detections, with significance ranging from 3σ − 7σ and flux level f 900 0.050 µJy (m AB [LyC] 27.2).
Thus, in total, 13 of 136 galaxy targets were removed from the final KLCS sample 14 . The 124 galaxies remaining in the final sample are listed (along with their coordinates, redshifts, and optical photometry) in Table 2. 6. MEASUREMENTS OF INDIVIDUAL KLCS SPECTRA 6.1. Residual Interloper Contamination As discussed in §1, foreground contamination leading to false-positive detection of LyC flux is a major concern for any putative detection of ionizing radiation from high redshift galaxies, and recent experience has shown that candidate LyC detections based on ground-based imaging surveys should be viewed with caution until additional observations -particularly multi-band HST imaging -can confirm the association of the detected flux with the known z ∼ 3 galaxy or is more likely to be due to an unrelated object at a different (lower) redshift. Three examples of targets identified as spectroscopic blends of two redshifts within the 1-D extraction footprint of the primary target. Each panel has line identifications marked with colors corresponding to the two redshifts indicated (lower redshift in red, the higher redshift in blue). Since the lowerredshift features in the top 2 panels are sufficiently low that non-ionizing flux would contaminate the rest-frame LyC region of the higher redshift source, they were removed from the analysis sample. The bottom panel, which shows a spectroscopic blend of galaxies with z = 2.756 and z = 2.822, was retained in the sample but assigned the lower of the two redshifts.
In this section, we estimate the likelihood that our triage of the 1-D KLCS spectra has successfully identified most of the contaminated sources that would lead to false positive detections of LyC signal.
Recently, imaging surveys for LyC (Siana et al. 2007;Inoue & Iwata 2008;Nestor et al. 2011;Vanzella et al. 2012;Mostardi et al. 2013) have used Monte Carlo simulations based on the surface density of objects in a deep U band images to estimate the probability that a random faint galaxy with (e.g.) z < 2.75 would fall close enough to the centroid position of a targeted galaxy to cause a false-positive LyC detection. The probability of a chance superposition in- creases as the sensitivity limit for LyC detections becomes more sensitive; any galaxy with UV continuum surface brightness exceeding the typical statistical detection threshold over the bandpass used for LyC sensing is a potential source of contamination.
For the particular case of UV-color-selected LBGs at z ∼ 3 -which require the presence of a photometric break between the observed U n (3550/600) and G (4730/1100) and a relatively blue color between G and R (6830/1000) to have been selected for observation in the first place -the most likely sources of contamination are flat-spectrum (i.e., zero color on the AB magnitude system) galaxies with apparent magnitudes bright enough to yield a photometric or spectroscopic detection at λ obs 3500 Å, but faint enough that the resulting perturbation of the U n GR photometry does not scatter the object out of the color selection window. The KLCS spectroscopically-observed sample has apparent U n ≥ 26.13 15 , median U n = 27.7, and typical spectroscopic detection limits in the LyC detection band (rest-wavelength range [880,910] Å) of σ 900 0.01µJy (10 −31 ergs s −1 cm −2 Hz −1 , or m AB [LyC] = 28.9; see Figure 4 and Table 2.) If we consider objects in the magnitude range 26 < m AB (3500) < 28.0 as the most likely potential contaminants of the KLCS sample, we can use our deepest available UV images obtained under seeing conditions similar to those of the KLCS spectroscopic observations 16 to estimate the probability of contamination of the KLCS slit apertures. The average surface density of detections in the range 26 ≤ m AB [3500]Å ≤ 28 is Σ avg = 88.7 ± 2.4 arcmin −2 , where the uncertainty represents the scatter in Σ avg among the 4 fields. Each KLCS spectroscopic extraction aperture subtends 1 . 2 by 1 . 35, or a solid angle of 1.6 arcsec 2 (4.4×10 −4 arcmin −2 ). The probability that the centroid of a m AB = 26 − 28 object (at any redshift) falls within the extraction aperture for any single KLCS target is P 88.7 × 4.4 × 10 −4 = 0.039; thus, in a sample of 128, we 15 Only 15% have Un < 27. 16 These include the NB3420 image in the HS1549 field (Mostardi et al. 2013), NB3640 in SSA22a (Nestor et al. 2011), and the Un images in the Q1422+23 field ) and Q0821+31 fields (KBSS; Rudie et al. 2012;Reddy et al. 2012;Steidel et al. 2014.) expect N c ∼ 128 × 0.039 5 will be affected by such contamination.
While this type of argument cannot exclude the possibility that there remain unidentified false-positive detections in the KLCS sample, the fact that the number of slit apertures predicted to be affected by chance superposition of UV-faint foreground galaxies is similar to the number of spectra identified as blends with foreground galaxies suggests that our spectroscopic "triage" has likely removed most of the contaminants that would lead to false positive LyC detections.

Galaxy Systemic Redshifts
Because our primary goal is to measure the residual flux averaged over a relatively broad window in rest-wavelength, precise systemic redshift measurements are not critical. However, since we will be combining individual spectra into composite stacks ( §8), we assigned our best estimate of z s for each galaxy based on the information in hand. Of the 124 sources in the KLCS analysis sample, 55 (44%) have measurements of nebular [O III] emission lines in the rest-frame optical obtained using Keck/MOSFIRE (McLean et al. 2012;Steidel et al. 2014.) In all such cases, the measured z neb is assumed to define the systemic redshift z sys , with an uncertainty of 15 − 20 km s −1 .
For objects lacking nebular emission line measurements, we used estimates of z sys based on the full KBSS-MOSFIRE sample Strom et al. 2017) with z > 2.0 and existing high-quality rest-UV spectra. These were used to calibrate relationships between z sys and redshifts measured from features in the rest-frame FUV spectra, strong interstellar absorption lines (z IS ) and/or the centroid of Lyman-α emission (z Lyα ). As for previous estimates of this kind (e.g., Adelberger et al. 2003;Steidel et al. 2010;Rudie et al. 2012), we adopt rules that depend on the particular combination of features available in each spectrum, where spectra fall into one of 3 categories: (a) those with measurements of both z Lyα and z IS ; (b) those with z IS but without z Lyα (generally because the Lyα feature appears in absorption); and (c) those with only z Lyα available.
For galaxies in category (a), comprising 60% of the KLCS sample, for category (b) ( 30% of the KLCS sample) : with ∆v IS = 110 km s −1 ; for category (c) ( 10% of the KLCS sample): with ∆v Lyα = 235 km s −1 . The residual uncertainties after applying these rules are σ 110, 130, and 130 km s −1 (rms) for spectra in categories (a), (b), and (c), respectively, estimated from the subset of KLCS with measurements of z neb . The resulting value of z sys is given for each galaxy in Table 2.

LyC Measurements
All quantitative measurements or limits on the flux density of residual LyC emission have been evaluated over a fixed bandpass in the source rest frame, placed just shortward of the Lyman limit: We also define a rest-frame bandpass to represent the FUV non-ionizing flux density, The choice of the interval [880,910] for the LyC measurement is a compromise based on two considerations. First, ideally the LyC measurement should be made as close as possible to the rest-frame Lyman limit of the galaxy, and integrated over the smallest wavelength range feasible given the noise level of the spectra, in order to minimize the effects of H I in the IGM along the line of sight. As discussed in detail in §7, the opacity of the IGM to LyC photons is the largest source of uncertainty in the measurement of the emergent ionizing photon flux from a galaxy. The mean free path of H-ionizing photons emitted at z s = 3.05 corresponds to a redshift interval of ∆z ∼ 0.25 (e.g., Rudie et al. 2013), or to ∼ 60 Å in the rest frame; i.e., the flux density at λ 0 850 Å is reduced by an average factor of e (∼ 2.72) relative to its emergent value. By using the [880,910] interval, to a first approximation the average IGM LyC optical depth due to intervening material is τ igm < ∼ 0.5. An additional consideration is the wavelength range over which the Keck/LRIS-B system sensitivity remains high. Although the UV sensitivity of Keck/LRIS-B is the highest among instruments of its kind , it decreases with wavelength shortward of 3500 Å, particularly when atmospheric opacity is included. The lowest redshift source included in KLCS has z s = 2.718 (Q0100-MD6), placing λ 0 = 880 Å at an observed wavelength of 3272 Å , where the LRIS-B end-to-end throughput with the d500 beamsplitter has dropped to 20% from > ∼ 50% near 4000 Å. Including the atmosphere, these numbers reduce to ∼ 11% and ∼ 40% at a typical airmass of 1.10 on Mauna Kea. At the mean redshift of the KLCS sample ( z s = 3.05), [880,910] corresponds to an observed-frame interval [3560,3690] Å, where the instrumental throughput averages 35%.
Thus, [880,910] is observable over the full range of KLCS, and is narrow enough to minimize the IGM opacity against which escaping flux must be detected, but sufficiently broad to allow for improved photon statistics while reducing the fluctuations due to discrete H I absorption lines arising from intervening material, superposed on the rest-frame LyC. Figure 3 presents the measured values of f 900 , with objects grouped and ordered according to slitmask and the relative physical position on the mask (along the slit direction) for each target. Targets having > 3σ detections of f 900 are labeled. The values of f 900 and their associated statistical error (σ 900 ) for KLCS galaxies are listed in Table 2. Also given in Table 2 are measurements of the non-ionizing UV continuum flux density f 1500 (equation 7) measured directly from the individual 1-D flux-calibrated spectra. Each slit mask used has its own panel, with the order of objects following their relative physical placement on the slitmask. The total exposure time in hours for each mask is indicated under the mask name. The errors shown are ±1σ statistical errors derived accounting for counting statistics and detector noise. A total of 15 out of 124 galaxies are nominal detections, with f 900 > 3σ 900 ; these are marked with square (blue) symbols. Note that galaxies Q0933-MD75 and Q0933-D20 were observed twice, on masks q0933_L2 and q0933_L3-points labeled in red represent the weighted average measurement from the two independent spectra.
The measurements and uncertainties for f 900 and f 1500 [and their ratio ( f 900 / f 1500 ) obs ] were obtained directly from the 1-D spectra and associated error arrays based on the noise model  Shapley et al. (2016). described in §A.5.1; we have made no attempt to apply aperture corrections to the spectroscopic flux density measurements (the photometric measurements in the U n GR system are also provided in Table 2) but we believe that the relative spectrophotometry of the 1D spectra has systematic errors < ∼ 5%; uncertainties in the calibration of LRIS-R relative to LRIS-B spectra may contribute additional < ∼ 10% systematic uncertainty in ( f 900 / f 1500 ) obs .
The flux error-bars in Figure 3 are statistical errors based on the noise model presented in §A.5.1, with typical values of σ 900 = 0.01 µJy; most objects on each mask have f 900 measurements consistent with zero to within 1-2σ 900 .
Grouping observations by slit-mask allows us to monitor the presence of systematic errors from mask to mask. It is also a useful way to inspect the data for correlations with object position on a given slit-mask. It is apparent from Figure 3 that mask B20902-L1 (and possibly others) has residual systematic errors manifesting as correlated behavior of f 900 as a function of position on the slitmask. Although the maximum deviation from zero in the b20902_L1 panel is only ∼ 2σ 900 , there appear to be systematic undulations relative to zero with amplitude comparable to the random errors. The systematics appear to have larger excursions to f 900 < 0, as might occur from over-subtraction of the background due to contamination of slit regions used to model it by either scattered light (judged to be the dominant factor for mask b20902_L1), or by contributions from unmasked sources falling along the slit.
Systematic over-subtraction of the background level was also noted for the sample of 14 sources presented by Shapley et al. (2006) (S06) -these authors found that the measured f 900 for objects without significant detections was centered around an unphysical negative value (see their Figure 5). As discussed in detail in §4.2 (see also Appendix A.6), considerable effort was invested in improving the flat fielding (with its tendency to exacerbate problems with non-uniform scattered light ( §A.1) and 2-D background subtraction. The tests summarized in appendix A.6 suggest that these were generally successful. We show below that the procedures have also reduced the residual systematic errors in extracted 1-D spectra to a level significantly smaller than the random errors.
To verify that this is the case, we excluded all galaxy measurements for which | f 900 | > 3σ 900 where σ 900 is the statis-tical error estimate. For the remaining sample of 106 measurements, f 900 /σ 900 = 0.31±1.33, with median f 900 /σ 900 = 0.44, and inter-quartile range [-0.50,+1.41]. Excluding the two masks (B20902-L1 and Q1009-L1) with the most obvious systematic issues, the results remain largely unchanged: f 900 /σ 900 = 0.28 ± 1.32, with median 0.44 and inter-quartile range [-0.66,+1.25]. Under the null hypothesis that f 900 = 0 and that systematic errors are small compared to normally distributed random errors, one expects f 900 /σ 900 = 0.0 ± 1.0. As we shall see, the true values of f 900 are likely greater than zero for some fraction of the sample with | f 900 /σ 900 | < 3; since the standard deviation of f 900 /σ 900 is only ∼ 30% larger than expected under the null hypothesis, we believe that residual systematic errors in f 900 are sub-dominant compared to statistical errors. In the remainder of the analysis below, we continue to assume that our model for random errors in the spectra ( §A.5.1) accurately describes the uncertainties. Thus, 15 of 124 galaxies have f 900 > 3σ 900 , which henceforth are referred to as "detected"; their properties are listed individually in Table 3. The same objects are marked using blue squares in Figure 3. Note that one of the objects in Table 3, Q1549-C25 [( f 900 / f 1500 ) obs = 0.08], has been discussed in detail by Shapley et al. (2016); in addition to the LyC flux measurement, it has also been observed with multi-band HST imaging (see Mostardi et al. 2015) which indicates no evidence for a contaminating foreground source that might have caused a false-positive LyC detection 17 . The only other confirmed LyC detection of a z ∼ 3 galaxy is "Ion-2" (z s = 3.22), which has ( f 900 / f 1500 ) obs 0.13 and rest-frame UV luminosity ∼ L * uv de Barros et al. 2016). Most of the objects in Table 3 have ( f 900 / f 1500 ) obs similar to the two confirmed LyC detections 18 .  Figure 4 shows all 124 KLCS targets on a plot of σ 900 versus redshift. Note that σ 900 shows a trend of increasing toward lower redshift; this is due to the decreasing instrumental sensitivity at rest-frame wavelengths [880,910] Å which (as we have seen) changes by a factor of a few over the range 2.75 < ∼ z < ∼ 3.4 (where [880,910] falls at observed wavelengths 3300 < ∼ (λ/Å) < ∼ 4000.) Figure 5 shows a combination of formal detections andfor the non-detections -the adopted 3σ lower limits on the ratio ( f 1500 / f 900 ) obs using the measurements and uncertainties given in Table 2. Note that the formally detected objects lie well within the distribution of 3σ upper limits of the nondetections. The implications of these results for the distribution of LyC flux within the observed sample requires a quantitative assessment of the extent to which H I gas outside of galaxies (but along the line of sight) has censored our ability to detect LyC flux when it is present; we address this issue in §7.
6.4. LyC Detection vs. Other Galaxy Properties In this section we briefly review the properties of galaxies with individual LyC detections versus the majority that do not have individual detections. We argue below ( §8) that unambiguous interpretation of LyC detection statistics requires the use of ensembles of galaxies. Figure 6 compares the distributions in 3 empirical galaxy properties for the LyC-detected and LyC-undetected subsamples of KLCS (with 15 and 109 objects, respectively.) The left-most panel of Fig. 6 shows the rest-UV luminosity distribution, relative to L * uv in the far-UV (rest-frame 1700 Å) luminosity function at z ∼ 3 (Reddy & Steidel 2009). As discussed in §5 above, most KLCS galaxies have L uv within a factor of a few of L * UV . The sub-sample with formal LyC detections occupies a slightly narrower range of luminosity, though a twohas been confirmed spectroscopically, with ( f 900 / f 1500 ) obs = 0.052 ± 0.011. After accounting for the factor of ∼ 1.7 lower 90th-percentile transmission at zs = 4 compared to zs = 3.05 (see Table B2), this value is roughly equivalent to a measurement of ( f 900 / f 1500 ) obs 0.088 ± 0.019 at zs = 3.05, entirely consistent with the typical KLCS detection listed in Table 2.
sample KS test shows that the two luminosity distributions are statistically consistent with being drawn from the same parent population 19 .
As one of the most easily-observed and measured spectroscopic characteristics of high redshift star-forming galaxies, the rest-frame equivalent width of the Lyman-α emission line, W λ (Lyα), is a useful diagnostic, and has been shown to correlate strongly with other more subtle spectroscopic features present in the spectra of LBGs Kornei et al. 2010). We measured W λ (Lyα) for the KLCS galaxy sample following the method described in Kornei et al. (2010); the values are listed in Table 2. The center panel of Figure 6 (see also Table 2) shows the distribution of W λ (Lyα) for the full KLCS galaxy sample, divided according to whether or not they are formally detected in the LyC band [880,910]. Although the detected sub-sample tends to have Lyα in emission, and the fraction of objects with LyC detected appears to be correlated with Lyα equivalent width, a two-sample KS test cannot reject the null hypothesis that the sub-samples are drawn from the same parent population.
As described in more detail elsewhere Steidel et al. 2003;Reddy & Steidel 2009) the z ∼ 3 LBG U n GR color selection imposes small systematic differences in the redshift selection function depending on the intrinsic galaxy properties, in the sense that intrinsically redder galaxies are less likely to satisfy the selection criteria at the high redshift end of the selection window. The main reason for this is increased line blanketing from the Lyα forest with redshift. Galaxies with very strong Lyα features (in either emission or absorption) affect the observed G − R color for similar reasons, since Lyα falls in the observed G band throughout the KLCS redshift range. However, since the full KLCS sample has spectroscopic measurements, we can correct the observed G − R color of individual galaxies for both effects, thereby yielding estimates of the intrinsic UV continuum color. Figure 7 shows the measurements for individual KLCS galaxies as a function of redshift in terms of (G − R) 0 , the proxy for continuum color after correction for the mean IGM absorption in the G band and the individual W λ (Lyα) . The KLCS galaxies with individual > 3σ LyC detections are circled, and their distribution is compared with the full sample in the rightmost panel of Figure 6.
As for L uv and W λ (Lyα), the sub-sample having direct individual detections is consistent with being drawn from the same distribution in (G − R) 0 as the full KLCS sample. We discuss the statistical connections between LyC leakage and galaxy properties in more detail in §8.

THE OPACITY OF THE INTERGALACTIC MEDIUM
The opacity of the IGM has been reasonably wellquantified in a statistical sense (e.g., Madau 1995;Faucher-Giguère et al. 2008;Prochaska et al. 2009;Rudie et al. 2013;Inoue et al. 2014) from high resolution spectroscopic surveys of relatively bright QSOs. Modeling the statistics of IGM absorption is essential for understanding the implications of any survey seeking to quantify the intensity of ionizing photons escaping from high redshift galaxies. In order to convert our observations of ( f 900 / f 1500 ) obs into the more relevant spectrum of emergent ionizing radiation from galaxies, we used a set of IGM attenuation models using a Monte Carlo technique described by Nestor et al. (2011) (see also Bershady et al. 1999;Figure 6. (Left:) Distribution of rest-frame UV luminosity relative to L * UV at z 3 for galaxies with individually detected f 900 (dark/blue) compared to that of the sample with f 900 < 3σ 900 (light/yellow). (Center:) As for the lefthand panel, comparing the distribution of rest-frame equivalent width of Lyα, W λ (Lyα). (Right:) As for the lefthand panel, comparing the distribution of rest-UV continuum color (see §6.4). 2-sample Kolmogorov-Smirnov (KS) tests applied to all 3 cases cannot significantly reject the null hypothesis that the distributions of formally detected and formally undetected objects are drawn from the same parent population. , with H I distribution function ( f (N HI , X)) parameters updated based on the KBSS QSO sightlines, including the effects of the CGM (Rudie et al. 2013.) The models produce full realizations of individual IGM sightlines toward a source with redshift z s by drawing from the empiricallycalibrated incidence of intervening H I as a function of N HI and redshift, over the range 12 ≤ (log (N HI /cm −2 ) ≤ 22.0 and 1.6 ≤ z ≤ z s , the details of which are presented in Appendix B. Each simulated spectrum includes both line blanketing from Lyman series absorption lines (most relevant for the low-N HI systems) and LyC opacity (dominated by systems having log N HI 16 − 18) as a function of observed wavelength over the range 3100 ≤ λ obs ≤ 1216(1 + z s ). By creating an ensemble of simulations with the same redshift distribution as the sources in the observed KLCS sample, one can make precise statistical statements about the effect of the IGM on the measurements of f 900 .
Monte Carlo simulations were run for two separate models based on the parametrization of f (N HI , X) ( Figure B1.) The first, which we call "IGM-Only", assumes that every sightline to a KLCS source is equivalent to randomly selecting H I absorbers from the distribution function in a way that depends only on redshift and N HI , i.e. the "Average IGM" parametrization in Figure B1. Rudie et al. (2012) showed that regions within 300 kpc (physical) of galaxies at z ∼ 2.4 give rise to a significantly higher incidence of logN HI > ∼ 14 absorption. To account for this, a second set of realizations, called "IGM+CGM", includes a model for regions of enhanced H I absorption arising in the circum-galactic medium (CGM) near the source galaxies. As discussed in Appendix B, the CGM H I absorber frequency distribution function f (N HI , X)[CGM] ( Figure B1) is based on the KBSS survey results (Rudie et al. , 2013 detailing the distribution of N HI along sightlines passing within 50-300 physical kpc (pkpc) and within 700 km s −1 in redshift (i.e., |∆z| ≤ 0.0023(1 + z s ) of spectroscopically identified galaxies with 2.0 < ∼ z < ∼ 2.8. For modeling lines of sight in the "IGM+CGM" simulations we used the CGM parameters for f (N HI , X) (Table B1) for redshifts z s − ∆z ≤ z ≤ z s and the "Average IGM" formulation z < (z s − ∆z), where ∆z = 0.0023(1 + z s ) (c∆z = 700 km s −1 .) For most purposes in this paper, the "IGM+CGM" is the more relevant of the two; results from both are included in Table B2. The differences between "IGM+CGM" and "IGMonly" opacity model are discussed below.
Since we have adopted [880,910] Å in the rest frame of the source as our LyC measurement band, of greatest interest is the prediction for the statistical reduction by the CGM+IGM of the flux density at observed wavelengths 880(1 + z s ) ≤ λ obs ≤ 910(1 + z s ). We define With this definition of t 900 , we can correct observed values of f 900 / f 1500 obs for IGM(+CGM) attenuation where f 900 / f 1500 out is the emergent flux density ratio that would be measured by an observer at z = 0 if there were no opacity contribution from the CGM and IGM along the line of sight; as discussed in more detail below ( §10, §11.4),  Table B2), the dashed black curve is the mean transmission, and the magenta solid curve is 1 − D B , the mean transmission between Lyβ and the Lyman limit in the rest frame of the source (see discussion in AppendixB). The KLCS redshift range is (lightly) shaded with a rectangular box. Right: As in the lefthand panel, for the "IGM+CGM" opacity model.
f 900 / f 1500 out is the quantity relevant to calculations of ionizing emissivity of galaxy populations. Figure 8 illustrates how the distributions of "IGM-Only" and "IGM+CGM" transmission t 900 depend on z s , with the range selected for KLCS shaded yellow. Table B2 summarizes the results of the IGM modeling in terms of the percentiles (10, 25, 50, 75, and 90th) of the IGM or IGM+CGM transmission at particular rest wavelength intervals of interest, all as a function of source redshift. Source redshift values were modeled using small (∆z = 0.05) increments over the KLCS redshift range, but we have also included values for sources with z s > 3.5 to provide some intuition about the rapid decline in IGM transmission as redshift increases 20 .
Even within the redshift range spanned by the KLCS sample, t 900 varies by almost a factor of 2; however, if we treat the KLCS sample as an ensemble of 124 sightlines with 2.75 < ∼ z s ≤ 3.55, the IGM+CGM simulations predict that the ensemble average t 900 = 0.371 with 68% confidence interval t 68 = [0.353, 0.392], i.e. an uncertainty of ∼ 5% on t 900 for the ensemble 21 . Fig. 9 shows the full distribution of t 900 expected for ensembles of 124 sightlines with the same redshift distribution as KLCS, for both the IGM Only and IGM+CGM opacity models.
One can see from Figure 9 that there are two main effects of including the opacity of the CGM: (1) it increases by a factor of ∼ 3 the number of sightlines with very low transmission (t 900 < 0.05); and (2) it significantly decreases the fraction of sightlines expected to have transmission near the maximum of t 900 0.6. The latter would be (all other factors being equal) the most likely to yield detectable LyC signal in the spectra of individual galaxies. This point illustrates the importance of accounting statistically for LyC attenuation by gas which is 20 It is likely that our Monte Carlo simulations underestimate the LyC opacity for z > 3.5, since our assumption about the evolutionary parameter γ = 1.0 describes the incidence of LLSs well over the range 2 < z < 3.5 but the slope appears to steepen to γ 2 by z ∼ 4 (Prochaska et al. 2009;Songaila & Cowie 2010). 21 The expected t 900 is very close to the average expected if all galaxies had zs = 3.05, the mean redshift of KLCS. Figure 9. Histogram of the probability density function of t 900 , the net transmission in the rest wavelength interval [880,910] for sources having identical zs distribution as the KLCS ensemble. The histograms show the average of 1000 sets of 124 sightlines: blue-hatched includes only IGM opacity while the red-hatched includes IGM+CGM opacity. The continuous curves, color coded in the same way, show the cumulative fraction of sightlines with lower t 900 (both curves refer to the righthand axis). Note that the main effect of including the CGM is to flatten the distribution for t 900 > ∼ 0.05 while increasing the expected number of opaque sightlines (t 900 < 0.05) by a factor of 3 and decreasing the relative number of sightlines with the highest t 900 . outside of the galaxies, but near enough to be observationally indistinguishable from a case of zero LyC photon escape from the galaxy ISM.
For some purposes (see §8), it is useful to form ensemble average transmission spectra covering a wider range in restwavelength than we have used to measure t 900 . Figure 10 illustrates the rest-wavelength dependence of ensemble average transmission spectra for z s = 3.05. These spectra were created b Given an ensemble of N realizations of t900 for source redshift zs, the ratio of the half-width of the 68% confidence interval of t900 and the true value of t900 .
from an ensemble of 10000 realizations of the full transmission spectra predicted by the Monte Carlo IGM+CGM model, sorting them by t 900 , and averaging in percentile bins. shows that for > ∼ 10% of sightlines at z s = 3.05 ("Top 10%" spectrum in the figure) there will be little or no attenuation over and above that produced by line blanketing. However, given typical detection limits for LyC flux, at least the "Bottom 25%" shown in Figure 10 will appear be be optically thick to their own LyC radiation, due entirely to attenuation by intervening H I at D > 50 kpc arising due to the opacity of the combined CGM+IGM. 7.1. IGM Transmission: Sampling Issues Assuming the "IGM+CGM" opacity model, at z ∼ 3, the 68% confidence interval for a single realization of t 900 at z s 3.05 is t 900 = [0.038, 0.589]; the corresponding 68% confidence intervals are t 900 = [0.076, 0.665] at z s = 2.70 , and t 900 = [0.017, 0.471] at z = 3.50. The implications of these broad distributions are worth stating explicitly: the LyCdetectability of a galaxy with a given intrinsic (i.e., emergent) ( f 900 / f 1500 ) out is to a great extent controlled by the statistics of the IGM transmission, which is uncertain by factors of between 5 and 10 even over the limited redshift range 2.70 < ∼ z s < ∼ 3.50. Conversely, a single measurement of ( f 900 / f 1500 ) obs cannot be converted into an intrinsic property of the source, for the same reason (see, e.g., Vanzella et al. 2015;Shapley et al. 2016.) However, t 900 (z s ) , the mean transmissivity of the IGM for a source at z = z s , and its uncertainty δt 900 (z s , N) can be quantified for an ensemble of N sightlines using the Monte Carlo models described in Appendix B. For example, assuming z s = 3.0, the number of independent IGM sightlines one must sample in order to reduce δt 900 (3.0, N) to less than 10 (5) percent of t 900 (3.0) is N = 36 (150). In other words, if the spectra of 36 sources at z s = 3.0 are averaged to produce a measurement of f 900 / f 1500 obs , then ( f 900 / f 1500 ) out = ( f 900 / f 1500 ) obs /(t 900 ± δt 900 ), and the IGM correction contributes a fractional uncertainty to the inferred emergent flux density ratio of 10%. Some example statistics relevant for the KLCS redshift range and sample size are given in Table 4.
This type of analysis is useful when one has observations of a particular class of object (grouped by known property, e.g. L uv , z s , W λ (Lyα), inferred extinction) forming a subset of the full sample: as long as the subset has sufficiently large N, the statistical knowledge of the IGM opacity can be used to derive the average intrinsic LyC properties of that class. The validity of this procedure depends on (1) the assumption that each line of sight in the sample is uncorrelated with any other line of sight in the same ensemble and (2) that the IGM+CGM opacity model is an accurate statistical description of the true opacity.
The first assumption -that lines of sight are independentis almost certainly not valid when a survey is conducted in a single contiguous field of angular size ∼ 10 (transverse scale of ∼ 5 pMpc at z ∼ 3), as has often been the case for reasons of practicality (e.g., Shapley et al. 2006;Nestor et al. 2011;Mostardi et al. 2013Mostardi et al. , 2015Vanzella et al. 2010;Siana et al. 2015). Correlated sightlines could be especially problematic in fields known to contain significant galaxy over-densities at or just below the source redshifts. As has been discussed by (e.g.) Shapley et al. (2006): if observed sources are located in regions containing more (or less) gas near N HI = 10 17 cm −2 than average, their "local" IGM could skew significantly away from expectations if an "average" line of sight is assumed. The effect is likely to be negligible for Lyα forest blanketing, but could strongly influence f 900 , since it relies on the statistics of small numbers through the incidence of relatively high column density H I over a small redshift path (∆z 0.14 for z s = 3.05 for our LyC region [880,910] Å.) The full KLCS sample is relatively insensitive to the effects of correlations between IGM sightlines by virtue of the fact that it is comprises 9 independent survey regions. Similarly, the IGM opacity model is based on the statistics of 15 independent QSO sightlines in the KBSS survey (Rudie et al. 2013), and the CGM corrections to the IGM model are based on regions near galaxies selected using essentially identical criteria to those used for KLCS (albeit at slightly lower redshift; see Rudie et al. 2013).
Nevertheless, it is worth pointing out caveats associated with the CGM+IGM opacity models possibly relevant even for KLCS. First, the adopted IGM+CGM opacity model probably under-estimates the CGM opacity, since it is based on lines of sight to background QSOs with 50 < D tran ≤ 300 pkpc (i.e., projected angular distances 6 < θ < 37 ), which by virtue of the cross-section weighting correspond to an average impact parameter of D tran > ∼ 200 pkpc. On the other hand, every line of sight to a (source) galaxy includes gas with physical distance from the source of 50-300 pkpc. If N HI continues to increase with decreasing galactocentric radius (as is likely), the transverse sightlines used to estimate the CGM contribution for the opacity model would systematically underestimate the total CGM opacity. While the CGM+IGM corrections applied to the galaxy spectra in the KLCS sample are likely to be appropriate for the range of galaxy properties we are sensitive to in this work (due to the similarity in the galaxy properties between KBSS and KLCS mentioned above), it could be dangerous to apply the corrections to sources selected using substantially different criteria -for example, if most of the ionizing radiation field is produced by objects with a different overall environment (e.g., low-density regions harboring fainter sources might be expected to be surrounded by less intergalactic and circumgalactic gas; see Rakic et al. 2012;Turner et al. 2014) then the CGM contribution to the net opacity toward those sources might be over-estimated.
Finally, the absence of a clear distinction between ISM, CGM, and IGM leads to an issue of semantics: in considering the "escape" of ionizing photons, one must also define what constitutes escape, i.e., how far must the ionizing photon travel before it is counted as having escaped, and what is the probability that it will be observable? For the remainder of this paper, we assume that ionizing photons absorbed within a galactocentric radius of r = 50 pkpc have by definition not escaped. We then use the IGM+CGM statistics outlined in Appendix B to correct the observations back to the r 50 pkpc "surface" -our working definition of a galaxy's "LyC photosphere". 7.2. LyC Detectability: Spectroscopy vs. Imaging There are distinct practical advantages in using a measurement band that samples a fixed bandwidth in the rest frame of each source, 880-910 Å, placed just shortward of the intrinsic Lyman limit. It is possible to use images taken through a comparably-narrow bandpass (i.e., 100 − 200 in the observed frame; e.g., Inoue & Iwata 2008;Nestor et al. 2011;Mostardi et al. 2013), but a narrow-band filter will be optimally-placed only for sources at a fixed redshift, which as discussed above can also be problematic in terms of sample variance due to gas-phase overdensities and/or correlated sightlines arising in the intra-protocluster medium.
Similar issues affect observations where a broad-band filter is used for the LyC detection band, such as surveys conducted using the F336W filter in the WFC3-UVIS camera on board HST (e.g., Mostardi et al. 2015;Siana et al. 2015;Vasei et al. 2016). In this case, since the broad bandpass must lie entirely shortward of the rest-frame Lyman limit of sources to provide unambiguous detections of ionizing photons, it is confined to source redshifts z s > 3.07; however, even for z s 3.1, the photon-number-weighted mean flux density measured in the filter bandpass is considerably reduced relative to the rest-frame [880,910] wavelength interval. An example assuming z s = 3.10 is shown in Figure 11, where t 900 = 0.372 (i.e., close to the mean value for random IGM+CGM sightlines to sources with z s = 3.10; Table B2), but the F336W Figure 11. An example of a single IGM sightline, simulated using the Monte Carlo models described in the appendix, for a source at z = 3.10. The IGM transmission vs. observed wavelength is shown, with the position of the restframe Lyman limit indicated with the vertical dashed red line. The light yellow shaded region illustrates the range of observed wavelength used to calculate t 900 , which in this realization has t 900 = 0.372 -close to the mean value expected for z ∼ 3.1 (see Table B2). The violet shaded region is the filter response function for the HST WFC3/UVIS F336W filter, used for many studies of LyC leakage at z > ∼ 3. The mean transmitted flux integrated over the F336W bandpass has t F336W = 0.045, i.e., a factor of > ∼ 8 lower than the spectroscopic measure of t 900 . Figure 12. Comparison of the cumulative distribution of IGM+CGM transmission evaluated for zs = 3.1 (orange), zs = 3.3 (magenta), and zs = 3.5 (green) using the t 900 bandpass (solid curves) vs. t F336W or t LBC−U (shortdashed curves), where F336W is the WFC3-UVIS filter and LBC-U is the U filter described by Grazian et al. (2016). Note that the largest differences appear at relatively high values of transmission.
band-averaged transmission t 336 = 0.045, ∼ 8.3 times smaller (2.3 mag). Figure 12 shows cumulative distribution functions of t 900 , t 336 , and t LBC (the broadband U filter used by Grazian et al. 2016) for a large ensemble of sightlines to z s = 3.1 and z s = 3.5. At z = 3.1, the median transmission in the relevant LyC de-  uated: [880,910] is the spectroscopic band used in this work; the LBC-U filter is close to SDSS u , and is described by Grazian et al. (2016); HST-F336W is the HST/WFC3-UVIS filter F336W.
Similarly, compared to the ground-based U LBC used to measure LyC emission from spectroscopically identified galaxies at z ∼ 3.3 by Grazian et al. (2016Grazian et al. ( , 2017, the spectroscopic [880,910] Å bandpass is > ∼ 2 times more likely to include sightlines with t LyC > 0.2, and 9 times more likely to have sightlines with t LyC > 0.4. These differences could be quite large if there is limited dynamic range for LyC detection from individual sources -as has been the case for all surveys to date -potentially producing large differences in the fraction of sources with significant LyC detections even for surveys that nominally reach the same flux density limit in the LyC band. Care must also be exercised when comparing the results of surveys that use different LyC detection bands and/or IGM opacity models: Table 5 compares our determination of the mean IGM+CGM transmission evaluated for 3 different LyC detection passbands. The ensemble of sightlines used at each z s was identical. Note that we find that the mean IGM+CGM transmission evaluated using the U LBC filter for z s = 3.30 (z s = 3.40) is t LBC 0.149 (0.099), to be compared with the value assumed by Grazian et al. (2016Grazian et al. ( , 2017, t LBC = 0.28. The latter value was based on the IGM opacity models of Inoue et al. (2014) 22 . Additional ambiguities relevant to the quantitative comparison of LyC results arise due to differences in the definition of "escape fraction", discussed in more detail in §9.1.
The main point here is that the probability of detecting LyC emission-or of setting interesting limits on f 900 / f 1500 -depends very sensitively the source redshift z s , the bandwidth and relative wavelength/redshift sensitivity of the LyC detection band, and the fidelity of the correction for the IGM+CGM opacity.

INFERENCES FROM COMPOSITE SPECTRA OF KLCS GALAXIES
As discussed in §7, it is potentially misleading to interpret individual measurements of the quantity ( f 900 / f 1500 ) obs because of the large expected variation in t 900 from sightline to sightline. There are significant advantages associated with considering only ensembles of galaxies (sharing particular properties) that are large enough that the uncertainty in the ensemble average IGM correction is reduced. Since we have relatively high quality spectra of 124 objects remaining after cleaning the sample of potential contamination or obvious systematic issues, in this section we combine various subsets of the KLCS spectra to form high S/N spectroscopic composites. The resulting spectra are then used to obtain sensitive measurements of f 900 / f 1500 obs for which t 900 is welldetermined.
The individual KLCS spectra have a range of S/N owing to differences in apparent magnitude, redshift, and observing conditions. One might be tempted to combine them so as to maximize the S/N of the resulting composite, but the necessary weighting involved would invalidate the approach we adopt to correct for intergalactic opacity, which assumes that every sightline through the IGM is contributing equally to the net suppression of the intrinsic spectrum. Thus, we used the following approach to forming composites.
All extracted one-dimensional spectra (and their associated error spectra) of KLCS galaxies were shifted to the rest frame using the values of z sys (Table 2), where each was resampled to a common wavelength grid using spline interpolation and scaled according to its observed continuum flux density f 1500 . In forming a composite spectrum, the mean value at each restwavelength dispersion point was calculated after outliers with values differing by more than 3σ from the median were rejected. This algorithm is very effective in removing residuals from imperfectly subtracted night sky lines and other defects present in individual spectra. Figure 13 shows the spectroscopic stack of all 124 galaxies in the KLCS analysis sample. Fig. 13 also addresses the extent to which the IGM opacity models described in §7 are consistent with the attenuation actually observed in the composite galaxy spectrum. The observed spectrum (black) has been corrected using the calculated mean IGM+CGM transmission spectrum appropriate for a sample having z s 3.05 ( Figure 10), normalized so that f 1500 = 1. The orange spectrum is the similarly-normalized intrinsic spectrum of the stellar+nebular continuum for the stellar population synthesis (SPS) model found to be most successful in reproducing the observed composite spectrum (see §9.2). The cyan spectrum in Figure 13 is the model SPS spectrum after reddening using the best-fitting attenuation relation and color excess E(B-V) (discussed in §9.2 below) normalized so that f 1500 = 1. Note that the continua of the model spectrum (cyan) and the IGM+CGM-corrected observed spectrum (black) are in excellent agreement over the rest wavelength range 950-1210 Å-the range sensitive to the accuracy of the IGM Lyα forest Monte Carlo modeling vs. redshift ( §7 and Appendix B.) Figure 14 shows a zoomed-in version of Figure 13, illustrating the full KLCS sample before (top panel) and after (bottom panel) applying the IGM+CGM correction to the observed spectrum. The model spectrum shown in each panel is identical, and all spectra have been normalized to unity at rest-frame 1500 Å.

KLCS Subsamples
We formed a number of subsets of the final KLCS statistical sample based on empirical criteria that could be measured easily from photometry or spectroscopy of individual objects; Figure 13. The composite spectrum of all 124 galaxies included in the KLCS statistical sample (black). The stacked spectrum has been corrected for the mean IGM+CGM opacity appropriate for an ensemble with zs = 3.05 as shown in Figure 10. Some of the easily-identified spectral features are indicated, where the labels have been color-coded according to whether they are primarily interstellar (dark green), nebular (blue), stellar (red), or excited fine structure emission (purple). The orange (turquoise) spectrum is the best-fitting stellar population synthesis model (see §9.2) before (after) applying reddening according to Reddy et al. (2016a)  134 ± 0.007 · · · · · · · · · · · · All, not detected g 106 3.044 + 10.6 1.07 0.50 0.010 ± 0.003 · · · · · · · · · · · · z (Q1) 31 2.    Table 6.) The values of each, and the wavelength interval over which they have been evaluated, are indicated by magenta bars in both panels. The model spectrum (cyan histogram) is identical in both panels, and also the same as that shown in Figure 13 using the same color coding.
these include L uv , W λ (Lyα), and rest-UV continuum color (G − R) 0 ( §6.4). In view of the results of §7.1, a minimum subsample size N > ∼ 30 was maintained so that the uncertainty in the IGM+CGM correction is < ∼ 10% (see Table 4). Thus, for each of the aforementioned parameters we split the sample of 124 into 4 independent quartiles consisting of 31 galaxies each.
Additional subsets were formed according to the following criteria: objects with L uv ≥ L * uv (48% of the sample, very similar to a combination of the L uv (Q1) and L uv (Q2) subsamples) and those with L uv < L * uv (52% of the sample); according to whether Lyα appears in net emission (W λ (Lyα) > 0; 60% of the sample) or net absorption (W λ (Lyα) ≤ 0; 40% of the sample); and, finally, grouping together the galaxies exceeding the often-used threshold W λ (Lyα) > 20 Å for "Lyman Alpha Emitters" (LAEs). The KLCS LAE sub-sample (28 galaxies, or 22.6% of the total) is based upon the spectroscopicallymeasured W λ (Lyα), and happens to be nearly identical to the W λ (Lyα) (Q4) sub-sample of 31 galaxies. Table 6 includes values of parameters measured directly from the composite spectra or the mean value among the objects comprising the subsample: the number of galaxies in the subsample (N), the mean redshift of the objects in the subsample ( z s ), W λ (Lyα) (measured directly from the composite spectrum), the mean UV luminosity relative to L * uv ( L uv /L * uv ), and the measured flux density ratio f 900 / f 1500 obs . Values of t 900 for the transmission in the rest wavelength range [880,910] Å are given for both the "IGM only" and "IGM+CGM" Monte Carlo models ( §7 and Appendix B), where each value and its uncertainty were calculated by drawing ensembles of sightlines with the same number and z s distribution as the subsample, averaging the values of t 900 , and repeating 1000 times. The tabulated uncertainties reflect the 68% confidence interval for the mean transmission of the 1000 ensembles for each subsample. The composite spectra of several of the sub-samples listed in Table 6 are shown in Figure 15.
The corrected (emergent) flux density ratio f 900 / f 1500 out (equation 10) is also listed in Table 6, where the quoted errors include uncertainties in both the measurement and in t 900 . Table 6 also includes entries for the full KLCS sample ("All"), and "All" divided into subsets according to whether or not the individual LyC measurement had | f 900 /σ 900 | > 3 23 .
Finally, Table 6 includes entries for subsamples of KLCS 23 Note that we have not included entries that require a value of t 900 because the Monte Carlo models assume t 900 is independent of any property used to select the sample -clearly invalid in the case of a known detected or undetected sub-sample. formed according to z s , where z s (Q1) is the lowest-redshift quartile, and z s (Q4) the highest-redshift quartile. These are intended to show that there is no strong dependence of the results on source redshift once the redshift-dependent IGM (or IGM+CGM) corrections have been applied; the composite spectra of these two redshift subsamples are among those shown in Figure 15 On the basis of the results summarized in Table 6, there are two easily-measured empirical characteristics that correlate most strongly with a propensity to "leak" measurable LyC radiation: L uv and W λ (Lyα). Figure 16 illustrates, for selected sub-samples indicated on the figure, the dependence on L uv and W λ (Lyα) of the measured f 900 / f 1500 obs and inferred f 900 / f 1500 out for both IGM transmission models. The top panel of Figure 16 suggests an almost bimodal dependence of f 900 / f 1500 out on L uv , where the 2 lowestluminosity quartiles (i.e., the lower luminosity half) of the KLCS sample have f 900 / f 1500 out 0.13, whereas the UV-brighter half of the sample has f 900 / f 1500 out consistent with zero (3σ upper limit f 900 / f 1500 out < ∼ 0.02.) The transition luminosity below which f 900 / f 1500 out appears to increase from zero to 13% is very close to L * uv , which also happens to lie close to the median L uv of the KLCS sample.
The trend of f 900 / f 1500 obs and f 900 / f 1500 out with W λ (Lyα) -illustrated in the bottom panel of Figure 16 -is similar, albeit perhaps exhibiting a more gradual dependence of f 900 / f 1500 out on W λ (Lyα) compared to the relatively abrupt luminosity dependence of the subsamples grouped according to L uv . Figure 17 shows the KLCS sample color-coded according to quartiles in L uv (top) and W λ (Lyα) (bottom). There is clearly an inverse correlation between W λ (Lyα) and L uv in the sense that there is a dearth of galaxies with bright L uv and large W λ (Lyα) -e.g., the brighter two quartiles in L uv (Q1 and Q2) are dominated by galaxies with W λ (Lyα) < ∼ 10 Å (median 0 Å), and most of the galaxies in the two largest-W λ (Lyα) quartiles (Q3 and Q4) have L uv < L * uv 24 . Because of the well-established correlation between L uv and W λ (Lyα) within all LBG samples (including KLCS; see ,e.g., Shapley et al. 2003;Gronwall et al. 2007;Kornei et al. 2010;Stark et al. 2010), it is difficult to say with certainty which is the more reliable indicator of high f 900 / f 1500 out ; however, for reasons discussed in detail in §9.3 and §11.1 below, we believe that W λ (Lyα) is more directly linked to whether or not a galaxy has significant LyC leakage, and that the L uv dependence of f 900 / f 1500 out results from the UV luminosity dependence of the gas-phase covering fraction/column density along the line of sight.

GALAXY PROPERTIES AND LYC ESCAPE
While a determination of the total ionizing emissivity contributed by an ensemble of galaxies can be calculated directly from f 900 / f 1500 out , without knowing anything about the galaxies except for their FUV luminosity function, one would like to understand the physical causes of the large differences inferred among the sample. More generally, one would like to place LyC escape in the context of other galaxy properties, including the massive star populations and the ISM conditions that modulate LyC "leakage". Understanding which galaxies contribute -and why -to the metagalactic ionizing radiation field at z ∼ 3 will certainly improve our ability to make quantitative sense of the reionization era, where direct measurements of f 900 / f 1500 out will not be possible. These questions also bear considerably on our understanding of galaxies, the CGM, and the IGM at z ∼ 3 where the impact of LyC leakage is direct and highly relevant. In this section, we extend our analysis of the KLCS results to include properties of the galaxies, their stellar populations, and the radiative transfer of ionizing and non-ionizing FUV light through the ISM.
We begin by briefly reviewing the commonly-used forms of the LyC "escape fraction" ( f esc ) in the context of the KLCS dataset.

Escape Fraction: Definitions
The most common definition of escape fraction, often used in theoretical studies of re-ionization, is the fraction of ionizing photons produced by stars in a galaxy that escape into the IGM without being absorbed by H I within the galaxy (e.g., Wyithe & Cen 2007;Wise & Cen 2009). This quantity has also been called the "absolute escape fraction" ( f esc,abs ) to convey the fact that the ionizing photon budget against which the leaking LyC is compared includes all of the ionizing UV photons whether or not they are evident in observations. Thus, in order to determine f esc,abs , a measure of the intrinsic production rate of ionizing photons is necessary; unfortunately, direct empirical estimates of the production rate of ionizing photons (N ion ) are relatively impractical for high redshift galaxies. Instead, one typically uses SED modeling of young stellar populations (including dust extinction, population age, IMF, etc.) to estimate the intrinsic ionizing photon production. Doing this accurately requires knowledge of the intrinsic SED of the massive stellar populations in the extreme UV (EUV), especially for photon energies in the range 1-4 Ryd, as well as a detailed understanding of the distribution and composition of dust grains that attenuate and redden the intrinsic spectrum of the stars in the galaxy.
An alternative definition intended to be closer to the most readily-available measurements of high redshift galaxy populations is the "relative escape fraction" ( f esc,rel ; Steidel et al. 2001), is the attenuation in magnitudes as a function of rest wavelength and k λ parametrizes the attenuation relation. The term (L 900 /L 1500 ) int is the intrinsic ratio of the unattenuated stellar population of the galaxy, prior to transfer through the ISM, the CGM, and the IGM. Note that there are inconsistencies in the definition of f esc,rel in the literature having to do with whether or not dust attenuation is included; some authors have implicitly assumed that A λ (900) − A λ (1500) = 0, i.e. that the attenuation by dust affects ionizing and non-ionizing UV equally, or that dust affects only the non-ionizing UV [i.e., A λ (900 Å) = 0]. The definition of f esc,rel used by Grazian et al. (2017), expressed using the notation we have adopted here, is f esc,rel = ( f 900 / f 1500 ) out /(L 900 /L 1500 ) int . In order to avoid ambiguity, in what follows below we have attempted to be clear about any assumptions made in mapping f 900 / f 1500 out to more model-dependent quantities whenever relevant; §9.4 and 9.5 consider ISM models both with and without dust attenuation of emergent LyC flux.

Stellar Population Synthesis Models
The intrinsic energy distribution of the stellar sources, parametrized as (L 900 /L 1500 ) int , is a crucial ingredient to understanding the relationship between ionizing photon production and escape into the IGM. The intrinsic break at the Lyman limit in the integrated spectrum of stars depends on the age, metallicity, initial mass function (both slope and upper mass limit), and the effects of binary evolution of massive stars. Depending on the assumptions, SPS models predict spectra falling in the range 0.15 < ∼ (L 900 /L 1500 ) int < ∼ 0.75 (see, e.g., Leitherer et al. 1999;Eldridge et al. 2017). In the absence of external constraints, the factor of ∼ 5 uncertainty on (L 900 /L 1500 ) int translates directly into uncertainties on estimates of escape fraction.
Fortunately, high quality far-UV composite spectra such as those of the KLCS sub-samples significantly constrain the likely range of stellar population parameters consistent with the observations. Using similar far-UV composites constructed from the spectra of z 2.4 galaxies that were also observed in the near-IR as part of KBSS-MOSFIRE, Steidel et al. (2016) (S16) showed that, among the SPS models considered, the only set which could simultaneously match the details of the far-UV stellar spectrum -and correctly match the observed ratios of strong nebular emission lines in the spectra of the same galaxies -were those with low stellar metallicity, (Fe/H) 0.1 (Fe/H) , and which include the evolution of massive stars in binary systems (BPASSv2; Eldridge et al. 2017).
Following the procedure detailed in S16, we fit each of the composite KLCS spectra in Table 6 with a range of SPS models similar to that described by S16, with the following additions: we tested 3 different attenuation relations for reddening the SPS model spectra, in addition to varying the metallicity of the stars and the upper mass limit and slope of the IMF. We found that the newest publicly released BPASSv2.1 (Eldridge et al. 2017) models with stellar metallicity Z * = 0.001 ( 0.07 Z , assuming the solar abundance scale of Asplund et al. 2009), IMF slope α = −2.35, and upper stellar mass limit of 300 M (hereafter referred to as BPASSv2.1-300bin-z001), provided the lowest χ 2 (i.e., the best fit) for every KLCS composite listed in Table 6. However, the attenuation relation accompanying the best-fitting model varied among the sub-samples; Table 7 summarizes the parameters of the best fit for each sub-sample listed in Table 6.
Although three different attenuation relations were tested in the fitting, none of the composite spectra was best-fit by the Calzetti et al. (2000) attenuation relation. As indicated in Table 7, 13 of the 23 composite spectra were best fit by the attenuation relation constructed by combining the results of Reddy et al. (2015) with the far-UV extension of Reddy et al. (2016a) (R16), while the remaining 10 spectra favored the steeper, SMC attenuation curve which combines the empirical line-of-sight SMC extinction curve from Gordon et al. (2003) (assuming R V = 2.74) with an extension to the FUV derived as in Reddy et al. (2016a).
The values of (L 900 /L 1500 ) mod listed in Table 7 assume (for the time being) that the same attenuation curve and E(B − V ) value affects both the ionizing and non-ionizing UV light, and that the far-UV attenuation curves can be extrapolated to λ 0 ≤ 950 Å, the shortest wavelength over which they are empirically constrained (see Reddy et al. 2016a). Here we remind the reader that the values of f esc,rel listed in Table 7 are as defined in equation 11. Although the fits to the stellar continuum exclude all wavelength pixels shortward of 1070 Å, the fitted spectra appear to be excellent representations of the observed spectra all the way down to ∼ 930 Å, where the confluence of the stellar and interstellar Lyman series lines in the wavelength interval 910 − 930 Å depresses the galaxy spectrum relative to the model; an example is shown in Figures 13 (see also Figure 14.) The BPASSv2.1-300bin-z001 continuous star formation SPS model, including the contribution of the nebular continuum as described by S16, has (L 900 /L 1500 ) int 0.28 ± 0.03 25 , prior to reddening according to the values in Table 7; our assumption that the entire spectrum is reddened by the same attenuation relation further reduces the intrinsic L 900 /L 1500 by an amount that depends on E(B−V ) and the attenuation curve. Thus, the values of (L 900 /L 1500 ) mod in Table 7 correspond to the predicted emergent spectrum that would be observed if there were reddening by dust, but no photoelectric absorption by H I in the ISM along the line of sight. Table 7 lists the inferred values of f esc,rel , equal to the ratio between f 900 / f 1500 out (Table 6) and (L 900 /L 1500 ) mod (see equation 11), the flux density ratio for the best-fit (reddened) SPS spectrum. Figure 18 illustrates the associated values for each of the subsample composite spectra.
At this point, the value of f esc,rel has accounted for the relative attenuation between 900 and 1500 Å, but not for the absolute attenuation affecting f 1500 . In the context of the modeled UV SEDs summarized in Table 7, this factor is simply where k λ (1500) = 8.91, 13.05 for the R16 and SMC attenuation relations, respectively. We then define the absolute escape fraction as f esc,abs = f esc,rel /C1500 ; both C1500 and f esc,abs are included in Table 7. Figure 19 compares the values of this inferred quantity with the ob- 25 The quantity (L 900 /L 1500 ) int is closely related to the quantity ξ ion (see, e.g., Robertson et al. 2015), the number of H-ionizing photons produced per unit non-ionizing UV specific luminosity at rest-frame 1500 Å. The BPASSv2.1-300bin-z001 model used in the present case has ξ ion 25.5±0.1 for continuous star-formation ages of 7.5 < ∼ log(t/yr) < ∼ 8.7.  Table 6), with the values of (L 900 /L 1500 ) mod based on the best-fit continuum-reddened stellar population synthesis model. The ratio between the values of the red (squares) and filled green dots is defined as f esc,rel (Table 7). The error bars on the IGMonly points (open squares) have been suppressed for clarity (see Table 6 for values). The dashed horizontal line shows the value of (L 900 /L 1500 ) int for the SPS model. ( f 900 / f 1500 ) obs (see Table 6). As for Figure 18, the error bars have been suppressed for f esc,abs values for the "IGM-only" opacity model; the values of the points are collected in Table 7.
served f 900 / f 1500 obs for each KLCS subset composite spectrum. Note that all of the KLCS subsets with f esc,abs > ∼ 0.20whether based on UV color, UV luminosity, or W λ (Lyα) -are best-fit using the "line-of-sight" SMC extinction curve. This can probably be attributed to the relatively blue UV color of these composites, as the bluest/youngest star forming galaxies  (Reddy et al. , 2016a. b Flux density ratio for stellar population synthesis model, after reddening according to the specified extinction relation. The intrinsic value is (L900/L1500)int 0.28. c Inferred extinction correction at λ0 = 1500 Å for the best-fit SPS model. d Inferred bolometric luminosity from the attenuation-corrected L1500, in solar luminosities. e Ratio between ( f900/ f1500)out and (L900/L1500)mod, for the best-fitting SPS model. f Inferred absolute escape fraction, where fesc,abs ≡ fesc,rel/C1500 (see Fig. 19). g Values of fesc,rel and fesc,abs for the subsample with individual detections assuming that the IGM corrections are drawn from the top 12% of the t900 distribution function at z = 3.093 (see Table 6).
are known to be most consistent with a steep UV attenuation curve , rather than grayer extinction curves such as those of Calzetti et al. (2000) or Reddy et al. (2016a), perhaps due to the prevalence of small dust grains (or the absence of large ones). We note that fitting reddened SPS models to the observed spectra to estimate f esc,abs (Table 7 and Figure 19) rather than the more empirical f 900 / f 1500 out reveals that, in addition to luminosity and W λ (Lyα) dependence noted in §8 for f 900 / f 1500 out (each of which becomes steeper relative to f esc,abs ), accounting for continuum reddening corrections suggests that f esc,abs also depends strongly on UV color ((G − R) 0 ), in qualitative agreement with the analysis presented by Reddy et al. (2016b).
We return to a discussion of the implications of the various measures of escape fraction in §12. 9.3. "Spectral Morphology" and LyC Escape As described in §9.2 above, the relative escape fraction is f esc,rel 0.5 for the most-leaking KLCS subsets, with corresponding f esc,abs 0.3 − 0.4; both values are modeldependent in the sense that they are evaluated relative to assumed SPS models (whereas ( f 900 / f 1500 ) out is SPS modelindependent). We have shown that there is strong correlation between W λ (Lyα) and f esc , and it has been known for some time ) that W λ (Lyα) is anti-correlated with the strength of low-ionization (e.g., C II, Si II, O I) interstellar absorption features observed in the same spectra. The presence of strong lines of neutral and singly-ionized species in gas seen in absorption against the UV continuum from massive stars suggests significant neutral H column densities, N HI > 10 17 cm −2 , and the relatively constant W λ /λ for lines of the same species but varying λ f , where λ is the rest wavelength of the transition and f is the oscillator strength (e.g., Si II λ1190, 1193, 1260, 1304, 1526) indicates the presence of saturation 26 .
The observation that strong IS lines often do not become black at maximum optical depth then suggests partial covering of the continuum source (hereinafter, f c < 1) by the lowion-bearing gas. If the same partial covering were to apply to the H I, there would follow an obvious relationship between gas covering fraction f c and the probability that significant LyC could "escape" via the uncovered regions. Conversely, if a galaxy spectrum's low-ionization IS lines are black over a considerable range in velocity, one would expect no significant LyC photon leakage in our direction (e.g., Steidel et al. 2001;Pettini et al. 2002;Shapley et al. 2003Shapley et al. , 2006Quider et al. 2009;Steidel et al. 2010;Heckman et al. 2011;Jones et al. 2012). However, as recently shown (Vasei et al. 2016;Reddy et al. 2016b), f c < 1 for low-ionization metal lines is a necessary, but not sufficient, predictor of significant LyC leakage from a galaxy. The metal lines would be relatively insensitive to low-metallicity H I, and possibly also to metalenriched gas with LyC optical depths 1 < ∼ τ LyC < ∼ 10 where neutral and singly-ionized metallic species might not be the dominant ionization stages in the gas.
To examine trends in the residual depth and the velocity extent of interstellar absorption features, we divided each of the KLCS composites by its best-fit SPS model, and formed continuum-normalized spectra with the strongest stellar features removed. Figure 20 shows the line profiles of Si II λ1260 and λ1526, C II λ1334, and the Lyβ λ1025.7 and Ly λ937.8 for several of the KLCS subsets. Also shown in each panel of Figure 20 is the profile of Lyα emission, after subtracting the stellar continuum, dividing by a factor of 15, and adding 1 (both of the latter for display purposes). The spectral resolution (FWHM) for the metal lines (labeled "LRIS-R") and for the Lyman lines (labeled "LRIS-B") are shown in the bottom right of each panel. Figure 20 is arranged according to the measured f esc,abs , from largest (top left) to smallest (bottom right).
In addition to the trends already noted between W λ (Lyα) and f 900 / f 1500 out , there is a similar trend between the apparent depth of the Lyman series and low-ionization metallic absorption lines and f esc . In this context, it is noteworthy that the metallic lines shown in Figure 20 were recorded using LRIS-R, which have (as discussed in §3) spectral resolving power R ∼ 1400 (σ inst 90 km s −1 ), whereas the Ly series absorption lines were measured using LRIS-B, with R ∼ 800 (σ inst 160 km s −1 ). In spite of this, the Lyman series absorption lines in all cases have equal or larger apparent optical depth than the metal lines. The similarity of the apparent depth of Lyman lines (at least, for v sys < ∼ 0) in all of the composites shown suggests log(N HI /cm 2 ) > ∼ 18 -in which case all of the observed Ly lines would be strongly saturatedcoupled with partial covering of the far-UV continuum. In general, the maximum depths of both H I and metallic absorption lines occur at −300 < ∼ v < ∼ − 100 km s −1 , consistent with the hypothesis that the bulk of the cool, metal-enriched gas along the sightline "down the barrel" is out-flowing. The Ly line profiles exhibit more varied behavior for velocities v > 0, where absorption is less-deep overall, and where Lyβ absorption is significantly stronger than Ly , suggesting that both f c and log(N HI ) are lower for gas with v > 0 physically located on the near-side of the stars comprising the UV continuum.
The shapes and depths of the metal lines at v > 0 are difficult to interpret in the same way, since all of the metallic species shown in Fig. 20 are resonance lines that may manifest "emission filling" from scattered line photons re-emitted in our direction (see, e.g., Prochaska et al. 2011;Erb et al. 2012;Martin et al. 2012;Scarlata & Panagia 2015). However, all 3 of the lines illustrated have nearby excited fine structure lines which can act as "relief valves" for the resonance lines, so that a substantial fraction of absorbed resonance photons may be re-emitted in the non-resonance excited fine structure (EFS) lines -in this case, Si II * λ1264.9, C II * λ1335.7, S II * λ1533.4. Each panel of Figure 20 marks the position of the galaxy systemic redshift relative to the profiles of the EFS emission lines -with the same color coding used for the absorption profiles of the nearby resonance line. The implications of the EFS emission lines for LyC escape are discussed in §11.3.
In any case, the Ly series absorption lines between the Lyman limit (911.75 Å) and Ly-β provide the most direct constraints on N HI and/or the continuum covering fraction ( f c ) of H I near the galaxy systemic redshift along the line of sight. They are particularly useful for determining N HI when log (N HI /cm −2 ) < 18 -where the LyC would be optically thin regardless of f c ) -or when log (N HI /cm 2 ) > 20, where prominent damping wings of Lyα can be recognized even in spectra of modest spectral resolution. In between these limits (i.e., for 18 < ∼ log(N HI /cm −2 ) < ∼ 20), N HI is not well-constrained by the spectra, but the fact that the observable Lyman series lines are strongly saturated provides redundant constraints on the covering fraction of optically-thick H I. If these saturated Lyman series lines are resolved in velocity space, the residual intensity at maximum line depth relative to the continuum should be constant and equal to 1 − f c , where f c is the fraction of far-UV continuum covered by optically thick material. At the resolution of the KLCS galaxy spectra shortward of rest-frame Lyα (σ inst 160 km s −1 ) the cores of the higher Lyman series lines with σ v < ∼ 100 km s −1 would not be resolved, in which case any apparent residual flux in the cores of the Lyman series lines would correspond to a lower limit on the fraction f c of the stellar FUV continuum covered by optically thick H I.
In §9.2 above, we described fitting SPS spectra to the far-UV spectra of the observed KLCS composites in order to estimate f esc,rel and f esc,abs in a self-consistent manner. These spectral fits assumed that both ionizing and non-ionizing components of the integrated UV light from stars are reddened by the same dust screen with attenuation that is a smooth function of rest-wavelength governed by one of three model attenuation relations (Calzetti, R16, and revised SMC). However, introducing the parameter f c to describe the fraction of the FUV stellar continuum covered by optically thick H I requires a decision about how to handle dust reddening -i.e., whether the "covered" and "uncovered" portions of the stellar continuum have been reddened by the same dust opacity. It is easy to imagine that the dust column density could be lower along lines of sight that do not intersect any ISM gas Figure 20. Profiles of normalized intensity vs. systemic velocity for selected IS absorption lines for selected KLCS sub-samples after normalizing by the best-fit SPS model as described in §9.2: C II λ1334.5 (orange), Si II λ1260.4 (magenta), Si II λ1526.7 (blue), Lyβ λ1025.7 (yellow), Ly λ937.8 (pale green), and Lyα emission (red). The panels are arranged (left to right and top to bottom) in order of decreasing f esc,abs (last column of Table 7.) Lyα emission profiles are based on the SPS continuum-subtracted spectrum; each has been scaled down by a factor of 15 and shifted by +1.0 for display purposes. The three additional emission profiles at velocities > 0 km s −1 are excited fine structure lines of C II * λ1335.7 (orange), Si II * λ1264.8 (magenta), and Si II * λ1533.4 (blue). Short vertical dashed lines, color coded in the same way, mark the position of vsys = 0. The legend in the lower right of each panel shows the approximate spectral resolution for LRIS-R (relevant for the metal lines) and LRIS-B (relevant for the Lyman series lines).
with log(N HI /cm −2 ) > ∼ 18; if LyC escape requires geometric "holes" that have been cleared of all gas (neutral and ionized) between the stars and the IGM, then it may be that any residual LyC flux has not been attenuated or reddened by dust at all (a counter-example in a LyC-emitting galaxy has been discussed by Borthakur et al. (2014).) In sections 9.4 and 9.5, we consider two simple geometric models of the ISM intended to bracket the range of possibilities in extinction/reddening. In the first, which we call the "screen model", all light is attenuated and reddened by the same foreground screen-the standard assumption implicit in estimating dust extinction from far-UV photometry or spectroscopy, and the one used in §9.2 for calculations summarized in Table 7. The second alternative, which we call the "holes model", assumes that the covered portion of the FUV continuum is reddened by dust, but that the uncovered portion passes through the ISM without significant attenuation from dust or photoelectric absorption. We discuss the two cases separately below. 9.4. The "Screen" Model The emergent spectrum after passing through the ISM and CGM in the screen model context can be expressed as where A λ = k λ E(B − V ), S ν,int is the intrinsic stellar spectrum prior to attenuation and absorption in the ISM, and e −τ (λ) is the transmission spectrum through the ISM due to line and continuum absorption which depends on N HI and (to a lesser extent) the kinematics of the absorbing gas.
Using an approach similar to that described by Reddy et al. (2016b), we modeled the ISM 27 for each of the KLCS composite spectra by starting with the best-fit SPS model (Table 7) and fitting for additional H I line and continuum absorption shortward of Lyα. The model fits are identical to those described in §9.2, where the SPS model is reddened using the attenuation relation that provides the best fit to the CGM+IGMcorrected spectrum in the rest wavelength range 1070-1740 Å, after masking the regions listed in the top portion of Table 8. Simultaneously, we fit for the values of f c and log(N HI ) by comparing the predicted model spectrum with the observed spectrum, where the additional parameters are constrained in the wavelength intervals listed in the bottom part of Table 8. The regions near Lyα (and, to a lesser extent, near Lyβ) are sensitive to H I damping wings for log(N HI /cm 2 ) > ∼ 20.0, while the depths of Lyman series lines and the residual LyC emission constrain f c . Assuming the same SPS model and attenuation relations used in Table 7, we recorded the values of E(B − V), f c , and log(N HI /cm 2 ) that minimized the residuals with respect to the data; these are listed in Table 9. We estimated the uncertainties on each parameter by repeatedly perturbing every pixel of the observed spectrum by an amount consistent with the error spectrum (which includes contributions from both sample variance and shot noise) and refitting. The typical uncertainties on E(B − V), f c and log(N HI /cm 2 ) were found to be ±0.002, ±0.03, and ±0.10, respectively. Figure 21. Example fits using the "screen" model for the galaxy ISM, which includes H I absorption with column density log(N HI /cm −2 ) covering a fraction fc of the stellar UV continuum. The values of log(N HI /cm −2 ) and fc (see Table 9) are constrained by the damping wings of Lyα and Lyβ the depth of the higher order Lyman series lines, and the residual LyC flux in the range λ 0 = 880 − 910 Å; the model spectrum is otherwise identical to that used for the fits listed in Table 7. Each of the observed spectra (black histograms) has been corrected using the CGM+IGM model appropriate for the source redshift distribution of the galaxies in the sub-sample (thus, adjusted so that f 900 / f 1500 = f 900 / f 1500 out (see Table 7.) The best-fit value of f esc,abs computed directly from the models, and the RMS scatter of the best-fit values obtained from repeated fits of the perturbed observed spectra, are listed in the last column of Table 9. Example fits for the full KLCS galaxy sample ("All"), the highest L uv quartile (LQ1), and the subset with W λ (Lyα) > 20 Å (LAEs), are shown in Figure 21. As noted by Reddy et al. (2016b), f c and N HI generally have little covariance because one ( f c ) is determined by the minima at the cores of Ly series lines and the LyC region, while logN HI is constrained by the broad damping wings of Lyα. However, as discussed above, column densities in the range 18 < ∼ log(N HI /cm −2 ) < ∼ 20 are relatively poorly constrained because there are no discernible Lyα damping wings, yet the resolved Lyman series lines have saturated cores. This affects the fit for the "All, detected" sample, where f c and N HI for a fixed f 900 / f 1500 out are less well-determined 28 .
The tabulated f esc,abs values are generally very similar to those listed in Table 7; the difference in the present case is that the fits include constraints provided by the the depth of the Ly series lines and the Lyα damping wings in addition to the measurement of ( f 900 / f 1500 ) out . Thus, the values listed in Table 9 constitute more stringent checks on the internal consistency of the geometrical model of the ISM: rather than calculating f esc,abs from the SPS fit to the stellar continuum and ( f 900 / f 1500 ) obs alone, we now require consistency between the diminution of LyC flux relative to the SPS model and the other observational constraints on the column density and geometric covering fraction of the gas-phase ISM. In this sense, the fits summarized in Table 9 provide additional support for our simplified model for the galaxy ensembles. We return to a 28 Formally, the uncertainties are ±10% and ±0.75 dex, respectively. discussion of the implications in §10 below. 9.5. The "Holes" Model An alternative model for the effect of the galaxy ISM on the far-UV spectrum is one in which leaking LyC photons escape via completely transparent "holes" in the ISM, as discussed by Reddy et al. (2016b) (see also Zackrisson et al. 2013). This type of model is qualitatively supported by the fact that the equivalent widths of interstellar lines of low-ionization species is correlated with inferred reddening by dust, but the same is not true for higher ionization species like C IV (e.g., Shapley et al. 2003.) In the "Holes" model the spectrum is decomposed into two distinct components, one of which is reddened and absorbed by foreground dust and gas, and the other unattenuated and unabsorbed -and is therefore a scaled version of the intrinsic SPS spectrum. The parameter f c is the fraction of the intrinsic stellar spectrum covered by opticallythick absorbing material along our line of sight, and the net observed spectrum is given by the linear combination where A λ = k λ E(B − V ) cov and E(B − V ) cov is the reddening applied to the covered portion of the galaxy The second term on the righthand side of equation 16 also accounts for H I line and continuum opacity. For each KLCS composite spectrum, we fitted for the model values of f c , log(N HI /cm −2 ), and E(B− V ) cov that minimized χ 2 ; the best-fit model parameters are given in Table 10. Figure 22 shows the decomposition of the best fitting holes models for the same composites shown in Figure 21.
The main difference between the screen and holes models is the interpretation of f c : for the holes case, the differen- Figure 22. Fits to the same observed spectra as in Figure 21, but using "holes" models for the galaxy ISM. The model spectra in this case are linear combinations of a component which escapes the galaxy without attenuation (with relative amplitude 1 − fc; blue) and a component (orange) that is both reddened by dust and subject to line and continuum opacity due to H I; the attenuated component has an initial (unabsorbed) relative amplitude fc. For this type of model (see Table 10), since the reddened portion of the spectrum must pass through gas with large N HI (with τ LyC > ∼ 1000 for typical values of log N HI ), f esc,abs (1 − fc). tial attenuation between the covered and uncovered components requires accounting for extinction as part of the calculation of f c ; this is easily seen in the top panel of Figure 22, which shows that the relative contributions of the attenuated and unattenuated components to the observed (non-ionizing) spectrum are in the ratio ∼ 3 : 1, while geometrically the ratio is f c /(1 − f c ) ∼ 10. The concept f esc,rel is probably less useful in the context of the holes models, since f 900 / f 1500 out now has a more complex relationship to the intrinsic SED of the SPS spectrum than in the screen case (because of the assumption of two independent components, one affected by dust and the other unaffected.) However, in the holes models, f esc,abs is very straightforward to estimate provided τ (LyC) 1 for the covered component, in which case Equation 17 applies to all of the models listed in Table 10 except for the composite of individually-detected galaxies ("All, detected"). As discussed previously, the appropriate CGM+IGM correction for this sub-sample is much more uncertain than for any of the other composites; with the assumptions adopted earlier -that the detected subsample has IGM+CGM transmission drawn from the highest 12% -the best-fit model spectrum has 1 − f c ≈ 0.20 but f esc,abs has an additional contribution from ionizing photons that have passed through the H I layer without being completely attenuated (because the LyC optical depth is of order unity). 9.6. Model Assessment Models fits for both the screen and the holes models can simultaneously reproduce the overall spectral shape, the depth of the Lyman series absorption lines, the wings of the interstellar Lyα absorption feature, and the residual LyC flux using the same underlying SPS model. There is a slight tendency for both sets of ISM models to over-estimate the residual LyC flux for composites having the lowest values of f 900 / f 1500 out : e.g., the L uv > L * uv sub-sample has f 900 / f 1500 out = 0.006 ± 0.008 (Table 6) but the best-fit model gives f 900 / f 1500 out ∼ 0.03 (Tables 9 and 10). Such discrepancies would be expected if the covering fraction of the optically thick ISM layer is slightly underestimated by the cores of the Lyman series lines, or if there is H I with 17.5 < ∼ log(N HI /cm 2 ) < ∼ 20 that does not contribute significantly to the Lyα damping wings and whose Lyman series lines are not fully resolved in the spectra (and thus their depth is underestimated). Since the parameter f c is constrained in our modeling by simultaneously matching the LyC residuals (which should not depend on spectral resolution) and the Lyman series line depths (which may depend on resolution), larger reduced χ 2 evaluated over those spectral regions would flag inconsistencies.
Another possibility of course is that f 900 / f 1500 out is underestimated due to residual systematics in the background subtraction, in this case over-subtraction; this possibility cannot be ruled out entirely, but seems less likely than the presence of unresolved intermediate column density H I. Broadly speaking, the success of the simple models presented above demonstrates that there are real trends within the KLCS sample among the depth of the Lyman series lines, the residual LyC flux, and the column density of the dominant component of ISM N HI absorption; the trends -increased LyC leakage in galaxies with lower N HI , lower inferred f c , and larger W λ (Lyα) -are in the direction one might have antici-pated. We consider a simple physical interpretation of the results in §11.1, where we conclude that the "Holes" model does surprisingly well in explaining the observed interconnections between the observable quantities.

IMPLICATIONS OF ESCAPE FRACTION RESULTS
The term "escape fraction" has come to represent several different quantities, depending on the context and the underlying assumptions. We have gone through the exercise in the previous sections of calculating a number of different quantities related to LyC escape, justifying underlying assumptions by demonstrating that they do not lead to internal inconsistencies. Here we briefly digress to review the different quantities and the extent to which they are model-or assumption-dependent. As above, we begin with the most model-independent quantities, and proceed in order of increasing model dependence: 10.1. f 900 / f 1500 obs Assuming that the observations are free of systematics that would invalidate the results, this quantity is entirely modelindependent; however, it is of only limited use for individual objects because it is so strongly modulated by the variations of CGM+IGM opacity among different lines of sight ( §7). With ensembles in which more than 30 objects with similar redshifts and common observed properties are considered, the mean value of ( f 900 / f 1500 ) obs can be corrected using empirically-calibrated IGM+CGM models to obtain f 900 / f 1500 out . For the KLCS sample, 0 ≤ f 900 / f 1500 obs ≤ 0.14 among the galaxy subsamples (Table 6), with a sample mean of f 900 / f 1500 obs = 0.021 ± 0.002. 10.2. f 900 / f 1500 out This parameter is the most important one, in practical terms, for estimating the contribution of a galaxy population of known far-UV (e.g., at rest-frame 1500 Å) luminosity function to the metagalactic ionizing radiation field (see §11.4.) By our definition, this ratio characterizes the flux density of ionizing photons averaged over the 880-910 Å window, relative to the non-ionizing UV flux density near 1500 Å, statistically corrected for IGM+CGM opacity back to the "galaxy photosphere" at r gal = 50 kpc. As discussed in §7, the statistical uncertainty on the IGM correction is drastically reduced for an ensemble of galaxies provided the lines of sight are independent and the sample selection does not depend significantly on the IGM properties at fixed redshift. Remaining systematic uncertainties in the values of t 900 for an ensemble will stem primarily from the treatment of the CGM and possible systematic differences in the large-scale IGM environment. Table 6 includes values of t 900 under two different assumptions for the H I column density distribution for distances r > 50 pkpc from a galaxy; t 900 values with and without an enhancement in N HI from the CGM typically differ by ∼ 20%, which likely over-estimates the true level of systematic uncertainty, but would be a conservative upper limit. As Table 6 indicates, such systematics, if present, would affect the inferred values of f esc,rel and f esc,abs by a similar factor. 10.3. f esc,rel As discussed in §9.1 above, f esc,rel is a useful concept if the intrinsic (L 900 /L 1500 ) int of a galaxy's integrated stellar population is known, since the latter ratio is the absolute maximum possible ( f 900 / f 1500 ) out ; i.e., f esc,rel approaches unity only when there is no H I with appreciable LyC optical depth between the observed stars and r gal ∼ 50 pkpc.
If f esc,rel ≤ 1, (L 900 /L 1500 ) int 0.28 ( §9.2; Figure 23), and t 900 0.37 (Table 6), then the maximum value of f 900 / f 1500 obs expected for an ensemble of galaxies with z s 3.05 is f 900 / f 1500 obs,max 0.11. For an individual galaxy with z s = 3.05, t 900 < 0.67 (the ∼ 99th percentile IGM+CGM transmission; see Table B2), ( f 900 / f 1500 ) obs,max 0.30×0.67 0.20. Accounting for dust 29 or assuming a more typical IGM sightline would immediately reduce the expected upper limits on observed values. These small numbers would be reduced by large factors if the LyC regions were not sampled optimally ( §7.2), or the source redshift were significantly higher ( §7). A corollary to these considerations is that one should be suspicious of individual measurements of LyC leakage if ( f 900 / f 1500 ) obs > ∼ 0.2, and of ensemble measurements with f 900 / f 1500 obs > ∼ 0.11. Clearly, assuming a different SPS model, with a different (L 900 /L 1500 ) int , would also change the inferred f esc,rel given a measurement of, or limit on, ( f 900 / f 1500 ) out . Figure 23 shows the metallicity and age dependence of (L 900 /L 1500 ) int 29 For the attenuation relation that provides the best fit to the bluest galaxies in the KLCS sample (SMC), k λ (900) − k λ (1500) = 10.14; e.g., for E(B − V ) = 0.05, A(900) − A(1500) 0.51 mag, reducing ( f 900 / f 1500 ) obs,max to ∼ 0.13 and f esc,rel < ∼ 0.19.
(including self-consistently the contribution of nebular continuum emission to the non-ionizing UV light) for constant SFR BPASSv2.1-300bin models. In the context of such models, typical z ∼ 3 galaxies have SED-inferred ages of ∼ 50 − 500 Myr (7.7 < ∼ log(t/yr) < ∼ 8.7); the Z = 0.001 model that best reproduces the details of the observed spectra has 0.24 < ∼ (L 900 /L 1500 ) int < ∼ 0.33 over the same age range. Assuming younger ages for the SPS models would, naively, lower the required value of f esc,abs to match an observed ( f 900 / f 1500 ) out ; however, younger SPS models are somewhat bluer in the non-ionizing FUV, so that matching an observed galaxy SED requires more reddening by dust -and all reasonable dust attenuation relations have k λ (900) > k λ (1500), meaning that additional reddening reduces f 900 relative to f 1500 . For the screen models, fits of BPASSv2.1 models to the full KLCS galaxy composite with the same IMF and stellar metallicity (Z * = 0.001) but varying the age over the range log(t/yr) = [7.5 → 9 However, the best fits to the spectraboth in terms of reproducing the details of the stellar continuum and matching the observed LyC residuals, Lyman series line depths, and Lyα damping wings -are obtained for 7.9 < ∼ log(t/yr) < ∼ 8.5, where the parameters are consistent with those listed in Table 9 for the fiducial log(t/yr) = 8.0 model (see the bottom panel of Figure 23.) 10.4. f esc,abs In this section we have tried to emphasize that the parameter often seen as the "holy grail" of LyC studies -the absolute escape fraction f esc,abs -is several steps removed from an observational measurement, even at z ∼ 3 where direct measurements of residual LyC emission are feasible. Each of the steps is model-dependent, beginning with uncertainty in the shape and absolute intensity of the unattenuated EUV and FUV light from stars, followed by the net effect of dust on both the ionizing and non-ionizing UV, and, finally, the rather poorly-defined notion of what constitutes "escape" of an ionizing photon from the galaxy in which it was produced.
However, we have also attempted to leverage recent progress in understanding the ionizing sources through joint analysis of the stellar and nebular spectra (e.g., Steidel et al. 2016;Strom et al. 2017) and to follow a self-consistent thread capable of reconciling the detectability of LyC radiation leakage with all other available empirical constraints. With these caveats, we have measured f esc,abs for a sample of L * uv galaxies at z ∼ 3. Focusing on the "holes" models, for reasons discussed in more detail in the next section, the bottom line is that some galaxies appear to leak a substantial fraction of the ionizing radiation they produce-up to 30% for those in the highest quartile of W λ (Lyα), while others -notably the brightest galaxies, and those which have Lyα strongly in absorption -leak at most imperceptibly (< 3% at 2σ.) Averaged over the KLCS sample, all of which has L uv within a factor of a few of L * uv at z ∼ 3, f esc,abs ∼ 9% given our assumed SPS and reddening models (Table 10).
It is important to keep in mind the degree to which f esc,abs depends on the latter assumptions -if, for example, we had assumed that the attenuation relation is that of Calzetti et al. (2000) and that the appropriate SPS model is a solarmetallicity Starburst99 model see Steidel et al. 2016 for comparisons with the BPASSv2.0 models), then the best-fitting model parameters give f esc,abs = 0.25 30 , a factor of more than two times higher than our fiducial model (Table 7), but still under-predict the observed LyC flux by a factor of > ∼ 3. Changing the attenuation relation to R16 but keeping the SPS model fixed, f esc,abs = 0.33, with similarly poor χ 2 relative to the observed spectrum (Figure 23.) With the wrong choice of SPS model, it is not possible to simultaneously match the shape of the FUV spectrum, the Lyman series line depth, and the residual LyC emission (see also Reddy et al. 2016b.) 10.5. Comparison to Other Recent Results Given the many quantities commonly used to describe LyC measurements, it is useful to compare our results to others from the recent literature. The most directly comparable, in that the results involve composite spectral stacks, have been discussed by Marchi et al. 2017a,b, using spectra obtained as part of the VUDS Survey (Le Fèvre et al. 2015). Due to the minimum wavelength of 3800 Å the samples include only galaxies having 3.5 ≤ z s ≤ 4.3 (median z s = 3.8.) Marchi et al. (2017a) use a sample of 33 galaxies which were conservatively culled from an original sample of 46 using multi-band HST imaging to eliminate any significant possibility of foreground contamination. Their spectral stack, using a weighted average, has f 900 / f 1500 obs = 0.008 ± 0.004 at median redshift z = 3.80. If we apply our own determination of the mean IGM+CGM transmission at z 3.8, t 900 0.25 ( 10% lower than the transmission assumed by Marchi et al. 2017a) to facilitate a direct direct comparison with LLS, one obtains f 900 / f 1500 out 0.032 ± 0.016, somewhat smaller than that of the full KLCS sample ( f 900 / f 1500 out = 0057 ± 0.006; Table 6.). However, the Marchi et al. (2017a) sample is on average 1 magnitude brighter (in terms of M uv ) than KLCS; in addition, only 12 of 33 galaxies (36%) have Lyα appearing in emission, compared with 74/124 (60%) in the KLCS sample, in spite of the significantly higher mean redshift of the former. As discussed above, the KLCS results would predict that a brighter sample, dominated by galaxies with Lyα observed in absorption, would have relatively low f 900 / f 1500 out (see Table 6.) Marchi et al. (2017b) use a different subset of VUDS to investigate the dependence of f 900 / f 1500 out on properties believed to be correlated with LyC leakage based on existing detections at low and high redshift. They conclude that galaxies with large W λ (Lyα) and compact physical sizes (in the UV continuum and/or in Lyα emission) show evidence for significantly higher f 900 / f 1500 out than the remainder of the sample. Again adopting the measurements of f 900 / f 1500 obs and applying our estimated mean IGM+CGM transmission at z 3.8 for consistency, we find that the subsample with W λ (Lyα) > 50 Å has f 900 / f 1500 out 0.3, a factor of 2 higher than the KLCS subsamples with the largest f 900 / f 1500 out measured at z = 3.05. However, all subsets of the z s = 3.8 spectra imply f 900 / f 1500 out > ∼ 0.17, and the authors are careful not to claim to have measured uncontaminated values of f 900 / f 1500 out . In fact, if one adopts lowestmeasured f 900 / f 1500 out as the zero point, the differential between their highest and lowest f 900 / f 1500 out subsamples would yield f 900 / f 1500 out 0.16, entirely consistent with the maximum f 900 / f 1500 out among the KLCS sub-samples (Table 6.). Thus, we conclude that the results of both studies are consistent with those of KLCS.

A MODEL FOR LYMAN CONTINUUM ESCAPE
11.1. Trends with W λ (Lyα) In §8, we observed a strong dependence of the residual LyC emissivity f 900 / f 1500 out on spectroscopically-measured W λ (Lyα), and noted ( Figure 16) that the relation could be captured approximately by a linear relation The parametrization of equation 18 was chosen in part because the fiducial BPASSv2.1-300bin-t8.0 continuous star formation models used to fit the sub-sample composite spectra in §9.4 and §9.5 predict that W λ (Lyα) 110 Å, assuming "Case B" (ionization bounded) recombination and no selective attenuation of Lyα photons relative to the FUV continuum near Lyα 31 . The term inside square brackets is therefore the ratio of observed W λ (Lyα) to the maximum value expected under the most favorable conditions for Lyα production, i.e., where every LyC photon has been absorbed and converted to Lyα or nebular continuum in accordance with Case B expectations. An observation of the maximum W λ (Lyα) would also require that, within the spectroscopic aperture, the intensity of the Lyα emission line relative to the continuum flux density has not been significantly altered by resonant scattering of Lyα (i.e, spatial redistribution and/or selective absorption of Lyα photons by dust grains; see, e.g., Steidel et al. 2011). The lefthand panel of Figure 24 illustrates the nearly-linear relationship between 1 − f c (Tables 9 and 10 for the Screen and Holes models, respectively) and W λ (Lyα) among the composite spectra. This behavior strongly supports the hypothesis of a direct causal link between W λ (Lyα) in emission and the depth/strength of interstellar absorption features observed in the same galaxy spectra (e.g., Shapley et al. 2003;Steidel et al. 2010Steidel et al. , 2011Jones et al. 2012); it also lends credence to the use of simple geometric models of the galaxy ISM to understand the interplay between residual LyC emission, the depth of Lyman series absorption lines, and Lyα escape. The sense of the relationship is that the fraction of the Case B W λ (Lyα) (110 Å in our model) measured within a spectroscopic aperture is directly related to the portion of the stellar UV continuum that is not covered by optically-thick H I (1 − f c ) along the same line of sight.
Of course, one does not expect the production rate of Lyα photons to be compatible with Case B assumptions if there is significant detection of residual LyC flux, since any detected ionizing photons will not have produced the corresponding Lyα photons expected for an ionization-bounded geometry. For a typical value 1 − f c 0.1, the reduction in the Lyα source function would be small; however, in the limit of very high 1 − f c ∼ 1, an observer would see Lyα emission along the same line of sight only if the photons have been scattered in our direction after (initially) being emitted in another. In all likelihood this suggests that the correlation between (1 − f c ) or f esc,abs and W λ (Lyα) would turn over f esc,abs > ∼ 0.5. Con- 31 The predicted W λ (Lyα) varies by only ∼ ±10% over the age range 7.7 < ∼ log(t/yr) < ∼ 8.7 for the same models.  Tables 9 and 10) (1-fc) for each of the composite spectra, where fc is the inferred fraction of the stellar continuum covered by optically thick H I based on the joint constraints provided by the depth of the Lyman series absorption and the residual LyC emission ( f 900 / f 1500 out). (Right:) As for the lefthand panel, where (1-fc) has been replaced by f esc,abs under the Screen and Holes models. Note that, because of the way fc is defined for Holes models, (1-fc) is essentially identical to f esc,abs except for the "All, Detected" subsample.
The relationships shown in Figure 24 are surprisingly tight, given the complexity of Lyα radiative transfer compared to that of ionizing photons; for example, Lyα can be scattered by H I with N HI 10 17 cm −2 , and Lyα resonant scattering may increase the probability of absorption by dust or of scattering into a diffuse Lyα halo, in which case an aperture tailored to the apparent angular scale of the FUV continuum light (as for KLCS spectra) might miss a large fraction of the emergent Lyα flux (see, e.g., Steidel et al. 2011;Wisotzki et al. 2016). However, a "down the barrel" spectrum which includes a LyC detection will also capture the fraction of Lyα photons that pass through the same low optical depth holes in the ISM. The tightness of the relationships shown in Figure 24 suggests that most Lyα photons detected within the spectroscopic aperture may have passed through the same holes in the ISM. This scenario is easier to explain in the context of the "holes" models, as shown by the small scatter in both panels of Figure 24 evidently, applying dust attenuation to the entire EUV/FUV continuum makes it more difficult to account for the tightness of the correlation between f 900 / f 1500 out and W λ (Lyα). We also expect that using a larger aperture to measure the same galaxies would increase the scatter between W λ (Lyα) and f 900 / f 1500 out as an increasing fraction of the observed Lyα photons will have been scattered to large galactocentric radii via less direct channels through the ISM than the LyC photons.
In the discussion that follows, we continue to adopt the "holes" model as the simplest picture that remains consistent with the current observations.

Lyα Emission Line Kinematics
Theoretical predictions concerning the relationship between the kinematics and escape of Lyα emission and the leakage of significant LyC emission have varied depending on the physical model assumed for the structure and kinematics of the ISM (e.g., Yajima et al. 2014;Verhamme et al. 2015;Dijkstra et al. 2016;Verhamme et al. 2017). In general, most models predict that if the ISM is porous enough to allow the escape of ionizing photons along a particular line of sight, then one should also observe Lyα emission with a significant component near v sys 0 and/or a double-peaked profile with small velocity separation, which might indicate relatively low N HI along the line of sight. However, the picture that seems indicated by the KLCS results is that, while there is an overall correlation between the inferred value of N HI associated with the covered portion of the stellar continuum and f esc,abs (see, e.g., Table 10), it is the covering fraction ( f c ) that appears most closely linked with LyC leakage (as well as with W λ (Lyα)). Since the measured values of N HI are far too high to have allowed any detectable LyC leakage if the entire stellar UV continuum were covered, the apparent connection between N HI and f esc,abs may simply indicate that the incidence of clear channels through the ISM is higher when the characteristic N HI (with covering fraction f c ) is lower. Similar results have recently been obtained for a sample of LyC-detected low-redshift galaxies with N HI measurements Chisholm et al. 2018).
Interestingly, among the KLCS subsample composites, none exhibits a Lyα peak near v ∼ 0 -including the composite formed entirely of individual LyC detections -as shown in Figure 25. Indeed, the overall The overall shapes of the Lyα profiles -including the subsample comprised of individual LyC detections -are remarkably similar, with a fraction f (Lyα, blue) 0.31 ± 0.02 of the total Lyα emission emerging with v sys < 0. The center panel of Figure 25 shows that the velocity relative to systemic of the Lyα red peak (v peak ) decreases with increasing W λ (Lyα) , as found by Erb et al. (2014); Trainor et al. (2015) for similar-sized samples of LBGs and LAES at z ∼ 2 − 3. Nevertheless, the bottom panel of Figure 25) shows that the KLCS subsamples with f esc,abs > ∼ 0.1 tend to have a nearly constant v peak ∼ 250 km s −1 , whereas those with f esc,abs < ∼ 0.1 exhibit a wide range in v pk . The faint LAE sample of Trainor et al. (2015) (whose LyC properties are unknown, but which have Lyman β line depths similar to the KLCS LAE sample) continues a trend of slightly smaller Lyα velocity shifts for intrinsi- Figure 25. (Top:) Stellar continuum subtracted Lyα emission line intensity profiles (i.e., the area under the profile is proportional to W λ (Lyα)) for selected KLCS sub-samples. The profiles are color-coded according to the legend, which also includes the value of W λ (Lyα)/Å as given in Table 6. The spectrum labeled "faint LAEs" is the stacked profile of 300 continuum-faint Lyα-selected LAEs at z ∼ 2.7 from Trainor et al. (2015) with Muv ∼ −18 (i.e., ∼ 10 times fainter than the LAEs in the KLCS sample). (Center:) The systemic velocity of the red peak of Lyα emission, v peak , color-coded as in the legend for the top panel. Although W λ (Lyα) varies by a factor of ∼ 15 and v peak is clearly correlated with W λ (Lyα) , the fraction of the total Lyα emission with v < 0 is remarkably constant among all of the composite profiles: f (Lyα, blue) = 0.31 ± 0.02. (Bottom:) As in the center panel, but with the measured values of W λ (Lyα) replaced by the inferred f esc,abs . cally fainter galaxies at fixed W λ (Lyα) (center panel of Figure 25); this, along with the accompanying smaller maximum blueshifted velocity of the low-ionization metals, is attributed to lower outflow velocities. It is therefore quite possible that the characteristic v pk at fixed f esc,abs would also be smaller for galaxies with M uv > −19.5.
Clearly, Lyα peaking at v = 0 is not a necessary condition for significant LyC leakage. The non-zero vales of v pk among LyC-detected subsamples suggests that most Lyα photons escaping along the same line of sight as LyC photons have been scattered at least once from gas moving with v sys > 0. This suggests that the holes through the ISM allowing LyC photons to escape are not generally optically thin to Lyα photons without the benefit of a velocity kick to move them off resonance of gas near v sys = 0. 11.3. Resonantly-Scattered Metal Lines As briefly mentioned in §9.3, complementary insight into the geometric/kinematic distribution of optically thick material can be obtained from observations of the non-resonance excited fine structure emission associated with UV resonance absorption lines. Since the interstellar Si IIλ1260 and Si IIλ1526 absorption lines are resonance transitions, modulo the selective absorption of line photons by dust grains, there should be no net emission or absorption provided that the measurement aperture includes the entire scattering medium (e.g., Erb et al. 2012). For this reason, it is potentially interesting to measure the equivalent widths of the resonance absorption and the associated non-resonance EFS emission within the (fixed) spectroscopic apertures used for KLCS. Figure 26 compares the mean equivalent widths of the resonance lines Si IIλ1260 and Si IIλ1526 with the average equivalent width of the EFS emission lines Si II * λ1265 and Si II * 1533 as a function of W λ (Lyα) for the KLCS subsamples. The lefthand panel of Figure 26 shows the wellestablished anti-correlation between W λ (Lyα) and W λ (LIS) (e.g., Shapley et al. 2003;Du et al. 2018) where LIS refers to low-ionization resonance absorption lines (red filled points). In addition, the ratio of Si II * emission (skeletal points) to Si II absorption increases with W λ (Lyα) -there is a weak trend of increasing strength of W λ (Si II * with W λ (Lyα) (magenta dashed line in Figure 26) but a strongly decreasing trend of W λ (Si II) with increasing W λ (Lyα) causes the ratio W λ (Si II * )/W λ (Si II) to vary from ∼ 0.3 to ∼ 0.9 across the sample.
Although a detailed discussion of the implications of the EFS emission is deferred to future work, one interpretation of the trends illustrated in Figure 26 is that the physical extent of the scattering medium producing resonance absorption lines in "down the barrel" spectra is more compact for galaxies with higher LyC escape fractions. In such a scenario, the KLCS spectroscopic aperture (projected physical size 9.2 × 10.3 kpc at z s ) contains a significantly larger fraction of the gas with log N HI > ∼ 10 17 cm −2 in the CGM of galaxies having the highest LyC escape fraction. In other words -not surprisingly -galaxies are more likely to leak along an observed line of sight when they exhibit more spatially compact distributions of optically-thick gas in the interface between the ISM and CGM.
An interesting question is whether the more compact distribution of optically thick material is the "cause" or the "effect" of higher LyC escape. It is quite possible that the lowionization CGM is reduced by a more intense local ionizing radiation field in the inner CGM. That is, the reduced spatial extent of optically-thick low-ionization material surrounding LyC leaking galaxies might be caused at least in part by self-ionization. For example, we consider the W λ (Lyα) (Q4) subsample, which has νL ν (1500Å) 1.6 × 10 11 L and f 900 / f 1500 out 0.17 (Table 6), thus producing an average escaping ionizing luminosity νL ν (LyC) 10 44 ergs s −1 . The average specific intensity at the Lyman limit as a function of galactocentric distance r would be I ν (r) 10 −19 (r/50 kpc) −2 ergs s −1 cm −2 Hz −1 : at r = 50 kpc, the local ionizing radiation field intensity would exceed that of the metagalactic background by a factor of > 10, and could dominate the ionization state of the CGM to r > ∼ 150 kpc around individual galaxies.
11.4. Ionizing Emissivity of z ∼ 3 LBGs Quantitative measurement of LyC flux from the KLCS sample of galaxies at z ∼ 3 is the first step in addressing the issue of the connection between directly-observed galaxy properties and the propensity to leak significant LyC radiation. We have shown that even this step has required combining N ∼ 30 independent galaxy spectra according to various easily-measured empirical criteria in order to accurately correct the ensemble for the effects of intergalactic and circumgalactic gas along our line of sight.
If the KLCS galaxy sample were a representative subset of the population of star-forming galaxies that contributes to the metagalactic ionizing radiation field at z ∼ 3, then the values of f 900 / f 1500 out in Table 6 are all that is needed to convert a measurement of the far-UV (non-ionizing) emissivity uv = L uv Φ(L uv )dL uv , where Φ(L uv ) is the galaxy luminosity function evaluated in the rest-frame FUV (typically 1500-1700 Å), into the ionizing emissivity contributed by the same population of galaxies, LyC Luv,max Luv,min f 900 / f 1500 out × L uv Φ(L uv )dL uv (19) as discussed by many authors (see, e.g., Steidel et al. 2001;Shapley et al. 2006;Nestor et al. 2011Nestor et al. , 2013Mostardi et al. 2013;Grazian et al. 2017.) This method is especially useful because Φ(L uv ) has been measured over a substantial range in M uv at redshifts z 2 − 9 (e.g., Reddy & Steidel 2009;Oesch et al. 2010;Bouwens et al. 2015;Finkelstein et al. 2015;Parsa et al. 2016.) Note that f 900 / f 1500 out appears inside the integral in equation 19; as discussed in the previous section, although KLCS samples only a limited range in M uv (−19 , evidently only the UV-fainter half of the sample has significant f 900 / f 1500 out . We argued in §11.1 that a surprisingly tight relationship between the covering fraction f c and the measured W λ (Lyα) in the spectroscopic aperture suggests that the most direct signature of leaking LyC emission appears to be W λ (Lyα), and that the apparent dependence on L uv is secondary. However, it is straightforward to obtain approximate numbers for the contribution of KLCS-like galaxies to LyC using the empirical dependence of f 900 / f 1500 out on either L uv or W λ (Lyα).
If we assume for the moment that f 900 / f 1500 out = 0 for galaxies with M uv ≤ −21.0 and f 900 / f 1500 out = 0.11 for −21 < M uv ≤ 19.5 (Table 6 and Figure 16), then the galaxies that contribute significantly to LyC at z ∼ 3 have uv 1.21 × 10 26 ergs s −1 Hz −1 Mpc −3 and thus contribute LyC 0.11 × 1.21 × 10 26 13.0 × 10 24 ergs s −1 Hz −1 Mpc −3 . The contribution would be slightly lower, LyC 11 × 10 24 ergs s −1 Hz −1 Mpc −3 , assuming f 900 / f 1500 out = 0.094 as for the "IGM Only" transmission model. 11.4.2. Weighting by W λ (Lyα) Alternatively, we can use our derived empirical relation between W λ (Lyα) and f 900 / f 1500 out to estimate the contribution of KLCS-like galaxies to LyC . In this case, we take f 900 / f 1500 out outside the integral in equation 19 and use the distribution function n[W λ (Lyα)] of UV continuum-selected z ∼ 3 galaxies to estimate contributions to LyC relative to the integral over the UVLF, uv .
We begin with equation 18, the best linear relation between f 900 / f 1500 out and W λ (Lyα); where again the statistic f 900 / f 1500 out is adopted as the least model-dependent link between the FUV luminosity of galaxies and their contribution to the ionizing emissivity, since it does not depend on assumptions regarding SPS models, ISM radiative transfer, or attenuation by dust. However, for the present purpose, we assume that W λ (Lyα, CaseB) 110 Å as discussed in §11.1 above, and a model for mapping between f 900 / f 1500 out and W λ (Lyα) given by equation 18. For W λ (Lyα) > 0.5W λ (Lyα, CaseB) = 55 Å, which is unconstrained by the KLCS composite spectra, we assume that f 900 / f 1500 out reaches a maximum when (W λ (Lyα)/110 Å) ∼ 0.5, since at some point the increasing escape of LyC photons would result in a decreasing Lyα source function compared to Case B expectations. A naive expectation is that the source function of Lyα emission along a particular line of sight would be reduced by a factor (1-f c ) relative to Case B, followed by a further reduction proportional to f c due to scattering of Lyα photons by the covered portion of the UV continuum. In such a simplified picture, one might expect W ≡ W λ (Lyα)/W λ (Lyα, CaseB) ≈ 2f c (1 − f c ), which has a maximum W ≈ 0.5 when f c 0.5 and goes to zero as f c → 0 and as f c → 1. Solving for f c in terms of W gives Using equation 18 and the holes model fit (1 − f c ) ≈ 0.75W ( Figure 24) to connect f 900 / f 1500 out to 1 − f c over the range constrained by the data results in a model which maps the expected f 900 / f 1500 out to W over the full range of each: Two models mapping f 900 / f 1500 out to W λ (Lyα) are illustrated in the top panel of Figure 27; one is the best linear relation from equation 18, reflected about W = 0.5, the assumed maximum. The second is based on equation 22. Also shown in the top panel of Figure 27 is the relative incidence of W λ (Lyα) parametrized as an exponential of the form n(W λ ) ∝ exp(−W λ /W 0 ), with W 0 = 23.5 Å, very close to the observed distribution function in the KLCS sample, as well as that obtained from previous spectroscopic samples at z ∼ 3 Kornei et al. 2010.) Reddy et al. (2008) showed that the observed distribution of W λ (Lyα) Figure 27. (Top:) Estimated f 900 / f 1500 out as a function of W λ (Lyα). The red and blue solid curves correspond to the distribution functions of the same color coding in the legend. The dark green curve shows the relative incidence of W λ (Lyα) for a continuum-selected spectroscopic sample, which weights the incidence of the values of f 900 / f 1500 out expected in an ensemble. The dashed curves show the expected distribution of W λ (Lyα) weighted by the relative contribution to the total LyC emissivity, i.e., proportional to n[W λ ] × f 900 / f 1500 out. (Bottom:) The probability distribution function for f 900 / f 1500 out for an ensemble of LBGs, accounting for the distribution functions in the top panel for W λ (Lyα) > 0 and assuming that 40% of the ensemble has spectroscopic W λ (Lyα) < 0 (i.e., Lyα in net absorption.) The solid curves refer to the righthand axis, showing the cumulative fraction with lower f 900 / f 1500 out. The ensemble averages for the two distributions are indicated with the same color coding. at z ∼ 3 must be close to the intrinsic distribution for the galaxies contributing to the FUV continuum luminosity function. We further assume that 40% of continuum-selected galaxies have spectroscopic W λ (Lyα) ≤ 0 (and therefore have f 900 / f 1500 out 0). It is then straightforward to compute the probability distribution of f 900 / f 1500 out expected for an ensemble of continuum-selected galaxies with M uv < −19.5 (L uv > 0.25L * uv ) as in the KLCS sample. The resulting PDFs of the quantity f 900 / f 1500 out for an ensemble of galaxies with the assumed n[W λ (Lyα)] are shown in the bottom panel of Figure 27. Note that just over 50% of the PDF falls in lowest bin of f 900 / f 1500 out but the ensemble average is f 900 / f 1500 out 0.035.
Returning to the estimate of the contribution to LyC , the W λ (Lyα)-based f 900 / f 1500 out = 0.035 implies that galaxies with M uv < −19.5 contribute LyC 0.035 × ×1.71 × 10 26 5.9 × 10 24 ergs s −1 Hz −1 Mpc −3 . Because of the strong empirical connection between f 900 / f 1500 out and W λ (Lyα), and because previous work has shown that spectroscopic samples are unbiased with respect to W λ (Lyα), we are more confident of the W λ (Lyα)-based estimate of LyC than of the luminosityweighted estimate in §11.4.1.
For comparison, recent indirect estimates of the total ionizing emissivity at z ∼ 3 inferred from modeling the transfer of ionizing radiation through the IGM coupled with measurements of the mean transmissivity of the Lyman α forest (Becker & Bolton 2013) find LyC 10 × 10 24 ergs s −1 Hz −1 Mpc −3 with a systematic uncertainty of ∼ 0.4 dex (factor of 2.5). 32 . Thus, without extrapolating our results to lower luminosities that are not directly constrained by the KLCS sample, the ionizing emissivity LyC [L uv > 0.25L * uv ] 5.9 × 10 24 ergs s −1 Hz −1 Mpc −3 comprises ∼ 60% of the estimated total. In addition, it exceeds the total ionizing emissivity of QSOs at z ∼ 3, estimated to lie in the range 1.6 − 5 × 10 24 ergs s −1 Hz −1 Mpc −3 (Hopkins et al. 2007;Cowie et al. 2009.) A more detailed discussion of the implications of our result for the ionizing emissivity LyC and photoionization rate Γ at z ∼ 3, as well as a comprehensive comparison to other observational constraints, is deferred to a separate paper.

SUMMARY AND DISCUSSION
We have presented the results of a deep spectroscopic survey of z = 3.05 ± 0.18 Lyman-break-selected sources obtained using LRIS-B+R on the Keck 1 telescope. The spectra cover the observed wavelength range 3200 − 7200 Å with spectral resolving power of R = 800 − 1400, enabling direct measurement of the rest-frame LyC region 880 ≤ λ 0 /Å ≤ 910 for sources with redshifts z > 2.75, as well as the far-UV spectrum over the range 930 < ∼ λ 0 /Å < ∼ 1800. Very careful attention was paid to minimizing systematic errors associated with flat-fielding, background subtraction, and spectral extraction. In addition, the initial sample of 136 spectra was carefully examined to remove those contaminated by foreground (lower-redshift) galaxies or affected by other potential systematic errors. The final KLCS sample of 124 galaxies spans a UV luminosity range −19.5 > ∼ M uv > ∼ − 22.0 (0.25 < ∼ (L uv /L * uv ) < ∼ 3), and is representative of the larger spectroscopic parent samples from which it was drawn.
Measurements of the emergent LyC emissivity of galaxies are limited both by the dynamic range of individual spectra and -more importantly -by large and stochastic variations in IGM transmission for individual lines of sight averaged over the LyC measurement band [880,910] Å. Because of the importance of an accurate statistical model of the opacity of gas lying outside of the galaxy ISM, but between the source and the observer, to measurement of the emergent ratio of ionizing to non-ionizing flux density ratio f 900 / f 1500 out , we used the statistics of H I measured in both the general IGM and in 32 The systematic uncertainties are caused primarily by degeneracies among the spectral index of the ionizing sources and the effective mean free path of ionizing photons, both of which affect the conversion between a measurement of LyC and the corresponding photoionization rate Γ; see Becker & Bolton (2013) for detailed discussion. the CGM of galaxies from the Keck Baryonic Structure Survey (KBSS; Rudie et al. 2013) to produce new Monte Carlo simulations of the IGM transmission as a function of source redshift. These simulations were used to model the transmission for individual sightlines at a given source redshift, and to predict the average transmission for ensembles of sources with the same redshift distribution as the real KLCS sample.
We focused our analysis on composite spectra formed from subsets sharing common attributes that are easily measured from individual spectra. The ensemble spectra have much higher S/N than individual measurements, thus substantially increasing the dynamic range for LyC detection; in addition, the sub-sample sizes have been chosen so that the marginalized IGM+CGM transmission expected for the ensemble is known accurately. The mean transmission appropriate to each ensemble was used to correct the observed flux density ratio f 900 / f 1500 obs to the more relevant emergent flux density ratio f 900 / f 1500 out without adding significant additional uncertainty.
Our main results may be summarized as follows: • Averaged over the full KLCS sample, we find f 900 / f 1500 out = 0.057 ± 0.006. This value is independent of the nature of the ionizing sources within the galaxies, and is model-independent except for small residual uncertainties in the IGM+CGM transmission; the same value, within statistical uncertainties, is found for galaxies in the lowest and highest quartile in redshift. However, we find two other galaxy characteristics that correlate strongly with the propensity to leak measurable LyC radiation: rest-frame UV luminosity (L uv ) and rest-frame Lymanα equivalent width [W λ (Lyα)]. In particular, galaxies belonging to the lowest quartile ( L uv /L * uv 0.5) have f 900 / f 1500 out = 0.138 ± 0.024, while the highest quartile ( L uv /L * uv 1.7) have f 900 / f 1500 out = 0.000 ± 0.003). Similarly, galaxies in the highest quartile of Lyα equivalent width ( W λ (Lyα) = 42 Å) have f 900 / f 1500 out = 0.166 ± 0.025, while the lowest quartile ( W λ (Lyα) = −18 Å) have f 900 / f 1500 out = 0.013 ± 0.011.
• We show that trends in LyC emission are strongly linked to the apparent depth, relative to the continuum, of low-ionization metal lines and, especially, the depth of Lyman series absorption lines (other than Lyα), which always indicate a higher covering fraction than the metal lines ( §9.3.) In most cases, there is evidence for strong saturation among these ISM absorption features, suggesting that the apparent depth of the features is directly related to the continuum covering fraction of neutral gas with appreciable N HI . The Lyα emission line equivalent width W λ (Lyα), long known to be inversely correlated (and probably causally connected) to the low-ion absorption line strength (e.g., Shapley et al. 2003;Steidel et al. 2010Steidel et al. , 2011) is a manifestation of the effect of the same gas on the emergent EUV/FUV spectrum.
• We use the results above to explore more modeldependent characterizations of the LyC emission from the galaxies, including the relative and absolute escape fraction of ionizing photons from the parent galaxies, f esc,rel , f esc,abs respectively. We show ( §9.2) that the non-ionizing (i.e., 920-1750 Å) composite FUV spectra can be very well matched by SPS models with low stellar abundance ([Fe/H]∼ −1.15; Z * /Z 0.07) which include binary evolution of massive stars reddened by plausible far-UV attenuation relations. These models have intrinsic L 900 /L 1500 0.28 ± 0.03 over the plausible range of star formation age 7.5 < ∼ log(t/yr) < ∼ 8.7, corresponding to ξ ion 25.5 ± 0.1.
• We show ( §9.4 and §9.5) that the same SPS and attenuation models can be used with simple geometrical models of the ISM to simultaneously account for the depth of the Lyman series absorption lines, the residual LyC emission, the damping wings of the Lyα absorption profile (which constrains N HI in the dominant absorbing gas) as well as the stellar FUV spectrum. The best-fit models provide constraints on N HI , the continuum covering fraction f c , and the intrinsic and escaping fraction of both ionizing and non-ionizing stellar photons.
• We examine in more detail ( §11.1) the tight relationship between W λ (Lyα) and the inferred uncovered portion (1 − f c ) measured from model fits to the composite spectra. For our favored model geometry, in which both LyC and a fraction of Lyα photons escape along the line of sight through the same "holes" in the ISM, the average LyC escape fraction for KLCS galaxies at z ∼ 3 is f esc,abs = 0.09 ± 0.01. Within the sample, f esc,abs 1 − f c varies linearly with the observed W λ (Lyα) ( §11.1), over the range 0 < ∼ f esc,abs < ∼ 0.29 and 0 ≤ W λ (Lyα)/Å < ∼ 50. We argue that the kinematics of Lyα emission are not uniquely connected with the propensity to leak LyC radiation into the IGM; rather, it is the presence of holes (parametrized by 1 − f c ) in the scattering/absorbing medium along our line of sight which apparently controls both W λ (Lyα) and f 900 / f 1500 out .
• Finally, we use the observed interrelation of L uv , f c , W λ (Lyα), and f 900 / f 1500 out to estimate the contribution of galaxies with M uv < −19.5 to the total ionizing emissivity LyC at z ∼ 3 ( §11.4). Because the correlation between the spectroscopically-measured W λ (Lyα) and inferred f esc,abs is both strong and linear, our most confident estimate comes from integrating the non-ionizing UVLF weighted by the product n[W λ (Lyα)] × f 900 / f 1500 out , where n[W λ (Lyα)] is the relative incidence of W λ (Lyα) and f 900 / f 1500 out is the empirically-measured ionizing to non-ionizing continuum ratio evaluated as a function of W λ (Lyα). We find a mean effective f 900 / f 1500 out ∼ 0.035 for galaxies with −19.5 ≥ M uv ≥ −22.1, accounting for > ∼ 50% of the total LyC estimated using indirect methods, and exceeding by a factor of ∼ 1.2 − 3.7 the total ionizing emissivity of QSOs at z ∼ 3.
Throughout this paper, we have deliberately concentrated on the least model-dependent, most easily-measurable parameters -UV luminosity and W λ (Lyα) -in seeking significant trends between galaxy properties and the propensity to leak significant LyC flux. There has recently been a great deal of interest in identifying observational signatures that strongly correlate with LyC escape but are more amenable to measurement at very high redshift where direct measurements of LyC flux become difficult (e.g., Steidel et al. 2011;Heckman et al. 2011;Zackrisson et al. 2013;Nakajima & Ouchi 2014;Erb et al. 2014;Trainor et al. 2015Trainor et al. , 2016Verhamme et al. 2017). Among the most promising indirect indicators discussed are anomalous nebular line ratios, possibly due to density-bounded H II regions, where nebular emission from lower-ionization species (such as [OII] or [SII]) may be suppressed. One proposed manifestation is a large value of the ratio O32 ≡ I([OIII]4960 + 5008)/I([OII]3727 + 3729), which has the advantage of being easily-observable by JWST in the redshift range z 6 − 9. The subset of the KLCS sample with z > 2.92, comprising 94 of 124 galaxies in the sample (76%), allows measurement of [O III] λλ4960, 5008 and [O II] λλ3727,3729 in the K and H atmospheric windows from the ground 33 . We are currently completing Keck/MOSFIRE observations of the relevant KLCS sub-sample to investigate correlations between these excitation/ionization-sensitive line ratios and direct LyC observations. The results will be presented in future work.
Deep ground-based near-IR photometry, Spitzer/IRAC coverage, and selected coverage with deep HST WFC3-IR and ACS imaging are also well underway -these observations, when completed, will provide stellar population parameters for all of KLCS. The HST observations will serve as a final verification that the sample has been cleaned of all contaminating sources. Additionally, in view of the apparently-strong anti-correlation with UV luminosity and UV reddening ( §8; Reddy et al. 2016b) and positive correlation with W λ (Lyα) of the LyC escape fraction, we are in the process of extending KLCS to fainter UV luminosities, bluer UV colors, and stronger W λ (Lyα) using Keck/LRIS-B+R.
Continued development of state-of-the-art stellar population synthesis models (e.g., Eldridge et al. 2017;Choi et al. 2017), particularly those focusing on the rest-frame UV, will clearly be important to modeling the sources of LyC emission in a fully-consistent treatment including the stars, gas, and dust in forming galaxies. "Rules of thumb" established at z 3 will be crucial to interpreting sources at much higher redshift where only indirect measurements are possible.

A. KLCS DATA REDUCTION DETAILS
A.1. Scattered Light and Flat-Fielding As discussed in §4, adequate illumination for spectroscopic flat-fields necessitates the use of twilight sky, observations. However, because of the rapidly increasing system sensitivity as wavelength increases from 3100 Å to 4000 Å, it is difficult in practice (particularly given the very short periods of evening and morning twilight on Mauna Kea) to obtain good exposures in the deep UV while controlling the amount of stray light reaching the LRIS-B detector.
Extensive tests performed during several LRIS observing runs established that stray light at a very low level was adversely affecting the LRIS-B observations, particularly at the shortest wavelengths where the sky background is faintest. We found that most of the problem could be attributed to stray light present in the twilight sky flat fields, which were invariably obtained just after sunset (or just before sunrise). Evidently scattering of sky light inside the telescope dome allowed low-level, undispersed stray light to reach the LRIS-B detector, which would then imprint diffuse undulations in the background level in flat-fielded science exposures. These subtle fluctuations had a significant effect on the precision of subsequent background subtraction.
Once the problem was identified, we adopted an improved strategy for obtaining twilight sky flat-fields: for each slitmask we obtained two sets of exposures, one of which had the LRIS-B U n = u filter (bandpass 3500/600 Å) inserted in the beam. The U n filter cuts off longward of ∼ 3900 Å and efficiently blocks the vast majority of stray light, since the LRIS-B filters are inserted just in front of the spectrograph camera. The presence of the filter made it possible to safely obtain flat field exposures with higher signal-to-noise ratio (S/N) shortward of 3900 Å. The U n twilight spectroscopic images were then combined with shorter twilight sky exposures taken without the filter in place. Since most of the stray light problem was confined to regions with λ < 3900 Å, the two sets of flats were "spliced" together to create composite twilight sky flats with high S/N throughout the 3200 Å−5000 Å wavelength range, with greatly reduced contamination by scattered light.

A.2. Image Rectification
After bias subtraction and CCD gain correction, the relevant "footprint" on the spectroscopic frames are identified by tracing the slit edges across the dispersion direction, fitting the edge positions with a 3rd order polynomial that is then used to extract each slit and transform it into a rectilinear image by applying nearest-integer pixel shifts in the spatial direction only. This produces a rectified slit image that requires no re-sampling of pixel intensities.
In the absence of differential atmospheric refraction, the previous step would make object traces run parallel to the CCD columns and thus simple linear traces could be used to extract one-dimensional (1D) spectra. For the data obtained prior to the ADC installation in July 2007, the wavelength dependence of objects within slits was accounted for, as described in §A.4.

A.3. Slit Illumination Correction
Imperfect milling or cleaning of LRIS slitmasks prior to insertion in the instrument can cause spatially-abrupt changes in slit illumination along the spatial direction, since (if everything is set up correctly) the slit plane is where both the telescope and the spectrograph camera are focused. This slit illumination function will also incorporate any intrinsic illumination gradients due to field-position-dependent vignetting within the telescope or instrument optics. Since our observations are attempting to detect signal at the level of < ∼ 1% of the night sky background, the typical 1-2% divots in the slit function need to be corrected before the sky subtraction stage. Left uncorrected, small systematic deviation from zero mean might remain in the background subtracted images; in short exposures this problem can easily go unnoticed, as other sources of random error are dominant. However, its systematic influence becomes clear in the average of many exposures.
Custom procedures were developed to correct for variations in slit function for both science and calibration exposures, accounting for small amounts of instrument flexure between each and the images used for illumination correction (i.e., shifts in the positions of the slit edges, and therefore of any slit defects, on the detector.) The slit function for each slit on a mask was obtained from rectified two-dimensional spectrograms ( §A.2) of the twilight sky by collapsing in the dispersion direction over the wavelength range deemed unaffected by residual scattered light, and averaging with iterative σ-rejection. The slit illumination correction was then applied to all rectified science and calibration exposures.

A.4. Atmospheric Dispersion Correction and
Object Tracing Since we are interested both in detecting LyC flux when present, and obtaining sensitive upper limits when it is not, it is advantageous to be confident about the location of an object's trace whether or not flux is evident. An often-used technique for tracing object positions in 2-D spectrograms is to determine spatial centroids regularly along the whole dispersion axis, which of course relies on the presence of signal over the same range. There is a tendency for traces extrapolated beyond the point where there is detectable flux to "wander" in search of upward fluctuations in the noise level. This is clearly undesirable in our particular application.
Prior to the deployment of the ADC on Keck 1, the trace of an object spectrum would differ from the trace of its slit edge by an amount determined by differential atmospheric refraction (DAR). Slitmask alignment for KLCS was always performed using direct imaging mode through the mask with a G-band filter in the LRIS-B beam, and guiding was maintained using an offset guide field CCD camera with peak sensitivity in V band. The exact value of the apparent spatial shift with time of source positions within slits depends on wavelength, airmass, slit position angle with respect to the elevation direction (the parallactic angle) and atmospheric conditions (temperature, pressure, water vapor pressure) at the time of observation. For reasons of practicality, KLCS observations were obtained (both before and after installation of the ADC) using a single mask and sky PA over a range of hour angle, albeit as close to the meridian as possible.
As a consequence, individual exposures of the same source could not be spatially registered at all wavelengths simultaneously. The difference in the apparent spatial position of the object over the 3200-5000 Å range of LRIS-B can be as large as ∼ 1 (∼ 7 − 8 LRIS-B pixels) if the slits have been oriented near the parallactic angle to minimize light losses due to DAR.
We developed an approach that deals with this problem in the earliest KLCS masks taken prior to August 2007. Using the atmospheric dispersion model of Filippenko (1982) we predict the value of the apparent spatial offset of the object position as a function of wavelength, and implement the shifts during the slit rectification step so that the object position in the rectified image is parallel to the slit edges, at a constant pixel location. Tests performed on observations of bright sources demonstrated that using the average atmospheric conditions for a given night on Mauna Kea provided sufficient precision in predicting the wavelength-dependent image shift from atmospheric dispersion. After August 2007, this procedure became unnecessary, as the ADC removes all but < ∼ 0.1 of differential shift over the full LRIS-B wavelength range.
A.5. Background Subtraction Once careful slit illumination and flat-field corrections have been applied ( §A.3 and §A.1), a low-order polynomial becomes an adequate description of the background intensity along the spatial direction. The accuracy of the background estimation depends to a large extent on the number of spatial pixels at each dispersion point that sample only "background" and not the outer parts of either the principal target or any part of nearby objects that happen to fall on the slit. Our experience has been that, at the level of sensitivity reached in the KLCS (m AB ∼ 28.0, or ∼ 2 × 10 −31 ergs s −1 cm −2 Hz −1 ,) the fraction of "clean" pixels remaining after source masking can be in the minority. We found it helpful, in terms of the number of spectra ultimately useful for LyC measurements, to enforce relatively long individual slits centered on high priority targets (slit lengths > 15 arcsec), sometimes at the expense of the number of targets that could be accommodated on a single slitmask.
We also deliberately lowered the priority for slit assignment of objects lying in obviously crowded regions on the deep images. In spite of this pre-screening of potential targets, the majority of the slits included serendipitous sources detected in the vicinity of the primary target. In fact, after combining ≥ 8 hours of LRIS spectroscopic observations the sensitivity is high enough to clearly detect the spectra of all sources seen in our deepest ground-based broad-band G or R images, which typically reach 1σ limits of G ∼ 29.5 and R ∼ 28.8 mag arcsec −2 ; . Using faint object positions predicted from the deep images we were in a position to exclude most or all pixels along the slit that might be contributing something other than sky background.
Background subtraction was performed in two stages: on the first pass we used object masks created from objects with significant detections in single 1800s LRIS-B exposures; in the second stage we compared the stack of the 2D background-subtracted spectrograms with broad band images and masked additional slit regions containing any source visible in either. An iterative fit with a polynomial function of order 2 to 4 was used to determine the final background model, which was subtracted separately for each 1800s exposure. The sky-subtracted spectrograms were then combined in 2-D, after shifting into spatial and spectral registration, by averaging, with bad pixel and outlier rejection.
The effectiveness of the above background subtraction procedures is discussed in more detail in §A.6.
A.5.1. Noise Model Throughout the reduction process, care was taken to preserve the original pixel sampling and number of photoelectrons per pixel in order to track formal uncertainties through all steps. A stack of all science frames made just prior to the 2-D background subtraction stage was made as a record of the number of sky+object photoelectrons recorded in every pixel. These images are used to estimate the variance each pixel [i, j] in the 2d spectrogram, where N exp is the number of 1800s exposures averaged to produce the value of pixel [i, j], C[i, j] is the number of photoelectrons recorded in the pixel, and ρ is the detector read noise, assumed to be ρ = 3.8e − pix −1 for LRIS-B.
A.6. Testing for Systematic Errors in 2D Spectrograms In this section, we outline tests performed to assess the quality of data reduced as in §4.2. Most importantly, we tested for the presence of residual systematic errors in the 2-D background subtraction by creating an average two-dimensional spectrum for an ensemble of final 2-D background-subtracted spectrograms for a set of KLCS targets which each led to > 3σ detections of residual flux in the rest-frame [880-910] Å. Fully processed, background-subtracted two-dimensional rectified spectrograms were wavelength calibrated and then put onto a common observed wavelength grid; We used the same object masks that had been created in the second pass of the sky subtraction phase for each slit ( §A.5), so that when the spectrograms were averaged, the masked pixels in each were excluded. The effective exposure time for the resulting stack [ Figure A1] is ∼ 90 hours, so that the random noise (due to photon statistics and detector read noise) is expected to be lower by a factor of 3 compared to any single observation in KLCS, thus more sensitive to possible systematic errors.
The rms noise values within the boxes shown in Figure A1 are higher in the 3450 − 3650 Å range due to lower system sensitivity (since the pixel intensity values are calibrated onto a flux density scale.) If the background subtraction is accurate, with systematics that are small compared to random errors, the distribution of pixel values is expected to well-described by a normal distribution with zero mean. The histograms to the right of the image show that this is the case for both wavelength regions. In addition, the measured rms noise level σ is entirely consistent with the noise model described in §A.5.1 (see below). Figure A1. (a) Two-dimensional composite spectrogram, combined using the object masks created during the background subtraction second pass, described in §A.5, for the 10 slits containing nominally detected LyC sources. Histograms of pixel fluxes, in units of 0.01 µJy (10 −31 ergs s −1 cm −2 Hz −1 ), are shown to the right for the regions indicated with colored boxes on the image: blue corresponds to the observed wavelength range 3450-3650 Å and red to the range 4000 − 4100 Å, and both windows are 80 spatial pixels ( 10.8 arcsec) wide. The values of the mean and standard error in the mean are shown for each histogram, along with the fitted value of σ assuming zero mean and a normal distribution, as expected for regions on the 2-D images without background subtraction systematics and with no residual flux from unmasked objects. (b) Similar to (a), except that in this case, the primary target galaxy was left unmasked, and the 2-D spectrograms were registered spatially to place the target at a constant Y-pixel location. The pixel statistics within the colored boxes (shown at right) should be sensitive to possible systematic sky subtraction errors related to the position of the primary target along the slit. As the histograms show [based on the same total number of pixels in each sample as in (a)], the residuals are consistent with Gaussian noise with zero mean, and σ pix values in the crucial 3450-3650 Å range are identical to within 0.3% to those in the same wavelength interval in (a).
As a second test for background subtraction systematics, we re-made the stacked composite spectrum using the same slits, but without masking the primary target galaxy (all other object masks remain the same). In this case, the spectrograms were spatially shifted before averaging so as to align the slit position of the primary target in the final image. The resulting two-dimensional stacked composite is presented in Figure A1(b); in this case, the stack would be sensitive to possible systematic errors that might be correlated with the position of the target object (generally the brightest object on the slit.) Unlike the previous case, which tested our ability to mask all significant sources in performing 2-D background subtraction, systematic errors in sky subtraction should become more prominent as more high-sensitivity spectrograms are averaged. To the right of Figure A1 are histograms of intensities drawn from regions on either side of the main target, beginning ±4 . 0 from the target centroid, evaluated in the same rest-wavelength intervals as in Figure A1(a) (with the same total number of pixels, divided into two separate spatial regions). The results are very similar -with zero mean and essentially identical σ within each wavelength region -to those obtained in Figure A1(a).
In evaluating the accuracy of our noise model, we examined each extracted 1-D spectrum after wavelength calibration, but prior to re-sampling in the dispersion direction or any other processing of the 1-D spectra and 1 σ error arrays. These spectra, examples of which are shown in Figure A2, remain in units of raw photoelectrons summed over the 10 (spatial) pixel extraction aperture of the final 2-d background-subtracted stack, so that the spectrum is in units of electrons per wavelength pixel (1 pix 2.14 Å) averaged over typically 18 1800s exposures. Each panel in Figure A2 shows the mean and standard error in the mean e − pix −1 evaluated over the rest-frame [880,910] interval, F 880−912 , the rms per pixel within that interval (σ 880−912 ), and the expected rms per pixel given by the background pixels and noise model described in §A.5.1 (σ sky ). The top two rows (i.e. the first 4 spectra shown) are formal non-detections of LyC flux, and the bottom row shows two examples of formally detected objects (Westphal/GWS-MM37 and Q1549-D3). These examples illustrate that the predicted and measured noise levels in the spectral regions used for LyC measurement are very close to one another. Note that the top-left panel shows the combined spectrum of Q0933-D20 after averaging the 2 independent 1-D spectra with inverse variance weighting.

B. IMPROVED IGM OPACITY MODELS FOR GALAXY SOURCES
For the purposes of generating Monte Carlo models of the IGM opacity, we adopted a combination of the measured N HI distribution functions f (N HI , X) IGM and f (N HI , X) CGM as shown in Figure B1, both as presented by Rudie et al. (2013) based on data from the Keck Baryonic Structure Survey (KBSS). Here f (N HI , X) is the distribution function, first introduced by Carswell et al. (1984), describing the incidence of H I absorption systems per unit co-moving path length X, as defined by Bahcall & Figure A2. Example extracted LRIS-B spectra, in units of mean photoelectrons per pixel prior to re-binning and flux calibration of the 1-D spectra. In each panel, the red curve is the 1σ noise level per pixel expected from the noise model. The vertical lines in each panel mark, moving left to right, rest wavelengths 880, 912, 1026 (Lyβ), and 1216 (Lyα). The top left panel shows one of the two galaxies observed on 2 separate masks, after combining the independently reduced spectra in 1-D, with inverse-variance weighting. The annotations above each panel show the object name, redshift, mean and error in the mean of the number of photoelectrons per pixel over the [880,912] rest wavelength interval (F 880−912 ), the rms per pixel over the same range (σ 880−912 ), and the noise level per pixel expected over the same pixels based on the noise model in §A.5.1. Peebles (1969). The parameter X is cosmology-dependent, but for spatially flat CDM models it is given by where H(z) = H 0 Ω Λ + Ω m (1 + z) 3 1/2 (B2) is the Hubble parameter. The CGM statistics (see Rudie et al. 2013) are derived based on measurements of f (N HI , X) for sightlines within 300 pkpc and ±700 km s −1 [∆z ≤ 0.0023(1 + z s )] of spectroscopically-identified KBSS galaxies in the redshift range 2.0 < ∼ z < ∼ 2.8. Aside from having somewhat lower redshift, the galaxies in KBSS were selected using very similar criteria to those of KLCS and occupy a nearly-identical range of stellar mass and star formation rate (see Adelberger et al. 2004;Steidel et al. 2004.) Figure B1 shows that within the CGM so-defined, the incidence of H I absorbers is ∼ 6 times higher than the IGM average at the high-N HI end of the f (N HI , X) distribution; for log N HI < ∼ 14 the CGM regions asymptotically approach the statistics of the general IGM. In general, f (N HI , X) is redshift-independent only if the physical properties of the absorbers do not evolve with redshift. Since the measurements in Rudie et al. (2013) are based on QSO sightlines with a path-length averaged redshift z = 2.37, we applied a redshift correction to the measured f (N HI , X) by expressing the number of absorbers of column density between N HI,min and N HI,max as N abs = NHI,max NHI,min z2 z1 N −β HI A(1 + z) γ dN HI dz where N abs is the number of absorbers with N HI,min < N HI < N HI,max expected over the redshift path interval [z1, z2], A is a constant chosen to match the observations at the effective redshift of the measured sample, β is the slope of f (N HI , X), and γ is a power Figure B1. Distribution functions f (N HI , X) inferred from the H I absorber catalog of Rudie et al. (2013). The best-fit power law slopes β for the "general IGM" (blue and black) and regions within 300 pkpc and δv < 700 km s −1 from a galaxy in the observed galaxy sample (such regions include ∼ 7% of the total effective survey volume; Rudie et al. 2012). The IGM has been divided into low-N HI and high-N HI subsamples at log (N HI /cm −2 ) = 15.2 with separate slopes and normalizations allowed (Table B1). Note that the incidence of absorbers with log(N HI /cm −2 ) > 16 is ∼ 6 times higher in the CGM of galaxies than in the general IGM (see text for discussion.) law exponent describing the redshift evolution of the number of absorbers per unit redshift, dN abs /dz ∝ (1 + z) γ . For reasons discussed in detail by Rudie et al. (2013), we divided the full range of N HI into two regimes: log N HI < 15.2 and log N HI > 15.2, with separate β, γ, and normalization A, as indicated in Figure B1. The values adopted for the MC simulations are summarized in Table B1 34 .
In the present case, our main interest is a quantitative estimate of the reduction in flux in a small wavelength range [880,910] Å in the rest-frame of a galaxy source. The effective flux decrement in this interval includes blanketing from Lyman series absorption lines -dominated by systems with log (N HI /cm −2 ) < 15.2) -and Lyman continuum opacity, dominated by higher N HI systems at redshifts close to z s . Most of the line blanketing component is contributed by Lyα and Lyβ falling at observed wavelengths near 900(1 + z s ) Å. For example, at z s = 3.10, the relevant redshifts are z(Lyα) 1.97 − 2.07 and z(Lyβ) 2.52 − 2.64. This implies that the measurements of Rudie et al. (2013) require no extrapolation in redshift to accurately estimate the line opacity contribution to the transmission in the wavelength interval [880,910] in the rest frame for a source at z s 3. Rudie et al. (2013) verified that the overall incidence of the low-column density lines is consistent with dN abs /dz ∝ (1 + z) 2.5 for the N HI range 12 < logN HI < 15.2, which has been assumed in our Monte Carlo realizations.
For Lyman continuum opacity, the most relevant absorption systems are those with 16 < log(N HI /cm −2 ) < 18 (Rudie et al. 2013) and redshifts within our chosen LyC interval of [880,910] in the source rest-frame, i.e., z s − ∆z ≤ z ≤ z s , where ∆z = 0.036(1 + z s ).  The LyC opacity depends on the higher N HI portion of the distribution measured by Rudie et al. (2013), including systems with log (N HI /cm −2 ) > 17.2 (Lyman limit systems), which are not well sampled in that data set. However, Rudie et al. (2013) showed that the power-law extrapolation of f (N HI , X) to log (N HI /cm −2 ) > 17.2 accurately reproduced the best available constraints on the incidence of Lyman limit systems (LLSs) at z ∼ 2.4, and, assuming γ = 1.0, at z ∼ 3 as well. The point representing the LLSs (at log (N HI /cm −2 ) = 17.7 in Figure B1 is the actual location of the data point for 17.2 ≤ logN HI ≤ 18.2 assuming the extrapolated slope and the observed total redshift path density of LLSs at z ∼ 3. For the MC models, LLSs were drawn from a distribution that assumes β = 1.473 over the full range 15.2 ≤ logN HI ≤ 21.0 35 . Monte Carlo realizations of IGM sightlines were made, using the same code as in Shapley et al. (2006) and Nestor et al. (2011) albeit with more detailed treatment of the higher-order Lyman series lines and with normalization parameters updated as in Table B1. The realizations include the full range 12 <logN HI < 21 drawn in a Poisson fashion according to the appropriate distribution function and redshift. Each realization for a given assumed source redshift z s is a synthetic spectrum covering the observed wavelength range 3160 − (1 + z s )1215.7 Å. A set of 10,000 realizations at each of 17 assumed source redshifts (2.70-3.50 in steps of 0.05) was performed using the "IGM-only" version of f (N HI , X). Table B2 summarizes the results in terms of percentiles in transmission t 900 , the mean unabsorbed fraction in the source rest-frame [880,910] interval. Also provided for each source redshift is the value of 1 − D B , where D B is the mean flux decrement from Lyman line blanketing in the rest-frame interval [920,1015] (Oke & Korycansky 1982). As shown in §7, 1 − D B is a close approximation to the maximum transmission expected in the [880,910] LyC interval at each value of z s , for sightlines without significant continuum opacity from high-N HI absorbers.
As in Rudie et al. (2013), a second set of Monte Carlo realizations was made by drawing absorbers from the same IGM distribution f (N HI , X) IGM , except within 700 km s −1 [∆z < 0.0023(1 + z s )]) of the source, where absorbers were instead drawn from the f (N HI , X) CGM distribution. This set of simulations is hereafter referred to as "IGM+CGM"; the statistics as a function of source redshift are summarized in the bottom half of Table B2.