Lyman Continuum Galaxy Candidates in COSMOS

Star-forming galaxies are the sources likely to have reionized the universe. As we cannot observe them directly due to the opacity of the intergalactic medium at $z\gtrsim5$, we study $z\sim3\text{--}5$ galaxies as proxies to place observational constraints on cosmic reionization. Using new deep \textit{Hubble Space Telescope} rest-frame UV F336W and F435W imaging (30-orbit, $\sim40$~arcmin$^2$, $\sim29\text{--}30$~mag depth at 5$\sigma$), we attempt to identify a sample of Lyman continuum galaxies (LCGs). These are individual sources that emit ionizing flux below the Lyman break ($<912~\text{\AA}$). This population would allow us to constrain cosmic reionization parameters such as the number density and escape fraction ($f_{\rm esc}$) of ionizing sources. We compile a comprehensive parent sample that does not rely on the Lyman-break technique for redshifts. We present three new spectroscopic candidates at $z\sim3.7\text{--}4.4$, and 32 new photometric candidates. The high-resolution multi-band HST imaging and new Keck/Low Resolution Imaging Spectrometer (LRIS) redshifts make these promising spectroscopic LCG candidates. Using both a traditional and probabilistic approach, we find the most likely $f_{\rm esc}$ values for the three spectroscopic LCG candidates are $>100\%$, and therefore not physical. We are unable to confirm the true nature of these sources with the best available imaging and direct blue Keck/LRIS spectroscopy. More spectra, especially from the new class of 30 m telescopes, will be required to build a statistical sample of LCGs to place firm observational constraints on cosmic reionization.


INTRODUCTION
During the Epoch of Reionization (EoR), the universe was transformed from almost completely neutral to nearly entirely ionized (e.g., Loeb & Barkana 2001). This transition was thought to occur at z 6 (e.g., Fan et al. 2006;Pentericci et al. 2011;Mason et al. 2018), but recent results show that it could have continued inhomogeneously in patches down to z 5 (e.g., observations: Becker et al. 2015;Bosman et al. 2018;Eilers et al. 2018;Yang et al. 2020;Kusakabe et al. 2020;theory: Kulkarni et al. 2019;Keating et al. 2020;Nasir & D'Aloisio 2020). The sources responsible for cosmic reionization are not yet known but it is thought that star-forming galaxies are the main contributors. This is because the number density of quasars has been found not to significantly contribute to the ultraviolet (UV) background beyond z ∼ 3 (e.g., Hopkins et al. 2007;Jiang et al. 2008;Fontanot et al. 2012).
Although star-forming galaxies remain the most likely source of ionizing flux during the EoR, the type responsible for the bulk of the UV background radiation remains elusive. Many studies of the EoR have focused on abundant, young, low-mass, star-forming galaxies as the main contributors (e.g., Ouchi et al. 2009;Wise & Cen 2009;Yajima et al. 2011;Finkelstein et al. 2012;Bouwens et al. 2015;Paardekooper et al. 2015;Robertson et al. 2015). However, recent modeling found that most of the ionizing photon budget ( 80%) could come from massive, UV-bright galaxies with high escape fractions of ionizing radiation (f esc ; Naidu et al. 2020). This theory is bolstered by observations of ultra-luminous Lyman-α (Lyα) emitters during the EoR (z ∼ 6.6; e.g., Hu et al. 2016;Matthee et al. 2018;Songaila et al. 2018;Meyer et al. 2021;Taylor et al. 2020). These powerful UV-flux emitters are thought to be located at the sites of massive reionization bubbles, which is supported by results from simulations (e.g., Gronke et al. 2020).
Given the opacity of the intergalactic medium (IGM) at z 5 to UV radiation, direct observations of the flux that ionized the universe are not possible (e.g., Inoue & Iwata 2008;Inoue et al. 2014). However, we require more observational constraints on f esc to understand reionization. While we cannot directly observe the photons that reionized the universe, we can look to z ∼ 3-5 galaxies that serve as close proxies to those during the EoR. Identifying high-redshift galaxies that emit Lyman continuum (LyC) flux, ionizing UV radiation below the Lyman break (< 912Å), has long been a focus of observational cosmology. We can calibrate observable properties of galaxies during the EoR, such as optical nebular emission lines, to the LyC flux emission. We can therefore determine this relation for z ∼ 3-5 galaxies and extrapolate it to z > 6 (e.g., Bassett et al. 2019).
One possible explanation for the lack of detections relative to the extensive LyC search campaigns is the methods used to select possible candidates. Identifying high-redshift LyC-emitter candidates requires secure redshifts and a method of detecting leaking LyC flux. An efficient way of measuring the redshifts of distant UV-bright star-forming galaxies is the Lymanbreak technique (Steidel et al. 1996). This method uses large area multi-broadband imaging to efficiently identify high-redshift galaxies, rather than expensive deep spectroscopy of individual sources. An important question to ask is if LyC flux is added below the Lyman break, is the galaxy always identifiable as a Lymanbreak galaxy (LBG)?
This was the hypothesis that Cooke et al. (2014) explored by artificially adding LyC flux of various strengths to LBG spectra and measuring the resulting colors. Meštrić et al. (2020) also tested this theory by searching in new CFHT Large Area U -band Deep Survey (CLAUDS; Sawicki et al. 2019) images for LyC flux for galaxies with publicly available spectroscopic redshifts. This search resulted in two new candidates. As found by Cooke et al. (2014); Meštrić et al. (2020), these "Lyman continuum galaxies" (LCGs) are typically outside of the LBG color-color selection boxes. A fraction of LBGs with LyC flux still reside in the traditional LBG selection region (see Cooke et al. 2014). However, their location is inconsistent with the model predictions for their redshift. All galaxies that have LyC flux are included in the LCG definition. However, understanding those outside of the LBG selection region, that may be missed by other searches, is a primary goal of this work.
Even spectroscopic redshifts in the literature can have selection biases that favor LBGs. Targets for spectroscopic surveys are usually pre-selected to have a representative distribution using either color selection or photometric redshifts. Galaxies are likely to have the securest photometric redshifts if they exhibit a clean Lyman break. Therefore, spectroscopic redshift selected samples can also be biased against strong LyC detection.
In this work, we present new deep HST imaging of the Cosmic Evolution Survey field (COSMOS; Scoville et al. 2007) to investigate LCGs. We use accurate photometric redshifts from the FourStar Galaxy Evolution Survey (ZFOURGE; Straatman et al. 2016), and 18 other public photometric and spectroscopic redshift catalogs, to identify candidates. ZFOURGE uses 30-band photometry in the COSMOS field from 0.5-8 µm, including medium infrared (IR) bands that sample the Balmer break at 1 z 4. This enables accurate photometric redshifts without relying on the Lyman break that may be less prominent due to escaping LyC. We attempted to spectroscopically follow up our LCG candidates using Keck/Low Resolution Imaging Spectrometer (LRIS) and other telescopes. We present our results from the photometric and spectroscopic candidate searches here.
We describe the HST observations, custom reduction, and photometry in Section 2. We outline the redshift catalogs and selection methods used to identify LCG candidates in Section 3. We detail our various spectroscopic follow-up campaigns and reduction in Section 4. The sample refinement, properties, and analysis are described in Section 5. We discuss these results in the context of the literature in Section 6. We summarize the paper and give our main conclusions in Section 7. Throughout this paper we use AB magnitudes (m AB ; Oke & Gunn 1983) and assume the latest Planck cosmological parameters (H 0 = 67.4 kms −1 Mpc −1 , Ω m = 0.315; Planck Collaboration et al. 2020).

Observations
High-resolution HST data are required for this work to identify any potential low-redshift contaminants along the line of sight. The new HST data presented here were taken in Cycle 25 (GO 15100; PI Cooke; 30 orbits) in three blocks between 2018 February 21 and 2018 April  Shapley et al. (2003). The LBG spectrum is redshifted to z = 3.087 (F336W z lim ; gray) and z = 4.391 (F435W z lim ; black). Four HST bands are shown (left to right): new WFC3/F336W (purple), new ACS/F435W (light blue), CANDELS-ACS/F606W (light green), COSMOS-ACS/F814W (gold), and three HSC bands (left to right): g (dark blue), r (dark green), i (orange). Key spectral features are highlighted on the black (z = 4.391) spectrum with blue vertical lines: Lyman break at rest ∼ 912Å (dashed), Lyα at rest ∼ 1216Å (dot-dashed), non-ionizing UV continuum rest ∼ 1500Å (dotted).
15. The data were taken using two filters in parallel: F336W on Wide Field Camera 3 (WFC3)/UVIS and F435W on Advanced Camera for Surveys (ACS)/WFC.
For galaxies at z = 3-5, LyC flux (< 912Å) is cleanly sampled in the WFC3/F336W band at z > 3.087 and in the ACS/F435W band at z > 4.391. These limits assume a < 0.3% red leak ( 8 mag fainter than the UV continuum) as used in Smith et al. (2018) and we refer to these limits as z lim . Figure 1 shows an LBG composite (∼ 200 spectra binned by Lyα strength) from Shapley et al. (2003). The LBG composite is redshifted to z lim = 3.087 (gray) and z lim = 4.391 (black). Transmission curves for photometric bands used in this work are overlaid. The HST bands: WFC3/F336W, ACS/F435W, ACS/F606W, ACS/F814W, and three ground-based Hyper-Suprime Cam (HSC) bands: g, r, i are shown. We also highlight key features on the LBG template at z = 4.391 (black) with vertical dark blue lines: the Lyman break at rest ∼ 912Å (dashed), Lyα at rest ∼ 1216Å (dot-dashed), and the non-ionizing UV continuum at rest ∼ 1500Å (dotted).
The 30 orbits of HST data were originally intended to be taken in two F336W pointings with two parallel F435W pointings. Given observation scheduling constraints, the first F336W pointing (F336W 1) was observed in two overlapping visits with different position  (Koekemoer et al. 2007) with survey footprints relevant to this work and our new LCG candidates overlaid. Two deep WFC3/F336W (purple) and three parallel ACS/F435W (light blue) pointings with their names and orbit depths shown. The CANDELS-ACS/F606W (light green; Grogin et al. 2011;Koekemoer et al. 2011) and ZFOURGE survey (yellow;Straatman et al. 2016) footprints are also shown. The LCG candidates are split into three categories: spectroscopic candidates/Group A (red stars and ZFOURGE IDs), and photometric candidates, both Group B (blue circles) and C (dark blue diamonds, lowest priority candidates). See Section 5.1 for more details on the LCG candidate group categorization. angles (PAs) for a total of 15 orbits. This means that the parallel F435W pointings are split up into F435W 1a (5 orbits) and F453W 1b (10 orbits). The second F336W pointing (F336W 2) and F435W parallel (F435 2) both include the full 15 orbits. The total area of all the pointings is ∼ 40 arcmin 2 at various depths. Figure  2 shows all five HST pointing footprints (two F336W and three F435W), with names and respective orbit depths, overlaid on the COSMOS-ACS/F814W background image cutout (Koekemoer et al. 2007). We also show the CANDELS-ACS/F606W (Grogin et al. 2011;Koekemoer et al. 2011) footprint (light green) overlaid.
We use full-orbit exposures to reduce the total background contribution caused by post-flashing each exposure to 12 e − (as recommended in Anderson et al. 2012). The longer exposures also reduce read noise and maximize the counts per exposure for faint sources. This enables improved charge transfer efficiency (CTE) correction which is better suited for higher count levels. We use a five-point dither pattern per visit to sub-sample the point-spread function (PSF). We employ a larger dither pattern between each visit to remove the chip gap and low-level, large-scale pattern noise (Rafelski et al. 2015). The total exposure time of both F336W pointings shown in Figure 2 is ∼ 44000 s. The exposure time of the deepest F435W pointing (F435W 2) is ∼ 40000 s, and this exposure time is split by a 1:2 ratio between F435W 1a and F435W 1b respectively. See Table 1 for a summary.

HST data reduction
Given the faint nature of potential LCG candidates in the rest-far UV (FUV), having the best quality HST images is essential. We therefore make several improvements to the standard Space Telescope Science Institute (STScI) WFC3/UVIS reduction pipeline, to clean and calibrate our images beyond what was previously available from the Mikulski Archive for Space Telescopes (MAST). We also developed and apply new additional calibrations to the CTE-corrected single-exposure calibrated images (called FLCs). Some of these improvements have been added to the standard WFC3/UVIS MAST pipeline (detailed in Section 2.2.1), while others are additional calibration improvements (outlined in Section 2.2.2). We provide further information and plots ( Figure A1) to demonstrate our improvements in Appendix A.

HST/WFC3 darks pipeline improvements
The standard darks pipeline for WFC3/UVIS on HST left some residual systematic and instrumental effects. Previously, Rafelski et al. (2015) identified the cause and corrected for some of the worst of these residual effects in the WFC3/UVIS darks. In this work, we adapt the STScI WFC3/UVIS darks Python pipeline (copied April 19, 2019; Bourque & Baggett 2016). The improvements include those outlined in Rafelski et al. (2015) that were previously implemented in custom IDL routines, and we add further new improved calibrations. Here we briefly outline the major changes we make to the WFC3/UVIS darks pipeline which have since been adopted by the WFC3 team as standard for their delivered data products (released May 2021). The WFC3/UVIS MAST data also include the new time-dependent photometry flux calibrations that are essential for accurate photometry (Bohlin et al. 2020;Bajaj et al. 2020;Calamida et al. 2021).
• A new charge transfer efficiency (CTE) correction (calwf3 v3.6.0; Anderson 2020; Anderson et al. 2021) is used. The new code uses a pixel-based algorithm that optimizes parameters for the data to mitigate read noise. This reduces background noise and improves the correction of bright sources such as cosmic rays which reduces residual gradients across the images.
• The detectors are periodically heated in an anneal to reduce hot pixels. The period between each heating is referred to as an "anneal cycle", within which the detector has a different "fingerprint". Darks are now processed within concurrent anneal cycles (i.e., not using the darks from the previous anneal cycle for pixel replacement). This results in cleaner corrections to the data, and reduced blotchy patterns and noise (see also appendix in Rafelski et al. 2015).

Additional calibrations
The improvements described so far have been folded into the official WFC3/UVIS pipeline coupled with the new time-dependent photometry. We therefore download the WFC3/F336W data presented in this paper from MAST (on May 28, 2021) as a starting point. We then apply additional corrections described below to these data to further improve their calibration. We make no changes to the ACS reduction pipeline, so the F435W FLCs are those from MAST. A Jupyter notebook of the additional calibrations and some steps required to easily apply them to WFC3 and ACS (if relevant) FLCs is available on GitHub 1 . We briefly summarize these corrections below and give more details in Appendix A.
• We developed a novel hot pixel flagging routine to apply a variable flagging threshold as a function of distance from the "readout" amplifier. The method used in the MAST data uses a constant threshold. Due to the imperfect CTE correction, this results in ∼ 30% of hot pixels being missed further from the readout.
• We flag the negative divots adjacent to readout cosmic rays (ROCRs) in the data quality (DQ) arrays of each FLC. These ROCRs land on the array while it is being read out and are overcorrected by the CTE code.
• We equalize the overall signal levels in the four different amplifier regions (quadrants) of each FLC. This produces smoother images without the background discontinuities.

Ancillary data
We collate publicly available ancillary HST imaging to aid with source detection and for escape fraction estimates. At z 3.1 (z lim for F336W), the non-ionizing UV continuum around rest ∼ 1500Å is partially covered by F606W. For z 4.4 (z lim for F435W), the UV continuum lands in F814W. See Figure 1 for the positions of the bands on spectra shifted to both z lim values. We therefore use the CANDELS-ACS/F606W (Grogin et al. 2011;Koekemoer et al. 2011) 2 and COSMOS-ACS/F814W (Koekemoer et al. 2007) 3 images in our search and analysis of LCGs. See Figure 2 for the F606W footprint (light green) and the F814W band (background image). The COSMOS-F814W image spans the entire COSMOS field (∼ 1.4 deg × 1.4 deg), and thus we have F814W spanning our new HST images.
We combine the F814W tiles spanning the new HST pointings and ZFOURGE footprint with the Montage v6.0 software 4 . We then reconstruct the header World Coordinate System (WCS) information of the resulting mosaic with the tweakreg routine from the Driz-zlePac package (Gonzaga et al. 2012;Hoffmann et al. 2021). This step makes the WCS compatible for use with astrodrizzle that is used for combining the new HST images. The F606W and F814W mosaics are cropped down to the area of interest covering the new HST pointings, shown in Figure 2. Both the F606W and F814W images have 0.03 pixel scales and are aligned to the COSMOS survey coordinates. The CANDELS-F606W images have an exposure time of 3300s and depth of 28.5 mag at 5σ (Grogin et al. 2011;Koekemoer et al. 2011). The COSMOS-F814W images have an exposure time of 2028s with a 27.2 mag point-source depth at 5σ (Koekemoer et al. 2007). See Table 1 for a summary of image depths.

Creating mosaics
The new F336W and F435W FLCs are drizzled into two separate images using tools from DrizzlePac (v3.1.6 for F435W, v3.2.1 for F336W; Gonzaga et al. 2012;Hoffmann et al. 2021). The updatewcs (with use db=False), tweakreg, and tweakback packages are used for image alignment. We developed a tool that we use for cleaning image edges to improve alignment with tweakreg (cleanedges 5 ). The new F336W and F435W FLCs are drizzled onto the COSMOS coordinates using the F814W band as a reference image for the same footprint. We create both an exposure time map ('EXP') and inverse variance weight map ('IVM') for each of the new filter drizzles. The final F336W and F435W mosaics have 0.03 pixel scales and are aligned with pixel orientation north. See Table 2 for a summary of the final drizzle parameters.
We measure the depth of the images by taking the sigma-clipped median root-mean-squared (RMS) of 1000 apertures placed on sky (avoiding sources) in the HST pointings with our hst sky rms code 6 . We then measure the flux required for a 5σ detection (relative to our median sky RMS) in a 0.2 radius aperture and convert this to a magnitude for each filter with timedependent zeropoints from STScI 7 . The resulting RMS image depths of F336W 1 (15 orbits, misaligned PAs) is m AB = 29.77 at 5σ and for F336W 2 (15 orbits) is m AB = 30.11 at 5σ. The RMS magnitude depths for the three F435W pointings with varying orbit depths are m AB = 28.52 at 5σ for F435W 1a (5 orbits), m AB = 29.11 at 5σ for F435W 1b (10 orbits), and m AB = 29.78 at 5σ for F435W 2 (15 orbits). See Table 1 for a summary of image properties. All our F336W and F435W data products are available as a High Level Science Product at MAST via 10.17909/t9-pe8e-s980 8 .
We use the make model star image tool from the psf tools.psfutils module in the wfc3 photometry package to plant empirical PSF models into copies of our individual FLCs and drizzle them together. The resulting image includes a grid of ∼ 1000 empirical PSFs that can be stacked like real stars. We stack the PSFs in the images using new PSF 9 https://github.com/spacetelescope/wfc3 photometry stacking tools 10 . We extract each PSF and interpolate them onto a sub-pixel grid. The PSFs are fitted with Moffat profiles, checked for quality, aligned, and the median is taken. We interpolate the sub-pixel sample average PSF profiles back to the native pixel scale of the images and determine a full-width-half-max (FWHM) value (see Table 1).
For the F336W and F435W images we stack the empirical PSFs. For F606W and F814W we stack 77 real stars. We create a convolution kernel using a Hanning Window to translate the PSF of each band to the F814W PSF with the largest FWHM. We convolve all images using their respective kernels to obtain PSF-matched images. For the new F336W and F435W images, we apply different PSF corrections for each pointing given their different depths and PAs.

Isophotal flux measurements
To identify sources in the PSF-matched images and extract their flux, we use tools from the photutils Python package. The benefit of isophotal photometry, i.e., extracting flux for each object within segmentation footprints, is that it cleanly extracts flux for the faint blue sources that we are targeting. We first create RMS maps from the inverse variance weight maps (IVM) output by astrodrizzle (RMS = 1/ (IVM)). We then create a noise-equalized (NEQ) image of all our PSF-matched HST bands. To do this, we divide the PSF-matched images by their RMS maps to create signal-to-noise (S/N ) maps, mask invalid values, and sum them together. We also create negatives of all the images (× − 1), including the NEQ image, that we use to aid with source detection.
We make two different segmentation maps from the NEQ images using the photutils.detect threshold and photutils.detect sources routines. We create a clean map with few spurious objects for source detection and flux extraction (i.e., isophotes), using a 2D Gaussian kernel with standard deviation of 3 pixels, threshold of 1.8, and 15 connected pixels for detection (npixels). We test values for our clean segmentation map by finding the threshold required to find few sources in our negative images while maximizing flux extraction of the real sources. We also make a low threshold map (0.8, npixels=10) with many spurious sources for rigorously masking objects when measuring a local sky value (sky masks).
We perform isophotal photometry with photutils by measuring flux within the source footprints determined from the NEQ image and giving this as a mask to the aperture photometry routine. We first determine a local sky value in an annulus around each source. We heavily mask sources with our rigorous sky mask and take a 3σ-clipped median of the "sky" pixels. We multiply this sky value by the isophotal area of the source and subtract it from the source flux.
To calculate the errors on the isophotal fluxes we use the photutils.calc total error function. To account for correlated noise that results from drizzling larger native pixels onto a finer grid, we determine a correlated-noise correction factor (as described in Fruchter & Hook 2002;Hoffmann et al. 2021). The photutils.calc total error function then takes the image data, the RMS map multiplied by the correlated noise factor, and the exposure-time map to determine an error map that includes shot noise. The photutils.aperture photometry routine takes the image data with the local background subtracted, the isophotal extraction footprint, and the total error array to generate accurate photometry and errors.

LyC contamination calculations
We measure the probability of LyC flux contamination for the new F336W and F435W bands using a similar method to determining the image depths (included in our hst sky rms code 11 ). We make a 3σ segmentation map for each pointing and randomly place 10,000 apertures onto it. We use an aperture width of 0.4 which is the average size of our LCG candidates (see Section 5.1). If there are 15 source pixels within any randomly placed aperture (our npixels value for source detection), it is considered contaminated. Taking a ratio of the clean and contaminated 10,000 randomly placed apertures gives us a contamination rate for each pointing. 11 https://github.com/lprichard/hst sky rms We find a 0.025, 0.018, and 0.013 probability of any F435W 1a-, F435W 1b-, F435W 2-selected source (respectively) being contaminated. We also find a 0.050 and 0.014 contamination rate estimate for F336W 1 and F336W 2 respectively (see Table 1 for a summary). Although both F336W pointings have 15-orbit depths, F336W 1 has misaligned PAs (see Figure 2) so the image edges are shallower and more noisy. Toward the image center at full depth, the image properties will be close to that of F336W 2. These are conservative estimates as the source flux, size, and redshift are not considered, only randomly overlapping sources in the segmentation map.

Ground-based photometry
As seen in Figure 2, not all our HST images are covered by CANDELS-F606W that we use in our analysis of LCGs. For F435W 2, we require ground-based data that span the field and a similar wavelength range as F606W for making color-color plots, for illustrative purposes rather than sample selection. We opt to use HSC on the 8. Due to the large size of the ground-based PSFs relative to the HST bands, when convolving the F435W band to the HSC-i PSF, the flux from the faint sources that we are interested in is lost altogether. We therefore derive a correction to convert the F435W isophotal fluxes to 1 -aperture-and HSC-PSF-matched fluxes. We achieve this by measuring 1 aperture fluxes on our F435W image convolved to the HSC-i PSF for all sources. We then plot the HST isophotal fluxes against all the valid PSF-and aperture-matched F435W fluxes and derive the following relation for converting between the two.
Here F F435W,1 ap is the F435W flux measured in 1 apertures from the HSC-i PSF-matched image. F F435W,iso is the measured F435W isophotal flux. We apply this correction to the F435W isophotal fluxes to convert them to 1 -aperture HSC-i PSF-convolved fluxes so they can be directly compared to the groundbased HSC bands.

Redshift catalogs
To select LCGs, we require a sample of high-redshift galaxies not selected with the Lyman-break technique. We use the ZFOURGE survey  as our primary source of redshifts as it comprises 30 bands of imaging, including medium IR FourStar bands (Persson et al. 2013) from the 6.5 m Magellan Baade telescope at Las Campanas Observatory. These medium IR bands accurately sample the region around the Balmer break for well-constrained photometric redshifts (δz ∼ 2% quoted accuracy; Straatman et al. 2016) of 1 < z < 4 galaxies (van Dokkum et al. 2009). ZFOURGE spans a roughly 11 × 11 (13 × 13 with dithering) region of COSMOS (yellow footprint in Figure 2). As shown in Figure 2, the ZFOURGE footprint does not span the full area of our new HST F435W images (blue footprints), with around ∼ 1 /3 of F435W 1b and F435W 2 outside the ZFOURGE area.
The ZFOURGE survey is useful in the search for LCGs as these medium bands sample the Balmer break. This means that we are less dependent on the Lyman break, which may be less prominent in these galaxies, for determining photometric redshifts. Strong emission lines can bias photometric redshift measurements. Narrow bands that are sensitive to emission were omitted in the ZFOURGE photometric redshift calculations, the medium bands are less sensitive to emission and were included, but may still bias the results. We therefore draw from multiple redshift sources where possible for selecting LCG candidates.
ZFOURGE sources are K s -band selected in images with 5σ depths of 26.2-26.5 mag at a typical seeing of ∼ 0.4 . As ZFOURGE galaxies are K s -band selected, the candidates must be sufficiently bright in this red band to be detected and we will miss some faint bluer sources. We therefore collate other available redshift information in the COSMOS field to construct a more complete parent sample. These include photometric redshifts from the COSMOS2015 catalog (Laigle et al. 2016) and the CANDELS catalog of COSMOS (Nayyeri et al. 2017).
We compiled a catalog of available spectroscopic redshifts in the field from 3D-HST (Brammer et al. 2012;Momcheva et al. 2016) Lilly et al. 2009). We note that many targets in spectroscopic surveys are pre-selected using color selection and photometric redshifts. Many of these are biased towards objects with the strongest Lyman breaks. Therefore, many of the publicly available spectroscopic catalogs may be biased against bright leaking LyC.
We first spatially crop all catalogs down to our region of interest (shown in Figure 2). We then cross-match the redshift catalogs within 1 apertures (0.5 search radii). We start by matching everything in COSMOS to ZFOURGE, then add any candidates from COSMOS that do not overlap. We then cross-match the other catalogs, starting with CANDELS, and then the spectroscopic redshift catalog. The CANDELS catalog contains photometric redshifts produced by six different authors and methods with the same data. We therefore calculate an average CANDELS photometric redshift (as shown to be most accurate in e.g., Dahlen et al. 2013;Barro et al. 2019) and standard deviation of the available redshifts and use both in our redshift constraints for selecting our samples.
Within our F336W and F435W footprints there are 8311 unique sources at all redshifts (our "All-z Sample"). This reference sample splits to 3065 objects in the F336W footprints and 6430 sources in the F435W footprints (with overlap). We create a comprehensive "Parent Sample" from the All-z Sample by selecting all objects from the catalogs with a redshift above the LyC z lim for the F336W or F435W filters. We remain as agnostic as possible at the candidate selection phase and include any galaxy in the Parent Sample that has at least one redshift to place it above z lim for F336W or F435W to probe LyC. The CANDELS redshift constraint we use for selecting our Parent Sample was our calculated average minus one standard deviation of the CANDELS photometric redshifts. There are 574 unique objects in the Parent Sample: 380 objects have at least one redshift z > 3.087 and are in the F336W footprints and 242 objects have at least one redshift z > 4.391 and are in the F435W footprints (with overlap).

Visual inspection
We visually inspect all objects in the Parent Sample with at least three different investigators examining each. We inspect the sources in all four HST bands described here, their NEQ and source-detection segmentation maps to evaluate if there is any visible LyC flux. If any of the investigators notes a detection of possible flux in the filter probing the rest-frame LyC (F336W and/or F435W), they are re-evaluated multiple times. We remove any obvious erroneous redshifts (large bright local galaxies, star spikes, and bright saturated stars) from the sample. We end up with a long list of 77 unique objects that have varying degrees of confidence both in their LyC detection (strength, offset, confusion over the source for close objects) and redshift. We rank these candidates based on their number of visual confirmations, LyC confidence, and redshift confidence. All 77 sources are used as potential candidates for spectroscopic follow up to remain as agnostic as possible and as a relatively high density on the sky is required for multi-object slit spectroscopy (MOS).

SPECTROSCOPIC DATA
To spectroscopically confirm the redshift and true nature of our LCG candidates, we require deep, blue spectroscopy. One of the few telescope and instrument combinations that satisfy our constraints is the Low-Resolution Imaging Spectrometer (LRIS; Oke et al. 1995;McCarthy et al. 1998;Rockosi et al. 2010) on Keck I. LRIS is a MOS instrument that allows the efficient observation of targets in the same field. The observing site on Maunakea provides enough elevation to reduce the atmospheric extinction. Thicker atmospheres at lower sites can impede the transmission of blue flux, which would require significantly more observing time to achieve our primary science objectives.
Our team has been attempting to spectroscopically confirm LCGs in COSMOS selected independently of the Lyman-break technique since 2014, with four dedicated runs prior to the HST proposal. This work provided the motivation and foundation for the HST program presented here (PID 15100). We have since been awarded time for four dedicated separate observing runs, including three NASA/Keck and one Large Binocular Telescope (LBT), that are part of the HST program. We also applied for time on two additional related projects in COSMOS with some overlapping targets and science goals to this work. Unfortunately, nearly all this time has been weathered out. In total, we have been awarded 12 dedicated nights (over 20 calendar nights) that are part of, or related to, the HST program presented here. In this paper, we present ∼ 1.5 nights of usable Keck/LRIS data taken in moderate observing conditions (see Table 3 for a summary).

Observations
We present Keck/LRIS data of our LCG candidates taken on 2016 February 9 and 10 (2016A W034LA, PI Cooke) and 2020 January 20-22 (2019B N010, PI Rafelski). We require good seeing (< 1 ), clear skies and dark nights for our observations. To spectroscopically detect the LyC, we require the objects to have optical magnitudes of ∼ 24 mag. We then need around a half night (∼ 6 hours) in our optimal conditions to achieve spectroscopic detection of LyC flux. This assumes there is no contamination from neighboring sources or LRIS detector issues. The conditions for the 2016 run were relatively poor, spread over two nights, but we were able to obtain some usable exposures for zf 9775. In 2020, around ∼ 80% of the time was usable but with a relatively poor average seeing of ∼ 1.2 (relative to what is possible at Keck). This poorer seeing disperses light over the edges of the LRIS slit and therefore increases the exposure time required to get sufficient signal-tonoise.
The LRIS instrument uses pre-cut metal masks with slits to disperse the light of our targets across blue and red spectroscopic arms. We use the 400/3400 grism on the blue side and the 400/8500 grating for the red side with the 560 dichroic. We create LRIS masks with the autoslit software. The slit widths are ∼ 1.2 and longer length slits are preferred to improve background subtraction. The masks are designed to encompass as many of our high-priority candidates as possible. We fill any remaining space with lower-priority LCG candidates and fillers that include extreme emission-line galaxies, Lyα emitters and LBGs. We then select at least four stars, ideally five to six, placed as near as possible to the edge of the LRIS masks for acquisition on the sky and mask alignment.
The highest priority objects are placed as close to the north of the detector (with PA=0) as possible due to amplifier issues with the red side. At the time of observations, both chips could not reliably be read out, leaving no coverage at redder wavelengths for half the mask. The best objects are placed in the central third of the detector when possible (in x and y) to accommodate full wavelength coverage and to minimize instrument distortion at the edges. With an optimal distribution on the sky, about ∼ 25 candidates fit on a mask, including ∼ 5 high-priority targets. However, when factoring in all the observational constraints, fewer than that are often in the prime detector and dispersion positions. We aim to observe each mask for a total of ∼ 6 hours in 1200 s exposures to reduce data loss in variable weather conditions and contamination by cosmic rays. Typically, we observed just a few hours per mask if at all.
The Keck/LRIS spectra are reduced using the Image Reduction and Analysis Facility software (IRAF; Tody 1986) following the standard reduction procedure for multi-slit spectroscopy. We also use twilight flats to aid with the blue flux calibrations. A 2×2 binning was used for the 2016 observations and a 1×1 binning was used for the 2020 data.

Determining redshifts
The Keck/LRIS observations yielded redshifts for three LCG candidates: zf 9775, zf 11754, and zf 14000. Only 14 (seven > 3σ in F336W/F435W) galaxies from our long candidate list were successfully observed due to the weather conditions limiting the number of masks we could observe and their spatial distribution. We present spectra for three of the > 3σ candidates, but were unable to get redshifts for the other four observed. A full analysis of the spectroscopic properties of the candidates (e.g., Lyα profiles) is planned to be presented in a future paper. Sources with confident spectroscopic LRIS redshifts and non-detections of LyC are the focus of another paper (Meštrić et al. 2021). Figure 3 (continues over three pages) shows 3-pixel boxcar-smoothed LRIS spectra (black) and a composite of ∼ 200 LBGs with similar Lyα strength to the galaxy overlaid (blue; Shapley et al. 2003). The postage stamps on the left show the LRIS mask slits (top, 30 across) and candidate orientation (bottom, 5 across). Below each spectrum we show zoom-ins on some spectroscopic features. We also highlight other spectroscopic features with shading: the UV continuum red-ward of Lyα ( 1216Å; yellow), the Lyα forest (green), the Lyβ forest (violet), and the LyC (< 912Å; blue panel).
To determine redshifts for the candidates, we compare each spectrum to an LBG composite from Shapley et al. (2003) and look at the fit to the profile of the UV continuum. We also investigate all the highlighted spectroscopic information available, including the flux decrement as a result of the Lyα forest, the presence of a Lyα emission/absorption feature, and the fit to ∼ 6-12 ISM absorption features (depending on the wavelength coverage of the spectrum). We summarize the properties of the three spectra shown in Figure 3 and redshifts we determine from each below.
• zf 9775: We show the blue and red halves of the 2016 LRIS spectrum. The profile of the LRIS spectrum matches the LBG composite, with distinctive breaks and well-aligned ∼ 10 common strong UV ISM absorption features, all assessed to determine the redshift (e.g., see zoom-ins below). This spectrum has a confident redshift of z = 4.365 (with δz ∼ 0.003 for all spectra). However, there is some ambiguity over the nature of this detection. The spectrum is dominated by the optically bright component (left), and we are unable to resolve the component emitting the blue flux (right; see components in Figure  4). Therefore, zf 9775 remains a "spectroscopic candidate" rather than a confirmed LyC detection.
• zf 11754: We show just the red side of the spectrum as the blue side is contaminated by the close bright source in the upper right of the postage stamp. The red half of the spectrum shows a distinct Lyα line (panel below) and break on either side from the UV continuum (yellow) to Lyα forest (green). We are confident in the redshift from the red side of the spectrum of z = 4.470 (using ∼ 6 common ISM lines). However, due to the blue side contamination and the fact that it is so close to the other source, there remains some ambiguity as to the nature of the LyC detection seen in the HST F435W images. This contamination could be increasing the blue flux through the Lyα forest (green) and even red-ward of Lyα. The red side of the spectrum also shows some deviation from the composite which could be real (i.e., a very young population of stars and clean line of sight), or a product of extraction with such a close neighboring source. zf 11754 therefore remains a spectroscopic candidate rather than a confirmed detection.
• zf 14000: There is a tentative broad feature that could be Lyα and this exists at a slight break in the spectrum between what would be the UV continuum (yellow) and Lyα forest (green). This smaller decrement, if indeed a physical Lyα forest, could be due to the source being observed along a cleaner line of sight. This could be the case if considering the strong LyC detection which would require a high IGM transmission sight line. There is contamination seen in the 2D spectrum (upper trace) due to slight slit drift throughout the observations, from alignment star catalog inconsistencies or length of observations. However, the upper trace does not affect the extraction of zf 14000 due to their sufficient separation. The tentative redshift we determine for this source is z = 3.727, and it remains a spectroscopic candidate rather than confirmation. We rule out the possibility that this object is a star as its blue temperature is consistent with it being a K star, but the spectrum does not support that. Its spectral energy distribution (SED) is also not consistent with a late M-type star. Most catalogs list this source as non-stellar, but they do have a variety of redshifts (z = 0.301-3.953, see Table 5). For it to be a K4-M0 star, it needs to be from 30-128 kpc away to be as faint as it is. Only the reddest M dwarf would be realistic (at ∼ 2-10 kpc) and they are too red for the SED. Three of the LCG candidates have low-quality ancillary spectra. Two of the spectra place the candidates above the redshift limit of their respective bands for LyC to be detected if present (z lim ), and one places it below z lim . Two spectra ( Figure B1) are from the DEIMOS 10K survey and are assigned their lowest-quality redshift flag (DM18 hereafter; Hasinger et al. 2018). The spectra are for zf 14000 (DM18 ID: L677152, z DM18 = 3.779) and zf 11423 (DM18 ID: L658665, z DM18 = 3.655). There is tentative feature that could add weight to the LRIS redshift for zf 14000, but with low confidence. The spectrum for zf 11423 shows two very tentative features, and conflicts with three photometric redshifts for this source, so we cannot be confident of its DM18 redshift. There is a 3D-HST G141 grism spectrum (Brammer et al. 2012;Momcheva et al. 2016) for zf 9775 (3D-HST ID: 9997, z 3DHST = 0.446) in Figure B2. There are no strong features to provide evidence for either the 3D-HST or LRIS redshift. The 3D-HST value also conflicts with the closely matching independent redshifts from ZFOURGE, COSMOS, average CANDELS, and our confident LRIS redshift (z ∼ 4.2-4.4). We therefore deem the 3D-HST redshift unreliable. See Appendix B for more details.

Categorizing the LCG candidates
We refine our LCG candidates and categorize them based on our confidence in their detection and redshift. We select only LCG candidates with 3σ flux detections in their LyC probing band (determined from the available redshifts). We then remove candidates with close neighboring sources (within 1.2 apertures), and those contaminated by bright saturated nearby stars. To help rule out the possibility that our candidates are active galactic nuclei we check the spectra (where available) for AGN FUV emission features. We then cross-match the candidate list against the public XMM-Newton X-ray (Cappelluti et al. 2009) and Chandra COSMOS legacy (Marchesi et al. 2016) catalogs within 0.5 radii and find no matches.
This refinement results in a final sample of 35 sources with 26 sources in the F336W pointings and 9 in F435W. With cross-matching and overlap, our final LCG candi-  (3). This is within just our four untargeted HST fields (∼ 40 arcmin 2 ).
The three LCG candidates for which we have Keck/LRIS spectra for, zf 9775, zf 11754, and zf 14000, are the spectroscopic candidates that we categorize as Group A. We then categorize the remaining sample of 32 sources as "photometric" LCG candidates. We define a redshift confidence parameter (z conf ) to categorize the LCG candidates, with z conf = 1 being most confident and z conf = 4 being least confident.
• z conf = 1 (most confident): Has at least one highquality spectroscopic redshift and at least two redshifts total, both above the redshift for F336W or F435W to probe LyC flux (z lim ). This only applies to the three Group A candidates with LRIS spectra.
• z conf = 2: Has more than two photometric redshifts > z lim OR has one or more redshifts from a reliable source > z lim . We consider a reliable source a tentative public spectroscopic redshift that we have in-spected (i.e., zf 11423 see Section 4.2, and Figure B1), photometric redshifts from ZFOURGE or the CAN-DELS average minus one standard deviation. This is Group B, which contains 10 LCG candidates.
• z conf = 3: Has one photometric redshift > z lim from a less reliable source (i.e., COSMOS2015, that we found to deviate more from the other available redshifts) and no other redshifts (i.e., not enough information to rule the redshift reliable or not). We place these candidates in Group C.
• z conf = 4 (least confident): Has one photometric redshift > z lim from a less reliable source (COS-MOS2015) and other conflicting redshifts that place the source at < z lim . We place these candidates in Group C along with z conf = 3 sources, with a total of 22 LCG candidates.
In Figure 4 we show 3 postage stamps of the three Group A spectroscopic LCG candidates (with z conf = 1), ordered by magnitude in F814W (m F814W ) from brightest to faintest. In Figure 5 we show the 10 Group B photometric LCG candidates ordered by m F814W (with z conf = 2; that includes the tentative spectro-  Table 5 for a full summary of all available redshift information and Table 6 for the photometry of all LCG candidates. See section 5.1 for details of their categorization and parameters. scopic redshift of zf 11423, see Section 4.2 and Figure  B1). In Figure 6 we show our 22 Group C photometric LCG candidates (with z conf = 3 and 4) also ordered by m F814W .
For each of the LCG candidate figures we show the available HST bands for the targets (new F336W and F435W, CANDELS-F606W and COSMOS-F814W), the NEQ combined image, and the source-detection segmentation map of the NEQ image (NEQ-Seg) used for isophotal fluxes (right stamp). We apply minor smoothing to the blue bands by convolving the images with a 2D Gaussian kernel (1.5-pixel standard deviation) to highlight the often-faint flux. We indicate which of the blue bands contain clean LyC flux if their > z lim redshifts are accurate. Above the row of stamps for each candidate, we show their IDs, any > z lim redshifts, S/N of their LyC flux, z conf , and m F814W values. We over plot a 1.2 aperture for scale. See Figure 2 for the spatial distribution of the candidates split by Groups A (with IDs shown), B, and C. See Table 5 for a summary of all the LCG candidate details and redshift information, and Table 6 for their photometry.

Color-color plots
Color selections are often used to identify high-redshift galaxies and LBGs (e.g. Steidel et al. 1996). We therefore plot our LCG candidates on color-color plots with the available bands to show their colors. However, these filter combinations are not the normal bands used for high-redshift LBG selection, nor are they optimal for it. Given that we do not select any galaxies based on colors, we use these purely for illustrative purposes. Figure 7 shows the LCG candidates split by category: Group A (red stars), Group B (blue circles), and Group C (dark blue diamonds). We show the Parent Sample of all galaxies > z lim (black points) and the All-z Sample within our HST footprints (gray points). The evolutionary tracks for four Shapley et al. (2003) LBG composites with increasing Lyα strength (red, yellow, green, turquoise) redshifted from z = 2.7 to 5 are shown. We show LBGs with no LyC flux (dotted curves) and increasing LyC flux as observed from Earth (rest < 912Å) as compared to observed rest ∼ 1500Å flux (i.e., R obs =1%, 2%, 5%, 10%, 20%; see Cooke et al. 2014) that peel the tracks off at different colors from the top. Note that these tracks are not at derived f esc values but observed flux ratios (see Section 5.4 for more details). The tracks are marked at increasing redshift (left to right) with '+'s at z = 3.5, 4 and 5.
We show the color-color plots for the HST bands (Figure 7 top row) with colors calculated from our isophotal magnitudes. The isophotal magnitudes are extracted within the same footprints on each band for each respective source (see Figures 4 to 6). The F336W iso −F606W iso vs. F606W iso −F814W iso plot (top left) shows all but one of our LCG candidates selected within the F336W band. The missing candidate is zf 13676 from Group B (a lower-priority candidate) as it is not covered by F606W. The F435W iso −F606W iso vs. F606W iso −F814W iso plot (top right) is missing one of the top Group A candidates (zf 11754) as it does not have F606W coverage. We therefore use our derived correction (see Section 2.4.4, Equation 1) to directly compare our HSC PSF-and aperture-matched F435W fluxes with ground-based HSC bands (bottom plot). The F435W 1 ap −HSC r 1 ap vs. HSC r 1 ap −HSC i 1 ap shows zf 11754 (red star). Some LCG candidates (mostly in Group C) are not shown on the axes ranges we present. This may indicate contamination of the lower confidence LCG candidates if they are that far away from the All-z cloud (gray points).
LBG selection boxes for z = 2-4 galaxies (higherredshift LBGs are off the plot) with these filter combinations are shown in gray in Figure 7. The boxes are comprehensive to span possible LBG regions. However, they are purely illustrative and should not be used for the clean selection of LBGs. We start with the boxes for a similar filter combination presented in Alavi et al. (2016) (F336W−F435W vs. F435W−F814W). The box is defined by the location of color-color evolution tracks of star-forming galaxies for these filters (not shown here, see Alavi et al. 2016). The star-forming galaxy tracks are model predictions from Bruzual & Charlot (2003) synthetic stellar populations with constant star formation, 0.2 Z , 100 Myr ages and color excesses E(B−V)=[0,0.1,0.2,0.3]. We apply the Madau (1995) cosmic opacity prescription and dust extinction law from Calzetti et al. (2000). We translate these model-based LBG boxes in color space to our filter combinations. We then expand them to encompass the Shapley et al. (2003) tracks (dotted lines shown here) as needed. Most of the LCG candidates, including all the top Group A targets, sit outside of the comprehensive illustrative LBG boxes.

Escape fraction estimates
In determining f esc , the escape fraction of ionizing UV flux relative to the non-ionizing continuum, many assumptions are required that have significant degeneracies. These include the SED of a galaxy, its starformation history (SFH), metallicity, and assumed line- Not all the targets have F606W available (see Figure 2), so we also use ground-based HSC bands (see Figure 1 for a comparison and Section 2.4.4). We show the Group A spectroscopic candidates (red stars): top left plot has zf 9775 (upper point) and zf 14000 (lower point), the bottom plot has zf 11754. We also show the photometric candidates in Group B (blue circles), and Group C (dark blue diamonds). For reference, we show the Parent Sample at a common redshift to the candidates (black points) and All-z Sample within the same HST footprints (gray points). Overlaid are evolutionary tracks of four LBG composites from Shapley et al. (2003) with increasing Lyα strength (red, yellow, green, turquoise) redshifted from z = 2.7 to 5. We show LBGs with no observed LyC flux (dotted curves) and increasing observed LyC to UV continuum flux ratio (R obs =1%, 2%, 5%, 10%, and 20%) that peel the tracks off at different levels from the top respectively. The tracks are marked at increasing redshift (left to right) with '+'s at z = 3.5, 4 and 5. Selection boxes for LBGs at z = 2-4 galaxies (higher-redshift LBGs are off the plot) with these filter combinations are shown in gray. Most of the LCG candidates, including all the top Group A targets, sit outside of the illustrative LBG boxes. See Section 5.3.

of-sight IGM distribution through which it is observed.
There are often a variety of escape-fraction related parameters quoted in the literature which can further complicate the picture.
The simplest observational parameter for measuring LyC emission is the flux ratio of the observed LyC to observed UV continuum (R obs (λ); see Cooke et al. 2014). We measure R obs = 4%, 3% and 9% for zf 9775, zf 11754 and zf 14000 respectively (see Table 4 for a summary). These values assume the LyC flux to be in F336W (for zf 9775 and zf 14000) or F435W (zf 11754), and the UV continuum to be F814W for the redshifts (z ∼ 3.7-4.4) of our sources (see Figure 1). The R obs values vary slightly from the color-color plot approximations (where they overlap with the tracks) given the different filter combinations and galaxy properties.
The parameter required to constrain cosmic reionization is the absolute f esc , or f abs esc . This is the fraction of escaping LyC photons, produced by O and B stars in star-forming galaxies, that reach the IGM without being absorbed by the interstellar or circumgalactic medium (CGM; e.g. Here (F LyC /F 1500 ) obs is the observed rest-frame LyC to UV flux ratio, (L LyC /L 1500 ) int is the ratio of the intrinsic ionizing to non-ionizing luminosity density, and τ LyC IGM is the redshift-dependent IGM attenuation along the line of sight. Both (L LyC /L 1500 ) int and τ LyC IGM are highly model and assumption dependent. This f rel esc parameter can be converted to an f abs esc estimate (proposed in Inoue et al. 2005;Siana et al. 2007) using where k λ is the reddening law (Calzetti et al. 2000) and E(B − V ) is the total dust attenuation (see e.g., Meštrić et al. 2020, for a summary). We determine f abs esc for the three Group A candidates using the equations above as is traditionally done in the literature. We first determine average IGM transmission of the LyC probing filter using the commonly adopted models from Inoue et al. (2014) (T LyC IGM,I14 ). These models consider only the IGM contributions along the line of sight and do not include attenuation by the CGM. We find average T LyC IGM,I14 = 0.0118, 0.0261, and 0.0429 for zf 9775, zf 11754, and zf 14000 respectively. The intrinsic ratio between the luminosity in the LyC band over the non-ionizing band is measured for the best fitting Binary Population and Spectral Synthesis (BPASSv2.1; Eldridge et al. 2017) model for each candidate. We determine f rel esc using Equation 2 and since our best fitting E(B − V ) = 0, f abs esc = f rel esc , so no further correction for dust is applied. With the attenuation values from Inoue et al. (2014), we find escape fractions of f abs esc,I14 = 8.75, 3.63, and 5.87 for zf 9775, zf 11754, and zf 14000 respectively.
We then calculate a new average IGM transmission of the LyC probing filter (T LyC IGM ) by generating 10,000 sight lines that include contributions from both the IGM and CGM (see Appendix C for more details). We find average T LyC IGM = 0.0015, 0.0141, and 0.0125 for zf 9775, zf 11754, and zf 14000 respectively. Using the same calculation as before with these new attenuation estimates, we find f abs esc = 111.3, 9.8, and 20.23 for zf 9775, zf 11754, and zf 14000 respectively. Note that these are not percentages but fractions so are well over 100% (i.e., 980% for zf 11754) and are not physical given these assumptions.
To further investigate f esc of the Group A candidates, we also adopt a statistical approach to generate realistic f esc probability distribution functions (PDFs), f PDF esc , that account for some of the sources of uncertainty. We use the probabilistic method described in Bassett et al. (2022) that we outline in more detail in Appendix C. Probabilistic approaches to estimating f esc have also been adopted in previous studies Vanzella et al. 2016). We show the f esc PDFs for one LCG candidate, zf 11754, in Figure 8 (left panel). The PDFs are scaled by the maximum probability (max(P)) assuming the model flux is equal to the observed flux. The PDFs are shown at each of the 25 BPASS model ages and for the best fitting model (zem5 = Z * = 1e −5 , log(age) = 6.0) and the weighted average at this fixed metallicity (red line). We also show the probabilityweighted average f esc PDF at each metallicity (Figure 8, middle panel). The red line shows the weighted average for all 275 BPASS metallicities and ages. See Figure C1 for grids of weights for each BPASS model for zf 11754.
The f esc estimates are constrained from 0 to 1, and the maximum is at 1. These plots show that the most consistent f PDF esc for zf 11754 is 100%. The right panel of Figure 8 shows a histogram of fluxes (blue) across all 10,000 sight lines for the redshift of zf 11754 (z = 4.470) and 10,000 f esc values between 0 and 1 (10 8 trials). The observational data for zf 11754 are represented by a Gaussian (black curve) with mean of its observed LyC flux in F435W (µ) and width of the flux error (σ). Where the blue histogram tails off on the right is f PDF esc = 100% and high transmission. The peak for zf 11754, i.e., its  highest probability f PDF esc , sits even further right of that tapered tail. The right panel therefore shows that the actual probability of this estimate being realistic is low and that it is likely > 100%. See Table 4 for summary of all f esc values and related parameters.
From the full range of f esc and IGM transmission combinations tested, ∼ 1/10, 000 sight lines would produce a high enough transmission to obtain the observed LyC flux provided f esc = 1. Given the very low likelihood of this happening, we deem that the more likely explanation is that the f esc > 1. If adopting a similar approach to Inoue et al. (2014) where the CGM is not included in the sight lines, the likelihood of this high-transmission and high-escape combination only slightly increases and the models become less physical. We are unable to provide absolute probabilities of these values being possible since they would depend on the intrinsic PDF of f esc for a population of star-forming galaxies at the redshifts of the Group A candidates (z = 3.7-4.4), which is unknown. All of the probabilities we present are relative, so we therefore quote only the f esc value that has the highest relative probability within the framework of our analysis.
Quantifying the errors is not possible as the PDF analysis indicates the most likely true value of f PDF esc is > 100%, and therefore not constrained. The main limitation here is that even in the 99 th percentile, IGM transmission is ∼ 0.2 in the F435W band at the redshift of zf 11754 (see Appendix C, Figure C2). Even for the most transparent sight lines with the largest intrinsic LyC fluxes and f esc = 1 as inputs, the resulting flux of the model is > 3 times the measured uncertainty below the measured flux (i.e., they are too faint). This is also the case for the other two spectroscopic candidates zf 9775 and zf 14000, so we only present the plots for zf 11754 to illustrate the method. An additional caveat for zf 14000 is that from the full array of SEDs, there are no models that fit the data well, whereas the bestfit models to zf 11754 and zf 9775 are realistic. We therefore have additional doubts about the validity of the f PDF esc measurement, and therefore LyC detection, of zf 14000.

Escape fraction comparison
The most likely value of f abs esc and f PDF esc for all three Group A spectroscopic candidates is > 100%, and therefore not physical with these models and assumptions. It is possible that some assumptions that go into the f PDF esc calculations may not be appropriate for our candidates. Many of the assumptions are those commonly adopted by others for the more traditional f esc estimates (e.g., f abs esc calculated with Equation 3). The models (BPASSv2. 1;Eldridge et al. 2017), metallicity ranges (Steidel et al. 2018;Fletcher et al. 2019;Bassett et al. 2021a), and IGM sight line calculations (Inoue et al. 2014;Steidel et al. 2018;Bassett et al. 2021a) are all based on established methods but do not include strong emission lines in the SEDs, only exponentially declining SFHs. Although as the best fitting model for all three Group A candidates is "zero-age" for BPASS (equivalent to a single, instantaneous burst with an age of 10 6 years), the choice of SFH will not affect the fact that f PDF esc is unconstrained at > 100% (see Appendix C), The statistical f PDF esc estimate approach presented here has been tested on a spectroscopically confirmed z = 3.64 galaxy (VUDS 511227001;Le Fèvre et al. 2015) in Bassett et al. (2022). CLAUDS (Sawicki et al. 2019) U -band shows a LyC detection for this galaxy. The resulting PDF has a peak within the 0 to 1 range, with f PDF esc = 0.45 +0.39 −0.28 . If we derive many more possible sight lines (10 7 rather than 10 4 ), it is possible that at least one would be transparent enough to match the observed flux for our candidates. However, all three would be unlikely to be within our four untargeted HST pointings.
The true form of f PDF esc could also vary by several factors, including galaxy mass, which further complicates its accurate estimation (e.g., Finkelstein et al. 2019;Naidu et al. 2020;Bassett et al. 2021a). The broad range of f esc estimates in the literature, and even for each individual source, makes it hard to infer meaningful interpretations for the ionizing flux contributions from highredshift galaxies. There only exists a relatively small sample of individually detected LyC candidates and they make up a diverse population of chance detections.
A key conclusion from Meštrić et al. (2020) was that only the cleanest lines of sight could work for observing LCGs. Therefore, these detections tell us as much about the random line of sight of the galaxies as it does about the galaxies themselves. It is likely that a population of high-redshift (z ∼ 3-5) galaxies may hold some of the most valuable clues about the sources that reionized the universe due to their proximity in time, and likely their physical properties. However, given the severe attenuation of LyC photons by the intervening IGM, this population remains incredibly elusive.
Only in large statistical numbers can we average over possible IGM sight lines. Even then, many of those that we observe in the LyC may exist along the cleanest lines of sight. Non-detections of LyC may be the best path forward for constraining the global ionizing photon budget. The KLCS clean non-detection (< 3σ) sample of 107 galaxies has a f abs esc = 0.05 ± 0.01 (Steidel et al. 2018;Pahl et al. 2021). Bian & Fan (2020) stack 54 Lyα emitters at z 3.1 and find no visible LyC but 3σ upper limits of f abs esc < 14-32%. LBG stacks from Ji et al. (2020) give upper limits of f abs esc 0.63%. The spectroscopic assessment in Cooke et al. (2014) gives f abs esc upper limits of 0.086 for z ∼ 3 galaxies. Using a simplified stacking method of f abs esc for galaxies without AGN at z ∼ 2.3-4.3, Smith et al. (2018Smith et al. ( , 2020 find limits of 0.08, 0.06, 0.57 for F225W, F275W and F336W respectively from their stacks. Meštrić et al. (2021) analyze quoted upper limits in the literature and find global f abs esc values of 0.086 at 2 z ≤ 3, 0.105 at 3 ≤ z ≤ 4 and 0.084 at 2 z < 6.

Do the LCG candidates have line-of-sight low-redshift interlopers?
Even with deep HST LyC imaging and secure LRIS redshifts for zf 9775 and zf 11754, the LCGs presented in this paper remain candidates. This is in part due to the quality of the spectra given consistently poor weather conditions for our observations (see Table 3) and slit positioning on our targets. Nevertheless, from the spectra and imaging alone, these appear to be promising candidates. However, the f abs esc values and statistical f PDF esc estimates (that explore a broader range of possible assumptions than are usually adopted), we find that the most likely f esc is > 100% and therefore not physical given these models and assumptions. The candidates shown on the color-color plots ( Figure  7) mostly sit outside the comprehensive LBG region. This was hypothesized by Cooke et al. (2014) and found with a small number of LCG candidates by Meštrić et al. (2020). The Group A candidates sit on the LBG tracks with observed LyC flux artificially added (R obs ∼ 5% for zf 9775, R obs ∼ 10% for zf 11754 and zf 14000, see Table 4 for calculated values). This could imply that these LCGs are similar to LBGs that exhibit LyC flux (as shown in KLCS; Steidel et al. 2018;Pahl et al. 2021) or they have flux contamination.
Are the Group A spectroscopic LCG candidates contaminated by low-redshift interlopers? Spatially resolving the LyC in high-resolution imaging is one of the best ways to identify low-redshift contaminants (e.g., Vanzella et al. 2012). We find that the chance of contamination is higher than the LCG candidate detection rates of sources in the field. We derive 1.3% to 5.0% probabilities of 3σ overlapping sources in our deep HST F336W and F435W images probing the LyC (see Section 2.4.3, Table 1). The redshift-independent contamination rates are higher than the detection rates relative to the Allz Sample from the combined redshift catalog (0.85% in F336W, 0.14% in F435W). Therefore, it is possible that our candidates are high-redshift sources contaminated by low-redshift interlopers and our selection method is biased towards finding these chance alignments.
Our blue HST images (∼ 29-30 mag depths at 5σ) are deeper than other searches for LyC emitters. Comparable surveys from Siana et al. (2015) reach ∼ 29.27 mag at 5σ and Smith et al. (2018Smith et al. ( , 2020 reach ∼ 29.11 mag at 5σ in F336W (with data from Oesch et al. 2018). Deeper bluer filter images could help to resolve the question of possible contamination. However, with deeper images, picking up faint foreground sources is more likely (e.g., Alavi et al. 2014). At the high redshift of our candidates (z ∼ 3.7-4.4) the sight line is long and therefore the chance of coincident sources is higher, compared with z ∼ 1 (e.g., Siana et al. 2007Siana et al. , 2010.
Higher signal-to-noise spectra with maximal spatial separation of components along the slit of zf 9775 might help us identify the true nature of the candidate. If detected, emission lines in the LyC part of the spectrum could confirm this double source as a low-redshift interloper rather than a physical merger. The LyC emission for zf 11754 is also slightly offset from the optical emission and could also be contaminated. zf 14000 is very centered but this is the Group A candidate with the largest uncertainties with respect to its nature, redshift, and f esc estimates. However, it may be that LyC flux does not trace non-ionizing flux, for example mergers can disrupt gas in galaxies, inhomogeneously exposing star-forming regions.
At higher redshifts, the IGM opacity increases significantly. This further reduces the chance of identifying strong LyC detection, particularly at the redshifts of our Group A spectroscopic candidates (z ∼ 3.7-4.4). For zf 9775 and zf 14000, the band cleanly sampling their LyC flux is F336W (zf 9775 also has F435W which contains some forest lines), which would be around ∼ 560-695Å. The expected LyC flux is low in F336W, e.g., intrinsic: F int = 0.034 µJy, attenuated: F att = 5.2e −5 ±3.2e −4 µJy, ≈ 34.6±6.6 mag. Although, Ion3 does show faint LyC down to 750Å . The results presented in this work are a cautionary tale that even with a blue spectroscopic redshift, detection in a "LyC" probing band, and HST imaging to remove contamination, the true nature of the sources may still not be clear.

Future prospects & need for spectroscopic follow-up
The Group A spectroscopic candidates are the three brightest galaxies of all our LCG candidates, with isophotal m F814W < 25 mag. Given the poor to moderate observing conditions (with poorer-than-optimal seeing), obtaining spectra for galaxies fainter than m F814W ∼ 25 proved challenging. Of the three top candidates, additional improved spectra are required to confirm them as genuine LCGs. A slit orientation to maximally spatially resolve the two components of zf 9775 is needed. A slit for zf 11754 placed to avoid its bright neighboring sources is also required. We have already attempted to target these three Group A candidates with these optimized slit positions (2020B N168, PI Prichard, weathered out). We were recently re-awarded this time for January 2022 (2021B N87, PI Prichard).
The remaining Group B and Group C photometric candidates are relatively faint with isophotal optical magnitudes ( 26 mag, down to 28.3 mag). Note that the isophotal magnitudes are a conservative clean flux estimate (best for accurate colors) that may not include all the flux of a source. In good conditions, with seeing < 1 , we expect to get redshifts in 4-6 hours down to m F814W ∼ 26. If a candidate has a strong line, assumed to be Lyα in emission, a redshift can be determined for a source fainter than m F814W ∼ 26.
High S/N spectroscopy is essential for identifying genuine LCGs. In a recent complimentary LCG study with MOSFIRE, Bassett et al. (2022) found that the ZFOURGE photometric redshifts for all seven LCG candidates were overestimated by δz ∼ 0.2, most likely due to template mismatch. These candidates were therefore ruled out as LCGs as their true spectroscopic redshift placed them at a redshift such that the U -band was not cleanly sampling the LyC. Although the ZFOURGE stated accuracy is ∼ 2% errors , Bassett et al. (2022) concluded that the search for LCGs is biased towards finding overestimated redshifts (with 94% confidence). The results of this paper further emphasize that HST data are required for removing lowredshift interlopers, and demonstrate the requirement of blue spectra to aid with the removal of contaminants.
Redder wavelengths can be used to confirm redshifts, but may miss blue features to aid with contaminant removal. We have had supporting programs at redder wavelengths using Keck/MOSFIRE (2018B W151, PI Bassett) that focused on the [OIII]/[OII] relation with LyC (Bassett et al. 2019(Bassett et al. , 2022. The combination of MOS observations spanning the bluest wavelengths does provide the best opportunity for getting redshifts and removing low redshift contaminants efficiently, along with the opportunity to spectroscopically detect the LyC in optimal observing conditions if bright enough. Keck/LRIS is not the only blue MOS option we have tried, we were awarded telescope time with the Large Binocular Telescope/MODS (LBT) that was also weathered out (see Table 3 for a summary). Integral field units (IFUs) could be valuable for exploring spectroscopic properties of individual LyC emission regions and for removing contaminants. However, this requires high pixel resolution (< 1 ) for these small sources and is less efficient for building up statistical samples.
Additional information in the blue part of the spectrum can help reveal the nature of the detection and if it is contaminated. This expensive spectroscopic campaign will be handled with superior efficiency for galaxies such as our candidates with the new class of blue MOS instruments on 30 m telescopes. Additionally, the proposed Keck Wide-Field Imager (KWFI; Gillingham et al. 2020) would provide photometry down to ∼ 29 mag in u in just ∼ 5 hours and ∼ 3 hours in g over an area ∼ 150× that of HST to rapidly identify photometric candidates for spectroscopic follow-up. This combination could efficiently detect and confirm statistical samples of LCGs that could help to constrain valuable cosmic reionization parameters.

SUMMARY & CONCLUSIONS
We present new deep WFC3/F336W and ACS/F435W HST images (with depths ∼ 29-30 mag at 5σ), that we use to identify a population of highredshift LCG candidates. We process the WFC3/UVIS images with our new custom pipeline that reduces noise, gradients, hot pixels and offsets in the images to produce the cleanest HST UV data possible. Some of our custom reduction steps for WFC3/UVIS have now been adopted by the WFC3 team and are available for download directly from MAST for all existing data. We apply additional corrections to the F336W images that are not included on MAST but are easily applied to data with a new public reduction notebook 13 .
We collate a comprehensive Parent Sample of possible candidates without relying on the Lyman-break technique. We use photometric redshifts from the ZFOURGE survey, COSMOS, CANDELS, and publicly available spectroscopic redshifts from 16 surveys within our new HST footprints. Although no catalog will be free from selection bias, this is one of the most comprehensive and agnostic LCG searches possible with the available data. There are 3065 objects at all redshifts in the two F336W footprints and 6430 in the three F435W footprints (our All-z Sample). 380 sources are at z > 3.087 in F336W and 242 at z > 4.391 in F435W, the redshift limits (z lim ) for the bands to cleanly probe LyC (Parent Sample). We visually inspect all sources in the Parent Sample multiple times to obtain a long list of possible candidates (77) for spectroscopic confirmation.
The bulk of awarded spectroscopic observing time was weathered out or had less-than-optimal conditions. We present the spectra obtained for some of our sources with Keck/LRIS (∼ 1.5 nights). We acquire spectroscopic redshifts for the three brightest LCG candidates (i.e., Group A spectroscopic candidates): zf 9775 at z = 4.365, zf 11754 at z = 4.470, and zf 14000 at z = 3.727. The list of candidates is cleaned and we select only those with > 3σ LyC detections. The remaining 32 sources are photometric LCG candidates that we split by redshift confidence. Group B has 10 LCG candidates with relatively confident photometric (and one low-quality public spectroscopic) redshift. Group C has 22 candidates with more unreliable redshifts. We summarize the main results of this work below.
• We find a 6.84% detection rate of our > 3σ LCG candidates in the F336W footprints and 3.72% in the F435W footprints relative to the Parent Sample. We find number densities of ∼ 1.8 arcmin −2 at z > 3.087 in the F336W pointings and ∼ 0.1-0.4 arcmin −2 at z > 4.391 in the F435W pointings (at varying depths).
• We determine conservative redshift-independent LyC contamination probabilities of 1.8%-3.3% depending on the field. These are higher than the LCG candidate to All-z Sample detection rates (0.85% for F336W, 0.14% for F435W). Most of the candidates are in Group C, whose redshifts we are least confident in and that are likely to be dominated by low-redshift interlopers.
• The HST images and Keck/LRIS spectroscopy for the three Group A spectroscopic candidates imply that these detections could be genuine. We determine f esc using the traditional (f abs esc ) and a statistical (f PDF esc ) approach. The results imply that the most likely escape fraction estimates for all three spectroscopic candidates for both methods is > 100%, and therefore not physical given the models and assumptions used.
• From the purely illustrative color-color plots, most candidates sit outside of the LBG box. They lie on tracks with observed LyC to UV continuum flux ratios (R obs ). Calculated directly from their LyC (F336W or F435W) to UV continuum (F814W) bands are R obs = 4%, 3% and 9% for zf 9775, zf 11754 and zf 14000 respectively. This could imply that the LCG candidates are similar to LBG galaxies with LyC emission, or that they have flux contamination.
• The three Group A spectroscopic LCG candidates and 32 Group B and C photometric LCG candidates require spectroscopic follow up to confirm their true nature. Even with spectroscopic redshifts, UV-HST imaging and multi-band photometry, more information is required to confirm these as genuine LCGs. The future KWFI instrument will be valuable for efficiently identifying new faint LCG candidates. The new class of 30 m telescopes will be able to spectroscopically confirm the faintest candidates presented here so that larger statistical samples can be built up for constraining cosmic reionization parameters.

ACKNOWLEDGMENTS
We thank the anonymous reviewer for a helpful report. This research is based on observations made with the NASA/ESA HST obtained from the STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program ID 15100 (PI Cooke) and we acknowledge financial support for this work from HST. We also acknowledge some additional support from HST program ID 15647 (PI Teplitz) for work related to the final data calibrations presented here. This work is also based on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST. Part of this research was funded by the Australian Research Council Centre of Excellence for All-sky Astrophysics in 3 Dimensions (ASTRO-3D), CE170100013, the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), CE110001020, and the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), CE170100004.
Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to NASA through the agency's scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.  This paper is based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at National Astronomical Observatory of Japan. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. We make several improvements to the WFC3/UVIS reduction pipeline (Bourque & Baggett 2016) for which we provide details here. The WFC3/UVIS channel has a charge coupled device (CCD) detector. Photons landing on the array are converted to charge which is measured by being transferred across the array to the "readout" amplifiers. Due to degradation of the detector over time from radiation, charge traps develop, reducing the charge transfer efficiency (CTE). This effect is corrected after the fact (e.g., Anderson 2020; Anderson et al. 2021) to recover as much of the original flux as possible.
The detector also has imperfections such as hot pixels. The WFC3/UVIS CCD undergoes anneals (brief heating of the detector every ∼month) to reduce the number of hot pixels. In each anneal cycle (period between anneals), the detector has a different "fingerprint". Darks and data are typically processed as soon as possible after being observed. This speed means that darks from the previous anneal cycle (using a different detector fingerprint than the observations) are used to calibrate the data which diminishes the quality of the previously available data products.
The new CTE code (calwf3 v3.6.0; Anderson 2020; Anderson et al. 2021) optimizes parameters for the data to mitigate read noise, which reduces background noise and gradients across the images. Figure A1 (top left panel) shows the RMS measured in 1000 randomly placed apertures on sky (avoiding sources) of a drizzled image processed with the old CTE code (red) and the new code (blue). The background noise is significantly reduced with the new code. The CTE corrections are applied to both the science data and the darks in the reduction. Figure A1 (top right panel) shows the signal gradients across one chip of a single FLC. We show the median signal of each row (dots) and the running median of the row values (lines) for the old (red) and new (blue) CTE code. The new CTE correction results in much flatter signal across the chips. The new CTE code along with the use of concurrent darks (from the same anneal cycle as the data are observed) are now included as standard by WFC3, and are available from MAST (as of May 2021).
We developed a novel flagging routine to identify hot pixels more effectively. The method used by MAST adopts a constant threshold, and due to the imperfect CTE correction, increasingly more hot pixels are missed further from the readout. To correct for this loss of hot pixels caused by CTE degradation, we derive a variable threshold for every chip to identify roughly a constant number of hot pixels across the chip as expected. The most accurate number of hot pixels is closest to the readout with the shortest distance to travel across the array. We therefore take this value to be the true number of hot pixels, and find the threshold that selects just below this value for each set of 50 rows. We then fit a function to those values to derive a variable threshold conversion for each chip. Figure A1 (bottom panel) shows the number of hot pixels as a function of row number with the readout locations at the edges of the plot. We show the hot pixels identified using a constant threshold (red) and the new variable threshold method (blue) for the old (faint lines) and new CTE code (bold lines). The new CTE code increases the number of hot pixels missed as the flux pedestal subtracted is larger, so more pixels are below the previous constant threshold. Figure A1. Top left: Comparison of RMS measured from images processed with the old CTE code (red) vs. the new CTE code (blue). The RMS of the background is measured from 1000 random apertures on sky (avoiding sources). The new CTE code significantly reduces the background with an average of 3.3 × 10 −4 e/s from 3.9 × 10 −4 e/s, valuable for the detection of faint sources. Top right: Signal gradients across one chip of an individual FLC. The median of each row (points) and the running medians (solid lines) are shown for the old (red) and new (blue) CTE codes. The new CTE reduces the gradients over the chips. Bottom: Number of hot pixels identified as a function of row number across chip 2 (left) and chip 1 (right). The previous method using a constant threshold (red lines) misses ∼30% of hot pixels. Using a new variable threshold (blue lines) identifies the correct number of hot pixels across each entire chip. Bold and faint lines show the new and old CTE corrections respectively.
The constant threshold misses around ∼ 30% of the expected hot pixels. For a detailed step-by-step run through of all the steps we use for the WFC3/UVIS reduction, see our pipeline and supplementary codes on GitHub 14 .
We flag negative divots adjacent to the readout cosmic rays (ROCRs) in the FLC data quality (DQ) arrays to prevent them from propagating through to the drizzled images. As these cosmic rays appear during readout, they fall closer to the amplifier than they appear in the FLCs, which causes over-corrections by the CTE code. Due to their nature, these flags increase with distance from the amplifier while the array is being readout, and affect a total of Figure A2. Comparison of a reduced CTE-corrected single-visit exposure (FLC) using the old and new reduction methods. An FLC previously available on MAST is shown on the left. On the right is an FLC reduced with the new darks (including the new hot pixel method) and new CTE (that reduces noise and gradients). An additional correction to equalize the bias offsets between the four quadrants is also applied (medsub).
∼ 0.2% detector pixels. We flag divots within five pixels of a CR hit (away from the readout direction) that are 3σ below the image mean as bad detector pixels (value of 4 in the DQ arrays).
We equalize the overall signal levels in the four amplifier quadrants of each FLC. Typically, small bias offsets exceeding one electron can exist between these regions, causing discontinuities in the background. To remove these differences, we calculate the offsets between the median level in each individual quadrant and the average of all four quadrants, multiply the difference by the flat field, and subtract these offsets from the original data. As the zeropoints are determined as an average over all four quadrants, the equalization of the amps does not affect our photometry as we have a roughly even distribution of sources in the images. Figure A2 shows a comparison of an FLC previously available on MAST (left) and one processed using the new concurrent darks, the new CTE, and the amplifier bias offsets equalized (medsub). Gradients, noise, offsets, and blotchy patterns are greatly reduced using the new reduction methods. The new CTE and concurrent darks are now applied as standard by STScI and available from MAST. The additional corrections can be easily applied to reduced FLCs from MAST using a new publicly available notebook 15 .

B. ANCILLARY SPECTRA
We show the DEIMOS 10K (DM18; Hasinger et al. 2018) spectra for two of the LCG candidates in Figure B1. Both the unsmoothed spectra (pale background spectra) and 15-pixel boxcar-smoothed spectra (dark green/black) are shown overlaid. We show the DM18 spectrum of a Group A spectroscopic candidate zf 14000 (DM18 ID: L677152, z DM18 = 3.779; top plot). We shift the spectrum to our LRIS redshift (green) and the DM18 redshift (black) to compare. We highlight one potentially real spectral feature in red (C IV) that may add weight to the LRIS redshift, although with low confidence.
The lower panel of Figure B1 shows the DM18 spectrum of one of the Group B photometric candidates zf 11423 (DM18 ID: L658665, z DM18 = 3.655). We highlight two potentially real spectral features in red. However, as the spectrum is noisy and is shown heavily smoothed here, we are not confident of the DM18 redshift. This galaxy therefore remains in Group B with the same confidence as our photometric candidates. Both spectra have assigned quality flags from DM18 of 1 (lowest quality), and a quality flag in the combined spectroscopic catalog of 2 (from 1 uncertain to 4 best) based on a comparison of available redshifts in each of the 16 catalogs. Figure B1. DEIMOS 10K spectra (DM18; Hasinger et al. 2018) of two LCG candidates. We show the unsmoothed spectra (pale background spectra), 15-pixel boxcar-smoothed spectra (dark green/black lines), LBG composite (blue) from Shapley et al. (2003), and spectral features (vertical lines). Top: DM18 spectrum of a Group A spectroscopic candidate, zf 14000 (DM18 ID: L677152). The redshift we determine from the LRIS spectrum ( Figure 3) is z = 3.727. We show the DM18 spectrum shifted to the redshift from LRIS (green) and the spectrum at the DM18 redshift of zDM18 = 3.779 (black) for comparison. We highlight a potentially real feature which adds weight to the LRIS redshift (red). Bottom: DM18 spectrum of a Group B candidate, zf 11423 (DM18 ID: L658665, zDM18 = 3.655). We highlight two potentially real spectral features (red) that might add weight to the DM18 redshift but with low confidence. Figure B2. 3D-HST spectrum (Brammer et al. 2012;Momcheva et al. 2016) for zf 9775 (3D-HST ID: 9997) redshifted to the LRIS redshift (zLRIS = 4.365, green, top panel), and that quoted in 3D-HST (z3DHST = 0.446, black, bottom panel). The 3D-HST spectrum does not provide sufficient evidence for either redshift, we therefore deem the 3D-HST redshift invalid due to the confident redshift we derive from the LRIS spectrum ( Figure 3).
We show the unsmoothed 3D-HST G141 grism spectrum (Brammer et al. 2012;Momcheva et al. 2016) of the Group A spectroscopic candidate zf 9775 in Figure B2. The same spectrum is shown at rest wavelengths assuming the redshift derived from our LRIS spectrum (z LRIS = 4.365, green, top panel) and that quoted in 3D-HST (z 3DHST = 0.446, black, bottom panel). There are no redshift quality flags provided by the 3D-HST survey. However, as the zf 9775 3D-HST redshift is derived from a single < 1σ S III line, it is low confidence. The 3D-HST grism spectra are not visually inspected. A recent paper performed visual inspection of a large number of 3D-HST spectra and found that a quarter of the published redshifts are not correct (Henry et al. 2021). Additionally, redshifts for zf 9775 from ZFOURGE, COSMOS, average CANDELS, and LRIS are all in close agreement with a stated redshift of z ∼ 4.2-4.4. The 3D-HST redshift (z = 0.446) deviates from these values significantly.

C. ESCAPE FRACTION SIMULATIONS
We summarize the key steps for estimating f esc PDFs for the three Group A spectroscopic LCG candidates. This method is detailed in full in Bassett et al. (2022). See Section 5.4 for a summary of the results.
1. We use a range of possible SED templates from BPASS (v2. 1;Eldridge et al. 2017). We combine these with a single, exponentially declining SFH with e-folding time 0.1 Gyr. We also use BPASS models with metallicity Z * = 0.001-0.020 (as in Steidel et al. 2018;Fletcher et al. 2019;Bassett et al. 2021a), initial mass function slope α = −2.35, and stellar-mass limit 300 M . We consider only SED models with ages 2.5 × 10 8 years ((L 900 /L 1500 ) int > 0.05) as fainter LyC is highly unlikely to be detected in galaxies at z 3.7-4.4.

2.
We estimate E(B−V) at each BPASS model age. We apply each E(B−V) value to a Reddy et al. (2016) attenuation curve, determine photometric fluxes of attenuated models using the ZFOURGE filters and minimize χ 2 . With the optimal E(B−V) for a given age, we calculate the likelihood that the model represents our LCG candidates. This Figure C2. IGM transmission curves for zf 11754 (z = 4.470), the highest redshift spectroscopic LCG candidate in Group A. Shown is the median (black line) of the 10,000 IGM curves generated. The 5 th (red), 25 th (orange), 75 th (yellow), and 95 th (green) percentiles are also shown. The Lyman break (∼ 912Å) is shown (red dashed line). Overlaid is the LyC probing filter for this candidate F435W (turquoise). The curves generated for all three of the Group A spectroscopic candidates are very similar due to their close redshifts. The median and percentile IGM curves do not represent individual lines of sight.
just one distribution for zf 11754 in Figure C2. The median and percentile IGM curves do not represent individual sight lines.
5. We then determine the probability that the photometric LyC flux and error is consistent with f esc considering the SED models and IGM transmission curves. We apply our ensemble of 10,000 IGM curves and 10,000 f esc values (between 0 and 1) to every BPASS model. We measure the probability (P obs ) that a given mock observation is consistent with the observations. We take the value at a model flux for a Gaussian distribution with center of the observed flux and FWHM of the error.
6. The probability of any f esc value is taken as the median across all 10,000 sight lines. The PDFs of f esc are for an individual SED age, metallicity, and E(B−V). We therefore use our best fitting models from Step 2 to generate the most appropriate f esc PDF (and corresponding f PDF esc estimate) for each of our sources.