neutral fraction decreases slowly from z ∼ 12 until z ∼ 6 (e.g. Finkelstein et al. 2019 ).

To draw conclusions about the evolution of the neutral fraction from observations of galaxies, one must attempt to understand the ionizing emissivity of galaxies as a function of cosmic time. This quantity is commonly parametrized as a function of three variables: the UV luminosity function ( ρ UV ), the ionizing photon production efficiency ( ξ ion ), and the fraction of ionizing luminosity that escapes the interstellar and circumgalactic medium (ISM and CGM) and proceed to ionize the IGM ( f esc ) (Robertson et al. 2015 ). While constraints are available for both ρ UV and ξ ion well into the epoch of reionization (Madau & Dickinson 2014 ;Stark et al. 2015Stark et al. , 2017 ) ), f esc is uniquely difficult to ascertain in the early Universe. Estimating f esc requires direct observations of the ionizing radiation from galaxies in the Lyman continuum (LyC) spectral region. The transmission through the general IGM of LyC photons escaping a galaxy depends sensitively on emission redshift and decreases rapidly beyond z ∼ 3.5 (Vanzella et al. 2012 ). This drop off is due to LyC absorption from trace amounts of H I and makes direct determinations of f esc impossible during the epoch of reionization itself.

MNRAS 521, 3247-3259 (2023) Models of reionization are distinguished by their assumptions about f esc . Finkelstein et al. ( 2019 ) present a 'democratic' model for reionization that assumes that the process is driven by faint sources with high f esc values. In contrast, the 'oligarchical' model of Naidu et al. ( 2020 ) concludes that massive luminous ( M UV < -18 and log(M * / M ) > 8) galaxies provide the bulk of ionizing photons during reionization. For testing assumptions of f esc during reionization, z ∼ 3-4 galaxies provide an essential laboratory, and can discern the fundamental properties that go v ern f esc at the highest redshifts these measurements can be made. Critically, these galaxies may be closer analogues to reionization era galaxies than those in the local Universe (e.g. Flury et al. 2022a , b ).

A number of LyC observational surv e ys at z ∼ 3-4 hav e attempted to measure average f esc values and potential correlations between f esc and the properties of galaxies. Success has been found by stacking deep observations of the LyC either photometrically (e.g. Be gle y et al. 2022 ) or spectroscopically (Marchi et al. 2017 ). Here, we focus on the Keck Lyman Spectroscopic (KLCS) surv e y, which included deep K eck/Lo w Resolution Imaging Spectrometer (LRIS; Oke et al. 1995 ;Steidel et al. 2004 ) spectra of Lyman break galaxies (LBGs) at z ∼ 3 (Steidel et al. 2018 ). Steidel et al. ( 2018 ) reported an average f esc = 0.09 ± 0.01 from 124 LBGs, estimated by stacking their rest-UV spectra with co v erage of the LyC region. After careful treatment of line-of-sight contamination using HST imaging, the average KLCS f esc was corrected to 0.06 ± 0.01 upon removal of four apparently contaminated galaxies from the sample (Pahl et al. 2021 ). This significant reduction in sample-averaged f esc highlights the importance of foreground decontamination in studies of LyC at high redshift. Galaxy properties correlated with inferred f esc in the KLCS are determined to be Ly α equi v alent width ( W λ (Ly α)) and UV luminosity ( L UV ), such that galaxies with stronger Ly α emission and lower L UV luminosities tend to have higher f esc (Steidel et al. 2018 ;Pahl et al. 2021 ). In the 2018 paper, we explained why W λ (Ly α) is more fundamental in its correlation with f esc , as W λ (Ly α) is modulated by the neutral gas co v ering fraction of a galaxy (Reddy et al. 2016b ;Gazagnes et al. 2020 ), which similarly modulates the escape of ionizing radiation. Trends between f esc , W λ (Ly α), and L UV have also been recovered in complementary LyC surveys at z ∼ 3 (Marchi et al. 2017(Marchi et al. , 2018 ; ;Be gle y et al. 2022 ), but L UV and W λ (Ly α) ultimately represent a limited parameter space from which to construct a comprehensive picture of LyC escape in star-forming galaxies.

Promising indirect indicators of f esc may surface from the feedback of star formation and its effect on the ISM and CGM of a galaxy. Cosmological zoom-in simulations coupled with radiative transfer calculations indicate that feedback from recent dense star-formation can induce fa v ourable channels in the ISM and CGM that allow ionizing photons to escape (Ma et al. 2020 ). Understanding f esc as a function of the surface density of star formation ( SFR ), stellar age, or specific star-formation rate would allow observational comparison to these simulations, and empirically connect the history of star formation to f esc . Additionally, dust attenuation is intricately linked to the neutral gas co v ering fraction in the ISM and CGM, but the relationship between f esc and E ( B -V ) at z ∼ 3 has thus far only been investigated using rest-UV observations (Reddy et al. 2016a , b ;Steidel et al. 2018 ). Thanks to multiband photometry available for the KLCS, which can constrain stellar population parameters, we can explore these relationships at high redshift, many for the first time.

In Pahl et al. ( 2022 ), we began by examining 35 galaxies from the KLCS that were co v ered by HST imaging, enabling measurements of rest-UV sizes. Together with SFR estimates from fits to multiband photometry, we measured SFR and attempted to constrain f esc versus SFR . We ultimately determined that the limited KLCS subsample with HST imaging was too small and unrepresentative to determine trends with f esc and galaxy property. In the present work, we extend the analysis of Pahl et al. ( 2022 ), by examining SED-modelled measurements of stellar mass (M * ), E ( B -V ), stellar age, and SFR instead, which allow nearly the entire KLCS sample to be utilized. By performing stacking of rest-UV spectra as a function of galaxy property, we investigate the dependence of f esc on these galaxy properties. Significant correlations will test existing reionization models and strongly inform future ones.

We organize the paper as follows: in Section 2 , we review the spectroscopic observations of the KLCS sample and its ancillary photometric measurements, and provide an overview of the SED and spectral fitting methodology. In Section 3 , we present the SED-modelled parameters for individual galaxies and estimates of ionizing escape from binned subsamples. In Section 4 , we explore similar observational analyses from the literature, connections to cosmological zoom-in simulations, and implications for reionization. We summarize our main conclusions in Section 4.3 .

Throughout this paper, we adopt a standard CDM cosmology with m = 0.3, = 0.7 and H 0 = 70 km s -1 Mpc -1 . The f esc values reported in this paper are absolute escape fractions, equi v alent to f esc, abs in Steidel et al. ( 2018 ), and defined as the fraction of all H-ionizing photons produced within a galaxy that escapes into the IGM. We also employ the AB magnitude system (Oke & Gunn 1983 ).

In order to understand how ionizing photon escape is tied to measurable characteristics of galaxy stellar populations or the spatial distribution of the interstellar gas, we require integrated photometric measurements that sample a wide wavelength baseline as well as direct constraints on the LyC emission. Both types of measurements are available for KLCS galaxies. In this section, we outline the data included in our analysis, featuring an o v erview of the KLCS sample, associated spectra, and the multiband photometry available for KLCS galaxies. We explain the methodology of SED fitting to determine galaxy properties and spectral modeling to estimate parametrizations of ionizing-photon escape such as f esc .

The primary goal of the KLCS was to examine the hydrogen-ionizing spectra of star-forming galaxies at z ∼ 3 (Steidel et al. 2018 ). To this end, 137 galaxies were observed with LRIS on the Keck I telescope on Mauna Kea, Hawai'i. Each object was observed for 8.2 h at minimum. These observations began in 2006 and were concluded in 2008. Of 137 targets, 13 galaxies were remo v ed due to either instrumental defects or spectroscopic evidence of contamination by foreground galaxies. The final sample presented in Steidel et al. ( 2018 ) numbered 124 galaxies. Of these, 15 galaxies apparently had significant flux density in the LyC spectral region, defined as having f 900 > 3 σ 900 , where f 900 is the average flux density between 880-910 Å in the rest frame, and σ 900 is the standard deviation of flux densities in the same spectral region. Objects meeting this criteria were defined as individual LyC 'detections,' with the remaining 109 galaxies labelled as LyC 'non-detections.' Despite the efforts of Steidel et al. ( 2018 ) to produce a clean sample of LyC leakers at z ∼ 3 by looking for evidence of spectral blending, foreground contamination remains a significant concern for surv e ys of LyC at high redshift (Vanzella et al. 2012 ;Mostardi et al. 2015 ). A galaxy along the line of sight to a z ∼ 3 source

Table 1. Photometric bands used in SED modelling.

Photometric bands

s , IRAC1, IRAC2, V 606 , J 125 , H 160 a Observ ed with Keck/LRIS. b Observ ed with F ourStar at the Magellan Baade 6.5m telescope. c Observed with the Multiobject Spectrometer for Infra-Red Exploration (MOSFIRE) on the Keck I telescope. d Observed with the COSMIC prime focus imager on the Palomar 5.08 m telescope (see Steidel et al. 2003 ). e Observed with the Prime Focus Imager on the William Herschel 4.2 m telescope (WHT) (see Steidel et al. 2003 ). f Observed with NIRC on the Keck I telescope (Shapley et al. 2001 ). g Observed with the Wide Field Infrared Camera (WIRC) on the Palomar 5.08m telescope. h From the Canada-France-Hawaii Telescope (CFHT) Legacy Survey. i Observed with CFHT/WIRCam as part of the WIRCam Deep Surv e y (Bielby et al. 2012 ).

can provide non-ionizing photons that masquerade as rest-frame LyC assuming a single redshift of z ∼ 3. As such, Pahl et al. ( 2021 ) presented new HST measurements of the 15 individual LyC detections in the KLCS sample, which were the objects most likely to be significantly contaminated by foreground light. These data were taken across five survey fields, including seven ACS/F606W ( V 606 ) pointings and 11 WFC3/F125W ( J 125 ), and WFC3/F160W ( H 160 ) pointings. Each pointing was observed for three orbits in each filter. Pahl et al. ( 2021 ) also utilized existing HST data for one object (Q1549-C25) from Mostardi et al. ( 2015 ) and Shapley et al. ( 2016 ). Based on the V 606 J 125 H 160 colours of the subcomponents in the high resolution, HST light profiles, two individual LyC detections were determined to be likely contaminated. An additional 24 LyC nondetections were included in the aforementioned HST pointings and were also analysed. Two of these were found to have associated subcomponents with colours consistent with foreground sources. In total, four galaxies were remo v ed from the KLCS, for a final sample size of 120, including 13 galaxies individually detected in LyC.

Several galaxy properties can be estimated from broad-band photometric measurements. In addition to the U n GR images used for original photometric selection of z ∼ 3 candidates (see Steidel et al. 2003 ), longer wavelength photometry of the KLCS has been obtained. We summarize the photometric information available for the objects in KLCS in Table 1 . A subset of these measurements are also summarized and analysed in Pahl et al. ( 2022 ). Specifically, optical, near-IR, and mid-IR data were available for the majority of galaxies in KLCS. We required at least one photometric measurement entirely redward of the Balmer break in order to accurately constrain the stellar populations. The filters that fulfilled this requirement for the KLCS were the H , K s , and Spitzer/IRAC bands. Thirteen objects do not have sufficient IR measurements (i.e. did not have any photometric measurements entirely redward of the Balmer break) and were remo v ed from our sample. In addition, we remo v ed two objects with significant scattered light in their ground-based light profiles from nearby objects, and one galaxy identified with multiple redshifts in the original KLCS spectrum.

HST photometry was included for the 35 objects observed in the HST V 606 J 125 H 160 pointings presented in Pahl et al. ( 2021 ). In addition to these 35 objects, 11 objects were co v ered by at least one HST filter, without the full V 606 J 125 H 160 data set required for contamination analysis. HST H 160 imaging was also available for four objects in the Q0100 field and two objects in Q1009 (Law et al. 2012 ). We remeasured integrated photometry for all objects in KLCS with V 606 J 125 H 160 HST data available largely following the methodology of Pahl et al. ( 2021Pahl et al. ( , 2022 ) ). In an effort to improve consistency between all photometric measurements, we adopted HST magnitudes from the same MAG AUTO parameter from SEXTRAC-TOR (Bertin & Arnouts 1996 ) that was employed by the groundbased measurements, rather than the isophotal HST magnitudes adopted by Pahl et al. ( 2021Pahl et al. ( , 2022 ) ). Finally, we note that HST H 140 measurements are available for three objects in Q0100 and five in Q1009, which were included in our analysis (Chen et al. 2021 ).

We attempted to correct the photometry from potential biases resulting from strong emission lines that lie in the bandpass of individual filters. Notably, we used existing Keck/MOSFIRE spectra with co v erage of [O II ] λλ3726, 3729, H β, and [O III ] λλ4959, 5007 rest-optical lines to correct broad-band H and K s flux measurements, depending on the wavelength of the observed line. We identified eight objects with neither Keck/MOSFIRE spectra nor additional photometry redward of the Balmer break aside from H or K s . We remo v ed these objects from the sample to ensure all galaxies in our analysis had at least one trustworthy photometric measurement redward of the Balmer break, free of potential emission-line bias. In addition, we used Ly α equi v alent widths presented in Steidel et al. ( 2018 ) to correct broadband G and V 606 flux measurements if the observed wavelength of the line was contained in the respective bandpass.

The final sample with sufficient multiband photometry for robust stellar population modelling consisted of 96 galaxies, of which 12 were individual LyC detections, which we define as the 'KLCS SED' sample. In Fig. 1 , we display the KLCS SED sample as a function of key observables from Steidel et al. ( 2018 ), including spectroscopic redshift z sys , Ly α equi v alent width W λ (Ly α), and UV luminosity ( L UV /L * UV , where the characteristic luminosity L * UV corresponds to M * UV = -21 . 0). We simultaneously present the characteristics of the parent KLCS sample of 120 galaxies. Median z sys , L UV , and W λ (Ly α) values for KLCS SED are consistent with those for the full KLCS sample.

In order to estimate stellar-population parameters such as stellar mass (M * ), star-formation rate (SFR), stellar age, and E ( B -V ) for the galaxies in the KLCS SED sample, we employed SED fits to the multiband photometry available for these objects. We broadly followed the fitting methodology of Reddy et al. ( 2022 ) (also see Pahl et al. 2022 ). In brief, we utilized BPASS stellar-population synthesis models (BPASS v2.2.1; Eldridge et al. 2017 ) assuming a Chabrier ( 2003 ) initial-mass function. We assumed a constant star-formation history (SFH) with stellar population ages greater than 50 Myr, such that stellar ages would not be less than typical dynamical-time-scales of z ∼ 3 star-forming galaxies (Reddy et al. 2012 ). We adopted constant SFHs as the y hav e been shown to reproduce independent measurements of SFR for galaxies z ≥ 1.5 (Reddy et al. 2012 ). Constant SFHs may also provide a better description of galaxies at the stellar masses of our sample ( ∼10 9 -10 10.5 M ), which may have less bursty SFHs than galaxies at lower masses (e.g. Dom ínguez et al. 2015 ). We adopted assumptions of metallicity of 0.14 times solar and an SMC dust attenuation curve (Gordon et al. 2003 ). We examined each SED fit individually for outlier photometric measurements, and dropped Spitzer/IRAC data with clear evidence of blending from nearby sources. Given our SED fitting methodology, we note that galaxies fit with larger masses tended to have higher SFRs and stellar ages. Over the mass range of our sample, we do not find a significant trend with M * and E ( B -V ), which may be expected considering the weak relationship that has been reco v ered between M * and E ( B -V ) when assuming the SMC dust extinction curve in the interpretation of galaxy colours vs. M * SED (McLure et al. 2018 ).

While f 900 can be measured for each object individually, constraining the LyC leaking in the vicinity of a galaxy requires an understanding of the attenuation of the signal from neutral hydrogen along the line of sight in the IGM and CGM. The transmission of LyC emission varies significantly from sightline to sightline at the redshifts of our sample, introducing large uncertainties on individual LyC measurements (Steidel et al. 2018 ). To circumvent this sightline to sightline variability, we used binned subsamples and composite spectra that reflect average effects of IGM and CGM attenuation on the LyC spectral region as in Steidel et al. ( 2018 ). In order to understand how ionizing-spectral properties vary with the properties produced by SED fits described in the previous section, we binned the KLCS SED sample as a function of M * , SFR, E ( B -V ), age, and specific starformation rate (sSFR; sSFR ≡ SFR/M * ). We created three bins for each property, each containing 32 galaxies, to ensure that the mean IGM + CGM transmission is known with 10 per cent uncertainty (Steidel et al. 2018 ) in subsequent composite spectra. We also created an 'all' sample, containing all 96 galaxies from KLCS SED, and binned subsamples for W λ (Ly α) and L UV .

For each subsample, we generated composite spectra representing the average spectral properties of the component galaxies. Following the methodology of Steidel et al. ( 2018 ) (also see Pahl et al. 2021Pahl et al. , 2022 ) ), each individual spectrum is first normalized to the average flux density in the non-ionizing UV spectral region, 1475-1525 Å in the rest frame. Using the set of normalized spectra for each binned sample, we then computed the sigma-clipped mean of the distribution of flux densities at each rest-frame wavelength increment, with clipping applied at 3 σ . We did not apply sigma clipping to the Ly α spectral region (1200 -1230 Å) in order to conserve the inferred composite Ly α profile. The error on the mean flux density at each w avelength w as propagated from the values of individual error spectra.

For each composite spectrum, we computed f 900 / f 1500 obs , which is the ratio between the average flux densities in the LyC region (880-910 Å, f 900 ) and the non-ionizing UV continuum (1475-1525 Å, f 1500 ). While this ratio is useful for discerning the average observed ionizing photon leakage relative to the non-ionizing ultraviolet luminosity density, as discussed abo v e, the spectra must be corrected for lowered transmission from the IGM in the LyC region in order to understand the average effect of LyC leakage has on its environment. We corrected the spectra using average 'IGM + CGM' transmission functions from Steidel et al. ( 2018 ), calculated at the mean redshift of each composite subsample, and based on the statistics of H I absorption systems along QSO sightlines presented by Rudie et al. ( 2012Rudie et al. ( , 2013 ) ). To demonstrate the characteristics of the composite spectra used in our analysis, we display the 'all' composite before and after the IGM + CGM transmission correction in the upper panel of Fig. 2 . Using corrected spectra, we repeated the measurement of the ratio of f 900 to f 1500 , defined as f 900 / f 1500 out , which applies to the Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 The 'all' composite alongside the same spectrum corrected from the average attenuation from the IGM and CGM at the mean redshift of the composite, z mean = 3.05. The uncorrected spectrum is shown with a thin, orange curve, while the corrected composite is shown with a thick, maroon curve. An inset is included to highlight the LyC spectral region. Bottom: IGM-and CGM-corrected composite spectrum alongside the best-fitting spectrum from the modelling process. The corrected composite is again shown with the thick, maroon curve. The best-fitting BPASS model is presented as a thin green curve. This model is summed from two component spectra, an attenuated portion displayed as a dashed pink line, and an unattenuated portion displayed as a dotted blue line, as per the 'holes' model of Steidel et al. ( 2018 ). An inset is included to highlight the LyC spectral region. The free parameters values of the fit that produced the model curves are f c = 0.96, E ( B -V ) = 0.093, and log (N HI /cm -2 ) = 20.85. ratio that would be observed at 50 proper kpc from galaxy centre (see Steidel et al. 2018 ).

While f 900 / f 1500 out is a useful empirical measurement of leaking LyC, f esc remains e xtensiv ely used in reionization modelling. In order to calculate the average f esc for each subsample, we require both an understanding of the intrinsic UV spectrum of the galaxies and the average effects from any intervening gas in the ISM. Thus, f esc is dependent on the assumed stellar population synthesis model, and we follow the well-motivated choices for such models discussed in Steidel et al. ( 2018 ). We introduced consistency between our multiwavelength and spectroscopic modeling by again using the BPASS stellar-population synthesis models of Eldridge et al. ( 2017 ). We coupled these models with an SMC e xtinction curv e (Gordon et al. 2003 ) and a range of E ( B -V ) from 0.0 to 0.6, and assumptions of metallicity of 0.07 times solar. This metallicity is similar to that assumed for the SED fitting and is consistent with the spectral modelling of Steidel et al. ( 2018 ) and Pahl et al. ( 2021 ). We model the ISM geometrically using the 'holes' approach, which assumes LyC light escapes through a patchy neutral-phase gas (Zackrisson, Inoue & Jensen 2013 ;Reddy et al. 2016bReddy et al. , 2022 ) ). The free parameters of the fit included the neutral gas co v ering fraction f c , the column density of neutral hydrogen N HI , and the dust attenuation from the foreground gas E ( B -V ) cov (i.e. the uncorrected portion is assumed to be dust free). In general, f esc is defined from f c , where f esc = 1f c . To demonstrate the fitting process, we display a fit to the corrected fullsample composite in the lower panel of Fig. 2 . Here, the modelled spectrum in green is split into an unattenuated (blue) and attenuated (pink) portion, representing the light that either escaped through clear sightlines in the ISM or was partially reduced by intervening material, respectively.

In order to estimate the uncertainty in average escape parameters for a given set of galaxies, we must understand the level of variability induced from sample construction, while simultaneously including the uncertainty on galaxy property measurements. We perturbed each Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 Figure 3. Galaxy property distributions of the KLCS SED sample. The five properties displayed here were inferred from SED fits performed for each galaxy. The median and standard deviations with respect to a given measurement are presented as large data points with capped error bars, while the typical (median) error on individual measurements are presented as smaller data points in the upper left with uncapped error bars. The edges of bins used for generation of composite spectra are shown as vertical dashed lines. The full sample was sorted according to each galaxy property and divided into three equal-sized bins (n = 32), which were then used to generate composite spectra.

individual measurement randomly by a Gaussian characterized by a width equal to its error, and constructed 100 modified parameter histograms for each galaxy property. We then binned each perturbed sample distribution into three subsamples of increasing galaxy property, with 32 galaxies in each subsample. Finally, with the goal of understanding how sample variance affects generated composite spectra, we used bootstrap resampling of the galaxies of each binned subsample. For each of the 100 sets of 32 galaxies (generated from the modified parameter distributions), we drew one new subsample with replacement. We subsequently created composite spectra and measured ionizing-photon escape for each random draw, using the process described earlier in this section. The mean and standard deviation of the f 900 / f 1500 obs , f 900 / f 1500 out , and f esc distributions generated from the 100 composite spectra were used as the fiducial value and error estimate for the corresponding binned sample. The errors determined from this bootstrap resampling were larger than those associated with measurement uncertainty, average IGM + CGM transmission variability, and errors from individual measurements.

Based on SED fits, we estimated M * , SFR, stellar age, E ( B -V ) and sSFR for each galaxy. We present the distribution of these SEDmodelled parameters for the KLCS SED sample in Fig. 3 . The median and standard deviation of each respective measurement distribution are displayed, respectively, as dark points and error bars. In the figure, we use dashed vertical lines to indicate the edges of the three equal-sized (n = 32) samples that comprise the bins for generating composite spectra.

A composite spectrum was generated for each binned sample detailed in Fig. 3 , and three estimates of ionizing-photon escape were measured, as described in Section 2.3 . The first, f 900 / f 1500 obs , is the ratio of ionizing to non-ionizing flux density directly observed in the composite generated from individual spectra. The second, f 900 / f 1500 out , is the same ratio instead measured from a composite corrected for mean line-of-sight attenuation from the IGM and CGM. Finally, f esc is a parameter estimated via stellar-population synthesis and ISM modeling of the full rest-UV composite. We display the

Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 (Steidel et al. 2018 ;Reddy et al. 2016bReddy et al. , 2022 ) ). Blue circles and purple stars are shifted left and right, respectively, for visual clarity.

three measurements of ionizing escape and their respective errors for each subsample binned as a function of galaxy property in Fig. 4 . These results are tabulated in Table 2 , which includes the median galaxy properties of each subsample.

To determine whether ionizing-photon escape is correlated with measured galaxy properties, we define a 'significant' correlation as fulfilling two criteria: the escape parameter varies monotonically across the three bins, and the difference between the escape parameter Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 in the highest and lowest bins was > 1 σ . Using f esc as an example, we define a > 1 σ difference as

where f esc, highest is measured from the third tertile of a given galaxy property, f esc, lowest is from the first tertile, and σ f esc , highest and σ f esc , lowest are their corresponding errors derived as described in Section 2.3 . As seen in the upper left-hand panel of Fig. 4 , both f 900 / f 1500 out and f esc are significantly, ne gativ ely correlated with M * in the KLCS SED sample, such that lower mass galaxies have higher escape fractions. The two measures of LyC escape are also significantly ne gativ ely correlated with SFR, shown in the upper central panel. While we find no significant correlation with f 900 / f 1500 out and E ( B -V ), we do find a significant ne gativ e correlation with f esc and E ( B -V ), as displayed in the middle left-hand panel. Discrepancies between these two parameters of LyC escape arise from the fact that the fitting process to determine f esc incorporates additional information from the composite spectrum, including the Lyman series absorption features and the UV spectral shape. Correlation between ionizing-escape parameters and E ( B -V ) is expected in our analysis considering that neutral gas and dust are spatially associated (Reddy et al. 2016b ;Du et al. 2018 ;Pahl et al. 2020 ). We also see an anticorrelation between f esc and E ( B -V ) when inferring E(B-V) from rest-UV spectral modeling. Finally, we find no significant correlation between f 900 / f 1500 out or f esc and stellar age or sSFR, seen in the middle left and center panels of Fig. 4 , respectively. We additionally computed these trends using a modified SED fitting process that assumed a Calzetti et al. ( 2000 ) dust extinction curve, and found qualitative agreement in the trends between f esc and galaxy properties, although we note that a stronger f esc -E ( B -V ) trend was reco v ered when using the SMC curve. We consider our SED modelling process that utilizes the SMC curve fiducial, considering greater consistency has been found between SED-fit SFR measurements and H α-based SFRs when assuming an SMC dust extinction curve rather than Calzetti et al. ( 2000 ) at z ∼ 2 (Reddy et al. 2022 ). Additionally, Reddy et al. ( 2018 ) demonstrated that z ∼ 2 star-forming galaxies have an IRX -β relation consistent with predictions only when assuming an SMC dust curve.

We also note that f esc and f 900 / f 1500 out are significantly correlated with W λ (Ly α) and L UV in this analysis, seen in the middle right-hand and bottom panels of Fig. 4 , mirroring the results of the full KLCS presented in Steidel et al. ( 2018 ) and Pahl et al. ( 2021 ). The positive trend between f esc and W λ (Ly α) has also been confirmed in additional z ∼ 3-4 LyC surv e ys (Marchi et al. 2017(Marchi et al. , 2018 ; ;Fletcher et al. 2019 ;Be gle y et al. 2022 ). In Pahl et al. ( 2022 ), we argued that reco v ering these well-established spectral trends is important for determining whether a sample is sufficiently large and representative for examining relationships between f esc and other galaxy properties. Considering the KLCS SED sample has both the size (n = 96) and dynamic range of galaxy properties to confidently reco v er trends between f esc and W λ (Ly α)/ L UV , we conclude that the KLCS SED sample is sufficient and representative, fulfilling the requirements for determining the trends between f esc and galaxy property presented in this section. We note that if subtle correlations do exist between f esc and age or sSFR, we may require a larger sample to discern these trends.

The connections between f esc and galaxy properties at z ∼ 3 provide key insights into the physics of ionizing-photon escape, and also indicate the most appropriate assumptions for f esc at even higher redshift, during the epoch of reionization. We find that galaxies with higher f esc tend to have lower E ( B -V ), which is consistent with a physical picture in which dust is spatially coincident with neutralphase gas in a galaxy, such that an ISM with a higher neutral-gas co v ering fraction will be both dustier and have lower associated f esc (Reddy et al. 2016b ;Du et al. 2018 ;Pahl et al. 2020 ). In addition, ionizing photons are more attenuated by dust than non-ionizing photons (Reddy et al. 2016a ). We find a ne gativ e trend between both f esc and f 900 / f 1500 out and M * , highlighting the fact that more massive galaxies at z ∼ 3 have conditions that are less conducive to LyC escape. This relationship is likely due to the fact that more massive galaxies tend to be dustier (e.g. Whitaker et al. 2017 ;McLure et al. 2018 ). Finally, the lack of trend between f esc and either stellar age or sSFR is in tension with the physical picture advanced in simulations that bursts of recent star formation induce favorable channels in the ISM and CGM for ionizing photons to escape (Ma et al. 2020 ). In this section, we introduce comparisons between our results and recent LyC surv e ys both at z ∼ 3 and in the local Universe. We also connect our ionizing-photon escape trends or lack thereof with radiative transfer modelling of simulated galaxies and the predictions from reionization models.

Direct comparisons can be made between our reported trends between f esc and galaxy property and those found in recent LyC surv e ys at z ∼ 3. Of particular note are the recent photometric LyC measurements of 148 galaxies from the VANDELS surv e y at 3.35 ≤ z ≤ 3.95 (Be gle y et al. 2022 ). These authors constrained the average f esc of the sample as f esc = 0.07 ± 0.02, consistent with f esc = 0.06 ± 0.01 measured from the uncontaminated KLCS (Pahl et al. 2021 ). The VANDELS LyC sample was binned in two as a function of a variety of galaxy properties. A positive correlation between f esc and W λ (Ly α) and a ne gativ e correlation between f esc and L UV reported in the VANDELS analysis aligns with the correlations found in the KLCS SED sample and the full KLCS (Steidel et al. 2018 ;Pahl et al. 2021 ). Best-fitting f esc values were also calculated for two bins of increasing M * in the VANDELS sample, which we display alongside our f esc versus M * measurements for the KLCS SED sample in Fig. 5 . A weak anticorrelation was observed between f esc and M * in the VANDELS analysis when utilizing maximumlikelihood estimation (MLE) to determine f esc , consistent with trends found using indirect measurements of f esc in VANDELS (Saldana-Lopez et al. 2022a ). No correlation was found when using a Bayesian estimate of f esc . The trend between the VANDELS MLE f esc values and M * is remarkably consistent with our results. In addition, the reported f esc values for both Bayesian and MLE methods from VANDELS are consistent with our f esc constraints at comparable M * ; ho we ver, we note that the VANDELS results use modelling that more closely resembles the 'screen' model of Steidel et al. ( 2018 ), rather than the 'holes' model used in this work. Using the 'screen' model results in ∼ 30 per cent higher f esc than using the 'holes' model in the KLCS (Steidel et al. 2018 ), which is still consistent with the VANDELS results. The VANDELS analysis also reco v ered a significant anticorrelation between f esc and UV dust attenuation, where dust attenuation was quantified in terms of the UV slope, β. These results are qualitatively consistent with the anticorrelation between f esc and E ( B -V ) we present in the central left-hand panel of Fig. 4 .

Stellar population parameters have also been explicitly correlated with W λ (Ly α) in galaxy surv e ys at z ∼ 2-5. These trends are informative for interpreting LyC escape considering that the strength Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 Figure 5. Inferred f esc as a function of stellar mass from this work alongside trends from observation and modelling. Estimates of f esc for three bins of increasing M * for the KLCS SED sample are presented as yellow boxes, and are identical to values presented in Fig. 4 . The f esc constraints from two bins of increasing M * from 148 z ∼ 3.5 galaxies from VANDELS (Be gle y et al. 2022 ) are displayed as purple circles. Solid, purple circles represent f esc fit by maximum-likelihood analysis, while skeletal purple circles represent f esc fit by Bayesian analysis. Predictions for f esc from the FIRE-2 cosmological simulations at a particle mass of m B ∼ 7000 M are displayed as horizontal bars (Ma et al. 2020 ). Predicted f esc as a function of M * at z ∼ 4 for the fiducial model of Naidu et al. ( 2020 ) are displayed as dark circles.

of Ly α emission is similarly modulated by the neutral gas co v ering fraction (e.g. Steidel et al. 2010Steidel et al. , 2011Steidel et al. , 2018 ; ;Verhamme et al. 2015 ;Reddy et al. 2016b ). Stacks of rest-UV spectra at z ∼ 2-5 have demonstrated anticorrelations between W λ (Ly α) and both M * and SFR (Du et al. 2018 ;Pahl et al. 2020 ), mimicking the anticorrelations between f esc and these parameters that we presented in Fig. 4 . Meanwhile, surv e ys at this redshift have shown either no strong correlation between W λ (Ly α) and age (Du et al. 2018 ;Pahl et al. 2020 ) or only a weak correlation (Reddy et al. 2022 ). These analyses of Ly α escape in combination with our f esc trends indicate that stellar age may not as closely linked to the configuration of neutral-phase gas in the ISM and CGM of a galaxy as much as other galaxy properties, such as M * , L UV , E ( B -V ), and SFR.

While the z ∼ 3-4 Universe is an excellent laboratory to test LyC escape physics in galactic environments more similar to those at z > 6, LyC surv e ys in the local Universe are afforded advantages such as the ability to examine the direct ionizing signals from intrinsically fainter galaxies in the dwarf galaxy regime, which may dominate the ionizing background during the epoch of reionization (Robertson et al. 2015 ;Finkelstein et al. 2019 ). In addition, local surv e ys a v oid the sightline variability of the IGM that necessitates binning at z ∼ 3 (Steidel et al. 2018 ), enabling constraints on f esc for individual objects.

The Low-redshift Lyman Continuum Surv e y (LzLCS) analysed 66 galaxies at z = 0.2-0.4 observed with the HST /COS, and reported 35 galaxies individually detected in LyC (Flury et al. 2022a , b ). The galaxies were indirectly selected to be strongly leaking using [O III ] λ5007/[O II ] λλ3726, 3729, SFR surface density, and UV spectral slope, in contrast the LBG-selected KLCS. The correlation between a number of galaxy properties and f esc were considered, where f esc was inferred from stellar-population synthesis fits to COS UV spectra, similar to our determinations of f esc for the KLCS SED sample. The LzLCS f esc values appear to decrease as a function of increasing M * , consistent with the ne gativ e trend, we present in the upper left-hand panel of Fig. 4 . Ho we ver, the correlation coefficient between f esc and M * was determined not to be significant, mirroring other local explorations of the two variables (Izotov et al. 2021 ). Augmenting this result, an examination of the LzLCS sample in tandem with archi v al observ ations (totaling 89 star-forming galaxies at z ∼ 0.3) found that galaxies at lower M * tend to have both bluer spectral slopes and higher f esc (Chisholm et al. 2022 ). This analysis focused on the strong inverse correlation found between f esc and the UV spectral slope at 1550 Å ( β 1500 ) in the expanded sample. While UV spectral slope can encapsulate both the intrinsic spectral slope of the stellar population and the degree of dust reddening in the UV, β 1500 was strongly correlated with E ( B -V ) and uncorrelated with stellar age, indicating that β 1500 for the LzLCS is primarily reflecting the degree of dust reddening (also see Saldana-Lopez et al. 2022b ). Apparent anticorrelations between f esc and β 1500 are supported by our anticorrelation of f esc and E(B-V) seen in the center left panel of Fig. 4 . None the less, Chisholm et al. ( 2022 ) make predictions for f esc vs. M UV at z ∼ 3 that are too low when compared to the f esc = 0.06 ± 0.01 of the KLCS (at M UV ∼ -21), despite reproducing our qualitative relationship between f esc and L UV shown in the center right-hand panel of Fig. 4 .

A weaker, but still significant trend of f esc and sSFR was also observed in the LzLCS, contrasting with the lack of trend between f esc and sSFR that we presented in the central panel of Fig. 4 . We note that the LzLCS modelling allows for arbitrarily-short stellar ages, in contrast with our SED fitting procedure, which ensures ages are greater than typical dynamical time-scales (50 Myr). Finally, f esc is strongly correlated with W λ (Ly α) in the LzLCS analysis, which is broadly consistent with the correlation found for both the KLCS SED sample in Fig. 4 and the full KLCS in Steidel et al. ( 2018 ) and Pahl et al. ( 2021 ).

We also note that strong correlations are found between f esc and Ly α peak separation ( v sep ) and star-formation rate surface density ( SFR ) in the LzLCS. These potential tracers of f esc remain unconfirmed at z ∼ 3, and, in particular, the number of galaxies with HST imaging in the KLCS remains insufficient for testing the connection between f esc and SFR (Pahl et al. 2022 ). Ho we ver, future work will examine potential connections between f esc and v sep at z ∼ 3 by leveraging higher-resolution spectroscopy of the Ly α profiles of KLCS galaxies and f esc and SFR in the VANDELS surv e y.

Theoretical predictions for f esc in a variety of galactic environments can help elucidate fundamental relationships between galactic physics and escaping ionizing radiation in the earliest galaxy populations, where direct LyC detections are impossible. The feedback in realistic environments (FIRE-2; Hopkins et al. 2018 ) project was coupled with radiative transfer in post-processing to examine f esc in cosmological zoom-in simulations of galaxies, evolved down to z = 5 (Ma et al. 2020 ). An increase of f esc with increasing mass was found up to log (M * / M ) ∼ 8, and a subsequent decrease in f esc was found at log (M * / M ) > 8. The increasing relationship between f esc and M * was determined to be due to an increasing efficiency of star formation and feedback, while the decrease at the high-mass end can be explained by increasing dust attenuation. We display this trend as dark horizontal bars in Fig. 5 , specifying the simulation resolution that extends to the stellar masses of our sample (baryonic particle mass m b ∼ 7000 M ). The ne gativ e trend between f esc and M * at log (M * / M ) > 8 found in Ma et al. ( 2020 ) is consistent with the ne gativ e trend we find in the KLCS SED sample, which has a median log (M * / M ) = 9.6. We do find o v erall lower f esc values than the FIRE-2 results at fixed M * . This discrepancy may be expected considering higher f esc values were found with increasing redshift at fixed M * in the simulated galaxies.

Evidence of a turno v er in the relationship between f esc and M * was also found in Kostyuk et al. ( 2022 ), which utilized the IllustrisTNG (Marinacci et al. 2018 ;Naiman et al. 2018 ;Nelson et al. 2018Nelson et al. , 2019 ; ;Pillepich et al. 2018 ;Springel et al. 2018 ) cosmological simulations coupled with the radiative transfer code CRASH (Graziani, Maselli & Ciardi 2013 ). These authors also found significant scatter in the relationship between f esc and M * , due to both differences in ionizing photon production rates and the distribution of stars within the neutral ISM. The SPHINX simulations have also similarly reported this relationship, finding a f esc trend that peaks at log (M * / M ) ∼ 7 and drops off strongly at lower and higher masses (Rosdahl et al. 2022 ).

Finally, Ma et al. ( 2020 ) also explored potential synchronization of periods of intense star formation and ele v ated f esc v alues. They find that feedback from star formation clears sightlines in the ISM and CGM of a galaxy, creating fa v ourable conditions for ionizing photons to escape. This process leads to a correlation between a b urst of star -formation and high f esc , albeit with a few Myr time delay. If true, one might expect a higher f esc in galaxies with shorter stellar ages and ele v ated sSFR. We find no correlation between f esc and these two properties in the KLCS SED analysis, as shown in Fig. 4 . The absence of an observed trend could be explained by a less bursty SFH than those found in Ma et al. ( 2020 ), which would reduce potential dependencies between f esc and stellar age.

Models of reionization and their predicted timelines are built upon assumptions regarding f esc , which are impossible to constrain directly in the reionization era. Some assume single values of f esc for all galaxies for simplicity (typically 10-20 per cent ; Robertson et al. 2015 ;Ishigaki et al. 2018 ), others assume that f esc depends on halo mass (Finkelstein et al. 2019 ), or that f esc depends on one particular galaxy property (Naidu et al. 2020 ;Matthee et al. 2022 ).

We compare our f esc versus M * trend to the predictions of the fiducial model of Naidu et al. ( 2020 ), which concludes that reionization is 'oligarchical,' such that the most luminous, massive [ M UV < -18 and log (M * / M ) > 8] galaxies at z > 6 contribute the bulk of the ionizing photon budget. In this model, a direct relationship between f esc and SFR is assumed such that f esc = 1.6 × 0 . 4 SFR . As massive and UV bright galaxies tend to have high SFR , the assumed connection between f esc and SFR results in a positive relationship between f esc and M * . We display this trend determined at z = 4 as dark circles in Fig. 5 . The trend we observe between f esc and M * at z ∼ 3 is inconsistent with the direction and magnitude of the model curve within the mass range where the model and observations o v erlap. Specifically, we show that f esc decreases with increasing M * at log (M * / M ) ∼ 9.5. We also find significantly lower f esc at fixed M * than is predicted by the model. Some of this offset at fixed mass is likely due to the average value of f esc = 0.09 ± 0.01 from Steidel et al. ( 2018 ), used as a constraint in the model, considering f esc of the KLCS was corrected to f esc = 0.06 ± 0.01 after removal of foreground contamination (Pahl et al. 2021 ). Additionally, the fiducial model of Naidu et al. ( 2020 ) does not explicitly consider dust, which we find is a significant factor modulating the escape fraction of galaxies in our sample. To conserv ati vely match observed relationships between f esc and M * at z ∼ 3, assumed f esc values of galaxies at M * > 10 9 M should be no higher than the f esc values of M * ∼ 10 8 . 5 M galaxies. Specifically, the f esc value for the most massive M * data point from Naidu et al. ( 2020 ) should shift to become lower than or equal to the f esc value for the second-most-massive data point. Meeting this requirement at z ∼ 4, which is the closest point of contact between the Naidu et al. ( 2020) model and our observations, would require a reduction in f esc for M * > 10 9 galaxies by a factor of two in the fiducial model of Naidu et al. ( 2020 ). Satisfying this criterion at z > 6 during the epoch of reionization would require a similar reduction. This adjustment would significantly shift the burden of reionization to lower mass galaxies (e.g. Finkelstein et al. 2019 ).

The rapidity of reionization depends strongly on the population of galaxies that dominates the ionizing emissivity o v er cosmic time.

Our results indicate that fainter, less massive galaxies with lower dust content have conditions favourable for escaping ionizing radiation, broadly consistent with other recent LyC observations at z ∼ 3 (Be gle y et al. 2022 ) and in the local Universe (Flury et al. 2022a , b ;Chisholm et al. 2022 ). If these trends were present within the epoch of reionization, the process of reionization may have started early and progressed gradually, such that the IGM neutral fraction is 20 per cent at z ∼ 7 (Finkelstein et al. 2019 ), in slight tension with neutral fraction constraints from Ly α damping wing measurements (Bolton et al. 2011 ;Greig et al. 2017 ;Ba ˜ nados et al. 2018 ). Chisholm et al. ( 2022 ) calculate the ionizing emissivity between z ∼ 4-8 using empirical relations between f esc and β 1500 found in the LzLCS, which are consistent with our results connecting f esc and E ( B -V ), and match constraints indicating that the ionizing emissivity flattens out at z < 5.5 (Becker & Bolton 2013 ;Becker et al. 2021 ). Ho we ver, as noted in Section 4.1 , average f esc values assumed by Chisholm et al. ( 2022 ) are too low at z ∼ 3 when compared to the KLCS. As an alternative scenario, Matthee et al. ( 2022 ) instead directly tie f esc to the strength of Ly α emission and build a model that produces rapid reionization and a flattened evolution of the ionizing emissivity at z < 6. Predicted trends between f esc and L UV from the model of Matthee et al. ( 2022 ) qualitatively match the KLCS anticorrelation, but underpredict the average f esc at z ∼ 3. None the less, such a prescription is promising considering that the relationship between f esc and W λ (Ly α) appears to be one of the most fundamental in our analysis. The Matthee et al. ( 2022 ) model assumes f esc = 50 per cent for half of Ly α emitters with L Ly α > 10 42 erg s -1 , based on the Ly α line profile shapes of z ∼ 2 Ly α emitters (Naidu et al. 2022 ). Our ongoing spectroscopic observing program to explore the connection between Ly α profile shape and LyC escape in the KLCS will test this formalism, which relies on a correlation between f esc and Ly α peak separation that currently lacks direct observational support at high redshift.

Both the Chisholm et al. ( 2022 ) and Matthee et al. ( 2022 ) models highlight important existing relationships found between f esc and galaxy property in our analyses, and present ionizing emissivities that both o v ercome recombination in the IGM at z ∼ 8 and a v oid o v erproducing ionizing photons at z < 6. The trends between f esc and galaxy properties presented in this work are vital for anchoring assumptions of f esc during the epoch of reionization, where direct constraints on f esc are impossible. Future reionization models can utilize these relationships to ensure consistency between reionization-era f esc prescriptions and our empirical results, particularly in comparable galaxy populations that have similar luminosities and masses to those of our sample. We will extend our analysis of f esc and galaxy properties to lower L UV in future work, which will elucidate most Downloaded from https://academic.oup.com/mnras/article/521/3/3247/7080157 by California Institute of Technology user on 31 July 2023 MNRAS 521, 3247-3259 (2023) fundamental predictors of f esc for a larger dynamic range of galaxy properties.

In this work, we examine the underlying processes behind the escape of ionizing radiation by exploring trends between f esc and galaxy properties at z ∼ 3. We accomplish this goal by leveraging multiband photometry of galaxies observed spectroscopically as part of KLCS. We examined a subsample of 96 KLCS galaxies with photometry suitable for SED fitting, and determined galaxy population parameters of M * , SFR, sSFR, E ( B -V ), and age from these stellar-population synthesis fits. For each galaxy property, we sorted the 96 galaxies and divided them into three equal-sized bins, constructing a rest-UV composite spectrum for each bin. The main results regarding the estimated Lyman-continuum escape parameters of f 900 / f 1500 out and f esc and their relationships with galaxy properties are as follows:

(i) We find significant correlations between f esc and W λ (Ly α) and anticorrelations between f esc and L UV in the KLCS SED subsample, indicating that our sample is representative of the full KLCS and appropriate for constraining f esc as a function of other galaxy properties (Pahl et al. 2022 ).

(ii) We find significant anticorrelation between f esc and E ( B -V ) across three bins of increasing E ( B -V ), although no correlation between f 900 / f 1500 out and E ( B -V ). The f esc result indicates that dust modulates escaping ionizing radiation at z ∼ 3. Such modulation naturally arises due to the spatial coincidence of neutral-phase gas and dust (Reddy et al. 2016b ;Du et al. 2018 ;Pahl et al. 2020 ), and the fact that dust directly absorbs LyC photons (Reddy et al. 2016a ). These results are broadly consistent with anticorrelations found between f esc and UV spectral slope at z ∼ 3.5 (Be gle y et al. 2022 ) and in the local Universe (Flury et al. 2022a , b ;Chisholm et al. 2022 ).

(iii) Both f 900 / f 1500 out and f esc are significantly correlated with M * and SFR. Trends between f esc and M * have also been suggested in other LyC surv e ys at high and low redshift (Be gle y et al. 2022 ;Flury et al. 2022a , b ). The sense of the relationships we observe is consistent with reco v ered anticorrelation between f esc and M * for log (M * / M ) > 8 galaxies in cosmological simulations (Ma et al. 2020 ;Kostyuk et al. 2022 ).

(iv) Some cosmological zoom-in simulations of reionization-era galaxies connect stellar feedback and fa v orable ISM/CGM conditions for LyC escape (e.g. Ma et al. 2020 ), which would plausibly manifest as ele v ated estimates of f esc in galaxies with shorter inferred stellar ages and higher sSFR. Ho we ver, we find no correlation between f 900 / f 1500 out or f esc and stellar age or sSFR, providing no direct observational support for the synchronization of recent bursts of star formation and the escape of ionizing photons at these masses. These trends are consistent with the absent or weak correlation found between W λ (Ly α) and stellar age in earlier work (Du et al. 2018 ;Pahl et al. 2020 ;Reddy et al. 2022 ).

These results represent a comprehensive exploration of f esc and SED-modelled properties at high redshift, grounding assumptions of f esc for galaxies in the reionization era. Significant unknowns still remain for f esc and its dependencies in galaxies less luminous than those in our sample, particularly at z ∼ 3. In future work, we will extend our examination of f esc and galaxy properties down to lower UV luminosities. Additional indirect diagnostics of f esc that hav e pro v en promising in the local Universe can also be tested at high redshift with the KLCS. Ongoing follow-up of the Ly α line profiles of KLCS galaxies will elucidate potential trends between f esc and Ly α peak separation, and Keck/MOSFIRE spectra in hand for a substantial subset of the KLCS will enable an examination of the relationships between nebular emission-line properties and f esc . Using the results summarized in this section in tandem with future analyses of the KLCS, we will attempt to offer a unified picture of escaping ionizing radiation at z ∼ 3. This picture is vital for understanding the contribution of star-forming galaxies to reionization at earlier times.

We acknowledge support from NSF AAG grants 0606912, 0908805, 1313472, 2009313, 2009085, and 2009278. Support for program HST-GO-15287.001 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under N ASA contract N AS5-26555. CS was supported in part by the Caltech/JPL President's and Director's program. Based in part on observations obtained with Me gaPrime/Me gaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at Terapix available at the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Surv e y, a collaborativ e project of NRC and CNRS. This paper includes data gathered with the 6.5 m Magellan Telescopes located at Las Campanas Observatory, Chile. We thank D. Kelson for the use of his FourCLift FourStar Reduction code and for his assistance with it. We wish to extend special thanks to those of Hawaiian ancestry on whose sacred mountain, we are privileged to be guests. Without their generous hospitality, most of the observations presented herein would not have been possible.