# The Low-Redshift Lyman Continuum Survey II:
New Insights into LyC Diagnostics

Sophia R. Flury Department of Astronomy, University of Massachusetts Amherst,
Amherst, MA 01002, United States  Anne E. Jaskot  Department of Astronomy,
Williams College, Williamstown, MA 01267, United States  Harry C. Ferguson
Space Telescope Science Institute, 3700 San Martin Drive Baltimore, MD 21218,
United States  Gábor Worseck Institut für Physik und Astronomie, Universität
Potsdam, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam, Germany  Kirill Makan
Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Str.
24/25, D-14476 Potsdam, Germany  John Chisholm  Department of Astronomy,
University of Texas at Austin, Austin, TX 78712, United States  Alberto
Saldana-Lopez  Observatoire de Genève, Université de Genève, 51 Ch. des
Mailletes, 1290 Versoix, Switzerland  Daniel Schaerer  Observatoire de Genève,
Université de Genève, 51 Ch. des Mailletes, 1290 Versoix, Switzerland  Stephan
R. McCandliss Center for Astrophysical Sciences
Department of Physics and Astronomy, Johns Hopkins University, 3400 North
Charles Street, Baltimore, MD 21218, USA  Bingjie Wang  Department of Physics
and Astronomy, Johns Hopkins University, Baltimore, MD 21218, United States
N. M. Ford  Department of Astronomy, Williams College, Williamstown, MA 01267,
United States  M. S. Oey  Department of Astronomy, University of Michigan, Ann
Arbor, MI 48109, United States  Timothy Heckman  Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, United States
Zhiyuan Ji  Department of Astronomy, University of Massachusetts Amherst,
Amherst, MA 01002, United States  Mauro Giavalisco  Department of Astronomy,
University of Massachusetts Amherst, Amherst, MA 01002, United States  Ricardo
Amorín Instituto de Investigación Multidisciplinar en Ciencia y Tecnología,
Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile  Departamento de
Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La
Serena, Chile  Hakim Atek  Institut d’astrophysique de Paris, CNRS UMR7095,
Sorbonne Université, 98bis Boulevard Arago, F-75014 Paris, France  Jeremy
Blaizot  Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France  Sanchayeeta
Borthakur  School of Earth & Space Exploration, Arizona State University,
Tempe, AZ 85287, United States  Cody Carr  Minnesota Institute for
Astrophysics, School of Physics and Astronomy, University of Minnesota, 316
Church Str SE, Minneapolis, MN 55455, United States  Marco Castellano  INAF,
Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone,
Italy  Stephane De Barros  Observatoire de Genève, Université de Genève, 51
Ch. des Mailletes, 1290 Versoix, Switzerland  Mark Dickinson  National
Optical-Infrared Astronomy Research Laboratory, Tucson, AZ, United States
Steven L. Finkelstein  Department of Astronomy, The University of Texas at
Austin, Austin, TX, United States  Brian Fleming  Laboratory for Atmospheric
and Space Physics, Boulder, Colorado, United States  Fabio Fontanot
INAF–Osservatorio Astronomico di Trieste, Via G.B. Tiepolo, 11, I-34143
Trieste, Italy  Thibault Garel  Observatoire de Genève, Université de Genève,
51 Ch. des Mailletes, 1290 Versoix, Switzerland  Andrea Grazian  INAF-
Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Matthew Hayes  The Oskar Klein Centre, Department of Astronomy,
Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden  Alaina Henry
Space Telescope Science Institute, 3700 San Martin Drive Baltimore, MD 21218,
United States  Valentin Mauerhofer  Observatoire de Genève, Université de
Genève, 51 Ch. des Mailletes, 1290 Versoix, Switzerland  Univ. Lyon, Univ.
Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574,
69230 Saint-Genis-Laval, France  Genoveva Micheva  Leibniz-Institute for
Astrophysics Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany  Goran
Ostlin  The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-10691 Stockholm, Sweden  Casey Papovich  George P. and Cynthia
Woods Mitchell Institute for Fundamental Physics and Astronomy, Department of
Physics and Astronomy, Texas A&M University, College Station, TX, United
States  Laura Pentericci  INAF, Osservatorio Astronomico di Roma, via Frascati
33, I-00078 Monteporzio Catone, Italy  Swara Ravindranath  Space Telescope
Science Institute, 3700 San Martin Drive Baltimore, MD 21218, United States
Joakim Rosdahl  Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France  Michael
Rutkowski  Department of Physics and Astronomy, Minnesota State University,
Mankato, MN, 56001, United States  Paola Santini  INAF, Osservatorio
Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy
Claudia Scarlata  Minnesota Institute for Astrophysics, School of Physics and
Astronomy, University of Minnesota, 316 Church str SE, Minneapolis, MN 55455,
United States  Harry Teplitz  Infrared Processing and analysis Center,
California Institute of Technology, Pasadena, CA 91125, United States  Trinh
Thuan  Astronomy Department, University of Virginia, Charlottesville, VA
22904, United States  Maxime Trebitsch  Astronomy, Kapteyn Astronomical
Institute, Landleven 12, 9747 AD Groningen, The Netherlands  Eros Vanzella
INAF, Osservatorio Astronomico di Bologna, via Gobetti 93/3 I-40129 Bologna,
Italy  Anne Verhamme  Observatoire de Genève, Université de Genève, 51 Ch. des
Mailletes, 1290 Versoix, Switzerland  Univ. Lyon, Univ. Lyon 1, ENS de Lyon,
CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-
Laval, France

###### Abstract

The Lyman continuum (LyC) cannot be observed at the epoch of reionization
($z\gtrsim 6$) due to intergalactic H I absorption. To identify Lyman
continuum emitters (LCEs) and infer the fraction of escaping LyC, astronomers
have developed various indirect diagnostics of LyC escape. Using measurements
of the LyC from the Low-redshift Lyman Continuum Survey (LzLCS), we present
the first statistical test of these diagnostics. While optical depth
indicators based on $\rm Ly\alpha$, such as peak velocity separation and
equivalent width, perform well, we also find that other diagnostics, such as
the [O III]/[O II] flux ratio and star formation rate surface density, predict
whether a galaxy is a LCE. The relationship between these galaxy properties
and the fraction of escaping LyC flux suggests that LyC escape depends
strongly on H I column density, ionization parameter, and stellar feedback. We
find LCEs occupy a range of stellar masses, metallicities, star formation
histories, and ionization parameters, which may indicate episodic and/or
different physical causes of LyC escape.

cosmology: reionization – ISM: interstellar absorption – galaxies: emission
line galaxies, intergalactic medium, star formation

††journal: ApJ††software:  astropy  (Astropy Collaboration et al., 2013,
2018),  matplotlib  (Hunter, 2007),  numpy  (van der Walt et al., 2011),
scipy  (P. Virtanen, 2020)

##  1 Introduction

While the reionization of the Universe likely completed by a redshift of
$z\sim 6$ (e.g., Becker et al., 2001; Paoletti et al., 2020; Pahl et al.,
2020; Yang et al., 2020), substantial contention persists concerning which
galaxies are responsible. Intense, concentrated star formation seems to be a
necessary characteristic of reionizing galaxies because of the associated
stellar feedback (e.g., Heckman et al., 2001; Clarke & Oey, 2002; Heckman et
al., 2011; Sharma et al., 2017; Naidu et al., 2020). Many predictions suggest
that dwarf galaxies are the primary Lyman continuum (LyC) photon donors
because the weak gravitational potentials of low mass galaxies exacerbate the
clearing effects of feedback (e.g., Razoumov & Sommer-Larsen, 2010; Wise et
al., 2014; Paardekooper et al., 2015; Izotov et al., 2021). Others indicate
that more massive galaxies dominate reionization because stronger
gravitational potentials better facilitate efficient star formation, resulting
in a higher LyC photon budget (e.g., Wyithe & Loeb, 2013). An increase in dust
extinction with galaxy mass could reduce LyC escape at the highest masses,
causing galaxies with masses in between these two regimes to dominate
reionization (e.g., Ma et al., 2020).

Central to this debate is the fraction of LyC photons that escape from a
galaxy to reionize the intergalactic medium (IGM). This “escape fraction”
$f_{esc}^{LyC}$ relates the cosmic star formation rate density and LyC
production efficiency to the cosmic ionization rate. Based on results from
stacked spectra at $z\sim 3$ (e.g., Steidel et al., 2018) and from
cosmological simulations (e.g., Finkelstein et al., 2012), galaxies that
exhibit $f_{esc}^{LyC}$$\gtrsim 0.05$ are the most relevant to reionization.
Due to the increase in intergalactic Lyman continuum absorption with redshift,
direct measurement of the LyC, let alone $f_{esc}^{LyC}$, is statistically
unlikely beyond $z\sim 4$ (e.g., Vanzella et al., 2012; Worseck et al., 2014)
to ascertain which galaxies are responsible for reionizing the universe.

Therefore, we must study lower redshift LyC-leaking galaxies to develop
indirect diagnostics for $f_{esc}^{LyC}$. These diagnostics may also provide
insight into the physical mechanism(s) and environmental conditions that
contribute to leaking of LyC photons. In the following, we will explore some
of the most promising diagnostics for $f_{esc}^{LyC}$, including the
equivalent widths (EWs) of nebular emission lines, emission line flux ratios,
the profile of the $\rm Ly\alpha$ emission line, the star formation rate
surface density ($\Sigma_{\rm SFR}$), and the UV continuum slope ($\beta$) and
magnitude ($M_{1500}$).

Of these properties, nebular emission line EWs are one of the easiest to
measure. Because hydrogen Balmer lines originate from recombination across the
entire extent of H II regions, measuring a smaller EW in these lines may
indicate lower optical depth: if gas does not extend out to the Strömgren
radius, not all LyC photons will photoionize H I, and thus, the nebula will
emit fewer Balmer photons (e.g., Bergvall et al., 2013). To date, few tests of
this diagnostic have been conducted. Moreover, many local LyC emitters (LCEs)
exhibiting high Balmer emission line EWs have been found (Izotov et al.,
2016b, 2018a, 2018b). In any case, Balmer line EWs are sensitive not only to
nebula size but also starburst age and star formation history (Zackrisson et
al., 2017; Binggeli et al., 2018; Alavi et al., 2020). Plagued by such
degeneracies, the EW diagnostic may be best when paired with another
indicator.

Similarly motivated by optical depth, the emission line flux ratio $O_{32}$$=$
[O III]$\lambda 5007$/[O II]$\lambda\lambda 3726,3729$ probes the relative
sizes of the O+2 and O+ ionized zones. When the nebula is optically thin to
LyC, the O+ zone will be smaller relative to the O+2 zone (Jaskot & Oey, 2013;
Nakajima & Ouchi, 2014). Albeit straightforward, $O_{32}$ also depends
strongly on ionization parameter to the extent that $O_{32}$ is commonly used
to measure it (e.g., Kobulnicky & Kewley, 2004; Stasińska et al., 2015; Strom
et al., 2018). It is unclear whether this dependence will obfuscate
$f_{esc}^{LyC}$ diagnostics. For instance, $f_{esc}^{LyC}$ itself may depend
on ionization parameter if the incident stellar radiation field is responsible
for a density-bounded scenario wherein gas does not persist out to the
Strömgren radius. $O_{32}$ also does not probe the outer extent of the nebula,
the region that is most sensitive to optical depth. Observationally, there is
substantial scatter in the $f_{esc}^{LyC}$ \- $O_{32}$ relation (e.g., Izotov
et al., 2018b), making this diagnostic tenuous.

Complementing the $O_{32}$ diagnostic, the
$[$O I$]\lambda 6300$ emission line traces neutral gas at the outer extent of
a nebula because O0 has an ionization potential similar to that of hydrogen
($13.62$ eV) and is coupled to hydrogen recombination by charge exchange
reactions (Osterbrock & Ferland, 2006). Comparing the $[$O I$]\lambda 6300$
flux to the flux of another emission line can thus indicate the prevalence of
neutral gas at the nebula boundary. Stasińska et al. (2015) consider the
$O_{31}=$[O III]$\lambda 5007$/[O I]$\lambda 6300$ flux ratio as a diagnostic
for LyC escape. They argue that, while sensitive to physical conditions and
other ionization sources, $O_{31}$ could indicate a density-bounded nebula as
it traces the size of the O0 zone relative to the O+2 zone. Similarly, the [O
I]/H$\beta$ ratio could indicate a deficiency in neutral gas relative to
hydrogen recombining from H II to H I across the nebula, comparable to the
effect described by Wang et al. (2019) for [S II]/H$\alpha$. The [O I] line is
not without complications, though. In models where LyC optical depth and
ionization parameter are not constant across the nebula (the so-called “picket
fence” scenario, cf. Heckman et al., 2001), Ramambason et al. (2020) predict
that $O_{31}$ and [O I]/H$\beta$ depend on conditions such as starburst age,
ionization parameter, and LyC escape conditions. Shocks can also produce [O I]
emission, further complicating the [O I] line as a diagnostic. Even so, these
[O I] flux ratios have promise as indicators of $f_{esc}^{LyC}$. Jaskot et al.
(2019) found that both flux ratios correlate with the $\rm Ly\alpha$ profile
peak velocity separation, itself an indicator of $f_{esc}^{LyC}$ (e.g., Izotov
et al., 2018b).

The $\rm Ly\alpha$ recombination emission line is sensitive to the neutral H I
opacity of a galaxy (e.g., Behrens et al., 2014; Verhamme et al., 2015).
First, the shape of the $\rm Ly\alpha$ profile depends on both the amount of
absorption and re-emission of $\rm Ly\alpha$ photons by H I in the host galaxy
(essentially a “scattering” effect, e.g., Verhamme et al., 2017; Gazagnes et
al., 2020). Second, the line-of-sight H I column density affects the fraction
of $\rm Ly\alpha$ photons that escape (e.g., Henry et al., 2015; Yang et al.,
2017b; Izotov et al., 2018b). The discovery of triple-peaked $\rm Ly\alpha$
profiles associated with LCEs (e.g., Rivera-Thorsen et al., 2019; Izotov et
al., 2018b; Vanzella et al., 2020) suggests low optical depth cavities in the
ISM facilitating both $\rm Ly\alpha$ and LyC photon escape (Byrohl & Gronke,
2020). While these results are promising, the $\rm Ly\alpha$ line can be
sensitive to IGM extinction, metallicity, and stellar populations in the host
galaxy (e.g., Henry et al., 2015; Yang et al., 2017b). Moreover, the
increasing IGM neutral covering can significantly extinguish $\rm Ly\alpha$ at
high redshifts (e.g., Pentericci et al., 2011; Schenker et al., 2014),
particularly the blue peak in the emission line profile (e.g., Laursen et al.,
2011; Yang et al., 2017a; Hayes et al., 2021). Focusing on $\rm Ly\alpha$
emitters may additionally introduce a sampling bias since a non-negligible
fraction of high-redshift Lyman break galaxies and some LCEs exhibit little to
no $\rm Ly\alpha$ (e.g., Vanzella et al., 2009; Ji et al., 2020).
Nevertheless, at reionization, $\rm Ly\alpha$ falls in an ideal ground-based
observing window (Bagley et al., 2017; Yajima et al., 2018), making $\rm
Ly\alpha$ more feasible a diagnostic of $f_{esc}^{LyC}$ than rest-frame
optical lines.

Star formation rate surface density, on the other hand, may serve
simultaneously as a probe of the physical mechanism of LyC leakage and a proxy
for $f_{esc}^{LyC}$. Feedback from star formation can blow bubbles/chimneys in
the host galaxy’s ISM (e.g., Heckman et al., 2001; Clarke & Oey, 2002; Sharma
et al., 2016; Naidu et al., 2020; Gazagnes et al., 2020), suggesting high
$\Sigma_{\rm SFR}$ will correspond to high $f_{esc}^{LyC}$. Previously
detected LCEs also appear to be compact (e.g., Izotov et al., 2018b; Marchi et
al., 2018). Compactness indicates concentrated star formation and may
correlate with LyC emission (e.g., Marchi et al., 2018), although it is yet
unclear whether this is a defining characteristic of all LCEs. High
$\Sigma_{\rm SFR}$ appears to correlate with weaker low-ioniztion metal
absorption lines (e.g., Alexandroff et al., 2015), which suggests a connection
with low optical depth.

All the diagnostics outlined above are plausible but remain, as of yet,
insufficiently explored due to small sample sizes. Future observations with
JWST and other observatories require reliable indirect diagnostics to infer
$f_{esc}^{LyC}$ at reionization. The recent Low-Redshift Lyman Continuum
Survey (LzLCS; Flury et al., 2022, hereafter Paper i) presents an
unprecedented opportunity to investigate the properties of nearby ($z\sim
0.3$) LCEs to test each diagnostic. In this paper, we summarize the survey in
§2 but refer the reader to Paper i for details. We provide assessments of the
success of indirect diagnostics for selecting LCEs and inferring
$f_{esc}^{LyC}$ (§3). In §4, we compare different parameters and consider the
implications both for diagnostics and for populations of LCEs. Finally, we
discuss our results in the contexts of different LyC escape scenarios,
cosmological simulations, and high redshift surveys in §5. As in Paper i, we
assume $H_{0}=70$ km s-1 Mpc-1, $\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$.

##  2 The LzLCS Sample

In Paper i, we outline the LzLCS sample selection, data processing, LyC
measurements and ancillary measurements. Here, we summarize the LzLCS but
refer the reader to Paper i for details. The LzLCS consists of _HST_ /COS
G140L observations of 66 new LCE candidates in the nearby ($z\sim 0.3-0.4$)
Universe. The significance of the LzLCS over previous surveys is not just its
size but also its scope: the LzLCS spans a much broader range in, among other
properties, $O_{32}$, UV $\beta$, $\Sigma_{\rm SFR}$, stellar mass,
metallicity, and burst age than prior studies.

From the COS spectra, we measure the LyC and $\rm Ly\alpha$ flux, the $\rm
Ly\alpha$ EW, and the spectral slope of the attenuated starlight continuum. To
assess the detection of the LyC, we define the probability $P(>N|B)$ given by
the Poisson survival function that the observed signal arises by chance from
background fluctuations (see Worseck et al., 2016; Makan et al., 2020). In
summary, we detect LyC emission from 35 galaxies at $>$97.725% confidence
($P(>N|B)<0.02275$, $>2\sigma$ significance) with 13/35 of these galaxies
having fairly-detected LyC flux ($P(>N|B)<0.00135$, $3-5\sigma$ significance)
and 12/35 of these galaxies having well-detected LyC flux
($P(>N|B)<2.867\times 10^{-7}$, $>5\sigma$ significance). Following Chisholm
et al. (2019), we also perform SED fitting using instantaneous burst model
spectra from Starburst99 (Leitherer et al., 2010) to infer the intrinsic and
observed UV magnitudes at 1500 Å, starlight attenuation, starburst age, and
intrinsic LyC flux (see Saldana-Lopez et al., 2022, for details). From the COS
acquisition images, we determine the UV half-light radius from the background-
subtracted empirical cumulative light distribution, which is consistent with
results from Sérsic profile fits.

Optical spectra from the SDSS provide additional information. We measure
nebular emission line fluxes and EWs from these spectra. To correct the fluxes
for reddening, we infer nebular extinction from Balmer lines assuming Case B
recombination with Storey & Hummer (1995) atomic data and the Cardelli et al.
(1989) reddening law. From the extinction-corrected fluxes, we obtain
diagnostic flux ratios like $O_{32}$. Since [O III]$\lambda 4363$ is detected
in the majority of the galaxies in our sample, we determine gas-phase
metallicities by the direct method, inferring the auroral line flux by the ff-
relation (Pilyugin et al., 2006) when not available.

Combining the SDSS and COS measurements with GALEX observations, we infer
additional properties of the sample. We obtain stellar masses by fitting
aperture-matched photometry from SDSS and GALEX using Prospector (Leja et al.,
2017; Johnson et al., 2019). From the Balmer-line star formation rate
(Kennicutt & Evans, 2012), UV half-light radius, and stellar mass, we
calculate the star formation rate surface density and specific star formation
rate. From the H$\beta$ and Ly$\alpha$ fluxes, we also determine the fraction
of $\rm Ly\alpha$ photons $f_{esc}^{Ly\alpha}$ which have escaped the host
galaxy.

As discussed in §1, the key objective of establishing LyC escape diagnostics
is to determine which properties most strongly correlate with the LyC escape
fraction $f_{esc}^{LyC}$. We measure $f_{esc}^{LyC}$ by three different
methods: (i) the empirical $F_{\lambda LyC}/F_{\lambda 1100}$ ratio, (ii)
using H$\beta$ to infer the intrinsic LyC from Starburst99 continuous
starburst model spectra ($f_{esc}^{LyC}$$({\rm H\beta})$), and (iii) using
spectral energy distribution (SED) fits to the COS FUV spectra to infer the
intrinsic LyC ($f_{esc}^{LyC}$$({\rm UV})$). We summarize these metrics below.

  1. i.

The $F_{\lambda LyC}/F_{\lambda 1100}$ ratio is simply a ratio of flux density
measurements and is therefore free of any assumptions. However, $F_{\lambda
LyC}/F_{\lambda 1100}$ is sensitive to a number of factors in addition to the
LyC escape fraction, including dust attenuation, star formation rate,
metallicity, and burst age, complicating the interpretation of this quantity
as a metric of $f_{esc}^{LyC}$.

  2. ii.

The H$\beta$ approach introduced by Izotov et al. (2016b) uses the H$\beta$
flux, a measure of the total number of absorbed ionizing photons, the H$\beta$
EW, a proxy for burst age, and the gas phase metallicity to determine which
Starburst99 model to use to predict the intrinsic LyC flux density
corresponding to the measured LyC. The observed LyC flux density relative to
the intrinsic LyC flux density then gives the $f_{esc}^{LyC}$. This method can
be sensitive to the assumed star formation history, particularly for older
burst ages or for galaxies with high stellar mass. In Paper i, we find that a
continuous star formation scenario provides the best agreement with the other
two $f_{esc}^{LyC}$ metrics.

  3. iii.

The UV approach relies on fitting the COS UV spectra over the rest-frame
wavelength range of 925-1345 Å. These SED fits consist of a uniform Reddy et
al. (2016) attenuation law applied to a weighted linear combination of single-
burst Starburst99 spectra and accompanying nebular continua (see Chisholm et
al., 2019; Saldana-Lopez et al., 2022). The best-fit SED allows us to predict
the intrinsic LyC flux density, which we then use to obtain the
$f_{esc}^{LyC}$. While the UV fits are still affected by assumptions about the
star formation history and dust attenuation, the COS UV $f_{esc}^{LyC}$ is
less sensitive to systematics than the H$\beta$ $f_{esc}^{LyC}$, making this
metric the most reliable of the three.

We discuss the systematics and caveats of each $f_{esc}^{LyC}$ measurement in
detail in Paper i. Since each metric has its caveats, we choose to consider
all three in our assessment of $f_{esc}^{LyC}$ diagnostics below.

In addition to the 66 LzLCS galaxies, we include 23 LCE candidates which have
been observed by _HST_ /COS and previously published in other studies (Izotov
et al., 2016a, b, 2018a, 2018b, 2021; Wang et al., 2019). To obtain LyC fluxes
and $f_{esc}^{LyC}$ values in a consistent manner, we re-process their HST/COS
observations. Our significance estimates are comparable to the published
values. However, our more stringent $P(>N|B)<0.02275$ detection criterion
includes two fewer LCEs than previous assessments, yielding a total of 15 LyC
detections out of 23 published observations. For consistency, we also re-
measure and infer the diagnostic properties from SDSS and COS observations in
the same manner as the LzLCS galaxies. A detailed comparison of these
published LCE candidates with the LzLCS is presented in Paper i.

##  3 Indirect ${f_{esc}^{lyc}}$ Diagnostics

The LzLCS aims to test various $f_{esc}^{LyC}$ diagnostics proposed in the
literature as discussed in §1. To that end, we compare our measured
$f_{esc}^{LyC}$ against each parameter in Figures 1-17. As discussed in §2, we
consider three different metrics of $f_{esc}^{LyC}$ since each is susceptible
to various systematic uncertainties. Another reason to use different
$f_{esc}^{LyC}$ indicators stems from the fact that one of the $f_{esc}^{LyC}$
measurements is sometimes obtained from the same data as the proposed indirect
diagnostics. Additional $f_{esc}^{LyC}$ metrics allow us to assess whether an
apparent correlation is real or simply introduced by using the same data in
the abscissa and ordinate variables.

To quantify possible correlations, we compute the Kendall $\tau$ rank
correlation coefficient following the Akritas & Siebert (1996) prescription
for censored data to account for upper limits on $f_{esc}^{LyC}$ when
$P(>N|B)>0.02275$ and any left- or right-censored indirect indicators of
$f_{esc}^{LyC}$. For non-emitters, we use $f_{esc}^{LyC}$ derived from the
$1\sigma$ upper limit on the LyC flux for all three $f_{esc}^{LyC}$
indicators. If we impose a $3\sigma$ detection requirement on the LyC instead
of the adopted $2\sigma$ criterion, we find no significant change in the value
of $\tau$. We report $\tau$ for each of the three $f_{esc}^{LyC}$ indicators
($F_{\lambda LyC}/F_{\lambda 1100}$, $f_{esc}^{LyC}$(H$\beta$), and
$f_{esc}^{LyC}$(UV)) in Table 1 for the combined LzLCS and Izotov et al.
(2016a, b, 2018a, 2018b, 2021); Wang et al. (2019) samples. We also report the
probability $p$ that the measured $\tau$ is consistent with the null
hypothesis that $f_{esc}^{LyC}$ is not correlated with the indirect
diagnostic. We consider correlations to be significant if $p<1.350\times
10^{-3}$ ($3\sigma$ confidence) and strong if $|\tau|>0.2$.

Table 1: Kendall $\tau$ correlation coefficients for proposed LyC $f_{esc}^{LyC}$ diagnostics for the combined LzLCS and Izotov et al. (2016a, b, 2018a, 2018b, 2021); Wang et al. (2019) sample, accounting for upper limits in $F_{\lambda LyC}/F_{\lambda 1100}$ and $f_{esc}^{LyC}$ following Akritas & Siebert (1996). $p$ values indicate the false-positive probability that the correlation is real. $\sigma$ values, given by ${\rm probit}(1-p)$, indicate the significance of the correlation, i.e., the number of standard deviations separating the correlation statistic from that of the null hypothesis. The characteristic $1\sigma$ uncertainties in $\tau$ estimated by bootstrapping values of the $f_{esc}^{LyC}$ metrics and indirect diagnostics are $\pm 0.05$. Each correlation is assessed using all 89 objects in the combined LzLCS sample except $v_{sep}$, which is only measured for 7 LzLCS and 20 published objects.  Diagnostic | $F_{\lambda LyC}/F_{\lambda 1100}$ |  $f_{esc}^{LyC}$(H$\beta$) |  $f_{esc}^{LyC}$(UV)
---|---|---|---
| $\tau$ | $p$ | $\sigma$ | $\tau$ | $p$ | $\sigma$ | $\tau$ | $p$ | $\sigma$
$f~{}_{esc}^{Ly\alpha}$ | $0.292$ | $5.186\times 10^{-5}$ | $3.882$ | $0.343$ | $1.942\times 10^{-6}$ | $4.618$ | $0.324$ | $6.774\times 10^{-6}$ | $4.351$
EW(Ly$\alpha$) | $0.320$ | $8.687\times 10^{-6}$ | $4.296$ | $0.234$ | $1.141\times 10^{-3}$ | $3.051$ | $0.342$ | $2.011\times 10^{-6}$ | $4.610$
$v_{sep}$ | $-0.493$ | $3.103\times 10^{-4}$ | $3.422$ | $-0.422$ | $2.033\times 10^{-3}$ | $2.873$ | $-0.530$ | $1.055\times 10^{-4}$ | $3.705$
$\log_{10}O_{31}$ | $-0.149$ | $0.039$ | $1.761$ | $-0.144$ | $0.045$ | $1.693$ | $-0.151$ | $0.036$ | $1.796$
$\log_{10}$[O I]/H$\beta$ | $-0.148$ | $0.041$ | $1.745$ | $-0.145$ | $0.044$ | $1.709$ | $-0.145$ | $0.044$ | $1.705$
$\log_{10}O_{32}$ | $0.290$ | $5.678\times 10^{-5}$ | $3.860$ | $0.198$ | $6.024\times 10^{-3}$ | $2.511$ | $0.347$ | $1.438\times 10^{-6}$ | $4.679$
EW(H$\beta$) | $0.223$ | $1.953\times 10^{-3}$ | $2.886$ | $0.109$ | $0.132$ | $1.117$ | $0.283$ | $8.366\times 10^{-5}$ | $3.764$
M1500,obs | $0.045$ | $0.533$ | $0.000$ | $-0.013$ | $0.857$ | $0.000$ | $0.098$ | $0.174$ | $0.940$
M1500,int | $0.228$ | $1.591\times 10^{-3}$ | $2.950$ | $0.157$ | $0.029$ | $1.895$ | $0.320$ | $8.978\times 10^{-6}$ | $4.289$
$\beta_{1200}$ | $-0.221$ | $2.200\times 10^{-3}$ | $2.848$ | $-0.261$ | $2.966\times 10^{-4}$ | $3.435$ | $-0.283$ | $8.366\times 10^{-5}$ | $3.764$
$\log_{10}M\star$ | $-0.089$ | $0.216$ | $0.785$ | $-0.074$ | $0.307$ | $0.503$ | $-0.167$ | $0.021$ | $2.040$
COS NUV $r_{50}$ | $-0.388$ | $7.179\times 10^{-8}$ | $5.261$ | $-0.301$ | $2.938\times 10^{-5}$ | $4.018$ | $-0.382$ | $1.193\times 10^{-7}$ | $5.166$
$\rm\log_{10}\Sigma_{SFR,H\beta}$ | $0.368$ | $3.884\times 10^{-7}$ | $4.941$ | $0.264$ | $2.650\times 10^{-4}$ | $3.465$ | $0.325$ | $7.099\times 10^{-6}$ | $4.341$
$\rm\log_{10}\Sigma_{SFR,F_{\lambda 1100}}$ | $0.070$ | $0.334$ | $0.429$ | $0.068$ | $0.347$ | $0.394$ | $-0.035$ | $0.632$ | $0.000$
$\rm\log_{10}{sSFR}$ | $0.110$ | $0.128$ | $1.138$ | $0.043$ | $0.554$ | $0.000$ | $0.181$ | $0.012$ | $2.254$
$\rm\log_{10}\Sigma_{sSFR,H\beta}$ | $0.290$ | $6.320\times 10^{-5}$ | $3.833$ | $0.208$ | $4.167\times 10^{-3}$ | $2.638$ | $0.346$ | $1.859\times 10^{-6}$ | $4.627$
$\rm 12+\log_{10}\left(\frac{O}{H}\right)$ | $-0.187$ | $9.484\times 10^{-3}$ | $2.346$ | $-0.130$ | $0.070$ | $1.475$ | $-0.211$ | $3.420\times 10^{-3}$ | $2.705$

Below, we examine diagnostics motivated by LyC optical depth and by properties
related to stellar feedback mechanisms facilitating LyC escape. Many proxies
exhibit an intrinsic, often significant correlation with $f_{esc}^{LyC}$ but
with substantial scatter in the data. This scatter often extends from a
parameter-dependent upper bound on $f_{esc}^{LyC}$ down to
$f_{esc}^{LyC}$$=0$. These upper bounds in $f_{esc}^{LyC}$ suggest that line
of sight effects (e.g., orientation, covering fraction, etc.) may be
responsible for obscuring otherwise clear trends, although relationships
between various properties may also contribute to the scatter.

Figure 1: Comparison of $F_{\lambda LyC}/F_{\lambda 1100}$ (left),
$f_{esc}^{LyC}$ from H$\beta$ (center), and $f_{esc}^{LyC}$ from fitting the
COS UV spectrum (right) with $f_{esc}^{Ly\alpha}$ for strong LCEs ($>5\sigma$
LyC detection and $f_{esc}^{LyC}$$>0.05$, filled black), weak LCEs ($>2\sigma$
LyC detection but not strong, filled grey), and non-emitters ($<2\sigma$ LyC
detection, open grey). For non-emitters, $f_{esc}^{LyC}$ is shown for the
$1\sigma$ upper limit on the LyC flux. LzLCS results are shown as squares
(strong LCEs), circles (weak LCEs), and triangles (non-detections) while the
Izotov et al. (2016a, b, 2018a, 2018b, 2021) galaxies are shown as stars and
the Wang et al. (2019) are shown as diamonds. Characteristic $1\sigma$
uncertainties are shown in the upper right. Dotted line indicates 1:1
equivalence.

To illustrate the combination of obscuring scatter and upper-bound effects, we
refer to Figure 1 comparing $f_{esc}^{LyC}$ and $f_{esc}^{Ly\alpha}$. For each
of the LyC metrics, the highest $f_{esc}^{LyC}$ values at any given
$f_{esc}^{Ly\alpha}$ closely follow the dotted line indicating where
$f_{esc}^{LyC}$=$f_{esc}^{Ly\alpha}$. However, the lowest $f_{esc}^{LyC}$
values at any given $f_{esc}^{Ly\alpha}$ are largely invariant with
$f_{esc}^{Ly\alpha}$, consistent with $f_{esc}^{LyC}$$\approx 10^{-3}$
estimated in the non-detections. Thus, in this example, the identity line
$f_{esc}^{LyC}$=$f_{esc}^{Ly\alpha}$ is an upper bound, or “envelope”, to the
distribution of $f_{esc}^{LyC}$values at a particular value of
$f_{esc}^{Ly\alpha}$. The envelope represents the intrinsic relationship
between $f_{esc}^{LyC}$ and $f_{esc}^{Ly\alpha}$, a trend obscured by the
scatter apparent in the observed distribution. Here, this upper bound could
originate from a line of sight effect. When optically thin regions align with
the aperture, $f_{esc}^{LyC}$=$f_{esc}^{Ly\alpha}$. When optically thin
regions are not aligned with the aperture, $f_{esc}^{LyC}$ decreases
significantly while resonant scattering of $\rm Ly\alpha$ photons into the
line of sight maintains a high $f_{esc}^{Ly\alpha}$. This upper bound effect
appears in the $f_{esc}^{Ly\alpha}$, [O I], $O_{32}$, EW H$\beta$, $M_{1500}$,
$\beta_{1200}$, $M_{\star}$, $r_{50}$, sSFR, and $\rm 12+\log_{10}(O/H)$
diagnostics for at least one of the LyC escape metrics, although the scatter
down to $f_{esc}^{LyC}$=0 tends to be more pronounced in $F_{\lambda
LyC}/F_{\lambda 1100}$ and H$\beta$ $f_{esc}^{LyC}$.

To highlight the intrinsic trends obscured by downward scatter in
$f_{esc}^{LyC}$, we calculate the fraction of galaxies that are prodigious
LCEs (LyC detected at $>5\sigma$ and either $f_{esc}^{LyC}$ or $F_{\lambda
LyC}/F_{\lambda 1100}$$>0.05$) in each diagnostic. This “LCE fraction”
indicates the affinity of strong LCEs for certain global properties. Such a
preference could be used to select LCEs in future studies.

We define the LCE fraction as the ratio of strong LCEs ($P(>N|B)<2.867\times
10^{-7}$ for $5\sigma$ significance as in §2, and $f_{esc}^{LyC}$ indicator
exceeds 5%) to total galaxies in the combined sample in a given bin of a
particular property. Bin counts are given as the sum of the probabilities for
the strong LCE and total LCE samples, i.e., the Poisson binomial expected
value for each bin. We define the final LCE detection fraction as the ratio of
the strong LCE counts to that of the total counts in a given bin. We determine
uncertainty in the LCE detection fraction using a Monte Carlo simulation of
the ratio of Poisson binomial distributions using $10^{4}$ independent
Bernoulli trials for each measurement to determine the distribution of counts
in each bin. In cases of upper limits where no uncertainty is reported for a
measurement (notably the [O I] line), we follow the root-finding minimization
method from Gehrels (1986) to define confidence intervals on the detection
fraction. The $F_{\lambda LyC}/F_{\lambda 1100}$ flux ratio tends to agree
more with $f_{esc}^{LyC}$(H$\beta$) than with $f_{esc}^{LyC}$(UV); however,
$f_{esc}^{LyC}$(UV) is less dependent on assumptions about stellar populations
and star formation history than $f_{esc}^{LyC}$(H$\beta$) and contains
corrections for dust that are absent in $F_{\lambda LyC}/F_{\lambda 1100}$.
Despite some disagreements in values of LCE fraction, $f_{esc}^{LyC}$, and
$F_{\lambda LyC}/F_{\lambda 1100}$, the trend in LCE fraction does not change
significantly depending on the $f_{esc}^{LyC}$ indicator.

###  3.1 $\rm Ly\alpha$

Figure 2: Fraction of galaxies in a given bin of $f_{esc}^{Ly\alpha}$ (left),
EW($\rm Ly\alpha$) (center), or $v_{sep}$ (right) that are strong LCEs
according to $F_{\lambda LyC}/F_{\lambda 1100}$ (grey diamonds),
$f_{esc}^{LyC}$(H$\beta$) (open squares), and $f_{esc}^{LyC}$(UV) (filled
black circles) for the combined LzLCS and published sampels.

As both $\rm Ly\alpha$ and the LyC are sensitive to the line of sight, we
anticipate the former to trace the LCE fraction rather well even as the
emission line experiences higher optical depths than the continuum for the
same column density. We demonstrate this relation is the case in Figure 2. The
LCE fraction increases nearly monotonically with both $f_{esc}^{Ly\alpha}$ and
the $\rm Ly\alpha$ EW. Deviations from this trend occur at high
$f_{esc}^{Ly\alpha}$ due to two galaxies having abnormally high $\rm
Ly\alpha$/H$\beta$ ratios corresponding to values of $f_{esc}^{Ly\alpha}$
close to $1$.

The overwhelming success of all three $\rm Ly\alpha$ diagnostics in
correlating with $f_{esc}^{LyC}$ demonstrates that both $\rm Ly\alpha$ and LyC
escape depend on the distribution of neutral hydrogen in a galaxy. While all
three are promising candidates for $f_{esc}^{LyC}$ indicators at high
redshift, $\rm Ly\alpha$ EW is perhaps the easiest to measure as it requires
neither rest-frame optical nor high resolution spectroscopy. The wealth of
$\rm Ly\alpha$ EW measurements of LCEs at higher redshifts ($z\sim 3-4$, e.g.,
Steidel et al., 2018; Marchi et al., 2018; Fletcher et al., 2019; Pahl et al.,
2021) suggests these diagnostics will likely remain applicable for galaxies at
the epoch of reionization. We discuss this possibility in detail in §5.3.

####  3.1.1 $\rm Ly\alpha$ Escape Fraction

In Paper i, we derived $f_{esc}^{Ly\alpha}$ from $\rm Ly\alpha$ and H$\beta$
emission line fluxes using Case B recombination coefficients determined by the
measured electron temperatures and densities. The $f_{esc}^{LyC}$exhibits a
strong correlation with the fraction of escaping $\rm Ly\alpha$ photons.
Consistent with previous studies (e.g., Dijkstra et al., 2016; Verhamme et
al., 2017; Izotov et al., 2020), we typically find that
$f_{esc}^{Ly\alpha}$$\gtrsim$$f_{esc}^{LyC}$ and that the envelope of
$f_{esc}^{LyC}$ increases with increasing $f_{esc}^{Ly\alpha}$(Figure 1). This
result illustrates that $\rm Ly\alpha$ has an advantage over the LyC in
escaping a galaxy because it can resonantly scatter. However, $\rm Ly\alpha$
and LyC escape have a clear physical connection: our results in Figure 1
suggest only $\rm Ly\alpha$ emitters (LAEs) can be strong LCEs. Indeed, from
Figure 1, $f_{esc}^{LyC}$$=$$f_{esc}^{Ly\alpha}$ serves as an approximate
envelope to the distribution of $f_{esc}^{LyC}$.

There remains a great deal of scatter in $f_{esc}^{LyC}$ values for any given
$f_{esc}^{Ly\alpha}$, which is to be expected. Observed $\rm Ly\alpha$ flux
depends on scattering effects from dust and gas along and outside the line of
sight, making $f_{esc}^{Ly\alpha}$ sensitive to dust and gas content, ISM
clumpiness, and gas dynamics (e.g., Dijkstra et al., 2016). While ISM content
and geometry also affect $f_{esc}^{LyC}$, the scatter evident in Figure 1
suggests the dependence of $f_{esc}^{LyC}$ on these properties may differ from
that of $f_{esc}^{Ly\alpha}$.

####  3.1.2 $\rm Ly\alpha$ Equivalent Width

Figure 3: Same as Figure 1 but for $\rm Ly\alpha$ EW. Dashed line is the
relation from Pahl et al. (2021) based on measurements of the LyC from stacks
of spectra from galaxies at $z\sim 3$.

In place of $f_{esc}^{Ly\alpha}$, we consider $\rm Ly\alpha$ EW as another
possible indicator of $f_{esc}^{LyC}$. While less directly tracing optical
depth than $f_{esc}^{Ly\alpha}$, $\rm Ly\alpha$ EW is far easier to measure at
high redshift because, unlike $f_{esc}^{Ly\alpha}$, it does not require
observing the Balmer lines. Furthermore, both Henry et al. (2015) and Yang et
al. (2017b) demonstrate a strong correlation between $f_{esc}^{Ly\alpha}$ and
$\rm Ly\alpha$ EW for Green Peas (GPs).

From Figure 3, $f_{esc}^{LyC}$ is bounded by an upper limit that increases
with $\rm Ly\alpha$ EW. We find that the relation between $\rm Ly\alpha$ EW
and $f_{esc}^{LyC}$ given by Pahl et al. (2021) provides a reasonable
description of this envelope in $f_{esc}^{LyC}$ values up to $\sim 150$ Å.
Interestingly, this trend extends beyond the 110 Å maximum predicted by
Steidel et al. (2018), perhaps due to their assumed continuous star formation
history.

The correlation between $f_{esc}^{LyC}$ and $\rm Ly\alpha$ EW is one of the
stronger, more significant of the diagnostics we consider, with $\tau\approx
0.3$ across all three $f_{esc}^{LyC}$ indicators. Ten of the sixteen galaxies
with high $\rm Ly\alpha$ EW are strong leakers, indicating a possible
transition in $\rm Ly\alpha$ EW at 100 Å. Galaxies with EWs above this value
are far more likely to be strong LCEs than those below it (see also Figure 2).
Further supporting this notion of a 100 Å transition, LCEs with EWs below this
value are not as prodigious emitters as above it, which may suggest LCEs with
high EW contribute more significantly to the cosmic LyC budget.

$\rm Ly\alpha$ is sensitive to both column density and $f_{esc}^{LyC}$, which
affects the distribution of $f_{esc}^{LyC}$ over the EWs. As the optical depth
decreases, more $\rm Ly\alpha$ and LyC photons can escape. The drop in column
density will then begin to limit the intrinsic emission of $\rm Ly\alpha$
photons, causing $\rm Ly\alpha$ EW to decrease with increasing $f_{esc}^{LyC}$
(Nakajima & Ouchi, 2014; Steidel et al., 2018). This effect may explain why
some strong LCEs persist at $\rm Ly\alpha$ EWs at or even below 100 Å.
However, $\rm Ly\alpha$ EW also depends very strongly on the continuum flux
and, by extension, the stellar population(s), making unclear how significantly
$f_{esc}^{LyC}$ affects the $\rm Ly\alpha$ EW. Nevertheless, $\rm Ly\alpha$ EW
is still one of the most promising diagnostics for $f_{esc}^{LyC}$, performing
comparably to $f_{esc}^{Ly\alpha}$.

####  3.1.3 $\rm Ly\alpha$ Peak Velocity Separation

$\rm Ly\alpha$ profiles are typically double-peaked with the velocity
separation of the two peaks directly depending on the H I column density
(e.g., Verhamme et al., 2015). Previously, Izotov et al. (2018b, 2021)
demonstrated that the velocity separation $v_{sep}$ of $\rm Ly\alpha$ peaks
strongly correlates with $f_{esc}^{LyC}$ as a result of the strong dependence
of both values on the H I column. We compare $f_{esc}^{LyC}$ with $v_{sep}$
for a subset of 7 LzLCS galaxies with G160M measurements (Henry et al., 2015;
Yang et al., 2017b) and a subset of 20 published LCEs (Verhamme et al., 2017;
Izotov et al., 2018b, 2021) in Figure 4. We find a strong relationship between
the values. Neither weak leakers nor non-detections persist below $v_{sep}\sim
250$ km s-1 regardless of $f_{esc}^{LyC}$ metric, which is consistent with
predictions by Verhamme et al. (2015). This result favors the H I column
density interpretation that motivated this diagnostic.

In Figure 4, we also compare $f_{esc}^{LyC}$ and $v_{sep}$ to the relation
from Izotov et al. (2018b) and find that their fit roughly describes the
envelope of $f_{esc}^{LyC}$ values. Moreover, we find $v_{sep}$ has one of the
most pronounced correlations of any indirect $f_{esc}^{LyC}$ indicator
regardless of $f_{esc}^{LyC}$ metric with $\tau\sim-0.4$. However, additional
measurements of $v_{sep}$ from $\rm Ly\alpha$ profiles are necessary to fully
test this $f_{esc}^{LyC}$ diagnostic.

Figure 4: Same as Figure 1 but for $\rm Ly\alpha$ peak velocity separation
$v_{sep}$. Results include only 7 galaxies from the LzLCS plus the Izotov et
al. (2016a, b, 2018a, 2018b, 2021) samples. Dashed line is the description by
Izotov et al. (2018b).

###  3.2 $[$O I$]\lambda 6300$

Figure 5: Same as Figure 1 but for $O_{31}$ and [O I]/H$\beta$.

Few of the LCEs have detected [O I]$\lambda 6300$, with 39 of the 50 LCEs in
the combined sample having only upper limits in the emission line flux,
indicating that [O I] may be particularly weak in LCEs (see, however, Plat et
al., 2019; Ramambason et al., 2020). The preponderance of [O I] upper limits
in LCEs over non-emitters prevents determination of the shape of the
distribution of $f_{esc}^{LyC}$ with respect to the [O I] flux ratios, as
evidence in Figure 5. These upper limits may contribute to the lack of
significant correlation (as shown in Table 1, $p>0.09$ for all three
$f_{esc}^{LyC}$ metrics for both [O I] flux ratios).

Figure 6: Same as Figure 2 but for $O_{31}$ (left) and [O I]/H$\beta$ (right).

We compute the [O I] LCE fraction using upper limits when no flux is detected.
The LCE fraction distribution (Fig. 6) appears relatively monotonic over the
[O I] flux ratios $O_{31}$ and [O I]/H$\beta$. As shown in Figure 6, there is
little difference between the three LyC escape indicators. Above
$\log_{10}O_{31}\sim 1.8$ and below [O I]/H$\beta\sim-1.3$, the LCE fraction
is a roughly constant 30%. The trends with both flux ratios suggest that LCE
detection fraction increases with decreasing [O I] flux.

###  3.3 $O_{32}$

Figure 7: Same as Figure 1 but for the extinction-corrected $O_{32}$ flux
ratio. Dashed line indicates the Izotov et al. (2018a) relation. Note that the
number of strong LCEs changes depending on the $f_{esc}^{LyC}$ metric. Figure
8: Same as Figure 2 but for $O_{32}$.

As shown in Figure 7, we see increasing average and maximum $f_{esc}^{LyC}$
with increasing $O_{32}$. The total combined sample exhibits a significant
($|\tau|>0.2$ at $>3\sigma$ confidence) correlation with $O_{32}$ in all three
$f_{esc}^{LyC}$ metrics (see Table 1). The $f_{esc}^{LyC}$ relation proposed
in Izotov et al. (2018a, dashed line in Figure 7) serves as a description of
the envelope in $f_{esc}^{LyC}$ values, particularly in the case of the UV-
derived escape fraction (right panel of Figure 7). For the other two
$f_{esc}^{LyC}$ metrics, several strong LCEs have $f_{esc}^{LyC}$ in excess of
the predicted relation, notably the two LCEs from Wang et al. (2019). These
outlying galaxies have higher stellar masses $>10^{9}$ M⊙, which makes values
life $F_{\lambda LyC}/F_{\lambda 1100}$ and $f_{esc}^{LyC}$(H$\beta$) more
sensitive to the star formation history. Because the majority of the outlying
objects exhibit agreement with this envelope when using the UV $f_{esc}^{LyC}$
metric (which assumes a non-parametric star formation history), we interpret
these outliers as having $f_{esc}^{LyC}$ systematically biased by this effect
rather than evidence that the Izotov et al. (2018b) prescription is not
appropriate at low $O_{32}$.

Scatter of $f_{esc}^{LyC}$ below the envelope could be caused by a number of
effects related to line-of-sight optical depth and other galaxy properties.
Galaxy-to-galaxy variations in opening angle, covering fraction, burst age,
ionization parameter, star-formation rate, metallicity, extinction law and
dust, and/or optical depth may indicate that additional physical properties
are relevant to determining $f_{esc}^{LyC}$. Unfortunately the degenerate
dependence of $O_{32}$ on metallicity and ionization parameter substantially
complicates any physical interpretation of this result (Sawant et al., 2021).
While the mass-metallicity relation appears to evolve with redshift (Sanders
et al., 2021), the coupling of $\rm 12+\log(O/H)$ to ionization parameter
observed at low redshift persists at $z\sim 2-4$ (Sanders et al., 2020).
Therefore, the underlying properties which affect $O_{32}$ remain related in
the same way regardless of epoch, supporting the extension of the $O_{32}$
diagnostic to higher redshifts. Regardless of the physical connection between
$O_{32}$ and $f_{esc}^{LyC}$, the scatter in this diagnostic indicates that an
isotropic density-bounded escape scenario alone is likely incompatible with
the observations. We discuss interpretations and caveats of $O_{32}$ as an
$f_{esc}^{LyC}$ diagnostic in greater detail in §4 and §5.

With these nuances to $O_{32}$ in mind, we do find that $f_{esc}^{LyC}$
correlates strongly with $O_{32}$, particularly for $f_{esc}^{LyC}$ derived
from the UV spectrum. Further evidence for this relationship is the LCE
fraction shown in Figure 8. The highest $O_{32}$ galaxies in the combined
sample ($O_{32}$$\gtrsim$10) have LCE fractions $\gtrsim 0.5$. However, LCE
fractions $>0.1$ persist across the full range of $O_{32}$ values regardless
of $f_{esc}^{LyC}$ indicator. As we demonstrate in Figure 7, while most
galaxies with $O_{32}$$>5$ in the combined sample are strong LCEs along the
line of sight, some galaxies with lower $O_{32}$ can still be strong LCEs. Two
scenarios might contribute to this trend with $O_{32}$: separate populations
of LCEs with LyC escaping under different conditions or an evolutionary
sequence in which LyC escape recurs at later times when fewer early type stars
persist.

###  3.4 H$\beta$ Equivalent Width

Figure 9: Same as Figure 1 but for the rest-frame H$\beta$ EW.

As evident in Figure 9, the average and maximum values of $f_{esc}^{LyC}$
increase with H$\beta$ EW. $f_{esc}^{LyC}$(UV) has a significant and strong
correlation with H$\beta$ EW, indicating at least a moderate relationship
between the two. The distribution of $f_{esc}^{LyC}$ appears to occupy two
distinct ranges in EW (see also Figure 10), one with high ($\gtrsim 150$ Å) EW
containing the most prodigious strong LCEs ($f_{esc}^{LyC}$$\gtrsim 0.1$) and
one with low ($\lesssim 150$ Å) EW containing additional strong LCEs and the
majority of weaker LCEs. Two thirds of the galaxies in our combined sample
with H$\beta$ EW $>150$ Å have detected LyC. Similarly, 28 of the 39 non-
detections appear concentrated at EWs below 150 Å.

Figure 10: Same as Figure 2 but for rest-frame EW(H$\beta$).

Quantifying these results, the fraction of galaxies that are strong LCEs
appears to increase with EW, which we demonstrate in Figure 10. This may
indicate a preference for young, strong bursts of star formation. As noted
above, though, strong LCEs have a wide range in EWs. Thus, the increase in
strong LCE fraction with EW demonstrates a dearth of weak LCEs and non-
emitters rather than an increase in the prevalence in LCEs at high EW. In
other words, high H$\beta$ EW indicates a strong LCE, but high EW galaxies do
not account for the full population of LCEs. Some scatter persists in LCE
fraction, particularly at lower EW, depending on the $f_{esc}^{LyC}$
indicator. Uncertainties in the LCE fraction are large enough to suggest any
change in LCE fraction between the various indicators may not be significant.
However, our results suggest a significant occurrence of LCEs for EWs below
150 Å, consistent with the $f_{esc}^{LyC}$ trends in Figure 9.

The prevalence of LCEs at high EW demonstrates that H$\beta$ EW is more
sensitive to star-formation history or sSFR than it is to optical depth.
Studies at higher redshifts find similar results in combined H$\beta$+[O III]
EWs (e.g., Castellano et al., 2017; Endsley et al., 2021; Saxena et al.,
2021), suggesting that this trend persists in earlier cosmological epochs. As
stellar mass and starburst age dominate the H$\beta$ EW, LCEs with high EWs
are likely young starbursts and/or low mass galaxies. The lack of weak LCEs
and non-emitters at high EWs may indicate that young ages and/or low masses
better facilitate high LyC escape fractions. However, a subset of between 6
and 10 strong LCEs have EWs below 150 Å, indicating that strong LCEs occupy a
wide range of ages, star-formation histories, and/or sSFRs. These low-EW LCEs
are likely more evolved galaxies with higher stellar mass and metallicity. We
discuss age and mass further in §4.

###  3.5 FUV Magnitude

Figure 11: Same as Figure 1 but for M1500.

The immense scatter in Figure 11 anticipates the weak, insignificant
correlation coefficients for observed $M_{1500}$ in all three $f_{esc}^{LyC}$
metrics. Strong leakers persist below $M_{1500}=-20.25$, implying a population
of strong LCEs at higher stellar mass or higher UV luminosity. Figure 12
illustrates this persistence of LCEs across a range of $M_{1500}$. The LCE
fraction is relatively constant with observed FUV magnitude with a slight
increase at the faint end of the distribution ($M_{1500}>-19.25$), especially
in the UV $f_{esc}^{LyC}$ where LCEs have $f_{esc}^{LyC}$$\gtrsim 0.02$ and
9/15 LCEs are strong LCEs. This mild increase may indicate that strong LCEs
are more likely fainter and/or less massive galaxies, particularly given the
slight trend in $f_{esc}^{LyC}$ upper bounds with $M_{1500}$.

Figure 12: Same as Figure 2 but for M1500.

The persistence of LCEs at higher UV luminosities could imply that LCEs
consist of multiple epochs of recent star formation and/or higher stellar
masses. Indeed, EW H$\beta$ declines with increasing UV luminosity such that,
for the UV $f_{esc}^{LyC}$, 5 of the 6 strong LCEs and 9 of the 12 weak LCEs
with $M_{1500}<-20.25$ have H$\beta$ EWs $<100$ Å. A major caveat to this
interpretation is the effect of dust on the observed UV magnitude. With the
LzLCS, we find $f_{esc}^{LyC}$ depends both on dust (Saldana-Lopez et al.,
2022) and age (§3.4) but less-so on mass (§3.7). Such effects complicate
direct interpretation of $M_{FUV}$ as an $f_{esc}^{LyC}$ diagnostic without
ancillary information. We discuss $M_{1500}$ in a high redshift context
further in §5.3.

###  3.6 UV $\beta_{1200}$

Figure 13: Same as Figure 1 but for the UV spectral index $\beta_{1200}$.

Flury et al. (2022) measured the UV slope $\beta_{1200}$ from the COS spectra
over a range of 1050-1350 Å. While we present results for this measurement
below, we note that $\beta_{1200}$ can be quite different from the commonly
used $\beta_{1500}$. Therefore, $\beta_{1200}$ provides only a relative sense
of the relationship between $\beta_{1500}$ and $f_{esc}^{LyC}$.

Figure 13 demonstrates that the highest $f_{esc}^{LyC}$ values occur for the
steepest values of $\beta_{1200}$. Indeed, the highest values of
$f_{esc}^{LyC}$ all occur where $\beta\leq-2$ while a distinct envelope in
$f_{esc}^{LyC}$ describes $f_{esc}^{LyC}$ values for $\beta\gtrsim-2$. Though
strong LCEs can persist across a wide range of $\beta\in[-3,-1]$, the UV slope
has a significant correlation coefficient for two of the three $f_{esc}^{LyC}$
indicators ($\tau\sim-0.25$, see Table 1). As implied by the correlation
coefficients and Figure 13, bluer galaxies tend to have higher
$f_{esc}^{LyC}$. The distribution of LCE fraction over UV $\beta_{1200}$ in
Figure 14 is equally informative. LCE fraction increases as the slope steepens
from $\sim 0$ at $\beta_{1200}>-1$ to nearly 0.6 at $\beta_{1200}\sim-2.5$,
suggesting that the strongest LCEs contain young stellar populations without
substantial dust attenuation.

Figure 14: Same as Figure 2 but for UV $\beta_{1200}$.

These trends in $\beta_{1200}$ make sense because early type stars produce the
majority of the LyC and low dust extinction better facilitates LyC escape.
However, some galaxies with $-2<\beta<-1$ are strong LCEs, with the Wang et
al. (2019) strong LCEs being even redder than -1. One or a combination of
scenarios could describe LCEs with $\beta_{1200}>-2$: (i) a patchy dust screen
that reddens the UV spectrum but still allows FUV and LyC photons to escape
(the so-called “picket fence” model, e.g., Heckman et al., 2001; Saldana-Lopez
et al., 2022) and/or (ii) multiple generations of stars that, through feedback
processes, facilitate LyC escape by clearing out gas while bolstering the
redder part of the UV spectrum with emission from older stellar populations
(see discussion in Micheva et al. 2017 of this effect in Mrk 71). We discuss
the relationship between $\beta_{1200}$, dust, and stellar population age
further in §4.1.

###  3.7 Stellar Mass

Figure 15: Same as Figure 1 but for stellar mass.

Simulations disagree as to whether dwarf or more massive star-forming galaxies
are the major source of LyC photons for reionization. Figure 15 suggests that
dwarf galaxies ($M_{\star}<10^{9}$ $M_{\odot}$) are more likely to be strong
LCEs. More than half (31/52) of the dwarf galaxies in the combined sample are
LCEs, 13 of which are strong LCEs.

With increasing mass, the number of strong LCEs and corresponding
$f_{esc}^{LyC}$(UV) values appear to decrease. Even so, galaxies with
$M_{\star}>10^{9}$ M☉ can still be prodigious LCEs, some even having
$f_{esc}^{LyC}$$>0.1$. The $\tau$ coefficient and corresponding $p$ value
suggest the correlation is neither strong nor significant for any
$f_{esc}^{LyC}$ metric. Furthermore, the other two $f_{esc}^{LyC}$ metrics do
not indicate much of an envelope in $f_{esc}^{LyC}$ with $M_{\star}$. As
evident in Figure 16, the LCE fraction distribution is also relatively flat
over the range of stellar masses. Like $M_{1500}$, the distribution of
$f_{esc}^{LyC}$ with respect to $M_{\star}$ suffers from substantial scatter.
Unlike $M_{1500}$, the stellar masses have relatively large uncertainties.
These large uncertainties may contribute to the apparent scatter and even to
the flat LCE fraction distribution.

Figure 16: Same as Figure 2 but for $M_{\star}$.

###  3.8 $\Sigma_{\rm SFR}$, ${\rm sSFR}$, and Half-light Radius

Figure 17: Same as Figure 1 but for UV half-light radius $r_{50}$,
$\Sigma_{\rm SFR}$, and ${\rm sSFR}$. Dashed line in $\Sigma_{\rm SFR}$ is the
$f_{esc}^{LyC}$$\propto\Sigma_{\rm SFR}^{0.4}$ relation from Naidu et al.
(2020).

$\Sigma_{\rm SFR}$ exhibits a threshold of $\Sigma_{\rm SFR}$$\sim 10$ M⊙ yr-1
kpc-2 above which nearly all strong LCEs appear, as evidenced by Figure 17. We
see scatter in $f_{esc}^{LyC}$ similar to other diagnostics examined here,
again implying the effects of orientation or variations in host galaxy
properties. For comparison, we also show $f_{esc}^{LyC}$ as it varies with
half-light radius in Figure 17. The prominent transition in $f_{esc}^{LyC}$ at
$r_{50}\approx 0.7$ kpc suggests the concentration of star formation dominates
$\Sigma_{\rm SFR}$ in strong LCEs.

This affinity of strong LCEs for high $\Sigma_{\rm SFR}$ and low $r_{50}$
suggests higher concentrations of star formation provide the feedback
necessary to clear LyC escape paths in the ISM. Consistent with this
interpretation, we find that the maximum observed value of $f_{esc}^{LyC}$
increases with ${\rm sSFR}$. Although less pronounced, this envelope in
$f_{esc}^{LyC}$ is similar to the $f_{esc}^{LyC}$$\propto$$\Sigma_{\rm
SFR}$0.4 relation described by (Naidu et al., 2020). These trends suggest
concentrated star formation may be a good indicator of $f_{esc}^{LyC}$ and
could be pivotal to understanding the origin of reionization.

Figure 18: Same as Figure 2 but for $r_{50}$, $\Sigma_{\rm SFR}$, and ${\rm
sSFR}$.

We show the LCE fraction with respect to both $\Sigma_{\rm SFR}$ and ${\rm
sSFR}$ in Figure 18. In both cases, LCE fraction clearly increases with the
concentration of star formation. The LCE fraction changes relatively little
with $\Sigma_{\rm SFR}$ until reaching $\approx 10$ M⊙ yr-1 kpc-2, at which
point LCE fraction increases from $\sim 10$% to 60% over half a decade in
$\Sigma_{\rm SFR}$. LCE fraction increases more gradually with ${\rm sSFR}$,
suggesting that the concentration of star formation is more important than the
effects of stellar mass.

Comparing $\Sigma_{\rm SFR}$, sSFR, and $\rm\Sigma_{sSFR}=$$\Sigma_{\rm
SFR}$$/M_{\star}$ highlights this point as $\Sigma_{\rm SFR}$ and
$\rm\Sigma_{sSFR}$ yield much higher $\tau$ values than sSFR. The $\tau$
values of $\Sigma_{\rm SFR}$ and $\rm\Sigma_{sSFR}$ are also quite comparable,
which demonstrates that factoring in stellar mass does not produce a better
diagnostic. This result comes with the caveat that the LzLCS stellar masses
have high uncertainties, which could mitigate any effects stellar mass may
have on the observed trends.

High $\Sigma_{\rm SFR}$, like those seen for GP galaxies, need not imply
unusually strong SNe feedback (e.g., Chisholm et al., 2017; Jaskot et al.,
2017). Because $\Sigma_{\rm SFR}$ is derived from H$\beta$, the apparent
envelope in $f_{esc}^{LyC}$ could instead arise from an increased compactness
of H II regions. Indeed, $\tau$ is higher and more significant for $r_{50}$
than for $\Sigma_{\rm SFR}$, sSFR, or $\rm\Sigma_{sSFR}$. Such concentrated
star forming regions could deplete the surrounding gas so that LyC photons can
leave the local H II region into the more diffuse host galaxy and potentially
escape into the IGM.

###  3.9 Metallicity

In Figure 19, we compare $f_{esc}^{LyC}$ to the gas phase metallicity $\rm
12+\log_{10}(O/H)$. The LCEs with $f_{esc}^{LyC}$$\gtrsim 10$% all reside
below $\rm 12+\log_{10}(O/H)$$\sim 8.1$. Much like the case of $M_{\star}$,
the strong LCEs with $f_{esc}^{LyC}$$\lesssim 10$% and the weak LCEs are more
evenly distributed from $\rm 12+\log_{10}(O/H)$$=7.7-8.5$. While substantial
scatter appears in all three diagnostics, roughly two thirds (28/39) of the
non-LCEs fall _above_ $\rm 12+\log_{10}(O/H)$$\sim 8.1$. With a correlation
coefficient of $\tau\sim-0.2$ at about 2$\sigma$ significance, $f_{esc}^{LyC}$
does not have a significantly strong correlation with $\rm 12+\log_{10}(O/H)$;
however, the preferences of the most prodigious LCEs for lower metallicities
and non-LCEs for higher metallicities is promising as a diagnostic.

Figure 19: Same as Figure 1 but for gas phase metallicity $\rm
12+\log_{10}(O/H)$.

The LCE fraction shown in Figure 20 also suggests $\rm 12+\log_{10}(O/H)$ is a
promising, albeit complicated, diagnostic. As with the $f_{esc}^{LyC}$ trends
shown in Figure 19, we see strong LCEs both above and below $\rm
12+\log_{10}(O/H)$$\sim 8.1$. Given the uncertainty in the LCE fraction, it is
unclear whether the lack of LCEs at intermediate values of $\rm
12+\log_{10}(O/H)$$\sim 8.1$ implies the existence of distinct LCE populations
separated by metallicity.

Figure 20: Same as Figure 2 but for $\rm 12+\log_{10}(O/H)$.

###  3.10 Summary of Correlations

Ten of the seventeen indirect $f_{esc}^{LyC}$ diagnostics considered above are
significantly strong ($|\tau|>0.2$ with $p<0.00135$, $3\sigma$ significance):
$f_{esc}^{Ly\alpha}$, EW $\rm Ly\alpha$, $v_{sep}$, $O_{32}$, EW H$\beta$,
$\beta_{1200}$, $M_{1500,int}$, NUV $r_{50}$, $\Sigma_{\rm SFR}$, and
$\rm\Sigma_{sSFR}$. The $f_{esc}^{Ly\alpha}$, EW $\rm Ly\alpha$,
$\beta_{1200}$, NUV $r_{50}$, and $\Sigma_{\rm SFR}$ diagnostics are
significantly strong across both $f_{esc}^{LyC}$ metrics. These five
diagnostics are thus the most compelling for identifying LCEs and inferring
$f_{esc}^{LyC}$. H$\beta$ EW and $M_{1500,int}$ are significantly strong
correlations only with the COS UV $f_{esc}^{LyC}$. In general, $f_{esc}^{LyC}$
derived from fitting the UV spectrum gives the strongest, most significant
correlations for nearly all the indirect diagnostics. The success of this
particular LyC escape indicator is particularly compelling because, as
discussed in §2, it is the most reliable and direct measure of
$f_{esc}^{LyC}$.

The indirect $f_{esc}^{LyC}$ diagnostics motivated by optical depth contain
the majority of the strong correlations. Likely owing to their unique
sensitivity to H I, the diagnostics based on $\rm
Ly\alpha$–$f_{esc}^{Ly\alpha}$, $\rm Ly\alpha$ EW, and profile peak
separation–have the strongest correlations as measured by Kendall’s $\tau$ for
censored data. However, the $v_{sep}$ correlations lack the significance of
the other $\rm Ly\alpha$ correlations, in part because of insufficient data.
While not as strongly correlated with $f_{esc}^{LyC}$ as the $\rm Ly\alpha$
properties, the $O_{32}$ flux ratio still exhibits a strong, significant
correlation. The [O I] flux relative to other emission lines’ is less
successful as an $f_{esc}^{LyC}$ diagnostic, likely owing to the faintness of
[O I].

Of the remaining diagnostics, which are motivated by the physical mechanism
behind LyC escape, the half-light radius and $\Sigma_{\rm SFR}$ have the
strongest, most significant correlations with $f_{esc}^{LyC}$. Indeed, strong
LCEs exhibit clear associations with highly concentrated star formation,
distinguishing which sorts of galaxies or regions within galaxies are emitting
LyC. Both H$\beta$ EW and UV $\beta_{1200}$ also correlate strongly and
significantly with $f_{esc}^{LyC}$. While stellar mass and attenuated UV
magnitude lack strong, significant correlations, these two diagnostics still
exhibit envelopes in $f_{esc}^{LyC}$ and suggest the strongest LCEs have
$M_{\star}<10^{9}$ M⊙. Such envelopes in $f_{esc}^{LyC}$ indicate that
multiple properties contribute to an ideal escape scenario, including
starburst age, ionization parameter, star formation rate, and orientation of
optically-thin channels.

##  4 Relationships Between Diagnostics

The results of these indirect diagnostics, particularly the envelopes in
$f_{esc}^{LyC}$, prompt further investigation of two-dimensional diagnostics
to determine trends in $f_{esc}^{LyC}$ with multiple parameters and explore
sources of the scatter apparent in many diagnostics. Below, we consider
relationships between diagnostics regarding stellar population age
($\beta_{1200}$ vs H$\beta$ EW, $O_{32}$ vs H$\beta$ EW), star formation
history ($O_{32}$ vs $M_{\star}$, $O_{32}$ vs $\rm 12+\log_{10}(O/H)$),
ionization structure ($O_{32}$ vs $O_{31}$), and ionization vs mechanical
feedback ($\Sigma_{\rm SFR}$ vs $O_{32}$). Given fewer systematic
uncertainties associated with the UV SED $f_{esc}^{LyC}$, we proceed with this
particular $f_{esc}^{LyC}$ metric in our subsequent analysis and discussion.

###  4.1 UV $\beta_{1200}$ vs H$\beta$ EW

Figure 21: UV $\beta_{1200}$ vs H$\beta$ EW for the LzLCS strong (squares) and
weak (dots) LCEs and non-emitters (triangles) color-coded by $f_{esc}^{LyC}$
derived from the UV spectrum. Solid lines are the Zackrisson et al. (2013)
$\sim 30$% solar metallicity picket fence model predictions for
$f_{esc}^{LyC}$$=0.0$ (magenta), $0.5$ (purple), and $0.9$ (cyan). There
appears to be no association of $f_{esc}^{LyC}$ with a particular sequence in
the UV $\beta_{1200}$ \- EW(H$\beta$) plane. While symbols are color-coded by
the UV-fit $f_{esc}^{LyC}$, the observed trends here and in subsequent figures
are similar if one of the other $f_{esc}^{LyC}$ metrics are used.

Although most galaxies with H$\beta$ EWs $<150$ Å do not appear associated
with LyC escape (§3.4), Zackrisson et al. (2013) propose the UV $\beta_{1200}$
– H$\beta$ EW plane could constrain $f_{esc}^{LyC}$. We compare their model
predictions for $f_{esc}^{LyC}$$=0.0$, $0.5$, and $0.9$ to the LzLCS results
in Figure 21 for $\sim 30$% solar-metallicity, picket fence nebulae over a
range of $10^{6}$ to $7\times 10^{8}$ year burst ages. Our combined results do
not agree with their predictions. We do see a mild gradient in $f_{esc}^{LyC}$
over UV $\beta_{1200}$ that appears to shift to steeper continua with
increasing EW, as their predictions suggest. However, the observed locus
appears shifted to higher spectral slopes: only one previously published LCE
falls into the $f_{esc}^{LyC}$$>0$ space predicted by Zackrisson et al.
(2013).

A major caveat to this apparent disagreement is the effect of dust, which may
substantially affect the predicted sequences in the UV $\beta_{1200}$ \-
H$\beta$ EW plane (Zackrisson et al., 2017). Moreover, the star formation
history also determines where a galaxy resides in Figure 21 as underlying
older stellar populations can strongly affect the optical continuum and thus
the H$\beta$ EW. LzLCS galaxies appear to be older and/or dustier than the
Zackrisson et al. (2013) models, which places our sample at lower H$\beta$ EW
and higher UV $\beta_{1200}$ than their predictions. Interpreting the
discrepancy comes with an additional complication that the UV $\beta$
predicted in the Zackrisson et al. (2013) is measured at redder wavelengths
than those accessible in the COS spectrum. Because an attenuated young stellar
SED peaks at $\sim 1100-1200$ Å, $\beta_{1200}$ is necessarily different from
the Zackrisson et al. (2013) UV $\beta$ predictions. That being said, we find
no distinguishable trends in $f_{esc}^{LyC}$ in the
$\beta_{1200}$-EW(H$\beta$) plane.

Whereas the Zackrisson et al. (2013) models suggest that strong LCEs should
have low H$\beta$ EWs and low UV $\beta_{1200}$, the strongest LCEs have high
($>100$ Å) H$\beta$ EWs while those LCEs with lower EWs have higher UV
$\beta_{1200}$. The latter group of LCEs may contain porous H II regions with
optically thin channels cleared by SN feedback or turbulence (see §5). In
particular, a two-stage starburst or dusty young starburst could be
responsible for the high UV $\beta_{1200}$, low H$\beta$ EW exhibited by some
of the strong LCEs.

###  4.2 $O_{32}$: Age, Mass, Metallicity

Figure 22: $O_{32}$ vs H$\beta$ EW. Symbols as in Figure 21. A strong
correlation is evident, owing to the coupling of ionization parameter with
starburst age.

Starburst age and ionization parameter are related since the earliest O stars
die off within the first 2-4 million years. Because these stars dominate the
LyC flux, properties sensitive to ionization parameter like $O_{32}$ will
consequently depend, at least in part, on the age of the starburst (e.g.,
Jaskot & Oey, 2013; Izotov et al., 2017) in addition to other properties.
Figure 22 demonstrates a tight relationship between $O_{32}$ and H$\beta$ EW.

The trend in Figure 22 could result from a variety of processes. Decreasing
H$\beta$ flux due to fewer ionizing photons would decrease the EW in tandem
with $O_{32}$ as the O+2 zone declines relative to O+. Due to the mass-
metallicity correlation for galaxies (e.g., Tremonti et al., 2004), increased
continuum flux due to higher galaxy mass would decrease the EW while
increasing the metallicity, the latter serving to decrease $O_{32}$. We find
that $O_{32}$ and EW decrease with increasing mass, but because mass
correlates with $\rm 12+\log_{10}(O/H)$, it is not clear whether mass
dominates this trend.

We find that LCEs occupy a wide range in $O_{32}$ and H$\beta$ EW. Roughly
half of the LCEs occur at $\log_{10}$$O_{32}$$>1$ or H$\beta$ EW$>100$ Å. Non-
emitters do not appear in this part of the diagram while all the
$f_{esc}^{LyC}$$>0.2$ leakers reside here. The remaining LCEs have lower
H$\beta$ EWs ($\approx 50-150$ Å) where many of the non-emitters also reside.
These separate sets of LCEs may imply different starburst ages or stellar mass
populations at which LyC escape occurs.

Escaping LyC photons should affect both H$\beta$ EW and $O_{32}$. The former
should decrease due to a drop in H$\beta$ flux while the latter should
increase due to a drop in [O II] flux. A lack of preference of $f_{esc}^{LyC}$
for any part of the distribution demonstrates that age and ionization
parameter mitigate any effects of escaping LyC photons on the measured
properties.

Figure 23: $O_{32}$ compared to metallicity indicator $\rm 12+\log_{10}(O/H)$
and SFH indicator $M_{\star}$. Symbols as in Figure 21.

Gas-phase and stellar metallicities also affect $O_{32}$. As collisionally
excited lines, [O III] and [O II] depend on electron temperature, which is
higher in low metallicity gas due to decreased cooling by metal emission
lines. Higher stellar metallicity changes the opacity of stellar interiors and
thus affects their photosphere temperatures. Secondarily, high stellar
metallicity increases line blanketing in the LyC of O and B stars. These
metallicity effects on stellar atmospheres serve to decrease the ionizing
photon budget. These dependencies on gas-phase and stellar metallicity
indicate that $O_{32}$ correlates with an abundance indicator like $\rm
12+\log_{10}(O/H)$, which we show in Figure 23.

While there is some scatter, a trend persists between $O_{32}$ and $\rm
12+\log_{10}(O/H)$ over nearly an entire decade. Although high $f_{esc}^{LyC}$
LCEs tend to have higher $O_{32}$ and lower $\rm 12+\log_{10}(O/H)$, LCEs
populate the majority of the distribution in Figure 23. Harder ionizing
spectra at lower metallicity could facilitate a density-bounded escape
scenario. Weaker feedback from low-metallicity stars may also generate a
dense, clumpy gas geometry that enables LyC escape (Jaskot et al., 2019),
perhaps owing to low covering fractions (e.g., Gazagnes et al., 2020; Saldana-
Lopez et al., 2022). Mechanical feedback from stars is stronger at higher
metallicity, which might explain the subpopulation of high LCEs at higher
metallicity. However, higher metallicity is also linked to higher UV
attenuation (e.g., Wang et al., 2019), thereby limiting LyC escape. The
reduced $f_{esc}^{LyC}$ of LCEs at higher metallicity may also be due to the
fact that, at higher galaxy masses, increased gravity may reduce the
effectiveness of feedback and thus limit LyC escape.

Following the known mass-metallicity relation, $O_{32}$ ought to anti-
correlate with $M_{\star}$ just as it does with $\rm 12+\log_{10}(O/H)$.
Indeed, in Figure 23, we see such a trend until $M_{\star}\lesssim 10^{9}$ M⊙
(corresponding to $O_{32}$$\approx 2-3$), at which point the two properties
appear to decouple. With weaker gravitational potentials, low $M_{\star}$
galaxies also undergo more bursty star formation (e.g., Weisz et al., 2012;
Guo et al., 2016). Such episodic star formation results in a range of stellar
feedback and ionization parameters, which would lead to a wide range in
$O_{32}$.

For a given metallicity or stellar mass, strong LCEs tend to have the highest
$O_{32}$ values. The leakers with $f_{esc}^{LyC}$$>0.2$ concentrate at low
$\rm 12+\log_{10}(O/H)$, low $M_{\star}$, and high $O_{32}$. Thus, the
strength of a starburst relative to the mass of a galaxy may determine its
likelihood of leaking LyC.

###  4.3 Oxygen Ionization Structure

Figure 24: $O_{32}$ vs $O_{31}$. Symbols as in Figure 21. Stasińska et al.
(2015) predictions for a 2 Myr starburst in a filled sphere (solid) and a
spherical shell (dashed) nebula bounded by radiation (magenta) and density
(purple). As a burst ages, $O_{32}$ decreases at a fixed $O_{31}$.

Photoionization models predict that $O_{31}$ should increase by 0.5 dex or
more as the LyC optical depth decreases from ionization to density boundary
conditions at a fixed $O_{32}$ (e.g., Stasińska et al., 2015; Plat et al.,
2019; Ramambason et al., 2020). We compare our results to the predicted
sequences for radiation- and density-bounded nebulae models from Stasińska et
al. (2015) in Figure 24. The combined samples reside on or above the
radiation-bounded sequence. We find no relationship between $f_{esc}^{LyC}$
and any offset from the locus. Because the strongest LCEs have only upper
limits in [O I], these galaxies may in fact be consistent with density-bounded
or picket fence scenarios. LCEs with detected [O I] reside along the
radiation-bounded model predictions, which may suggest that an isotropic
density-bounded escape scenario is not realized. Indeed, Ramambason et al.
(2020) predict that density-bounded channels within an otherwise optically-
thick medium can describe the distribution of $O_{32}$ and $O_{31}$ flux
ratios we observe for LCEs with $f_{esc}^{LyC}$$\lesssim 0.1$. The [O I] flux
from LCEs with high $O_{32}$ and disproportionately low $O_{31}$ may also be
contaminated by shocks, causing $O_{31}$ to shift to lower values than other
LCEs (Stasińska et al., 2015; Plat et al., 2019).

###  4.4 $\Sigma_{\rm SFR}$: Ionization & Age

Figure 25: $\Sigma_{\rm SFR}$ vs $O_{32}$. Symbols as in Figure 21.

$\Sigma_{\rm SFR}$ describes the concentration of O and B stars. As a result,
$\Sigma_{\rm SFR}$ influences the ionization parameter, which depends strongly
on the number density of early type stars. At least in the case of density-
bounded LyC escape, we anticipate higher $O_{32}$ for LCEs than for non-
emitters of the same $\Sigma_{\rm SFR}$ due to a deficit in [O II] flux.

We find LCEs prefer high $\Sigma_{\rm SFR}$ and high $O_{32}$, although only
one of the two is necessary to demonstrate the likelihood of LyC escape.
Figure 25 demonstrates that, with few exceptions, LCEs have high $O_{32}$,
high $\Sigma_{\rm SFR}$, or a combination thereof. The LCE population with
$\log_{10}\Sigma_{SFR}<1$ $M_{\odot}$ yr-1 kpc-2 is sparse. Thus, high
$O_{32}$ and high $\Sigma_{\rm SFR}$ are each indicative of LCEs, in
particular strong LCEs, with the most extreme LCEs having both high $O_{32}$
and high $\Sigma_{\rm SFR}$.

Combined with the results of UV $\beta_{1200}$ and H$\beta$ EW, the
combination of $\Sigma_{\rm SFR}$ and $O_{32}$ could suggest two different
types of strongly star-forming LCEs. One population is younger with high
$O_{32}$, high H$\beta$ EW, and low metallicity. The other population has
lower $O_{32}$, lower H$\beta$ EW, and higher metallicity, and resides in the
upper left of Figure 25. The latter suggests mechanical feedback, while the
former population may point to the additional role of ionization.

##  5 Discussion

###  5.1 Insight Into Lyman Continuum Escape

The substantial scatter evident in all considered diagnostics could be
attributed to one or both of two physical explanations. First, orientation
could cause the observed scatter since $f_{esc}^{LyC}$ likely depends on line
of sight in an anisotropic, radiation bounded escape scenario (e.g., Zastrow
et al., 2013; Cen & Kimm, 2015). Alternatively, the observed scatter might be
due to a time delay between the onset of star formation and the escape of the
LyC since feedback requires sufficient time to clear out escape channels, thus
decoupling properties from $f_{esc}^{LyC}$ to some extent (e.g., Trebitsch et
al., 2017; Barrow et al., 2020; Naidu et al., 2022). Whatever might cause this
scatter, it is likely physical because an upper envelope in $f_{esc}^{LyC}$
persists across several diagnostics. These envelopes suggest genuine trends do
exist between $f_{esc}^{LyC}$ and different parameters, but orientation, time
delays, and/or starburst properties like age and ionization parameter might
influence the values of $f_{esc}^{LyC}$, causing the observed scatter.

As implied by the envelopes and correlations in the $f_{esc}^{LyC}$
diagnostics, the strongest correlations are $\rm Ly\alpha$ EW,
$f_{esc}^{Ly\alpha}$, and $v_{sep}$ for the diagnostics motivated primarily by
optical depth. Other diagnostics with strong correlations include $O_{32}$,
H$\beta$ EW, $\Sigma_{\rm SFR}$, and $r_{50}$, which may give insights into
the physical mechanisms responsible for LyC escape. Some properties, notably
$O_{32}$ and H$\beta$ EW, exhibit a broad range in galaxies with
$f_{esc}^{LyC}$$>0.05$. Different LyC escape scenarios (e.g., Kakiichi &
Gronke, 2019; Gazagnes et al., 2020) could populate different parts of these
diagnostics because the corresponding conditions (e.g., optical depth,
porosity of the ISM, or feedback mechanism driving LyC escape) correspond to
different ranges of the observable parameter space of $O_{32}$ and H$\beta$
EW.

Investigation of two-dimensional diagnostics sheds light on these possible
explanations. At higher mass and higher metallicity, the lower $O_{32}$, lower
H$\beta$ EW LCEs likely clear out their birth cloud by mechanical feedback
from early type stars. This serves to more uniformly distribute the dense gas
and dust leftover from star formation into an optically thick bubble around
the O and B stars. SNe and turbulent feedback subsequently clear optically
thin channels in this bubble.

The younger, lower mass, lower metallicity nature of the high $O_{32}$, high
H$\beta$ EW LCEs indicates this second class of LyC leakers are like the GP
galaxies. Such galaxies, while strongly star forming, likely lack the high
velocity outflows necessary to expel the remnant birth cloud material and form
a uniform obscuring medium (Jaskot et al., 2019). However, they may have
sufficient radiative feedback to evacuate cavities in the ISM (Kakiichi &
Gronke, 2019). As a result, these GP LCEs might contain an ISM that is
optically thin except where dense clouds cause optically thick patches of
relatively small covering fractions. We speculate that the range of parameter
spaces occupied by LCEs indicates a shift in conditions and escape mechanism
from young GP starbursts with large optically thin cavities to older,
optically thick star forming galaxies with a porous ISM.

While these putative differences in LCEs spurred our comparison of
$f_{esc}^{LyC}$ proxies, the two-dimensional diagnostics also demonstrate how
different parameters depend on one another. These comparisons show that
metallicity, mass, $\Sigma_{\rm SFR}$, age, and ionization parameter are
inextricably linked. Furthermore, we find that age and ionization parameter
dominate any optical depth effects on properties like $O_{32}$ and H$\beta$
EW. In other words, young ages and strong radiation boost $O_{32}$ and
H$\beta$ EW in the strongest LCEs and may be critical for high levels of LyC
escape. It is also plausible that a feedback-induced time delay introduces LCE
“stragglers” into the distributions of $f_{esc}^{LyC}$ (e.g., Trebitsch et
al., 2017; Secunda et al., 2020).

###  5.2 Comparison with Simulations

Cosmological simulations provide a critical link between the low redshift LCEs
and high redshift surveys. Predictions from simulations connect parameters
such as stellar mass and $\Sigma_{\rm SFR}$ to $f_{esc}^{LyC}$ at the epoch of
reionization. Because LzLCs galaxies and other GP-like galaxies are reasonable
analogs of high-redshift galaxies, comparing our results to these predictions
from simulations can be useful both for translating indirect diagnostics to
high redshift and for interpreting the distributions we observe.

One of the most contentious issues regarding reionization is the mass regime
of galaxies that provide the ionizing photons. Some simulations predict that
$f_{esc}^{LyC}$ decreases with stellar mass (e.g., Kimm & Cen, 2014) and UV
luminosity (e.g., Trebitsch et al., 2017; Secunda et al., 2020) while others
predict that $f_{esc}^{LyC}$ increases with stellar mass and UV luminosity
(e.g., Ma et al., 2020). The strongest LCEs in the combined sample of local
LCEs are smaller, fainter galaxies with $M_{\star}\lesssim 10^{9}$ M⊙ and
$M_{FUV}\gtrsim-19$. Reionization is likely driven by galaxies with similar
properties.

The LzLCS samples the upper end of the simulated $M_{\star}$-$\Sigma_{\rm
SFR}$ distribution of LCEs in Naidu et al. (2020). We find that the combined
sample of LCEs is largely consistent with their predictions, including the
range in $f_{esc}^{LyC}$. However, our average $f_{esc}^{LyC}$ toward higher
stellar mass bins is about 0.1 dex lower than theirs even if considering only
the strongest LCEs. Moreover, although we find qualitative agreement between
the envelope in $f_{esc}^{LyC}$ with respect to $\Sigma_{\rm SFR}$ and the
$f_{esc}^{LyC}$$\propto\Sigma_{SFR}^{0.42}$ relation prescribed by Naidu et
al. (2020), the $M_{1500}<-20$ LCEs in the LzLCS are not strong leakers at the
same rate as their $M_{1500}>-19$ counterparts. While seemingly contrary to
predictions by Naidu et al. (2020), the fact that brighter, higher mass
galaxies are less likely to be prodigious LCEs in our sample may be consistent
with their simulations as intermediate-mass galaxies at $z>4$ tend to have
emission-line properties similar to those of the $z\sim 0.3$ low-mass galaxies
of this study. The $f_{esc}^{LyC}$-$\Sigma_{\rm SFR}$ envelope also agrees
with the sharp increase in $f_{esc}^{LyC}$ predicted by Sharma et al. (2017),
although we do not see the 20% $f_{esc}^{LyC}$ values they suggest for
$M_{\star}>10^{9}$ $M_{\odot}$. This might suggest that properties in addition
to $\Sigma_{\rm SFR}$ (e.g., metallicity, dust) affect LyC escape.

$O_{32}$ is a more empirical measurement than $M_{\star}$ or $\Sigma_{\rm
SFR}$ and directly traces the ionization parameter. Together, these attributes
have favored $O_{32}$ as an $f_{esc}^{LyC}$ diagnostic in the past.
Unfortunately, properties like $\Sigma_{\rm SFR}$ and metallicity can also
contribute to $O_{32}$, thus making interpretation of $O_{32}$ unclear (see
Nakajima & Ouchi, 2014). As with the LzLCS, recent studies have found
substantial scatter in the $f_{esc}^{LyC}$-$O_{32}$ diagnostic (e.g., Izotov
et al., 2018b; Bassett et al., 2019), which Bassett et al. (2019) attribute to
a combination of orientation and opening angle effects. With these physical
and empirical caveats in mind, recent cosmological hydrodynamic simulations
have predicted where high redshift LCEs should appear in the $O_{32}$-$R_{32}$
diagnostic (e.g., Katz et al., 2020). We show the distribution of LzLCS,
Izotov et al. (2016a, b, 2018a, 2018b, 2021); Wang et al. (2019), and
published LCEs (de Barros et al., 2016; Nakajima et al., 2020) and LAEs (Reddy
et al., 2021) at $z\sim 2-3$ in Figure 26 along with $z\in[0.2,0.4]$ star-
forming galaxies taken from the Thomas et al. (2013) SDSS catalog for
comparison. Simulations by Barrow et al. (2020) suggest the entire
distribution of LCEs shifts 0.5 dex lower in $R_{23}$ at high redshift due to
decreases in metallicity while maintaining high $O_{32}$ through high
ionization parameter. On the other hand, simulations by Katz et al. (2020)
suggest LCEs at high redshift shift to a locus $\sim 1$ dex higher in $O_{32}$
than that of the SDSS galaxies, intersecting the low-metallicity “tail” where
the combined low-redshift LCEs reside.

In contrast to these theoretical predictions, LCEs at moderate redshift
($z\sim 3$, Nakajima et al. 2020) still reside in this tail, despite having
fesc values in excess of 10%. Thus, even in earlier epochs, properties like
age, ionization parameter, and metallicity may continue to outweigh the
effects of LyC escape in setting $O_{32}$. Star-forming galaxies from the
MOSDEF survey at $z\sim 2-3$ have $O_{32}$ values more consistent with non-
LCEs and typical SDSS star-forming galaxies. These galaxies are likely
analogous to Lyman break galaxies (LBGs), suggesting that LBGs are only weakly
leaking, if at all. Reddy et al. (2021) find that these galaxies have older
burst ages (continuous star formation ages of $\sim 100$ Myr), moderate
ionization parameters ($\log U\sim-3$), and relatively higher metallicities
($\sim$30-40% solar), which could readily account for the observed $O_{32}$
values. These same properties may be associated with low levels of LyC escape.
The LzLCS results suggest that $O_{32}$ and $f_{esc}^{LyC}$ may be indirectly
related, as both $O_{32}$ and $f_{esc}^{LyC}$ may depend on the same physical
properties such as ionization parameter and metallicity. As such, these
underlying properties will predominantly determine the observed $O_{32}$ flux
ratio at $z\sim 3$ and at $z\sim 6$.

Figure 26: The $O_{32}$ vs $R_{23}$ diagnostic for the combined sample.
Symbols as in Figures 1-19 with the addition of LCEs from Nakajima et al.
(2020) with measured $O_{32}$ and $R_{23}$ (green hexagons), the $z=3.2$ LCE
from de Barros et al. (2016, blue triangle), and the $z=1.85-3.49$ LBGs from
Reddy et al. (2021, red crosses). For reference, we show SDSS star-forming
galaxies (grey dots) selected from the value-added catalog of Thomas et al.
(2013) with requirements that $z\in[0.2,0.4]$ and that strong emission lines
have $F/\sigma>3$. We show the locus of high-redshift LCEs predicted from
cosmological simulations by Katz et al. (2020, dash-dotted) and Barrow et al.
(2020, dashed).

###  5.3 Implications for High Redshift

One of the principal objectives of the LzLCS is to ascertain which diagnostics
are best suited to identifying LCEs and inferring their $f_{esc}^{LyC}$. While
$f_{esc}^{Ly\alpha}$, EW $\rm Ly\alpha$, $v_{sep}$, $O_{32}$, $\beta_{1200}$,
NUV $r_{50}$, and $\Sigma_{\rm SFR}$ all exhibit strong significant
correlations with $f_{esc}^{LyC}$ at $z\sim 0.3$, comparisons with
observations of these properties at higher redshifts are necessary to extend
our results to the epoch of reionization.

$\rm Ly\alpha$ offers one of the most promising indicators of LyC escape: both
$f_{esc}^{Ly\alpha}$ and $\rm Ly\alpha$ EW correlate well with
$f_{esc}^{LyC}$, at least in part because $\rm Ly\alpha$ is also sensitive to
the line-of-sight neutral hydrogen column density. However, this sensitivity
complicates the detection of $\rm Ly\alpha$ because the neutral hydrogen
fraction in the IGM increases dramatically with redshift (Tang et al., 2019;
Thomas et al., 2021), accounting for the steep drop in LAEs observed at
$z\gtrsim 6$ (e.g., Pentericci et al., 2011; Treu et al., 2013; Schenker et
al., 2014). Although $\rm Ly\alpha$ falls in the optical observing window at
redshifts of $z\sim 6$, few $\rm Ly\alpha$ emitters may be observable (an
effect exacerbated by the shift of $\rm Ly\alpha$ into the near infrared at
$z>7$). But it does suggest that LCEs are frequently LAEs. The few LAEs
detected at $z>6$ may be leaking LyC. Indeed, the size of ionized bubbles
necessary to produce observed double-peaked $\rm Ly\alpha$ at $z\sim 6$ imply
LyC escape (Meyer et al., 2021).

Star-forming galaxies observed at higher redshifts allow a more direct
comparison with the LCEs from our sample. The $z\sim 2-3$ LBG-like star-
forming galaxies from Reddy et al. (2021) have lower $f_{esc}^{Ly\alpha}$ and
EW $\rm Ly\alpha$ than the $z\sim 0.3$ LCEs, suggesting that LBGs are either
extremely weak LCEs or non-emitters. The $z\sim 3$ individual LCEs from
Nakajima et al. (2020) and $z\sim 2-3$ stacked spectra LCEs with
$f_{esc}^{LyC}$$\geq 0.2$ from Steidel et al. (2018), Pahl et al. (2021), and
Naidu et al. (2022) exhibit $\rm Ly\alpha$ EWs of $30-110$ Å and LyC escape
fractions slightly larger than those of $z\sim 0.3$ LCEs at similar EWs. As
shown in Figure 3, the Pahl et al. (2021) relation between $f_{esc}^{LyC}$ and
EW $\rm Ly\alpha$ (based on the $z\sim 3$ stacks from Steidel et al., 2018) is
qualitatively consistent with the $f_{esc}^{LyC}$-$\rm Ly\alpha$ envelope in
the LzLCS galaxies.

At higher redshifts, the Marchi et al. (2018) analysis of stacked spectra of
$z=3.5-4.3$ LAEs indicates that galaxies with $\rm Ly\alpha$ EWs $\geq 50$ Å
have relative $f_{esc}^{LyC}$ significantly higher than 0.05 while those with
$\rm Ly\alpha$ EWs $<50$ Å do not. The individual LyC detection in _Ion3_ at
$z=4$ exhibits quadruply-peaked $\rm Ly\alpha$ with an EW of $\sim 40$ Å and a
relative escape fraction of 0.6 (Vanzella et al., 2018). These $z\sim 2-4$
results are consistent with the fact that 22/26 of the LzLCS strong LCEs have
EW $\rm Ly\alpha$ $>40\ \AA$. Thus, EW $\rm Ly\alpha$ as a LCE identifier and
$f_{esc}^{LyC}$ diagnostic should be readily extensible to the epoch of
reionization.

Studies by Steidel et al. (2018) and Pahl et al. (2021) at $z\sim 3$ suggest
that fainter galaxies have higher $f_{esc}^{LyC}$. As shown in §3.5, we only
find weak correlations with $M_{1500}$, although we do see a possible
dependence of maximal $f_{esc}^{LyC}$ on $M_{1500}$. That being said, our
sample probes fainter magnitudes than studies at higher redshifts. Saldana-
Lopez et al. (2022) found that individual LCE detections at higher redshifts
are consistent with LzLCS results for $F_{\lambda\rm LyC}/F_{\lambda 1500}$
(see their Figure 16). LyC measurements using stacked spectra from Steidel et
al. (2018) suggest that a possible envelope in $F_{\lambda\rm LyC}/F_{\lambda
1500}$ persists out to $M_{1500}\sim-22$. Despite scatter, the increase in max
$F_{\lambda\rm LyC}/F_{\lambda 1500}$ with magnitude suggests that $M_{1500}$
may still be a useful tool to estimate galaxy contributions to reionization.
However, additional observations of UV-faint LCEs at $z\sim 3$ are necessary
to confirm that the dependence of maximal $F_{\lambda LyC}/F_{\lambda 1100}$
on $M_{1500}$ at $z\sim 3$ is similar to that observed in the combined LzLCS
at $z\sim 0.3$.

Several indirect diagnostics driven by physical mechanism also show promise at
high redshift. Star-forming galaxies at $z\sim 6-8$ have UV slopes and star
formation rates comparable to those of LCEs in this study (e.g., Bhatawdekar &
Conselice, 2021), which may suggest that the $\beta_{1200}$ and $\Sigma_{\rm
SFR}$ diagnostics at $z\sim 0.3$ can be used to identify LCEs and even infer
$f_{esc}^{LyC}$ at the epoch of reionization. Shibuya et al. (2015, 2019) find
that star-forming galaxies become increasingly compact, approaching UV half-
light radii of $0.5$ to $0.7$ kpc at $z\sim 6$. Such compactness is typical
for the $z\sim 0.3$ LCEs in this study and LAEs in other studies (e.g., Kim et
al., 2021). Moreover, this UV concentration corresponds to an average increase
in $\Sigma_{\rm SFR}$ up to values $\gtrsim 10$ M⊙ yr-1 kpc-2 at $z\sim 6$
(Shibuya et al., 2015), which is comparable to the $\Sigma_{\rm SFR}$ of
$z\sim 0.3$ LCEs. Such concentrated star formation at the epoch of
reionization may be indicative of LyC escape. Star-forming galaxies at $z\sim
5-6$ could have higher-velocity outflows than galaxies of the same star-
formation rates and stellar masses at $z\sim 0$ (Sugahara et al., 2019). Such
extreme outflows might indicate augmented stellar feedback, in turn improving
LyC escape indicators like $\Sigma_{\rm SFR}$, $M_{\star}$, and UV $\beta$ at
$z\sim 6$ relative to $z\sim 0.3$.

The effect of gravity on $f_{esc}^{LyC}$ diagnostics is not straightforward.
In §3.7, we found no correlation of $f_{esc}^{LyC}$ with $M_{\star}$. We did
find a significant trend between $f_{esc}^{LyC}$ and ${\rm sSFR}$ in §3.8,
which may indicate that strong feedback better facilitates LyC escape in
weaker gravitational potentials. At $z\sim 2.3$ (Sanders et al., 2016), $z\sim
3$ (Reddy et al., 2021), and $z=7$ (Endsley et al., 2021), galaxies tend to
have ${\rm sSFR}$ in the 1-10 Gyr-1 range, which corresponds to weaker LCEs
and non-emitters at $z\sim 0.3$. Taking the §3.8 results at face value, these
higher redshift galaxies might not be prodigious LCEs. However, the analysis
by Saxena et al. (2021) of individual LCEs at $z\sim 3$ suggests galaxies with
lower ${\rm sSFR}$ values of 1-10 Gyr-1 have $f_{esc}^{LyC}$ in excess of 10%
at higher redshift. No $z\sim 0.3$ LCEs with $f_{esc}^{LyC}$$>10$% have sSFR
$<10$ Gyr-1. This striking difference at $z\sim 3$ may point to an increase in
burst-enabled LyC escape with redshift.

Other properties may readily persist at higher redshifts (e.g., H$\beta$ EW,
[O III]$\lambda 5007$ EW, $O_{32}$, $R_{23}$, cf. Nakajima et al. 2020; Boyett
et al. 2021; Endsley et al. 2021; Reddy et al. 2021; Saxena et al. 2021) but
are difficult to observe at $z\gtrsim 6$ without space-based facilities like
JWST. However, as we enter the JWST era, indirect indicators of
$f_{esc}^{LyC}$ will prove immensely useful via the trends presented in §3.
Indeed, unlike observing the LyC, the reionization of the IGM is relatively
insensitive to the anisotropic escape of LyC photons. The broad range of LCE
properties may still complicate interpreting some diagnostics, as might the
interdependent nature of parameters, since it is not clear how either will
differ at $z\sim 6$.

LCEs at the epoch of reionization could be intrinsically unlike those in the
LzLCS. Dust and the corresponding attenuation laws might alter the optical
depth of the ISM at the epoch of reionization. Low metallicities in the early
universe could suppress fragmentation during star formation, resulting in a
top-heavy IMF with higher concentrations of O and B stars than at low redshift
for the same $\Sigma_{\rm SFR}$. Thus, a more physically complete
understanding of the ISM and star formation at high redshift are paramount to
successfully interpreting the signatures of LyC escape. Nevertheless, we find
that young ages, concentrated star formation, and high ionization parameter
are important for LyC escape.

##  6 Conclusion

Using the results for the 66 galaxies at $z\sim 0.3$ in the LzLCS, we have
conducted the first thorough statistical test of indirect diagnostics of
$f_{esc}^{LyC}$ and LyC escape. We find that $f_{esc}^{Ly\alpha}$, EW $\rm
Ly\alpha$, $v_{sep}$, $O_{32}$, EW H$\beta$, $\beta_{1200}$, $M_{1500,int}$,
NUV $r_{50}$, $\Sigma_{\rm SFR}$, and $\rm\Sigma_{sSFR}$ provide the best
correlations with $f_{esc}^{LyC}$. However, many of these diagnostics exhibit
an upper limit in the values of $f_{esc}^{LyC}$ that varies with the property
in question. We interpret these results as underlying trends that are obscured
by a combination of line-of-sight effects and variations in other physical
properties. Indirect diagnostics arising from stellar feedback mechanisms
(e.g., stellar winds, SNe) tend to suffer from this effect more so than those
based on optical depth. Such mechanism-based trends are still useful to infer
$f_{esc}^{LyC}$ while also providing insight into the nature of LyC escape. We
note that the apparent success of each diagnostic may depend on the
$f_{esc}^{LyC}$ metric used, with templates fit to the COS spectra giving the
strongest, most significant correlations of properties with $f_{esc}^{LyC}$.
Since many of the trends show significant scatter for any given diagnostic
parameter, we also perform a quantitative assessment of the fraction of
galaxies which are LCEs.

Our statistical tests include the following diagnostics: fraction of $\rm
Ly\alpha$ photons ($f_{esc}^{Ly\alpha}$) which escape the host galaxy, EW($\rm
Ly\alpha$), the velocity separation ($v_{sep}$) between the $\rm Ly\alpha$
profile peaks, the [O I]$\lambda 6300$/H$\beta$ flux ratio, the [O
III]$\lambda 5007$/[O I]$\lambda 6300=O_{31}$ flux ratio, the [O III]$\lambda
5007$/[O II]$\lambda 3726,3729=O_{32}$ flux ratio, EW(H$\beta$), $M_{1500}$,
UV power law slope $\beta_{1200}$, host galaxy stellar mass $M_{\star}$, half-
light radius $r_{50}$, $\Sigma_{\rm SFR}$, and sSFR=SFR$/M_{\star}$.

Of these diagnostics, $f_{esc}^{Ly\alpha}$, EW $\rm Ly\alpha$, $v_{sep}$,
$O_{32}$, $\beta_{1200}$, NUV $r_{50}$, and $\Sigma_{\rm SFR}$ exhibit some of
the strongest, most significant correlations with $f_{esc}^{LyC}$. Trends in
$f_{esc}^{LyC}$ with $\rm Ly\alpha$ properties like $f_{esc}^{Ly\alpha}$, EW
$\rm Ly\alpha$, and $v_{sep}$ suggest line-of-sight optical depth plays a key
role in identifying LCEs and inferring $f_{esc}^{LyC}$. Other diagnostics,
including $O_{32}$, $r_{50}$, $\Sigma_{\rm SFR}$, and $\beta_{1200}$, also
demonstrate strong, significant correlations with $f_{esc}^{LyC}$, indicating
concentrated star-formation, young stellar populations, and high ionization
play pivotal roles in LyC escape.

From these diagnostics, we have obtained deeper insight into the nature of
LCEs and the conditions that may best facilitate LyC escape. LCEs tend to have
compact ($r_{50}<0.5$ kpc) regions of intense star formation ($\Sigma_{\rm
SFR}$$\gtrsim 10$ M⊙ yr-1 kpc-2), indicating stellar feedback may be crucial
for LyC escape.

The distribution of $f_{esc}^{LyC}$ over other properties like ionization
parameter (measured by $O_{32}$), burst age (measured by H$\beta$ EW), stellar
mass, and metallicity suggests a broad range in LCE galaxy properties that
persists across many diagnostics. This diversity suggests that different
physical conditions can lead to LyC escape or that a time delay between star
formation and LyC escape can decouple $f_{esc}^{LyC}$ from certain galaxy
properties. With that being said, galaxies with high $O_{32}$, high H$\beta$
EW, and low $\beta$ are more likely to be strong LCEs ($f_{esc}^{LyC}$$\gtrsim
0.05$), suggesting that galaxies like the Green Peas are the most substantial
LyC leakers.

To better understand how certain properties can affect one another and the
inferred $f_{esc}^{LyC}$, we juxtapose $f_{esc}^{LyC}$ diagnostics. Zackrisson
et al. (2013) proposed a comparison of $\beta_{1200}$ with EW(H$\beta$). We
find that, while not following the predictions from Zackrisson et al. (2013),
LCEs with higher $f_{esc}^{LyC}$ tend to have higher EW(H$\beta$) ($>100$ Å)
and steeper UV continua ($<-2$). Moreover, LCEs distinguish themselves as
having either higher EW(H$\beta$) and steeper UV continua or lower
EW(H$\beta$) ($<100$ Å) and more shallow UV continua ($>-2$).

We also compare $O_{32}$ and EW(H$\beta$), finding the two are strongly
correlated. Through comparisons with mass and metallicity, we find higher
$O_{32}$ ($>5$), higher $f_{esc}^{LyC}$ LCEs have low masses and low
metallicities while the lower $O_{32}$ ($<5$), lower $f_{esc}^{LyC}$ LCEs have
higher masses and higher metallicities. These results point to a distinct
difference in LyC escape mechanism, as discussed above. The oxygen ionization
structure indicated by $O_{32}$ and $O_{31}$ (e.g., Stasińska et al., 2015)
demonstrates that LCEs cannot be readily described by a simple density-
bounded, isotropic escape scenario. Even so, all LCEs display concentrated
star formation, indicating that locally intense star formation is necessary
for LyC escape regardless of mechanism.

For further insight into LCEs, we compare our results with cosmological and
galaxy evolution simulation predictions. Consistent with some simulations
(e.g., Trebitsch et al., 2017; Secunda et al., 2020), we demonstrate that
compact, UV-faint, low-mass galaxies are far more likely to be significant LyC
leakers than other galaxies. While we do find that concentrated star formation
is a significant indicator of LyC escape, simulations that predict this
relationship suggest LCEs ought to have higher stellar masses than we find in
our sample (e.g., Sharma et al., 2017; Nakajima et al., 2020).

While we demonstrate that $\rm Ly\alpha$ is the best $f_{esc}^{LyC}$ indicator
at low redshift and that our LCEs are comparable to LyC-leaking LAEs at $z\sim
2-3$, the increasing IGM opacity at $z\gtrsim 4$ renders $\rm Ly\alpha$, at
best, difficult to detect from galaxies at the epoch of reionization. From a
holistic interpretation of other indirect diagnostics, a galaxy with high EW
$\rm Ly\alpha$, high $f_{esc}^{Ly\alpha}$, high $O_{32}$, high EW(H$\beta$),
low $\beta_{1200}$, low NUV $r_{50}$, and high $\Sigma_{\rm SFR}$ is likely an
LCE. Thus, observing rest-frame properties like $O_{32}$, EW(H$\beta$), and
$\Sigma_{\rm SFR}$ with _JWST_ is critical to identifying and understanding
the galaxies responsible for cosmic reionization.

Support for this work was provided by NASA through grant number _HST_
-GO-15626 from the Space Telescope Science Institute. Additional work was
based on observations made with the NASA/ESA Hubble Space Telescope, obtained
from the data archive at the Space Telescope Science Institute from _HST_
proposals 13744, 14635, 15341, and 15639. STScI is operated by the Association
of Universities for Research in Astronomy, Inc. under NASA contract NAS
5-26555. Funding for the Sloan Digital Sky Survey IV has been provided by the
Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science,
and the Participating Institutions. SDSS-IV acknowledges support and resources
from the Center for High Performance Computing at the University of Utah. The
SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research
Consortium for the Participating Institutions of the SDSS Collaboration. RA
acknowledges support from ANID Fondecyt Regular 1202007.

## References

  * Akritas & Siebert (1996) Akritas, M. G., & Siebert, J. 1996, MNRAS, 278, 919, doi: 10.1093/mnras/278.4.919
  * Alavi et al. (2020) Alavi, A., Colbert, J., Teplitz, H. I., et al. 2020, arXiv e-prints, arXiv:2007.05519.  https://arxiv.org/abs/2007.05519
  * Alexandroff et al. (2015) Alexandroff, R. M., Heckman, T. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 810, 104, doi: 10.1088/0004-637X/810/2/104
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., SipHocz, B. M., et al. 2018, aj, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bagley et al. (2017) Bagley, M. B., Scarlata, C., Henry, A., et al. 2017, ApJ, 837, 11, doi: 10.3847/1538-4357/837/1/11
  * Barrow et al. (2020) Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, arXiv e-prints, arXiv:2010.00592.  https://arxiv.org/abs/2010.00592
  * Bassett et al. (2019) Bassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2019, MNRAS, 483, 5223, doi: 10.1093/mnras/sty3320
  * Becker et al. (2001) Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850, doi: 10.1086/324231
  * Behrens et al. (2014) Behrens, C., Dijkstra, M., & Niemeyer, J. C. 2014, A&A, 563, A77, doi: 10.1051/0004-6361/201322949
  * Bergvall et al. (2013) Bergvall, N., Leitet, E., Zackrisson, E., & Marquart, T. 2013, A&A, 554, A38, doi: 10.1051/0004-6361/201118433
  * Bhatawdekar & Conselice (2021) Bhatawdekar, R., & Conselice, C. J. 2021, ApJ, 909, 144, doi: 10.3847/1538-4357/abdd3f
  * Binggeli et al. (2018) Binggeli, C., Zackrisson, E., Pelckmans, K., et al. 2018, MNRAS, 479, 368, doi: 10.1093/mnras/sty1061
  * Boyett et al. (2021) Boyett, K. N. K., Stark, D. P., Bunker, A. J., Tang, M., & Maseda, M. V. 2021, arXiv e-prints, arXiv:2110.15858.  https://arxiv.org/abs/2110.15858
  * Byrohl & Gronke (2020) Byrohl, C., & Gronke, M. 2020, arXiv e-prints, arXiv:2006.10041.  https://arxiv.org/abs/2006.10041
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Castellano et al. (2017) Castellano, M., Pentericci, L., Fontana, A., et al. 2017, ApJ, 839, 73, doi: 10.3847/1538-4357/aa696e
  * Cen & Kimm (2015) Cen, R., & Kimm, T. 2015, ApJ, 801, L25, doi: 10.1088/2041-8205/801/2/L25
  * Chisholm et al. (2017) Chisholm, J., Orlitová, I., Schaerer, D., et al. 2017, A&A, 605, A67, doi: 10.1051/0004-6361/201730610
  * Chisholm et al. (2019) Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182, doi: 10.3847/1538-4357/ab3104
  * Clarke & Oey (2002) Clarke, C., & Oey, M. S. 2002, MNRAS, 337, 1299, doi: 10.1046/j.1365-8711.2002.05976.x
  * de Barros et al. (2016) de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51, doi: 10.1051/0004-6361/201527046
  * Dijkstra et al. (2016) Dijkstra, M., Gronke, M., & Venkatesan, A. 2016, ApJ, 828, 71, doi: 10.3847/0004-637X/828/2/71
  * Endsley et al. (2021) Endsley, R., Stark, D. P., Charlot, S., et al. 2021, MNRAS, doi: 10.1093/mnras/stab432
  * Finkelstein et al. (2012) Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012, ApJ, 758, 93, doi: 10.1088/0004-637X/758/2/93
  * Fletcher et al. (2019) Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87, doi: 10.3847/1538-4357/ab2045
  * Flury et al. (2022) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, arXiv e-prints, arXiv:2201.11716.  https://arxiv.org/abs/2201.11716
  * Gazagnes et al. (2020) Gazagnes, S., Chisholm, J., Schaerer, D., Verhamme, A., & Izotov, Y. 2020, A&A, 639, A85, doi: 10.1051/0004-6361/202038096
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
  * Guo et al. (2016) Guo, Y., Rafelski, M., Faber, S. M., et al. 2016, ApJ, 833, 37, doi: 10.3847/1538-4357/833/1/37
  * Hayes et al. (2021) Hayes, M. J., Runnholm, A., Gronke, M., & Scarlata, C. 2021, ApJ, 908, 36, doi: 10.3847/1538-4357/abd246
  * Heckman et al. (2001) Heckman, T. M., Sembach, K. R., Meurer, G. R., et al. 2001, ApJ, 558, 56, doi: 10.1086/322475
  * Heckman et al. (2011) Heckman, T. M., Borthakur, S., Overzier, R., et al. 2011, ApJ, 730, 5, doi: 10.1088/0004-637X/730/1/5
  * Henry et al. (2015) Henry, A., Scarlata, C., Martin, C. L., & Erb, D. 2015, ApJ, 809, 19, doi: 10.1088/0004-637X/809/1/19
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Izotov et al. (2016a) Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016a, Nature, 529, 178, doi: 10.1038/nature16456
  * Izotov et al. (2016b) Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683, doi: 10.1093/mnras/stw1205
  * Izotov et al. (2018a) Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2018a, MNRAS, 474, 4514, doi: 10.1093/mnras/stx3115
  * Izotov et al. (2020) —. 2020, MNRAS, 491, 468, doi: 10.1093/mnras/stz3041
  * Izotov et al. (2017) Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2017, MNRAS, 471, 548, doi: 10.1093/mnras/stx1629
  * Izotov et al. (2021) Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2021, arXiv e-prints, arXiv:2103.01514.  https://arxiv.org/abs/2103.01514
  * Izotov et al. (2018b) —. 2018b, MNRAS, 478, 4851, doi: 10.1093/mnras/sty1378
  * Jaskot et al. (2019) Jaskot, A. E., Dowd, T., Oey, M. S., Scarlata, C., & McKinney, J. 2019, ApJ, 885, 96, doi: 10.3847/1538-4357/ab3d3b
  * Jaskot & Oey (2013) Jaskot, A. E., & Oey, M. S. 2013, ApJ, 766, 91, doi: 10.1088/0004-637X/766/2/91
  * Jaskot et al. (2017) Jaskot, A. E., Oey, M. S., Scarlata, C., & Dowd, T. 2017, ApJ, 851, L9, doi: 10.3847/2041-8213/aa9d83
  * Ji et al. (2020) Ji, Z., Giavalisco, M., Vanzella, E., et al. 2020, ApJ, 888, 109, doi: 10.3847/1538-4357/ab5fdc
  * Johnson et al. (2019) Johnson, B. D., Leja, J. L., Conroy, C., & Speagle, J. S. 2019, Prospector: Stellar population inference from spectra and SEDs.  http://ascl.net/1905.025
  * Kakiichi & Gronke (2019) Kakiichi, K., & Gronke, M. 2019, arXiv e-prints, arXiv:1905.02480.  https://arxiv.org/abs/1905.02480
  * Katz et al. (2020) Katz, H., Ďurovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164, doi: 10.1093/mnras/staa2355
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, doi: 10.1146/annurev-astro-081811-125610
  * Kim et al. (2021) Kim, K. J., Malhotra, S., Rhoads, J. E., & Yang, H. 2021, ApJ, 914, 2, doi: 10.3847/1538-4357/abf833
  * Kimm & Cen (2014) Kimm, T., & Cen, R. 2014, ApJ, 788, 121, doi: 10.1088/0004-637X/788/2/121
  * Kobulnicky & Kewley (2004) Kobulnicky, H. A., & Kewley, L. J. 2004, ApJ, 617, 240, doi: 10.1086/425299
  * Laursen et al. (2011) Laursen, P., Sommer-Larsen, J., & Razoumov, A. O. 2011, ApJ, 728, 52, doi: 10.1088/0004-637X/728/1/52
  * Leitherer et al. (2010) Leitherer, C., Ortiz Otálvaro, P. A., Bresolin, F., et al. 2010, ApJS, 189, 309, doi: 10.1088/0067-0049/189/2/309
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Ma et al. (2020) Ma, X., Quataert, E., Wetzel, A., et al. 2020, arXiv e-prints, arXiv:2003.05945.  https://arxiv.org/abs/2003.05945
  * Makan et al. (2020) Makan, K., Worseck, G., Davies, F. B., et al. 2020, arXiv e-prints, arXiv:2012.07876.  https://arxiv.org/abs/2012.07876
  * Marchi et al. (2018) Marchi, F., Pentericci, L., Guaita, L., et al. 2018, A&A, 614, A11, doi: 10.1051/0004-6361/201732133
  * Meyer et al. (2021) Meyer, R. A., Laporte, N., Ellis, R. S., Verhamme, A., & Garel, T. 2021, MNRAS, 500, 558, doi: 10.1093/mnras/staa3216
  * Micheva et al. (2017) Micheva, G., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 845, 165, doi: 10.3847/1538-4357/aa830b
  * Naidu et al. (2020) Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109, doi: 10.3847/1538-4357/ab7cc9
  * Naidu et al. (2022) Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2022, MNRAS, 510, 4582, doi: 10.1093/mnras/stab3601
  * Nakajima et al. (2020) Nakajima, K., Ellis, R. S., Robertson, B. E., Tang, M., & Stark, D. P. 2020, ApJ, 889, 161, doi: 10.3847/1538-4357/ab6604
  * Nakajima & Ouchi (2014) Nakajima, K., & Ouchi, M. 2014, MNRAS, 442, 900, doi: 10.1093/mnras/stu902
  * Osterbrock & Ferland (2006) Osterbrock, D., & Ferland, G. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (University Science Books)
  * P. Virtanen (2020) P. Virtanen, R. Gommers, T. O. 2020, Nature Methods, 17, 261–272, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Paardekooper et al. (2015) Paardekooper, J.-P., Khochfar, S., & Dalla Vecchia, C. 2015, MNRAS, 451, 2544, doi: 10.1093/mnras/stv1114
  * Pahl et al. (2020) Pahl, A. J., Shapley, A., Faisst, A. L., et al. 2020, MNRAS, 493, 3194, doi: 10.1093/mnras/staa355
  * Pahl et al. (2021) Pahl, A. J., Shapley, A., Steidel, C. C., Chen, Y., & Reddy, N. A. 2021, MNRAS, 505, 2447, doi: 10.1093/mnras/stab1374
  * Paoletti et al. (2020) Paoletti, D., Hazra, D. K., Finelli, F., & Smoot, G. F. 2020, J. Cosmology Astropart. Phys, 2020, 005, doi: 10.1088/1475-7516/2020/09/005
  * Pentericci et al. (2011) Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, ApJ, 743, 132, doi: 10.1088/0004-637X/743/2/132
  * Pilyugin et al. (2006) Pilyugin, L. S., Thuan, T. X., & Vílchez, J. M. 2006, MNRAS, 367, 1139, doi: 10.1111/j.1365-2966.2006.10033.x
  * Plat et al. (2019) Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978, doi: 10.1093/mnras/stz2616
  * Ramambason et al. (2020) Ramambason, L., Schaerer, D., Stasińska, G., et al. 2020, A&A, 644, A21, doi: 10.1051/0004-6361/202038634
  * Razoumov & Sommer-Larsen (2010) Razoumov, A. O., & Sommer-Larsen, J. 2010, ApJ, 710, 1239, doi: 10.1088/0004-637X/710/2/1239
  * Reddy et al. (2016) Reddy, N. A., Steidel, C. C., Pettini, M., & Bogosavljević, M. 2016, ApJ, 828, 107, doi: 10.3847/0004-637X/828/2/107
  * Reddy et al. (2021) Reddy, N. A., Topping, M. W., Shapley, A. E., et al. 2021, arXiv e-prints, arXiv:2108.05363.  https://arxiv.org/abs/2108.05363
  * Rivera-Thorsen et al. (2019) Rivera-Thorsen, T. E., Dahle, H., Chisholm, J., et al. 2019, Science, 366, 738, doi: 10.1126/science.aaw0978
  * Saldana-Lopez et al. (2022) Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2022, arXiv e-prints, arXiv:2201.11800.  https://arxiv.org/abs/2201.11800
  * Sanders et al. (2016) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23, doi: 10.3847/0004-637X/816/1/23
  * Sanders et al. (2020) Sanders, R. L., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 491, 1427, doi: 10.1093/mnras/stz3032
  * Sanders et al. (2021) Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19, doi: 10.3847/1538-4357/abf4c1
  * Sawant et al. (2021) Sawant, A. N., Pellegrini, E. W., Oey, M. S., López-Hernández, J., & Micheva, G. 2021, ApJ, submitted
  * Saxena et al. (2021) Saxena, A., Pentericci, L., Ellis, R. S., et al. 2021, MNRAS, doi: 10.1093/mnras/stab3728
  * Schenker et al. (2014) Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, ApJ, 795, 20, doi: 10.1088/0004-637X/795/1/20
  * Secunda et al. (2020) Secunda, A., Cen, R., Kimm, T., Gotberg, Y., & de Mink, S. E. 2020, arXiv e-prints, arXiv:2007.15012.  https://arxiv.org/abs/2007.15012
  * Sharma et al. (2016) Sharma, M., Theuns, T., Frenk, C., et al. 2016, MNRAS, 458, L94, doi: 10.1093/mnrasl/slw021
  * Sharma et al. (2017) —. 2017, MNRAS, 468, 2176, doi: 10.1093/mnras/stx578
  * Shibuya et al. (2015) Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15, doi: 10.1088/0067-0049/219/2/15
  * Shibuya et al. (2019) Shibuya, T., Ouchi, M., Harikane, Y., & Nakajima, K. 2019, ApJ, 871, 164, doi: 10.3847/1538-4357/aaf64b
  * Stasińska et al. (2015) Stasińska, G., Izotov, Y., Morisset, C., & Guseva, N. 2015, A&A, 576, A83, doi: 10.1051/0004-6361/201425389
  * Steidel et al. (2018) Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123, doi: 10.3847/1538-4357/aaed28
  * Storey & Hummer (1995) Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41, doi: 10.1093/mnras/272.1.41
  * Strom et al. (2018) Strom, A. L., Steidel, C. C., Rudie, G. C., Trainor, R. F., & Pettini, M. 2018, ApJ, 868, 117, doi: 10.3847/1538-4357/aae1a5
  * Sugahara et al. (2019) Sugahara, Y., Ouchi, M., Harikane, Y., et al. 2019, ApJ, 886, 29, doi: 10.3847/1538-4357/ab49fe
  * Tang et al. (2019) Tang, M., Stark, D. P., Chevallard, J., & Charlot, S. 2019, MNRAS, 489, 2572, doi: 10.1093/mnras/stz2236
  * Thomas et al. (2013) Thomas, D., Steele, O., Maraston, C., et al. 2013, MNRAS, 431, 1383, doi: 10.1093/mnras/stt261
  * Thomas et al. (2021) Thomas, R., Pentericci, L., Le Fèvre, O., et al. 2021, A&A, 650, A63, doi: 10.1051/0004-6361/202038438
  * Trebitsch et al. (2017) Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224, doi: 10.1093/mnras/stx1060
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Treu et al. (2013) Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJ, 775, L29, doi: 10.1088/2041-8205/775/1/L29
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vanzella et al. (2009) Vanzella, E., Giavalisco, M., Dickinson, M., et al. 2009, ApJ, 695, 1163, doi: 10.1088/0004-637X/695/2/1163
  * Vanzella et al. (2012) Vanzella, E., Guo, Y., Giavalisco, M., et al. 2012, ApJ, 751, 70, doi: 10.1088/0004-637X/751/1/70
  * Vanzella et al. (2018) Vanzella, E., Nonino, M., Cupani, G., et al. 2018, MNRAS, 476, L15, doi: 10.1093/mnrasl/sly023
  * Vanzella et al. (2020) Vanzella, E., Caminha, G. B., Calura, F., et al. 2020, MNRAS, 491, 1093, doi: 10.1093/mnras/stz2286
  * Verhamme et al. (2015) Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7, doi: 10.1051/0004-6361/201423978
  * Verhamme et al. (2017) Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13, doi: 10.1051/0004-6361/201629264
  * Wang et al. (2019) Wang, B., Heckman, T. M., Leitherer, C., et al. 2019, ApJ, 885, 57, doi: 10.3847/1538-4357/ab418f
  * Weisz et al. (2012) Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, ApJ, 744, 44, doi: 10.1088/0004-637X/744/1/44
  * Wise et al. (2014) Wise, J. H., Demchenko, V. G., Halicek, M. T., et al. 2014, MNRAS, 442, 2560, doi: 10.1093/mnras/stu979
  * Worseck et al. (2016) Worseck, G., Prochaska, J. X., Hennawi, J. F., & McQuinn, M. 2016, ApJ, 825, 144, doi: 10.3847/0004-637X/825/2/144
  * Worseck et al. (2014) Worseck, G., Prochaska, J. X., O’Meara, J. M., et al. 2014, MNRAS, 445, 1745, doi: 10.1093/mnras/stu1827
  * Wyithe & Loeb (2013) Wyithe, J. S. B., & Loeb, A. 2013, MNRAS, 428, 2741, doi: 10.1093/mnras/sts242
  * Yajima et al. (2018) Yajima, H., Sugimura, K., & Hasegawa, K. 2018, MNRAS, 477, 5406, doi: 10.1093/mnras/sty997
  * Yang et al. (2017a) Yang, H., Malhotra, S., Rhoads, J. E., et al. 2017a, ApJ, 838, 4, doi: 10.3847/1538-4357/aa6337
  * Yang et al. (2017b) Yang, H., Malhotra, S., Gronke, M., et al. 2017b, ApJ, 844, 171, doi: 10.3847/1538-4357/aa7d4d
  * Yang et al. (2020) Yang, J., Wang, F., Fan, X., et al. 2020, arXiv e-prints, arXiv:2009.13544.  https://arxiv.org/abs/2009.13544
  * Zackrisson et al. (2013) Zackrisson, E., Inoue, A. K., & Jensen, H. 2013, ApJ, 777, 39, doi: 10.1088/0004-637X/777/1/39
  * Zackrisson et al. (2017) Zackrisson, E., Binggeli, C., Finlator, K., et al. 2017, ApJ, 836, 78, doi: 10.3847/1538-4357/836/1/78
  * Zastrow et al. (2013) Zastrow, J., Oey, M. S., Veilleux, S., & McDonald, M. 2013, ApJ, 779, 76, doi: 10.1088/0004-637X/779/1/76