# The HDUV Survey: Six Lyman Continuum Emitter Candidates at $z \sim 2$ Revealed by HST UV Imaging

R. P. Naidu, P. A. Oesch, N. Reddy, B. Holden, C. C. Steidel, M. Montes, H.

Atek, R. J. Bouwens, C. M. Carollo, A. Cibinel, et al.

► To cite this version:

R. P. Naidu, P. A. Oesch, N. Reddy, B. Holden, C. C. Steidel, et al., The HDUV Survey: Six Lyman Continuum Emitter Candidates at $z \sim 2$ Revealed by HST UV Imaging. The Astrophysical Journal, 2017, 847, 10.3847/1538-4357/aa8863. insu-03747442

HAL Id: insu-03747442

https://insu.hal.science/insu-03747442v1

Submitted on 4 Sep 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The HDUV Survey: Six Lyman Continuum Emitter Candidates at $z \sim 2$ Revealed by HST UV Imaging*

R. P. Naidu<sup>1,2</sup>, P. A. Oesch<sup>2,3,4</sup>, N. Reddy<sup>5,14</sup>, B. Holden<sup>6</sup>, C. C. Steidel<sup>7</sup>, M. Montes<sup>2</sup>,

H. Atek<sup>8</sup>, R. J. Bouwens<sup>9</sup>, C. M. Carollo<sup>10</sup>, A. Cibinel<sup>11</sup>, G. D. Illingworth<sup>6</sup>, I. Labbe<sup>9</sup>, D. Magee<sup>6</sup>,

L. Morselli<sup>12</sup>, E. J. Nelson<sup>13</sup>, P. G. van Dokkum<sup>2,4</sup>, and S. Wilkins<sup>11</sup>

## Abstract

We present six galaxies at $z \sim 2$ that show evidence of Lyman continuum (LyC) emission based on the newly acquired UV imaging of the Hubble Deep UV legacy survey (HDUV) conducted with the WFC3/UVIS camera on the *Hubble Space Telescope (HST)*. At the redshift of these sources, the HDUV F275W images partially probe the ionizing continuum. By exploiting the *HST* multiwavelength data available in the HDUV/GOODS fields, models of the UV spectral energy distributions, and detailed Monte Carlo simulations of the intergalactic medium absorption, we estimate the absolute ionizing photon escape fractions of these galaxies to be very high—typically >60% (>13% for all sources at 90% likelihood). Our findings are in broad agreement with previous studies that found only a small fraction of galaxies with high escape fraction. These six galaxies compose the largest sample yet of LyC leaking candidates at $z \sim 2$ whose inferred LyC flux has been observed at *HST* resolution. While three of our six candidates show evidence of hosting an active galactic nucleus, two of these are heavily obscured and their LyC emission appears to originate from star-forming regions rather than the central nucleus. Extensive multiwavelength data in the GOODS fields, especially the near-IR grism spectra from the 3D-*HST* survey, enable us to study the candidates in detail and tentatively test some recently proposed indirect methods to probe LyC leakage. High-resolution spectroscopic follow-up of our candidates will help constrain such indirect methods, which are our only hope of studying $f_{esc}$ at $z \sim 5-9$ in the *JWST* era.

**Key words:** dark ages, reionization, first stars – galaxies: evolution – galaxies: high-redshift

## 1. Introduction

Identifying the sources that dominated cosmic reionization in the first 1 Gyr of cosmic time is still one of the key open questions of observational extragalactic cosmology. Recent advances in tracing the buildup of galaxies during the epoch of reionization (EoR) at $z > 6$ indicate that ultrafaint galaxies are very abundant in the early universe and that they dominate the UV luminosity density. This has led several authors to speculate that the faint galaxy population is the main driver for reionization, a scenario that can reconcile several independent measurements of the reionization history (e.g., Bouwens et al. 2006, 2012, 2015; Oesch et al. 2009; Ouchi et al. 2009; Bunker et al. 2010; McLure et al. 2010; Finkelstein et al. 2012; Grazian et al. 2012; Duncan & Conselice 2015; Robertson et al. 2015).

The main unknown in these studies is the fraction of ionizing photons that escape galaxies into the intergalactic medium (IGM), the so-called escape fraction, $f_{esc}$. This remains unconstrained observationally during the EoR. The typical conclusion of reionization calculations is that the escape fraction of galaxies has to be $\gtrsim 10\%$. Otherwise, their ionizing photon production falls short of the required value to complete reionization by galaxies, and other sources such as active galactic nuclei (AGNs) are needed to contribute significantly (e.g., Giallongo et al. 2015; Madau & Haardt 2015; Mitra et al. 2015, 2016; Feng et al. 2016; Price et al. 2016).

Direct observational constraints on $f_{esc}$ are effectively impossible to obtain at $z \gtrsim 4.5$ and into the EoR owing to the high opacity of the intervening IGM absorption (e.g., Prochaska et al. 2010; Inoue et al. 2014). However, at lower redshifts such direct studies of Lyman continuum (LyC) photons (at rest wavelength $\lambda < 912$ Å) are possible. Until recently, the few constraints on $f_{esc}$ that existed in the local universe were only upper limits indicating very low values of only a few percent at most (e.g., Leitherer et al. 1995; Deharveng et al. 2001; Grimes et al. 2009; Siana et al. 2010;

<footnote n="1"><sup>1</sup> Yale-NUS College, 12 College Avenue West, 138614, Singapore; rohan.naidu@u.yale-nus.edu.sg</footnote>

<footnote n="2"><sup>2</sup> Astronomy Department, Yale University, New Haven, CT 06511, USA</footnote>

<footnote n="3"><sup>3</sup> Geneva Observatory, Université de Genève, Chemin des Maillettes 51, 1290 Versoix, Switzerland</footnote>

<footnote n="4"><sup>4</sup> Yale Center for Astronomy and Astrophysics, Yale University, New Haven, CT 06511, USA</footnote>

<footnote n="5"><sup>5</sup> University of California, Riverside, 900 University Avenue, Riverside, CA 92507, USA</footnote>

<footnote n="6"><sup>6</sup> UCO/Lick Observatory, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="7"><sup>7</sup> Cahill Center for Astronomy and Astrophysics, California Institute of Technology, MS 249-17, Pasadena, CA 91125, USA</footnote>

<footnote n="8"><sup>8</sup> Institut d'astrophysique de Paris, 98bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="9"><sup>9</sup> Leiden Observatory, Leiden University, NL-2300 RA Leiden, The Netherlands</footnote>

<footnote n="10"><sup>10</sup> Institute for Astronomy, ETH Zurich, 8092 Zurich, Switzerland</footnote>

<footnote n="11"><sup>11</sup> Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH, UK</footnote>

<footnote n="12"><sup>12</sup> Excellence Cluster Universe, Boltzmannstr. 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="13"><sup>13</sup> Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85741 Garching bei München, Germany</footnote>

<footnote>*Received 2016 November 22; revised 2017 July 27; accepted 2017 August 22; published 2017 September 14*</footnote>

<footnote>* Based on observations made with the NASA/ESA *Hubble Space Telescope*, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.</footnote>

<footnote n="14"><sup>14</sup> Alfred P. Sloan Foundation Fellow.</footnote>
Rutkowski et al. 2016)—far too small compared to the ≳10% required for cosmic reionization. However, recent work with the COS spectrograph on the Hubble Space Telescope (HST) has identified a subsample of highly star-forming galaxies (SFGs) in the local universe that appear to show significant and detectable LyC emission (Borthakur et al. 2014; Izotov et al. 2016a, 2016b; Leitherer et al. 2016).

At higher redshift, the situation is similar. At $z \sim 2-3$, the LyC shifts into the observed ~2000–3500 Å range, allowing UV-sensitive instruments to directly detect ionizing photons. Early observations resulted in confusing results, with many of the direct detections being attributed to contamination by foreground sources (Vanzella et al. 2010, 2012; Nestor et al. 2011; Mostardi et al. 2015; Siana et al. 2015; Grazian et al. 2016). However, recently a small sample of three galaxies with confirmed direct detections of their LyC emission has emerged (Mostardi et al. 2015; Shapley et al. 2016; Vanzella et al. 2016). One of these is Ion2 at $z = 3.2$, which was originally identified in Vanzella et al. (2015, henceforth V15) using a method similar to the one we adopt in this paper. In particular, V15 simulated the UV flux of sources with secure spectroscopic redshifts to identify candidate LyC sources in the GOODS-S broadband imaging data. This resulted in two candidates, one of which, Ion2, has been confirmed as an LyC leaker through direct follow-up imaging (de Barros et al. 2016; Vanzella et al. 2016). Building up the sample size of such confirmed sources is crucial to aid our understanding of LyC photon escape from galaxies and of cosmic reionization.

In this paper we exploit the newly obtained UV imaging by the WFC3/UVIS camera on HST from the Hubble Deep UV (HDUV) imaging survey over the two GOODS/CANDELS-Deep fields (Oesch et al. 2016, submitted), along with data from the UVUDF survey (Teplitz et al. 2013; Rafelski et al. 2015). The HDUV filter set covers ~2500–3700 Å and directly images LyC photons in $z \gtrsim 2$ galaxies. These UV data are combined with archival HST imaging at longer wavelengths and spectroscopic redshifts from the literature to search for potential LyC candidates over the full redshift range $z = 1.9$ to $z = 4$, using a technique similar to the one presented in V15. This search also provides the basis for a future paper in which we will constrain the average escape fraction of SFGs at $z \sim 2-3$.

This paper is structured as follows. In Section 2 we describe the imaging and the spectroscopic data used for the analysis. The methodology of our candidate search is outlined in Section 3. We present six candidate LyC emitters in Section 4, calculate their $f_{esc}$ and situate them in the context of other efforts to understand LyC leakage in Section 5, and finally summarize our findings while looking toward the future in Section 6.

Throughout this paper, we adopt $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, $H_0 = 70$ km s⁻¹ Mpc⁻¹, i.e., $h = 0.7$, largely consistent with the most recent measurements from Planck (Planck Collaboration et al. 2016). Magnitudes are given in the AB system (Oke & Gunn 1983).

## 2. Data

### 2.1. Photometry

As shown previously, any study of LyC emission at high redshift requires data at excellent spatial resolution in order to avoid contaminating flux from foreground sources that lie close

in projection along the line of sight (e.g., Vanzella et al. 2010). Hence, in this paper we only analyze objects for which HST images are available. The novel data that let us search for LyC candidates in the relatively unexplored redshift range of $z \sim 2-3$ come from deep UV imaging (down to 27.5–28.0 mag at $5\sigma$) of the GOODS-North and GOODS-South fields in the F275W and F336W bands acquired by the HDUV legacy survey (GO-13871; see Oesch et al., 2016, submitted), including all the F275W data taken by the CANDELS survey (Koekemoer et al. 2011). Additionally, we include the previously released version 2 of the UVUDF images<fnref n="15" /> (Teplitz et al. 2013; Rafelski et al. 2015).

At $z \sim 2-3$, the LyC is probed by the HDUV bands, and we show in Section 3 how LyC leakage may be inferred from these photometric data. Crucially, the HDUV survey’s coverage area is a subset of that of the 3D-HST (Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016) and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) surveys, as well as the previous GOODS ACS imaging (Giavalisco et al. 2004). This complementarity provides continuous multiwavelength imaging (using ACS and WFC3) from the UV to the near-IR, along with an abundance of grism redshifts (from 3D-HST), and facilitates the reliable calculation of UV continuum slopes ($\beta$) and the derivation of physical properties of galaxies through spectral energy distribution (SED) fitting.

### 2.2. 3D-HST Grism Spectra and Other Spec-Z

For a reliable LyC emitter search, we require accurate spectroscopic redshifts. Fortunately, the HDUV/GOODS fields have extended spectroscopic coverage from several surveys conducted over many years (e.g., Dawson et al. 2001; Cowie et al. 2004; Reddy et al. 2006; Wuyts et al. 2008; Yoshikawa et al. 2010). Most of the secure redshifts from these surveys are already compiled in the 3D-HST catalogs from Skelton et al. (2014), which we use for our analysis. We also harvest newly available spectroscopic redshifts from the VUDS (DR1; Le Fèvre et al. 2015) and MOSDEF (Kriek et al. 2015) surveys. We ensure that only high-probability redshifts are included in our analysis (e.g., confidence class 3+ in the VUDS release).

Additionally, grism spectra from the 3D-HST survey are available for most of the sources that we analyze, and they are particularly reliable when prominent emission lines are detected (Momcheva et al. 2016). The WFC3/G141 grism used in the 3D-HST survey spans 1.1–1.7 µm and is perfectly situated to capture the distinctive [O III] $\lambda\lambda4959, 5007$ doublet at $z \sim 1.9–2.4$. Thus, the grism redshifts derived in this redshift window are anchored to well-detected emission lines, and this gives us a sizable sample of sources for which the LyC is reliably located in the HDUV/F275W filter. Furthermore, the near-IR spectra enable the analysis of emission-line ratios. In a narrow window around $z \sim 1.9–2.0$, both [O II] and [O III] fall in the G141 grism’s spectral range, and almost always $H\beta$ is available along with [O III], though blended with [O III] given the very low spectral resolution of the grism.

## 3. Methodology

The HDUV filter set probes exclusively ionizing photons for galaxies at $z > 2.4$ in F275W and at $z > 3.1$ in F336W. For

<footnote n="15">¹⁵ https://archive.stsci.edu/prepds/uvudf/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Summary of the selection technique adopted in this paper illustrated with the first LyC candidate GS 30668 ($z = 2.172$). Top left: at GS 30668's redshift, >50% of the F275W band is covered by the LyC (shown as the shaded blue region), and hence LyC leakage may be inferred from the flux observed in that band. The SED shown is a BPASSv2 model SED with $\beta_{UV}$ and age consistent with GS 30668, and $E(B - V) = 0.03$ (using the dust extinction curve from Reddy et al. 2016a). Top right: four realizations of the IGM transmission curves toward a $z = 2$ galaxy from the Inoue et al. (2014) Monte Carlo sampling, with $\lambda = 912$ Å indicated by blue dotted lines. A total of 10,000 such curves are convolved with a Gaussian distribution of SEDs ($\mu = \beta_{Obs}$, $\sigma = \sigma_{\beta_{Obs}}$) in our method to compute the expected color distribution in the bottom left panel. Bottom left: the observed F275W–F336W color (red dot with error bars) lies almost entirely blueward of the color distribution generated under the assumption that $f_{esc} = 0$ (shown in green). The probability for $f_{esc} = 0$ for GS 30668 is $P(f_{esc} = 0) = 2\%$, which makes it an LyC candidate. Further, the actually measured color perfectly lines up with the distribution generated under the assumption that $f_{esc} = 1$ (shown in blue), strongly suggesting a high value of $f_{esc}$ for this galaxy. Bottom right: same color distributions as on the left, but following the original V15 method, in which a single highly ionizing SED is used (selected from our BPASSv2 SED grid that uses the Reddy et al. 2016a dust curve; see the text), and there is no Monte Carlo treatment of the observed color. This method results in a tighter color distribution, but it also calculates a very low $P(f_{esc} = 0)$ for GS 30668. All six candidates found using our method of Gaussian SED distributions are also candidates according to the V15 method.</figcaption></figure>

sources at somewhat lower redshift, the filters cover both ionizing and nonionizing wavelengths (see Figure 1 for an example of a $z = 2$ galaxy). However, even at these lower redshifts, LyC emitters with non-negligible $f_{esc}$ can be identified by modeling the UV SED and estimating the contribution of nonionizing photons to the filter flux of a given galaxy. In particular, V15 developed such a technique, which was successful in selecting Ion2 as a highly probable LyC candidate. Ion2 was subsequently followed up and confirmed with HST imaging and until recently was the only spectroscopically confirmed LyC leaker at high redshift ($z = 3.212$; see de Barros et al. 2016; Vanzella et al. 2016, for more details on Ion2).

In this paper, we build on the V15 technique and adapt it to use a new dust curve and different SED models before applying it to the HDUV data set. In brief, we fit the UV continuum

slope of a galaxy to identify representative UV model SEDs, which are then used to derive the expected color of that galaxy in the HST filters straddling the LyC edge. The latter is done via a Monte Carlo simulation of the IGM transmission representing 10,000 lines of sight. We then identify galaxies whose measured HST colors are inconsistent with an escape fraction of $f_{esc} = 0$ but indicate $f_{esc} > 0$.

Explicit details of the selection procedure of sources with non-negligible $f_{esc}$ are described below:

1. *Input Galaxy Sample:* In order to obtain reliable results, we only apply our method to galaxies that have a secure spectroscopic redshift from the literature, as well as a reliable flux measurement ($S/N > 3$) available in a filter that probes >50% of LyC photons (i.e., such that the central wavelength of the filter lies below rest-frame 912 Å).
The idea here is that if the galaxy is an LyC leaker, the measured filter flux will contain a contribution from LyC photons that we can attempt to infer. A secure redshift is thus important, since it ensures that the LyC falls within a particular filter.

2. SED grid: We assemble our SED grid using the latest Binary Population and Spectral Synthesis models (BPASSv2; J. J. Eldridge et al. 2017, in preparation) rather than the Bruzual & Charlot (2003, henceforth BC03) models used in V15. This choice is motivated by the better performance of BPASS models in matching the spectral properties of SFGs at $z \sim 2-3$ (e.g., Steidel et al. 2016; Strom et al. 2017), as well as those of young, massive star clusters, which dominate the rest-frame UV region (Wofford et al. 2016). The BPASS models are also more consistent with the intrinsic 900–1500 Å flux density ratio inferred for $\sim L^*$ galaxies at $z \sim 3$ (Reddy et al. 2016b).

We use the fiducial BPASSv2 galaxy templates (initial mass function [IMF] with a slope of $-1.30$ between 0.1 and $0.5M_{\odot}$ and a slope of $-2.35$ between 0.5 and $100M_{\odot}$), to which we self-consistently added nebular continuum and line emission to build an extremely fine template grid that uniformly spans various ages, metallicities, and magnitudes of dust extinction. The following parameter space is covered in our grid: $Z/Z_{\odot} = 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.7, 1, 1.5, 2; E(B-V)$ from 0 to 0.6 in linear steps of 0.03; age from 1 Myr to 10 Gyr in steps of 0.1 dex; and a constant star formation history. To account for dust attenuation, we use the Calzetti attenuation curve (Calzetti et al. 2000) for wavelengths greater than 1500 Å. In the UV region blueward of 1500 Å that is critical to our study, instead of simply extrapolating the Calzetti curve (as was done in previous studies, like V15), we use the newly derived dust curve from Reddy et al. (2016a). This new curve predicts a factor of $\sim 2$ lower dust attenuation of LyC photons than the Calzetti curve for $E(B-V) \sim 0.15$, typical for $\sim L^*$ galaxies. The difference in dust attenuation is smaller for bluer $E(B-V)$.

3. $\beta_{UV}$—IGM Monte Carlo color simulation and candidate selection: In order to select SEDs that best represent a given galaxy, we compare the observed UV continuum slopes ($\beta_{obs}$) with the ones from the SEDs ($\beta_{SED}$). $\beta_{obs}$ is calculated as described in Castellano et al. (2012) by performing multiband fitting assuming $f_{\lambda} \propto \lambda^{\beta}$. The flux measured in the HST filters that span rest-frame 1300–3300 Å depending on the object’s redshift is used to measure its $\beta_{obs}$. Theil–Sen regression, which is robust against outliers (Wilcox 1998), is used to fit the slope and calculate confidence intervals ($\sigma_{\beta_{obs}}$).

In the original V15 method, every galaxy is only matched to a single SED: from the grid described above, one selects the SEDs that satisfy $\beta_{Obs} - \sigma_{\beta_{obs}} < \beta_{SED} < \beta_{Obs} + \sigma_{\beta_{obs}}$ and $\text{Age}_{SED} < \text{Age}_{\text{Universe}}$ at $z_{spec}$, and from these, one picks the SED with the maximum $LLyC/L1500$ (LLyC is calculated over 850–900 Å). This results in the bluest simulated color when $f_{\text{esc}} = 0$ is assumed (calculated for the LyC-containing filter and an adjacent filter). For a galaxy to qualify as an LyC leaker, the galaxy must display a color bluer than the bluest simulated color. This method is thus the most conservative to identify galaxies that are inconsistent with a zero escape fraction.

In the method used in this work, instead of relying on a single extremely ionizing template, we account for the variance in potential ionizing fluxes of a given galaxy by performing weighted sampling (10,000 times) from the SED library. Each SED in the library is given a weight based on its UV continuum slope as per a Gaussian distribution centered at $\beta_{obs}$ and with $\sigma = \sigma_{\beta_{obs}}$. For each of the 10,000 selected SED templates, we apply a realization of IGM attenuation drawn from 10,000 possible sightlines at the $z_{spec}$ of the given galaxy (IGM transmission functions from Inoue et al. 2014). The IGM-applied SEDs are used for a Monte Carlo sampling of the expected color distribution of the galaxy in our HST filters by assuming $f_{\text{esc}} = 0$. In the V15 method, the color distribution arises from the convolution of the single extremely ionizing template with the 10,000 IGM sightlines and filters.

If the color distribution, which has no LyC photons contributing to it (i.e., assuming $f_{\text{esc}} = 0$), falls largely redward of the observed color, the galaxy is an LyC candidate since even in transparent IGM sightlines LyC photons would be required to reproduce the observed color.

This idea is formalized in terms of a probability (Equation (1) from V15):

$$P(f_{\text{esc}} = 0) = N_{\text{color}}/N_{\text{total}}. \quad (1)$$

1. $N_{\text{color}}$ is the number of IGM-$\beta$-Color$_{obs}$ realizations for which Color$_{simulated} > Color_{obs}$.

2. $N_{\text{total}} = 10^9$. We treat Color$_{obs}$ as a Gaussian distribution with width corresponding to the photometric scatter and sample 100,000 times from this distribution. Each of the 100,000 Color$_{obs}$ samples is compared to the expected color from each of the 10,000 IGM-applied SEDs, which have been described earlier. Thus, $N_{\text{total}} = 100,000$ Color$_{obs}$ samples $\times 10,000$ IGM-applied SEDs $= 10^9$.

LyC candidates are the galaxies for which this probability $P(f_{\text{esc}} = 0)$ is low. We repeat the calculation of color histograms under the assumption of $f_{\text{esc}} = 1$ for a consistency check (see blue histograms in Figures 1 and 2).

While both the V15 method and our method can be used to select LyC leakers, our method can also be used to constrain $f_{\text{esc}}$ since it accounts for the entire diversity in the UV SED of the galaxy, as well as the variation in the IGM transmission (see Section 5.1 for more details). In general, the two methods identify the same candidate LyC emitters. However, there are some marked differences. For instance, in Figure 1, we show the expected color distributions computed with both methods for the first candidate selected in the HDUV data set (GS 30668). This source is a high-probability candidate with $P(f_{\text{esc}} = 0) = 2\%$ (our method) and $P(f_{\text{esc}} = 0) = 0.02\%$ (V15). However, the color distributions computed using our method are much wider, since we adopt not just a single SED but a distribution of SEDs according to the variance in the $\beta_{obs}$ of the galaxy. In the case of V15, the same SED is used for all realizations, and the color distribution is only a reflection of the variance of IGM transmission functions.

## 4. Candidate LyC Leakers

The procedure outlined in Section 3 allows us to identify likely LyC emitter sources in the redshift range $z = 1.9–4$. The upper limit in the redshift range is primarily imposed by the decreasing IGM transmission to higher redshift, which makes it
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Selection of LyC candidates from Monte Carlo color distributions (see also Figure 1, bottom left panel). The galaxies whose simulated color distributions are shown here (all at $z \sim 2$) are selected as LyC candidates since their observed F275W–F336W colors (indicated in red with error bar) are inconsistent with the $f_{esc} = 0$ distributions (shown in green), resulting in low $P(f_{esc} = 0)$ values.</figcaption></figure>

virtually impossible to directly probe LyC photons at $z > 4$ based on broadband images. Below the lower redshift limit, our bluest filter (F275W) contains <50% of the LyC. We thus applied our procedure to all sources listed in the 3D-HST GOODS catalogs with a secure redshift in the range $z = 1.9–4$ for which the relevant photometry was available.

For a source to qualify as an LyC emitter candidate, we set the following criteria: (1) $3\sigma$ detection in the LyC-containing band, $5\sigma$ detection in the adjacent redward band used to calculate the color distribution, and S/N greater than 5 for the calculated color; (2) flux measurements from at least three bands to fit the $\beta_{UV}$ slope; (3) clean morphology in all available HST images to rule out contamination from low-$z$ interlopers with chance projections; and (4) $P(f_{esc} = 0) < 15\%$ (set

arbitrarily). By “clean morphology” we mean that the source has no near neighbors coincident with the hypothesized LyC flux and that the source has uniform color in its *VJH* red, green, and blue (RGB) stamps (F606W, F125W, and F160W cutouts form the three channels of the RGB stamps).

Using these criteria, we identify six candidates in the HDUV +UVUDF survey area. Their color histograms are shown in Figure 1 for the first source and Figure 2 for the remaining five, while their basic properties are listed in Table 1. Multiwavelength HST stamps of the six candidates are shown in Figure 3. In principle, we extended our search up to $z \sim 4$; however, all six candidates lie at $z \sim 2$, and they were selected based on the F275W–F336W color. This is not necessarily surprising, given that most of our input spectroscopic redshifts from the literature,
<table><caption>**Table 1**
Summary of Lyc Candidates</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">Star-forming Galaxies</th><th colspan="3">Active Galactic Nuclei<sup>a</sup></th></tr><tr><th>GS 30668</th><th>GS 35257</th><th>GS 14633</th><th>GN 21231</th><th>GN 19591</th><th>GN 19913</th></tr></thead><tbody><tr><td>R.A.</td><td>3:32:35.47</td><td>3:32:24.93</td><td>3:32:46.46</td><td>12:36:46.74</td><td>12:36:48.31</td><td>12:36:35.6</td></tr><tr><td>Decl.</td><td>-27:46:16.89</td><td>-27:44:51.61</td><td>-27:50:36.64</td><td>+62:14:45.97</td><td>+62:14:16.64</td><td>+62:14:24.0</td></tr><tr><td>Redshift<sup>b</sup></td><td>2.172<sup>+0.001</sup><sub>-0.003</sub></td><td>2.107<sup>+0.002</sup><sub>-0.003</sub></td><td>2.003<sup>+0.001</sup><sub>-0.002</sub></td><td>2.004<sup>+0.002</sup><sub>-0.002</sub></td><td>1.998<sup>+0.001</sup><sub>-0.003</sub></td><td>2.012<sup>+0.001</sup><sub>-0.002</sub></td></tr><tr><td>f<sub>esc</sub> [%]<sup>c</sup></td><td>60<sup>+40</sup><sub>-38</sub></td><td>72<sup>+28</sup><sub>-48</sub></td><td>62<sup>+38</sup><sub>-31</sub></td><td>&gt;71</td><td>&gt;13</td><td>(&sim;100)<sup>d</sup></td></tr><tr><td>UV slope &beta; (f<sub>&lambda;</sub> &propto; &lambda;<sup>&beta;</sup>)</td><td>-2.23 &plusmn; 0.08</td><td>-1.93 &plusmn; 0.07</td><td>-1.92 &plusmn; 0.01</td><td>-1.86 &plusmn; 0.07</td><td>-1.27 &plusmn; 0.07</td><td>-1.11 &plusmn; 0.41</td></tr><tr><td>log M<sub>gal</sub>/M<sub>&odot;</sub><sup>e</sup></td><td>9.07<sup>+0.05</sup><sub>-0.06</sub></td><td>9.37<sup>+0.00</sup><sub>-0.03</sub></td><td>9.23<sup>+0.00</sup><sub>-0.01</sub></td><td>10.44<sup>+0.00</sup><sub>-0.03</sub></td><td>9.99<sup>+0.12</sup><sub>-0.03</sub></td><td>11.3<sup>+0.00</sup><sub>-0.06</sub></td></tr><tr><td>log SFR [M<sub>&odot;</sub> yr<sup>-1</sup>]<sup>e</sup></td><td>0.24<sup>+0.02</sup><sub>-0.04</sub></td><td>0.38<sup>+0.0</sup><sub>-0.33</sub></td><td>-0.56<sup>+0.00</sup><sub>-1.39</sub></td><td>1.06<sup>+0.01</sup><sub>-0.00</sub></td><td>0.92<sup>+1.29</sup><sub>-0.00</sub></td><td>1.55<sup>+0.00</sup><sub>-0.41</sub></td></tr><tr><td>log sSFR [yr<sup>-1</sup>]<sup>e</sup></td><td>-8.83<sup>+0.07</sup><sub>-0.04</sub></td><td>-8.99<sup>+0.00</sup><sub>-1.31</sub></td><td>-9.79<sup>+0.00</sup><sub>-1.38</sub></td><td>-9.38<sup>+0.04</sup><sub>-0.00</sub></td><td>-9.07<sup>+1.30</sup><sub>-0.00</sub></td><td>-9.75<sup>+0.00</sup><sub>-0.35</sub></td></tr><tr><td>log Age[yr]<sup>e</sup></td><td>8.9<sup>+0.1</sup><sub>-0.2</sub></td><td>8.0<sup>+0.1</sup><sub>-0.1</sub></td><td>8.0<sup>+0.0</sup><sub>-0.0</sub></td><td>9.5<sup>+0.0</sup><sub>-0.0</sub></td><td>7.7<sup>+0.6</sup><sub>-0.1</sub></td><td>9.5<sup>+0.0</sup><sub>-0.0</sub></td></tr><tr><td>E(B &minus; V)<sup>e</sup></td><td>0.0<sup>+0.0</sup><sub>-0.0</sub></td><td>0.07<sup>+0.01</sup><sub>-0.03</sub></td><td>0.02<sup>+0.00</sup><sub>-0.02</sub></td><td>0.0<sup>+0.01</sup><sub>-0.00</sub></td><td>0.22<sup>+0.10</sup><sub>-0.00</sub></td><td>0.17<sup>+0.01</sup><sub>-0.10</sub></td></tr><tr><td>[O III]/[O II]<sup>f</sup></td><td>9.47 &plusmn; 3.81</td><td>&gt;5.5</td><td>3.23 &plusmn; 1.41</td><td>0.81 &plusmn; 0.12</td><td>2.23 &plusmn; 0.17</td><td>10.33 &plusmn; 5.43</td></tr><tr><td>EW<sub>rest</sub>([O III]) [Å]<sup>f</sup></td><td>1211 &plusmn; 55<sup>g</sup></td><td>168 &plusmn; 31</td><td>501 &plusmn; 89</td><td>102 &plusmn; 10</td><td>384 &plusmn; 22</td><td>94 &plusmn; 4</td></tr><tr><td>EW<sub>rest</sub>(H&beta;)[Å]</td><td>130 &plusmn; 45<sup>g</sup></td><td>&lt;19</td><td>95 &plusmn; 38</td><td>49 &plusmn; 9</td><td>97 &plusmn; 10</td><td>13 &plusmn; 3</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Note that these fits do not use AGN templates.

<sup>b</sup> Redshift is the 3D-HST z<sub>grism</sub> for all sources.

<sup>c</sup> Median of f<sub>esc</sub> distribution with 16th and 84th percentile error bars. When >50% of the distribution is truncated, we state the 10th percentile as a lower limit (see Section 5.1).

<sup>d</sup> For GN 19913, the only candidate emitter with LyC emission consistent with an active nucleus, we use two sets of AGN templates to estimate f<sub>esc</sub> as 23<sup>+21</sup><sub>-15</sub>% and 14<sup>+21</sup><sub>-14</sub>% (see Appendix and Section 5.1).

<sup>e</sup> Derived using FAST (Kriek et al. 2009) and BC03.

<sup>f</sup> [O III] refers to the merged doublet, i.e., [O III] λ4959 + [O III] λ5007.

<sup>g</sup> The 3D-HST pipeline overestimates equivalent widths when the continuum detected by the grism is very faint, like in the case of GS 30668 (see Figure 4). In such a case we use the grism line flux along with a continuum extrapolation from EAZY (Brammer et al. 2008) to calculate the equivalent width.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** 3" HST postage stamp images of LyC candidates. The F275W and F336W images (first two columns) are the new data acquired by the HDUV survey, while the rest are from HST GOODS archival imaging (see Section 2.1). These high-resolution, multiwavelength images allow us to conclusively rule out flux contamination from neighboring sources. For instance, in the case of GN 19591 we investigate a potential interloping clump toward the bottom right of the central source (visible in F435W, F606W, F775W), but we conclude that this does not pose a problem (see Section 4.2)—such a check would not be possible at ground-based resolution. In terms of morphology, our candidates are in general compact. Two curiosities are the AGNs GN 21231 and GN 19591, whose F275W detections (which include hypothesized LyC leakage) appear to be extended and not concentrated on a central point source (see the AGN GN 19913 for comparison, which is barely resolved in F275W), which may implicate stars instead of the active nucleus as the origin of LyC flux (see Section 4.2 for details).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** 3D-HST grism spectra of LyC candidates. In each panel the top strip contains the 2D spectrum of the source and the bottom part shows the extracted 1D spectrum (black), along with a best-fit model (red) (see Momcheva et al. 2016, for details about the 3D-HST pipeline). In all the spectra, the distinctive merged [O III] doublet (O IIIλ4959 + O III λ5007) is unambiguously detected. On the other hand, the O II and Hβ detections are in general tentative. Deeper, high-resolution spectra are required to better constrain [O III]/[O II] and Hβ/β<sub>UV</sub>, two promising indirect methods to infer the escape fraction of galaxies.</figcaption></figure>

and in particular from the 3D-HST grism data, lie at $z \sim 2$. In total, our input galaxy sample with reliable spectroscopic redshifts contained 1124 galaxies.

We pay careful attention to redshift quality, since our entire selection procedure hinges on secure redshifts. All our candidates have 3D-HST grism spectra with well-detected emission lines (shown in Figure 4), and we use the associated $z_{\text{grism}}$ measurements in our analysis. Additionally, spectroscopic redshifts were already available from the literature for four of the sources that corroborate the $z_{\text{grism}}$ (Reddy et al. 2006, for GN 19591, GN 19913, and GN 21231; Kirkpatrick et al. 2012, for GN 21231; Trump et al. 2011, for GS 30668).

We investigated the sources in detail and tested whether they show any sign of AGN activity that could contribute to the ionizing photons, or whether they had any nearby neighbors that could contaminate the UV color measurements. We split

the sample into two classes: (1) SFG candidates and (2) likely AGNs. These classes are discussed separately in the following sections.

### 4.1. Star-forming Galaxy Candidates

The first three of our candidates are the most convincing, as they show no signs of AGN activity and no nearby, potentially contaminating galaxies. These are GS 30668, GS 35257, and GS 14633. These sources are not detected in the deepest available Chandra GOODS images and associated CANDELS-matched catalogs that include even extremely faint X-ray sources down to a flux of $1(8) \times 10^{17}$ erg cm$^{-2}$ s$^{-1}$ in the [0.5–2] keV ([0.5–10] keV) energy band (Cappelluti et al. 2016). For two of these sources (GS 30668 and GS 14633) we can further examine the Mass–Excitation diagram (Juneau et al. 2014), since [O III] ($>5\sigma$) and Hβ ($>2\sigma$) measurements are available.
GS 30668 and GS 14633 are within the star-forming region of the diagram (Figure 7, Appendix). The compact and isolated objects in this section are likely galaxies where ionizing radiation escapes, similar to the source Ion2 identified in V15 using a technique similar to the one we have adopted here.

## 4.2. AGN

Two sources in our sample are known AGNs: GN 21231 and GN 19913, which have been well studied in the literature (GN 19913: e.g., Smail et al. 2004; Bluck et al. 2011; GN 21231: e.g., Evans et al. 2010; Kirkpatrick et al. 2012). In addition to them, we classify GN 19591 as a likely AGN based on a strong Spitzer/MIPS and even a Herschel detection, as well as the hint of a power-law SED across the Spitzer/IRAC bands. However, we note that GN 19591 is not detected in the Chandra Deep Field North (2 Ms) images and has spectral features that could easily belong to an SFG (e.g., [O III]/[O II] ~ 2.2, [O III] λ5007/Hβ ~ 3.2).

GN 21231 was identified as an AGN by Kirkpatrick et al. (2012) and was found to display the distinctive 9.7 µm silicate line in its IR spectrum, which indicates a large column of dust along the line of sight. This massive amount of dust is perhaps why GN 21231 is heavily obscured and not detected in Chandra images. The dust probed by the 9.7 µm line is concentrated in the innermost, central 2 pc of the galaxy (Köhler & Li 2010), however. It is thus fair to conjecture that the ionizing flux from the active nucleus in GN 21231 is largely suppressed by the dust torus around it. Thus, any LyC flux we detect probably originates from the star-forming parts of the galaxy. This hypothesis is supported by our FAST SED fit, which yields $E(B-V) = 0$, indicating that while the central nucleus may be dusty, the rest of the galaxy is not. The morphology of the ionizing radiation in the F275W filter, which is extended and not concentrated on a central point source, is another indication that the LyC flux originates from stars.

Similar to GN 21231, GN 19591, which is not detected in Chandra 2 Ms images, might host a heavily obscured AGN. Based on the extended, resolved LyC radiation in the F275W image, it is likely that star-forming regions are responsible for the ionizing flux, however. The similarity of morphologies in the F336W and F275W images is further evidence of the star-forming regions being correlated with LyC emission.

In order to quantitatively verify that the hypothesized LyC emission observed in the F275W images of GN 21231 and GN 19591 may have its origins in the star-forming parts of these galaxies, we measure the FWHM of the flux and construct radial profiles (see Figure 8, Appendix). Indeed, while the FWHM of the point-spread function (PSF) in the F275W image is 0.11, the FWHM of GN 21231 and GN 19591 is 0.17 ± 0.01 and 0.26 ± 0.08, respectively. The F275W flux being >50% wider than the PSF is a strong indication that the ionizing radiation is not dominated by a point source, i.e., the nucleus of the AGN. Viewed together, these two obscured AGN candidates hint that AGNs may play an important role in clearing the interstellar medium (ISM) and facilitating the escape of LyC photons emitted by stars in these galaxies.

In contrast, the situation in the AGN GN 19913 is clearer, and indicates that it is not useful for a study of LyC escape from stars. The UV morphology of GN 19913 is concentrated on a central point source, consistent with the ionizing photons originating from the AGN.

# 5. Discussion

## 5.1. Probability Distributions of the LyC Escape Fraction

We can use the color histograms presented in Figures 1 and 2 to estimate probability distributions of the absolute $f_{esc}$ for each source. In particular, for each realization of the IGM transmission we can compute the required escape fraction to bring the simulated F275W-F336W color into agreement with the observed color. For a single IGM line of sight, $f_{esc}$ is calculated as follows:

$$f_{esc} = \frac{10^{-0.4 \text{ Color}_{\text{Obs}}} - 10^{-0.4 \text{ Color}_{f_{esc}=0}}}{10^{-0.4 \text{ Color}_{f_{esc}=1}} - 10^{-0.4 \text{ Color}_{f_{esc}=0}}}. \quad (2)$$

1. Color$_{Obs}$ is the F275W-F336W color of the source. This quantity is sampled $10^5$ times from a Gaussian distribution that accounts for the photometric scatter (as described in Section 3).

2. Color$_{f_{esc}=0}$ and Color$_{f_{esc}=1}$ are the simulated colors under the assumption of $f_{esc} = 0$ and $f_{esc} = 1$, respectively.

We limit $0 \le f_{esc} \le 1$, and so for a particular sightline if the observed color is bluer (redder) than the simulated $f_{esc} = 1$ ($f_{esc} = 0$) color, we truncate $f_{esc}$ to 1 (0). For every candidate we have a billion estimates of $f_{esc}$ ($10^5$ samplings from the observed color Gaussian for each of the 10,000 IGM realizations), based on which we compute the cumulative distribution function of $f_{esc}$ (see Figure 5). In Table 1 we summarize these measurements. We state the median $f_{esc}$ with the 16th and 84th percentiles as error bars. When >50% of the $f_{esc}$ estimates for a source are truncated at $f_{esc} = 1$, we state the 10th percentile as a lower bound.

The estimated $f_{esc}$ values for our candidates range from ~60% to ~100%.<fnref n="16" /> Such high escape fractions are expected, given our selection procedure in which we set the stringent criterion of $P(f_{esc} = 0) < 15\%$. This is also in agreement with the findings of Ion2, which was selected as an LyC candidate by the V15 method and was later confirmed to have an $f_{esc} \gtrsim 50\%$ through follow-up imaging with HST (Vanzella et al. 2016). Note that it is likely that many sources in our parent input sample show significant escape fractions, but they are missed in our selection since we cannot reliably separate the two color histograms. A future paper, currently in preparation, will address the average escape fraction of galaxies in the HDUV fields. In general, our findings are in broad agreement with previous studies that found only a small fraction of galaxies to show high $f_{esc}$. This could be explained by generally high gas covering fractions with few clear sightlines out of galaxies or similarly by the fact that the escape of ionizing photons is a stochastic process, with short periods of time of high $f_{esc}$ (e.g., Wise et al. 2014; Cen & Kimm 2015).

We have not included AGN templates in our estimates of $f_{esc}$. For GN 21231 and GN 19591 this is because, as discussed in Section 4.2, the LyC flux seems to be originating from the star-forming parts of the galaxy. And since we calculate $f_{esc}$ using LyC-UV colors, it would be inappropriate to use AGN SEDs. This line of reasoning is shored up by the excellent

<footnote n="16">¹⁶ Note that $f_{esc} = 100\%$ or higher either is unphysical or can be excluded owing to our detection of strong emission lines in these objects (Figure 4). Such high inferred values rather reflect our limited knowledge of the intrinsic UV SEDs below the Lyman limit.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example $f_{\text{esc}}$ distributions for two of our candidates. The cumulative $f_{\text{esc}}$ distributions based on 100,000 $f_{\text{esc}}$ realizations (10 observed color samplings per IGM line of sight) are shown, with percentiles for the distribution indicated in blue text. In the left panel, for SFG GS 30668, we state the median $f_{\text{esc}}$. On the right, more than half of the $f_{\text{esc}}$ values calculated for the source GN 21231 are truncated to 1. In this case, we state a lower limit equal to the 10th percentile. See Table 1 for similarly calculated $f_{\text{esc}}$ values for all our candidates.</figcaption></figure>

agreement of the observed color for these two sources with the simulated color from the BPASS templates (see second row of Figure 2), as well as the divergence of their observed color from that predicted by pure-AGN SED templates (Stevans et al. 2014; Siebenmorgen et al. 2015). For GN 19913, the LyC flux appears to be due to AGN activity, and our estimate of $f_{\text{esc}}$ from the BPASS SEDs is an unphysical ~100%. We calculate and quote this number, however, since it may be instructive to future users of our technique who wish to find pure-AGN LyC candidates or simply select AGNs based on UV color excesses from large data sets that have not been pre-classified into AGNs and SFGs (as discussed and envisioned in Vanzella et al. 2015). We have also calculated $f_{\text{esc}}$ for GN 19913 using the AGN SEDs described in Stevans et al. (2014) and Siebenmorgen et al. (2015) as $23_{-15}^{+21}\%$ and $14_{-14}^{+21}\%$, respectively, and confirmed its LyC candidature. The simulated color distributions for GN 19913 (shown in Figure 9, Appendix) generated with AGN SEDs agree well with its observed color, which was too blue for almost all the sightlines simulated with BPASS templates (bottom panel of Figure 2).

## 5.2. Comparison with Confirmed High-z LyC Leakers

Here we provide a short comparison of our LyC emitter candidates with the three previous, confirmed LyC sources at $z > 2$: Ion2 ($z = 3.21$), Q1549-C25 ($z = 3.15$), and MD5b ($z = 3.14$). In particular, we estimated several physical parameters for our LyC leaker candidates based on the broadband photometry and SED fitting using FAST (Kriek et al. 2009). We use the same models as the 3D-HST survey (Skelton et al. 2014), except with a metallicity of $0.2Z_{\odot}$. Specifically, we use BC03 models with a Chabrier (2003) IMF, reddened by the Calzetti et al. (2000) dust curve, and exponentially declining SFHs, with $\tau$ in the range $\log\tau \text{ yr}^{-1} = 7-10$. We allow for ages in the range $\log A \text{ yr}^{-1} = 7.6$ up to the age of the universe at the given redshift. The filters are the full filter set of the Skelton et al. (2014) catalog, plus the two HDUV filters. No AGN templates are used in SED fitting.

The derived parameters are summarized in Table 1. We confirmed that the strong emission lines do not affect the physical parameters significantly by excluding the filters that were most affected for each galaxy (H band and K band). The parameters calculated by excluding these bands are essentially

identical to those in Table 1, except with slightly larger uncertainties. In general, our galaxies have $\sim L_{\text{UV}}$ at $z \sim 2$ (Reddy & Steidel 2009).

Ion2 and GS 30668 share many remarkable similarities. Both display an extreme EW([O III]$_{4959+5007}$ + H$\beta$) (Ion2: ~1600 Å; GS 30668: ~1350 Å) and [O III]/[O II] $\gtrsim 10$. Even the multiband fitted $\beta_{\text{UV}}$ (Ion2: $-2.2 \pm 0.2$; GS 30668: $-2.2 \pm 0.1$), EW(H$\beta$) (Ion2: $112 \pm 60$ Å; GS 30668: $130 \pm 45$ Å) and $E(B-V) = 0$ for these sources resemble each other.

Broadly speaking, our candidates display little to no dust extinction, consistent with Ion2, Q1549-C25, and MD5b. The exception is GN 19591 with $E(B-V) = 0.22_{-0.00}^{+0.10}$. This trend supports the idea that dust attenuation is not conducive to the escape of LyC radiation. Furthermore, similar to previous LyC leakers (in particular MD5b), the SFGs in our sample are generally very young (~50–160 Myr). GS 30668 is the exception, with an age of ~800 Myr (once again, consistent with Ion2’s reported age).

Our sources thus provide some evidence that LyC emission generally occurs in young galaxies with little dust extinction and thus blue SEDs, or in sources whose ISM is highly excited and that show very strong rest-frame optical emission lines (such as Ion2) (consistent with Jones et al. 2013; Wofford et al. 2013; Borthakur et al. 2014; Alexandroff et al. 2015; Rivera-Thorsen et al. 2015; Trainor et al. 2015; Dijkstra et al. 2016; Nakajima et al. 2016; Reddy et al. 2016b). However, a more systematic study of the average escape fractions of galaxies in the HDUV field will have to be performed to better connect the physical properties that lead to LyC emission with significant $f_{\text{esc}}$ from galaxies.

## 5.3. The Peak Escape Fraction as a Function of Redshift

It is interesting to put the SFGs identified in this paper in a broader context. Figure 6 shows a compilation of absolute $f_{\text{esc}}$ measurements for star-forming sources reported in the literature. These include recent direct detections for individual galaxies at low redshift (Leitet et al. 2011; Borthakur et al. 2014; Izotov et al. 2016a, 2016b; Leitherer et al. 2016), as well as detections or limits from individual $z > 2$ sources (Mostardi et al. 2015; Shapley et al. 2016; Vanzella et al. 2016; Vasei et al. 2016) or averages from small subsamples of galaxies at $z > 2$
<figure><img src="image_8.png" /></figure>

**Figure 6.** Compilation of absolute $f_{esc}$ measurements for star-forming sources reported in the literature. The red symbols include direct detections from galaxies at low redshift (Leitet et al. 2011; Borthakur et al. 2014; Izotov et al. 2016a, 2016b; Leitherer et al. 2016), as well as detections or limits from $z > 2$ sources using different methods (Mostardi et al. 2015; Leethochawalit et al. 2016; Matthee et al. 2017; Shapley et al. 2016; Vanzella et al. 2016; Vasei et al. 2016). The plot does not show population-average escape fractions and depicts $f_{esc}$ measured from individual galaxies and small subsamples. The candidate sources studied in this paper (shown in purple) occupy the relatively unexplored $z \sim 2$ region in redshift space, and they double the number of direct high-$z$ $f_{esc}$ detections. The redshift of some sources was slightly offset for clarity. The shaded area in the upper half of the graph represents $f_{esc} > 10\%$, a necessary condition for SFGs to drive reionization. While the population-average escape fraction at $z > 2$ still has to be measured reliably, it is clear that at least some of the few individually detected sources at high redshift satisfy this criterion. Only one such source is currently known at $z < 0.5$, hinting at a possible evolution of the maximally achievable escape fraction as a function of cosmic time. Note, however, that the $z > 2$ sources were selected based on their high $f_{esc}$, and that the absolute value of the $f_{esc}$ measurements depends on the exact assumptions made (see the text).

(Leethochawalit et al. 2016; Matthee et al. 2017). The plot does not show population-average escape fractions, which still have to be measured reliably based on large samples of galaxies with HST imaging in the future (e.g., Siana et al. 2015). In addition to our star-forming candidates, we also show the two AGNs for which we have good evidence that the ionizing photons we detect are emitted by star-forming regions (GN 21231 and GN 19591).

While the lower-redshift sources that are now being detected as LyC emitters typically still show a relatively low escape fraction of <15%, a significant fraction of the high-redshift detections reach $f_{esc} \gtrsim 50\%$. Given that most of the high-redshift points included in Figure 6 were selected based on their high $f_{esc}$, it is clear that they are not likely representative of the average galaxy at these redshifts. However, the fact that several galaxies with likely $f_{esc} > 50\%$ at $z > 2$ have been found, while no such sources have (so far) been seen at $z < 2$, hints at a possible evolution of the maximally achievable escape fraction from galaxies as a function of cosmic time (see also Inoue et al. 2006). We refer to the maximum observed $f_{esc}$ at a particular redshift as the “peak escape fraction” for that redshift.

Note that the various derivations of $f_{esc}$ in the literature use different assumptions (e.g., SED frameworks, IMF, median/mean stacking, etc.), which will affect the absolute values that are reported and shown in Figure 6. For instance, SED models that include binary stellar populations (like BPASSv2, used in this work) produce a larger number of ionizing photons and thus lead to lower $f_{esc}$ values compared with models like BC03 that have often been used in the past literature. The magnitude of this effect depends on the exact assumptions, but it is of order $\sim 2-3 \times$. This

is still smaller than the order-of-magnitude difference seen between the reported $f_{esc}$ values found for low- and high-redshift sources.

Another caveat for the above conclusion of an evolving peak escape fraction is that we do not have a complete sampling of lower-redshift LyC sources. While LyC photons can be directly observed at $z \gtrsim 1$ through UV imaging surveys, the current lower-redshift LyC emitters are all obtained through targeted, individual follow-up observations with UV spectrographs. Even though the most likely LyC candidate sources are typically followed up, it is not guaranteed that no sources with $f_{esc} > 20\%$ exist, and it will be important to continue to search for these with future observations.

### 5.4. Linking $f_{esc}$ to $z > 5$ Observables

The opacity of the IGM prevents any direct measurement of $f_{esc}$ beyond $z \gtrsim 4.5$. So in order to study $f_{esc}$ in the EoR, we need to link it to quantities that may be measured at very high redshifts. Several such indirect indicators of $f_{esc}$ have been discussed in the literature, including (1) the line ratio [O III]/[O II], which potentially traces density-bounded H II regions (e.g., Jaskot & Oey 2013; Nakajima & Ouchi 2014; Faisst 2016); (2) the strengths of nebular emission lines such as $H\beta$ compared with the total star formation rate (Zackrisson et al. 2013, 2017); (3) the shape of the Ly$\alpha$ line profile (Verhamme et al. 2015, 2017); or (4) the absorption strength of low-ionization lines and Lyman series lines, which are related to the covering fraction of absorbing gas (e.g., Heckman et al. 2011; Leethochawalit et al. 2016; Reddy et al. 2016b). With the limited data we already have on our candidates, we can discuss the first two indicators, which we do in the following sections.

#### 5.4.1. [O III]/[O II]

It has been shown that $f_{esc}$ can correlate with the oxygen line ratio [O III]/[O II], due to a higher expected [O III] flux at a given [O II] flux in density-bounded nebulae (e.g., Nakajima & Ouchi 2014). Faisst (2016) used a compilation of eight detections and four upper limits of $f_{esc}$ to show a tentative positive correlation with [O III]/[O II]. Out of these sources, Ion2 was the sole representative of the $z > 0$ universe. It is thus interesting to test whether our sources agree with this correlation.

For an escape fraction of 0.6, the relationship derived in Faisst (2016) predicts $[O III]/[O II] \sim 11$. GS 30668 displays extreme [O III] emission and may satisfy the Faisst (2016) prediction, but its [O II] flux still needs to be reliably measured. GS 14633 has a significantly lower value of $[O III]/[O II] \sim 3$, albeit with large uncertainty.

It is worthwhile to turn to Stasińska et al. (2015), who used a large sample of galaxies with extreme [O III]/[O II] and a careful analysis of photoionization models to conclude that $[O III]/[O II]$ on its own is an insufficient diagnostic tool for the leakage of LyC photons and must be used along with other lines like [Ar III], [O I], and He II and considerations of the gas covering fraction (Reddy et al. 2016b). Follow-up observations of our candidates to obtain high-resolution spectra will help make definitive statements about the $[O III]/[O II]$ approach toward constraining $f_{esc}$.

#### 5.4.2. EW($H\beta$)-$\beta_{UV}$

Zackrisson et al. (2013) show via simulations that the EW($H\beta$)-$\beta_{UV}$ diagram is an effective selector of high $f_{esc}$ at
$z > 6$, and in a follow-up study Zackrisson et al. (2017) conclude that a rest-frame EW(Hβ) < 30 Å is sufficient to select for $f_{esc} > 0.5$ at $z \sim 7-9$. Since their conclusions and diagrams only apply to dust-free SEDs with $\beta < -2.3$ (typical of $z > 6$ galaxies), we can only discuss GS 30668 ($\beta = -2.23 \pm 0.08$; $E(B-V) = 0$). If GS 30668’s Hβ flux were detected with higher certainty, we could qualitatively verify the Zackrisson et al. (2017) prediction, since both Ion2 and GS 30668 have similar metallicity, age, $\beta_{UV}$, and $f_{esc}$ and should occupy essentially the same point on the EW(Hβ)-β$_{UV}$ diagram (we have discussed the resemblance of Ion2 and GS 30668 earlier in Section 5.2). It will be important to get high-quality spectra to test such indirect methods with larger samples of directly detected LyC emitters in the future.

## 6. Summary and Outlook

In this paper we presented six galaxies that likely exhibit a large fraction of escaping ionizing photons at $z \sim 2$. These are among the first sources detected in ionizing photons at significant redshift.

The novel data that made the discovery of these candidates possible came from the HDUV survey (Oesch et al. 2016, submitted), the deepest large-area UV survey undertaken by HST to date. The F275W and F336W measurements from HDUV in combination with multiwavelength archival GOODS+CANDELS imaging provide continuous, high-resolution HST photometry from the UV to near-IR. Building on a selection method first described in Vanzella et al. (2015), we developed an SED-modeling Monte Carlo method to detect flux excesses in the UV photometry that imply LyC leakage. In our analysis we use BPASS SEDs, a newly derived dust law for the LyC region (Reddy et al. 2016a), and well-tested realizations of the IGM transmission (Inoue et al. 2014). Based on this method, we discovered six sources that have a high probability to be leaking LyC flux into the IGM, and we estimate their absolute $f_{esc}$—all very high, but ranging from >0.13 to unity (at 90% likelihood). A future paper, currently in preparation, will address the average escape fraction of galaxies in the HDUV fields. In general, our findings are in broad agreement with previous studies that found only a small fraction of galaxies to show high $f_{esc}$.

While our sources are clearly not representative of the average galaxy at these redshifts, we are finding evidence that the maximally achievable $f_{esc}$ is evolving with cosmic time. Currently, no source with $f_{esc} > 13\%$ has been found at low redshift, while several of the individual detections at $z > 2$ (including our galaxies) are consistent with $f_{esc} > 50\%$ (Figure 6).

Thanks to how richly studied the GOODS fields are, we are able to draw from existing literature and ancillary data (chiefly the 3D-HST grism survey, whose redshifts also helped in the selection) to investigate these sources in some detail. We use very deep Chandra X-ray data, Spitzer fluxes, and a Herschel study to identify three of our sources as AGNs. In two of the AGN sources, it is likely that the LyC flux nevertheless predominantly originates from star-forming regions, aided by the clearing out of the ISM by the active nucleus. In the remaining three sources the ionizing radiation is likely to originate purely from stars.

A comparison of the SFGs in our sample with the three previously known high-z LyC sources proves to be quite revealing. GS 30668 and Ion2, extreme O III emitters with two of the largest EW([O III]4959+5007 + [Hβ]) recorded at high z, resemble each other in many aspects. In general, our candidates

are ~$L_{UV}^*$ galaxies, and they show relatively young stellar population ages of ~100 Myr and little dust extinction, as has been found for previous LyC emitters.

Looking to the future, it will be important to use candidates like the ones presented here to calibrate indirect methods of estimating $f_{esc}$, since LyC photons are effectively impossible to observe beyond $z \gtrsim 4.5$. This includes relationships between $f_{esc}$ and parameters like O III/O II and Hβ/β$_{UV}$. In this work, we show tentatively that these relations have promise. High-quality, high-resolution spectra that capture features like [O II], Hβ, Si, and Lyα are required for our candidates. Studies like these will allow us to infer $f_{esc}$ for galaxies directly in the EoR with James Webb Space Telescope (JWST) NIRSPEC observations in the future to derive a self-consistent picture of cosmic reionization.

We thank the referee for a critical appraisal of this paper. R.N. was supported by Yale Astronomy's Dorrit Hoffleit Undergraduate Research Scholarship, Alice & Peter Tan, and Yale-NUS College's Summer Independent Research Program (SIRP) while working on this research. P.O. acknowledges support by the Swiss National Science Foundation through the SNSF Professorship grant 157567 "Galaxy Build-up at Cosmic Dawn."

We are grateful to Akio Inoue for providing Monte Carlo realizations of the IGM transmission at high redshift. The primary data for this work were obtained with the Hubble Space Telescope operated by AURA, Inc., for NASA under contract NAS 5-26555. Support for this work was provided by NASA through grant HST-GO-13871 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST).

Facilities: HST (ACS, WFC3).

## Appendix
### SED Fits without Emission-line-containing Bands

Figure 7 depicts the Mass-Excitation diagram described in Section 4.1, while Figures 8 and 9 support our arguments for the origins of the Lyman Continuum flux in the three AGN candidates.

<figure><img src="image_9.png" /></figure>

**Figure 7.** Mass-Excitation (MEx) diagram (Juneau et al. 2014) for star-forming LyC candidates. GS 30668 and GS 35257 have [O III] and Hβ measurements available from the 3D-HST survey. As per the MEx diagram, which separates AGNs from SFGs (separating boundary shown in red), these two sources fall in the SFG region. However, we note that the Hβ flux for these sources is <3σ detected, and future spectroscopic follow-up will help locate them on this diagram with higher certainty.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Radial profiles of LyC candidates GN 19591 (red) and GN 21231 (green). The PSF of the F275W filter is shown in blue, and $1\sigma$ errors on the normalized flux are shown by dashed lines for both the sources. These two sources are known AGNs, but their radial profiles show that their LyC flux in the F275W filter is extended and inconsistent with a pure point-source origin (i.e., the active nucleus). The FWHM of both these sources is >50% broader than that of the PSF.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as Figure 2, except the simulated color is based on AGN SED templates. GN 19913, for which we posit the LyC flux to be escaping from the active nucleus, has very low $P(f_{\text{esc}} = 0) = 5\%$, which makes it an LyC candidate with $f_{\text{esc}} = 23^{+21}_{-15}\%$. Further, the observed color for this source is consistent with the color simulations, unlike in Figure 2, which did not include AGN SEDs. The distributions shown here arise from the Stevans et al. (2014) double power law. The Siebenmorgen et al. (2015) models are also in agreement with these findings: $P(f_{\text{esc}} = 0) = 17\%$, $f_{\text{esc}} = 14^{+21}_{-21}\%$.</figcaption></figure>

## ORCID iDs

R. P. Naidu https://orcid.org/0000-0003-3997-5705

P. A. Oesch https://orcid.org/0000-0001-5851-6649

N. Reddy https://orcid.org/0000-0001-9687-4973

B. Holden https://orcid.org/0000-0002-6153-3076

C. C. Steidel https://orcid.org/0000-0002-4834-7260

H. Atek https://orcid.org/0000-0002-7570-0824

R. J. Bouwens https://orcid.org/0000-0002-4989-2471

C. M. Carollo https://orcid.org/0000-0003-1624-7609

A. Cibinel https://orcid.org/0000-0003-4578-514X

I. Labbé https://orcid.org/0000-0002-2057-5376

E. J. Nelson https://orcid.org/0000-0002-7524-374X

P. G. van Dokkum https://orcid.org/0000-0002-8282-9888

## References

Alexandroff, R. M., Heckman, T. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 810, 104

Bluck, A. F. L., Conselice, C. J., Almaini, O., et al. 2011, MNRAS, 410, 1174

Borthakur, S., Heckman, T. M., Leitherer, C., & Overzier, R. A. 2014, Sci, 346, 216

Bouwens, R. J., Illingworth, G. D., Blakeslee, J. P., & Franx, M. 2006, ApJ, 653, 53

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJL, 752, L5

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Bunker, A., Wilkins, S., Ellis, R., et al. 2010, MNRAS, 409, 855

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cappelluti, N., Comastri, A., Fontana, A., et al. 2016, ApJ, 823, 95

Castellano, M., Fontana, A., Grazian, A., et al. 2012, A&A, 540, A39

Cen, R., & Kimm, T. 2015, ApJL, 801, L25

Chabrier, G. 2003, PASP, 115, 763

Cowie, L. L., Barger, A. J., Hu, E. M., Capak, P., & Songaila, A. 2004, AJ, 127, 3137

Dawson, S., Stern, D., Bunker, A. J., Spinrad, H., & Dey, A. 2001, AJ, 122, 598

de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51

Deharveng, J.-M., Buat, V., Le Brun, V., et al. 2001, A&A, 375, 805

Dijkstra, M., Gronke, M., & Venkatesan, A. 2016, ApJ, 828, 71

Duncan, K., & Conselice, C. J. 2015, MNRAS, 451, 2030

Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37

Faisst, R. L. 2016, ApJ, 829, 99

Feng, Y., Di-Matteo, T., Croft, R. A., et al. 2016, MNRAS, 455, 2778

Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012, ApJ, 758, 93

Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83

Giavalisco, M., Dickinson, M., Ferguson, H. C., et al. 2004, ApJL, 600, L103

Grazian, A., Castellano, M., Fontana, A., et al. 2012, A&A, 547, A51

Grazian, A., Giallongo, E., Gerbasi, R., et al. 2016, A&A, 585, A48

Grimes, J. P., Heckman, T., Aloisi, A., et al. 2009, ApJS, 181, 272

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Heckman, T. M., Borthakur, S., Overzier, R., et al. 2011, ApJ, 730, 5

Inoue, A. K., Iwata, I., & Deharveng, J.-M. 2006, MNRAS, 371, L1

Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805

Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016a, Natur, 529, 178

Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683

Jaskot, A. E., & Oey, M. S. 2013, ApJ, 766, 91

Jones, T. A., Ellis, R. S., Schenker, M. A., & Stark, D. P. 2013, ApJ, 779, 52

Juneau, S., Bournaud, F., Charlot, S., et al. 2014, ApJ, 788, 88

Kirkpatrick, A., Pope, A., Alexander, D. M., et al. 2012, ApJ, 759, 139

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Köhler, M., & Li, A. 2010, MNRAS, 406, L6

Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Leethochawalit, N., Jones, T. A., Ellis, R. S., Stark, D. P., & Zitrin, A. 2016, ApJ, 831, 152

Leitet, E., Bergvall, N., Piskunov, N., & Andersson, B.-G. 2011, A&A, 532, A107

Leitherer, C., Ferguson, H. C., Heckman, T. M., & Lowenthal, J. D. 1995, ApJL, 454, L19

Leitherer, C., Hernandez, S., Lee, J. C., & Oey, M. S. 2016, ApJ, 823, 64

Madau, P., & Haardt, F. 2015, ApJL, 813, L8

Matthee, J., Sobral, D., Best, P., et al. 2017, MNRAS, 465, 3637

McLure, R. J., Dunlop, J. S., Cirasuolo, M., et al. 2010, MNRAS, 403, 960

Mitra, S., Choudhury, T. R., & Ferrara, A. 2015, MNRAS, 454, L76

Mitra, S., Choudhury, T. R., & Ferrara, A. 2016, arXiv:1606.02719

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Mostardi, R. E., Shapley, A. E., Steidel, C. C., et al. 2015, ApJ, 810, 107

Nakajima, K., Ellis, R. S., Iwata, I., et al. 2016, ApJL, 831, L9

Nakajima, K., & Ouchi, M. 2014, MNRAS, 442, 900

Nestor, D. B., Shapley, A. E., Steidel, C. C., & Siana, B. 2011, ApJ, 736, 18

Oesch, P. A., Carollo, C. M., Stiavelli, M., et al. 2009, ApJ, 690, 1350

Oesch, P. A., Montes, M., Reddy, N., et al. 2016, ApJ, submitted

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Ouchi, M., Mobasher, B., Shimasaku, K., et al. 2009, ApJ, 706, 1136

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Price, L. C., Trac, H., & Cen, R. 2016, arXiv:1605.03970
Prochaska, J. X., O'Meara, J. M., & Worseck, G. 2010, ApJ, 718, 392

Rafelski, M., Teplitz, H. I., Gardner, J. P., et al. 2015, AJ, 150, 31

Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778

Reddy, N. A., Steidel, C. C., Erb, D. K., Shapley, A. E., & Pettini, M. 2006, ApJ, 653, 1004

Reddy, N. A., Steidel, C. C., Pettini, M., & Bogosavljevic, M. 2016a, ApJ, 828, 107

Reddy, N. A., Steidel, C. C., Pettini, M., Bogosavljevic, M., & Shapley, A. 2016b, ApJ, 828, 108

Rivera-Thorsen, T. E., Hayes, M., Östlin, G., et al. 2015, ApJ, 805, 14

Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJL, 802, L19

Rutkowski, M. J., Scarlata, C., Haardt, F., et al. 2016, ApJ, 819, 81

Shapley, A. E., Steidel, C. C., Strom, A. L., et al. 2016, ApJL, 826, L24

Siana, B., Teplitz, H. I., Ferguson, H. C., et al. 2010, ApJ, 723, 241

Siana, B., Shapley, A. E., Kulas, K. R., et al. 2015, ApJ, 804, 17

Siebenmorgen, R., Heymann, F., & Efstathiou, A. 2015, A&A, 583, A120

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Smail, I., Chapman, S. C., Blain, A. W., & Ivison, R. J. 2004, ApJ, 616, 71

Stasińska, G., Izotov, Y., Morisset, C., & Guseva, N. 2015, A&A, 576, A83

Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159

Stevans, M. L., Shull, J. M., Danforth, C. W., & Tilton, E. M. 2014, ApJ, 794, 75

Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164

Teplitz, H. I., Rafelski, M., Kurczynski, P., et al. 2013, AJ, 146, 159

Trainor, R. F., Steidel, C. C., Strom, A. L., & Rudie, G. C. 2015, ApJ, 809, 89

Trump, J. R., Weiner, B. J., Scarlata, C., et al. 2011, ApJ, 743, 144

Vanzella, E., de Barros, S., Castellano, M., et al. 2015, A&A, 576, A116

Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41

Vanzella, E., Guo, Y., Giavalisco, M., et al. 2012, ApJ, 751, 70

Vanzella, E., Siana, B., Cristiani, S., & Nonino, M. 2010, MNRAS, 404, 1672

Vasei, K., Siana, B., Shapley, A. E., et al. 2016, ApJ, 831, 38

Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7

Verhamme, A., Orlitova, I., Schaerer, D., et al. 2017, A&A, 597, A13

Wilcox, R. 1998, Biometrical Journal, 40, 261

Wise, J. H., Demchenko, V. G., Halicek, M. T., et al. 2014, MNRAS, 442, 2560

Wofford, A., Charlot, S., Bruzual, G., et al. 2016, MNRAS, 457, 4296

Wofford, A., Leitherer, C., & Salzer, J. 2013, ApJ, 765, 118

Wuyts, S., Labbé, I., Förster Schreiber, N. M., et al. 2008, ApJ, 682, 985

Yoshikawa, T., Akiyama, M., Kajisawa, M., et al. 2010, ApJ, 718, 112

Zackrisson, E., Binggeli, C., Finlator, K., et al. 2017, ApJ, 836, 78

Zackrisson, E., Inoue, A. K., & Jensen, H. 2013, ApJ, 777, 39