11institutetext: Kapteyn Astronomical Institute, University of Groningen, P.O.
Box 800, 9700 AV Groningen, The Netherlands 22institutetext: Department of
Astronomy, Université de Genève, 51 Ch. Pegasi, 1290 Versoix, Switzerland
33institutetext: Cosmic Dawn Center (DAWN), Niels Bohr Institute, University
of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark 44institutetext:
Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482
Potsdam, Germany 55institutetext: d’Astrophysique de Paris, CNRS, Sorbonne
Université, 98bis Boulevard Arago, 75014, Paris, France 66institutetext:
Inter-University Centre for Astronomy and Astrophysics, Pune 411 007, India
77institutetext: Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95064, USA 88institutetext: Department of Physics,
ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich, 8093, Switzerland
99institutetext: Instituto de Astrofísica de Canarias, c/ Vía Láctea s/n,
E-38205 - La Laguna, Tenerife, Spain
Departamento de Astrofísica, Universidad de La Laguna, E-38205 - La Laguna,
Tenerife, Spain 1010institutetext: MIT Kavli Institute for Astrophysics and
Space Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA
1111institutetext: Department for Astrophysical and Planetary Science,
University of Colorado, Boulder, CO 80309, USA 1212institutetext: Department
of Physics and Astronomy, University of California, Riverside, 900 University
Avenue, Riverside, CA 92521, USA 1313institutetext: Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK.
1414institutetext: Cavendish Laboratory, University of Cambridge, 19 JJ
Thomson Avenue, Cambridge, CB3 0HE, UK 1515institutetext: European Southern
Observatory, Av. Alonso de Córdova 3107, 763 0355 Vitacura, Santiago, Chile

# Lyman Continuum Leaker Candidates at $z\sim 3-4$ in the HDUV Based on a
Spectroscopic Sample of MUSE LAEs

J. Kerutt<EMAIL_ADDRESS>1122 P. A. Oesch , 2233 L. Wisotzki 44 A.
Verhamme 22 H. Atek 55 E. C. Herenz 66 G. D. Illingworth 77 H. Kusakabe 22 J.
Matthee 88 V. Mauerhofer 11 M. Montes 99 R. P. Naidu 1010 E. Nelson 1111 N.
Reddy 1212 J. Schaye 66 C. Simmonds , 13131414 T. Urrutia 44 E. Vitte , 221515

(Received ; accepted )

###### Abstract

Context. In recent years, a number of Lyman continuum (LyC) leaker candidates
at intermediate redshifts have been found, providing insight into how the
Universe was reionised at early cosmic times. Around a hundred LyC leakers at
all redshifts are known by now, which enables us to analyse their properties
statistically.

Aims. Here we identify new LyC leaker candidates at $z\approx 3-4.5$ and
compare them to objects from the literature to get an overview of the
different observed escape fractions and their relation to the properties of
the Lyman $\alpha$ (Ly$\alpha$) emission line. The aim of this work is to test
indicators (or proxies) for LyC leakage suggested in the literature and to
improve our understanding of the kind of galaxies from which LyC radiation can
escape.

Methods. We use data from the Hubble Deep Ultraviolet (HDUV) legacy survey to
search for LyC emission based on a sample of $\approx 2000$ Ly$\alpha$
emitters (LAEs) detected previously in two surveys with the Multi-Unit
Spectroscopic Explorer (MUSE), MUSE-Deep and MUSE-Wide. Based on their known
redshifts and positions, we look for potential LyC leakage in the WFC3/UVIS
F336W band of the HDUV. The escape fractions are measured and compared in
different ways, including based on spectral energy distribution (SED) fitting
performed using the CIGALE software.

Results. We add twelve objects to the sample of known LyC leaker candidates
(five highly likely leakers and seven potential ones), one of which was
previously known, and compare their Ly$\alpha$ properties to their escape
fractions. We find escape fractions between $\sim 20\%$ and $\sim 90\%$,
assuming a high transmission in the intergalactic medium (IGM). We show a
method to use the number of LyC leaker candidates we find to infer the
underlying average escape fraction of galaxies, which is $\approx 12\%$.

Conclusions. Based on their Ly$\alpha$ properties, we conclude that LyC
leakers are not very different from other high-z LAEs and suggest that most
LAEs could be leaking LyC even if this can not always be detected due to the
direction of emission and the transmission properties of the IGM.

###### Key Words.:

galaxies: high redshift – galaxies: formation – galaxies: evolution –
cosmology: observations

##  1 Introduction

The Epoch of Reionisation (EoR) is the last phase transition of the universe,
where the intergalactic medium (IGM) went from a neutral to a mostly ionised
state (see review Wise, 2019 and e.g. Fan et al., 2006; Bouwens et al., 2007;
Ouchi et al., 2009a, b; Bouwens et al., 2015; Haiman, 2016). This mostly
coincides with the formation of the first galaxies. While recent years have
brought new insights, there are still aspects we need to understand about this
crucial phase in the development of the universe.

What we can constrain relatively well is the time of the EoR ($z\approx
6\\--8$, Fan et al., 2006), for example using the drop in the fraction of
Ly$\alpha$ emitters (LAEs, e.g. Ouchi et al., 2008; Kashikawa et al., 2011;
Caruana et al., 2012; Kusakabe et al., 2020), which is also seen in
simulations (e.g. Garel et al., 2021). This could either be due to the
intrinsic evolution of the LAE population or an effect of the IGM. While Lyman
continuum (LyC) emission is absorbed by the neutral hydrogen in the IGM at the
EoR, Ly$\alpha$ (Ly$\alpha$) is scattered out of the line of sight and is,
therefore, an indicator of the neutral fraction of the IGM.

In recent years, observational results pointed toward star-forming galaxies
such as LAEs being the best candidates for providing the ionising emission at
and after the EoR (Fontanot et al., 2014; Onoue et al., 2017; Matsuoka et al.,
2018; Japelj et al., 2017; Naidu et al., 2020; Matthee et al., 2022; Naidu et
al., 2022; Begley et al., 2022; Matsuoka et al., 2022), with only a minimal
contribution by active galactic nuclei (AGN, Cowie et al., 2009; Madau &
Haardt, 2015; Smith et al., 2020; Trebitsch et al., 2021). Since the number
density of AGN decreases rapidly at $z>3$ (e.g. Masters et al., 2012) and the
escape fraction of ionising photons from AGN is not as high as needed (e.g.
Micheva et al., 2017), the only possibility for AGN to contribute
significantly to the EoR would be through a large number of low luminosity
AGN, as has been claimed for example by Giallongo et al. (2015) and Grazian et
al. (2018). Recent James Webb Space Telescope (JWST) observations (Kocevski et
al., 2023; Übler et al., 2023; Harikane et al., 2023) seem to confirm this
possibility, while Parsa et al. (2018) found significantly fewer faint AGN at
redshifts $z>4$ (see also recent results using Subaru/Suprime-Cam and Hyper
Suprime-Cam in Onoue et al., 2017; Matsuoka et al., 2018, 2022 at redshifts
$z\approx 6$). It is also possible, that AGN contribute significant amounts of
LyC emission starting at $z\approx 2-3$, after the EoR (e.g. Becker et al.,
2015; Madau & Haardt, 2015; Faucher-Giguère, 2020; Trebitsch et al., 2021),
which means the question of the AGN contribution to the (re)ionisation of the
universe is not fully answered yet.

Having thus mostly ruled out AGN as the sources of reionisation, the focus has
shifted towards star-forming galaxies. An escape fraction of ionising photons
of $\approx 10\%$ is needed to explain the EoR (Pawlik et al., 2009; Vanzella
et al., 2012; Mitra et al., 2015; Giallongo et al., 2015; Robertson et al.,
2015; Madau & Haardt, 2015; Feng et al., 2016 but also Faucher-Giguère et al.,
2008; Matthee et al., 2017, who find lower necessary escape fractions under
certain assumptions).

Since the neutral fraction of hydrogen in the IGM is rising towards the EoR
(Madau, 1995; Inoue et al., 2014), it is not possible to directly observe the
LyC radiation responsible for the ionisation of the IGM. However, assuming the
properties of the sources of the EoR do not evolve much with redshift, we can
study them at lower redshifts, where the neutral fraction of the IGM allows
the LyC emission to get through (Hu et al., 1998; Steidel et al., 2001;
Shapley et al., 2006; Ouchi et al., 2008; Blanc et al., 2011; Vanzella et al.,
2012).

There have been several discoveries of LyC leaker candidates at $z\approx
3-4$, for example the well-studied Ion 1 Vanzella et al. (2010b); Ji et al.
(2020), Ion 2 Vanzella et al. (2015, 2020) and Ion 3 Vanzella et al. (2018)
(see table 5 in the appendix). A few LyC leakers at those redshifts have very
high escape fractions (e.g. $f_{\mathrm{esc}}^{\mathrm{LyC}}>50\%$ in Vanzella
et al., 2015; de Barros et al., 2016, $f_{\mathrm{esc}}^{\mathrm{LyC}}=52\%$
in Saxena et al., 2022, $f_{\mathrm{esc}}^{\mathrm{LyC}}=90\%$ in Marques-
Chaves et al., 2022, $f_{\mathrm{esc}}^{\mathrm{LyC}}=100\%$ in Rivera-Thorsen
et al., 2022). These high escape fractions could be the result of orientation
effects, as we might see them at the right angle where their LyC escapes,
since in general, escape fractions from individual objects (and stacks) tend
to be low. For example Grimes et al. (2009) find no LyC from a sample of 16
local starburst galaxies, Rutkowski et al. (2016) find non-detections with
upper limits of $f_{\mathrm{esc}}^{\mathrm{LyC}}<2.1\%$ for objects at
$z\approx 1$, and Grazian et al. (2016) only find two LyC leakers in 37
galaxies at $z\approx 3.3$. There also seems to be a trend with redshift
(Mitra et al., 2013; Fontanot et al., 2014; Faisst, 2016; Khaire et al., 2016;
Japelj et al., 2017). However, even at $z\approx 4$, the possible absorption
by the IGM has to be taken into account when measuring escape fractions (e.g.
Bassett et al., 2021), which is not easy considering the stochastic nature of
the different lines of sight (Madau, 1995; Inoue et al., 2014) and makes
measurements of LyC escape fractions rather uncertain.

At slightly lower redshifts, $z\approx 2-3$, where the IGM transmission is
reasonably high, there have been many successful detections of LyC leakers
(e.g. Vanzella et al., 2010a; Mostardi et al., 2015; Grazian et al., 2016;
Shapley et al., 2016; Steidel et al., 2018; Fletcher et al., 2019; Rivera-
Thorsen et al., 2019; Saha et al., 2020). However, it is important to exclude
lower-redshift interlopers, as they can be responsible for several assumed LyC
leaker candidates (Shapley et al., 2006; Vanzella et al., 2010b; Siana et al.,
2015; Mostardi et al., 2015), which is why reliable techniques to exclude such
interlopers are important (see e.g. Pahl et al., 2021, who suggest a colour
selection technique on resolved photometry).

Going to even lower redshifts can be a solution to avoid interlopers, and many
low-redshift analogues of high-redshift star-forming galaxies (in particular
LAEs and Lyman Break Galaxies, LBGs) have been studied in recent years.
Examples of such analogues are Green Peas (Cardamone et al., 2009), which
sometimes exhibit a high relative escape fraction of ionising emission, as
high as $\mathrm{f}_{\mathrm{esc}}^{\mathrm{LyC}}=73\%$ (Izotov et al., 2016a,
Izotov et al., 2016b, Izotov et al., 2018a, Izotov et al., 2018b), while
Izotov et al. (2021) find 15 LyC leakers among a sample of 20 Green Peas.
Other low-z galaxies that are leaking LyC emission are for example Haro 11
(Bergvall et al., 2006; Hayes et al., 2007; Leitet et al., 2011, 2013; Keenan
et al., 2017; Rivera-Thorsen et al., 2017b) with
$f_{\mathrm{esc}}^{\mathrm{LyC,rel}}=16.6^{+7.4}_{-6.5}\%$ (Leitet et al.,
2013), Tol 1247 with $f_{\mathrm{esc}}^{\mathrm{LyC}}=2.4\%$ (Leitet et al.,
2013), Mrk54 with $f_{\mathrm{esc}}^{\mathrm{LyC}}=6.2\%$ (Deharveng et al.,
2001), and J0921 with $f_{\mathrm{esc}}^{\mathrm{LyC}}\sim 1\%$ (Borthakur et
al., 2014, but see table 7 in the appendix). Recently, the number of known LyC
emitters at low-redshifts has been significantly increased with the Low-
Redshift Lyman Continuum Survey (LzLCs, Wang et al., 2021; Flury et al.,
2022a, b; Saldana-Lopez et al., 2022), which includes 35 LyC emitters at
$z=0.2-0.4$ with escape fractions of up to $50\%$ and notably 12 objects with
$f_{\mathrm{esc}}^{\mathrm{LyC}}>5\%$. A caveat to many of the studies of low-
redshift LyC leakers is that they are often selected to have specific
properties that make LyC leakage likely. Therefore, it is not certain whether
they are representative of the general galaxy population.

Individual detections of LyC leakers at high redshifts could introduce a
selection bias in the interpretation of the escape fraction since they could
be extreme and rare objects. However, since the IGM is stochastic and highly
variable along different lines of sight, it is also possible that they have
particularly low foreground opacities and are not different from the general
population. Several studies, therefore, used stacking to obtain the mean
escape fraction, with the first direct detection of LyC emission from star-
forming galaxies at high redshift from stacked spectra of 29 LBGs at $z=3.4$
(Steidel et al., 2001) and an escape fraction of $f_{\mathrm{esc}}\gtrsim
50\%$. Other studies followed, usually finding much lower escape fractions
(e.g. Siana et al., 2010; Grazian et al., 2017; Fletcher et al., 2019; Steidel
et al., 2018), often using LBGs and LAEs to find LyC emission (e.g. Grazian et
al., 2016; Japelj et al., 2017; Grazian et al., 2017; Naidu et al., 2018;
Marchi et al., 2018 and Rutkowski et al., 2016, 2017 at lower redshifts).
Results from these stacking analyses are not conclusive yet, as for example
Steidel et al. (2018) find $f_{\mathrm{esc}}^{\mathrm{LyC}}=9\%$ using LBGs
without prior selection for LyC leakage, barely enough to explain
reionisation. However, they find that fainter galaxies could have higher
escape fractions, which highlights the importance of the sample selection. In
a sample targeting preferentially [O III] emitters (which are thought to be
good candidates for LyC leakage), Fletcher et al. (2019) use 61 LAEs and find
$20\%$ individual LyC leakers, but when excluding the individually-detected
objects from the stack, they find an upper limit in the LyC escape fraction of
$0.3\%$. This could mean that only some objects have a high escape fraction,
but most have close to 0, or that most IGM transmission lines are not
transparent for LyC. Indeed, in a similar study, Bian & Fan (2020) do not find
any individual LyC leaker candidates among a sample of 54 LAEs at $z\approx
3.1$ in the GOODS-S field and no detection in a stack of those objects either.
In contrast, Begley et al. (2022) find an escape fraction of
$f_{\mathrm{esc}}=0.07\pm 0.02$ in a stack of 148 star-forming galaxies at
$z\simeq 3.5$. Individual detections are of course dependent on the depth of
the data and thus somewhat arbitrary unless the same data and detection
criteria are compared. Another aspect is that from simulations we expect the
escape fraction of ionising photons to vary over time (by more than six orders
of magnitude, Trebitsch et al., 2017) and between different sight-lines and
especially with inclination (e.g. Trebitsch et al., 2017; Smith et al., 2022).

The varying results on escape fractions from the literature can also in part
be attributed to the different methods and assumptions made to derive the
escape fractions. Not only is the absorption in the IGM uncertain for
individual objects, but also the stellar population properties such as
metallicity, age, star formation history, and initial mass function as well as
the dust absorption. In addition, only some studies take the effects of the
circumgalactic medium (CGM) into account (e.g. Reddy et al., 2016b; Steidel et
al., 2018; Pahl et al., 2021), while most group them with those of the IGM,
although significant ${H\textsc{I}}$ absorption has been found up to $\approx
500\,\mathrm{km/s}$ for LAEs at $z=2.9-3.8$ (Muzahid et al., 2021).

It would, therefore, be ideal to infer the LyC emission from other
observables, which are less affected by the IGM and model assumptions. Several
such proxies for LyC have been suggested (e.g. Dijkstra, 2014; Verhamme et
al., 2015, 2017; Marchi et al., 2018; Izotov et al., 2018b), including a
connection with Ly$\alpha$ emission, e.g. through the equivalent widths
(Micheva et al., 2017; Marchi et al., 2018; Steidel et al., 2018; Fletcher et
al., 2019; Pahl et al., 2021; Reddy et al., 2022; Flury et al., 2022b) or
double peaks (Verhamme et al., 2017; Vanzella et al., 2020; Izotov et al.,
2021). Contrary to such intuitive trends, one of the highest redshift LyC
leakers, Ion1 at $z=3.8$, has Ly$\alpha$ in absorption (Ji et al., 2020). This
can happen for moderate neutral hydrogen column densities (between
log[N${}_{{H\textsc{I}}}\,\mathrm{cm}^{-2}$] $\approx 13-17$), where the gas
would be optically thin for LyC, but Ly$\alpha$ is already scattered.
Similarly, the strong LyC leaker from Marques-Chaves et al. (2022) has a
surprisingly wide peak separation of $\Delta v(\mathrm{Ly}\alpha)=680\pm
70\,\mathrm{km}\,\mathrm{s}^{-1}$, indicating a lot of scattering for
Ly$\alpha$ and thus a rather large neutral hydrogen column density.

Nevertheless, some degree of connection between Ly$\alpha$ and LyC emission is
expected through their linked production mechanisms and similar liability to
interact with neutral hydrogen. A prerequisite for Ly$\alpha$ emission is the
production of LyC, which ionises neutral hydrogen, and then recombines to emit
a Ly$\alpha$ photon (in $\approx$ 2/3 of the cases for a temperature of
$10^{4}\,\mathrm{K}$, e.g. Dijkstra, 2014). Simplified, more LyC photons could
therefore mean more Ly$\alpha$ photons. While Ly$\alpha$ then scatters in
neutral hydrogen, LyC will be absorbed, which means that if there is LyC
emission, they both benefit from a relatively free path (meaning a low neutral
hydrogen column density) through the interstellar medium (ISM) to easily
escape the galaxy. Not only the creation of Ly$\alpha$ thus depends on the
presence of LyC, but due to its resonant nature, the shape of the Ly$\alpha$
line traces the neutral hydrogen column density in the ISM and CGM. A high
neutral hydrogen column density will prevent LyC from escaping and imprint
itself on the Ly$\alpha$ line shape through the FWHM and peak separation.
Therefore, both properties have been proposed as tracers of LyC emission.
Thus, LAEs could be ideal candidates to look for LyC emission (see e.g.
Rivera-Thorsen et al., 2017b; Vanzella et al., 2018; Gazagnes et al., 2020;
Izotov et al., 2021; Matthee et al., 2022), which is what we test in this
study.

Based on a sample of LAEs (see Kerutt et al., 2022) found with the Multi-Unit
Spectroscopic Explorer (MUSE, Bacon et al., 2010), we look for individual LyC
emitters in the Hubble Deep Ultraviolet (UV) Legacy Survey (HDUV, Oesch et
al., 2018) in the GOODS-South and North (Giavalisco et al., 2004). Naidu et
al. (2017) already found six candidates in the HDUV (three of which might
contain an AGN) using the two HST filter bands WFC3/UVIS F336W and F275W at
$z\approx 2$, but since the Lyman break is in the middle of F275W, computing
the escape fraction is challenging. Jones et al. (2018) use F275W in GOODS-
North for objects at $z\approx 2.4$, which guarantees that only LyC will fall
in the filter band, and find six candidates, four of which are interlopers. A
more recent study with the same filter band WFC3/UVIS F275W focused on the
GOODS-South and resulted in five LyC leaker candidates, two of which might be
interlopers (Jones et al., 2021). This demonstrates the need for carefully
analysing the potential leakers to estimate the contamination, which is what
we strive for in this paper. This allows us to better understand if the
proposed connection between the LyC and Ly$\alpha$ holds up at higher
redshifts. If not, this might be an indication that the Ly$\alpha$ and LyC
emission is not necessarily coming from the same region in the galaxy.

The paper is structured as follows: In Sect. 2 we discuss the data we use for
this study, from the HDUV survey and the MUSE-Wide and -Deep surveys. In Sect.
3 we search for LyC emission from our selected sample of LAEs and discuss
potential contamination. Having defined a sample of LyC leaker candidates, we
then measure their spectral energy distributions (SEDs) and derive LyC escape
fractions, which we discuss in Sect. 4. We compare the escape fractions to
Ly$\alpha$ properties in Sect. 5 in order to find correlations. We discuss our
results in Sect. 6 and give a summary and conclusion in Sect. 7. Throughout
this paper, we use AB magnitudes and physical distances and assume flat
$\Lambda\mathrm{CDM}$ cosmology with
$\mathrm{H}_{0}=70\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$,
$\Omega_{\mathrm{m}}=0.3$ and $\Omega_{\Lambda}=0.7$.

##  2 Data

We base our search for LyC emission on known LAEs from the MUSE-Wide survey
(Urrutia et al., 2019; Kerutt et al., 2022), using their positions and
redshifts as selection criteria to look for possible LyC leakage in the HST
filter bands WFC3/UVIS F275W and F336W, which we take from the HDUV survey
(Oesch et al., 2015, 2018). Both surveys are located in the Chandra Deep Field
South (CDFS, Giacconi et al., 2001) region, which is complemented by the Great
Observatories Origins Deep Survey (GOODS; Giavalisco et al., 2004) consisting
of multiple deep observations of the Hubble Space Telescope (HST) and Spitzer.
In Fig. 1 we show the footprints of the HDUV, the MUSE-Deep and -Wide surveys,
overlaid on an inverted HST image. As can be seen, the HDUV survey and the
MUSE-Wide and -Deep surveys overlap significantly, with most of the HDUV being
covered by MUSE observations.

Figure 1: Footprint of MUSE and HDUV pointings in the CDFS region. The
background is a V-band image taken from the Garching-Bonn Deep Survey (GaBoDS,
Hildebrandt et al., 2006). The 60 individual MUSE-Wide pointings relevant here
are shown in bright green, and the nine MUSE-Deep fields (Bacon et al., 2017;
Inami et al., 2017, also including the UDF 10) are marked blue. The HDUV
footprint is shown as the blue shaded area and covers most of the MUSE-Wide
fields as well as MUSE-Deep.

###  2.1 Spectroscopy from MUSE

There are two main observational techniques for detecting star-forming
galaxies at redshifts $z=3-6$: Using photometry to detect the Lyman break,
resulting in the detection of LBGs and using spectroscopy to detect emission
lines, at those redshifts mostly the Ly$\alpha$ line, resulting in LAEs. The
former method is potentially biased against finding strong LyC leaker
candidates, which would make the Lyman break less pronounced. A high LyC
escape fraction could result in an unusual SED and a less steep Lyman break,
making it hard to identify LyC leakers (e.g. Cooke et al., 2014; Vanzella et
al., 2016, but see also Steidel et al., 2018 who do find LyC emission in
LBGs). Furthermore, the search for LyC emission is often complicated by the
lack of precise redshift information (e.g. when using the Lyman break in
photometry). Therefore, in order to make use of the redshift information from
spectroscopy gained from the Ly$\alpha$ line, and because this line is
expected to correlate with LyC emission, we use a pre-selection of LAEs from
MUSE to start our search.

Integral field spectroscopy is ideal for finding emission line objects such as
LAEs without proxies (such as HST detections in the rest-UV), allowing for a
relatively unbiased sample selection (when it comes to galaxy properties such
as star formation rate, other emission lines or UV properties, which are often
used at low redshifts to select LyC leaker candidates). Thanks to the
Ly$\alpha$ line, we also have a direct estimate of the redshift, which can be
slightly shifted with respect to the systemic redshift (e.g. Verhamme et al.,
2018) but is good enough to search for LyC emission in broadband data.

Therefore, we use as the basis for our search a sample of high-redshift LAEs
that were found with MUSE (Bacon et al., 2010) installed at the ESO Very Large
Telescope (VLT) in Chile. MUSE has a spectral range covering
$4750\,\mbox{\AA}$ to $9350\,\mbox{\AA}$, allowing the detection of Ly$\alpha$
at $1215.67\,\mbox{\AA}$ in the redshift range of $2.9<z<6.7$.

Several studies within the MUSE consortium are concerned with finding emission
line galaxies, in particular LAEs, in MUSE data (e.g. Herenz et al., 2017;
Claeyssens et al., 2019). Here we use data from the MUSE-Wide survey111For the
first 44 fields, the data and data products such as cut-outs, mini-cubes and
extracted spectra as well as emission line catalogues are publicly available
and can be found at https://musewide.aip.de/project/. (Urrutia et al., 2019;
Kerutt et al., 2022), which consists of 100 MUSE pointings of one hour
exposure time, including 60 pointings in CDFS region (see Fig. 1), with a
field of view of one arcmin2 each, overlapping by $\approx 4\arcsec$. Among
the 100 pointings, there are also observations from the MUSE-Deep222The
catalogue of objects in the MUSE-Deep fields used here (and presented by Inami
et al., 2017) is available at http://cdsarc.u-strasbg.fr/viz-
bin/qcat?J/A+A/vol/page, the second data release can be found in Bacon et al.
(2023) and at https://amused.univ-lyon1.fr/. survey (Bacon et al., 2017; Inami
et al., 2017; Bacon et al., 2023), namely 9 mosaic fields with 10 hours
exposure time each and one field of 31 hours (the MUSE UDF-10, see Bacon et
al., 2017, 2023).

The emission line catalogue for MUSE-Wide (which here also includes the MUSE-
Deep observations) is constructed in a consistent way, using the two spatial
and one spectral dimension of integral field spectroscopy, assuring the
inclusion of potentially extended Ly$\alpha$ halos (which are ubiquitous, see
Wisotzki et al., 2016; Leclercq et al., 2017; Wisotzki et al., 2018). The
emission lines are detected using a matched filtering approach (Herenz &
Wisotzki, 2017) and each object is classified by at least three people using a
graphical user interface (Kerutt, 2017) that allows accessing all information
from the MUSE datacube and additional HST broadband images. Criteria for
classifying an object as an LAE were a typical asymmetric or double-peaked
line (e.g. Verhamme et al., 2006), no other emission lines visible in the MUSE
spectrum (except if they match the redshift), and no strong emission in a band
bluewards of the Lyman break. Especially the last criterion might potentially
bias the sample against strong LyC leakers, which is why we only used it as an
indication in case there are multiple potential UV continuum counterparts
visible in the HST data. In those cases, we chose the one with the more
pronounced Lyman break. This assures the exclusion of low-redshift
interlopers, which is especially important for our search for LyC emission
(see Sect. 3 where we describe further measures to reduce interlopers, such as
confidence flags). However, there is a remaining possibility that a few,
particularly strong LyC leakers could potentially have been falsely assigned
lower redshifts. The $15\%$ completeness limit for the detection of LAEs in
MUSE-Wide is
$\mathrm{log_{10}(L}_{\mathrm{Ly}\alpha}\,\mathrm{[erg/s]})\approx 42.2$
($\mathrm{log_{10}(L}_{\mathrm{Ly}\alpha}\,\mathrm{[erg/s]})\approx 42.7$) at
a redshift of $z\approx 3$ ($z\approx 6.5$), as discussed in Herenz et al.
(2019), who find a characteristic Ly$\alpha$ luminosity for their luminosity
function of $\mathrm{log}L^{*}\mathrm{[erg/s]}=42.2^{+0.22}_{-0.16}$.

This catalogue of emission line objects was the basis for the selection of
LAEs from Kerutt et al. (2022) used in this paper, which consists of 1920 LAEs
identified in the MUSE-Wide and -Deep data (excluding known AGN). We also use
the properties of the Ly$\alpha$ emission lines of those LAEs, which were
measured from spectra extracted from the MUSE data, constructed by summing
each spectral layer, and weighted by the wavelength-dependent Moffat (1969)
profile that best describes the MUSE point spread function (PSF, see Urrutia
et al., 2019). The positions for the spectral extractions were the highest SN
peak of the Ly$\alpha$ emission in MUSE. From these spectra we can gain
information on the full width at half maximum (FWHM), the peak separation (in
case of a double peak) and, in combination with broadband HST data, the rest-
frame equivalent width ($\mathrm{EW}_{0}$) of Ly$\alpha$, all of which we take
from Kerutt et al. (2022).

Figure 2: HST bands and MUSE throughput used in this study. The $x$-axis shows
the observed wavelength in Å, the left $y$-axis shows the redshift and
corresponds to the black dashed and solid lines, and the right $y$-axis shows
the filter throughput and corresponds to the coloured filter curves. The solid
line shows where the Ly$\alpha$ line would fall depending on redshift, and the
dashed line shows the Lyman limit at $912\,\mbox{\AA}$, delimiting the LyC,
which is bluewards of this line. Only the redshift range where Ly$\alpha$ is
visible for MUSE is shown. The two bands with the lowest wavelengths
(WFC3/UVIS F275W and F336W) are from the HDUV legacy survey and are used in
this study to identify LyC emission.

###  2.2 The HDUV survey

To identify possible LyC emission in our sample of MUSE LAEs, we use data from
the HDUV legacy survey, (Oesch et al., 2015, 2018), which includes the UVUDF
(Teplitz et al., 2013; Rafelski et al., 2015) and covers an area of $\approx
100\,\mathrm{arcmin}^{2}$ in the two GOODS/CANDELS-Deep (Giavalisco et al.,
2004; Koekemoer et al., 2011) fields and an area of
$43.4\,\mathrm{arcmin}^{2}$ in GOODS-S (Oesch et al., 2018). The HDUV survey
provides the two HST filters WFC3/UVIS F275W and F336W (see Fig. 2), reaching
a limiting magnitude of $27.5$ to $28\,\mathrm{mag}_{\mathrm{AB}}$ at
$5\,\sigma$ in $0\aas@@fstack{\prime\prime}4$ apertures (Oesch et al., 2018)
and covering a wavelength range from $2435-3032\,\mbox{\AA}$ and
$3096-3639\,\mbox{\AA}$, respectively. Thus, the Lyman break at
$912\,\mbox{\AA}$ falls in the HST filter band WFC3/UVIS F275W at a redshift
range of $z=1.67-2.3$ and in the band WFC3/UVIS F336W at a redshift range of
$z=2.4-3.0$, in other words, at redshifts $z>2.3$ the band F275W ($z>3.0$ for
F336W) is uncontaminated by non-LyC emission (see Fig. 2). We therefore use
the HST filter band F336W to search for LyC emission from the MUSE LAEs, which
have redshifts $z>2.9$.

###  2.3 Additional HST data

To measure the LyC escape fraction, we need measurements of the LyC emission
and the UV continuum (see Sect. 4). Since we will base our $f_{\mathrm{esc}}$
measurements on SED fitting (see Sect. 4.4), we need flux measurements in
various wavelength bands, which we get through aperture photometry. Even
though most of our objects already have measurements in e.g. the 3DHST
catalogue (Skelton et al., 2014), we want to keep the flux measurement method
consistent in all bands, including the LyC band WFC3/UVIS F336W. Therefore, we
perform aperture photometry (using an aperture of
$0\aas@@fstack{\prime\prime}35$ radius, see Sect. 3.2) not only for LyC but
also for other HST bands, including WFC3/UVIS F275W, ACS F435W, F606W, F775W,
F814W, F850LP, WFC3/IR F105W, F125W, F140W and F160W. Except for WFC3/UVIS
F275W and F336W, the data was taken from the Hubble Legacy Fields (HLF333The
HST images for all bands except the two HDUV ones WFC3/UVIS F275W and F336W
were taken from https://archive.stsci.edu/prepds/hlf/., Illingworth et al.,
2016) GOODS-S data release 2.0. For the aperture photometry, we need to
account for the different PSFs in the bands, which increase in size with
increasing wavelength. Thus, the filter band with the largest PSF is the
WFC3/IR F160W band. For measuring the photometry, we use images of each filter
band that are matched to the PSF in WFC3/IR F160W band444The PSF matched
images can be found here:
https://archive.stsci.edu/hlsps/hlf/v2.0/60mas_conv/. To create these images,
Whitaker et al. (2019) use a linear combination of Gaussian-weighted Hermite
polynomials based on a stack of isolated stars to determine the PSF in WFC3/IR
F160W. The images of the other bands were then convolved with the matching
kernel. We take the positions of the UV continuum counterparts from Kerutt et
al. (2022), which were determined by considering all detections in the band
ACS F814W within a radius of $0\aas@@fstack{\prime\prime}5$ from the maximum
SN of the Ly$\alpha$ emission in MUSE. Then at least two people examined the
potential counterparts, taking other bands into account, and decided on the
correct counterparts (see Kerutt et al., 2022 for details).

##  3 Sample Selection

In this section, we describe how we select candidates to detect LyC emission.
We then divide the found candidates into two groups, a gold
($\mathrm{SN}_{\mathrm{F336W}}>3$) and a silver
($\mathrm{SN}_{\mathrm{F336W}}>2$) sample, indicating their quality as
potential LyC leakers.

###  3.1 Redshift and Quality Cuts

To construct the sample, we use the redshift information from the MUSE
catalogue (Kerutt et al., 2022) in the MUSE-Wide fields as a basis to look for
LyC emission in the HDUV data. For the LAEs we find in MUSE, which are in a
redshift range $2.9<z<6.7$, the Lyman break at $912\,\mbox{\AA}$ (shortwards
of which is the LyC) lies in the observed wavelength range
$3563-7014\,\mbox{\AA}$, which means the LyC falls into the wavelength range
of the filter WFC3/UVIS F336W (see the dashed line in Fig. 2). For redshifts
$z<3$, there is an overlap of non-LyC emission in the WFC3/UVIS F336W band,
which is why we apply a redshift cut of $z>3$ for our sample selection.
Another restriction to keep in mind is the IGM transmission, which decreases
drastically towards higher redshifts (e.g. Inoue et al., 2014), from a mean
IGM transmission at $900\,\AA$ of $\tau_{\mathrm{IGM}}=0.56$ at $z=3$ to
$\tau_{\mathrm{IGM}}=0.17$ at $z=4.5$, which is why we do not consider any
objects beyond $z>4.5$ for our search for LyC leaker candidates. In this
redshift range ($3<z<4.5$), there are $743$ LAEs in the MUSE-Wide survey that
overlap with the GOODS-S part of the HDUV survey. The redshifts for most of
those objects are based on a single line, the Ly$\alpha$ line. Due to
radiative transfer processes in the neutral hydrogen in the ISM, it can be
asymmetric or double-peaked, sometimes mimicking the O II doublet. Therefore,
low SN Ly$\alpha$ lines are not easy to identify, which is why we add a cut in
the confidence parameter given in the MUSE-Wide survey. This subjective
parameter represents the certainty of the redshift classification, ranging
from zero (not certain at all) to three (very sure). We exclude objects with a
confidence below two, which leaves us with a final sample of $621$ LAEs at
$z=3.0-4.5$.

###  3.2 Measuring the LyC emission and signal-to-noise

To determine possible LyC leakers, we measure the flux in the WFC3/UVIS F336W
band at the same position as the UV continuum of each object. Since we perform
SED fitting later (see Sect. 4.4), we apply the same procedure for all other
bands as well.

We use an aperture of $0\aas@@fstack{\prime\prime}35$ radius to perform
aperture photometry on the images convolved to the WFC3/IR F160W PSF, after
subtracting the local background (using the median of the 3 sigma clipped
values in a cutout of $2^{\prime\prime}\times 2^{\prime\prime}$ arcseconds,
excluding the aperture itself). Following Skelton et al. (2014), we apply an
aperture correction $21\%$ to total fluxes to account for an aperture with a
radius of $0\aas@@fstack{\prime\prime}35$.

For the position of the apertures, we use the same location as the UV
continuum emission, which is where the LyC emission is expected as well.
Unlike e.g. Ly$\alpha$, which is resonantly scattered, it will be detectable
in the line of sight where it was first emitted. Therefore, we do not use the
Ly$\alpha$ positions from the MUSE detections, but rather the UV continuum
counterpart positions, which were determined using the Galfit software (Peng
et al., 2002, 2010) in the ACS F814W band (Giavalisco et al., 2004; Koekemoer
et al., 2011; Grogin et al., 2011), which goes down to a magnitude of 28.92
$\mathrm{mag}_{\mathrm{AB}}$ in MUSE-Wide and 29.16
$\mathrm{mag}_{\mathrm{AB}}$ in MUSE-Deep (for a $2\,\sigma$ detection, see
Kerutt et al., 2022 for more details). This implies that we only consider
objects that have a detection in the UV continuum, but since the WFC3/UVIS
F336W filter band is a bit shallower and only reaches a magnitude of 28
$\mathrm{mag}_{\mathrm{AB}}$, this criterion will not exclude any potential
LyC leakers. In case the LyC emission originates from a single star-forming
region in the galaxy, its position can still be slightly offset with respect
to the maximum of the UV continuum, but still be captured by our
$0\aas@@fstack{\prime\prime}35$ aperture photometry.

To measure the signal-to-noise (SN) of the LyC flux, we need to take into
account that the noise in the HDUV data is correlated. Therefore, we determine
the local noise properties by applying 100 random, non-overlapping apertures
of the same size as for the flux measurement to estimate the noise in an area
of $8^{\prime\prime}\times 8^{\prime\prime}$ around the object, which we do in
the original, unconvolved data. To not be influenced by a real signal, we mask
objects and bright/hot pixels. For the latter, we use a SN cut of three for
individual pixels. To mask objects, we use the segmentation map created in the
WFC3/IR F160W band in the 3DHST survey (Skelton et al., 2014). This band has a
wider PSF than WFC3/UVIS F336W and lower redshift interlopers tend to be
brighter and more extended in this wavelength range. To be sure not to include
any emission that might be more extended in WFC3/UVIS F336W than in the
WFC3/IR F160W band the segmentation map is based on, we grow the mask by one
pixel in each direction. The standard deviation of the 100 apertures provides
us with an estimate of the noise and the SN is the flux measurement in the
original, unconvolved data over this standard deviation.

We set the SN cut for our pre-selection to the relatively low value of SN=2,
since we have the prior knowledge from the MUSE catalogues that there is an
object at the correct redshift at the position we are examining in the LyC
band WFC3/UVIS F336W (see also below in Sect. 3.3 for a detailed discussion).
Including the redshift cut $3<z<4.5$ mentioned above, this leaves us with a
total number of 42 objects to take into consideration).

###  3.3 Identification of LyC Counterparts and Contamination

Having measured the SN ratio in the LyC, we use this as an indicator of the
possible presence of a LyC leaker. However, even with a SN cut of two, the
signal could still be created by various other influences such as random
noise, close neighbours, or low-redshift interlopers. Therefore, we use
several criteria to narrow down our sample of 42 possible LyC leakers and also
to separate our remaining candidates into likely leakers, which we call the
gold sample, and possible leakers, the silver sample, as done for example in
Fletcher et al. (2019). The difference between the gold and silver samples is
that for a silver object, an SN cut of two in the WFC3/UVIS F336W band is
applied, while for a gold object, we require $\mathrm{SN}>3$.

We use several criteria to make sure our sample of LyC leakers is not
contaminated by low-redshift interlopers or noise, as discussed in many
pioneering studies searching for individual LyC leakers (see e.g. Shapley et
al., 2006; Iwata et al., 2009; Vanzella et al., 2010b; Nestor et al., 2011,
2013; Mostardi et al., 2015). For our study, we have the ideal situation to be
able to use integral field spectroscopy to obtain reliable redshifts, in
addition to the high resolution of the space-based HST data.

We visually vet all 42 potential LyC leakers using red-green-blue (rgb) images
(see Figs. 3 and 4) of three HST filter bands (ACS F606W, WFC3/IR F125W and
F160W), overlaid by SN contours in the WFC3/UVIS F336W band. In these rgb
images, we can see if the object has a similar morphology in the infrared as
in the UV and if it matches the detection in the LyC band. If there are
several clumps with different colours, this is an indication of a potential
chance alignment of different objects (which was the case for six objects).
Different areas of the same object can have different colours depending on its
properties. However, we try to distinguish between this case and clearly
different clumps. If they do have similar colours, this could indicate an
ongoing merger, which could trigger star formation and thus also the
production of LyC photons.

Figure 3: RGB image for one of the LyC leaker candidates. The three images
that were combined for the main image are shown to the right of each pane and
are the HST bands WFC3/IR F160W (red), F125W (green) and ACS F606W (blue). The
white contours in the main panels show the signal in HST band WFC3/UVIS F336W,
where the LyC emission can be found. The contours show SN of 3, 2 and 1
(strongest to lightest) for an image smoothed with a Gaussian kernel of
$\sigma=1$ pixel. The yellow dotted contours show the extent of the Ly$\alpha$
emission in the MUSE data, also smoothed with a Gaussian kernel of $\sigma=1$
pixel. It has to be kept in mind, that the PSF of MUSE is larger than that of
HST, however, it has been shown that Ly$\alpha$ emission is usually ten times
more extended than the UV continuum (e.g. Wisotzki et al., 2016; Leclercq et
al., 2017; Wisotzki et al., 2018). The red cross in the main panel shows the
position of the pixel with the highest S/N in the Ly$\alpha$ emission found
with MUSE.

Figure 4: Same as Fig. 3 for the rest of the gold sample LyC leaker
candidates.

In addition to using rgb images, we look at cutouts of the objects in the
WFC3/UVIS F336W band and use the MUSE catalogue from Kerutt et al. (2022) as
well as the catalogue from Skelton et al. (2014) to highlight all other known
objects in the same cutout (see Figs. 13 and 14). This allows us to exclude
contamination by close, projected neighbours, that might only be visible in
the MUSE catalogue through their emission lines (which was the case for 2
objects). Another indication of low-redshift contamination is a strong signal
in the band WFC3/UVIS F275W with SN¿3. Even though LyC can be stronger at
lower wavelengths, here we expect the LyC to be significantly reduced by the
IGM (which is true for all but seven objects, three of which were obvious
contaminants).

Another criterion is the overlap between the signal found in LyC in WFC3/UVIS
F336W and the UV continuum, in this case, ACS F775W. For this, we use
segmentation maps made with SExtractor (Bertin & Arnouts, 1996) and determine
the overlap between the two bands by the percentage of pixels belonging to the
LyC emission which are also above the segmentation map threshold of $1\sigma$
in the UV continuum. A non-zero overlap is required for us to consider an
object a LyC leaker candidate.

This leaves us with a final sample of seven silver candidates and five gold
candidates. An overview of the final sample with IDs from different existing
catalogues and redshifts can be found in table 1.

###  3.4 Comparison to previous Studies

There have been previous studies of LyC leakage in the same region, one by
Saxena et al. (2022), who found 11 LyC leaker candidates in the CDFS using
narrow-band imaging data from Guaita et al. (2016) and the U-band from Nonino
et al. (2009). Of their 11 candidates, five are in the field of view of the
HDUV with MUSE-Deep IDs 1087, 8035, 6666, and 7820 (one is not in the MUSE
catalogue, see table 1 for the IDs). Of those, only two have potential
detections in the HDUV WFC3/UVIS F336W band, IDs 1087 and 8035. The former is
among our gold sample, the latter (ID 13385 in Saxena et al., 2022 at
$z=3.431$ with a SN in WFC3/UVIS F336W of $\approx 1.8$) was assigned a lower
redshift of $z=0.523$ in the MUSE-Deep catalogue (and was not detected in
MUSE-Wide), although with a low confidence of 1, indicating that the redshift
determination from MUSE might not be reliable. However, the MUSE redshift of
the latter agrees quite well with the photometric redshift of Skelton et al.
(2014) at $z=0.586$ (ID 25631).

Another recent study in the same field was done by Rivera-Thorsen et al.
(2022), also using data from HST and MUSE. Unlike previous studies, they first
run a source detection software on three HST UV filters (F225W, F275W and
F336W) and then determine whether the detected signal could be LyC emission.
They find six new LyC leaker candidates and confirm the object from Saxena et
al. (2022) that is also among our candidates (ID 3052076 with a
SN${}_{336}=4.07$). Of the other six, three are below our redshift cut of
$z>3$ imposed by the wavelength range of the WFC3/UVIS F336W filter band that
we use for detecting LyC emission. Another of their candidates is above our
redshift upper limit of $z<4.5$ in our catalogue but might be a different,
close-by object. The last candidate does not have a close enough counterpart
in our sample based on Kerutt et al. (2022).

Table 1: Overview of candidates

IDMW | IDMD | IDSkel | IDRaf | IDGuo | RA | Dec | zMW | zMD | zSkel | cMW | cMD | SN336
---|---|---|---|---|---|---|---|---|---|---|---|---
gold candidates |  |  |  |  |  |  |  |  |  |  |  |
1181371 | 7193 | 24480 | 6470 | - | 53.1358 | -27.7955 | 3.084 | 3.085 | 3.441 | 2 | 2 | 5.01
3052076 | 1087 | 24193 | 3506 | 12448∗ | 53.1679 | -27.798 | 3.457 | 3.462 | 1.711 | 3 | 3 | 4.07
109004028 | - | 15601 | - | 7570 | 53.0994 | -27.8392 | 3.267 | - | 1.492 | 2 | - | 3.03
122021111 | - | 16523 | - | 8005 | 53.1389 | -27.8354 | 3.794 | - | 0.724 | 3 | - | 3.57
126049137 | - | 20189 | - | 10131 | 53.2042 | -27.8172 | 4.426 | - | 1.222 | 2 | - | 4.39
silver candidates |  |  |  |  |  |  |  |  |  |  |  |
1521589 | 7169 | 26130 | 4873 | 13707 | 53.1283 | -27.7887 | 3.152 | 3.155 | 1.537 | 3 | 2 | 2.69
3452147 | 2134 | 23769 | 3431 | 12256 | 53.1541 | -27.7988 | 3.521 | 3.524 | 4.099 | 3 | 2 | 2.6
4062373 | 1360 | 27696 | 37765 | 14691 | 53.1792 | -27.7829 | 3.663 | 3.666 | 1.014 | 3 | 2 | 2.21
4172404 | 7121 | 27408 | 6974 | 14469 | 53.1851 | -27.7839 | 3.672 | 3.675 | 3.759 | 3 | 2 | 2.88
5622786 | 7137 | (20161?) | (781?) | - | 53.1604 | -27.8174 | 4.005 | 4.007 | - | 3 | 2 | 2.05
119004004 | - | 16269 | - | 7896 | 53.1891 | -27.8363 | 3.314 | - | 0.798 | 3 | - | 2.59
122032127 | - | 16000 | - | 7760 | 53.1326 | -27.8374 | 4.348 | - | 0.892 | 2 | - | 2.33

555IDMW: Identifier in the MUSE-Wide catalogue (Urrutia et al., 2019; Kerutt
et al., 2022 and Urrutia et al. in prep.). The colours indicate if the object
is placed in the gold or silver category, as explained in Sect. 3.3. IDMD:
Identifier in the MUSE-Deep catalogue (Inami et al., 2017). IDSkel: Identifier
from the 3D-HST CANDELS catalogue in Skelton et al. (2014). IDRaf: Identifier
from the UDF catalogue in Rafelski et al. (2015). IDGuo: Identifier from the
CANDELS GOODS-S catalogue in Guo et al. (2013). RA, DEC: Right ascension and
declination (Kerutt et al., 2022). zMW: redshift from MUSE-Wide, based on the
Ly$\alpha$ line, corrected to systemic using the full-width at half maximum
(FWHM) or peak separation (if the line has a double peak), based on Verhamme
et al. (2018). zMD: redshift from Inami et al. (2017). zSkel: redshift from
Skelton et al. (2014). cMW: Classification confidence (0 lowest and 3 highest)
of the redshift of the LAE based on three different investigators and a
consolidation (see Urrutia et al., 2019; Schmidt et al., 2021; Kerutt et al.,
2022). cMD: confidence from Inami et al. (2017). SN336: Signal-to-noise in the
LyC band WFC3/UVIS F336W. The object marked with ∗ is already found in Saxena
et al. (2022) and further discussed in Rivera-Thorsen et al. (2022).

##  4 Measuring Escape Fractions

In this section we will measure escape fractions in three different ways:
assuming a fixed ratio between the UV continuum and the LyC, using a Binary
Population and Spectral Synthesis (BPASS, Eldridge et al., 2017; Stanway &
Eldridge, 2018; Byrne et al., 2022, version 2.3), and fitting the SEDs using
the Code Investigating GALaxy Emission (CIGALE, Burgarella et al., 2005; Noll
et al., 2009; Boquien et al., 2019).

###  4.1 Definition and ingredients of the escape fraction

A common definition for the relative escape fraction $f_{\mathrm{esc,rel}}$
was suggested by Steidel et al. (2001) (but see also e.g. Siana et al., 2007;
Bian et al., 2017; Grazian et al., 2018; Steidel et al., 2018):

$f_{\mathrm{esc,rel}}=\left(\frac{f_{\mathrm{LyC}}}{f_{\mathrm{UVC}}}\right)_{\mathrm{obs}}\times\left(\frac{L_{\mathrm{UVC}}}{L_{\mathrm{LyC}}}\right)_{\mathrm{int}}\times\exp(\tau^{\mathrm{IGM}}_{\mathrm{LyC}})$
(1)

This relative escape fraction is a comparison between the observed flux ratio
of the LyC and UV continuum (UVC) fluxes
$(f_{\mathrm{LyC}}/f_{\mathrm{UVC}})_{\mathrm{obs}}$ and the intrinsic
luminosity ratio $(L_{\mathrm{UVC}}/L_{\mathrm{LyC}})_{\mathrm{int}}$,
corrected for the absorption of LyC photons by neutral hydrogen in the IGM,
$\exp(\tau^{\mathrm{IGM}}_{\mathrm{LyC}})$. This definition is independent of
dust and can be related to the absolute escape fraction by taking into account
the dust attenuation of the non-ionising UV continuum, $A_{\mathrm{UVC}}$:

$f_{\mathrm{esc,abs}}=f_{\mathrm{esc,rel}}\times 10^{-0.4\,A_{\mathrm{UVC}}}$
(2)

This definition of the escape fraction depends on knowledge of the dust
attenuation, for example, obtained through SED fitting or several photometric
bands. If there is no significant dust attenuation, the relative and absolute
escape fractions are the same (see for example the galaxy in Shapley et al.,
2016). Typically, the escape fraction is defined as the fraction between the
flux values at $900\,\mbox{\AA}$ for the LyC and $1500\,\mbox{\AA}$ for the UV
continuum.

For our escape fraction estimates, we use three different approaches, two of
which use fixed intrinsic flux ratios
$(L_{\mathrm{UVC}}/L_{\mathrm{LyC}})_{\mathrm{int}}$ (see Sect. 4.3) and the
last one based on SED fitting (see Sect. 4.4).

###  4.2 Observed flux ratio

The observed flux ratio $(f_{\mathrm{LyC}}/f_{\mathrm{UVC}})_{\mathrm{obs}}$
can be obtained from the measurements in the LyC band WFC3/UVIS F336W, as
described in Sect. 3.2. Depending on the redshift of the individual
candidates, the two bands used for the LyC and UV continuum respectively cover
different wavelength ranges and their effective wavelengths are not
necessarily at $900\,\mbox{\AA}$ and $1500\,\mbox{\AA}$. In the third column
in Table 2, we show the fraction
$(f_{\mathrm{LyC}}/f_{\mathrm{UV}})_{\mathrm{obs}}$ based directly on the
measurements in the respective bands used for detecting the LyC and UV
continuum. This fraction is free of any model assumptions and only contains
measured values. It is used to compute the escape fraction for the fixed LyC
to UV continuum ratio, while the escape fractions based on the BPASS model are
corrected to $900\,\mbox{\AA}$ and $1500\,\mbox{\AA}$ (see Sect. 4.3 below).
It has to be noted that except for the first object, all fractions are below
one, meaning the LyC flux is lower than the UV continuum, as expected. The
first object presents an interesting case and needs a special explanation. It
could be that the nebular continuum at the Lyman break is responsible for
creating a Lyman bump (see Inoue2010), resulting in an excess of LyC emission.

Table 2: LyC leaker candidates, various measurements of fesc

ID | z | $\left(\frac{f_{\mathrm{LyC}}}{f_{\mathrm{UV}}}\right)_{\mathrm{obs}}$ | TIGM | AV | method | $\left(\frac{L_{1500}}{L_{900}}\right)_{\mathrm{int}}$ | f${}_{\mathrm{esc,rel}}[\%]$ | f${}_{\mathrm{esc,rel}}^{\mathrm{IGM}}[\%]$ | f${}_{\mathrm{esc,abs}}^{\mathrm{IGM}}[\%]$
---|---|---|---|---|---|---|---|---|---
gold candidates |  |  |  |  |  |  |  |  |
1181371 | 3.08 | 1.17$\pm$0.22 | 0.66 | 0.42$\pm$0.11 | fixed | 3 | 350$\pm$66 | 530$\pm$100 | 360$\pm$77
|  |  |  |  | BPASS | 1.72 | 102$\pm$19 | 155$\pm$28 | 105$\pm$22
|  |  |  |  | CIGALE | 1.36 |  |  |  88$\pm$7 (90$\pm$6)
3052076 | 3.46 | 0.10$\pm$0.02 | 0.46 | 0.78$\pm$0.01 | fixed | 3 |  30$\pm$6 |  66$\pm$13 |  32$\pm$6
|  |  |  |  | BPASS | 1.72 |  17$\pm$3 |  37$\pm$7 |  18$\pm$3
|  |  |  |  | CIGALE | 4.32 |  |  |  23$\pm$5 (40$\pm$0)
109004028 | 3.27 | 0.14$\pm$0.04 | 0.57 | 0.58$\pm$0.13 | fixed | 3 |  43$\pm$12 |  76$\pm$21 |  45$\pm$14
|  |  |  |  | BPASS | 1.72 |  13$\pm$4 |  23$\pm$7 |  13$\pm$4
|  |  |  |  | CIGALE | 4.64 |  |  |  34$\pm$10 (47$\pm$6)
122021111 | 3.79 | 0.06$\pm$0.01 | 0.30 | 0.03$\pm$0.04 | fixed | 3 |  19$\pm$4 |  64$\pm$13 |  62$\pm$13
|  |  |  |  | BPASS | 1.72 |  12$\pm$3 |  40$\pm$10 |  39$\pm$10
|  |  |  |  | CIGALE | 3.98 |  |  |  66$\pm$11 (77$\pm$12)
126049137 | 4.43 | 0.27$\pm$0.12 | 0.05 | 0.89$\pm$0.24 | fixed | 3 |  81$\pm$37 | 1473$\pm$673 | 649$\pm$329
|  |  |  |  | BPASS | 1.72 |  49$\pm$22 | 891$\pm$400 | 393$\pm$196
|  |  |  |  | CIGALE | 1.39 |  |  |  69$\pm$10 (76$\pm$9)
silver candidates |  |  |  |  |  |  |  |  |
1521589 | 3.15 | 0.37$\pm$0.11 | 0.61 | 0.10$\pm$0.07 | fixed | 3 | 110$\pm$34 | 180$\pm$56 | 164$\pm$52
|  |  |  |  | BPASS | 1.72 |  26$\pm$8 |  42$\pm$13 |  39$\pm$12
|  |  |  |  | CIGALE | 1.55 |  |  |  79$\pm$15 (86$\pm$10)
3452147 | 3.52 | 0.48$\pm$0.15 | 0.46 | 0.70$\pm$0.18 | fixed | 3 | 143$\pm$46 | 314$\pm$101 | 165$\pm$60
|  |  |  |  | BPASS | 1.72 |  69$\pm$22 | 152$\pm$48 |  80$\pm$28
|  |  |  |  | CIGALE | 2.70 |  |  |  47$\pm$14 (56$\pm$14)
4062373 | 3.66 | 0.35$\pm$0.13 | 0.35 | 0.92$\pm$0.11 | fixed | 3 | 106$\pm$40 | 301$\pm$113 | 129$\pm$50
|  |  |  |  | BPASS | 1.72 |  61$\pm$23 | 173$\pm$65 |  74$\pm$28
|  |  |  |  | CIGALE | 4.93 |  |  |  74$\pm$13 (82$\pm$13)
4172404 | 3.67 | 0.52$\pm$0.15 | 0.35 | 0.49$\pm$0.22 | fixed | 3 | 156$\pm$45 | 443$\pm$128 | 282$\pm$100
|  |  |  |  | BPASS | 1.72 |  94$\pm$27 | 267$\pm$77 | 170$\pm$60
|  |  |  |  | CIGALE | 1.75 |  |  |  31$\pm$15 (47$\pm$15)
5622786 | 4.00 | 0.36$\pm$0.15 | 0.19 | 0.00$\pm$0.00 | fixed | 3 | 108$\pm$44 | 555$\pm$225 | 555$\pm$225
|  |  |  |  | BPASS | 1.72 |  62$\pm$25 | 319$\pm$128 | 319$\pm$128
|  |  |  |  | CIGALE | 3.86 |  |  |  53$\pm$11 (71$\pm$10)
119004004 | 3.31 | 0.18$\pm$0.06 | 0.57 | 0.81$\pm$0.26 | fixed | 3 |  54$\pm$19 |  96$\pm$34 |  45$\pm$19
|  |  |  |  | BPASS | 1.72 |  28$\pm$10 |  51$\pm$18 |  24$\pm$10
|  |  |  |  | CIGALE | 2.81 |  |  |  26$\pm$14 (50$\pm$14)
122032127 | 4.35 | 0.17$\pm$0.06 | 0.08 | 0.45$\pm$0.28 | fixed | 3 |  52$\pm$19 | 674$\pm$246 | 445$\pm$199
|  |  |  |  | BPASS | 1.72 |  23$\pm$9 | 298$\pm$117 | 197$\pm$92
|  |  |  |  | CIGALE | 1.25 |  |  |  77$\pm$14 (84$\pm$12)

666Overview of different escape fraction values based on three different
methods. (fLyC/f${}_{\mathrm{UV}})_{\mathrm{obs}}$: observed flux ratio
between the LyC band WFC3/UVIS F336W and the UV continuum in ACS F775W, given
in frequency space. TIGM: LyC transmission in the IGM at the redshift of the
object, using the mean of the $5\%$ lines of sight with the highest
transmission. AV: dust extinction from the SED fit of CIGALE.
(L1500/L${}_{900})_{\mathrm{int}}$: intrinsic luminosity ratio between the UV
continuum at $1500\,\mbox{\AA}$ and the LyC at $900\,\mbox{\AA}$. The three
different methods (in the method column) correspond to a fixed value of 3,
which is often used in the literature, a value of 1.72 based on a BPASS model
with a young age of 3 Myr and a metallicity of $\mathrm{Z}=0.002$, and the
third value is the ratio measured from the SED model of CIGALE, using the
stellar unattenuated model without IGM absorption. fesc,rel: relative escape
fraction for two different methods. The first corresponds to the intrinsic
flux ratio of 3 and uses the fluxes directly measured from the respective
bands (the third column in this table). The second escape fraction is based on
the BPASS model and uses flux measurements which are corrected to the
wavelengths $1500\,\mbox{\AA}$ and $900\,\mbox{\AA}$ using the model. There is
no relative escape fraction without IGM absorption from the SED fits using
CIGALE. f${}_{\mathrm{esc,rel}}^{\mathrm{IGM}}$: relative escape fractions
taking into account the IGM absorption shown in the fourth column of this
table. The three escape fractions correspond again to the different methods,
with the third one being the sightline-dependent escape fraction measured by
CIGALE. f${}_{\mathrm{esc,abs}}^{\mathrm{IGM}}$: absolute escape fractions
taking the dust extinction AV (fifth column of this table) into account, again
for the three different methods. Note that since CIGALE takes the dust
attenuation in the UVC already into account, we directly get the absolute
escape fraction. The value in brackets denotes the absolute escape fraction
from CIGALE when the dust description is extended below $912\AA$ but only at a
$10\%$ level.

###  4.3 Intrinsic ratio of ionising to non-ionising UV luminosities

The intrinsic ratio $(L_{\mathrm{UVC}}/L_{\mathrm{LyC}})_{\mathrm{int}}$
between the LyC and the UV continuum at $900\,\mbox{\AA}$ and
$1500\,\mbox{\AA}$ is not known as it cannot be measured directly but must be
inferred from stellar population models. Common values in the literature range
from $(L_{\mathrm{1500}}/L_{\mathrm{900}})_{\mathrm{int}}=3$ (frequency space
$\nu$, with luminosity given in erg/s/Hz, see e.g. Steidel et al., 2001;
Grazian et al., 2016, 2017; Japelj et al., 2017; Marchi et al., 2018) to $5$
(see e.g. Naidu et al., 2018) and even as high as $6\\--8$ (see e.g. Siana et
al., 2007), which depends on the underlying assumptions on the star-formation
history, metallicity, initial mass function (IMF) and age.

We use two fixed values for the intrinsic ratio. The first is chosen to
compare to results from the literature as
$(L_{\mathrm{1500}}/L_{\mathrm{900}})_{\mathrm{int}}=3$, which is what the
first relative escape fraction fesc,rel for each object in table 2 is based
on. For the second fixed value, we use an extreme BPASS model with a young age
of 3 Myr, binary evolution, an IMF upper mass cutoff of
$\mathrm{M}_{\mathrm{max}}=300\,\mathrm{M}_{\odot}$ and a metallicity of
$\mathrm{Z}=0.002$, which gives
$(L_{\mathrm{1500}}/L_{\mathrm{900}})_{\mathrm{int}}=1.72$. The young age of
the model is chosen such that we expect a large Lyman continuum luminosity
(see e.g. Rosdahl et al., 2022), but not too young in order to have enough
time to form escape channels for LyC. This way, the resulting escape fraction
measurements (the second fesc,rel in table 2) can be seen as a lower limit for
our objects. The subsolar metallicity we choose here is in agreement with
observations of high redshift, star-forming galaxies, see e.g. Steidel et al.,
2016; Cullen et al., 2020; Reddy et al., 2022, and matches reasonably well
with the results of the SED fitting for most objects (see table 4).

Using this model, we can correct the measured flux values in the LyC and UV
continuum bands to the values at $900\,\mbox{\AA}$ and $1500\,\mbox{\AA}$. For
this, we convolve the model spectrum with the filter band throughput (de-
redshifted to rest-frame wavelength) and get the flux ratio between the
effective wavelength corrected to rest-frame and the desired wavelength. This
ratio is used to correct the measured flux values, which are then used for the
escape fraction. The resulting escape fraction $f_{\mathrm{esc,rel}}$ is the
second value in the eighth column in the table 2.

###  4.4 SED Fitting

To better understand the properties of our objects and to obtain our third
measurement of the escape fraction, we model their SEDs using CIGALE
(Burgarella et al., 2005; Noll et al., 2009; Boquien et al., 2019). In this
way, we can take the additional information from other photometric bands into
account when modelling the underlying stellar population of the galaxy,
instead of assuming the same model for all objects as we did in Sect. 4.3. For
the SED fitting, we use the HST bands F275W, F336W, F435W, F606W, F775W,
F814W, F850LP, F105W, F125W, F140W, F160W as well as the HAWKI K-band and the
IRAC channels 1, 2, 3 and 4.

The results for the five best LyC leaker candidates (gold sample) can be seen
in Figs. 5 and 6. For the dust attenuation law, we use the CIGALE module
dustatt_modified_starburst, based on the Calzetti et al. (2000) starburst
attenuation curve. The model allows for a more flexible approach to the dust
attenuation (Boquien et al., 2019, 2022), especially by changing the slope
using the $\delta$ parameter and including options for the extinction of
emission lines using an extinction curve based on the Milky Way (MW, Cardelli
et al., 1989; O’Donnell, 1994), and Small and Large Magellanic Cloud (SMC,
LMC, based on Pei, 1992). It has been shown for example in Reddy et al., 2018,
that for high-redshift LAEs, an SMC-like dust curve would be more realistic
for the stellar continuum. To get as close as possible to the SMC curve, we
use initial values for the $\delta$ parameter of $\delta\approx-0.5$ (as shown
in Boquien et al., 2022), but leave the parameter-free for CIGALE to fit. When
it comes to the extinction of emission lines, a MW-like dust curve is adapted
here. However, we have verified that the choice of absorption curve for
emission lines does not change the escape fraction results significantly. For
the star formation history, we use the module sfh2exp, which consists of two
exponentials, one for the long-term star formation and one for the more recent
starburst. We use a Chabrier initial mass function with the module bc03 (based
on Bruzual & Charlot, 2003), leaving the metallicity as a free parameter (see
table 4).

The results of the SED fitting can be found in the appendix in table 4. We
note that our candidates have quite high specific star formation rates, which
could contribute to their high escape fractions.

###  4.5 IGM absorption

The last factor in the escape fraction calculation is the correction for the
IGM absorption of LyC photons, using the IGM transmissivity at
$900\,\mbox{\AA}$, $\tau_{\mathrm{900,\mathrm{IGM}}}$. Here we use the
modelled line of sight transmissions from Inoue et al. (2011) and Inoue et al.
(2014) for the different redshifts of our candidates. However, at redshifts
above 3, a significant fraction of lines of sight through the IGM are already
completely opaque for LyC photons. Using the median of the distribution of
possible IGM transmission lines would result in a value of zero transmission.
We argue that since we do detect LyC emission, this introduces a selection
effect in the lines of sight, favouring ones with high transmission of LyC
emission. For all escape fraction measurements we, therefore, choose an
arbitrary cut in the transmission values, using the mean of only the highest
$5\%$ of in the $10,000$ lines of sight from Inoue et al. (2014) at the
specific redshift of the objects. The IGM transmission values
$T_{\mathrm{IGM}}=\exp(-\tau^{\mathrm{IGM}}_{\mathrm{900}})$ are given in the
fourth column in table 2. However, CIGALE assumes that there is no IGM
transmission below the Lyman break, which is why we adjust the IGM
transmission as explained above, using only the highest $5\%$.

Thus, we obtain IGM corrected escape fractions
$f_{\mathrm{esc,rel}}^{\mathrm{IGM}}$, see the second to last column in table
2, for the assumption that the IGM transmission is high for all of our LyC
leaker candidates. Following a similar thought, we assume that our LyC leaker
candidates are not affected by absorption in the CGM. In the literature, the
effects of IGM and CGM are often not treated separately, however, it has been
shown by Muzahid et al. (2021) that LAEs at similar redshifts to ours, live in
extended neutral hydrogen environments, up to a distance of $\approx
500\,\mathrm{km/s}$ or $\approx 7$ virial radii. This would make it harder for
LyC emission to escape into the IGM. However, their results show that LAEs in
groups are more likely to have extended neutral hydrogen around them, hinting
at large-scale structures. While we cannot exclude that our candidates could
be affected by CGM absorption, we apply the same argument as for the IGM
absorption: The absorption is stochastic and depends on the line of sight (see
Fig. 5 in Muzahid et al., 2021. The fact that we do see LyC emission means
that we are looking at a relatively free line of sight or possibly a hole or
low column density in the CGM. Just as with the IGM, the CGM coverage is
difficult to treat for individual objects. We therefore assume a similarly
unusual CGM configuration for our LyC leaker candidates as for the IGM and
interpret the correction for the IGM absorption as including the CGM
absorption as well.

###  4.6 Dust discussion

As can be seen in Figs. 5 and 6, the LyC part of the spectrum is not affected
by dust in our models, since there are no descriptions for that part. Although
it would be possible to extrapolate the dust extinction law to smaller
wavelengths, we decided against such a procedure for our main escape fraction
measurements and keep the LyC unaffected by dust for several reasons. It is
not guaranteed, that the stars that emit most of the UV continuum live in the
same environment as the ones that produce the bulk of the LyC. It has been
speculated that even a single star cluster containing a handful of massive,
hot stars such as Wolf-Rayet stars or O-type stars, could be responsible for
most of the ionising emission that we see (e.g. Vanzella et al., 2020, 2022;
Meštrić et al., 2023). While it is true that our spatial resolution does not
allow us to study this assumption, the fact that we detect LyC emission shows
that there must be a clear path with a low neutral hydrogen column density
($\mathrm{N}_{{H\textsc{I}}}\lesssim 10^{17}\,\mathrm{cm}^{-2}$) through which
LyC could escape (see e.g. Reddy et al., 2016b; Steidel et al., 2018; Saldana-
Lopez et al., 2022; Reddy et al., 2022). Such a path could be blown in the ISM
for example via supernovae (SN), however for young starbursts $\sim 3$ Myr,
there might not have been enough time yet for SN to explode. Instead, paths
could be cleared by turbulence or radiative/stellar feedback (e.g. Kakiichi &
Gronke, 2021). In addition, a highly energetic ionising background would also
destroy dust, leaving the dust density low in such an escape channel. Since
dust and neutral hydrogen often coexist, we only need to take dust into
account in more neutral regions. However, the cross-section of dust is four
orders of magnitude smaller than that for photoionisation (e.g. Kakiichi &
Gronke, 2021), meaning the optical depth of dust is always smaller than that
of neutral hydrogen for LyC photons (see Reddy et al., 2016a; Ma et al., 2020;
Rosdahl et al., 2022), which is why the influence of dust on the LyC, even in
neutral regions, is negligible. An exception could be photoionised channels in
which the extinction by dust can become comparable to the absorption by
photoionisation for higher metallicities (significantly above
Z=$0.2\,Z_{\odot}$), however since LyC leakers at high redshifts are expected
to be low metallicity systems and indeed most of our candidates have
metallicities Z$\sim 0.2\,Z_{\odot}$, we argue that this effect is not
significant (Kakiichi & Gronke, 2021). It should also be kept in mind that
using a unity dust covering fraction might not be ideal, since the covering
fraction of neutral gas is not unity and the same processes that clear
channels in the ISM would also clear out the dust, given that the neutral
hydrogen and dust are dynamically coupled.

Several theoretical studies have also shown that indeed the LyC escape
fraction is not strongly influenced by the presence of dust. Kimm et al.
(2019) show on the scale of individual star-forming clouds of
$10^{6}\mathrm{M}_{\odot}$ that the absorption by neutral hydrogen is the
dominant factor, while Mauerhofer et al. (2021) confirm this for a simulated
galaxy of $2.3\times 10^{9}\,\mathrm{M}_{\odot}$. Yoo et al. (2020) study disc
galaxies with a stellar mass of $\sim 10^{9}\mathrm{M}_{\odot}$, which is the
same mass range as our galaxies (see table 4 in the appendix). They find that
dust changes the LyC escape fraction by only $\sim 10^{-3}-10^{-4}\%$, as they
assume that only $1\%$ of dust can survive in ionised regions. Without this
assumption, the escape fraction is reduced by $17\%$ ($37\%$ for solar
metallicities). However, Ma et al. (2020) find that the decrease in escape
fraction at higher masses of $\sim 10^{9}\,\mathrm{M}_{\odot}$ could be due to
dust absorption. For those studies, the presence of dust at high temperatures
plays an important role, as Ma et al. (2020) assume dust up to a temperature
of $10^{6}$ K, while for other studies like Mauerhofer et al. (2021) and Kimm
et al. (2019), dust can only survive up to several $10^{4}$ K.

This is why, for our main escape fraction measurements, we do not include any
dust extinction of the LyC in the SED fitting by CIGALE, which is why the
stellar unattenuated and attenuated models in Figs. 5 and 6 are the same for
wavelengths below $912\,\mbox{\AA}$. We are therefore treating the LyC and the
rest of the spectrum differently because we know that for wavelengths
$\lambda>912\,\mbox{\AA}$ we see the whole stellar population of the galaxy,
where dust has an influence, whereas, for $\lambda<912\,\mbox{\AA}$, we only
see emission that was unaffected by dust for the reasons explained above. This
results in a discontinuity in the fitted spectra in Figs. 5 and 6 at
$912\,\mbox{\AA}$, which however does not affect the measured properties. The
escape fraction from CIGALE is, therefore, a highly sightline-dependent
measure, since we assume that the ionising photons escape through ionised
channels which are low in dust. This is in contrast to the other measures of
escape fraction, assuming a fixed ionising to non-ionising ratio for the whole
galaxy. We can use the measured UV dust extinction from the CIGALE models to
correct the relative escape fractions measured from the two other methods to
absolute escape fractions (shown in table 2 and described above in Sect. 4.3),
using fixed values for the intrinsic ratio, shown in the last column of table
2. In order to estimate the potential influence of dust on our escape fraction
measurements from CIGALE, we provide a second CIGALE measurement in table 2.
Here we extend the dust treatment of CIGALE below $912\AA$, but reduce it by
$90\%$, thus accounting for the fact that since we do detect LyC, we can
assume a relatively low dust extinction in the direction of escape, in
addition to the arguments above. As expected, the LyC escape fraction
increases when dust is taken into account, by factors between $2\%$ and
$92\%$. Keeping in mind the uncertainties concerning the influence of dust on
the LyC, we can assume that the true escape fraction lies between the relative
and the absolute ones given in table 2. In this sense, the absolute escape
fraction without dust absorption in the LyC can be seen as a lower limit.

###  4.7 Escape fraction results

Having now measured the different LyC escape fractions with different methods,
we see that, depending on the assumptions made, escape fractions for the same
object can vary widely. For example, the first gold candidate has a relative
and even absolute escape fraction above $100\%$ if we assume an intrinsic UV
continuum to LyC ratio of three. Using the less steep intrinsic ratio of 1.72
from the BPASS model (see Sect. 4.3), we find more reasonable escape fractions
(below $100\%$), while the CIGALE value is again high at $\approx 90\%$.
Indeed, the escape fractions based on the BPASS model usually result in the
smallest derived escape fractions.

##  5 Lyman $\alpha$ properties of LyC leaker candidates

We can now look at the Ly$\alpha$ properties of the LyC leaker candidates and
compare them to the general population of LAEs. We consider Ly$\alpha$
properties such as FWHM, peak separation and equivalent width, which we take
from Kerutt et al. (2022).

Figure 5: Example of an SED fit using CIGALE (modified CIGALE plot using the
output files). The x-axis shows the logarithmic rest-frame wavelength in
Angstrom, the y-axis shows the logarithmic flux in micro Jansky. The purple
dashed line shows the stellar unattenuated emission and the yellow line shows
the stellar emission taking dust into account. The light blue line shows
nebular emission and the solid black line shows the composite CIGALE model
fit. The empty green dots are the measured values in the individual filter
bands (shown as transparent areas in different colours from purple to yellow)
and the dark green dots show the model values in the same filter bands. The
vertical dashed line indicates the LyC break at $912\mbox{\AA}$. The object
ID, redshift, and reduced $\chi^{2}$ are written on the top.

Figure 6: Same as Fig. 5 for the rest of the gold sample LyC leaker
candidates.

It has been suggested in the literature that Ly$\alpha$ emission could be a
good indicator of LyC leakage (see e.g. Micheva et al., 2017; Verhamme et al.,
2017; Marchi et al., 2018; Steidel et al., 2018; Fletcher et al., 2019;
Vanzella et al., 2020; Pahl et al., 2021; Reddy et al., 2022). Since LyC
emission is needed to ionise the neutral hydrogen in the ISM of a galaxy,
which recombines to emit a Ly$\alpha$ photon in $\approx$ 2/3 of the cases
(e.g. Osterbrock, 1989; Dijkstra, 2014), the two types of emission are closely
linked. The more ionising radiation is absorbed by the neutral hydrogen in the
ISM, the stronger the intrinsic Ly$\alpha$ emission. In the same vein, the
escape channels of Ly$\alpha$ and LyC could be similar, as Ly$\alpha$ will be
scattered and potentially absorbed by dust in a high neutral hydrogen column
density environment, while the LyC would be absorbed by the neutral hydrogen
as well. A low neutral hydrogen column density would thus permit an easier
escape for both Ly$\alpha$ and the LyC. Therefore, a high Ly$\alpha$
equivalent width, as well as a narrow line and a narrow peak separation, have
been proposed as good indicators of potential LyC leakage (e.g. Verhamme et
al., 2017; Vanzella et al., 2020; Izotov et al., 2021, see Sect. 6 for a
discussion).

###  5.1 Lyman $\alpha$ Properties of the LyC leaker candidates

In Table 3 we provide the rest-frame equivalent width, luminosity, flux, FWHM,
and the peak separation (for objects with a double peak) of Ly$\alpha$ as well
as absolute UV magnitudes for our candidates, all values taken from Kerutt et
al. (2022). It is worth mentioning that none of our candidates has a
Ly$\alpha$ $\mathrm{EW}_{0}>240\,\mbox{\AA}$ ($11\%$ of our final sample have
such high $\mathrm{EW}_{0}$), which is usually cited as the approximate upper
limit for normal stellar populations. This is interesting, as it was expected
that such unusual stellar populations (with low metallicities, potentially
containing population III stars, high star formation rates, and young ages)
would be more likely to produce noticeable amounts of LyC emission (see Sect.
6 for a discussion). Likewise, the FWHM and peak separations are both above
the average of the full sample (see Kerutt et al., 2022), contrary to
theoretical expectations which expect small FWHM and peak separations
(Verhamme et al., 2017). This can be seen in Fig. 7, where we show the
distributions of FWHM and peak separation against the Ly$\alpha$ luminosity of
the whole sample of LAEs, highlighting the respective values of our LyC leaker
candidates. Both the FWHM and the peak separation values of the candidates are
above the mean of the full sample, which is $\mathrm{FWHM}=218\pm
102\,\mathrm{km}\,\mathrm{s}^{-1}$ and peak sep. $=481\pm
244\,\mathrm{km}\,\mathrm{s}^{-1}$ (Kerutt et al., 2022), compared to
$\mathrm{FWHM}=270\pm 67\,\mathrm{km}\,\mathrm{s}^{-1}$ and peak sep. $=602\pm
69\,\mathrm{km}\,\mathrm{s}^{-1}$ for our sample of LyC leaker candidates. We
show the MUSE spectra of our candidates in the appendix in Fig. 12. As can be
seen from those spectra, some of the double peaks have a low S/N and are
rather uncertain. However, assuming the correlation of Verhamme et al. (2017)
between FWHM and velocity offset, as well as peak separation and velocity
offset, there is also a correlation between FWHM and peak separation, where
the peak separation is expected to be roughly twice as large as the FWHM,
which is the case for our double peaks, except for ID 1521589, which has a
rather high peak separation, but a convincing blue bump. Therefore we can
assume that our results are not influenced by a falsely identified blue peak.

Table 3: Overview of candidates - Ly$\alpha$ properties.

IDMW | EW0 | MUV | $\mathrm{log_{10}(L}_{\mathrm{Ly}\alpha})$ |  $10^{-18}$ fLyα | FWHM | peak sep.
---|---|---|---|---|---|---
| [Å] | [magAB] | $\mathrm{[erg/s]}$ | [erg/s/cm${}^{2}/\mbox{\AA}$] | [km/s] | [km/s]
gold candidates |  |  |  |  |  |
1181371 |  1$60\pm 13$ | $-18.25\pm 0.17$ | $42.0\pm 0.4$ | $12\pm 2$ | $277\pm 43$ |
3052076 |  1$18\pm 1$ | $-20.89\pm 0.04$ | $42.6\pm 0.1$ | $33\pm 1$ | $304\pm 33$ | $677\pm 32$
109004028 |  1$33\pm 7$ | $-19.65\pm 0.05$ | $42.3\pm 0.5$ | $22\pm 4$ | $318\pm 179$ |
122021111 |  1$45\pm 2$ | $-21.05\pm 0.02$ | $43.0\pm 0.1$ | $77\pm 3$ | $255\pm 13$ | $595\pm 11$
126049137 | $174\pm 28$ | $-19.61\pm 0.07$ | $43.0\pm 0.3$ | $54\pm 8$ | $342\pm 107$ | $661\pm 121$
silver candidates |  |  |  |  |  |
1521589 | $161\pm 14$ | $-17.64\pm 0.02$ | $42.2\pm 0.2$ | $18\pm 1$ | $188\pm 30$ | $556\pm 62$
3452147 | $104\pm 26$ | $-18.49\pm 0.24$ | $42.4\pm 0.3$ | $20\pm 2$ | $282\pm 45$ | $565\pm 30$
4062373 |  1$18\pm 2$ | $-20.01\pm 0.04$ | $42.2\pm 0.3$ | $13\pm 1$ | $264\pm 78$ |
4172404 |  1$94\pm 17$ | $-18.56\pm 0.18$ | $42.3\pm 0.2$ | $17\pm 1$ | $249\pm 22$ |
5622786 | $133\pm 33$ | $-17.80\pm 0.24$ | $42.2\pm 0.3$ | $10\pm 1$ | $203\pm 29$ | $459\pm 103$
119004004 | $118\pm 15$ | $-19.07\pm 0.11$ | $42.6\pm 0.2$ | $44\pm 3$ | $318\pm 37$ | $637\pm 87$
122032127 |  1$26\pm 8$ | $-19.20\pm 0.16$ | $42.0\pm 0.7$ |  1$5\pm 1$ | $123\pm 62$ |
mean |  1$78\pm 53$ | $-19.33\pm 1.14$ | $42.4\pm 0.3$ | $28\pm 20$ | $270\pm 67$ | $602\pm 69$

777IDMW: Identifier in Urrutia et al. (2019). EW0: Ly$\alpha$ rest-frame
equivalent width from MUSE-Wide in Å. MUV: Absolute UV continuum AB magnitude
at $1500\,\mbox{\AA}$, measured from HST bands.
$\mathrm{log_{10}(L}_{\mathrm{Ly}\alpha})$: Logarithmic Ly$\alpha$ luminosity
in $\mathrm{erg/s}$. $10^{-18}$ fLyα: Ly$\alpha$ flux in $10^{-18}$
erg/s/cm${}^{2}/\,\mbox{\AA}$. Ly$\alpha$ flux and luminosity were measured
from MUSE-Wide data within 3 Kron radii and include the blue bump (for double-
peaked lines). FWHM: Full width at half maximum of the main (red) peak of the
Ly$\alpha$ line in km/s. Peak sep.: Peak separation for lines with double
peaks in km/s. Both the FWHM and the peak separation are corrected for the
line spread function of MUSE. The Ly$\alpha$ measurements are taken from
Kerutt et al. (2022). The values in the last line give the mean values for all
candidates with the standard deviation of the measurements.

Figure 7: Left: Peak separation of the Ly$\alpha$ line as a function of the
logarithmic Ly$\alpha$ luminosity. The contours (each showing a $10\%$ number
difference) contain the full sample of LAEs from Kerutt et al. (2022) and show
the density of objects in the peak separation and Ly$\alpha$ luminosity plane.
The gold and silver dots show the individual values for the LyC leaker
candidates, with the sizes indicating the escape fraction values based on the
SED models by CIGALE. The gold candidate with a green circle is the object
already discovered in Saxena et al. (2022) and discussed in Rivera-Thorsen et
al. (2022). Objects without double peaks are placed at peak sep. = 0 at their
respective luminosities. The blue dots show data of low redshift LyC leakers
from Gazagnes et al. (2020), also featured in Izotov et al. (2016a, b, 2018a,
2018b) and Maji et al. (2022). Right: FWHM of the Ly$\alpha$ line as a
function of the logarithmic Ly$\alpha$ luminosity, again with contours showing
the distribution of objects in these values.

###  5.2 Connecting Lyman $\alpha$ Properties to Leakage

The Ly$\alpha$ properties of the LyC leakers are somewhat surprising, as it
would have been expected from theoretical models and previous observations of
low-redshift analogues (see e.g. Verhamme et al., 2017; Izotov et al., 2018b)
that a small peak separation and narrow FWHM, as well as a high equivalent
width, should be good indicators of LyC leakage. In the left panel of Fig. 8
we show the measured escape fraction from the SED models from CIGALE over the
peak separation of the Ly$\alpha$ line, with a comparison to literature values
from Izotov et al. (2018b) and Verhamme et al. (2017) for low-redshift LyC
leakers. To compare to LyC leakers at a similar redshift as ours, we use the
velocity offsets given in Fletcher et al. (2019) for LyC leakers at $z\approx
3.1$ and multiply them by two, assuming that there is a correlation between
the velocity offset and the peak separation (as shown in Verhamme et al.,
2018).

It has to be kept in mind that here we show the sightline-dependent escape
fraction measurement from CIGALE, where we assume no dust absorption of the
Lyman continuum since they escape along ionised channels with high gas
temperatures. The escape fraction estimates we get are not the global LyC
escape fractions of the galaxies. The same is true for Ly$\alpha$, which can
vary substantially with viewing angle (see Blaizot et al., 2023).

The results of these comparisons are further discussed below in Sect. 6 below.

Figure 8: Left: LyC escape fraction as a function of Ly$\alpha$ peak
separation. The silver and gold-filled dots are from our LyC leaker candidate
sample. The gold candidate with a green circle around it is again the object
already discovered in Saxena et al. (2022) and discussed in Rivera-Thorsen et
al. (2022). The black dash-dotted line is the relation in Izotov et al.
(2018b) for low-redshift LyC leakers. Open circles show low redshift objects
and filled dots show high redshift objects. The dark green data points are for
LyC leakers at $z\approx 3.1$ from Fletcher et al. (2019), where they give
velocity offsets with respect to systemic redshift, which was multiplied here
by two to estimate the potential peak separation, following the correlation
found between peak separation and the shift of the red peak with respect to
systemic velocity from Verhamme et al. (2018). The dark purple filled circle
is from Marques-Chaves et al. (2022) for $z=3.6$. The light green open circles
are taken from Verhamme et al. (2017), also for low-redshift analogues. The
black dots are taken from Flury et al. (2022b), showing the escape fractions
based on the COS UV spectra for only their strongest LyC leakers at low
redshifts ($z\approx 0.3-0.4$). Right: LyC escape fraction as a function of
the Ly$\alpha$ rest-frame equivalent width. The black dashed line is from
Steidel et al. (2018), just as the light green dots for objects at $z\approx
3$. The dark green dots are from Fletcher et al. (2019) for objects at
$z\approx 3$. The open black dots are again taken from Flury et al. (2022b)
and the open blue dots are from Pahl et al. (2021) at $z\approx 3$.

##  6 Discussion

Due to the need to infer LyC emission at the epoch of reionisation from other,
measurable properties of high-redshift galaxies, we have looked at the
Ly$\alpha$ line as a possible indicator. In this section, we discuss its
potential and also suggest a way to use the number of detected LyC leakers as
a possible way to predict the global escape fraction.

###  6.1 The Dispersion of LyC escape fractions

In table 2 we have seen that the escape fraction is sensitive to the method
and assumptions applied. This makes a comparison with literature results
difficult, as different studies use different values for the intrinsic ratio
(L1500/L${}_{900})_{\mathrm{int}}$ between the UV continuum and the LyC,
usually in the range of 3, 5 or 7, as well as for the IGM transmission. In
addition, dust extinction is also treated in different ways in the computation
of the absolute escape fraction. What is more, the limiting magnitudes of
different studies could result in biases, if the LyC escape is correlated with
UV luminosity and/or stellar masses. Keeping these caveats in mind, we
nevertheless compile a collection of literature values for the absolute escape
fraction as a function of redshift in Fig. 9 based in part on a collection of
LyC leakers from the literature in the appendix (5, 6 and 7). Simulations
predict an increase in the escape fraction with higher redshifts between
$z=4-8$ (e.g. Trebitsch et al., 2021). There are two more densely populated
areas in this plot, one is the region at $z<0.5$, containing low-redshift LyC
leakers as analogues of high-redshift ones. The other region is roughly
between $z=2-4$, where LyC emission is observable with the HST using the
filters WFC3/UVIS F275W and F336W, as done in this study. The region between
these intermediate redshift objects and the low-redshift analogues is scarcely
populated, mostly by stacking studies. In the range $z=3-4.5$, we show the
escape fraction measurements of our own LyC leaker candidates in gold and
silver (both for a fixed intrinsic flux ratio of 3 and based on the CIGALE SED
fitting), as well as upper limits for our non-detections, based on the CIGALE
SEDs. We use $2\sigma$ as the upper flux limit (since this is also our cut for
a detection) and the mean IGM transmission at the redshift bin of each object
(in steps of $\Delta z=0.1$), as well as the median dust extinction
$\mathrm{A}_{\mathrm{V}}$ of the 12 LyC leaker candidates to convert to
absolute escape fractions. Since the escape fraction measurements can vary
widely based on the assumptions and parameters used for the computation, but
mostly because of the different sample selections, potential biases, and
completeness issues, we cannot draw conclusions about a potential trend with
redshift from this plot and it should be seen as a compilation of literature
values.

Figure 9: Logarithmic absolute escape fraction of LyC over redshift for a
collection of literature values and our sample of gold and silver objects, in
their respective colours. Here, circles show measurements from SED fitting
using CIGALE, while squares are escape fractions using an intrinsic ratio of 3
to compare to the literature. The small blue dots represent individual
detections of LyC emission, for low-redshift analogues and high redshifts. The
small purple dots are from the LzLCS survey (Flury et al., 2022a, b). The
large green symbols indicate results from stacking analyses. Stacks, where
only upper or lower limits could be estimated, are shown with downwards-facing
and upwards-facing triangles, respectively. The large transparent circles are
the mean values of the respective sub-sample, with the grey containing both
the gold and silver samples. The light grey dots between $z=3-4.5$ are the
upper limits of our sample (based on an intrinsic ratio of 3), and the grey
downward triangle is the median of those upper limits. The magenta dot
corresponds to the inferred average escape fraction from the simulations
described in Sect. 6.3. The data shown is based on tables 5, 6 and 7 in the
appendix.

Most of the known low-redshift LyC leakers were found in the LzLCS (Flury et
al., 2022a, b; Saldana-Lopez et al., 2022; Chisholm et al., 2022) with a range
of escape fractions between 0 and $\approx 50\%$. In their studies, Flury et
al. (2022b) test among other things the assumed relation between Ly$\alpha$
properties and LyC leakage and find that the Ly$\alpha$ escape fraction,
equivalent width and peak separation correlates best with LyC escape
fractions. This connection is expected, as explained in the introduction and
Sect. 5. An optical depth of hydrogen ionising photons of one corresponds to a
neutral hydrogen column density of $\mathrm{N}_{\mathrm{HI}}=1.6\times
10^{17}\,\mathrm{cm}^{-2}$ (explained e.g. in Verhamme et al., 2015), however,
due to the larger cross-section of Ly$\alpha$, neutral hydrogen is optically
thick to Ly$\alpha$ down to $\mathrm{N}_{\mathrm{HI}}\sim
10^{13}\,\mathrm{cm}^{-2}$. In this density range the LyC emission can escape
more easily, while Ly$\alpha$ will be scattered, which results in a double
peak with increasing peak separation for increasing
$\mathrm{N}_{\mathrm{HI}}$.

From a theoretical point of view, it would therefore be expected that a narrow
Ly$\alpha$ line ($\mathrm{FWHM}\approx 200\,\mathrm{km/s}$) or a double peak
with a small peak separation ($\mathrm{peak\,sep.}<300\,\mathrm{km/s}$,
Verhamme et al., 2015), would indicate a low neutral hydrogen column density,
corresponding to the limit where the optical depth of LyC is around one.
However, if we look at Fig. 7, we see that our LyC leaker candidates are
typically not narrow and do not have narrow peak separations (our smallest
peak separation being $459\pm 103\,\mathrm{km/s}$). We also do not find a
connection between the peak separation and the measured LyC escape fraction in
the left panel of Fig. 8, although this has been shown to be a reliable
indicator for low-redshift analogues (Verhamme et al., 2017; Izotov et al.,
2018b; Flury et al., 2022b). For higher-redshift LyC leakers as in Fletcher et
al. (2019) at $z\approx 3.1$, the velocity offset with respect to systemic
also does not seem to be a good enough indicator for LyC escape. Similarly, a
high Ly$\alpha$ equivalent width (shown in the right panel of Fig. 7) is not a
requirement for a high LyC escape fraction. Still, Steidel et al. (2018) find
a correlation for restframe Ly$\alpha$ equivalent widths of up to $\approx
50\,\mbox{\AA}$, Marchi et al. (2018) find a LyC escape fraction of $33\%$ for
their stack of high Ly$\alpha$ EWs, and Begley et al. (2022) find a higher
escape fraction as well for their stack of the half of their sample with
higher Ly$\alpha$ equivalent widths (although it has to be noted that this
half of the sample only has a median $\mathrm{EW}_{0}=4.9\,\mbox{\AA}$). In
contrast, we find several of our LyC leaker candidates have high escape
fractions even with smaller equivalent widths. However, the highest Ly$\alpha$
equivalent width objects among our candidates have indeed high escape
fractions.

Another possible connection between LyC and Ly$\alpha$ has recently been
proposed by Maji et al. (2022), who analyse simulated galaxies at high
redshifts and find a strong correlation between the luminosities of the two
due to their production being related to the star formation rate of the
galaxy, as massive stars produce LyC photons, which in turn create Ly$\alpha$
photons. Here, we look at a possible connection between the LyC escape
fraction and the Ly$\alpha$ luminosity in Fig. 10. We do not detect any trend,
although there seems to be one for low-redshift leakers (see the figure
caption and also the LzLCS, Flury et al., 2022b).

Figure 10: LyC escape fraction versus the logarithmic Ly$\alpha$ luminosity
for the gold and silver sample of our LyC leaker candidates. In green, we show
the selection of low-redshift LyC leakers shown in Maji et al. (2022),
consisting of objects from Leitet et al. (2011, 2013); Borthakur et al.
(2014); Pardy et al. (2016); Izotov et al. (2016a, b); Verhamme et al. (2017);
Puschnig et al. (2017); Izotov et al. (2018a, b); Micheva et al. (2020);
Izotov et al. (2021); Gazagnes et al. (2021). The blue dots show low-redshift
leakers from Izotov et al. (2022).

###  6.2 Potential limitations of using Ly$\alpha$

The question remains why these correlations do not seem to work well for LyC
leakers at higher redshifts while they are sufficiently reliable at lower
redshifts (e.g. Verhamme et al., 2017; Flury et al., 2022b). One caveat is
that although there has been a lot of progress in the field in recent years,
the number of confirmed and robust LyC leakers at higher redshifts is still
small and we often need to rely on stacking analyses. An advantage of stacking
is that we can statistically correct for the IGM absorption, which is one of
the largest uncertainties for studies of individual LyC leakers.

The assumed IGM absorption is only one aspect of determining the LyC escape
fraction, however. When comparing escape fractions from different papers, like
in Figs. 8, 9 and 10, it has to be kept in mind that the assumptions that went
into the escape fraction can vary widely, not only when it comes to the IGM
absorption, but also concerning the intrinsic ratio between UV continuum and
LyC, as well as the age, metallicity and dust attenuation. Furthermore, the
samples might not be representative of the global population of LyC leakers,
as often only the brightest ones can be detected. Another caveat is
differences in resolution limits, which make it difficult to compare trends
with peak separation. In the MUSE data, we can go down to a peak separation of
$\sim 100\,\mathrm{km/s}$, however, the smaller the peak separation, the
harder it is to distinguish a double peak from a single peak. Keeping these
difficulties in mind, we can now think about physical reasons why the trends
with escape fraction observed at low redshifts do not seem to hold for our
sample.

The galaxies at redshifts around $z\approx 3-4.5$ are likely to have different
properties than the ones at low redshift. The discrepancy between our study
and the results from low-redshift LyC leakers might indicate that we are not
comparing the same kinds of objects after all. For the current study, we use a
catalogue of LAEs that was produced from a blind emission-line search without
pre-selection for any galaxy properties, unlike studies like for example the
LzLCs (Wang et al., 2021; Flury et al., 2022a, b; Saldana-Lopez et al., 2022),
which are aimed specifically at detecting LyC. One known difference, for
example, is that the morphological properties of high-redshift LAEs are more
irregular, which could mean that the kinematics and ISM properties are more
complex (e.g. Guaita et al., 2015) than the often compact low-redshift
leakers. A prominent example of a local LyC leaker with a complex morphology
is Haro 11 (e.g. Bergvall et al., 2006; Keenan et al., 2017; Rivera-Thorsen et
al., 2017b), which has three distinct clumps, two of which show no LyC escape.
However, all three show strong Ly$\alpha$ emission, with the highest
Ly$\alpha$ escape coming from the clump which is likely to be the one leaking
the LyC emission as well. There have been indications that this object also
has a complex velocity field, indicating a possible merger (e.g. Östlin et
al., 2015). This shows that without detailed knowledge of the morphology and
the origin of the Ly$\alpha$ and LyC emission, it is not easy to draw a
connection between the two. At higher redshift, Ly$\alpha$ emission does not
always come from one single object, but there could be clumps or two objects
close by. Sometimes this can even create the impression of a double peak,
although both lines come from different clumps (see Vitte et al. in prep,
using data from MUSE as well). For objects that do not show a clear connection
between Ly$\alpha$ and LyC, this could mean that the two are escaping from
different locations in the galaxy. Indeed, when looking at the positions of
Ly$\alpha$ and the LyC contours in Figs. 3 and 4, there seems to be an offset
for some of the gold candidates (roughly around
$0\aas@@fstack{\prime\prime}3-0\aas@@fstack{\prime\prime}5$ between the
respective SN peaks) and the LyC emission is typically not centred on the
Ly$\alpha$ emission. Another indicator for this is, that even for the objects
where we do not find a double-peaked Ly$\alpha$ line, the lines nevertheless
look asymmetric, indicating scattering or expanding gas, which disfavours an
ionised channel and again hints at different origins of the LyC and
Ly$\alpha$. It has to be kept in mind, that we do not see the overall LyC
escape fraction of a galaxy, but only the escape in our line of sight. It is
possible that only a few sightlines have visible LyC emission, which would
also explain why LyC leakers are so rare. Ly$\alpha$ emission however is
scattered and there can be a combination of different sources which form the
Ly$\alpha$ line, as for example Ji et al. (2020) argue that the fact they see
Ly$\alpha$ in absorption in the LyC leaker Ion 1 could be caused by a
combination of emission originating from the same place as the LyC photons and
absorption.

If however there is not enough neutral hydrogen to produce Ly$\alpha$ in the
first place, this could also explain low Ly$\alpha$ equivalent widths for high
LyC escape fractions. The reverse, a high Ly$\alpha$ equivalent width but no
LyC leakage, could occur in the case of a higher neutral hydrogen density in a
clumpy medium, allowing Ly$\alpha$ to scatter to escape the galaxy, while LyC
cannot and will be absorbed in the clumps (see Neufeld, 1991, but brought into
question by Gronke & Dijkstra, 2014; Duval et al., 2014). Another reason might
be that the LyC emission could already have been transformed into Ly$\alpha$
in dense neutral hydrogen regions. Indeed, all Ly$\alpha$ emission we see is
LyC emission that has been lost already, usually in the first
$10\,\mathrm{pc}$ around the stellar population where it was produced (e.g.
Paardekooper et al., 2015).

A high LyC escape fraction in objects with small Ly$\alpha$ equivalent width
can also be explained without morphological arguments. It is often assumed
that extreme stellar populations are needed for a high production of LyC and a
high Ly$\alpha$ equivalent width ($>240\,\mbox{\AA}$). However, it has to be
kept in mind that the LyC escape depends just as much on the foreground
opacity as on the stellar population properties. Indeed, Reddy et al. (2022)
argue that changes in the stellar population alone cannot explain high escape
fractions and only lead to moderate changes in the intrinsic LyC photon
production. Instead, a high LyC escape fraction with small Ly$\alpha$
equivalent width could indicate a neutral hydrogen column density low enough
($\mathrm{N}_{\mathrm{HI}}<10^{17}\,\mathrm{cm}^{-2}$) for LyC to escape more
easily than Ly$\alpha$, which could be scattered out of the line of sight. A
prominent example of a LyC leaker which even has Ly$\alpha$ in absorption is
Ion 1 (e.g. Vanzella et al., 2010b; Ji et al., 2020). Ji et al. (2020) argue
that since the optical thickness of Ly$\alpha$ photons is $\approx 10^{4}$
larger than that of LyC photons, the latter can escape more easily depending
on the neutral hydrogen column density.

Therefore, there might be holes with low neutral hydrogen column densities in
the ISM through which LyC could escape, even if most of the produced
Ly$\alpha$ emission is scattered in the rest of the galaxy. LyC could have
carved out ionised paths (as shown e.g. in Mainali et al., 2022 for the
sunburst arc) or they could be created by supernova-induced outflows and holes
(e.g. Clarke & Oey, 2002; Fujita et al., 2003). Ji et al. (2020) even suggest
the presence of escape channels with no neutral hydrogen, which would hinder
the production of Ly$\alpha$ emission. Another possible example of this is the
low-redshift LyC leaker Tololo 1247-232 (e.g. Leitet et al., 2013; Leitherer
et al., 2016), which seems to consist of multiple stellar populations with
different ages. The older population of $12\,\mathrm{Myr}$ could have cleared
a path for the LyC emission to escape, while the younger component of
$2\\--4\,\mathrm{Myr}$ could provide most of the ionising emission (Micheva et
al., 2018).

Another promising indicator of LyC escape is therefore the Ly$\alpha$ emission
at the systemic redshift (measured through non-resonant lines), as shown by
Naidu et al. (2022). Unfortunately, we do not have other spectral lines for
most of our sample to determine systemic redshifts. Another way in which
emission at systemic would manifest itself, even if some or most of the
Ly$\alpha$ emission was scattered in a high neutral hydrogen column density,
is by showing an additional peak in the middle between the other two, thus
creating a triple-peaked line. However, the resolution and depth of the
spectroscopic data for our current sample of LAEs are not sufficient to find
triple peaks and also result in rather large error bars on Ly$\alpha$
properties.

Another aspect that should be mentioned, is the possibility of interlopers
contaminating our sample. As discussed in Sect. 3.3, we take several measures
to avoid low-redshift interlopers and contamination from close neighbours, but
there is a remaining possibility that some of our objects could be
contaminated. Our candidate with the highest LyC over UV continuum ratio (ID
1181371) for example has a neighbour in the filter bands ACS F606W and F775W
(see Fig. 13), that could be the source of the Ly$\alpha$ emission. The
potential contaminant is at a distance of $\approx
0\aas@@fstack{\prime\prime}5$, though, which makes us confident that it is not
connected to the Ly$\alpha$ emission. Having otherwise excluded objects with
close neighbours, objects where the rgb images indicate different colours for
different parts of the objects (hinting at a chance alignment), objects that
have a too strong signal in the shorter wavelength band WFC3/UVIS F275W, as
well as objects where the potential LyC emission and the UV emission do not
overlap significantly, we are confident that most if not all of our candidates
are robust.

###  6.3 Expected LyC leaker fraction

Instead of deriving the LyC escape fraction from the individual objects, we
now want to turn the question around and derive the expected number of LyC
leakers that we could potentially find in our survey area, given a certain
assumed global escape fraction. This way, we can compare the actual number of
LyC leaker candidates to predictions for different escape fractions and thus
derive a global escape fraction.

For this approach, we first construct a LyC luminosity function. The basis for
this is a UV luminosity function of the typical Schechter (1976) shape:

$\phi(L)\,\mathrm{d}L=\phi^{*}\,\left(\frac{L}{L^{*}}\right)^{\alpha}\,\mathrm{e}^{-L/L^{*}}\,\frac{\mathrm{d}L}{L^{*}}$
(3)

With $\phi^{*}$ being the normalisation in comoving Mpc-3, $L^{*}$ the
characteristic luminosity in $\mathrm{erg}\,\mathrm{s}^{-1}$ and $\alpha$ the
faint-end slope. Here we use the Schechter function in terms of absolute
magnitude $\phi(L)\,\mathrm{d}L=\phi(M_{\mathrm{UV}})dM_{\mathrm{UV}}$, which
can be rewritten using
$\frac{L}{L^{*}}=10^{-0.4(M_{\mathrm{UV}}-M_{\mathrm{UV}}^{*})}$ as:

$\phi(M_{\mathrm{UV}})=\frac{\mathrm{ln}10}{2.5}\phi^{*}10^{-0.4(M_{\mathrm{UV}}-M_{\mathrm{UV}}^{*})(\alpha+1)}\exp\left(-10^{-0.4(M_{\mathrm{UV}}-M_{\mathrm{UV}}^{*})}\right)$
(4)

where $M_{\mathrm{UV}}^{*}$ is the characteristic absolute magnitude. For the
UV LF (thick lines in Fig. 11) we use values from Bouwens et al. (2021) in
steps of $\Delta z=0.1$ (starting at $z=3.05$) to make a redshift grid (shown
in different colours in Fig. 11).

To convert the UV LF to the LyC LF, we assume an intrinsic ratio of UV
continuum to LyC luminosities of $L_{\mathrm{UV}}/L_{\mathrm{LyC}}=3$ (in
frequency space), see the thin lines in Fig. 11.

Including the effects of IGM absorption using the transmission lines from
Inoue et al. (2014) shifts the luminosity function further down (see dashed
lines in Fig. 11). With this function, we can predict the number of observable
LyC leakers for a given survey size and escape fraction at different redshifts
by integrating the luminosity function (thick pink dashed line and pink dot in
Fig. 11).

Taking the information that we find 5 LyC leaker candidates up to a redshift
of $z=4.5$, for a survey area of 43 armin2, and a detection limit for a
$3\sigma$ detection in the HST band WFC3/UVIS F336W of
$28.75\,\mathrm{mag}_{\mathrm{AB}}$, we adjust the assumed escape fraction
until we match the number of detected LyC leakers. We find that for a global
escape fraction of $12\%$, we can reproduce the number of $\sim 4.9$ LyC
leaker candidates, which is close enough to our five candidates. From this
simple model and assuming the intrinsic UV to LyC continuum flux ratio is
accurate on average, we can thus conclude that the average escape fraction of
our sample is at least $\approx 12\%$ (assuming a distribution of line of
sight IGM transmissions based on Inoue et al., 2014). Incidentally, this
matches well with the $1\sigma$ limit found in a stack of LyC non-detections
of LBGs by Saxena et al. (2022).

For our silver sample, we set the detection limit to $2\sigma$ and our
faintest object has an AB magnitude of $\mathrm{m}_{\mathrm{AB}}=29.2$. If we
use this magnitude as a limit in Fig. 11, we get a number of $11.9$ expected
LyC leakers, meaning seven additional ones, which again matches well with our
sample of seven silver candidates. In a similar vein, Begley et al. (2022)
predict the expected number of LyC leaker candidates among their sample from
the average escape fraction they find in their stack star-forming galaxies at
$z\simeq 3.5$, which matches well with their individual detections for
$f_{\mathrm{esc}}=0.07\pm 0.02$. One caveat to mention is the influence of
cosmic variance on our results. Since we are looking at a relatively small
field of view and low number statistics of LyC leakers, the resulting escape
fraction from such a calculation could be different for different fields.
However, our calculation shows that it is plausible that galaxies at and after
the EoR could be leaking LyC emission with an average escape fraction of
$12\%$. However, most of them are not visible to us in LyC due to the
absorption in the IGM. Incidentally, recently Mascia et al. (2023) find an
inferred average escape fraction of $12\%$ (based on correlations with the O
III/O II ratio and $\beta$ parameter) in their sample of gravitationally
lensed galaxies at $z=4.5-8$ in the first results from the GLASS-JWST program
(Treu et al., 2022). The distribution of escape fractions for individual
galaxies however is not clear. Some studies favour reionisation by the
majority of faint galaxies (e.g. Finkelstein et al., 2019), while others
suggest that only a small minority of bright star-forming galaxies (such as
the one discovered in Marques-Chaves et al., 2021 or Marques-Chaves et al.,
2022) could be responsible for the entire ionising photon emission from star-
forming galaxies. However, Naidu et al. (2022), the authors claim that half of
the population of LAEs at $z\approx 2$ has an escape fraction of $25-50\%$,
which would mean that not only the brightest ones are responsible for the
ionising emission.

Coming back to Fig. 9, where we show the result from our exercise as the pink
dot, we can see that the estimated $12\%$ escape fraction is below most of the
individual measurements at $z=2-4.5$, but above or similar to values measured
from stacking and for lower redshifts of $z=0-1$. It is also close to the
median value of the upper limits for the non-detections. This suggests that
individual measurements are biased towards higher escape fraction values.

Figure 11: Expected number of LyC leakers (logarithmic, y-axis) versus
observed LyC AB magnitude. The thick, solid lines show the UV LF (the number
of objects expected in the HDUV area at the given redshift), while the thin,
solid lines show the predicted LyC LF, both colour-coded for different
redshifts (see colour bar). The dashed lines indicate the expected number at
different redshifts, using IGM transmission lines from Inoue et al. (2014) and
assuming an escape fraction of $12\%$. The colours show redshifts from 3.0 to
4.5 in redshift bins of $\Delta z=0.1$, which are all added up for the pink
dashed line. The pink dot marks the expected number of LyC leaker candidates
($4.9$) for the HDUV area in GOODS-S (with a size of $43$ arcmin2). This is
given for the depth in the HST band F336W of $\approx
28.75\,\mathrm{mag}_{\mathrm{AB}}$, which corresponds roughly to a $3\sigma$
detection.

##  7 Summary and Conclusions

We have identified five very likely LyC leaker candidates (that we have called
our ”gold” sample), four of which are new, among a sample of LAEs previously
detected by MUSE in the CDFS region, with an additional seven potential
candidates (that we have called our ”silver” sample). We measured their LyC
emission in the WFC3/UVIS F336W filter from the HDUV survey (Oesch et al.,
2015, 2018) and performed our own photometry in several filters and SED
fitting with the CIGALE software (Burgarella et al., 2005; Noll et al., 2009;
Boquien et al., 2019).

We assume that the ionising photons escape preferably through ionised channels
that are low in dust due to the high gas temperature (making our escape
fraction measurements from CIGALE sightline-dependent). In addition, assuming
a high transmission through the IGM (using the highest $5\%$ of the
distribution of Inoue et al., 2014), we find escape fractions ranging from
$22\%$ to $88\%$ in the redshift range $z=3.08$ to $z=4.43$, one of the
highest redshifts for a LyC leaker found so far. We used our knowledge of the
Ly$\alpha$ line properties from the MUSE spectra to check several previously
proposed indications of LyC leakage, such as FWHM, peak separation, equivalent
width and Ly$\alpha$ luminosity. We found no reliable correlation with any of
these values, even though especially the peak separation has proven to be a
good indicator for low-redshift LyC leakers. We argue that this could be
explained by large uncertainties on the measured escape fractions, as well as
the dependence on the underlying assumptions concerning dust extinction in the
ISM, absorption in the IGM and the intrinsic ratio between UV and LyC. It is
also possible that LyC leakers at higher redshifts have different properties
or production mechanisms of LyC photons than lower-redshift analogues.

We show a method to infer the global escape fraction from the number of found
LyC leaker candidates, by integrating the IGM corrected LyC LF down to the
depth of our data, for which we use the UV LF (Bouwens et al., 2021) and
assume an intrinsic UV to LyC ratio. Taking our five LyC leaker candidates
from the gold sample, we thus derive a global escape fraction of
$f_{\mathrm{esc}}=12\%$, which would be enough to reionise the universe. From
this exercise, we also see that in order to detect a significant number of
reliable LyC emitters, we need to go an order of magnitude deeper in the
WFC3/UVIS F336W band.

In this paper, we have focused on the detection of individual LyC leaker
candidates, while also comparing to literature results from stacking (see Fig.
9 and table 6 in the appendix). A possible next step to determine the global
LyC escape fraction from our sample of LAEs is stacking analysis, either by
stacking the WFC3/UVIS F336W cutouts to detect LyC from the photometry (see
Oesch et al. in prep.) or by stacking the spectra of the LAEs from MUSE. This
would reduce the possible redshift range for direct detection since the LyC
enters the MUSE wavelength only at a redshift of $z>4$. However, a stacking
analysis of LAE spectra can still give us information on UV lines that could
indicate both the ionising photon production as well as the ISM properties,
which are both connected to the production and escape of LyC emission (see
Kramarenko et al. in prep.). For example, recent studies have suggested that
the [O III]/[O II] ratio between [O III] line(s) at $5007$ and
$4959\,\mbox{\AA}$ to [O II] at $3727\,\mbox{\AA}$ could be a good indicator
of a strong ionisation field and low optical depths (proposed by Jaskot & Oey,
2013, see also e.g. Nakajima & Ouchi, 2014; Nakajima et al., 2016; Faisst,
2016; Izotov et al., 2018a; Katz et al., 2020 but Stasińska et al., 2015;
Naidu et al., 2018; Bassett et al., 2019 for counter-examples). Another
potential line would be Mg II $\lambda\lambda 2796,2803$ (Henry et al., 2018;
Katz et al., 2022), which correlates strongly with Lyman $\alpha$ emission and
could be used as a substitute if Ly$\alpha$ is not observable. Similarly,
${C\textsc{III}}]\lambda\lambda 1907,1909$ (Jaskot & Ravindranath, 2016) and
${C\textsc{IV}}\lambda 1550$ (Schaerer et al., 2022) trace low metallicities
and density-bounded regions of ionised hydrogen. Not only emission but also
residual flux in saturated low-ionisation interstellar absorption lines (e.g.
Borthakur et al., 2014; Mauerhofer et al., 2021) can indicate the state of the
ISM and thus potential LyC emission.

While LyC emission remains impossible to detect at the epoch of reionisation,
we are gaining new insights into this era through various JWST programs.
Through resolving Ly$\alpha$ emission and studying the morphology, kinematics
and offset to the systemic redshift, analysing H$\alpha$ emission and
obtaining Ly$\alpha$ escape fractions, and better knowledge of the spectra of
galaxies at redshifts $z>6$, surveys such as the ”First Reionization Epoch
Spectroscopically Complete Observations” (FRESCO, Oesch et al., 2023), the
”Cosmic Evolution Early Release Science” (CEERS, Finkelstein et al., 2023), or
GLASS-JWST-ERS (Treu et al., 2022) are already well underway to shedding more
light on this critical phase in the history of the universe.

##  8 Acknowledgements

We would like to thank our colleagues for many interesting discussions and
comments that helped greatly to improve this paper, among which are Rui
Marques-Chaves, Ivan Kramarenko, Moupiya Maji, Laia Barrufet, Sophie van
Mierlo, Vasily Kokorev, Edoardo Iani and Emil Rivera-Thorsen. The plots in
this paper were created using Matplotlib (Hunter, 2007) and we also used the
Python packages Numpy (van der Walt et al., 2011) and Astropy (Astropy
Collaboration et al., 2013). This work received support from the Swiss
National Science Foundation through grant PP00P2_190079\. This work was
supported by the Swiss State Secretariat for Education, Research and
Innovation (SERI) under contract number MB22.00072. The Cosmic Dawn Center
(DAWN) is funded by the Danish National Research Foundation under grant No.
140. GI acknowledges support for this work provided by NASA through grant HST-
GO-13872 awarded by the Space Telescope Science Institute, which is operated
by the Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-26555. HA is supported by CNES, focused on the HST mission, and
the Programme National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and
IN2P3, co-funded by CEA and CNES. Support for this work was provided by NASA
through the NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555. MM acknowledges support from the Project PCI2021-122072-2B,
financed by MICIN/AEI/10.13039/501100011033, and the European Union
“NextGenerationEU”/RTRP and IAC project P/302302. VM acknowledges support from
the NWO grant 016.VIDI.189.162 (“ODIN”).

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 8
  * Bacon et al. (2023) Bacon, R., Brinchmann, J., Conseil, S., et al. 2023, A&A, 670, A4
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1
  * Bassett et al. (2019) Bassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2019, MNRAS, 483, 5223
  * Bassett et al. (2021) Bassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2021, MNRAS, 502, 108
  * Becker et al. (2015) Becker, G. D., Bolton, J. S., & Lidz, A. 2015, PASA, 32, e045
  * Begley et al. (2022) Begley, R., Cullen, F., McLure, R. J., et al. 2022, MNRAS, 513, 3510
  * Bergvall et al. (2006) Bergvall, N., Zackrisson, E., Andersson, B. G., et al. 2006, A&A, 448, 513
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * 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
  * Blaizot et al. (2023) Blaizot, J., Garel, T., Verhamme, A., et al. 2023, MNRAS, 523, 3749
  * Blanc et al. (2011) Blanc, G. A., Adams, J. J., Gebhardt, K., et al. 2011, ApJ, 736, 31
  * Boquien et al. (2022) Boquien, M., Buat, V., Burgarella, D., et al. 2022, A&A, 663, A50
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Borthakur et al. (2014) Borthakur, S., Heckman, T. M., Leitherer, C., & Overzier, R. A. 2014, Science, 346, 216
  * Boutsia et al. (2011) Boutsia, K., Grazian, A., Giallongo, E., et al. 2011, ApJ, 736, 41
  * Bouwens et al. (2007) Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, ApJ, 670, 928
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Burgarella et al. (2005) Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413
  * Byrne et al. (2022) Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Cardamone et al. (2009) Cardamone, C., Schawinski, K., Sarzi, M., et al. 2009, MNRAS, 399, 1191
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, in Interstellar Dust, ed. L. J. Allamandola & A. G. G. M. Tielens, Vol. 135, 5–10
  * Caruana et al. (2012) Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2012, MNRAS, 427, 3055
  * Chisholm et al. (2017) Chisholm, J., Orlitová, I., Schaerer, D., et al. 2017, A&A, 605, A67
  * Chisholm et al. (2020) Chisholm, J., Prochaska, J. X., Schaerer, D., Gazagnes, S., & Henry, A. 2020, MNRAS, 498, 2554
  * Chisholm et al. (2022) Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104
  * Claeyssens et al. (2019) Claeyssens, A., Richard, J., Blaizot, J., et al. 2019, MNRAS, 489, 5022
  * Clarke & Oey (2002) Clarke, C. & Oey, M. S. 2002, MNRAS, 337, 1299
  * Cooke et al. (2014) Cooke, J., Ryan-Weber, E. V., Garel, T., & Díaz, C. G. 2014, MNRAS, 441, 837
  * Cowie et al. (2009) Cowie, L. L., Barger, A. J., & Trouille, L. 2009, ApJ, 692, 1476
  * Cullen et al. (2020) Cullen, F., McLure, R. J., Dunlop, J. S., et al. 2020, MNRAS, 495, 1501
  * de Barros et al. (2016) de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51
  * Deharveng et al. (2001) Deharveng, J.-M., Buat, V., Le Brun, V., et al. 2001, A&A, 375, 805
  * Dijkstra (2014) Dijkstra, M. 2014, PASA, 31, e040
  * Duval et al. (2014) Duval, F., Schaerer, D., Östlin, G., & Laursen, P. 2014, A&A, 562, A52
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058
  * Faisst (2016) Faisst, A. L. 2016, ApJ, 829, 99
  * Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614
  * Faucher-Giguère et al. (2008) Faucher-Giguère, C.-A., Lidz, A., Hernquist, L., & Zaldarriaga, M. 2008, ApJ, 688, 85
  * Feng et al. (2016) Feng, Y., Di-Matteo, T., Croft, R. A., et al. 2016, MNRAS, 455, 2778
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13
  * Finkelstein et al. (2019) Finkelstein, S. L., D’Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36
  * Fletcher et al. (2019) Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87
  * Flury et al. (2022a) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJS, 260, 1
  * Flury et al. (2022b) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022b, ApJ, 930, 126
  * Fontanot et al. (2014) Fontanot, F., Cristiani, S., Pfrommer, C., Cupani, G., & Vanzella, E. 2014, MNRAS, 438, 2097
  * Fujita et al. (2003) Fujita, A., Martin, C. L., Mac Low, M.-M., & Abel, T. 2003, ApJ, 599, 50
  * Garel et al. (2021) Garel, T., Blaizot, J., Rosdahl, J., et al. 2021, MNRAS, 504, 1902
  * Gazagnes et al. (2020) Gazagnes, S., Chisholm, J., Schaerer, D., Verhamme, A., & Izotov, Y. 2020, A&A, 639, A85
  * Gazagnes et al. (2021) Gazagnes, S., Koopmans, L. V. E., & Wilkinson, M. H. F. 2021, MNRAS, 502, 1816
  * Giacconi et al. (2001) Giacconi, R., Rosati, P., Tozzi, P., et al. 2001, ApJ, 551, 624
  * Giallongo et al. (2015) Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83
  * Giavalisco et al. (2004) Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93
  * Grazian et al. (2018) Grazian, A., Giallongo, E., Boutsia, K., et al. 2018, A&A, 613, A44
  * Grazian et al. (2016) Grazian, A., Giallongo, E., Gerbasi, R., et al. 2016, A&A, 585, A48
  * Grazian et al. (2017) Grazian, A., Giallongo, E., Paris, D., et al. 2017, A&A, 602, A18
  * Grimes et al. (2009) Grimes, J. P., Heckman, T., Aloisi, A., et al. 2009, ApJS, 181, 272
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35
  * Gronke & Dijkstra (2014) Gronke, M. & Dijkstra, M. 2014, MNRAS, 444, 1095
  * Guaita et al. (2015) Guaita, L., Melinder, J., Hayes, M., et al. 2015, A&A, 576, A51
  * Guaita et al. (2016) Guaita, L., Pentericci, L., Grazian, A., et al. 2016, A&A, 587, A133
  * Guo et al. (2013) Guo, Y., Ferguson, H. C., Giavalisco, M., et al. 2013, ApJS, 207, 24
  * Haiman (2016) Haiman, Z. 2016, in Astrophysics and Space Science Library, Vol. 423, Understanding the Epoch of Cosmic Reionization: Challenges and Progress, ed. A. Mesinger, 1
  * Harikane et al. (2023) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2303.11946
  * Hayes et al. (2007) Hayes, M., Östlin, G., Atek, H., et al. 2007, MNRAS, 382, 1465
  * Henry et al. (2018) Henry, A., Berg, D. A., Scarlata, C., Verhamme, A., & Erb, D. 2018, ApJ, 855, 96
  * Herenz et al. (2017) Herenz, E. C., Urrutia, T., Wisotzki, L., et al. 2017, A&A, 606, A12
  * Herenz & Wisotzki (2017) Herenz, E. C. & Wisotzki, L. 2017, A&A, 602, A111
  * Herenz et al. (2019) Herenz, E. C., Wisotzki, L., Saust, R., et al. 2019, A&A, 621, A107
  * Hildebrandt et al. (2006) Hildebrandt, H., Erben, T., Dietrich, J. P., et al. 2006, A&A, 452, 1121
  * Hu et al. (1998) Hu, E. M., Cowie, L. L., & McMahon, R. G. 1998, ApJ, 502, L99
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Illingworth et al. (2016) Illingworth, G., Magee, D., Bouwens, R., et al. 2016, arXiv e-prints, arXiv:1606.00841
  * Inami et al. (2017) Inami, H., Bacon, R., Brinchmann, J., et al. 2017, A&A, 608, A2
  * Inoue et al. (2011) Inoue, A. K., Kousai, K., Iwata, I., et al. 2011, MNRAS, 411, 2336
  * Inoue et al. (2014) Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805
  * Iwata et al. (2009) Iwata, I., Inoue, A. K., Matsuda, Y., et al. 2009, ApJ, 692, 1287
  * Izotov et al. (2022) Izotov, Y. I., Chisholm, J., Worseck, G., et al. 2022, MNRAS, 515, 2864
  * Izotov et al. (2016a) Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016a, Nature, 529, 178
  * Izotov et al. (2016b) Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683
  * Izotov et al. (2018a) Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2018a, MNRAS, 474, 4514
  * Izotov et al. (2021) Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2021, MNRAS, 503, 1734
  * Izotov et al. (2018b) Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018b, MNRAS, 478, 4851
  * Japelj et al. (2017) Japelj, J., Vanzella, E., Fontanot, F., et al. 2017, MNRAS, 468, 389
  * Jaskot & Oey (2013) Jaskot, A. E. & Oey, M. S. 2013, ApJ, 766, 91
  * Jaskot & Ravindranath (2016) Jaskot, A. E. & Ravindranath, S. 2016, ApJ, 833, 136
  * Ji et al. (2020) Ji, Z., Giavalisco, M., Vanzella, E., et al. 2020, ApJ, 888, 109
  * Jones et al. (2021) Jones, L. H., Barger, A. J., & Cowie, L. L. 2021, ApJ, 908, 222
  * Jones et al. (2018) Jones, L. H., Barger, A. J., Cowie, L. L., et al. 2018, ApJ, 862, 142
  * Kakiichi & Gronke (2021) Kakiichi, K. & Gronke, M. 2021, ApJ, 908, 30
  * Kashikawa et al. (2011) Kashikawa, N., Shimasaku, K., Matsuda, Y., et al. 2011, ApJ, 734, 119
  * Katz et al. (2022) Katz, H., Garel, T., Rosdahl, J., et al. 2022, MNRAS, 515, 4265
  * Katz et al. (2020) Katz, H., Ďurovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164
  * Keenan et al. (2017) Keenan, R. P., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 848, 12
  * Kerutt (2017) Kerutt, J. 2017, QtClassify: IFS data emission line candidates classifier
  * Kerutt et al. (2022) Kerutt, J., Wisotzki, L., Verhamme, A., et al. 2022, A&A, 659, A183
  * Khaire et al. (2016) Khaire, V., Srianand, R., Choudhury, T. R., & Gaikwad, P. 2016, MNRAS, 457, 4051
  * Kimm et al. (2019) Kimm, T., Blaizot, J., Garel, T., et al. 2019, MNRAS, 486, 2215
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJ, 954, L4
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36
  * Kusakabe et al. (2020) Kusakabe, H., Blaizot, J., Garel, T., et al. 2020, A&A, 638, A12
  * Leclercq et al. (2017) Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8
  * Leethochawalit et al. (2016) Leethochawalit, N., Jones, T. A., Ellis, R. S., Stark, D. P., & Zitrin, A. 2016, ApJ, 831, 152
  * Leitet et al. (2013) Leitet, E., Bergvall, N., Hayes, M., Linné, S., & Zackrisson, E. 2013, A&A, 553, A106
  * Leitet et al. (2011) Leitet, E., Bergvall, N., Piskunov, N., & Andersson, B.-G. 2011, A&A, 532, A107
  * Leitherer et al. (1995) Leitherer, C., Ferguson, H. C., Heckman, T. M., & Lowenthal, J. D. 1995, ApJ, 454, L19
  * Leitherer et al. (2016) Leitherer, C., Hernandez, S., Lee, J. C., & Oey, M. S. 2016, ApJ, 823, 64
  * Ma et al. (2020) Ma, X., Quataert, E., Wetzel, A., et al. 2020, MNRAS, 498, 2001
  * Madau (1995) Madau, P. 1995, ApJ, 441, 18
  * Madau & Haardt (2015) Madau, P. & Haardt, F. 2015, ApJ, 813, L8
  * Mainali et al. (2022) Mainali, R., Rigby, J. R., Chisholm, J., et al. 2022, ApJ, 940, 160
  * Maji et al. (2022) Maji, M., Verhamme, A., Rosdahl, J., et al. 2022, A&A, 663, A66
  * Marchi et al. (2018) Marchi, F., Pentericci, L., Guaita, L., et al. 2018, A&A, 614, A11
  * Marques-Chaves et al. (2021) Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2021, MNRAS, 507, 524
  * Marques-Chaves et al. (2022) Marques-Chaves, R., Schaerer, D., Amorín, R. O., et al. 2022, A&A, 663, L1
  * Mascia et al. (2023) Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155
  * Masters et al. (2012) Masters, D., Capak, P., Salvato, M., et al. 2012, ApJ, 755, 169
  * Matsuoka et al. (2022) Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2022, ApJS, 259, 18
  * Matsuoka et al. (2018) Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, ApJ, 869, 150