# No strong dependence of Lyman continuum leakage on physical properties of
star-forming galaxies at $\mathbf{3.1\lesssim z\lesssim 3.5}$

A. Saxena,1 L. Pentericci,2 R. S. Ellis,1 L. Guaita,3 A. Calabrò,2 D.
Schaerer,4,5 E. Vanzella,6 R. Amorín,7,8 M. Bolzonella,9 M. Castellano,2 F.
Fontanot,10,11 N. P. Hathi,12 P. Hibon,13 M. Llerena,8 F. Mannucci,14 A.
Saldana-Lopez,4 M. Talia,15 G. Zamorani9

1Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
2INAF – Osservatorio Astronomico di Roma, via Frascati 33, 00078, Monteporzio
Catone, Italy
3Departamento de Ciencias Fisicas, Universidad Andres Bello, Fernandez Concha
700, Las Condes, Santiago, Chile
4Observatoire de Genève, Département d’Astronomie, Université de Genève, 51
Chemin Pegasi, 1290 Versoix, Switzerland
5CNRS, IRAP, 14 avenue E. Belin, 31400 Toulouse, France
6INAF – Osservatorio Astronomico di Bologna, Via P. Gobetti 93/3, 40129
Bologna, Italy
7Instituto de Investigación Multidisciplinar en Ciencia y Tecnología,
Universidad de La Serena, Raúl Bitrán, 1305 La Serena, Chile
8Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan
Cisternas 1200 Norte, La Serena, Chile
9INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via
Piero Gobetti 93/3, 40129 Bologna, Italy
10INAF –- Astronomical Observatory of Trieste, Via G.B. Tiepolo 11, 34143
Trieste, Italy
11IFPU –- Institute for Fundamental Physics of the Universe, Via Beirut 2,
34151 Trieste, Italy
12Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218, USA
13European Southern Observatory (ESO), Vitacura, Chile
14INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze,
Italy
15University of Bologna – Department of Physics and Astronomy “Augusto Righi”
(DIFA), Via Gobetti 93/2, 40129 Bologna, Italy E-mail<EMAIL_ADDRESS>
(Accepted 2021 December 19. Received 2021 November 23; in original form 2021
September 09)

###### Abstract

We present Lyman continuum (LyC) radiation escape fraction
($f_{\textrm{esc}}$) measurements for 183 spectroscopically confirmed star-
forming galaxies in the redshift range $3.11<z<3.53$ in the Chandra Deep Field
South. We use ground-based imaging to measure $f_{\textrm{esc}}$, and use
ground- and space-based photometry to derive galaxy physical properties using
spectral energy distribution (SED) fitting. We additionally derive [O iii] +
H$\beta$ equivalent widths (that fall in the observed $K$ band) by including
nebular emission in SED fitting. After removing foreground contaminants, we
report the discovery of 11 new candidate LyC leakers at $\gtrsim 2\sigma$
level, with $f_{\textrm{esc}}$ in the range $0.14-0.85$. From non-detections,
we place $1\sigma$ upper limits of $f_{\textrm{esc}}$ $<0.12$, where the
Lyman-break selected galaxies have $f_{\textrm{esc}}$ $<0.11$ and ‘blindly’
discovered galaxies with no prior photometric selection have
$f_{\textrm{esc}}$ $<0.13$. We find a slightly higher $1\sigma$ limit of
$f_{\textrm{esc}}$ $<0.20$ from extreme emission line galaxies with rest-frame
[O iii] + H$\beta$ equivalent widths $>300$ Å. For candidate LyC leakers, we
find a weak negative correlation between $f_{\textrm{esc}}$ and galaxy stellar
masses, no correlation between $f_{\textrm{esc}}$ and specific star-formation
rates (sSFRs) and a positive correlation between $f_{\textrm{esc}}$ and EW0([O
iii] + H$\beta$). The weak/no correlations between stellar mass and sSFRs may
be explained by misaligned viewing angles and/or non-coincident timescales of
starburst activity and periods of high $f_{\textrm{esc}}$. Alternatively,
escaping radiation may predominantly occur in highly localised star-forming
regions, or $f_{\textrm{esc}}$ measurements may be impacted by stochasticity
of the intervening neutral medium, obscuring any global trends with galaxy
properties. These hypotheses have important consequences for models of
reionisation.

###### keywords:

cosmology: dark ages, reionisation, first stars – cosmology: early Universe –
galaxies: high-redshift – galaxies: evolution

††pubyear: 2021††pagerange: No strong dependence of Lyman continuum leakage on
physical properties of star-forming galaxies at $\mathbf{3.1\lesssim z\lesssim
3.5}$–12

##  1 Introduction

Constraining when and how the intergalactic medium (IGM) in the Universe made
a phase transition from neutral to completely ionised, i.e. the epoch of
reionisation (EoR) is an exciting challenge in observational cosmology. An
important period in the evolutionary history of our Universe, EoR was likely
triggered by UV photon production by the first stars, galaxies and
supermassive black holes in the Universe (Bromm & Yoshida, 2011) and
observations suggest that the reionisation of the Universe unfolded over the
redshift range $6\lesssim z\lesssim 10$ (e.g. Planck Collaboration et al.,
2016b). To understand how reionisation progressed, knowledge of the intrinsic
ionising spectra of sources in the early Universe is needed, in addition to
the fraction of Hydrogen-ionising Lyman continuum (LyC; $\lambda_{0}<912$ Å)
photons that actually escape into the neutral IGM driving reionisation (see
Dayal & Ferrara, 2018, for a recent review). The rapidly increasing neutrality
of the IGM at high redshifts (e.g. Inoue et al., 2014) makes it impossible to
directly observe LyC photons from reionisation era galaxies even with
facilities like the _James Webb Space Telescope (JWST)_ and therefore, the
only possibility to obtain constraints on the escape fraction of LyC photons
($f_{\textrm{esc}}$) is by studying lower redshift analogues of reionisation
era galaxies, where escaping LyC photons are not completely absorbed by the
comparatively less neutral IGM.

Assuming a constant LyC escape fraction for galaxies across the full range of
UV luminosities or stellar masses, various works (Robertson et al. 2013;
Robertson et al. 2015, but see also Naidu et al. 2020) have determined that
$f_{\textrm{esc}}$ $\approx 10-20\%$ is required to both reionise the Universe
by $z\sim 6$ and match the Thomson optical depth of electron scattering
observed in the cosmic microwave background (Planck Collaboration et al.,
2016b). High-resolution zoom-in simulations of reionisation era galaxies have
suggested that $f_{\textrm{esc}}$ could be higher in lower mass systems, as
the energy and momentum ejected by supernova explosions are able to clear out
channels of neutral gas more efficiently in lower mass environments,
facilitating LyC photon escape (Wise et al., 2014; Trebitsch et al., 2017).
Other simulations, however, have suggested a non-monotonic dependence of
$f_{\textrm{esc}}$ on stellar/halo mass, where $f_{\textrm{esc}}$ is lower at
the low mass end ($M_{\star}<10^{8}\,M_{\odot}$) due to inefficient star
formation, peaks at $M_{\star}\sim 10^{8}\,M_{\odot}$ and decreases again at
higher masses due to dust attenuation (e.g. Ma et al., 2020).

Alternatively, a strong dependence of $f_{\textrm{esc}}$ on the star-formation
rate surface density may imply that higher mass galaxies
($>10^{8}~{}M_{\odot}$) at $z>6$ may contribute $>80\%$ of the reionisation
photon budget and also explain the rapid conclusion of reionisation (Naidu et
al., 2020). Active galactic nuclei (AGN) at $z>6$, with a considerably higher
$f_{\textrm{esc}}$ on average compared to star-forming galaxies (Cristiani et
al., 2016; Grazian et al., 2018; Romano et al., 2019), may also contribute to
the observed rapid evolution of the IGM neutral fraction at $z>6$ (Fontanot et
al., 2012; Giallongo et al., 2015; Trebitsch et al., 2018; Finkelstein et al.,
2019; Dayal et al., 2020), and contributions to the reionisation budget may be
expected from low-metallicity X-ray binaries binary populations in high
redshift galaxies (Madau & Fragos, 2017; Schaerer et al., 2019; Saxena et al.,
2021).

However, no conclusive correlations between $f_{\textrm{esc}}$ and stellar
masses or star-formation rates have yet been identified observationally. LyC
measurements presented by Izotov et al. (2021) for a sample of very low mass
($M_{\star}<10^{8}\,M_{\odot}$) galaxies at $0.30<z<0.45$ showed that
$f_{\textrm{esc}}$ does not strongly depend on the stellar mass over a broad
mass range of $10^{7}-10^{10}\,M_{\odot}$. Radiative transfer modelling has
suggested that density-bounded H ii regions within star-forming galaxies,
traced by high [O iii]/[O ii] ratios (O32; Zackrisson et al. 2013), offer
ideal environments to achieve high $f_{\textrm{esc}}$ values (Jaskot & Oey,
2013; Nakajima & Ouchi, 2014; Faisst, 2016; Izotov et al., 2018a) and although
spectroscopic follow-up of candidate LyC leakers at $z\sim 3$ has shown that
galaxies with high $f_{\textrm{esc}}$ can also show high O32 (e.g. Vanzella et
al., 2016b; Nakajima et al., 2020), several other physical phenomena that do
not require high $f_{\textrm{esc}}$ can also lead to high O32 ratios within
galaxies (e.g. Stasińska et al., 2015; Izotov et al., 2017; Paalvast et al.,
2018; Jaskot et al., 2019). Therefore, a high O32 does not also necessarily
guarantee high $f_{\textrm{esc}}$ (Izotov et al., 2018a; Naidu et al., 2018;
Bassett et al., 2019; Nakajima et al., 2020).

Further, studies that have attempted to measure $f_{\textrm{esc}}$ from
intermediate redshift galaxies using both ground-based (e.g. Grazian et al.,
2016; Guaita et al., 2016; Marchi et al., 2017a; Japelj et al., 2017; Naidu et
al., 2018) and space-based (e.g. Siana et al., 2015; Mostardi et al., 2015;
Rutkowski et al., 2017; Fletcher et al., 2019) UV imaging have generally found
low $f_{\textrm{esc}}$ values ($<10-20\%$), with no clear differences in
$f_{\textrm{esc}}$ measured from photometrically selected Lyman break galaxies
(LBGs) or galaxies selected based on strong emission lines. LyC
$f_{\textrm{esc}}$ measured from compact Lyman alpha emitting galaxies (LAEs),
considered to be direct analogues of reionisation era galaxies, have reported
a higher incidence of individual LyC leakers (with $f_{\textrm{esc}}$ $\gtrsim
10\%$; e.g. Fletcher et al. 2019), but average $f_{\textrm{esc}}$ measured
through stacking of non-leakers remains low (Mostardi et al., 2013; Fletcher
et al., 2019; Bian & Fan, 2020; Smith et al., 2020) and consistent with that
derived for other types of galaxies. The candidate LyC leaking galaxies
identified at $z\sim 2-4$, with $f_{\textrm{esc}}$ $>20\%$ (e.g. de Barros et
al., 2016; Vanzella et al., 2016a, 2018; Shapley et al., 2016; Bian et al.,
2017; Rivera-Thorsen et al., 2019; Fletcher et al., 2019; Meštrić et al.,
2020; Marques-Chaves et al., 2021), also generally occupy a broad range of
stellar masses, star-formation rates and rest-optical nebular line strengths
(e.g. [O iii]), showing no clear preference or dependence on either of these
parameters.

To fully understand the nature of galaxies that drove reionisation, thereby
dictating both the spatial and temporal evolution of the neutral fraction of
the IGM at $z>6$, it is necessary to uniformly probe the dependence of
$f_{\textrm{esc}}$ on galaxy properties using large samples of star-forming
galaxies at intermediate redshifts. To do this, a combination of accurate
spectroscopic redshifts, deep rest-frame LyC data and multi-band photometry
that enables reliable constraints on the spectral energy distributions (SEDs)
is needed. This combination is crucial to overcome uncertainties in the
measurement of both LyC $f_{\textrm{esc}}$ as well as galaxy properties such
as stellar masses, star-formation rates (SFRs), histories, ages,
metallicities, dust attenuation, etc., which may also impact the inferred
$f_{\textrm{esc}}$ (Dayal & Ferrara, 2018). Accurate redshifts can also
minimise the uncertainties in the IGM transmission along the line-of-sight
that can further affect the $f_{\textrm{esc}}$ measurement (e.g. Bassett et
al., 2021).

Therefore, in this study we construct a sample of star-forming galaxies with
high quality spectroscopic redshifts at $z=3.11-3.53$ for which deep ground-
based imaging probing rest-frame LyC emission is available. Our aim is to
uniformly measure LyC $f_{\textrm{esc}}$ and key galaxy properties using SED
fitting to explore potential correlations between them. The galaxy sample used
in this study probes a range of UV luminosities, stellar masses and star-
formation rates, containing galaxies that were photometrically selected for
spectroscopic follow-up based on their Lyman break, as well as blindly
discovered galaxies from deep imaging-spectroscopy surveys of the sky. Each
galaxy benefits from deep Hubble Space Telescope (HST) and Spitzer/IRAC
imaging of their non-ionising wavelengths, which enables highly reliable SED
fitting. The redshift range of our sample is further advantageous due to the
redshifted [O iii] + H$\beta$ emission lines falling in the observed $K$
broadband, enabling the exploration of the connection between
$f_{\textrm{esc}}$ and [O iii] + H$\beta$ line strengths.

The layout of this paper is as follows. We describe the LyC imaging and the
spectroscopic data sets used in this study, along with the sample selection
and foreground contamination identification in §2. In §3 we describe the SED
fitting procedure used to derive key physical properties for our galaxies as
well as estimate [O iii] + H$\beta$ line strengths. In §4 we describe our LyC
$f_{\textrm{esc}}$ measurements and present the results for both individual
galaxies as well as $f_{\textrm{esc}}$ measured through stacking. We also
present stacked $f_{\textrm{esc}}$ measurements for sub-samples of galaxies,
including those for galaxies that show extremely strong [O iii] + H$\beta$
emission lines. In §5 we explore the dependence of $f_{\textrm{esc}}$ on key
galaxy properties such as stellar masses, specific star-formation rates and [O
iii] + H$\beta$ strengths for both individual $f_{\textrm{esc}}$ measurements
and stacks. Finally in §6 we discuss the implications of our results and
discuss mechanisms that likely regulate $f_{\textrm{esc}}$.

Throughout this paper, we assume a $\Lambda$CDM cosmology with
$\Omega_{\textrm{m}}=0.3$ and H${}_{0}=67.7$ km s-1 Mpc-1 taken from Planck
Collaboration et al. (2016a), and use the AB magnitude system (Oke & Gunn,
1983). All logarithms are in base 10, unless otherwise specified.

##  2 Data and sample selection

In this work, we focus on the Great Observatories Origins Deep Survey-South
(GOODS-S) extragalactic field111https://www.stsci.edu/science/goods/, which
benefits from deep photometric and spectroscopic observations of galaxies from
ground- and space-based observatories across a wide wavelength range, for
example through the Cosmic Assembly Near-infrared Deep Extragalactic Legacy
Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). This field also
benefits from deep narrow-band (Guaita et al., 2016) and broad-band imaging
data (Nonino et al., 2009) probing rest-frame LyC wavelengths at $z\sim
3-3.5$, and deep _HST_ and _Spitzer_ data that enables reliable SED fitting,
crucially at wavelengths red-wards of the Balmer break at $z\sim 3.5$ (around
$4000$ Å rest-frame). We first describe the available LyC filters in this
field below, followed by the spectroscopic data used in this study, possible
foreground contaminant identification and final galaxy sample creation.

###  2.1 LyC filters – narrow bands and U-band

To measure LyC signal, typically measured in the rest-frame wavelength range
$860<\lambda<912$ Å at $z>3$, we use narrow-band imaging data from Guaita et
al. (2016) and the publicly available $U$-band imaging data (Nonino et al.,
2009). We briefly describe the properties of these LyC filters below.

The images from the narrow bands (NB) _NB3727_ and _NB396_ were first
presented in Guaita et al. (2016) and we refer the readers to this paper for
full details of the observations, data reduction, calibration and other
properties associated with these images. Briefly, the _NB3727_ image was
obtained using the MOSAIC II instrument at the Cerro Tololo Interamerican
Observatory (CTIO), with the exposure time totalling to 36 hours, a total
coverage area of $31.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$ $\times$
$31.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$ on the sky and $1\sigma$
detection limit of 27.38 AB. _NB3727_ offers LyC coverage in the wavelength
range mentioned above for redshifts $3.11<z<3.34$.

The _NB396_ image was obtained using the WFI instrument mounted on the ESO
2.2m Telescope at La Silla, with a total exposure time of 13 hours, covering a
subset of the sky area encompassed by the _NB3727_ image with a $1\sigma$
detection limit at 27.72 AB. _NB396_ samples LyC signal in the range
$3.41<z<3.53$.

We also use the deep $U$-band image in the GOODS-S field taken by Nonino et
al. (2009) using VIMOS on the Very Large Telescope (VLT) covering an area of
$\approx 630$ arcmin2 and reaching a $1\sigma$ magnitude limit of 29.78 AB.
The $U$-band filter covers rest-frame LyC signal in the redshift range
$3.34<z<3.49$. Although flux bluewards of the Lyman break is sampled for
redshifts above $3.49$, the declining IGM transmissivity would strongly
suppress the expected flux at these shorter wavelengths (e.g. Inoue et al.,
2014) and may bias the measured LyC fluxes. The transmission curves and the
redshift range over which the rest-frame LyC flux is sampled by each filter
are shown in Figure 1.

Figure 1: Filter transmission profiles of the _NB3727_ , _NB396_ and the _U_
band filters used to target LyC emission in this study. The x-axis shows the
redshift at the rest-frame wavelength of 886 Å, which is the mid-point
wavelength of the LyC emission sampled between $860<\lambda<912$ Å in this
study (shaded regions for each filter). More details about the filter
properties and LyC redshift ranges covered are given in §2.1.

Therefore, owing to the rest-frame LyC wavelength coverage offered by imaging
data, the redshift range probed in this study is $3.11<z<3.53$. We now present
the spectroscopic data within the GOODS-S/CDFS field used to construct a
sample of star-forming galaxies in this redshift range.

###  2.2 Spectroscopic data

####  2.2.1 VANDELS

The majority of spectra used in this study are from VANDELS – a deep VIMOS
survey of the CANDELS UDS and CDFS fields – which is a recently completed ESO
public spectroscopic survey carried out using the VLT. The VANDELS survey
covers two well-studied extragalactic fields, the UKIDSS Ultra Deep Survey
(UDS) centred around RA = 02:17:38, Dec = $-$05:11:55, and the Chandra Deep
Field South (CDFS/GOODS-S) centred around RA = 03:32:30, Dec = $-$27:48:28. As
mentioned earlier, we use VANDELS spectra in the CDFS/GOODS-S field in this
study. The survey description and target selection can be found in McLure et
al. (2018), and detailed data reduction and redshift determination steps can
be found in Pentericci et al. (2018).

VANDELS spectroscopic targets were selected based on photometric redshifts (or
the presence of Lyman breaks in their broadband photometry) derived from deep
ground and space-based multi-band photometry. The parent spectroscopic sample
included bright star-forming galaxies ($i\leq 25.0$ AB) in the range $2.4\leq
z\leq 5.5$ with 99% of this subsample having specific star-formation rates
(sSFR) $>0.6$ Gyr-1, fainter Lyman-break galaxies (LBGs; $i\leq 27.5$ AB) in
the redshift range $3.0\leq z\leq 7.0$ with sSFR $>0.3$ Gyr-1 and passive
galaxies with sSFR $<0.1$ Gyr-1 (McLure et al., 2018). Over the redshift range
of interest for this study, the VANDELS galaxies are largely sampled from the
star-forming main sequence.

The latest VANDELS data release, DR4 (Garilli et al., 2021)
222http://vandels.inaf.it/dr4.html contains spectra of $\sim 2100$ galaxies in
the redshift range $1.0<z<7.0$, with on-source integration times ranging from
20 to 80 hours, where $>70\%$ of the targets have at least 40 hours of
integration time. The survey depths have been chosen so that the final spectra
have high signal-to-noise ratios (SNRs), enabling detailed absorption and
emission line studies, the study of ionising sources using emission line
ratios, the derivation of accurate metallicities and better constraints on
physical parameters such as stellar masses and star-formation rates. The
spectral resolution of the final spectra is $R\sim 600$.

The spectra used in this study were reduced using the Easylife data reduction
pipeline (Garilli et al., 2012) and the final data products delivered by
VANDELS include the extracted 1D spectra, the 2D re-sampled and sky-subtracted
spectra and catalogues with physical parameters, including spectroscopic
redshifts. The reliability of redshifts in the VANDELS database is recorded
using the following flags: 0 – no redshift could be assigned, 1 – 50%
probability to be correct, 2 – 70-80% probability to be correct, 3 – 95-100%
probability to be correct, 4 – 100% probability to be correct and 9 – spectrum
shows a single emission line. The typical accuracy of spectroscopic redshift
measurements is $\sim 150$ km s-1. Details about the final VANDELS data
products, redshift determination processes and redshift reliability flags
assignment can be found in Pentericci et al. (2018).

In the redshift range $3.11<z<3.53$, we identified 90 galaxies with
spectroscopic flags $3$ and $4$, which guarantees a $>95\%$ probability that
the spectroscopic redshift is correct and reliable. The VANDELS team has
performed aperture matched, PSF-homogenised photometry for every source
identified in the CANDELS footprint in GOODS-S across available imaging data
ranging from UV to NIR wavelengths (McLure et al., 2018), which we refer to as
the ‘master’ catalogue. Accurate and homogenised photometry is essential for
robust spectral energy distribution (SED) fitting to determine key galaxy
properties and we rely on photometry from this master catalogue for this
study.

Therefore, moving forward we require every source selected from other publicly
available spectroscopic data that are described below to have a matched entry
in the master catalogue and below we list the other publicly available data
sets from which sources with reliable spectroscopic redshifts were queried.

####  2.2.2 VUDS

The VIMOS Ultra-Deep Survey (VUDS; Le Fèvre et al.
2015)333http://cesam.lam.fr/vuds/ is a survey of $\sim 10000$ faint galaxies
in the redshift range $2<z\leq 6$ with $i$-band magnitudes down to $i\approx
27.0$ AB. The spectroscopic targets were also selected based on either
reliable photometric redshifts or the presence of a Lyman break (LBGs) at
higher redshifts, similar to VANDELS.

From the VUDS survey, we select sources in the GOODS-S/CDFS field with
redshift reliability greater than 95%, where the redshift reliability flags
are determined in a similar manner to VANDELS. We refer the reader to Le Fèvre
et al. (2015) and Tasca et al. (2017) for more information about source
selection, data reduction and redshift measurements. In total we select 20
sources with reliable spectroscopic redshifts from VUDS.

####  2.2.3 VIMOS10

The publicly available VLT/VIMOS spectroscopic survey in GOODS-S (Balestra et
al. 2010; VIMOS10) contains spectra obtained using both the Low Resolution and
Medium Resolution grisms for galaxies up to $z\sim 3.5$ with $B$ and $R$
magnitudes $\sim 24-25$ AB. The spectroscopic targets were predominantly LBGs,
similar to VANDELS and VUDS. We only use sources that have ‘secure’ redshifts,
i.e., $>99\%$ confidence of being correct owing to the relatively lower
integration times of this survey, and find that 43 sources satisfied our
selection criteria.

####  2.2.4 MUSE-WIDE Survey

In addition to photo-z or LBG selected spectroscopic targets, we query
reliable spectroscopic redshifts from the publicly available MUSE-WIDE survey
in the GOODS-S/CDFS field444https://musewide.aip.de/project/, which is a
‘blind’ integral-field unit (IFU) survey carried out without photometric pre-
selection (Urrutia et al., 2019). The unique aspect of surveys enabled by
MUSE555https://www.eso.org/sci/facilities/develop/instruments/muse.html on the
VLT is that no photometric pre-selection is required prior to obtaining a
spectroscopic redshift. Therefore, spectroscopically confirmed sources are
often blindly selected, owing to their emission lines alone and not
necessarily their photometric redshifts or Lyman (or Balmer) breaks, which may
bias the target selection. The MUSE-WIDE Survey is capable of obtaining
spectroscopic redshifts for $98\%$ of photometrically selected galaxies
brighter than $24.0$ AB in HST/F775W. This blind galaxy sample is therefore
highly complementary to the photo-$z$ or LBG selection implemented by other
data sets used in this study.

From this survey, we select sources with redshifts $3.11<z<3.53$ and reliable
spectroscopic redshifts, i.e. a ‘confidence’ value of 2 or 3, corresponding to
an expected error rate of $\lesssim 10\,\%$ in the correct identification of
emission and absorption features (Urrutia et al., 2019). We find that 44
galaxies in the MUSE-WIDE data release satisfy our selection criteria.

####  2.2.5 MUSE-HUDF Survey

We also use the published spectroscopic redshifts from the MUSE deep (coverage
area of $3^{\prime}\times 3^{\prime}$) and ultra deep (coverage area of
$1^{\prime}\times 1^{\prime}$) surveys in the Hubble Ultra Deep Field (HUDF;
Inami et al. 2017)666http://muse-vlt.eu/science/udf/, which is part of the
GOODS-S field. This survey is deeper than MUSE-WIDE and is capable of blindly
detecting sources as faint as $\approx 30.0$ AB in HST/F775W (Inami et al.,
2017).

We find that 42 sources from MUSE-HUDF lie in the redshift range of interest
and have reliable spectroscopic redshifts, represented by confidence levels 2
and 3 and these sources are added to our sample.

####  2.2.6 3D-HST

Additional redshifts were queried from the 3D-HST survey, which is a
spectroscopic survey conducted using the WFC3/G141 as well as the ACS/G800L
grisms on the _HST_ covering $\sim 600$ square arcminutes of well-studied
extragalactic fields, including GOODS-S (Brammer et al., 2012; Momcheva et
al., 2016). The 3D-HST survey obtained blind spectroscopy of galaxies, with
the redshift determination primarily relying on the [O ii] emission line in
the redshift range of interest to this study. Since there was no photometric
pre-selection required, the 3D-HST survey is similar to the MUSE surveys.

We add a further 9 sources to our study from 3D-HST that have “good” redshifts
(use_subj = 1) from the the publicly available
catalogues777https://3dhst.research.yale.edu/Data.php.

####  2.2.7 Source duplication and final sample

In total we shortlist 248 sources across the above-mentioned surveys for this
study, targeting galaxies based on photometric redshifts of Lyman-break
selection, which we refer to as the Lyman break or LBG-selected, and those
that were ‘blindly’ detected without any photometric pre-selection, which we
refer to as the ‘blind selected’ sample. We note here again that although the
spectroscopic surveys have different target selection strategies, we
eventually require a detection in the PSF-homogenised CANDELS master
photometric catalogue in GOODS-S for every source going forward.

Due to overlapping areas covered by the various spectroscopic surveys, a low
fraction sources have been observed by multiple surveys. In such cases we give
the highest precedence to the redshift obtained by VANDELS. In cases where the
multiply detected source was not part of VANDELS, we give preference to the
redshift recorded by the survey that is most sensitive out of the overlapping
surveys.

Removing duplicate observations, 219 unique sources remain in our final sample
and we show the breakdown of sources from various spectroscopic surveys in
Table 1. We now turn our attention to filtering out sources that may suffer
contamination to their LyC flux by emission from foreground sources along the
line-of-sight using high-resolution HST photometry, as well as identifying
possible AGN in the sample using deep X-ray data.

Table 1: Number of sources with spectroscopic redshifts from data sets used in
this study, at every stage of the selection described in §2.

Data set | $N$ | Reference
---|---|---
_Lyman-break (LBG) selected_ |  |
VANDELS DR4 | 90 | (1, 2)
VUDS | 20 | (3, 4)
VIMOS10 | 43 | (5)
_Blind selected_ |  |
MUSE-WIDE | 44 | (6)
MUSE HUDF | 42 | (7)
3D-HST | 9 | (8)
Total | 248 | $-$
Total unique | 219 | $-$
Possible contaminants | 34 | $-$
Possible AGN | 2 | $-$
Final sample | 183 | $-$
Total LBG selected | 111 | $-$
Total blind selected | 72 | $-$

References: (1) McLure et al. (2018); (2) Pentericci et al. (2018); (3) Le
Fèvre et al. (2015); (4) Tasca et al. (2017); (5) Balestra et al. (2010) (6)
Urrutia et al. (2019); (7) Inami et al. (2017); (8) Momcheva et al. (2016).

###  2.3 Identifying possible foreground contaminants

Figure 2: Ground-based LyC and HST postage stamps ($5^{\prime\prime}\times
5^{\prime\prime}$) for a target galaxy that is isolated (top) and a target
galaxy that has a possible foreground contaminant within $1\text{\,}\arcsec$
(bottom). For galaxies that are likely contaminated by foreground objects,
examining the colours of both the target galaxy and the contaminant aids
considerably in identifying photometric contamination in the LyC bands by a
foreground galaxy (see Figure 3) that is often mistakenly interpreted as LyC
signal (e.g. Mostardi et al., 2015; Grazian et al., 2016; Guaita et al.,
2016). Figure 3: Colour image ($5^{\prime\prime}\times 5^{\prime\prime}$)
using the F606W (blue), F775W (green) and F850LP (red) _HST_ filters of the
central and companion galaxy of the example shown in the bottom panel of
Figure 2. The colour of the companion is clearly different from the central
galaxy that is of interest, which indicates that the companion source is
likely a foreground object.

One of the biggest sources of error in the measurement of LyC signal from
intermediate redshift galaxies is contamination within the aperture from
emission from foreground galaxies along the line-of-sight to the source of
interest (e.g. Vanzella et al., 2010). Therefore, before proceeding to measure
LyC fluxes it is important to identify and remove sources that suffer from
possible contamination to their LyC signal in the form of image artefacts,
bright neighbouring sources or low-redshift interlopers along the line-of-
sight. We adopt a series of steps, similar to Guaita et al. (2016), to filter
the sample in order to maximise the accuracy of the LyC signal measurement.

Since all sources in our sample lie within the footprint of the GOODS-S field,
they benefit from deep _HST_ photometry at optical and near-infrared
wavelengths from the _HST_ Treasury Program (Dickinson et al., 2003) and
through CANDELS (Grogin et al., 2011; Koekemoer et al., 2011), offering
unparalleled resolution. This resolution is critical in the identification of
possible foreground neighbouring sources in the vicinity of the galaxies with
spectroscopic redshifts (e.g. Mostardi et al., 2015) that may strongly
contaminate any LyC signal observed from the relatively lower resolution
ground-based images.

To identify potential contamination, we create $20\text{\,}\arcsec$ $\times$
$20\text{\,}\arcsec$ postage stamps from images in the F435W, F606W, F775W and
F850LP _HST_ filters using publicly available
images888https://archive.stsci.edu/pub/hlsp/goods/v2/h_goods_v2.0_rdm.html. We
also create postage stamps of the same dimensions in the LyC images. Sources
that appear to be isolated in the _HST_ images, without any other component or
a nearby source within $1\text{\,}\arcsec$ of their centroid coordinates are
shortlisted for LyC signal measurement without any further filtration. An
example of such a source (with postage stamp sizes of $5^{\prime\prime}\times
5^{\prime\prime}$ for clarity) is shown in the top panel of Figure 2.

For sources that have one or more components within a radius of
$1\text{\,}\arcsec$ of the centroid in the _HST_ images, we must determine
whether these are a part of the galaxy in question (which may have a ‘clumpy’
morphology) or whether they are possibly foreground objects (e.g. Grazian et
al., 2016). In the bottom panel of Figure 2 we show an example (with
$5^{\prime\prime}\times 5^{\prime\prime}$ postage stamps for clarity) where
the spectroscopically confirmed source is surrounded by a nearby source that
lies within $1\text{\,}\arcsec$. Following Guaita et al. (2016), we examine
the colours of the central galaxy and the surrounding components to identify
contaminants in the following way. Using the F435W ($B$), F606W ($V$), F775W
($i$) and F850LP ($z$) images, we generate BGR colour images using the
make_lupton_rgb function of astropy.visualization in $BVi$, $BVz$, $BIz$ and
$VIz$. Additionally, we also produce a $V-z$ image for all sources. By then
comparing the colours of the spectroscopically confirmed target and its
companion(s) within $1\text{\,}\arcsec$, we determine whether the component(s)
close to the centroid of the source is (are) part of the system or not. For
example, Figure 3 shows the $Viz$ colour image for the source shown in the
bottom panel of Figure 2, which is indicative that the companion source is
likely to be a foreground contaminant.

Guaita et al. (2016) showed that the above-mentioned simpler and more
computationally efficient colour-based foreground identification method is
highly effective at identifying interlopers. Such an analysis also returns a
comparable contamination identification fraction compared to a full SED-
fitting based analysis of the individual components surrounding the
spectroscopically confirmed source. We also stress that removing sources that
show clumps of different colours is a highly conservative approach, as it
might result in the removal of galaxies with clumpy morphologies where the
individual clumps may harbour slightly different stellar populations.

Through this process, we remove 34 sources from our initial sample of 220
unique sources, representing $\approx 15\%$ of the sample and bringing the
sample size of unique sources that are unlikely to suffer from foreground
contamination to 186. The fraction of sources removed due to possible
contamination by foreground objects is highly consistent with an expected
fraction of $12.6\%$ for $z\gtrsim 3$ galaxies determined by Vanzella et al.
(2010), down to the photometric depths considered in this study assuming a
$1\text{\,}\arcsec$ radius for contaminant identification.

###  2.4 Removing possible AGN

Finally, we remove potential AGN from our sample by cross-matching the
coordinates of our sources with the _Chandra_ 7 Ms source catalogue compiled
by Luo et al. (2017) using a matching radius of $1.0\text{\,}\arcsec$. The
sources detected in the Luo et al. (2017) catalogue at $z\gtrsim 3$ are very
likely to be AGN, owing to X-ray luminosities of $>10^{42}$ erg s-1 at these
redshifts (e.g. Magliocchetti et al., 2020; Saxena et al., 2020a, 2021). We
identify and remove two sources with counterparts in the X-ray catalogue that
are also classified as AGN in the catalogue and therefore, the total number of
sources in our final sample comes to 183.

###  2.5 Final sample

After removing duplicates, sources with foreground contamination and AGN from
our sample, the total number of unique galaxies in the redshift range
$3.11<z<3.53$ with reliable spectroscopic redshifts is 183.

The median spectroscopic redshift of our sample is $z=3.38$ and our sample
includes the known over-density of star-forming galaxies at $z\sim 3.5$ in
GOODS-S reported by Forrest et al. (2017). The median HST/F606W band magnitude
of the galaxies in the final sample is $25.69\pm 0.84$. We now proceed to fit
spectral energy distribution (SED) models to the photometric data available to
derive key galaxy physical parameters.

##  3 SED fitting and physical properties

Since the main aim of this study is to explore the dependence of LyC
$f_{\textrm{esc}}$ on key galaxy properties such as stellar masses, star-
formation rates, etc. in addition to statistically exploring the correlations
between $f_{\textrm{esc}}$ and the strength of rest-frame nebular emission
lines tracing intense star-formation such as [O iii] + H$\beta$, obtaining
robust and reliable spectral energy distribution (SED) fits using high quality
photometry of galaxies in our sample is crucial.

These key galaxy properties are potentially important factors that govern
$f_{\textrm{esc}}$ from reionisation era galaxies, and in this section we use
deep _HST_ and _Spitzer_ /IRAC photometry in the observed frame optical and
near infrared bands, complemented by deep ground based $K$ band photometry in
combination with reliable spectroscopic redshifts to accurately constrain the
stellar population properties and dust extinction for the 183 sources in our
final sample. We additionally compare the derived physical properties of
galaxies in our sample that were LBG selected, and those that were blind
selected, exploring any inherent differences in galaxies that these two
selection techniques may introduce, thereby impacting the inferred
$f_{\textrm{esc}}$.

###  3.1 SED derived physical properties

To determine both the properties of stellar populations as well as the [O iii]
+ H$\beta$ line strengths using $K$ band photometry, we simultaneously fit the
stellar continuum and nebular line emission to the observed broadband
photometry of galaxies in our sample. We note that since both stellar
continuum and nebular emission are being simultaneously fitted, we include the
line-contaminated $K$ band photometry in the fitting.

Figure 4: Distribution of SED derived stellar masses and star-formation rates
for galaxies in the final sample, colour coded by metallicity. The dashed line
shows the star-forming main sequence derived by Santini et al. (2017) at
redshifts $3<z<4$. We find that most of our galaxies lie on the main sequence,
and span a range of masses and star-formation rates, which makes them
representative of the general galaxy population at $z\sim 3.3$

The SED fitting is performed using the Bayesian Analysis of Galaxies for
Physical Inference and Parameter EStimation (Bagpipes) code (Carnall et al.,
2018). Bagpipes is a python based tool to generate complex galaxy model
spectra, allowing flexible fitting to photometric data to derive key galaxy
physical properties. As mentioned earlier, we use a combination of stellar
continuum and nebular emission when fitting the SEDs. The stellar models are
based on Bruzual & Charlot (2003), constructed using the Kroupa & Boily (2002)
initial mass function (IMF). The nebular emission is computed using Cloudy
grids (Ferland et al., 2013), with the only free parameter being the log of
the ionisation parameter ($\log U$). When determining the strengths of the
emission lines from these grids, the metallicity of the line-emitting gas is
assumed to be the same as that of the stars.

We assume an exponentially declining star-formation history ($\tau$-model),
with stellar ages allowed to vary between 10 Myr and the age of the Universe
at the observed redshift of the galaxy. The value of $\tau$ is allowed to vary
between 100 Myr to 10 Gyr, with stellar metallicities allowed to vary in the
range $0-1\,Z_{\odot}$, where $Z_{\odot}$ is the solar metallicity. The
Calzetti et al. (2000) dust attenuation law is applied and the dust
extinction, $A_{1400}$ is allowed to vary between 0 and 2 mag. We use flat
posterior distributions for each of the parameters, which is the default
choice in Bagpipes. The nebular emission line strengths are set by the $\log
U$ parameter, which we allow to vary from $-3$ to $-2$, which are the
typically inferred values for $z\gtrsim 3$ star-forming galaxies (see Nakajima
et al., 2018; Saxena et al., 2020b, for example).

The chi-squared values of the best-fitting SEDs are evaluated, and individual
SEDs are then visually inspected to ensure that the broadband photometry, and
crucially the $K$ band flux, are well reproduced by the best-fitting model.
The best-fitting SEDs are then used to infer key physical properties for each
galaxy, such as stellar age, dust attenuation, stellar mass, metallicity,
star-formation rates, etc., in addition to the fluxes and rest-frame
equivalent widths of the [O iii] $\lambda 5007$ and the H$\beta$ lines derived
from the nebular component of the best-fitting SED. Thanks to the Bayesian
inference methodology of Bagpipes, posterior distributions for all physical
parameters are derived, leading to accurate uncertainty estimates.

The strength of the nebular emission that is added to the stellar continuum
using Bagpipes crucially depends on the best-fitting metallicity as well as
the ionisation parameter (see Carnall et al., 2018). Although the ionisation
parameter is a variable that can be controlled, the best-fitting metallicities
are inferred from the broadband SED fitting using both stellar and nebular
components.

The median stellar mass (with $1\sigma$ standard deviation) derived from SED
fitting of all 183 galaxies is $\log(M_{\star}/M_{\odot})=9.2\pm 0.5$, the
median star-formation rate is $8\pm 24$ M${}_{\odot}\,\textrm{yr}^{-1}$ , the
median stellar age is $1.5\pm 1.4$ Gyr, and the median stellar metallicity is
$0.6\pm 0.3$ $Z_{\odot}$, where $Z_{\odot}$ is the solar metallicity (Asplund
et al., 2009). Overall, the population of galaxies considered in this study
are consistent with low-mass, star-formation objects with young stellar ages
and low stellar metallicities.

In Figure 4 we show the distribution of the derived stellar masses and star-
formation rates of all the galaxies in the final sample. We mark the LBG-
selected galaxies in our sample using blue points, and the emission line
selected galaxies using orange squares. The histograms of both stellar masses
and SFRs for the LBG and blind selected samples are also shown. We show the
star-forming main sequence at $3<z<4$ determined by Santini et al. (2017)
using dashed lines.

The majority of galaxies in our sample lie on or around the main sequence, and
that the stellar masses and star-formation rates vary over a large range. We
also find that as expected, the blind selected galaxies tend to occupy the
lower stellar mass end of the main-sequence, consequently also exhibiting
lower star-formation rates compared to LBG-selected galaxies.

We note here, however, that owing to the magnitude limited nature of our
galaxy sample, a larger fraction of galaxies with lower stellar masses
detected in the spectroscopic surveys will cluster towards higher SFRs, as is
evident from Figure 4, thereby resulting in an overall higher sSFR probed at
lower stellar masses. Such an effect can potentially bias any possible
dependence of LyC $f_{\textrm{esc}}$ on key galaxy properties such as stellar
mass and/or sSFR.

###  3.2 [O iii] + H$\beta$ line strengths

The best-fitting SED includes a stellar continuum and a nebular emission line
component, from which the equivalent width of the emission lines can be
calculated. Particularly for the [O iii] and H$\beta$ lines that fall in the
observed $K$ broad band, we measure a median EW0([O iii] + H$\beta$) of
$\approx 230$ Å with standard deviation $225$ Å, with the LBG selected sample
having a median EW0([O iii] + H$\beta$) $=236.2$ Å and the blind selected
sample having a slightly lower median EW0([O iii] + H$\beta$) $=207.2$ Å. The
distribution of EW0([O iii] + H$\beta$) measured for LBGs and blind selected
galaxies is shown in Figure 5.

Figure 5: Histogram of the rest-frame [O iii] + H$\beta$ equivalent widths
measured for LBG-selected galaxies (blue) and blind selected galaxies (orange)
in our sample using SED fitting. The median EW0([O iii] + H$\beta$) of the
full sample is $228.6$ Å, and the derived equivalent widths follow a log-
normal distribution (similar to that reported by Stark et al. 2013. The EW
distributions for both LBGs and emission line galaxies are comparable, with
the median values (dashed lines and text) highly consistent with each other.

The most immediate comparison for our [O iii] + H$\beta$ equivalent widths can
be made with the study of Forrest et al. (2017), where using $K$ band
photometry for galaxies at $2.95<z<3.65$ they report an average EW0([O iii] +
H$\beta$) $=230\pm 90$ Å, which is highly comparable to the median value we
find for our sample. In Figure 6 we compare our SED derived line measurements
with the best-fitting relations derived between spectroscopic measurements of
the [O iii] $\lambda 5007$ line and the stellar masses of a large sample of
star-forming galaxies at $z\sim 2$ as part of the MOSDEF survey (Reddy et al.,
2018). We additionally show the stellar mass dependence of the [O iii] line
derived for a sample of extreme emission line galaxies (defined as EW0([O
iii]) $>225$Å by Tang et al. 2019) at $z\sim 2$ by Tang et al. (2019), which
also exhibit lower stellar masses.

We find that the distribution of EW0([O iii]) as a function of stellar mass
for galaxies in our sample broadly agrees with the best-fitting relation
derived from MOSDEF at higher stellar masses. The scatter around the best-
fitting relation observed in our sample is also consistent with what was
reported by Reddy et al. (2018).

We identify 59 galaxies in our sample that can be classified as extreme
emission line galaxies (EELGs), with EW0([O iii] + H$\beta$) $>300$ Å (e.g.
Stark, 2016)). The median EW0([O iii] + H$\beta$) in the EELG sample is
$\approx 500$ Å, which is slightly lower than the range of $\approx 660-800$ Å
reported by Forrest et al. (2017) for EELGs in their sample.

To further investigate the reliability of our inferred [O iii] + H$\beta$ line
strengths, we compare our rest-frame [O iii] + H$\beta$ equivalent width
measurements as a function of galaxy stellar to results from the MOSDEF
spectroscopic survey (Kriek et al., 2015), which spectroscopically observed a
sample of $\sim 1200$ star-forming galaxies at redshifts $1.4\lesssim
z\lesssim 3.8$ (Reddy et al., 2018) in Figure 6. We find that our measurements
at $z\sim 3-3.5$ broadly agree with the best-fitting relation obtained by
Reddy et al. (2018), with comparable scatter at stellar masses $>10^{9}$
$M_{\odot}$.

We also compare EWs of galaxies in our sample that may be considered as EELGs
following the Tang et al. (2019) classification (EW0([O iii]) $>225$Å) with
spectroscopic measurements at $z\sim 2$ from Tang et al. (2019), finding broad
agreement with their best-fitting relation (Tang et al., private
communication), which is also shown in Figure 6.

Figure 6: Distribution of EW0([O iii]) and galaxy stellar mass derived for our
sample of galaxies using SED fitting. Shown for comparison are the best-
fitting relations derived between the two quantities for relatively massive
star-forming galaxies at $z\sim 2$ by Reddy et al. (2018) (red solid line)
from the MOSDEF survey and extreme emission line galaxies (EELGs) at $z\sim 2$
by Tang et al. (2019) (blue dot-dashed line), where the black dashed line
marks EW0([O iii]) $>225$Å, used by Tang et al. (2019) to classify galaxies as
EELGs. The normalisation of the best-fitting relations have not been altered.
Galaxies from our sample that would be considered as normal emission line
galaxies (ELGs) are shown as red points, and those considered as EELGs are
shown using blue points. The distribution of our SED derived measurements (and
the scatter) are broadly consistent with spectroscopic measurements at lower
redshifts.

Our EELG [O iii] + H$\beta$ EWs are also comparable with photometrically
selected extreme line emitters found at $z>6$ – Endsley et al. (2020) measured
a median value of EW0([O iii] + H$\beta$) $\sim 700$ Å for star-forming
galaxies at $z\sim 7$ and Labbé et al. (2013) found an average rest frame [O
iii] + H$\beta$ EW of $670\pm 200$ Å at $z\sim 8$. Additionally, Smit et al.
(2015) identified galaxies with extreme Spitzer/IRAC colours at $z\sim 6.8$
and found that 50% of the strongest line emitters have EW0([O iii] + H$\beta$)
$\sim 1085$ Å. Using a similar technique, Roberts-Borsani et al. (2016)
identified four bright galaxies at $z\sim 7$ with EW0([O iii] + H$\beta$)
$\sim 1500$ Å, all of which have been confirmed to be strong Ly$\alpha$
emitters. Further, Castellano et al. (2017) measured an average EW([O iii] +
H$\beta$) of $\sim 1500$ Åfor Ly$\alpha$ emitting galaxies at $z\sim 6.8$.
Therefore, our sample of EELGs at $z\sim 3-3.5$ offers an excellent
opportunity to study the escaping ionising radiation from analogues of
reionisation-era galaxies.

We note here, however, that the line strengths of some of the brightest [O
iii] + H$\beta$ emitting galaxies, particularly those with fainter continua,
may be underestimated due to the lower limit on the ionisation parameter of
$-2.0$ permitted by Bagpipes as well as the depth of photometric data in the
$K$ band. We further note that the flux (and consequently stellar mass)
limited nature of our galaxy sample naturally results in incompleteness of
EELGs down to lower stellar masses.

###  3.3 Comparison with NIR spectroscopic observations

To test the accuracy of the inferred [O iii] + H$\beta$ strengths from SED
fitting, we compare our measurements with near-IR observations of a subset of
VANDELS galaxies in the CDFS and UDS fields in the redshift range $z\sim 3-4$
using MOSFIRE on Keck Telescope, as part of the recently completed NIRVANDELS
survey (Cullen et al., 2021). We note that only a fraction of NIRVANDELS
galaxies in the CDFS field are also part of the final sample in this study.

To do this, we fit SEDs to all NIRVANDELS galaxies using the same SED fitting
parameters and techniques that we used for the main galaxy sample described in
§3 and compare the observed and SED-derived [O iii] line fluxes. The
comparison between the SED derived and observed [O iii] line fluxes for
NIRVANDELS galaxies is shown in Figure 7, with the dashed line marking the
one-to-one relation. The error bars shown on the SED derived line fluxes are
at the 20% level. We find that the SED derived line fluxes are remarkably
close to the observed values (within $\sim 0.3$ dex) over a broad flux range.
We therefore conclude that simultaneously fitting the stellar continuum and
nebular line emission with the ionisation parameter in the range $-3\leq\log
U\leq-2$ gives realistic estimates of the [O iii] (and H$\beta$) line fluxes
for galaxies in our sample.

Figure 7: Comparison of [O iii] + H$\beta$ line fluxes measured from the best-
fitting SED and $K$-band excess and the observed line fluxes from NIR spectra.
From the limited number of galaxies in our sample that have observed line
fluxes from Cullen et al. (2021), we can conclude that our method of
estimating line fluxes from $K$-band excess works reasonably well.

Having fitted SEDs to all 183 unique galaxies in our final sample and measured
key galaxy properties as well as [O iii] and H$\beta$ line strengths, we now
proceed to measure LyC fluxes for galaxies using the available LyC images.

##  4 Lyman continuum escape fraction

###  4.1 Measuring $f_{\textrm{esc}}^{\textrm{rel}}$ and $f_{\textrm{esc}}$

In this section we describe our methodology for measuring the LyC flux from
the available LyC images. Since these filters cover slightly different rest-
frame LyC wavelengths at different redshifts as mentioned in §2.1, we split
our sample into the relevant redshift bins corresponding to LyC coverage
offered by the filters. NB3727 includes 79 galaxies, NB396 includes 83
galaxies, and U includes 91 galaxies. We show the number of galaxies and the
redshift ranges corresponding to each LyC filters along with $1\sigma$ depths
in Table 2. We note that owing to the overlapping wavelength coverage offered
by the the NB396 and the U-band filter, some galaxies are included in both
subsets.

Table 2: LyC filters used in this study along with $1\sigma$ depths (see §2.1), LyC redshift range covered and the number of sources. Note that some sources have overlapping LyC coverage in NB396 and U images. Filter | Tel./Instr. |  $1\sigma$ depth (AB) |  $z$ range | N | Reference
---|---|---|---|---|---
_NB3727_ | CTIO 4m/MOSAICII | 27.38 | $3.11-3.34$ | 79 | Guaita et al. (2016)
_NB396_ | La Silla 2.2m/WFI | 27.72 | $3.41-3.53$ | 83 | Guaita et al. (2016)
$U$ band | VLT/VIMOS | 29.78 | $3.34-3.49$ | 91 | Nonino et al. (2009)

An important parameter quantifying the escape of LyC photons from a galaxy in
a dust-independent manner is the relative $f_{\textrm{esc}}$, or
$f_{\textrm{esc}}^{\textrm{rel}}$ (Steidel et al., 2001), comparing the
observed ratio of ionising to non-ionising UV flux to the intrinsically
produced ratio by stellar populations within the galaxy, corrected for
absorption by the intervening neutral IGM. If the intrinsic ratio of
luminosities at $\sim 1400$ Å and $\sim 900$ Å,
$(L_{1400}/L_{900})_{\textrm{int}}$ is known (e.g. from stellar population
modelling), then $f_{\textrm{esc}}^{\textrm{rel}}$ can be calculated as

$f_{\textrm{esc}}^{\textrm{rel}}=\left(\frac{F_{900}}{F_{1400}}\right)_{\textrm{obs}}\times\left(\frac{L_{1400}}{L_{900}}\right)_{\textrm{int}}\times\frac{1}{\exp\left(-\tau_{900}\right)}$
(1)

where ${(F_{900}/F_{1400})_{\textrm{obs}}}$ is the observed ratio of the
ionising LyC flux to non-ionising UV flux at $\sim 1400$ Å, and
$\tau^{\textrm{IGM}}_{900}$ is the transmissivity of ionising photons (at
$900$ Å) through the intervening IGM (Siana et al., 2007). The absolute escape
fraction, $f_{\textrm{esc}}$, can then be calculated from
$f_{\textrm{esc}}^{\textrm{rel}}$ by applying a dust correction as follows

$f_{\textrm{esc}}=f_{\textrm{esc}}^{\textrm{rel}}\times 10^{(-0.4\times
A_{1400})}$ (2)

where $A_{1400}=10.33\times E(B-V)$ according to the Calzetti et al. (2000)
dust reddening law. Therefore, to calculate $f_{\textrm{esc}}$ we must
constrain each of these ingredients in equations (1) and (2), and below we
highlight how these are calculated for our galaxies.

####  4.1.1 LyC flux measurement

The observed flux ratios in the equations above require a measurement of the
LyC flux, which we measure from the available LyC images using aperture
photometry. Photometry in the LyC images is performed by placing circular
apertures at the coordinates in the CANDELS catalogues using the python
package photutils (Bradley et al., 2020). The average PSF of the _NB3727_
image is $1.4\text{\,}\arcsec$, of the _NB396_ is $0.94\text{\,}\arcsec$
(Guaita et al., 2016) and of the _U_ image is $0.8\text{\,}\arcsec$ (Nonino et
al., 2009). The radius of the circular aperture is set to either match or be
slightly larger than the PSF to capture small-scale PSF variations across the
LyC images. Therefore, we use an aperture with radius $0.7\text{\,}\arcsec$
for NB3727 and $0.5\text{\,}\arcsec$ for NB396 and U band images.

We measure the local background noise around the apertures by placing a
circular annulus (at the same positions as the aperture) with an inner radius
of $2\text{\,}\arcsec$ and an outer radius of $3\text{\,}\arcsec$. For
galaxies in which the flux measured in the aperture was lower than the
background flux measured in the annulus, we assigned a limiting LyC flux value
equal to the local background noise with a 100% error.

The non-ionising UV flux is calculated from the F606W HST image, which covers
the rest-frame $\sim 1400$ Å wavelengths for the redshift range of galaxies in
this study. Therefore, the ratio of the flux determined from the LyC bands and
from the F606W band gives us a measure of the term
$(F_{900}/F_{1400})_{\textrm{obs}}$.

####  4.1.2 Intrinsic non-ionising to ionising luminosity ratio

The intrinsic ratio between UV luminosity at non-ionising wavelengths (e.g.
$\sim 1400$ Å) and LyC wavelengths ($\sim 900$ Å) is calculated individually
for each source using the best-fitting SED. We create synthetic SEDs using the
best-fitting stellar population parameters determined and measure the
intrinsic ratios, $L_{1400}/L_{900}$, for each galaxy.

We find an average intrinsic ratio of $4.5$, with a standard deviation of
$0.7$. Individual intrinsic ratio measurements for each galaxy enabled by
accurate spectroscopic redshifts lead to a more accurate $f_{\textrm{esc}}$
calculation for individual galaxies compared assuming an average intrinsic
ratio for all galaxies in the sample, as done by previous studies measuring
LyC $f_{\textrm{esc}}$ at $z\sim 3$, where the assumed intrinsic ratio took
values of 3 (Steidel et al., 2001; Grazian et al., 2016, 2017; Marchi et al.,
2017b), 5 (Guaita et al., 2016; Vanzella et al., 2016a) or 7 (Fletcher et al.,
2019).

####  4.1.3 IGM transmissivity

We calculate the IGM transmissivity of LyC photons for each galaxy using the
prescriptions of Steidel et al. (2018). In addition to the IGM transmission,
the Steidel et al. (2018) model also takes into account the effects of
transmission through the circumgalactic medium (CGM) surrounding each galaxy
in the so-called “IGM+CGM” model. This results in an average LyC transmission,
$\exp(-\tau_{900})=0.371$ at $z=3.0$ and a much lower
$\exp(-\tau_{900})=0.264$ at $z=3.5$.

An advantage of having spectroscopic redshifts is that the IGM transmission
can be accurately calculated for each source at the rest-frame wavelength of
$900$ Å, the $\tau_{900}$. Hence, we calculate the LyC transmission through
both the IGM and the CGM, or $\exp(-\tau_{900})$ for each galaxy in our final
sample.

####  4.1.4 Dust attenuation

Finally, the dust attenuation (or $A_{1400}$) for each source is calculated
from the best-fitting SED. During the SED fitting process we allowed
$A_{1400}$ to vary between 0 and 2, which naturally limits the derived values.
This means that our SED fits will miss any significantly dusty galaxies, but
given the high-redshift star-forming nature of these galaxies, heavy dust
obscuration is unlikely to be the case. We find an average dust attenuation of
$A_{1400}=0.40$ mag in our sample with a standard deviation of $0.26$ mag.

Having defined the necessary ingredients to calculate $f_{\textrm{esc}}$ for
each source and described the measurements of each of these parameters, in the
following section we present the our measurements of the relative and the
absolute escape fractions of LyC photons from each galaxy.

###  4.2 Individual LyC $f_{\textrm{esc}}$ measurements

We measure non-zero LyC $f_{\textrm{esc}}$ at $\gtrsim 2\sigma$ significance
levels in 11 out of 183 galaxies in our sample. The LyC $f_{\textrm{esc}}$ is
in the range of $0.14-0.85$, and the galaxies showing non-zero
$f_{\textrm{esc}}$ span a wide range of stellar masses
($M_{\star}=10^{8.3}-10^{9.7}\,M_{\odot}$), star-formation rates (SFR
$=2.3-37.4\,M_{\odot}$ yr-1) and [O iii] + H$\beta$ equivalent widths (EW0([O
iii] + H$\beta$) $=111.8-1160.3$ Å).

We note that 8 out of 11 candidate LyC leakers were LBG selected, and 3 were
blindly selected for spectroscopic follow-up. The galaxy with the highest
measured $f_{\textrm{esc}}$ $=0.85\pm 0.33$ (CDFS 13385) was blindly selected
from the MUSE-HUDF survey and has a relatively low stellar mass of
$10^{8.3}\,M_{\odot}$, showing extreme [O iii] + H$\beta$ line emission
through its $K$ band photometry.

The best-fitting SED derived physical parameters for these 11 candidate LyC
leaking galaxies are given in Table 3, and the measured quantities (both
through SED fitting and aperture photometry) used to calculate
$f_{\textrm{esc}}^{\textrm{rel}}$ and $f_{\textrm{esc}}$ are given in Table 4.
The $5^{\prime\prime}\times 5^{\prime\prime}$ cutout galleries showing the LyC
image and the observed-frame optical HST photometry, along with the best-
fitting SEDs for the 11 LyC leaker candidates are shown in Appendix A. A full
spectroscopic analysis of these 11 LyC leaking candidates is beyond the scope
of the present paper and will be explored in greater detail in forthcoming
publications.

Table 3: Nature of selection for each galaxy, with LBG representing Lyman-break selected and BLD representing blind selection, as well as best-fitting SED derived physical properties for candidate LyC leaking galaxies. Also given is the LyC detection filter for each source. CANDELS ID | Sample | Survey | $z$ | LyC filter | log($M_{\star}$/$M_{\odot}$) | SFR ($M_{\odot}$ yr-1) | log(sSFR/yr-1) | EW0 ([O iii] + H$\beta$) (Å)
---|---|---|---|---|---|---|---|---
CDFS 5161 | LBG | VANDELS | 3.420 | NB396 | 9.3 | 16.3 | $-8.1$ | 268.2
CDFS 9358 | LBG | VIMOS10 | 3.229 | NB3727 | 9.6 | 37.0 | $-8.1$ | 241.3
CDFS 9692 | LBG | VANDELS | 3.470 | NB396 | 8.3 | 2.3 | $-8.0$ | 1160.3
CDFS 12448 | LBG | VANDELS | 3.460 | U | 9.6 | 12.5 | $-8.5$ | 119.3
CDFS 13385 | BLD | MUSE HUDF | 3.431 | U | 8.6 | 3.6 | $-8.1$ | 334.6
CDFS 15718 | BLD | MUSE HUDF | 3.439 | NB396 | 8.3 | 2.0 | $-8.0$ | 986.8
CDFS 16444 | BLD | MUSE HUDF | 3.128 | NB3727 | 9.3 | 23.1 | $-7.9$ | 456.4
CDFS 18454 | LBG | VUDS | 3.163 | NB3727 | 9.2 | 16.7 | $-8.0$ | 339.8
CDFS 19872 | LBG | VANDELS | 3.452 | NB396 | 9.6 | 5.6 | $-8.8$ | 111.8
CDFS 20745 | LBG | VANDELS | 3.495 | NB396 | 9.0 | 5.5 | $-8.3$ | 400.6
CDFS 24975 | LBG | VANDELS | 3.187 | NB3727 | 9.7 | 19.8 | $-8.5$ | 214.9
Table 4: Parameters required for the calculation of $f_{\textrm{esc}}^{\textrm{rel}}$ and $f_{\textrm{esc}}$ for candidate LyC leakers identified in this study. CANDELS ID | $z$ | $(L_{1400}/L_{900})_{\rm{int}}$ |  $A_{V}$ (mag) |  $F_{900}$ ($\mu$Jy) | $(F_{900}/F_{1400})_{\rm{obs}}$ | <exp($-\tau_{900}$)> | $f_{\textrm{esc}}^{\textrm{rel}}$ | $f_{\textrm{esc}}$
---|---|---|---|---|---|---|---|---
CDFS 5161 | 3.420 | 3.7 | $0.6\pm 0.1$ | $0.036\pm 0.015$ | $0.07\pm 0.02$ | 0.28 | $0.91\pm 0.39$ | $0.53\pm 0.24$
CDFS 9358 | 3.229 | 3.4 | $0.5\pm 0.1$ | $0.048\pm 0.020$ | $0.02\pm 0.01$ | 0.32 | $0.21\pm 0.11$ | $0.14\pm 0.07$
CDFS 9692 | 3.470 | 2.3 | $0.4\pm 0.1$ | $0.039\pm 0.020$ | $0.06\pm 0.03$ | 0.27 | $0.55\pm 0.25$ | $0.38\pm 0.19$
CDFS 12448 | 3.460 | 4.1 | $0.1\pm 0.1$ | $0.006\pm 0.001$ | $0.02\pm 0.01$ | 0.27 | $0.24\pm 0.10$ | $0.20\pm 0.12$
CDFS 13385 | 3.431 | 3.5 | $0.3\pm 0.1$ | $0.009\pm 0.002$ | $0.09\pm 0.02$ | 0.28 | $1.12\pm 0.24$ | $0.85\pm 0.33$
CDFS 15718 | 3.439 | 2.0 | $0.4\pm 0.1$ | $0.034\pm 0.009$ | $0.06\pm 0.02$ | 0.28 | $0.42\pm 0.12$ | $0.28\pm 0.10$
CDFS 16444 | 3.128 | 3.0 | $0.7\pm 0.1$ | $0.068\pm 0.014$ | $0.06\pm 0.01$ | 0.34 | $0.55\pm 0.11$ | $0.30\pm 0.07$
CDFS 18454 | 3.163 | 2.2 | $0.8\pm 0.1$ | $0.045\pm 0.021$ | $0.04\pm 0.02$ | 0.34 | $0.42\pm 0.21$ | $0.20\pm 0.10$
CDFS 19872 | 3.452 | 4.0 | $0.1\pm 0.1$ | $0.040\pm 0.020$ | $0.02\pm 0.01$ | 0.28 | $0.32\pm 0.16$ | $0.29\pm 0.30$
CDFS 20745 | 3.495 | 4.0 | $0.2\pm 0.1$ | $0.033\pm 0.015$ | $0.03\pm 0.02$ | 0.27 | $0.45\pm 0.23$ | $0.38\pm 0.30$
CDFS 24975 | 3.187 | 3.4 | $0.7\pm 0.1$ | $0.043\pm 0.021$ | $0.05\pm 0.02$ | 0.33 | $0.53\pm 0.21$ | $0.27\pm 0.11$

It is important to note that given the stochasticity of the IGM attenuation at
these redshifts (see for example Bassett et al., 2021), it is hard to
accurately determine the absolute LyC $f_{\textrm{esc}}$ for sources at
$z\gtrsim 3$. Therefore, the IGM stochasticity introduces an unavoidable
element of scatter in the $f_{\textrm{esc}}$ measurements, which in turn may
be an important source of scatter in the correlations between
$f_{\textrm{esc}}$ and key galaxy properties that will be investigated in the
sections that follow. We will discuss the effects of IGM stochasticity on the
inferred LyC $f_{\textrm{esc}}$ measurements in §6.

Given the relatively large sample size of 183 sources, we expect $\approx 4$
sources to appear as noise fluctuations at the $2\sigma$ level. From Table 4
we note that 4 out of 11 sources identified as candidate LyC leakers are close
to the $2\sigma$ confidence level in their respective LyC images and could be
consistent with noise fluctuations. Deeper photometric data at LyC wavelengths
is needed to confirm the LyC leaking nature of these sources.

Apart from the 11 candidate LyC leakers, all other galaxies in our sample
remain undetected at $\gtrsim 2\sigma$ levels in the available LyC images,
indicating that the global $f_{\textrm{esc}}$ from star-forming galaxies at
$z=3-3.5$ must be low. In the sections that follow we derive the global
$f_{\textrm{esc}}$ from different populations of galaxies in our sample
through stacking the low-significance LyC fluxes measured from the majority of
galaxies, and use these stacked measurements to set additional constraints on
any global dependence of $f_{\textrm{esc}}$ on key galaxy properties traced by
sub-populations of galaxies within our sample.

###  4.3 Stacked LyC $f_{\textrm{esc}}$ measurements

Only $\approx 6\%$ galaxies in our sample have a $\gtrsim 2\sigma$ detection
in the LyC images. Therefore, to constrain the global $f_{\textrm{esc}}$ we
must co-add the LyC signal from galaxies that were not classified as candidate
LyC leakers. To do this, we first remove the 11 individual LyC leakers from
the stacking sample and employ both the stacking of LyC images, as well as co-
adding the LyC flux measured through aperture photometry as described below.

We first investigate whether any LyC signal can be measured through stacking
the LyC images of individual galaxies falling within the redshift coverage of
each LyC filter. To do this, we produce $10^{\prime\prime}\times
10^{\prime\prime}$ cutouts of both the LyC image and the corresponding root-
mean-squared (RMS) error image centred on the CANDELS coordinates for galaxies
that are not individual LyC leakers. We then apply $5\sigma$ clipping to the
cutouts to masks bright pixels that could bias the stacked image. We then co-
add the sigma-clipped cutouts using a weighted averaging method, where the
corresponding RMS of each pixel of the galaxy cutout is used to assign a
weight to the pixels of each galaxy cutout. In Figure 8 we show the results
from stacking the LyC cutouts for each of the filters and as can be seen, we
do not find any statistically significant LyC signal even in the stacked
_NB3727_ , _NB396_ and _U_ band images, consistent with that reported by
Guaita et al. (2016) from the two narrow bands and by Naidu et al. (2018) from
the $U$ broad band.

Figure 8: Stacked image from all sources in _NB3727_ (left), _NB396_ (middle)
and _U_ band (right). We do not see any statistically significant flux in any
of the stacked images.

Instead of extracting a flux from these stacked LyC images and relying on the
sample averaged values for the intrinsic non-ionising to ionising UV ratio,
IGM transmissivity and dust attenuation to calculate $f_{\textrm{esc}}$, we
use a slightly different approach and derive global $f_{\textrm{esc}}$ using
the standard deviation of the weighted average of the individual LyC fluxes
measured through aperture photometry across the available LyC images as
described in §4.1.1 (see also Guaita et al., 2016). The weight assigned to the
$f_{\textrm{esc}}$ measurement for each source is the inverse square of the
error on $f_{\textrm{esc}}$, or $1/\Delta f_{\rm{esc}}^{2}$.

After background subtraction, the weighted average of sources that are not
classified as candidate LyC leakers would be expected to be close to zero,
with the standard deviation providing a robust $1\sigma$ upper limit on the
LyC signal from the sample (see also Guaita et al., 2016). This method takes
full advantage of the availability of spectroscopic redshifts and accurate SED
fitting that enable reliable measurements of both the values, errors and
limits on individual $f_{\textrm{esc}}^{\textrm{rel}}$ and $f_{\textrm{esc}}$
for each galaxy.

Galaxies in the redshift range covered by NB3727 do not overlap with the other
two LyC filters. However, galaxies with redshifts between $3.41-3.49$ have
overlapping LyC coverage in both NB396 and U images. Therefore, to assign a
LyC flux and $f_{\textrm{esc}}$ to each unique galaxy in our sample, we
prioritise the flux measured from the LyC image with a lower error, which
translates to a lower error on $f_{\textrm{esc}}$. We note here that
$f_{\textrm{esc}}$ measurements/limits recorded using both LyC filters are
broadly in agreement with each other.

Using the standard deviation of the weighted average of LyC flux from the
unique sample of galaxies that are not individually identified as LyC leakers,
we obtain a $1\sigma$ upper limit of $f_{\textrm{esc}}$ $<0.12$. We
additionally calculate the stacked $f_{\textrm{esc}}$ limit obtained from all
galaxies in the sample (i.e. both LyC detections and non-detections) finding a
slightly higher $1\sigma$ limit of $f_{\textrm{esc}}$ $<0.16$, which is not
surprising. This stacked measurement from the full sample potentially helps
minimise the effects of stochastic IGM transmission along the line of sight,
which we discuss further in §6.

####  4.3.1 LyC $f_{\textrm{esc}}$ from LBG and blind selected galaxies

Out of the 172 galaxies that were not identified as individual LyC detections,
103 are LBG selected and 69 are blindly selected. Using weighted averaging, we
derive $1\sigma$ $f_{\textrm{esc}}$ limits of $<0.11$ for the LBG selected and
a slightly higher $1\sigma$ limit of $<0.13$ for the blind selected sample. We
note here that since the blind selected sample is also fainter in UV
magnitudes, the flux limited nature of the available LyC images may bias the
inferred $f_{\textrm{esc}}$ to higher values.

A LyC $f_{\textrm{esc}}$ $<0.11$ for our LBG selected sample is comparable
with $f_{\textrm{esc}}$ $\approx 0.09$ reported by Steidel et al. (2018) by
stacking deep UV spectra of $z\sim 3$ LBGs and $f_{\textrm{esc}}$ $\approx
0.06$ obtained by a later reanalysis of the Steidel et al. (2018) sample by
Pahl et al. (2021), and with other measurements from LBGs at $z\sim 3$ (e.g.
Grazian et al., 2016; Guaita et al., 2016; Marchi et al., 2017a; Naidu et al.,
2018; Meštrić et al., 2021). The slightly higher $1\sigma$ $f_{\textrm{esc}}$
limit of $f_{\textrm{esc}}$ $<0.13$ measured for blind-selected galaxies is
also consistent with measurements of emission-line selected galaxies at
$z=2-4$ (e.g. Mostardi et al., 2015; Japelj et al., 2017; Ji et al., 2020;
Smith et al., 2020).

####  4.3.2 LyC $f_{\textrm{esc}}$ from extreme emission line galaxies (EELGs)

In this section we investigate the stacked $f_{\textrm{esc}}$ measurements of
the previously identified extreme emission line galaxies (EELGs; EW0([O iii] +
H$\beta$) $\geq 300$ Å) as well as normal emission line galaxies (ELGs; EW0([O
iii] + H$\beta$) $<300$ Å) in our sample (see §3.2). For this analysis we
include the 11 candidate LyC leakers to better understand the dependence of
$f_{\textrm{esc}}$ on [O iii] + H$\beta$ EWs.

There are 59 EELGs and 124 ELGs in our sample. The average EW0([O iii] +
H$\beta$) for EELGs is $421.2$ Å and the average EW0([O iii] + H$\beta$) for
ELGs is $183.5$ Å. By once again calculating the weighted average of
$f_{\textrm{esc}}$ for ELGs and EELGs, we determine $1\sigma$ upper limits of
$f_{\textrm{esc}}$$<0.14$ for ELGs and $f_{\textrm{esc}}$$<0.20$ for EELGs. We
note here that we obtain slightly higher $f_{\textrm{esc}}$ limits from bins
created based on emission line strengths compared to the LBG and blind-
selected samples because the individually identified LyC leakers are included
in the ELG and EELG bins.

Qualitatively speaking, the distribution of [O iii] + H$\beta$ EWs measured
for individual LyC leakers is similar to the distribution of the parent
sample, implying that the individual LyC leakers seem to have been randomly
drawn from the parent sample at least in terms of [O iii] + H$\beta$ EWs. This
is consistent with what was reported by Nakajima et al. (2020) for a sample of
LAEs at $z\sim 3$, where the [O iii] + H$\beta$ EWs of LAEs that are LyC
leakers and non-leakers are not significantly different.

In Table 5 we give the median galaxy physical properties and the measured
$f_{\textrm{esc}}$ for subsets created based on the initial spectroscopic
selection method, as well as on the strength of [O iii] + H$\beta$ emission
lines. No clear differences in $f_{\textrm{esc}}$ are seen between these
subsets, posing questions surrounding the origin of leaking LyC radiation from
only a low fraction of galaxies at intermediate redshifts. In the following
section we investigate the dependence of $f_{\textrm{esc}}$ on key galaxy
properties for both candidate LyC leakers and stacks of galaxies.

Table 5: Physical properties and $f_{\textrm{esc}}$ for LBG or blind selected galaxies, ELGs and EELGs, and the full sample. Note that the LBG and blind selected subsets are created only for those galaxies that are undetected in the LyC images at $\gtrsim 2\sigma$ level, but the ELG and EELG subsets include the individual LyC leakers as well. Subset | $N$ | log($M_{\star}$/$M_{\odot}$) | log(sSFR/yr-1) | EW0 ([O iii] + H$\beta$) (Å) | $f_{\textrm{esc}}$
---|---|---|---|---|---
LBG selected, LyC non-detections | 103 | $9.4$ | $-8.2$ | $236.2$ | $<0.11$
Blind selected, LyC non-detections | 69 | $9.0$ | $-8.3$ | $207.2$ | $<0.13$
ELGs (all) – EW0([O iii] + H$\beta$) $<300$ Å | 124 | $9.3$ | $-8.3$ | 183.5 | $<0.14$
EELGs (all) – EW0([O iii] + H$\beta$) $\geq 300$ Å | 59 | $8.9$ | $-8.0$ | 421.2 | $<0.20$
Full sample | 183 | $9.2$ | $-8.3$ | 230.9 | $<0.16$

##  5 Results

Having measured LyC fluxes and $f_{\textrm{esc}}$ values for every source in
the sample, and having used stacking to establish $f_{\textrm{esc}}$ limits
from the sample of galaxies that do not show any LyC leakage, in this section
we explore potential correlations between $f_{\textrm{esc}}$ and the
properties of the underlying stellar population in the galaxies. Constraining
the dependence of $f_{\textrm{esc}}$ on these key galaxy properties is
important to understand the key drivers of reionisation at $z>6$, and here we
rely on robust SED-fitting based measurements to explore these.

###  5.1 Dependence of $f_{\textrm{esc}}$ on stellar mass

We first explore the dependence of $f_{\textrm{esc}}$ on the stellar masses of
galaxies. In Figure 9 we show the $f_{\textrm{esc}}$ measured for the
candidate LyC leaking candidates as well as the stacked limits obtained for
the sub samples of LBG and blind selected galaxies in our sample. As reported
in §3.1, the median stellar mass of LBG selected galaxies is
$10^{9.3}\,M_{\odot}$ and that of the blind selected galaxies is
$10^{9.0}\,M_{\odot}$.

We find that the candidate LyC leaking galaxies occupy a large range of
stellar masses, ranging from $10^{8.3}-10^{9.7}\,M_{\odot}$. The galaxy with
the highest $f_{\textrm{esc}}$ in our sample, CDFS 13385 at $z=3.431$ with
$f_{\textrm{esc}}$ $=0.85\pm 0.33$ has a stellar mass of $\approx
10^{8.3}\,M_{\odot}$, which is lower than the stellar masses of known LyC
leaking galaxies with $f_{\textrm{esc}}$$>0.5$ at $z\gtrsim 3$ such as _Ion2_
(de Barros et al., 2016; Vanzella et al., 2016a), _Ion3_ (Vanzella et al.,
2018) and Q1549$-$C25 (Shapley et al., 2016), and comparable to the stellar
masses of LyC leaking LAEs identified by Fletcher et al. (2019). The LyC
leaker with the highest stellar masses in our sample has $M_{\star}\approx
10^{9.7}$ – CDFS 24975 at $z=3.187$ with $f_{\textrm{esc}}$ $=0.27\pm 0.11$.

We use least-squares fitting to obtain the best-fitting parameters for the
function $f_{\textrm{esc}}$ $=A\log(M_{\star}/M_{\odot})+B$, finding
$A=-0.10\pm 0.02$ and $B=1.17\pm 0.40$ (solid green line in Figure 9). The
slope of the best-fitting relation suggests a very weak dependence of
$f_{\textrm{esc}}$ on the stellar mass of galaxies in our sample. Although not
very strong, the mass dependence of $f_{\textrm{esc}}$ that we report from our
candidate LyC leakers is stronger than the no apparent mass dependence
reported by Izotov et al. (2021) for a sample of very low mass
($M_{\star}<10^{8}$ $M_{\odot}$) galaxies in the redshift range $z=0.31-0.45$.

We also find that there is no clear difference in the upper limits of
$f_{\textrm{esc}}$ obtained by stacking samples of LBG-selected and blind-
selected galaxies. Blind-selected galaxies on average have lower stellar
masses, and a slightly higher $f_{\textrm{esc}}$ upper limit, which is
consistent with the mass dependence obtained by fitting measurements for the
candidate LyC leakers in our sample.

Figure 9: LyC $f_{\textrm{esc}}$ as a function of stellar mass for candidate
LyC leakers in our sample (green points), and stacked measurements from LBG
selected (blue arrow), blind-selected galaxies (orange arrow) and the sample
average (black arrow) in this study. The number of galaxies used in the
stacking measurement are also given. We find a very weak correlation between
$f_{\textrm{esc}}$ and stellar mass, with a best-fitting slope $-0.10\pm 0.02$
(solid green line, uncertainty in shaded region). This dependence on stellar
mass is stronger that the recently observed lack of correlation between the
two quantities for a sample of low mass ($<10^{8}\,M_{\star}$) galaxies at
$0.3\lesssim z\lesssim 0.45$ by Izotov et al. (2021), however the scatter is
very large. We note that the stacked limits are not included when obtaining
the best-fitting relations in this study.

###  5.2 Dependence of $f_{\textrm{esc}}$ on specific SFR

We now investigate correlations between $f_{\textrm{esc}}$ and specific SFR
(sSFR), which is the star-formation rate per unit stellar mass of a galaxy.
The distribution of $f_{\textrm{esc}}$ and sSFR for both the candidate LyC
leakers and stacks of LBG and blind selected galaxies is shown in Figure 10.
As reported in §3.1, the blind selected galaxies in our sample have a slightly
lower median specific SFR of log(sSFR/yr-1) $=-8.3$ compared to log(sSFR/yr-1)
$=-8.2$ for LBG selected galaxies.

The candidate LyC leakers cover a large range of sSFRs, from log(sSFR/yr-1)
$\approx-8.8$ to log(sSFR/yr-1) $\approx-7.9$. We note that the majority of
LyC leaking candidates (seven out of 11) have relatively high sSFR values,
with log(sSFR/yr-1) $\gtrsim-8.1$, comparable to LyC leaking LAEs from
Fletcher et al. (2019). In particular, the strongest LyC leaker in our sample,
CDFS 13385, has log(sSFR/yr-1) $=-8.1$, which is lower than _Ion2_ and _Ion3_
but more than an order of magnitude higher than Q1549$-$C25 (Shapley et al.,
2016).

To quantify the dependence of $f_{\textrm{esc}}$ on sSFR, we again use least-
squares fitting and obtain best-fitting correlations to the function
$f_{\textrm{esc}}$ $=C\log(\rm{sSFR}/\rm{yr}^{-1})+D$ with a slope of
$C=-0.028\pm 0.015$ (solid green line in Figure 10). This best-fitting
suggests a very weak decrease of $f_{\textrm{esc}}$ with increasing sSFR, but
is consistent with no dependence of $f_{\textrm{esc}}$ on sSFR at the $\approx
2\sigma$ confidence level, consistent with what has been reported from disc
galaxy simulations at high redshifts (e.g. Yoo et al., 2020).

Assuming that the stellar mass of a galaxy is closely linked to its size at UV
wavelengths, no strong dependence of $f_{\textrm{esc}}$ on sSFR may also
translate to no strong dependence of $f_{\textrm{esc}}$ on the star-formation
surface density of galaxies (e.g. Model II of Naidu et al., 2020). This may
have bearing on the dominant contribution of the UV brightest sources towards
reionisation, and accurate size measurements of our candidates LyC leakers
that will be explored in a future study will help robustly test the dependence
of $f_{\textrm{esc}}$ on SFR surface densities.

Figure 10: Distribution of $f_{\textrm{esc}}$ and sSFR of candidate LyC
leakers, LBG and blind-selected galaxies as well as the sample average (same
legend as Figure 9). We find little to no correlation between
$f_{\textrm{esc}}$ and sSFR of individual LyC leakers at the $2\sigma$
confidence level, with a best-fitting slope $-0.028\pm 0.015$ (solid green
line with uncertainty in shaded region). The stacked $f_{\textrm{esc}}$ values
derived from LBG and blind-selected samples are consistent with the best-
fitting relation. The lack of dependence of $f_{\textrm{esc}}$ on sSFR has
also been seen in simulations (e.g. Yoo et al., 2020).

###  5.3 Dependence of $f_{\textrm{esc}}$ on [O iii] + H$\beta$ strengths

In this section we investigate the dependence of $f_{\textrm{esc}}$ on the [O
iii] + H$\beta$ line strengths for galaxies in our sample. In Figure 11 we
show the $f_{\textrm{esc}}$ and the log of EW0([O iii] + H$\beta$) for
individual LyC leakers (green points) as well limits derived for normal
emission line galaxies (ELGs; red arrow) and and extreme emission line
galaxies (EELGs; brown arrow), which were defined in §4.3.2. We note that 6
out of the 11 individual LyC leakers can be classified as EELGs, with EW0([O
iii] + H$\beta$) $>300$ Å.

For the individual LyC leaking candidates, CDFS 9692 at $z=3.470$ with
$f_{\textrm{esc}}$ $=0.38\pm 0.19$ has the highest EW0([O iii] + H$\beta$)
value of $1160.3$ Å. The strongest LyC leaker, CDFS 13385, also shows extreme
[O iii] + H$\beta$ line fluxes, with EW0([O iii] + H$\beta$) $=986.8$ Å, in
line with expectations from galaxies with high $f_{\textrm{esc}}$. However,
overall we do not find consistently high [O iii] + H$\beta$ EWs across our
candidate LyC leaking galaxies. Similarly, as reported in §4.3.2 we do not
find a considerably higher $f_{\textrm{esc}}$ using stacking for EELGs
compared to the sample of ELGs, but do find that the $f_{\textrm{esc}}$ limit
for EELGs is higher than that for ELGs.

A best-fitting relation between $f_{\textrm{esc}}$ and [O iii] + H$\beta$ line
strengths shows a positive correlation, indicating that galaxies with higher
$f_{\textrm{esc}}$ also tend to show stronger [O iii] + H$\beta$ lines, shown
in Figure 11. The best-fitting relation has a positive slope of $0.15\pm
0.04$, indicative of a mild dependence of $f_{\textrm{esc}}$ on rest-frame [O
iii] + H$\beta$ EW, with the general trend of increasing $f_{\textrm{esc}}$ at
higher [O iii] + H$\beta$ strengths consistent with what has been seen in the
local Universe (e.g. Izotov et al., 2018a). The higher $f_{\textrm{esc}}$
upper limit of $<0.20$ seen for EELGs compared to $f_{\textrm{esc}}$ $<0.14$
seem for ELGs is also consistent with the best-fitting relation obtained.

NIR spectroscopy of LyC leaking LAEs from the Fletcher et al. (2019) sample
presented by Nakajima et al. (2020), however, did not find uniformly high [O
iii] + H$\beta$ line strengths for LyC leakers, which is consistent with our
findings. Izotov et al. (2018a) reported a large spread in the
$f_{\textrm{esc}}$ values of local LyC leaking galaxies with high [O iii] line
strengths. Therefore, there is increasing evidence that suggests that high [O
iii] + H$\beta$ line strengths may be a necessary condition for
$f_{\textrm{esc}}$, but are by no means sufficient to guarantee high
$f_{\textrm{esc}}$. This is reflected in our galaxy sample as well, where
galaxies with extremely high [O iii] + H$\beta$ strengths do not all
necessarily appear to be strong LyC leakers.

Figure 11: Distribution of $f_{\textrm{esc}}$ and log(EW0([O iii] +
H$\beta$)/Å) for individual sources (green points) and limits derived for ELGs
(red arrow) and EELGs (brown arrow) in our sample. We find a positive
correlation between $f_{\textrm{esc}}$ and the strength of the [O iii] +
H$\beta$ lines for the candidate LyC leaker candidates, with best-fitting
slope of $0.15\pm 0.04$. High [O iii] + H$\beta$ line fluxes are seen when the
physical conditions within a galaxy favour low metallicities or intense star-
formation activity, which are important conditions to both produce a large
number of LyC photons as well as facilitate a high escape fraction of these
photons.

We also note here that given the volume limited nature of the spectroscopic
galaxy survey data used in this study, and the rarity on the sky of extreme
line emitting galaxies, we are likely to miss several extreme [O iii] +
H$\beta$ line emitters in our sample. Studies investigating correlations
between $f_{\textrm{esc}}$ and strength of [O iii] + H$\beta$ emission (or [O
iii]/[O ii] ratios) have largely focused on galaxy samples at relatively low
redshifts, and the lack of completeness of extreme line emitters that are
likely to exhibit high $f_{\textrm{esc}}$ is a caveat of any attempt to
investigate the dependence of $f_{\textrm{esc}}$ on [O iii] + H$\beta$
strength at $z>3$ from our data set.

Therefore, our analysis has found a weak negative correlation between
$f_{\textrm{esc}}$ and stellar mass, but contrary to expectations of models of
reionisation, we do not see any significant correlation between
$f_{\textrm{esc}}$ and specific SFR of galaxies. We have also found that only
a small fraction ($\approx 6\%$) of galaxies in our sample can be reliably
classified as candidate LyC leakers, although the uncertainties on the
measured $f_{\textrm{esc}}$ are large owing to low significance of the
measured LyC flux. The question that naturally arises following these
findings, and is paramount to understanding the nature of sources in the early
Universe that drove the bulk of reionisation, is: What regulates
$f_{\textrm{esc}}$ in a galaxy? We discuss the relevance and implications of
our findings in the following section.

##  6 Discussion – What regulates $f_{\textrm{esc}}$?

The escape fraction of ionising photons from star-forming galaxies crucially
depends on the fraction of ionising radiation absorbed/scattered by the
neutral gas and dust around the photon emitting regions. To have high
$f_{\textrm{esc}}$, a large amount of neutral gas and dust must have been
expelled from the star-forming regions through either supernova or AGN
feedback (e.g. Wise et al., 2014; Trebitsch et al., 2017, 2018; Barrow et al.,
2020; Katz et al., 2020). Given the dependence of high $f_{\textrm{esc}}$ from
star-forming galaxies on the presence of ionised channels through which LyC
photons can escape into the IGM, the orientation of the galaxies and potential
(mis)alignment of these channels with the observer will play a key role in
whether or not a galaxy appears to be leaking LyC radiation (e.g. Borthakur et
al., 2014).

Therefore, by selection only a fraction of intermediate redshifts reionisation
era analogue galaxies in any given sample will be identified as LyC leakers,
assuming that these galaxies uniformly leak non-zero fractions of LyC photons.
This geometrical misalignment of the viewing angles of a majority of sources
may explain the relatively low numbers of individual galaxies identified in
our sample, and across other samples that appear to be leaking significant
fractions of LyC photons (e.g. Fletcher et al., 2019). A further explanation
to the relatively low $f_{\textrm{esc}}$ values seen in the more relatively
massive star-forming galaxies at intermediate redshifts could be the higher
metallicities compared to the ultra-low metallicities expected from galaxies
in the reionisation era (Yoo et al., 2020).

Another important factor that has an impact on whether or not a galaxy
selected from a large sample appears to be a LyC leaker is the stochasticity
of the IGM along the line-of-sight. Using ensembles of simulated IGM
transmission functions along multiple lines-of-sight at $z\gtrsim 3$, Bassett
et al. (2021) showed that for a given underlying distribution of LyC
$f_{\textrm{esc}}$ in a galaxy sample, only those galaxies that lie along
sightlines with higher IGM transmission appear to be leaking LyC radiation,
exhibiting high $f_{\textrm{esc}}$. Inclusion of a possible IGM transmission
bias for both candidate LyC leakers and non leakers in this work may result in
reduced $f_{\textrm{esc}}$ values for the leakers, with an increase in the
upper limits of $f_{\textrm{esc}}$ inferred from stacking the non-leakers
(e.g. Bassett et al., 2021).

In addition to observational effects such as orientation or stochastic IGM
transmission, theoretical modelling has also revealed complex relationships
between stellar population properties such as metallicity, ionisation
parameter, etc. and $f_{\textrm{esc}}$, which impact the observed
$f_{\textrm{esc}}$ and its dependence on physical properties of stars in a
star-forming galaxy (e.g. Bassett et al., 2019). Further complexities are
introduced by geometrical and viewing angle effects that further impact
correlations between $f_{\textrm{esc}}$ and the strength of the ionising
radiation field and ionising photon production rates from young stars traced
by strong [O iii] + H$\beta$ line strengths (e.g. Katz et al., 2020), among
other galaxy properties. Barrow et al. (2020) showed that the periods during
which a galaxy exhibits intense star-formation activity and ionising photon
production, and periods corresponding to high $f_{\textrm{esc}}$ may be non-
coincident, with high periods of $f_{\textrm{esc}}$ often following starburst
activity (see also Kimm & Cen, 2014; Trebitsch et al., 2017; Kakiichi &
Gronke, 2021).

However, it is important to note that conditions enabling high
$f_{\textrm{esc}}$ following a starburst event may not necessarily coincide
with an equally high production rate or efficiency of ionising photons
($\xi_{\rm{ion}}$), as also seen in the simulations of Barrow et al. (2020).
Therefore, there is an element of fortuitousness introduced in whether or not
a star-forming galaxy selected from a large sample at any given epoch appears
to be leaking LyC photons, and crucially depends on periods of both high
$\xi_{\rm{ion}}$ as well as high $f_{\textrm{esc}}$ being coincident (see also
Fletcher et al., 2019; Nakajima et al., 2020).

Further, a spatially resolved analysis of the lensed LyC leaker dubbed
‘Sunburst’ at $z=2.37$ (Dahle et al., 2016; Rivera-Thorsen et al., 2019)
recently showed that a young massive star cluster dominates the leaking LyC
radiation from the entire galaxy ($f_{\textrm{esc}}^{\textrm{rel}}$ $\approx
0.43-0.93$), whereas older clusters show considerably lower LyC
$f_{\textrm{esc}}$ values (Vanzella et al., 2021). Without magnification due
to lensing, the spatially unresolved Sunburst arc would exhibit a much lower
LyC leakage overall ($f_{\textrm{esc}}^{\textrm{rel}}$ $\approx 0.06-0.12$) as
the LyC radiation escaping along other lines of sight would be missed due to
insufficient LyC imaging depths (see also Vanzella et al., 2020).

Therefore, in a scenario where a large fraction of LyC photons produced by
only a limited number of individual young star clusters within galaxies manage
to escape into the IGM through favourably aligned channels, there would be
increased scatter in the observed dependence of $f_{\textrm{esc}}$ and galaxy
properties such as stellar mass or specific SFR for any given sample of
galaxies. Although this would explain the lack of significant correlations
observed between $f_{\textrm{esc}}$ and the physical properties of candidate
LyC leakers in this study, most importantly with the sSFR, such an effect
complicates the observational understanding of what regulates the escape of
ionising photons and consequently, which sources were the key drivers of
reionisation. This makes studying indicators that trace the presence of LyC
photons escape channels and their dependence on key galaxy properties crucial
to understanding the physical and geometrical conditions that enable a high
apparent $f_{\textrm{esc}}$ in the epoch of reionisation (e.g. Steidel et al.,
2018).

One such probe is the Ly$\alpha$ line, whose shape, width and intensity can
offer insights into LyC and Ly$\alpha$ escape channels from clusters of young
stars within high redshift galaxies (e.g. Vanzella et al., 2018, 2020;
Verhamme et al., 2017; Izotov et al., 2018a; Izotov et al., 2018b; Matthee et
al., 2018; Meyer et al., 2021; Naidu et al., 2021). Additional high ionisation
emission lines like He ii $\lambda 1640$, C iv $\lambda\lambda 1548,1550$ and
C iii] $\lambda\lambda 1907,1909$ in the rest-UV may help trace stellar
populations that produce copious amounts of ionising photons, possibly leading
to higher $f_{\textrm{esc}}$ (e.g. Schaerer et al., 2018; Berg et al., 2019;
Senchyna et al., 2019; Llerena et al., 2021). Saxena et al. (2020b)
investigated He ii emitting star-forming galaxies in the redshift range $z\sim
2.5-5$ and did not find any strong dependence of the presence or the strength
of He ii emission on the stellar mass or star-formation rates of galaxies (see
also Nanayakkara et al., 2019), which may be consistent with factors that lead
to a lack of observed dependence of $f_{\textrm{esc}}$ on these properties as
discussed above.

A key advantage of our galaxy sample at $z\approx 3-3.5$ is the availability
of rest-frame UV spectra from various spectroscopic surveys, crucially also
covering the Ly$\alpha$ emission line, and in a future study we will explore
the correlations between $f_{\textrm{esc}}$, Ly$\alpha$ and other rest-UV
spectroscopic indicators for these newly discovered candidate LyC leakers (and
non-leakers). _JWST_ will give access, for the first time, to the full suite
of rest-frame optical emission lines in $z>3$ galaxies, enabling direct
measurements of the ionising radiation field strengths produced by young stars
(through [O iii] line strengths and [O iii]/[O ii] ratios) as well as
measurements of the ionising photon production efficiencies using H$\alpha$
measurements, which will considerably aid in identifying physical conditions
within star-forming galaxies that lend themselves to high $f_{\textrm{esc}}$.
The correlations between $f_{\textrm{esc}}$ and rest-frame optical emission
established at these intermediate redshifts will be instrumental in
interpreting the ionising photon escape of galaxies at $z=6-9$, for which JWST
will enable rest-optical spectroscopy, in a bid to understand the nature of
the main drivers of reionisation.

Therefore, going forward it is essential to obtain deeper LyC imaging for
intermediate redshift analogues of reionisation era galaxies to properly
account for the viewing angle and orientation effects when measuring
$f_{\textrm{esc}}$. Further, detailed studies of the rest-frame UV and optical
spectra of known LyC leakers are needed to investigate the presence of LyC
escape channels and study the dependence of these features on galaxy
properties. Better insights into the timescales and the incidence of high
$f_{\textrm{esc}}$ within star-forming galaxies may be needed from high-
resolution simulations (e.g. Mauerhofer et al., 2021) of a larger number of
reionisation era galaxies (with sufficient time resolution) to temporally tie
periods of intense star-formation to $f_{\textrm{esc}}$ and other galaxy
properties.

##  7 Summary and conclusions

We have measured the Lyman continuum (LyC) escape fraction
($f_{\textrm{esc}}$) using ground based LyC images available for 183
spectroscopically confirmed galaxies in the CDFS/GOODS-S field in the redshift
range $3.11<z<3.53$. The spectroscopic redshifts are taken from publicly
available databases, that primarily targeted the higher mass Lyman break (LBG)
selected galaxies or blindly discovered lower mass galaxies thanks to emission
line surveys.

We use available deep HST photometry and colours at non-ionising UV
wavelengths to remove possible foreground contaminants that may strongly bias
the inferred LyC flux and $f_{\textrm{esc}}$ as well as accurately constrain
the spectral energy distributions (SEDs) of our galaxies giving robust
measurements of key galaxy properties. Using a combination of stellar and
nebular emission, we also derive rest-frame optical emission line strengths,
primarily the [O iii] and H$\beta$ line fluxes (that fall in the observed $K$
band) for our galaxies.

We then measure LyC flux and $f_{\textrm{esc}}$ using aperture photometry and
detect LyC flux at $\geq 2\sigma$ significance levels in 11 out of 183
galaxies, with $f_{\textrm{esc}}$ values ranging from $0.14-0.85$. Around
$94\%$ of the galaxies in our sample remain undetected at the $\gtrsim\sigma$
level, and using weighted averaging of their LyC $f_{\textrm{esc}}$ limits, we
place a $1\sigma$ upper limit of $f_{\textrm{esc}}$ $<0.12$ from the galaxies
with no significant detections in the available LyC images. We additionally
compare the weighted average $f_{\textrm{esc}}$ limits from sub-samples of
galaxies that were LBG or blind-selected, finding $1\sigma$ limits of
$f_{\textrm{esc}}$ $<0.11$ for LBG-selected and $f_{\textrm{esc}}$ $<0.13$ for
blind-selected galaxies. We further place $f_{\textrm{esc}}$ limits from
normal emission line galaxies (ELGs; [O iii] + H$\beta$ $<300$ Å) and extreme
emission line galaxies (EELGs; [O iii] + H$\beta$ $>300$ Å), finding $1\sigma$
limits of $f_{\textrm{esc}}$ $<0.14$ for normal and $f_{\textrm{esc}}$ $<0.20$
for extreme emission line galaxies.

Although the uncertainties on the measured $f_{\textrm{esc}}$ values for the
11 candidate LyC leaking galaxies are high and a handful of sources have low
signal-to-noise ratios in the LyC images ($\approx 2\sigma)$, we find a weak
negative correlation between $f_{\textrm{esc}}$ and stellar mass at the
$2\sigma$ level, with a shallow slope of $-0.10\pm 0.02$. We do not find any
significant correlation between $f_{\textrm{esc}}$ and specific star-formation
rate of the individual LyC leakers at the $2\sigma$ level, with a slope of
$-0.028\pm 0.015$. Lastly, we find a positive correlation between
$f_{\textrm{esc}}$ and EW0([O iii] + H$\beta$) for the candidate LyC leakers
in our sample, suggesting that the [O iii] + H$\beta$ line strengths may be
important indicators of $f_{\textrm{esc}}$ as has been previously reported in
the literature.

From our results, it remains unclear what the dominant galaxy property is when
it comes to regulating $f_{\textrm{esc}}$ in LyC leaking galaxies. The low
detection fraction of LyC leakers across our sample is consistent with a
scenario where LyC photons escape through discrete ionised channels within
star-forming galaxies, which are seldom oriented along the line-of-sight of
the observer. The time-dependent nature of high $f_{\textrm{esc}}$ that often
follows periods of intense star-formation also suggests that only a low
fraction of galaxies at any given epoch would be actively leaking LyC photons,
consistent with our observations.

The stochasticity of the intervening neutral intergalactic medium, which may
play an important role in absorbing leaking LyC photons along the line-of-
sight is another factor that may impact LyC detections. Future work tying LyC
escape to the presence of ionised channels within galaxies using observations,
exploration of the time-dependent nature of LyC leakage from larger samples of
simulated reionisation era galaxies as well as investigation of IGM
stochasticity across redshifts are needed to fully understand these effects in
the context of contribution of star-forming galaxies towards the cosmic
reionisation photon budget.

In conclusion, our work has discovered a sample of 11 previously undiscovered
candidate LyC leakers at $z\sim 3-3.5$ adding crucial statistics that are
needed to understand the escape of ionising photons from intermediate redshift
analogues of $z>6$ galaxies. Studying lower redshift analogues of reionisation
era galaxies, where direct detection of escaping LyC radiation is possible,
remains the only way to constrain LyC escape from star-forming galaxies that
are expected to drive the bulk of reionisation. An increased sample of LyC
leaking galaxies complemented by deeper LyC/UV imaging, higher resolution
spectra across wavelengths and zoom-in simulations will help understand the
unique physical conditions that drive high $f_{\textrm{esc}}$, bringing us a
step closer towards unravelling the nature of galaxies and their ionising
radiation that drove cosmic reionisation.

## Acknowledgements

We thank the referee for constructive feedback that improved the quality of
this work. AS thanks Ross McLure and Mario Nonino for useful input and
discussions, Fergus Cullen for sharing NIRVANDELS measurements and Mengtao
Tang for sharing the best-fitting relation for emission line strengths from
their analysis. AS and LP thank Alice Shapley for useful comments and
discussions that improved the analysis presented.

AS and RSE acknowledge financial support from European Research Council
Advanced Grant FP7/669253. DS and ASL acknowledge support from Swiss National
Science Foundation. RA acknowledges support from ANID FONDECYT Regular
1202007. MLl acknowledges support from the National Agency for Research and
Development (ANID)/Scholarship Program/Doctorado Nacional/2019-21191036.

This work makes use of VANDELS data. The VANDELS collaboration acknowledges
the invaluable role played by ESO staff for successfully carrying out the
survey. This work is based on VUDS data, obtained with the European Southern
Observatory Very Large Telescope, Paranal, Chile, under Large Program
185.A-0791, and made available by the VUDS team at the CESAM data center,
Laboratoire d’Astrophysique de Marseille, France. This work is based on
observations taken by the MUSE-Wide and MUSE-HUDF Surveys as part of the MUSE
Consortium. This work is based on observations taken by the 3D-HST Treasury
Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS5-26555. MLl acknowledges support from the National Agency for
Research and Development (ANID)/Scholarship Program/Doctorado
Nacional/2019-21191036.

This work has made extensive use of Jupyter and IPython (Pérez & Granger,
2007) notebooks, Astropy (Astropy Collaboration et al., 2013, 2018),
matplotlib (Hunter, 2007), seaborn (Waskom & the seaborn development team,
2020), pandas (pandas development team, 2020) and topcat (Taylor, 2005). This
work would not have been possible without the countless hours put in by
members of the open-source developing community all around the world, with
special thanks to the Stackoverflow community.

## Data availability

The VANDELS data can be accessed using the VANDELS database at
http://vandels.inaf.it/dr4.html, or through the ESO archives. The VUDS data
can be accessed from the VUDS webpage: http://cesam.lam.fr/vuds/ or ESO
archives. MUSE data can be accessed from https://musewide.aip.de/project/ and
http://muse-vlt.eu/science/udf/ and the 3D-HST data can be accessed from
https://archive.stsci.edu/prepds/3d-hst/. The VIMOS10 data are available from
Balestra et al. (2010). All the code used to perform the analysis in this
paper was written in python and will be shared upon reasonable request to the
corresponding author.

## References

  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Balestra et al. (2010) Balestra I., et al., 2010, A&A, 512, A12
  * Barrow et al. (2020) Barrow K. S. S., Robertson B. E., Ellis R. S., Nakajima K., Saxena A., Stark D. P., Tang M., 2020, ApJ, 902, L39
  * Bassett et al. (2019) Bassett R., et al., 2019, MNRAS, 483, 5223
  * Bassett et al. (2021) Bassett R., Ryan-Weber E. V., Cooke J., Meštrić U., Kakiichi K., Prichard L., Rafelski M., 2021, MNRAS, 502, 108
  * Berg et al. (2019) Berg D. A., Chisholm J., Erb D. K., Pogge R., Henry A., Olivier G. M., 2019, ApJ, 878, L3
  * Bian & Fan (2020) Bian F., Fan X., 2020, MNRAS, 493, L65
  * Bian et al. (2017) Bian F., Fan X., McGreer I., Cai Z., Jiang L., 2017, ApJ, 837, L12
  * Borthakur et al. (2014) Borthakur S., Heckman T. M., Leitherer C., Overzier R. A., 2014, Science, 346, 216
  * Bradley et al. (2020) Bradley L., et al., 2020, astropy/photutils: 1.0.0, doi:10.5281/zenodo.4044744, https://doi.org/10.5281/zenodo.4044744
  * Brammer et al. (2012) Brammer G. B., et al., 2012, ApJS, 200, 13
  * Bromm & Yoshida (2011) Bromm V., Yoshida N., 2011, ARA&A, 49, 373
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Castellano et al. (2017) Castellano M., et al., 2017, ApJ, 839, 73
  * Cristiani et al. (2016) Cristiani S., Serrano L. M., Fontanot F., Vanzella E., Monaco P., 2016, MNRAS, 462, 2478
  * Cullen et al. (2021) Cullen F., et al., 2021, arXiv e-prints, p. arXiv:2103.06300
  * Dahle et al. (2016) Dahle H., et al., 2016, A&A, 590, L4
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1
  * Dayal et al. (2020) Dayal P., et al., 2020, MNRAS, 495, 3065
  * Dickinson et al. (2003) Dickinson M., Giavalisco M., GOODS Team 2003, in Bender R., Renzini A., eds, The Mass of Galaxies at Low and High Redshift. p. 324 (arXiv:astro-ph/0204213), doi:10.1007/10899892_78
  * Endsley et al. (2020) Endsley R., Stark D. P., Chevallard J., Charlot S., 2020, arXiv e-prints, p. arXiv:2005.02402
  * Faisst (2016) Faisst A. L., 2016, ApJ, 829, 99
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Finkelstein et al. (2019) Finkelstein S. L., et al., 2019, ApJ, 879, 36
  * Fletcher et al. (2019) Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87
  * Fontanot et al. (2012) Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413
  * Forrest et al. (2017) Forrest B., et al., 2017, ApJ, 838, L12
  * Garilli et al. (2012) Garilli B., Paioro L., Scodeggio M., Franzetti P., Fumana M., Guzzo L., 2012, PASP, 124, 1232
  * Garilli et al. (2021) Garilli B., et al., 2021, A&A, 647, A150
  * Giallongo et al. (2015) Giallongo E., et al., 2015, A&A, 578, A83
  * Grazian et al. (2016) Grazian A., et al., 2016, A&A, 585, A48
  * Grazian et al. (2017) Grazian A., et al., 2017, A&A, 602, A18
  * Grazian et al. (2018) Grazian A., et al., 2018, A&A, 613, A44
  * Grogin et al. (2011) Grogin N. A., et al., 2011, ApJS, 197, 35
  * Guaita et al. (2016) Guaita L., et al., 2016, A&A, 587, A133
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Inami et al. (2017) Inami H., et al., 2017, A&A, 608, A2
  * Inoue et al. (2014) Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805
  * Izotov et al. (2017) Izotov Y. I., Thuan T. X., Guseva N. G., 2017, MNRAS, 471, 548
  * Izotov et al. (2018a) Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., Fricke K. J., 2018a, MNRAS, 474, 4514
  * Izotov et al. (2018b) Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, MNRAS, 478, 4851
  * Izotov et al. (2021) Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Chisholm J., Thuan T. X., Fricke K. J., Verhamme A., 2021, MNRAS, 503, 1734
  * Japelj et al. (2017) Japelj J., et al., 2017, MNRAS, 468, 389
  * Jaskot & Oey (2013) Jaskot A. E., Oey M. S., 2013, ApJ, 766, 91
  * Jaskot et al. (2019) Jaskot A. E., Dowd T., Oey M. S., Scarlata C., McKinney J., 2019, ApJ, 885, 96
  * Ji et al. (2020) Ji Z., et al., 2020, ApJ, 888, 109
  * Kakiichi & Gronke (2021) Kakiichi K., Gronke M., 2021, ApJ, 908, 30
  * Katz et al. (2020) Katz H., et al., 2020, MNRAS, 498, 164
  * Kimm & Cen (2014) Kimm T., Cen R., 2014, ApJ, 788, 121
  * Koekemoer et al. (2011) Koekemoer A. M., et al., 2011, ApJS, 197, 36
  * Kriek et al. (2015) Kriek M., et al., 2015, ApJS, 218, 15
  * Kroupa & Boily (2002) Kroupa P., Boily C. M., 2002, MNRAS, 336, 1188
  * Labbé et al. (2013) Labbé I., et al., 2013, ApJ, 777, L19
  * Le Fèvre et al. (2015) Le Fèvre O., et al., 2015, A&A, 576, A79
  * Llerena et al. (2021) Llerena M., et al., 2021, arXiv e-prints, p. arXiv:2107.00660
  * Luo et al. (2017) Luo B., et al., 2017, ApJS, 228, 2
  * Ma et al. (2020) Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001
  * Madau & Fragos (2017) Madau P., Fragos T., 2017, ApJ, 840, 39
  * Magliocchetti et al. (2020) Magliocchetti M., et al., 2020, MNRAS, 493, 3838
  * Marchi et al. (2017a) Marchi F., et al., 2017a, A&A, 601, A73
  * Marchi et al. (2017b) Marchi F., et al., 2017b, A&A, 601, A73
  * Marques-Chaves et al. (2021) Marques-Chaves R., Schaerer D., Álvarez-Márquez J., Colina L., Dessauges-Zavadsky M., Pérez-Fournon I., Saldana-Lopez A., Verhamme A., 2021, MNRAS, 507, 524
  * Matthee et al. (2018) Matthee J., Sobral D., Gronke M., Paulino-Afonso A., Stefanon M., Röttgering H., 2018, A&A, 619, A136
  * Mauerhofer et al. (2021) Mauerhofer V., Verhamme A., Blaizot J., Garel T., Kimm T., Michel-Dansac L., Rosdahl J., 2021, A&A, 646, A80
  * McLure et al. (2018) McLure R. J., et al., 2018, MNRAS, 479, 25
  * Meštrić et al. (2020) Meštrić U., et al., 2020, MNRAS, 494, 4986
  * Meštrić et al. (2021) Meštrić U., Ryan-Weber E. V., Cooke J., Bassett R., Prichard L. J., Rafelski M., 2021, MNRAS, 508, 4443
  * Meyer et al. (2021) Meyer R. A., Laporte N., Ellis R. S., Verhamme A., Garel T., 2021, MNRAS, 500, 558
  * Momcheva et al. (2016) Momcheva I. G., et al., 2016, ApJS, 225, 27
  * Mostardi et al. (2013) Mostardi R. E., Shapley A. E., Nestor D. B., Steidel C. C., Reddy N. A., Trainor R. F., 2013, ApJ, 779, 65
  * Mostardi et al. (2015) Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A., Siana B., 2015, ApJ, 810, 107
  * Naidu et al. (2018) Naidu R. P., Forrest B., Oesch P. A., Tran K.-V. H., Holden B. P., 2018, MNRAS, 478, 791
  * Naidu et al. (2020) Naidu R. P., Tacchella S., Mason C. A., Bose S., Oesch P. A., Conroy C., 2020, ApJ, 892, 109
  * Naidu et al. (2021) Naidu R. P., et al., 2021, arXiv e-prints, p. arXiv:2110.11961
  * Nakajima & Ouchi (2014) Nakajima K., Ouchi M., 2014, MNRAS, 442, 900
  * Nakajima et al. (2018) Nakajima K., Fletcher T., Ellis R. S., Robertson B. E., Iwata I., 2018, MNRAS, 477, 2098
  * Nakajima et al. (2020) Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ, 889, 161
  * Nanayakkara et al. (2019) Nanayakkara T., et al., 2019, A&A, 624, A89
  * Nonino et al. (2009) Nonino M., et al., 2009, ApJS, 183, 244
  * Oke & Gunn (1983) Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
  * Paalvast et al. (2018) Paalvast M., et al., 2018, A&A, 618, A40
  * Pahl et al. (2021) Pahl A. J., Shapley A., Steidel C. C., Chen Y., Reddy N. A., 2021, MNRAS,
  * Pentericci et al. (2018) Pentericci L., et al., 2018, A&A, 616, A174
  * Pérez & Granger (2007) Pérez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21
  * Planck Collaboration et al. (2016a) Planck Collaboration et al., 2016a, A&A, 594, A13
  * Planck Collaboration et al. (2016b) Planck Collaboration et al., 2016b, A&A, 596, A108
  * Reddy et al. (2018) Reddy N. A., et al., 2018, ApJ, 869, 92
  * Rivera-Thorsen et al. (2019) Rivera-Thorsen T. E., et al., 2019, Science, 366, 738
  * Roberts-Borsani et al. (2016) Roberts-Borsani G. W., et al., 2016, ApJ, 823, 143
  * Robertson et al. (2013) Robertson B. E., et al., 2013, ApJ, 768, 71
  * Robertson et al. (2015) Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19
  * Romano et al. (2019) Romano M., Grazian A., Giallongo E., Cristiani S., Fontanot F., Boutsia K., Fiore F., Menci N., 2019, A&A, 632, A45
  * Rutkowski et al. (2017) Rutkowski M. J., et al., 2017, ApJ, 841, L27
  * Santini et al. (2017) Santini P., et al., 2017, ApJ, 847, 76
  * Saxena et al. (2020a) Saxena A., et al., 2020a, MNRAS, 496, 3796
  * Saxena et al. (2020b) Saxena A., et al., 2020b, A&A, 636, A47
  * Saxena et al. (2021) Saxena A., et al., 2021, MNRAS, 505, 4798
  * Schaerer et al. (2018) Schaerer D., Izotov Y. I., Nakajima K., Worseck G., Chisholm J., Verhamme A., Thuan T. X., de Barros S., 2018, A&A, 616, L14
  * Schaerer et al. (2019) Schaerer D., Fragos T., Izotov Y. I., 2019, A&A, 622, L10
  * Senchyna et al. (2019) Senchyna P., Stark D. P., Chevallard J., Charlot S., Jones T., Vidal-García A., 2019, MNRAS, 488, 3492
  * Shapley et al. (2016) Shapley A. E., Steidel C. C., Strom A. L., Bogosavljević M., Reddy N. A., Siana B., Mostardi R. E., Rudie G. C., 2016, ApJ, 826, L24
  * Siana et al. (2007) Siana B., et al., 2007, ApJ, 668, 62
  * Siana et al. (2015) Siana B., et al., 2015, ApJ, 804, 17
  * Smit et al. (2015) Smit R., et al., 2015, ApJ, 801, 122
  * Smith et al. (2020) Smith B. M., et al., 2020, ApJ, 897, 41
  * Stark (2016) Stark D. P., 2016, ARA&A, 54, 761
  * Stark et al. (2013) Stark D. P., Schenker M. A., Ellis R., Robertson B., McLure R., Dunlop J., 2013, ApJ, 763, 129
  * Stasińska et al. (2015) Stasińska G., Izotov Y., Morisset C., Guseva N., 2015, A&A, 576, A83
  * Steidel et al. (2001) Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665
  * Steidel et al. (2018) Steidel C. C., Bogosavljević M., Shapley A. E., Reddy N. A., Rudie G. C., Pettini M., Trainor R. F., Strom A. L., 2018, ApJ, 869, 123
  * Tang et al. (2019) Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572
  * Tasca et al. (2017) Tasca L. A. M., et al., 2017, A&A, 600, A110
  * Taylor (2005) Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R., eds, Astronomical Society of the Pacific Conference Series Vol. 347, Astronomical Data Analysis Software and Systems XIV. p. 29
  * Trebitsch et al. (2017) Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224
  * Trebitsch et al. (2018) Trebitsch M., Volonteri M., Dubois Y., Madau P., 2018, MNRAS, 478, 5607
  * Urrutia et al. (2019) Urrutia T., et al., 2019, A&A, 624, A141
  * Vanzella et al. (2010) Vanzella E., Siana B., Cristiani S., Nonino M., 2010, MNRAS, 404, 1672
  * Vanzella et al. (2016a) Vanzella E., et al., 2016a, ApJ, 821, L27
  * Vanzella et al. (2016b) Vanzella E., et al., 2016b, ApJ, 825, 41
  * Vanzella et al. (2018) Vanzella E., et al., 2018, MNRAS, 476, L15
  * Vanzella et al. (2020) Vanzella E., et al., 2020, MNRAS, 491, 1093
  * Vanzella et al. (2021) Vanzella E., et al., 2021, arXiv e-prints, p. arXiv:2106.10280
  * Verhamme et al. (2017) Verhamme A., Orlitová I., Schaerer D., Izotov Y., Worseck G., Thuan T. X., Guseva N., 2017, A&A, 597, A13
  * Waskom & the seaborn development team (2020) Waskom M., the seaborn development team 2020, mwaskom/seaborn, doi:10.5281/zenodo.592845, https://doi.org/10.5281/zenodo.592845
  * Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  * Yoo et al. (2020) Yoo T., Kimm T., Rosdahl J., 2020, MNRAS, 499, 5175
  * Zackrisson et al. (2013) Zackrisson E., Inoue A. K., Jensen H., 2013, ApJ, 777, 39
  * de Barros et al. (2016) de Barros S., et al., 2016, A&A, 585, A51
  * pandas development team (2020) pandas development team T., 2020, pandas-dev/pandas: Pandas, doi:10.5281/zenodo.3509134, https://doi.org/10.5281/zenodo.3509134

##  Appendix A LyC and HST Cutout galleries with best-fitting SEDs of
candidate LyC leakers

Figure 12: LyC and HST photometry (top), and best-fitting SED (orange) with
observed (blue points) and SED derived (orange ovals) for candidate LyC
leaking galaxies. We show the observed and SED predicted flux at LyC
wavelengths only for demonstration purposes – the SED derived LyC fluxes were
not used to infer $f_{\textrm{esc}}$ and the observed LyC fluxes were not used
during SED fitting to obtain galaxy physical properties.

Figure 12: Continued.

Figure 12: Continued.