# Multivariate Predictors of LyC Escape II: Predicting LyC Escape Fractions
for High-Redshift Galaxies111Based on observations made with the NASA/ESA
Hubble Space Telescope, obtained at the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5-26555. These observations are
associated with programs GO-15626, GO-13744, GO-14635, GO-15341, and GO-15639.

Anne E. Jaskot Department of Astronomy, Williams College, Williamstown, MA
01267, USA  Anneliese C. Silveyra  Department of Astronomy, Williams College,
Williamstown, MA 01267, USA  Department of Physics, University of Nevada,
Reno, NV 89557, USA  Anna Plantinga Department of Mathematics & Statistics,
Williams College, Williamstown, MA 01267, USA  Sophia R. Flury Department of
Astronomy, University of Massachusetts Amherst, Amherst, MA 01002, USA
Matthew Hayes Department of Astronomy, Oskar Klein Centre, Stockholm
University, AlbaNova, SE-10691 Stockholm, Sweden  John Chisholm Department of
Astronomy, University of Texas at Austin, Austin, TX 78712, USA  Timothy
Heckman Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  Laura Pentericci INAF - Osservatorio Astronomico di
Roma, via Frascati 33, 00078, Monteporzio Catone, Italy  Daniel Schaerer
Observatoire de Genève, Université de Genève, Chemin Pegasi 51, 1290 Versoix,
Switzerland  Maxime Trebitsch Astronomy, Kapteyn Astronomical Institute,
Landleven 12, 9747 AD Groningen, The Netherlands  Anne Verhamme Observatoire
de Genève, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland
Univ. Lyon, Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, 69230 Saint-Genis-Laval, France  Cody Carr Center for
Cosmology and Computational Astrophysics, Institute for Advanced Study in
Physics, Zhejiang University, Hangzhou 310058, China  Institute of Astronomy,
School of Physics, Zhejiang University, Hangzhou 310058, China  Henry C.
Ferguson Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD
21218, USA  Zhiyuan Ji Steward Observatory, University of Arizona, Tucson, AZ
85721, USA  Mauro Giavalisco Department of Astronomy, University of
Massachusetts Amherst, Amherst, MA 01002, USA  Alaina Henry Space Telescope
Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Rui Marques-
Chaves Observatoire de Genève, Université de Genève, Chemin Pegasi 51, 1290
Versoix, Switzerland  Göran Östlin Department of Astronomy, Oskar Klein
Centre, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden  Alberto
Saldana-Lopez Department of Astronomy, Oskar Klein Centre, Stockholm
University, AlbaNova, SE-10691 Stockholm, Sweden  Claudia Scarlata Minnesota
Institute for Astrophysics, School of Physics and Astronomy, University of
Minnesota, 316 Church St. SE, Minneapolis, MN 55455, USA  Gábor Worseck
VDI/VDE Innovation+Technik, Berlin, Germany  Xinfeng Xu Center for
Interdisciplinary Exploration and Research in Astrophysics, Northwestern
University, Evanston, IL 60201, USA

###### Abstract

JWST is uncovering the properties of ever increasing numbers of galaxies at
$z>6$, during the epoch of reionization. Connecting these observed populations
to the process of reionization requires understanding how efficiently they
produce Lyman continuum (LyC) photons and what fraction ($f_{\rm esc}$) of
these photons escape into the intergalactic medium. By applying the Cox
proportional hazards model, a survival analysis technique, to the Low-redshift
Lyman Continuum Survey (LzLCS), we develop new, empirical, multivariate
predictions for $f_{\rm esc}$. The models developed from the LzLCS reproduce
the observed $f_{\rm esc}$ for $z\sim 3$ samples, which suggests that LyC
emitters may share similar properties at low and high redshift. Our best-
performing models for the $z\sim 3$ galaxies include information about dust
attenuation, ionization, and/or morphology. We then apply these models to
$z\gtrsim 6$ galaxies. For large photometric samples, we find a median
predicted $f_{\rm esc}$=0.047-0.14. For smaller spectroscopic samples, which
may include stronger emission line galaxies, we find that $\geq 33$% of the
galaxies have $f_{\rm esc}$ $>0.2$, and we identify several candidate extreme
leakers with $f_{\rm esc}$ $\geq 0.5$. The current samples show no strong
trend between predicted $f_{\rm esc}$ and UV magnitude, but limited
spectroscopic information makes this result uncertain. Multivariate
predictions can give significantly different results from single variable
predictions, and the predicted $f_{\rm esc}$ for high-redshift galaxies can
differ significantly depending on whether star formation rate surface density
or radius is used as a measure of galaxy morphology. We provide all parameters
necessary to predict $f_{\rm esc}$ for additional samples of high-redshift
galaxies using these models.

##  1 Introduction

Reionization represents a fundamental transformation of the universe’s
hydrogen gas and dramatically illustrates the effect galaxies can have on
their surroundings. The presence of Ly$\alpha$ absorption in quasar spectra
suggests that reionization ended sometime near $z\sim 6$ (e.g., Fan et al.,
2006; McGreer et al., 2015; Robertson, 2022). The electron scattering optical
depth of the cosmic microwave background is consistent with this general
picture, constraining the midpoint of reionization to $z\sim 7.8$ (Planck
Collaboration et al., 2020). However, the exact timeline of reionization
remains uncertain, in part because of the inhomogeneous nature of the
reionization process (e.g., Eilers et al., 2018; Jung et al., 2020; Becker et
al., 2021).

The recently launched James Webb Space Telescope (JWST) is rapidly expanding
our knowledge of galaxy properties during the epoch of reionization. Surveys
are detecting more bright galaxies than anticipated, challenging our
understanding of early galaxy evolution (e.g., Castellano et al., 2022; Naidu
et al., 2022; Finkelstein et al., 2023; Donnan et al., 2023; Harikane et al.,
2023; McLeod et al., 2024). JWST images suggest that galaxies in the
reionization era are morphologically compact (e.g., Robertson et al., 2023;
Ormerod et al., 2024; Morishita et al., 2023), and spectroscopic observations
are tracing the evolution of nebular ionization and metallicity to redshifts
well above 7 (e.g., Schaerer et al., 2022; Sanders et al., 2023; Tang et al.,
2023; Fujimoto et al., 2023; Curti et al., 2023; Backhaus et al., 2023).
Measurements of nebular emission lines in early galaxies are also constraining
their ionizing, Lyman continuum (LyC) photon production rate, $\xi_{\rm ion}$,
a key quantity needed to understand reionization. Such observations suggest
that faint galaxies with bursty star formation histories may have elevated
$\xi_{\rm ion}$ (e.g., Atek et al., 2023; Simmonds et al., 2024; Saxena et
al., 2023a).

Several complementary approaches can give insights as to which galaxy
populations dominated the reionization process. A galaxy’s contribution to
reionization depends on its $\xi_{\rm ion}$ and the fraction $f_{\rm esc}$ of
these LyC photons that escape into the intergalactic medium (IGM). Using
constraints on galaxy number densities at $z>6$, some studies explore
different distributions of $\xi_{\rm ion}$ and $f_{\rm esc}$ among the galaxy
population and seek to match the observed constraints on the reionization
timeline. Depending on the model assumptions, either moderately bright or very
faint galaxies may provide the majority of the ionizing photons (e.g.,
Finkelstein et al., 2019; Naidu et al., 2020).

Cosmological hydrodynamical simulations can also model the progress of
reionization over time and can explore how and why LyC photons escape
simulated galaxies. Simulations consistently find that feedback plays an
essential role in clearing out obscuring material (e.g., Wise & Cen, 2009; Cen
& Kimm, 2015; Paardekooper et al., 2015; Trebitsch et al., 2017). LyC input
may peak shortly following a burst of star formation after supernovae explode
(e.g., Ma et al., 2015; Trebitsch et al., 2017), although other factors such
as binary star evolution or an extended star formation history may affect this
timeline (e.g., Ma et al., 2015; Barrow et al., 2020; Katz et al., 2023). The
geometry and mechanism of LyC escape may also vary across the galaxy
population (e.g., Katz et al., 2023; Bremer & Dayal, 2023). Despite general
agreement on the importance of feedback, however, simulation predictions have
yet to be confirmed observationally and the galaxies that are the primary
drivers of reionization have not yet been conclusively identified.

Directly detecting LyC and constraining $f_{\rm esc}$ observationally becomes
difficult above $z>4$ (e.g., Inoue et al., 2014) due to high IGM attenuation.
Fortunately, LyC observations at lower redshifts can help test simulation
predictions regarding which galaxy properties regulate LyC escape. Over the
past decade, the number of LyC detections has grown rapidly, with dozens of
LyC Emitters (LCEs) now known at $z\sim 2-3$ (e.g., Mostardi et al., 2015;
Shapley et al., 2016; Vanzella et al., 2016, 2018; Bassett et al., 2019;
Fletcher et al., 2019; Rivera-Thorsen et al., 2019; Ji et al., 2020; Saxena et
al., 2022), multiple detections of LyC emission in stacked samples at $z\sim
2-4$ (e.g., Marchi et al., 2018; Steidel et al., 2018; Bian & Fan, 2020;
Nakajima et al., 2020), and more than 50 known LCEs at $z<0.5$ (e.g., Leitet
et al., 2011, 2013; Borthakur et al., 2014; Izotov et al., 2016b, 2018b, 2021;
Leitherer et al., 2016; Wang et al., 2019; Flury et al., 2022a). At both low
and intermediate redshift, LCEs appear deficient in absorbing gas and dust
(e.g., Gazagnes et al., 2018; Chisholm et al., 2018; Saldana-Lopez et al.,
2022; Steidel et al., 2018; Ji et al., 2020), morphologically compact (e.g.,
Borthakur et al., 2014; Izotov et al., 2018b; Flury et al., 2022b; Vanzella et
al., 2016; Marchi et al., 2018; Rivera-Thorsen et al., 2019), and bright in
higher ionization emission lines such as [O iii] $\lambda$5007 (e.g., Izotov
et al., 2018b; Flury et al., 2022b; Vanzella et al., 2016; Fletcher et al.,
2019; Nakajima et al., 2020).

Building on these observational and theoretical efforts, several studies have
proposed diagnostics to predict $f_{\rm esc}$ at $z>6$ based on observable
properties. Using simulated galaxies from the SPHINX cosmological radiation
hydrodynamics simulation, Choustikov et al. (2023) develop a method to predict
$f_{\rm esc}$ from a linear combination of observables, including the UV
slope, E(B-V), H$\beta$ luminosity, UV magnitude, and nebular line ratios.
Most other studies have taken an empirical approach, constructing $f_{\rm
esc}$ diagnostics based on the large observational sample of LCEs at $z\sim
0.3$ (e.g., Verhamme et al., 2017; Wang et al., 2021; Flury et al., 2022b; Xu
et al., 2023). The largest of these $z\sim 0.3$ samples is the Low-redshift
Lyman Continuum Survey (LzLCS; Flury et al. 2022a), a set of 66 galaxies with
LyC measurements from the Hubble Space Telescope (HST) and ancillary
ultraviolet and optical data from HST and the Sloan Digital Sky Survey (SDSS;
Blanton et al. 2017). Combined with archival datasets (Izotov et al., 2016a,
b, 2018a, 2018b, 2021; Wang et al., 2019), this $z\sim 0.3$ sample, hereafter
the LzLCS+, consists of 89 galaxies with measured LyC or stringent upper
limits. Based on an analysis of the LzLCS+, Chisholm et al. (2022) find that
the UV slope $\beta_{\rm 1550}$ shows one of the strongest correlations with
$f_{\rm esc}$ and propose that this single variable can serve as a predictor
of $f_{\rm esc}$ at high redshift. Other recent studies of the LzLCS+ consider
$f_{\rm esc}$ diagnostics that incorporate information from multiple
variables. Saldana-Lopez et al. (2022) generate an equation to predict $f_{\rm
esc}$ from E(B-V) and low-ionization UV absorption lines. By tracing the gas
and dust that destroy LyC photons, these parameters closely track $f_{\rm
esc}$. Employing this method at high redshift may be a challenge, however, as
measuring weak absorption lines requires high signal-to-noise observations of
galaxy continua (but see Saldana-Lopez et al. 2023). Emission lines and
photometry offer a simpler, if less direct, means of predicting $f_{\rm esc}$.
Mascia et al. (2023a) propose a new $f_{\rm esc}$ diagnostic using $\beta_{\rm
1550}$, [O iii] $\lambda$5007/[O ii] $\lambda$3727=O32, and half-light radius,
three of the variables that correlate strongly with $f_{\rm esc}$ in the
LzLCS+. By applying this diagnostic at high redshift using JWST observations
from the GLASS and CEERS surveys, they find predicted $f_{\rm esc}$ of $\sim
0.1$ for galaxies at $z>6$ (Mascia et al., 2023a, b). With a similar
combination of variables ($\beta_{\rm 1550}$, O32, and UV magnitude), Lin et
al. (2024) develop a regression model for the probability of LyC escape. Their
model suggests that $f_{\rm esc}$ may be high ($\sim 0.2$) in brighter
galaxies in the epoch of reionization (Lin et al., 2024). These studies all
agree that plausible LCE candidates exist at high redshift, although they
differ in their $f_{\rm esc}$ prediction methods.

In Jaskot et al. (2024), hereafter Paper I, we use the LzLCS+ to develop new
empirical multivariate models for predicting $f_{\rm esc}$. Because the LzLCS+
contains both LyC detections and upper limits, we adopt the statistical
techniques of survival analysis, which are suitable for such censored data.
Specifically, we employ the semiparametric Cox proportional hazards model to
generate $f_{\rm esc}$ predictions based on a desired set of input
observables. We show that a model limited to observables accessible at $z>6$
can reproduce the observed $f_{\rm esc}$ in the LzLCS+ with a root-mean-square
scatter of 0.46 dex. Of these observables, three variables (O32, $\beta_{\rm
1550}$, and the star formation rate surface density) are statistically
significant in the fit, and a model limited to these three input observables
predicts $f_{\rm esc}$ as well as the full model. The Cox model technique can
be customized to include any combination of variables that are available for
most of the LzLCS+ galaxies and hence offers a flexible tool for predicting
$f_{\rm esc}$ at high redshift.

In this paper, we apply Cox models for $f_{\rm esc}$ to samples of high-
redshift galaxies. Following the techniques in Paper I, which we summarize in
§2, we generate new models for the sets of variables in published high-
redshift samples. In §3, we first test the models’ performance at $z\sim 3$
using samples with published LyC measurements. With published samples at
$z\gtrsim 6$, we then generate $f_{\rm esc}$ predictions for galaxies in the
epoch of reionization in §4. In §5, we compare our models with alternative
proposed $f_{\rm esc}$ prediction methods from the literature, and we discuss
the implications of our results for studies of reionization in §6. We
summarize our conclusions in §7. In the Appendix, we provide parameters for
all models as well as examples on how to apply these models to future samples.
We adopt a cosmology of $H_{0}=70$ km s-1 Mpc-1, $\Omega_{m}=0.3$, and
$\Omega_{\Lambda}=0.7$.

##  2 Methods

###  2.1 Sample: The Low-redshift Lyman Continuum Survey

We derive our empirical predictions of $f_{\rm esc}$ from the LzLCS+, a
combined, homogeneously processed dataset consisting of the Low-redshift Lyman
Continuum Survey (Flury et al., 2022a) and archival samples with HST Cosmic
Origins Spectrograph (COS) LyC observations (Izotov et al., 2016a, b, 2018a,
2018b, 2021; Wang et al., 2019). A full description of the data processing and
UV and optical measurements appears in Flury et al. (2022a), Saldana-Lopez et
al. (2022), and Paper I, but we summarize key points here.

The LzLCS+ contains 89 galaxies at $z\sim 0.3$, a distance where the LyC
redshifts into a sensitive wavelength range for the COS detector. The LzLCS
targets were selected on properties proposed to indicate LyC escape: high O32
ratios (O32 $\geq 3$), high star formation rate surface densities
($\Sigma_{\rm SFR}$$>0.1$ M☉ yr-1 kpc-2 ), and/or blue UV slopes ($\beta<-2$).
The LzLCS+ covers a wide range of parameter space, spanning $\sim 2$ dex in
O32, $\sim 2$ dex in $\Sigma_{\rm SFR}$, and an observed (not dust-corrected)
UV absolute magnitude range of $M_{\rm 1500}=-18.3$ to -21.5. In Paper I and
this study, we exclude one galaxy from the LzLCS+ sample: J1333+6246 (Izotov
et al., 2016b). This galaxy has visibly truncated emission lines in its SDSS
spectrum and unphysical Balmer line ratios, which together suggest that its
nebular line flux measurements may be inaccurate. For this work, our total
sample therefore includes 88 $z\sim 0.3$ galaxies, 49 of which have detected
LyC.

To measure the LyC, we use COS G140L observations, processed using the calcos
pipeline (v3.3.9) and the FaintCOS software routines (Worseck et al., 2016;
Makan et al., 2021). We measure the LyC in a 20 Å-wide wavelength bin near
rest-frame 900 Å, while excluding any wavelengths above 1180 Å in the observed
frame because of telluric emission. We follow the definitions of LyC
detections and upper limits from Flury et al. (2022a), where detections are
observations with a probability $<0.02275$ of originating from background
counts and the upper limit for non-detections represents the the 84th
percentile of the background count distribution. To correct for Milky Way
attenuation, we adopt the Green et al. (2018) dust maps and Fitzpatrick (1999)
attenuation law.

Flury et al. (2022a) investigate several different measures of $f_{\rm esc}$.
Here, as in Paper I, we adopt the absolute $f_{\rm esc}$, with Starburst99
(Leitherer et al., 2011, 2014) spectral energy distribution (SED) fits to the
UV continuum providing the estimate of the initial intrinsic LyC (Chisholm et
al., 2019; Saldana-Lopez et al., 2022). Alternative estimates of the intrinsic
LyC from H$\beta$ neglect LyC photons absorbed by dust, assume an isotropic
geometry, and require an assumed simplistic star formation history (e.g.,
Flury et al., 2022a). Another alternative measure of LyC escape is the
$F_{\lambda{\rm LyC}}/F_{\lambda{\rm 1100}}$ flux ratio. As a ratio of two
observed fluxes, this method requires no model assumptions. However, we found
in Paper I that this quantity was more difficult to predict from easily
accessible observables with the survival analysis models. Hence, for the
$f_{\rm esc}$ models in this paper, we proceed with the absolute $f_{\rm esc}$
estimates from the UV SED fits for the LzLCS+.

Optical and UV observations supply information about numerous other properties
of the LzLCS+ galaxies (see Flury et al. 2022a and Saldana-Lopez et al. (2022)
for details). We obtain stellar mass ($M_{*}$) estimates from Prospector (Leja
et al., 2017; Johnson et al., 2019) fits to the SDSS and Galaxy Evolution
Explorer (GALEX; Martin et al. 2003) photometry (Flury et al. 2022a, Ji et al.
in prep.). With multi-Gaussian fits to nebular lines, we measure nebular line
fluxes and equivalent widths (EWs), and we estimate the nebular dust
attenuation, E(B-V)neb, from Balmer line ratios. We calculate the oxygen
abundance using the direct method and the pyneb package (Luridiana et al.,
2015), adopting $n_{e}=100$ cm-3 and the estimated [O iii] $\lambda$4363 flux
from the Pilyugin et al. (2006) “ff-relation” in cases where the [S ii]
$\lambda\lambda$6716,6731 doublet or [O iii] $\lambda$4363 auroral line are
undetected. We estimate star formation rates (SFRs) from the Kennicutt & Evans
(2012) SFR calibration using the dust-corrected H$\beta$ luminosities and Case
B H$\alpha$/H$\beta$ ratio (Storey & Hummer, 1995). The COS near-UV
acquisition images allow us to measure the UV half-light radius $r_{\rm
50,NUV}$, which we then use to calculate $\Sigma_{\rm SFR}$ as

$\Sigma_{\rm SFR}=\frac{\rm SFR}{2\pi r^{2}_{\rm 50,NUV}}.$ (1)

In addition to $f_{\rm esc}$, we derive several other parameters from the HST
COS UV spectra. We obtain estimates of the dust attenuation E(B-V) from the UV
spectrum Starburst99 SED fits; we label this parameter as E(B-V)UV to
distinguish it from the nebular dust attenuation E(B-V)neb derived from the
Balmer lines. Although the G140L spectra do not extend to rest-frame 1500 Å,
we can estimate the “observed” (non-extinction corrected) absolute magnitude
at 1500 Å ($M_{1500}$) and the power law index slope at 1550 Å
($\beta_{1550}$) by extrapolating the SED fits to longer wavelengths (Saldana-
Lopez et al., 2022). The inferred $\beta_{\rm 1550}$ values do match the
observed values for the few galaxies with existing longer wavelength UV
spectra (Chisholm et al., 2022). We also measure Ly$\alpha$ properties from
the COS spectra. We derive Ly$\alpha$ fluxes and EWs by linearly fitting the
continuum within 100 Å of Ly$\alpha$, excluding regions affected by nebular or
stellar features, and integrating all emission above the continuum. We then
derive the Ly$\alpha$ escape fraction ($f_{{\rm esc,Ly}\alpha}$) using the
dust-corrected H$\beta$ flux and intrinsic Case B Ly$\alpha$/H$\beta$ ratio
(Storey & Hummer, 1995) appropriate for the galaxies’ measured electron
temperatures and densities. We note that the Ly$\alpha$ measurements represent
the net sum of both underlying absorption along the line of sight and
scattered Ly$\alpha$ emission within the aperture. Nine galaxies have
detectable Ly$\alpha$ absorption troughs that overlap with the Si iii
$\lambda$1206 absorption feature. We increase our uncertainties to account for
the change in flux from omitting wavelengths within 500 km s-1 of the Si iii
line. The inclusion or exclusion of this region has only minor effects,
changing the Ly$\alpha$ EW by $<3$Å and with typical changes to $f_{{\rm
esc,Ly}\alpha}$ of only 0.001.

###  2.2 Multivariate Survival Analysis: The Cox Proportional Hazards Model

In order to generate multivariate diagnostics for $f_{\rm esc}$, we need to
incorporate information from both the $f_{\rm esc}$ detections and the upper
limits. Within the field of statistics, survival analysis techniques are
appropriate for censored datasets that contain limits. One such survival
analysis method is the Cox proportional hazards model (Cox 1972; see Clark et
al. 2003 and Bradburn et al. 2003 for reviews and Feigelson & Nelson 1985 and
Isobe et al. 1986 for examples in astronomy). Here we describe the basic form
of the Cox model and its assumptions. We refer the reader to Paper I for a
more thorough discussion of this model and its application to the LzLCS+. To
implement the Cox model, we use the CoxPHFitter routine in the lifelines
python package (Davidson-Pilon, 2019).

The Cox proportional hazards model predicts the probability of a particular
$f_{\rm esc}$ value given a set of input variables. Like many other survival
analysis techniques, implementations of the Cox model typically assume the
dataset contains measurements and lower limits (so-called “right-censored”
data). In contrast, our dependent variable data consist of $f_{\rm esc}$
measurements and associated upper limits. Consequently, we transform our
$f_{\rm esc}$ values to the absorbed fraction of LyC, $f_{\rm abs}$ = 1 -
$f_{\rm esc}$, for use in the Cox model. However, for ease of interpretation,
we put the $f_{\rm abs}$ results back into the form of $f_{\rm esc}$ in all
figures and in the tabulated results in the Appendix.

As applied to the LzLCS+, the Cox proportional hazards regression model fits
for the probability of a LyC detection in an infinitesimally small increment
of $f_{\rm abs}$, for a given set of independent variables and assuming no
detection at a lower value of $f_{\rm abs}$ (higher value of $f_{\rm esc}$).
The Cox model assumes a particular functional form for this $f_{\rm esc}$
probability, which is known as the “hazard function”:

$h(f_{\rm abs}|x)=h_{0}(f_{\rm
abs})\exp[\sum_{i=1}^{n}b_{i}(x_{i}-{\bar{x}_{i}})].$ (2)

In this equation, $b_{i}$ are the best-fit coefficients for each input
variable $x_{i}$, and $\bar{x}_{i}$ is the mean value of each input variable
in the reference LzLCS+ dataset. The term $h_{0}$($f_{\rm abs}$) is the
baseline hazard function, the probability of having $f_{\rm abs}$ in the case
where all input variables match their average values within the LZLCS+. The
Cox model is semi-parametric, in that the dependence on the input variables
has a fixed exponential functional form, but the baseline hazard function,
$h_{0}$, is estimated non-parametrically. Despite its fixed functional form
for the input variable dependence, the Cox model allows the input variables to
take any form. For instance, we could define an input variable $x_{i}$ as a
measured value $j$, the logarithm of that measurement $\log_{\rm 10}$($j$),
the square of that measurement $j^{2}$, or any other functional form of our
choosing. In our case, we opt to use logarithmic forms of our input variables
where possible (e.g, $\log_{\rm 10}$(O32), $M_{\rm 1500}$, $\log_{\rm
10}$($M_{*}$)). In this logarithmic form, most variables have a similar order
of magnitude and scale in the same manner, and an order of magnitude increase
in a particular input variable simply translates into increasing the
probability $h(f_{\rm abs}|x)$ by a factor of $e^{b_{i}}$. In this paper, the
only variables that we do not use in a logarithmic fashion are EW(Ly$\alpha$)
and $f_{{\rm esc,Ly}\alpha}$, because they range from negative (Ly$\alpha$
absorption) to positive (Ly$\alpha$ emission) within the LzLCS+.

Although Equation 2 gives the probability of $f_{\rm abs}$, we would like to
predict the expected values of $f_{\rm abs}$ and $f_{\rm esc}$ and their
associated uncertainty. We adopt the median of the probability distribution as
this expected value (e.g., Bradburn et al., 2003; Davidson-Pilon, 2019); the
model predicts that $f_{\rm esc}$ will be above this value 50% of the time and
below it 50% of the time. To determine this median $f_{\rm abs}$
mathematically, we first calculate the survival function, $S$($f_{\rm abs}$),
the probability that we do not detect LyC at $f_{\rm abs,detect}<$ $f_{\rm
abs}$ (i.e., at $f_{\rm esc,detect}>$ $f_{\rm esc}$). In the Cox model, the
survival function is

$S(f_{\rm abs})=\exp[-{\rm HF}_{0}(f_{\rm abs})\cdot{\rm ph}(x)],$ (3)

where HF0 is the baseline cumulative hazard function

${\rm HF}_{0}(f_{\rm abs})=\int_{0}^{f_{\rm abs}}h_{0}(f)df$ (4)

and ph($x$) is the partial hazards function, which describes how the
probability scales with the set of input variables $x$:

${\rm ph}(x)=\exp[\sum_{i=1}^{n}b_{i}(x_{i}-{\bar{x}_{i}})].$ (5)

(e.g., Cox, 1972; Bradburn et al., 2003; Davidson-Pilon, 2019; McLernon et
al., 2023). The median $f_{\rm abs}$, which corresponds to our predicted
$f_{\rm esc}$ value, is the $f_{\rm abs}$ value where $S$($f_{\rm
abs}$)$=0.5$. In some cases, $S$($f_{\rm abs}$)$>0.5$, even for our largest
tabulated value of $f_{\rm abs}$, which implies that $f_{\rm abs}$ is closer
to 1 than we can determine and $f_{\rm esc}$ $\sim 0$. In this circumstance,
we infer that $f_{\rm esc}$ is arbitrarily small and report a predicted
$f_{\rm esc}$ $=0$.

The lifelines CoxPhFitter returns the best-fit coefficients $b_{i}$ and the
cumulative baseline hazard HF0($f_{\rm abs}$) for each value of $f_{\rm abs}$
corresponding to a LyC detection in the LzLCS+ (see Paper I for more details
about the lifelines methodology). In the Appendix, we provide these best-fit
parameters for the models in Paper I and the models in this paper, and we give
examples of how to use these models to calculate the predicted $f_{\rm esc}$
for a set of observed input variables. The coefficients and hazards can be
used to predict the expected $f_{\rm esc}$ for any galaxy, as long as it has
estimated values for each of the independent variables used in the model.

####  2.2.1 Uncertainty in the Predicted $f_{\rm esc}$

Because the survival function represents a probability distribution, we can
also use it to calculate the expected uncertainty in our predicted $f_{\rm
esc}$ estimate. The $f_{\rm abs}$ values where $S$($f_{\rm abs}$) = 0.159 and
$S$($f_{\rm abs}$) = 0.841 represent the range in the predicted $f_{\rm esc}$
corresponding to the Normal theory 1-$\sigma$ uncertainty. This uncertainty
range reflects the inherent scatter of the relationships between $f_{\rm esc}$
and the input variables, where this scatter can come both from measurement
uncertainties and from genuine variation among the galaxy population. In Paper
I, we tested the effect of measurement uncertainty by performing a Monte Carlo
(MC) resampling of each independent and dependent variable according to its
observational uncertainty. We found that the distribution of predicted $f_{\rm
esc}$ from the resampled inputs was nearly always smaller than the 1-$\sigma$
bounds inferred from the survival function. In other words, the inherent
scatter in the correlations, not the measurement error, is the dominant source
of uncertainty in the predicted $f_{\rm esc}$, and the survival function
bounds serve as a reasonable estimate of the uncertainty in the $f_{\rm esc}$
predictions.

####  2.2.2 Goodness-of-Fit Metrics

As in Paper I, we evaluate the Cox models’ performance using several
complementary metrics. The concordance index, $C$, is particularly useful, as
it is appropriate for censored data. The concordance index assesses whether
the model $f_{\rm esc}$ predictions correctly sort the dataset in the order of
its observed $f_{\rm esc}$. To evaluate this sorting, the concordance index
calculation compares each possible pair of data points. Concordant pairs are
those where the galaxy with higher observed $f_{\rm esc}$ also has higher
predicted $f_{\rm esc}$, discordant pairs are the opposite, and tied pairs
have identical predicted $f_{\rm esc}$. Some pairs with upper limits in
$f_{\rm esc}$ lead to ambiguous rankings and do not appear in the concordance
index calculation. With all pairs evaluated, the concordance index is
calculated as

$C=\frac{n_{c}+0.5n_{t}}{n_{c}+n_{d}+n_{t}},$ (6)

where $n_{c}$, $n_{d}$, and $n_{t}$ are the number of concordant, discordant,
and tied pairs, respectively. A value of $C=1.0$ indicates a perfect rank
ordering, 0.5 is perfectly random ordering, and 0 is perfect disagreement. The
major advantage of $C$ is that it includes both LyC detections and non-
detections. However, it only assesses the relative order of the predicted
$f_{\rm esc}$ values, not their accuracy. To measure the latter, we turn to
alternative quantities. One such metric is the $R^{2}$ statistic

$R^{2}=1-\frac{\sum_{i}(y_{i}-f_{i})^{2}}{\sum_{i}(y_{i}-{\bar{y}})^{2}},$ (7)

where $y_{i}$ are the observed values of $\log_{\rm 10}$($f_{\rm esc}$),
${\bar{y}}$ is their mean value, and $f_{i}$ are the predicted $f_{\rm esc}$
values from the Cox model. As in Paper I and Maji et al. (2022), we also
calculate a variant of $R^{2}$, the adjusted $R^{2}$, which accounts for the
number of free parameters $p$ in the model and number of data points $n$:

$R^{2}_{\rm adj}=1-(1-R^{2})\frac{n-1}{n-p-1}.$ (8)

When increasing the number of variables, $R^{2}_{\rm adj}$ increases only if a
variable improves the predictions more than expected by chance. Finally, we
report the root-mean-square (RMS) dispersion

${\rm RMS}=\sqrt{\frac{\sum_{i}(y_{i}-f_{i})^{2}}{n}}.$ (9)

As explained in Paper I, we evaluate these three quantities ($R^{2}$,
$R^{2}_{\rm adj}$, and RMS) using $\log_{\rm 10}$($f_{\rm esc}$). The scatter
in the $\log_{\rm 10}$($f_{\rm esc}$) predictions is relatively consistent
across the full range of observed $f_{\rm esc}$, whereas the scatter in the
linear $f_{\rm esc}$ changes systematically across this range. We can only
calculate the $R^{2}$, $R^{2}_{\rm adj}$, and RMS metrics for the LzLCS+
galaxies with both detected LyC and non-zero predicted $f_{\rm esc}$, since
these metrics require observed values of $\log_{\rm 10}$($f_{\rm esc}$) and
finite predicted $\log_{\rm 10}$($f_{\rm esc}$). Hence, these metrics only
indicate the accuracy of the model predictions for LCEs. However, the most
successful models according to the $R^{2}$, $R^{2}_{\rm adj}$, and/or RMS
metrics also tend to have high $C$ and vice versa (see Paper I).

####  2.2.3 Choice of Input Variables

In generating Cox models from the LzLCS+, we consider a limited set of input
variables. Although the dependent variable in the Cox model can contain
measurement limits, the independent variables cannot. We therefore only use
independent variables that have measurements for nearly all the LzLCS+
galaxies, which requires us to exclude fainter emission line measurements
(e.g., [O i] $\lambda$6300 and [S ii] $\lambda\lambda$6716,6731) or
measurements not widely available for the full sample (e.g., Ly$\alpha$
velocity peak separation). In this paper, all 88 galaxies in the LzLCS+ have
measurements for our chosen input variables with one exception. The non-leaker
J1046+5827 does not have a reported $\Sigma_{\rm SFR}$ or $r_{\rm 50,NUV}$,
and we do not include it in deriving Cox models that use these variables. The
Cox model may also fail to converge if we include highly collinear variables,
which trace nearly identical properties. Given this limitation, we choose to
use only one measure of UV dust attenuation (E(B-V)UV or $\beta_{\rm 1550}$)
per model, and we do not include both $\Sigma_{\rm SFR}$ and $r_{\rm 50,NUV}$
in a single model.

We explored a variety of variable combinations in Paper I, and our best model
attained $R^{2}=0.69$, $R^{2}_{\rm adj}=0.60$, RMS $=0.31$ dex, and $C=0.91$.
We provide the parameters for this model in the Appendix. However, this model
included the average EW of Lyman-series absorption lines (EW(H i,abs)) and
$f_{{\rm esc,Ly}\alpha}$, both of which will be affected by the IGM at $z>6$.
Models in Paper I without UV absorption lines or Ly$\alpha$ showed higher
scatter but an overall ability to identify LCEs, with $R^{2}=0.29-0.40$,
$R^{2}_{\rm adj}=0.14-0.35$, RMS $=0.44-0.47$ dex, and $C=0.83$ (Paper I; see
Appendix for these model parameters). In this paper, we are specifically
concerned with variables that are easily observable at $z>6$ or measured in
large samples at $z\sim 3$. Hence, we likewise omit UV absorption line
measurements, and we avoid Ly$\alpha$ measurements for most of the models
tailored to $z>6$ galaxies.

##  3 Testing the Cox Models at $z\sim 3$

We apply the Cox models developed on the LzLCS+ sample to galaxies at $z\sim
3$ and $z\gtrsim 6$. We use the $z\sim 3$ galaxies to test whether the models
based on the low-redshift LzLCS+ galaxies can successfully predict $f_{\rm
esc}$ for galaxies at high redshift. We then apply the Cox models to $z\gtrsim
6$ samples to predict $f_{\rm esc}$ for galaxies in the epoch of reionization.

###  3.1 High-redshift Datasets and Models

First, we compile samples of $z\sim 3$ galaxies with reported global absolute
LyC $f_{\rm esc}$ and at least three of the input variables from Paper I,
which include stellar mass $M_{*}$, $M_{\rm 1500}$, nebular EWs, metallicity,
optical nebular line ratios, $\Sigma_{\rm SFR}$, $r_{\rm 50,NUV}$, E(B-V)UV,
$\beta_{\rm 1550}$, and Ly$\alpha$ measurements. The $z\sim 3$ LyC
measurements include individual detections, reported LyC upper limits, and
stacked samples that average over variations in IGM attenuation (Hainline et
al., 2009; Vasei et al., 2016; James et al., 2018; de Barros et al., 2016;
Vanzella et al., 2016; Steidel et al., 2018; Pahl et al., 2021; Bassett et
al., 2019; Fletcher et al., 2019; Nakajima et al., 2020; Bian & Fan, 2020;
Marques-Chaves et al., 2021, 2022; Liu et al., 2023; Kerutt et al., 2023). We
exclude one AGN from the Fletcher et al. (2019) sample and one shock-dominated
galaxy from the Bassett et al. (2019) sample. Some well-known $z\sim 3$ LCEs
do not appear in our sample, because they lack published global absolute
$f_{\rm esc}$ measurements (e.g., Ion3 and the Sunburst Arc). We also note
that some of the high-redshift samples fall outside of the parameter space
probed by the LzLCS+, such that our model predictions will extrapolate for
these galaxies. The $z\sim 3$ data has limitations as well; the compiled
$z\sim 3$ samples differ in the methods used for $f_{\rm esc}$ and input
variable measurements, and individual $z\sim 3$ $f_{\rm esc}$ measurements
have high uncertainty due to unknown variations in IGM attenuation. We discuss
the limitations of the $z\sim 3$ comparison further in the following sections.

In Paper I, we introduced a JWST Cox Model, which could be applied to an ideal
$z>6$ sample, with eight relevant variables from the LzLCS+ included. We also
found that a “limited JWST Model”, fit using only the three top-ranked
variables ($\beta_{\rm 1550}$, $\log_{\rm 10}$(O32), and $\log_{\rm
10}$($\Sigma_{\rm SFR}$)) performed equally well for predicting $f_{\rm esc}$
in the LzLCS+ sample. Unfortunately, many of the required variables for both
the full and limited JWST models have not been measured for large samples of
$z\sim 3$ LCEs or $z>6$ galaxies. Consequently, we run new models limited to
the variables available for the $z\sim 3$ and for $z\gtrsim 6$ samples. In
Table 1 and Table 2, we list these models, the variables they include, and the
samples to which they apply.

For the $z\gtrsim 6$ samples, we prioritize models that can apply to large
samples (e.g., Endsley et al., 2021, 2023; Morishita et al., 2023), models
that can apply to faint galaxies (Atek et al., 2023), and models that have at
least two of the most statistically significant variables in Paper I (a
measure of dust attenuation plus a measurement of ionization or morphology).
The “TopThree” model includes only the three top-ranked variables from Paper
I: dust, O32, and $\Sigma_{\rm SFR}$. In the TopThree model, we use E(B-V)UV
rather than $\beta_{\rm 1550}$, since it enables us to compare with a larger
$z\sim 3$ sample. As noted in Chisholm et al. (2022), the $\beta_{\rm 1550}$
and E(B-V)UV measurements for the LzLCS+ track each other almost perfectly and
provide equivalent information. Most models in Table 1 use E(B-V)UV; models
with “$\beta$” in their name use $\beta_{\rm 1550}$ instead, and one model,
LAE-O32-nodust, has no dust attenuation measurement. All models except the
TopThree model include $M_{\rm 1500}$, and most models include $M_{*}$ as
well. The only models without $M_{*}$ are the LAE, LAE-O32-nodust, and
ELG-O32-$\beta$ models.

We provide the best-fit coefficients and cumulative baseline hazards for each
model in the Appendix, which can be used to derive the median predicted
$f_{\rm esc}$ for a given galaxy and the uncertainty from the 16-84th
percentiles of the $f_{\rm esc}$ probability distribution. As with the LzLCS+
galaxies, we use an MC method to sample the variable uncertainties for the
$z\sim 3$ and $z\gtrsim 6$ observations and re-generate the predicted $f_{\rm
esc}$. We again find that the uncertainty estimated from the Cox model
survival function dominates over the uncertainty from sampling the input
variables in nearly all cases.

Table 1: Cox Models for High-Redshift Predictions Model | Variables
---|---
| Dust | Ly$\alpha$ | Nebular | Luminosity | Morphology
TopThree | E(B-V)UV | — |  $\log_{\rm 10}$(O32) | — |  $\log_{\rm 10}$($\Sigma_{\rm SFR}$)
LAE | E(B-V)UV | EW(Ly$\alpha$) | — | $M_{\rm 1500}$ | —
LAE-O32 | E(B-V)UV | EW(Ly$\alpha$) |  $\log_{\rm 10}$(O32) |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) | —
LAE-O32-nodust | — | EW(Ly$\alpha$) |  $\log_{\rm 10}$(O32) | $M_{\rm 1500}$ | —
ELG-EW | E(B-V)UV | — |  $\log_{\rm 10}$(EW([O III]+H$\beta$)) |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) | —
ELG-O32 | E(B-V)UV | — |  $\log_{\rm 10}$(O32) |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) | —
ELG-O32-$\beta$ | $\beta_{\rm 1550}$ | — |  $\log_{\rm 10}$(O32) | $M_{\rm 1500}$ | —
ELG-O32-$\beta$-Ly$\alpha$ | $\beta_{\rm 1550}$ | $f_{{\rm esc,Ly}\alpha}$ |  $\log_{\rm 10}$(O32) |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) | —
R50-$\beta$ | $\beta_{\rm 1550}$ | — | — |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) |  $\log_{\rm 10}$($r_{\rm 50,NUV}$))
$\beta$-Metals | $\beta_{\rm 1550}$ | — | 12+$\log_{\rm 10}$(O/H) |  $M_{\rm 1500}$, $\log_{\rm 10}$($M_{*}$) | —

=2.5in

Table 2: Cox Model High-Redshift Samples

Model |  $z\sim 3$ Samples |  $z\sim 6$ Samples
---|---|---
TopThree | 1 Detection: | Ion2aaFor this variable, a positive EW denotes net absorption. |
| 1 Non-Detection: | Cosmic Horseshoebbfootnotemark:  |
| 1 with Input Limits: | J1316+2614ccfootnotemark:  |
LAE | 28 Detections+13 Stacks: | Ion2aaFor this variable, a positive EW denotes net absorption., Bassett et al. (2019), Fletcher et al. (2019)ddfootnotemark: , Liu et al. (2023)eefootnotemark:  | Saxena et al. (2023b)
|  |  Kerutt et al. (2023), J0121+0025fffootnotemark: , J1316+2614ccfootnotemark: , Steidel et al. (2018)ggfootnotemark:  |
| 6 Non-Detections+5 Stacks: |  Bassett et al. (2019), Steidel et al. (2018)ggfootnotemark: , Bian & Fan (2020) |
| 3 with Input Limits: |  Fletcher et al. (2019)ddfootnotemark:  |
LAE-O32 | 4 Detections | Ion2aaFor this variable, a positive EW denotes net absorption., Bassett et al. (2019), Fletcher et al. (2019)ddfootnotemark:  | Saxena et al. (2023b)
| 5 Non-Detections: | Bassett et al. (2019) |
| 5 with Input Limits: |  Bassett et al. (2019), Fletcher et al. (2019)ddfootnotemark: , J1316+2614ccfootnotemark:  |
LAE-O32-nodust | 5 Detections+1 Stack: | Ion2aaFor this variable, a positive EW denotes net absorption., Bassett et al. (2019), Nakajima et al. (2020), J1316+2614ccfootnotemark:  | Saxena et al. (2023b)
| 14 Non-Detections+2 Stacks: |  Bassett et al. (2019), Nakajima et al. (2020) |
| 19 with Input Limits: |  Bassett et al. (2019), Nakajima et al. (2020) |
ELG-EW | 4 Detections: | Ion2aaFor this variable, a positive EW denotes net absorption., Fletcher et al. (2019)ddfootnotemark:  |  Endsley et al. (2021, 2023),
| 4 with Input Limits: |  Fletcher et al. (2019)ddfootnotemark: , J1316+2614ccfootnotemark:  |  Bouwens et al. (2023), Tang et al. (2023)hhfootnotemark: ,
|  |  |  Fujimoto et al. (2023), Saxena et al. (2023b)
ELG-O32 | 4 Detections: | Ion2aaFor this variable, a positive EW denotes net absorption., Bassett et al. (2019), Fletcher et al. (2019)ddfootnotemark:  |  Williams et al. (2023), Tang et al. (2023)hhfootnotemark: ,
| 5 Non-Detections: | Bassett et al. (2019) |  Fujimoto et al. (2023), Saxena et al. (2023b)
| 5 with Input Limits: |  Bassett et al. (2019), Fletcher et al. (2019)ddfootnotemark: , J1316+2614ccfootnotemark:  |
ELG-O32-$\beta$ | 2 Detections: | Ion2aaFor this variable, a positive EW denotes net absorption., J1316+2614ccfootnotemark:  |  Williams et al. (2023), Schaerer et al. (2022)iifootnotemark: ,
|  |  |  Mascia et al. (2023a), Saxena et al. (2023b, a)
ELG-O32-$\beta$-Ly$\alpha$ |  |  | Saxena et al. (2023b)
R50-$\beta$ | 1 Detection: | Ion2aaFor this variable, a positive EW denotes net absorption. |  Mascia et al. (2023a), Morishita et al. (2023)jjfootnotemark:
$\beta$-Metals | 1 Detection: | Ion2aaFor this variable, a positive EW denotes net absorption. |  Schaerer et al. (2022)iifootnotemark: ,
| 1 with Input Limits: | J1316+2614ccfootnotemark:  |  Williams et al. (2023), Atek et al. (2023)kkfootnotemark:

Note. — (a) From de Barros et al. (2016); Vanzella et al. (2016, 2020). (b)
From Hainline et al. (2009); Vasei et al. (2016); James et al. (2018). We
adopt the average dust-corrected O32 ratio of the two similar-flux regions in
Hainline et al. (2009). The E(B-V)UV and $\Sigma_{\rm SFR}$ values are flux-
weighted averages of the regions in James et al. (2018). (c) From Marques-
Chaves et al. (2022).(d) From Fletcher et al. (2019); Nakajima et al.
(2020).(e) We adopt the $f_{\rm esc}$ determined from the SED-fitting method
in Liu et al. (2023).(f) From Marques-Chaves et al. (2021).(g) From Steidel et
al. (2018); Pahl et al. (2021).(h) We adopt $M_{*}$ calculated using the
bursty non-parametric star formation history, but include the alternative
values in the uncertainties. We use the direct method metallicities for the
two galaxies with secure [O iii] $\lambda$4363 detections and the modeled
metallicities for the others.(i) We correct $M_{*}$ for magnification, using
the values in Schaerer et al. (2022). For the source without a direct method
metallicity, we adopt the strong-line method 12+$\log_{\rm 10}$(O/H).(j) We
exclude sources flagged as having significant residuals in the morphological
fits.(k) Metallicities calculated using strong-line methods.

Table 3: Metrics for High-Redshift Cox Models

LzLCS+ Sample
---
Model | $N_{\rm Gal}$ | $N_{\rm Detect}$ | $R^{2}$ | $R^{2}_{\rm adj}$ | RMS | $C$
TopThree | 87 | 45 | 0.38 | 0.34 | 0.44 | 0.82
LAE | 88 | 46 | 0.21 | 0.15 | 0.50 | 0.81
LAE-O32 | 88 | 44 | 0.37 | 0.29 | 0.45 | 0.82
LAE-O32-nodust | 88 | 43 | 0.02 | -0.06 | 0.59 | 0.77
ELG-EW | 88 | 43 | 0.14 | 0.05 | 0.53 | 0.79
ELG-O32 | 88 | 42 | 0.42 | 0.36 | 0.44 | 0.79
ELG-O32-$\beta$ | 88 | 42 | 0.30 | 0.25 | 0.48 | 0.79
ELG-O32-$\beta$-Ly$\alpha$ | 88 | 43 | 0.40 | 0.32 | 0.45 | 0.81
R50-$\beta$ | 87 | 45 | 0.36 | 0.30 | 0.44 | 0.84
$\beta$-Metals | 88 | 44 | 0.11 | 0.02 | 0.55 | 0.77
$z\sim 3$ Sample
Model | $N_{\rm Gal}$ | $N_{\rm Detect}$ | $R^{2}$ | $R^{2}_{\rm adj}$ | RMS | $C$
TopThree | 2 | 1 | — | — | 0.15 | 1.00
LAE | 52 | 35 | -4.18 | — | 0.92 | 0.49
LAE-O32 | 9 | 4 | -0.35 | — | 0.32 | 0.71
LAE-O32-nodust | 22 | 6 | -0.12 | — | 0.31 | 0.80
ELG-EW | 4 | 4 | -4.81 | — | 0.65 | 0.50
ELG-O32 | 9 | 4 | 0.09 | — | 0.27 | 0.71
ELG-O32-$\beta$ | 2 | 2 | -0.91 | — | 0.04 | 1.00
ELG-O32-$\beta$-Ly$\alpha$ | — | — | — | — | — |
R50-$\beta$ | 1 | 1 | — | — | 0.09 | —
$\beta$-Metals | 1 | 1 | — | — | 0.43 | —
Combined Sample
Model | $N_{\rm Gal}$ | $N_{\rm Detect}$ | $R^{2}$ | $R^{2}_{\rm adj}$ | RMS | $C$
TopThree | 89 | 46 | 0.43 | 0.39 | 0.44 | 0.83
LAE | 140 | 81 | -0.47 | -0.53 | 0.71 | 0.70
LAE-O32 | 97 | 48 | 0.42 | 0.35 | 0.44 | 0.82
LAE-O32-nodust | 110 | 49 | 0.19 | 0.13 | 0.56 | 0.78
ELG-EW | 92 | 47 | 0.16 | 0.08 | 0.54 | 0.79
ELG-O32 | 97 | 46 | 0.48 | 0.42 | 0.43 | 0.80
ELG-O32-$\beta$ | 90 | 44 | 0.40 | 0.36 | 0.47 | 0.80
ELG-O32-$\beta$-Ly$\alpha$ | — | — | — | — | — | —
R50-$\beta$ | 88 | 46 | 0.41 | 0.35 | 0.43 | 0.84
$\beta$-Metals | 89 | 45 | 0.17 | 0.08 | 0.55 | 0.77

Note. — $N_{\rm Gal}$ is the number of galaxies assessed in the fit and does
not include any galaxies with limits for input variables. $N_{\rm Detect}$ is
the number of galaxies with LyC detections, finite $f_{\rm esc}$ predictions,
and no upper or lower limits for input variables; the $R^{2}$ and RMS metrics
use only these galaxies. The $R^{2}$ statistic measures how well the
predictive model explains the observed variance in the $\log_{\rm 10}$($f_{\rm
esc}$) data, with higher $R^{2}$ values corresponding to more accurate models.
The $R^{2}_{\rm adj}$ metric accounts for the number of parameters used in the
model and increases only if a variable improves the fit more than expected by
chance. RMS is the root-mean-square dispersion of the predicted vs. observed
$\log_{\rm 10}$($f_{\rm esc}$) for the LyC detections. Higher $C$ values
indicate that the model more accurately sorts the observations in the correct
order of increasing $f_{\rm esc}$. $C$ includes both detections and galaxies
with $f_{\rm esc}$ upper limits. See §2.2.2 for a full description of these
statistics.

We first test the Cox model predictions on $z\sim 3$ samples to see whether
the high-redshift LCEs behave similarly to the LzLCS+ sample. We derive Cox
models for the variable sets in Table 1 using only the LzLCS+ data; the high-
redshift samples are not included in the fitting process. We then use the
models to predict $f_{\rm esc}$ for each set of $z\sim 3$ measurements, and we
calculate goodness-of-fit metrics for the LzLCS+ sample alone, the high-
redshift sample alone, and for the combined sample of low- and high-redshift
galaxies. We list these goodness-of-fit metrics in Table 3.

###  3.2 Model Performance

Figures 1 and 2 show the predicted vs. observed $f_{\rm esc}$ for the LzLCS+
and $z\sim 3$ samples. As seen from the plots and the goodness-of-fit
statistics in Table 3, these models do not perform as well as the fiducial
model from Paper I, but several (TopThree, LAE-O32, ELG-O32, ELG-O32-$\beta$,
ELG-O32-$\beta$-Ly$\alpha$, R50-$\beta$) are comparable to or better than the
Paper I JWST model (see Tables 11-13 for the list of input variables used in
the Paper I models). The metrics for the LzLCS+ sample are $R^{2}=0.02-0.42$
(vs. 0.59 for the fiducial model and 0.29 for JWST), $R^{2}_{\rm adj}=-0.06$
to 0.36 (vs. 0.48 for the fiducial model and 0.14 for JWST), RMS$=0.44-0.59$
(vs. 0.37 for the fiducial model and 0.47 for JWST), and $C=0.77$ to 0.84 (vs.
0.88 for the fiducial model and 0.83 for JWST).

Figure 1: The $f_{\rm esc}$ predictions from the TopThree (a), LAE (b),
LAE-O32 (c), and LAE-O32-nodust (d) models. See Table 1 for model
descriptions. Red circles represent LzLCS+ LyC detections, and blue triangles
represent upper limits. We plot $z\sim 3$ galaxies in green, with stars
denoting LyC detections and triangles denoting upper limits. We identify the
two strongest high-redshift LCEs, Ion2 and J1316+2614, by a teal and a green
diamond, respectively. Light green symbols indicate high-redshift galaxies
that have a limit for one or more input variables. We draw circles around data
points representing high-redshift galaxy stacks, which are less subject to
uncertainty in IGM attenuation. The error bar in the upper left corner
indicates the median size of the uncertainties in the observed and predicted
$f_{\rm esc}$ for the combined sample of low- and high-redshift galaxies. The
dashed line shows a one-to-one correspondence. Several of the Cox models
predict $f_{\rm esc}$ for both low- and high-redshift galaxies with comparable
accuracy.  Figure 2: The $f_{\rm esc}$ predictions from the ELG-EW (a),
ELG-O32 (b), ELG-O32-$\beta$ (c), ELG-O32-$\beta$-Ly$\alpha$ (d), R50-$\beta$
(e), and $\beta$-Metals (f) models. See Table 1 for model descriptions. Ion2
is shown by an open symbol in panel d, because its prediction is based on a
limit in $f_{{\rm esc,Ly}\alpha}$. Other symbols are the same as in Figure 1.

The models developed on the LzLCS+ sample generally reproduce the observed
$f_{\rm esc}$ values of the $z\sim 3$ LCEs. In fact, the high-redshift samples
often have a lower RMS scatter than the predictions for the LzLCS+, and in
several models, the $R^{2}$ and $C$ values for the combined high- and low-
redshift sample are comparable to or higher than the $R^{2}$ and $C$ values
for the LzLCS+ sample alone. Notably, predictions for the two strongest high-
redshift LCEs, Ion2 and J1316+2614, match the observations nearly perfectly in
most cases (Figures 1 and 2), demonstrating the success of the model in
identifying extreme LCEs.

The goodness-of-fit metrics give clues as to which variables are most
important in predicting $f_{\rm esc}$. As seen in Table 3, the models with the
highest $R^{2}$ and lowest RMS for the LzLCS+ sample are the TopThree,
LAE-O32, ELG-O32 and variants, and the R50-$\beta$ models. These same models
are the only ones that include both a measure of UV dust attenuation and
either O32 or radius as variables. Turning to the $C$ metric, which includes
non-detections, the best-performing models for the LzLCS+ sample are the
R50-$\beta$, TopThree, LAE-O32, ELG-O32-$\beta$-Ly$\alpha$, and LAE models. To
properly estimate $f_{\rm esc}$ for weak and non-LCEs, UV dust attenuation
again appears important, but each of the best-performing models also includes
either morphological information or Ly$\alpha$ measurements. These additional
variables may help distinguish non-LCEs from LCEs.

By all metrics, the LAE-O32-nodust and $\beta$-Metals models perform the worst
for the LzLCS+ galaxies. The LAE-O32-nodust model is the only one that lacks
E(B-V)UV or $\beta_{\rm 1550}$, which highlights the crucial role of dust
extinction in LyC escape. Conversely, the $\beta$-Metals model relies almost
exclusively on dust attenuation to infer $f_{\rm esc}$ and lacks information
on O32, Ly$\alpha$, or morphology, which may better constrain the LyC
absorption due to H i.

####  3.2.1 Model Performance for $z\sim 3$ Samples

Comparing the models’ performance for the high-redshift and combined low- and
high-redshift samples is difficult, because each high-redshift model applies
to a different set of galaxies. For example, the ELG-O32-$\beta$ model has the
lowest RMS for the high-redshift galaxies, but this RMS is based on a single
galaxy: Ion2. Ion2 is the only galaxy that is included in most of the models
(all but ELG-O32-$\beta$-Ly$\alpha$, where it has an input limit). For Ion2,
we adopt an observed absolute $f_{\rm esc}$=0.75, in the middle of the 0.5-1
bounds reported by Vanzella et al. (2016). Using this $f_{\rm esc}$, Ion2’s
LyC escape is predicted most accurately by the ELG-O32-$\beta$ (RMS=0.02 dex)
and ELG-O32 (RMS=0.04 dex) models. For Ion2, dust attenuation, luminosity, and
ionization appear to be key factors in predicting its $f_{\rm esc}$, although
most other models are still consistent with its reported range of $f_{\rm
esc}$ (RMS=0.09-0.15 dex for the next six models). The $\beta$-Metals model,
which lacks information about O32 and morphology, is the only inaccurate model
(RMS=0.43 dex). More surprisingly, however, the TopThree model gives the
lowest predicted $f_{\rm esc}$=0.53 (RMS=0.15 dex), which highlights the role
of UV magnitude in the $f_{\rm esc}$ predictions, a point we discuss further
in §4.

Nine galaxies, including four LCEs, appear in each of the LAE, LAE-O32,
LAE-O32-nodust, and ELG-O32 model test samples. Although the LAE-O32-nodust
model performed worst for the LzLCS+, for this common subsample of nine high-
redshift galaxies, the LAE-O32-nodust model has the lowest RMS (0.09 dex) and
highest $C$ (0.76), compared to RMS values of 0.27-0.32 dex and $C$ values of
0.65-0.71 for these galaxies in the other three models. Given the generally
successful performance of all four of these models, models with O32 and/or
EW(Ly$\alpha$) appear able to predict $f_{\rm esc}$ for high-redshift galaxies
equally well. Hence, while models with E(B-V)UV and O32 or morphology work
best for the LzLCS+ sample, we cannot rule out the possibility that
EW(Ly$\alpha$) may also be important in predicting $f_{\rm esc}$ at $z\sim 3$.

Nevertheless, using EW(Ly$\alpha$) as the sole nebular measurement is not
sufficient to accurately predict $f_{\rm esc}$. The LAE model (Figure 1b),
which applies to the largest set of $z\sim 3$ LyC measurements, is one of the
least successful models at predicting $f_{\rm esc}$ in the high-redshift
galaxies. As seen in Table 3, this model has the largest RMS scatter for the
$z\sim 3$ sample (0.9 dex), lowest $C$ (0.49), and the second worst $R^{2}$
(-4.2). The other model with comparably poor $R^{2}$ and $C$ for $z\sim 3$
galaxies is the ELG-EW model, which is the only other model that exclusively
relies on an emission line EW for its input nebular information. Both models
also have some of the lowest $R^{2}$ values for the LzLCS+ sample, which
reflects the fact that these models substantially under-predict the $f_{\rm
esc}$ of the strongest LCEs in the LzLCS+. The LAE model shows this same
tendency towards under-prediction for the $z\sim 3$ LCEs as well; 13 $z\sim 3$
galaxies with observed $f_{\rm esc}$ $\geq 0.2$ are predicted to have $f_{\rm
esc}$ of only $\leq 0.03$ by the model.

The poor performance of the LAE model may result from limitations in both the
model itself and in the $z\sim 3$ data. The fact that both the LAE and ELG-EW
models fare poorly for the LzLCS+ and for the high-redshift galaxies,
especially for strong LCEs, suggests that nebular EW is a flawed tracer of
$f_{\rm esc}$. Indeed, [O iii] EW and Ly$\alpha$ EW should be low both for the
weakest and for the strongest LCEs. In non-leakers, a high optical depth and
corresponding lack of Ly$\alpha$ escape should lead to a low Ly$\alpha$ EW.
Low [O iii] and Ly$\alpha$ EWs could also indicate a weak current burst of
star formation, without significant feedback or LyC production, which might
result in undetected LyC. At lower optical depths and moderate $f_{\rm esc}$,
Ly$\alpha$ EW may increase due to enhanced escape, and [O iii] EWs may
likewise increase if these lower optical depths are preferentially associated
with stronger starbursts. However, the [O iii] and Ly$\alpha$ EWs will
decrease again in the strongest LCEs, as a general lack of nebular gas results
in the limited production of nebular emission lines (e.g., Zackrisson et al.,
2013; Nakajima & Ouchi, 2014). Thus, EW may be an ambiguous indicator of LyC
escape.

The $z\sim 3$ data also have limitations. The strong LCEs that are most
severely under-predicted by the LAE model all come from the Liu et al. (2023)
and Kerutt et al. (2023) samples of individual LAEs with LyC detections. Both
papers note that their LyC detections represent the extreme of the population
and may not be representative of the average galaxy population with these
parameters, whose $f_{\rm esc}$ may be much lower. In Paper I, we found a
similar result for the strongest LCEs in the LzLCS+, whose only distinguishing
feature compared to more moderate LCEs is their low line-of-sight H i content,
suggestive of a favorable orientation. Another limitation of the high-redshift
data is the uncertainty in the IGM attenuation. At $z\sim 3$, this attenuation
is significant and varies along different lines of sight (e.g., Rudie et al.,
2013; Inoue et al., 2014; Vanzella et al., 2016; Steidel et al., 2018). This
unknown sightline dependence leads to additional uncertainty in the reported
$f_{\rm esc}$ for individual galaxies and could cause some genuine LCEs to
appear as non-detections or to have overestimated $f_{\rm esc}$.

The stacked samples (Steidel et al., 2018; Nakajima et al., 2020; Bian & Fan,
2020) shown circled in Figure 1, average over these variations and are less
subject to this uncertainty. In the LAE model in Figure 1b, we see that the
stacked samples do indeed show less scatter than the individual high-redshift
detections. However, the model does not perfectly predict $f_{\rm esc}$ for
the stacked samples, and their scatter is still comparable to the scatter in
the individual LzLCS+ galaxies. In addition to averaging over IGM attenuation
variations, the stacks also average over any galaxy-to-galaxy variation in
physical properties. Hence, the input variables may not represent the true set
of properties of an individual system, leading to some uncertainty in the
predicted $f_{\rm esc}$.

Unfortunately, the small sample sizes in most models and the lack of a common
high-redshift sample across the models make it difficult to discern which
model parameters are most important for predicting $f_{\rm esc}$ at high
redshift. Our results suggest that if high-redshift galaxies behave like their
lower redshift counterparts in the LzLCS+ sample, E(B-V)UV or $\beta_{\rm
1550}$, O32, and $\Sigma_{\rm SFR}$ or half-light radius are essential
variables to include. For the existing $z\sim 3$ samples, models with at least
two of these parameters do successfully predict the $f_{\rm esc}$ of the
available high-redshift galaxies to within $\sim$0.3 dex.

This performance also shows that the LzLCS+ galaxies may indeed be reasonable
analogs for high-redshift galaxies and that the same observable and physical
properties may correlate with LyC escape at both low and high redshift (e.g.,
Saldana-Lopez et al., 2023; Schaerer et al., 2022; Mascia et al., 2023a). This
agreement is not trivial, as in principle, the $z\sim 0.3$ and $z\sim 3$
galaxies could have different star formation histories, dust properties,
morphologies, or other properties, any of which could affect $f_{\rm esc}$.
The LzLCS+ parameter space does cover the properties of the $z\sim 3$ samples,
with some exceptions. Seventeen $z\sim 3$ targets (8 LCEs and 9 non-LCEs) are
brighter than the LzLCS+ UV magnitude range by 0.1-3.18 magnitudes, including
the strong LCEs Ion2 and J1316+2614. An additional two LCEs in the LAE model
are fainter than the LzLCS+ by 0.5-0.7 magnitudes. The predictions for most
models, all those with $M_{\rm 1500}$ as a variable, therefore extrapolate to
an unobserved part of parameter space, yet still perform well.

The data in Figures 1 and 2 and metrics in Table 3 include $z\sim 3$ galaxies
that have measured values for all of the input independent variables in the
models222The reported E(B-V) for Ion2 is a stringent upper limit of
E(B-V)$<0.04$. For consistency with Ion2’s other reported measurements, we
follow Vanzella et al. (2016) in adopting E(B-V)=0.0, and we choose to include
Ion2 in our goodness-of-fit calculations.. However, additional galaxies in
these samples have upper or lower limits for some independent variables, which
can still provide potentially useful constraints on $f_{\rm esc}$. For these
galaxies, we adopt the limit as the value of the independent variable and
predict their $f_{\rm esc}$ using the LzLCS+ models. We plot these approximate
estimates in Figures 1 and 2 as light green symbols and an open symbol for
Ion2, but we do not include these galaxies in the goodness-of-fit metrics.

The galaxies with limits have one or more of the following: upper limits in
$M_{*}$, lower limits in $\Sigma_{\rm SFR}$, lower limits in EW([O
iii]+H$\beta$), lower limits in O32, lower limits in EW(Ly$\alpha$), lower
limits in $M_{\rm 1500}$, and lower limits in $f_{{\rm esc,Ly}\alpha}$. One
galaxy has an upper limit in EW([O iii]+H$\beta$) rather than a lower limit.
Given the coefficients for these variables in the models, the mass and
magnitude limits would cause the model to overestimate $f_{\rm esc}$, while
the lower limits in the nebular emission lines should generally cause the
model to underestimate $f_{\rm esc}$. For example, a galaxy with a lower limit
of EW(Ly$\alpha$) $>100$Å is an even stronger Ly$\alpha$ emitter than we
assume, and its predicted $f_{\rm esc}$ will be an underestimate.

Most of the $z\sim 3$ LCEs with input limits in the LAE-O32 and LAE-O32-nodust
models only have lower limits in O32, such that their predicted $f_{\rm esc}$
should be an underestimate. For galaxies with lower limits in both a nebular
property and $M_{\rm 1500}$, we cannot easily interpret the predicted $f_{\rm
esc}$ as a lower or upper limit. However, the best-fit coefficients in the
Appendix show that the models have a steeper dependence on the nebular lines
than they do on $M_{\rm 1500}$. Figures 1 and 2 show that for most of the
galaxies with input limits, the measured limits still constrain $f_{\rm esc}$
as accurately as for the rest of the sample. However, a few galaxies (90675,
101846, 105937 from Fletcher et al. 2019) have strongly under-predicted
$f_{\rm esc}$, which suggests that their nebular lines may be much stronger
than the reported limit.

Finally, we note that our model predictions do not account for systematic
uncertainties, including differences in methodology. We adopt the published
values for all data. However, each paper makes different assumptions, which
could affect the tabulated values. For instance, papers differ in their
adopted IGM transmission models, the models used in SED fits, and adopted dust
attenuation laws. These assumptions could lead to systematic differences in
properties like $f_{\rm esc}$, $M_{*}$, and E(B-V)UV among the different
samples. In addition, LyC measurements may be photometric or spectroscopic and
may cover different wavelength ranges. Despite these systematics, the models
work well in predicting $f_{\rm esc}$. Adopting consistent methodologies could
potentially result in better predictions; however, as shown by the scatter in
the LzLCS+ sample, which has a consistent methodology, the inherent
uncertainty in the model itself also limits the possible accuracy of
predictions.

####  3.2.2 Summary of $z\sim 3$ Results

In conclusion, Cox models derived using the $z\sim 0.3$ LzLCS+ sample can
successfully predict $f_{\rm esc}$ in $z\sim 3$ LCEs. This agreement suggests
that LCEs may have similar physical properties at both low and high redshift
(e.g. Saxena et al., 2023b). We find that the most accurate models include
E(B-V)UV and O32 or morphology measurements as variables. However, larger
high-redshift samples with a full suite of measurements are required to test
this result. Future observations of $z\sim 3$ LCEs with JWST will further
clarify whether the relationship between $f_{\rm esc}$ and physical properties
evolves with redshift or remains constant. Preliminary JWST observations
suggest that $z\gtrsim 6$ galaxies may indeed share numerous physical
properties with low-redshift analog samples like the LzLCS (e.g., Schaerer et
al., 2022; Endsley et al., 2023; Mascia et al., 2023a; Lin et al., 2024). The
models developed on the LzLCS+ sample and tested on $z\sim 3$ LCEs may
therefore apply equally well at the reionization epoch.

##  4 $f_{\rm esc}$ Predictions at $z\gtrsim 6$

###  4.1 Model Predictions for $z\gtrsim 6$ Samples

Given the success of our models in predicting $f_{\rm esc}$ at $z\sim 3$, we
now use the Cox models to estimate $f_{\rm esc}$ for galaxies at $z\gtrsim 6$,
where LyC is not detectable because of the IGM opacity. We consider several
models: the “ELG” models, R50-$\beta$, and $\beta$-Metals. The $\beta$-Metals
model is the most limited, as it only includes measurements of luminosity,
mass, UV slope, and metallicity and has no information on nebular line
strength or morphology. Based on the LzLCS+ galaxies, these models have an RMS
scatter of 0.44-0.55 dex in $f_{\rm esc}$ (see Section 3 and Table 3), with
the ELG-O32 or R50-$\beta$ models giving the most accurate predictions,
depending on the metric considered.

We can apply these models to several samples of reionization-era galaxies,
which span $z\sim 6-14$: photometric samples from Endsley et al. (2021),
Endsley et al. (2023), Bouwens et al. (2023), and Morishita et al. (2023) and
spectroscopic samples from Williams et al. (2023), Schaerer et al. (2022),
Tang et al. (2023), Fujimoto et al. (2023), Saxena et al. (2023b, a), Mascia
et al. (2023a), and Atek et al. (2023). The Endsley et al. (2021) sample is
based on Spitzer, HST, and ground-based observations; all other samples
incorporate JWST NIRCam photometry and several (Williams et al., 2023;
Schaerer et al., 2022; Tang et al., 2023; Fujimoto et al., 2023; Saxena et
al., 2023b, a; Mascia et al., 2023a; Atek et al., 2023) use JWST NIRSpec
spectroscopy as well. We select all galaxies at $z>5.9$ from these samples to
investigate $f_{\rm esc}$ in the epoch of reionization.

####  4.1.1 Photometric Samples

We list the predicted $f_{\rm esc}$ values and their associated uncertainties
for the $z\gtrsim 6$ galaxies in Tables 4-8. The ELG-EW and R50-$\beta$ models
(Tables 4 and 5) apply to the largest compilations of $z\gtrsim 6$ galaxies,
with 183 measurements corresponding to 180 unique galaxies for the ELG-EW
samples and 278 measurements at $z>5.9$ for the R50-$\beta$ samples, although
most of the galaxies for both models only have photometric redshifts. Most of
the galaxies in the ELG-EW samples are from Endsley et al. (2021, 2023), and
most of the R50-$\beta$ sample galaxies come from Morishita et al. (2023).

We show histograms of the ELG-EW and R50-$\beta$ model $f_{\rm esc}$
predictions in Figure 3. The median of the $f_{\rm esc}$ predictions for the
ELG-EW model samples is fairly low, $f_{\rm esc}$ $=0.047$, which is lower
than the average $f_{\rm esc}$ value of 0.1-0.2 required to reionize the
universe assuming a canonical ionizing photon production efficiency $\log_{\rm
10}$($\xi_{\rm ion}$)$\sim 25.3$(e.g., Finkelstein et al., 2015, 2019;
Robertson et al., 2015; Naidu et al., 2020). Of these $z\gtrsim 6$ galaxies,
27% have $f_{\rm esc}$ $\geq 0.1$ and only 6% have $f_{\rm esc}$ $\geq 0.2$.
However, as shown in Figure 2a, the ELG-EW model is one of the less accurate
models, with a greater tendency to underpredict the true $f_{\rm esc}$ of the
LzLCS+ sample. This underprediction implies that the ELG-EW variable set,
which includes EW([O iii]+H$\beta$) but not O32, does not distinguish strong
LCEs from weaker LCEs. The more accurate R50-$\beta$ model (see Table 3 and
Figure 2) implies a higher fraction of LCEs, with a median value of $f_{\rm
esc}$ $=0.14$, and 56% and 39% of the galaxies with $f_{\rm esc}$ $>0.1$ and
$f_{\rm esc}$ $>0.2$, respectively. Both models find a substantial fraction of
weak or non-leakers, and 25% of the galaxies in both models have $f_{\rm esc}$
$<0.03-0.04$. Taken at face value, our preliminary results would suggest that
$f_{\rm esc}$ values $>0.1$ are common but not ubiquitous among moderate to
bright galaxies ($M_{\rm 1500}<-18$) in the epoch of reionization (see also
Mascia et al. 2023a, b).

Table 4: Predictions for $z\gtrsim 6$ Galaxies from the ELG-EW Model

Source ID | $z_{\rm phot}$ | $z_{\rm spec}$ | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$ | Reference
---|---|---|---|---|---|---
COS-83688 | 6.70 | — | 0.190 | 0.056 | 0.602 | Endsley et al. (2021)
COS-87259 | 6.66 | — | 0.024 | 0.005 | 0.119 | Endsley et al. (2021)
COS-237729 | 6.87 | — | 0.119 | 0.033 | 0.488 | Endsley et al. (2021)
COS-312533 | 6.85 | — | 0.043 | 0.016 | 0.265 | Endsley et al. (2021)
COS-400019 | 6.88 | — | 0.109 | 0.031 | 0.463 | Endsley et al. (2021)

Note. — Predicted $f_{\rm esc}$ values from the ELG-EW model for the $z\gtrsim
6$ galaxies from Endsley et al. (2021, 2023), Bouwens et al. (2023), Tang et
al. (2023), Fujimoto et al. (2023), and Saxena et al. (2023b). $z_{\rm phot}$
and $z_{\rm spec}$ are the photometric and spectroscopic redshifts. $f_{\rm
esc,min}$ and $f_{\rm esc,max}$ represent the 15.9 and 84.1 percentiles of the
model $f_{\rm esc}$ predictions. For galaxies with upper limits on EW([O
iii]+H$\beta$), the $f_{\rm esc}$ predictions are also upper limits and are
marked accordingly. The Reference column lists the publication used for the
model input variables. The full, machine-readable version of this table is
available online. We show the first five rows as an example here.

Table 5: Predictions for $z\gtrsim 6$ Galaxies from the R50-$\beta$ Model

Source ID | $z_{\rm phot}$ | $z_{\rm spec}$ | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$ | Reference
---|---|---|---|---|---|---
J1235-15534 | 5.9 | — | 0.09 | 0.03 | 0.41 | Morishita et al. (2023)
JADESGDS-6734 | — | 5.92 | 0.11 | 0.03 | 0.44 | Morishita et al. (2023)
JADESGDS-11449 | — | 5.94 | 0.51 | 0.20 | 0.88 | Morishita et al. (2023)
JADESGDS-18169 | — | 5.94 | 1.00 | 0.61 | 1.00 | Morishita et al. (2023)
JADESGDS-33803 | — | 5.97 | 0.18 | 0.05 | 0.56 | Morishita et al. (2023)

Note. — Predicted $f_{\rm esc}$ values from the R50-$\beta$ model. $z_{\rm
phot}$ and $z_{\rm spec}$ are the photometric and spectroscopic redshifts.
$f_{\rm esc,min}$ and $f_{\rm esc,max}$ represent the 15.9 and 84.1
percentiles of the model $f_{\rm esc}$ predictions. The Reference column lists
the publication used for the model input variables. The full, machine-readable
version of this table is available online. We show the first five rows as an
example here.

We caution that some of the $z\gtrsim 6$ galaxies in the ELG-EW sample and
many in the R50-$\beta$ sample fall outside of the parameter space probed by
the LzLCS+, although a similar extrapolation did not seem to adversely affect
the $f_{\rm esc}$ predictions for the $z\sim 3$ galaxies (Section 3). Of the
183 measurements in the ELG-EW model sample, 14 galaxies have brighter UV
magnitudes than the LzLCS+ galaxies, one is fainter, and 17 have higher EW([O
iii]+H$\beta$). Likewise, the R50-$\beta$ sample contains 8 galaxies brighter
than and 30 galaxies fainter than the LzLCS+ sample. More concerningly, 114 of
the 278 galaxies in the R50-$\beta$ compilation are more compact than the
LzLCS+ galaxies, which could mean that their high inferred $f_{\rm esc}$
values result from an incorrect extrapolation into this compact regime.

####  4.1.2 Spectroscopic Samples

The ELG-O32 and ELG-O32-$\beta$ model predictions (Tables 6 and 7) apply to
fewer $z\gtrsim 6$ galaxies, only 17 and 27 galaxies, respectively, but all of
these galaxies have spectroscopic redshifts. Half the galaxies in the ELG-O32
model have predicted $f_{\rm esc}$ $<0.05$, although most of these are only
lower limits in $f_{\rm esc}$, and 25% of the ELG-O32-$\beta$ sample galaxies
have $f_{\rm esc}$ $\leq 0.01$. Like the R50-$\beta$ model predictions, low
$f_{\rm esc}$ is common, but the overall distribution also extends to very
high $f_{\rm esc}$. Both models include some galaxies with lower limits in
O32, whose predicted $f_{\rm esc}$ values therefore also correspond to lower
limits. Depending on how high the true $f_{\rm esc}$ for these galaxies are,
35-82% of the ELG-O32 model galaxies and 33-41% of the ELG-O32-$\beta$
galaxies have $f_{\rm esc}$ $\geq 0.2$.

Because these galaxies are all spectroscopically confirmed, these samples
could be biased toward galaxies with stronger emission lines and hence higher
$f_{\rm esc}$; the median EW([O iii]+H$\beta$) of the spectroscopic ELG-O32
sample is 1790 Å (Fujimoto et al., 2023; Tang et al., 2023; Saxena et al.,
2023b), compared with 690 Å for the photometric samples in the ELG-EW model,
for example (Endsley et al., 2021, 2023; Bouwens et al., 2023). In addition,
these models may not underpredict $f_{\rm esc}$ for strong LCEs to the same
extent as the ELG-EW model (see Figure 2 and Table 3). For the six galaxies
with predicted $f_{\rm esc}$ $>0.2$ in the ELG-O32 model, the ELG-EW
predictions are indeed lower by $\Delta$$f_{\rm esc}$=0.26 on average, which
suggests that model differences also account for some of the higher $f_{\rm
esc}$ values compared to the ELG-EW distribution.

Table 6: Predictions for $z\gtrsim 6$ Galaxies from the ELG-O32 Model

Source ID | Alternate Name | $z_{\rm spec}$ | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$ | Reference
---|---|---|---|---|---|---
CEERS-2 | CEERS1-3858 | 8.807 |  $>$0.178 |  $>$0.053 |  $>$0.567 | Fujimoto et al. (2023)
CEERS-3 | CEERS1-3908 | 8.005 |  $>$0.049 |  $>$0.023 |  $>$0.267 | Fujimoto et al. (2023)
CEERS-3 | CEERS1-3908 | 8.00 |  $>$0.006 |  $>$0 |  $>$0.032 | Tang et al. (2023)
CEERS-7 | CEERS1-6059 | 7.993 |  $>$0.119 |  $>$0.038 |  $>$0.453 | Fujimoto et al. (2023)
CEERS-20 | CEERS3-1748 | 7.769 |  $>$0 |  $>$0 |  $>$0 | Fujimoto et al. (2023)
CEERS-23 | CEERS6-7603 | 8.881 |  $>$0.008 |  $>$0 |  $>$0.041 | Fujimoto et al. (2023)
CEERS-23 | CEERS6-7603 | 8.881 |  $>$0.007 |  $>$0 |  $>$0.038 | Tang et al. (2023)
CEERS-24 | CEERS6-7641 | 8.998 |  $>$0.012 |  $>$0 |  $>$0.043 | Fujimoto et al. (2023)
CEERS-24 | CEERS6-7641 | 8.999 |  $>$0.013 |  $>$0 |  $>$0.049 | Tang et al. (2023)
CEERS-44 |  | 7.10 |  $>$0.402 |  $>$0.159 |  $>$0.731 | Tang et al. (2023)
CEERS-407 |  | 7.028 | 0.023 | 0.005 | 0.093 | Tang et al. (2023)
CEERS-498 |  | 7.18 |  $>$0.499 |  $>$0.201 |  $>$0.869 | Tang et al. (2023)
CEERS-499 |  | 7.168 |  $>$0 |  $>$0 |  $>$0 | Tang et al. (2023)
CEERS-698 |  | 7.47 | 0.722 | 0.498 | 1 | Tang et al. (2023)
CEERS-1019 |  | 8.678 | 1 | 0.644 | 1 | Tang et al. (2023)
CEERS-1025 |  | 8.715 | 0.301 | 0.110 | 0.622 | Tang et al. (2023)
CEERS-1027 |  | 7.819 |  $>$0.650 |  $>$0.453 | 1 | Tang et al. (2023)
CEERS-1038 |  | 7.194 |  $>$0.007 |  $>$0 |  $>$0.038 | Tang et al. (2023)
JADES-GS-z7-LA |  | 7.278 | 0 | 0 | 0.012 | Saxena et al. (2023b)
11027 |  | 9.51 | 0.021 | 0.005 | 0.071 | Williams et al. (2023)

Note. — Predicted $f_{\rm esc}$ values from the ELG-O32 model. $z_{\rm spec}$
is the spectroscopic redshift. $f_{\rm esc,min}$ and $f_{\rm esc,max}$
represent the 15.9 and 84.1 percentiles of the model $f_{\rm esc}$
predictions. For galaxies with lower limits on O32, the $f_{\rm esc}$
predictions are also lower limits and are marked accordingly. The Reference
column lists the publication used for the model input variables.

Table 7: Predictions for $z\gtrsim 6$ Galaxies from the ELG-O32-$\beta$ Model

Source ID | $z_{\rm spec}$ | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$ | Reference
---|---|---|---|---|---
21842 | 7.982 | 0.09 | 0.03 | 0.41 | Saxena et al. (2023a)
10013682 | 7.276 | 0.009 | 0.00 | 0.04 | Saxena et al. (2023a)
16625 | 6.631 | 0.37 | 0.12 | 0.68 | Saxena et al. (2023a)
18846 | 6.336 | 0.60 | 0.37 | 1.00 | Saxena et al. (2023a)
19342 | 5.974 | 0.47 | 0.18 | 0.82 | Saxena et al. (2023a)
9422 | 5.937 | 0.85 | 0.59 | 1.00 | Saxena et al. (2023a)
6002 | 5.937 | 0.11 | 0.03 | 0.44 | Saxena et al. (2023a)
12637 | 7.66 | 0.21 | 0.07 | 0.60 | Saxena et al. (2023a)
15362 | 6.794 | 0.00 | 0.00 | 0.03 | Saxena et al. (2023a)
13607 | 6.622 | 0.00 | 0.00 | 0.00 | Saxena et al. (2023a)
14123 | 6.327 | 0.07 | 0.03 | 0.34 | Saxena et al. (2023a)
58850 | 6.263 | 0.35 | 0.12 | 0.66 | Saxena et al. (2023a)
17138 | 6.204 | 0.00 | 0.00 | 0.02 | Saxena et al. (2023a)
9365 | 5.917 | 0.28 | 0.09 | 0.61 | Saxena et al. (2023a)
11027 | 9.51 | 0.01 | 0.00 | 0.05 | Williams et al. (2023)
4590 | 8.495 | 0.01 | 0.00 | 0.05 | Schaerer et al. (2022)
6355 | 7.664 | 0.06 | 0.03 | 0.35 | Schaerer et al. (2022)
10612 | 7.66 | 0.17 | 0.05 | 0.53 | Schaerer et al. (2022)
JADES-GS-z7-LA | 7.278 | 0.00 | 0.00 | 0.02 | Saxena et al. (2023b)
10025 | 7.875 | 0.03 | 0.005 | 0.11 | Mascia et al. (2023a)
100004 | 7.884 |  $>$0.02 |  $>$0.005 |  $>$0.11 | Mascia et al. (2023a)
10000 | 7.884 | 0.09 | 0.03 | 0.40 | Mascia et al. (2023a)
10021 | 7.288 | 0.53 | 0.24 | 0.88 | Mascia et al. (2023a)
100001 | 7.875 | 0.01 | 0.00 | 0.05 | Mascia et al. (2023a)
100003 | 7.88 | 0.58 | 0.30 | 1.00 | Mascia et al. (2023a)
100005 | 7.883 | 0.03 | 0.008 | 0.16 | Mascia et al. (2023a)
150008 | 6.23 |  $>$0.18 |  $>$0.05 |  $>$0.58 | Mascia et al. (2023a)

Note. — Predicted $f_{\rm esc}$ values from the ELG-O32-$\beta$ model. $z_{\rm
spec}$ is the spectroscopic redshift. $f_{\rm esc,min}$ and $f_{\rm esc,max}$
represent the 15.9 and 84.1 percentiles of the model $f_{\rm esc}$
predictions. For galaxies with lower limits on O32, the $f_{\rm esc}$
predictions are also lower limits and are marked accordingly. The Reference
column lists the publication used for the model input variables.

While we have generated initial predictions for $z\gtrsim 6$ galaxies using
our models, predictions for the largest samples are limited by a lack of
spectroscopic information. We do find evidence of high $f_{\rm esc}$ $>0.2$
among smaller spectroscopic samples at $z\gtrsim 6$, but these samples may be
biased toward stronger emission-line galaxies. To obtain more accurate
predictions, high-redshift studies should prioritize observations of nebular
line ratios and galaxy sizes for larger samples. As JWST surveys continue,
larger, more representative spectroscopic samples will improve estimates of
the average $f_{\rm esc}$ at $z>6$.

###  4.2 Trends with UV Magnitude

We examine the predicted $f_{\rm esc}$ values as a function of magnitude in
Figure 4. We plot predictions from the three models that apply to the largest
sample sizes (ELG-EW, ELG-O32-$\beta$, and R50-$\beta$) and from the model
that applies to the faintest galaxies ($\beta$-Metals; Table 8). We see no
strong trend with magnitude, although the few brightest sources with $M_{\rm
1500}\lesssim-21.5$ do appear to have higher $f_{\rm esc}$ on average, with a
median $f_{\rm esc}$ $=0.18$ vs. 0.07 for fainter sources. Aside from these
rare sources, the ELG-EW, ELG-O32-$\beta$, and $\beta$-Metals samples tend to
show $f_{\rm esc}$ within $\sim 0.02-0.3$ across a wide range of magnitudes,
from -16 to -22. The R50-$\beta$ model (Figure 4b) also shows a flat
distribution of $f_{\rm esc}$ with magnitude, but it predicts many more
galaxies at higher $f_{\rm esc}$. However, as previously noted, 41% of the
galaxies in this model require extrapolating predictions to galaxies with
extremely compact morphologies, which makes these predictions uncertain.

Table 8: Predictions for $z\gtrsim 6$ Galaxies from the $\beta$-Metals Model

Source ID | $z_{\rm spec}$ | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$ | Reference
---|---|---|---|---|---
18924 | 7.7 | 0.05 | 0.02 | 0.30 | Atek et al. (2023)
16155 | 6.87 | 0.03 | 0.006 | 0.18 | Atek et al. (2023)
23920 | 6.00 | 0.10 | 0.03 | 0.46 | Atek et al. (2023)
12899 | 6.88 | 0.17 | 0.04 | 0.55 | Atek et al. (2023)
8613 | 6.38 | 0.18 | 0.04 | 0.57 | Atek et al. (2023)
23619 | 6.72 | 0.08 | 0.03 | 0.42 | Atek et al. (2023)
38335 | 6.23 | 0.02 | 0.00 | 0.11 | Atek et al. (2023)
27335 | 6.76 | 0.08 | 0.03 | 0.42 | Atek et al. (2023)
11027 | 9.51 | 0.03 | 0.005 | 0.15 | Williams et al. (2023)
4590 | 8.495 | 0.05 | 0.02 | 0.35 | Schaerer et al. (2022)
6355 | 7.664 | 0.03 | 0.006 | 0.18 | Schaerer et al. (2022)
10612 | 7.66 | 0.05 | 0.02 | 0.36 | Schaerer et al. (2022)

Note. — Predicted $f_{\rm esc}$ values from the $\beta$-Metals model. $z_{\rm
spec}$ is the spectroscopic redshift. $f_{\rm esc,min}$ and $f_{\rm esc,max}$
represent the 15.9 and 84.1 percentiles of the model $f_{\rm esc}$
predictions. The Reference column lists the publication used for the model
input variables.

Any trends with magnitude or lack thereof in the other models should also be
taken with caution. All the models have a tendency to under-predict $f_{\rm
esc}$, but this tendency is most pronounced in the ELG-EW and $\beta$-Metals
models (Figure 2), which trace the extremes of the magnitude range in Figure
4. The more reliable ELG-O32-$\beta$ model, which better reproduces the LzLCS+
$f_{\rm esc}$ values, suggests that some of the galaxies in the $M_{\rm
1500}=-19$ to $-21$ magnitude range may have significantly high $f_{\rm esc}$.
The ELG-O32-$\beta$ sample does not include enough faint galaxies to reveal
whether or not a trend between $f_{\rm esc}$ and magnitude exists, however.

At the faintest magnitudes, the Atek et al. (2023) sample, included in the
$\beta$-Metals model includes several extremely low-luminosity lensed galaxies
observed as part of the JWST UNCOVER survey (Bezanson et al., 2022; Weaver et
al., 2024). For these galaxies, the $\beta$-Metals model therefore
extrapolates the LzLCS+ sample to an unobserved parameter space. The only
statistically significant predictor variable in the $\beta$-Metals model is
$\beta_{\rm 1550}$, and Figure 4 illustrates that these galaxies’ blue slopes
suggest at least moderate $f_{\rm esc}$ (0.02-0.18; Table 8). Crucially, even
these moderate $f_{\rm esc}$ values, which may be under-predicted, are still
greater than what is required for such faint galaxies to dominate reionization
given their high LyC photon production rates (Atek et al., 2023).

While moderate or high $f_{\rm esc}$ for faint, blue galaxies matches past
predictions (e.g., Chisholm et al., 2022), the high $f_{\rm esc}$ in the
brighter sources in Figure 4 may seem surprising at first glance. In all four
models, the fitted coefficients (see the Appendix) show that a brighter
$M_{\rm 1500}$ results in a higher $f_{\rm esc}$ at fixed values of the other
inputs (fixed E(B-V), fixed O32, etc). This effect does not come from the
association of higher observed UV luminosity with lower dust content. The data
do not show such a correlation, and, because the models already include dust
measurements, the model implies that brighter UV luminosities increase $f_{\rm
esc}$ at fixed dust attenuation. In addition, we see a similar relationship
with $M_{*}$, where higher $M_{*}$ is also associated with higher $f_{\rm
esc}$ at fixed values of the other parameters. Lin et al. (2024) find a
similar result using a different technique. Using the LzLCS+ data, they fit a
logistic regression model for the probability of having LyC escape as a
function of $M_{\rm 1500}$, $\beta_{\rm 1550}$, and O32 and likewise find a
brighter $M_{\rm 1500}$ increases the probability of LyC escape for fixed
values of the other variables. This scaling does not imply that higher UV
luminosities correlate with $f_{\rm esc}$ in the LzLCS+ sample; if anything,
the LzLCS+ indicates a slight trend in the opposite direction (Flury et al.,
2022b). Instead, this scaling implies that a bright galaxy with LCE-like
properties, such as high O32 and blue UV slope, is a more extreme object than
a faint galaxy with identical properties. For instance, at fixed O32, the
LzLCS+ galaxies with higher luminosities tend to have higher $f_{\rm esc}$,
but high O32 values are also uncommon among the bright galaxies as a whole,
leading to an overall trend of decreasing $f_{\rm esc}$ with luminosity (e.g.,
Flury et al., 2022b).

The enhanced $f_{\rm esc}$ for bright LCE candidates affects the predictions
in Figure 4 by boosting the $f_{\rm esc}$ of a bright, dust-poor emission-line
galaxy or compact galaxy relative to an otherwise similar faint galaxy. For
instance, for the ELG-EW model, the brightest galaxies should have higher
predicted $f_{\rm esc}$ only if they have equally low E(B-V) and equally high
EW as the fainter galaxies. Indeed, the ELG-EW samples, mostly galaxies from
Endsley et al. (2021, 2023), show no trend between E(B-V)UV and $M_{\rm 1500}$
or between EW([O iii]+H$\beta$) and $M_{\rm 1500}$; with comparable dust and
emission line strengths, the brighter galaxies therefore end up with higher
predicted $f_{\rm esc}$. Endsley et al. (2023) note that lower metallicities
or higher $f_{\rm esc}$ among the fainter galaxies could suppress their EWs
and account for the lack of an observed trend with $M_{\rm 1500}$. Hence,
other properties, such as O32, might be necessary to distinguish faint
galaxies with high $f_{\rm esc}$. The EWs in the Endsley et al. (2021, 2023)
samples also come from photometry, which may be more uncertain than
spectroscopic measurements (e.g., Duan et al., 2023).

The ELG-O32-$\beta$ samples have O32 measurements, but they similarly show no
strong trend between $\beta_{\rm 1550}$ and $M_{\rm 1500}$ or between O32 and
$M_{\rm 1500}$. Likewise, the Morishita et al. (2023) sample, which
constitutes the bulk of the R50-$\beta$ model sample, exhibits a relatively
flat trend between radius and $M_{\rm 1500}$. Only the $\beta$-Metals samples,
mostly from Atek et al. (2023), show any trend between physical properties and
$M_{\rm 1500}$, with the lowest metallicities and bluest slopes appearing
among the faintest galaxies.

The limitations of the models and the available $z\gtrsim 6$ datasets make it
difficult to discern any trends between $M_{\rm 1500}$ and $f_{\rm esc}$.
Limited measurements for galaxies at the brightest and faintest magnitudes
restrict us to some of the less accurate Cox models. The models also require
extrapolating to make predictions for galaxies outside the magnitude or radius
range of the LzLCS+ sample. Better constraints on $f_{\rm esc}$ within the
$z>6$ population will require larger spectroscopic samples across a wide range
of galaxy magnitudes as well as better observational constraints on $f_{\rm
esc}$ in a more diverse set of galaxies at lower redshifts.

Figure 3: The distribution of predicted $f_{\rm esc}$ from the ELG-EW model
(a) for $z\gtrsim 6$ galaxies from Endsley et al. (2021, 2023), Bouwens et al.
(2023), Tang et al. (2023), Fujimoto et al. (2023), and Saxena et al. (2023b)
and from the R50-$\beta$ model (b) for $z\gtrsim 6$ galaxies from Morishita et
al. (2023) and Mascia et al. (2023a). The vertical dotted lines show the
median and quartiles of the distributions ($f_{\rm esc}$ $=0.027,0.047,0.105$
for ELG-EW and $f_{\rm esc}$ $=0.039,0.141,0.467$ for R50-$\beta$). Both
models have some limitations but suggest that many, but not all, high-redshift
galaxies may have $f_{\rm esc}$ $>0.1$.  Figure 4: Predicted $f_{\rm esc}$ as
a function of $M_{\rm 1500}$ for $z\gtrsim 6$ galaxies. Panel a shows the
predicted $f_{\rm esc}$ from the ELG-EW model (circles) for galaxies with
photometric redshifts (gray) and spectroscopic redshifts (black), from the
ELG-O32-$\beta$ model (blue stars), and from the $\beta$-Metals model (green
diamonds). Panel b shows the predicted $f_{\rm esc}$ from the R50-$\beta$
model (squares) with photometric (gray) and spectroscopic (black) redshifts.
The crosses in the upper corners show the median uncertainties for each panel.
Dark gray dotted lines connect the predictions for the same galaxy in
different models. We see no strong trend between $f_{\rm esc}$ and $M_{\rm
1500}$, but many of the predictions suffer from a limited set of variables
(ELG-EW and $\beta$-Metals models) or require extrapolation outside the LzLCS+
parameter space (R50-$\beta$ model).

###  4.3 Notes on Individual Sources

####  4.3.1 Strong LCES

Several galaxies in the ELG-O32 and ELG-O32-$\beta$ models appear to be
extremely strong LCEs (Tables 6 and 7). One unusual galaxy, CEERS-1019, has a
predicted $f_{\rm esc}$ $=1$. Such a high $f_{\rm esc}$ would seemingly
conflict with the presence of strong emission lines, since it should
correspond to a complete absence of absorbing gas. However, given the
uncertainty in the model predictions, the $f_{\rm esc}$ of CEERS-1019 could be
as low as 0.6. Two other galaxies in the ELG-O32 model sample, CEERS-698 and
CEERS-1027, have $f_{\rm esc}$ $\geq 0.5$, and $f_{\rm esc}$ as high as 1 is
consistent with their uncertainties. Interestingly, all three of these strong
LCE candidates are strong Ly$\alpha$ emitters, which suggests that they reside
within ionized bubbles (Tang et al., 2023). Similarly, four galaxies in the
ELG-O32-$\beta$ samples have predicted $f_{\rm esc}$ $>0.5$: GLASS-JWST 10021
and 100003 from Mascia et al. (2023a) and 18846 and 9422 from JADES (Saxena et
al., 2023a). With O32$=70.6$, the ELG-O32-$\beta$ model infers $f_{\rm esc}$
$=0.85$ for 9422, much higher than the $f_{\rm esc}$ $=0.01$ value derived by
Saxena et al. (2023a) using the multivariate prediction method of Choustikov
et al. (2023); we discuss the differences between our predictions and those of
other models in §5.

The true $f_{\rm esc}$ of some of these extreme galaxies may not be quite as
high as the models predict, however. Because two of these strong LCEs,
CEERS-1019 ($f_{\rm esc}$ $=1$) and CEERS-698 ($f_{\rm esc}$ $\sim 0.7$), have
brighter UV magnitudes than the LzLCS+ sample by 0.4-0.9, their high predicted
$f_{\rm esc}$ could be unreliable. However, these same models successfully
predicted the high $f_{\rm esc}$ of the $z\sim 3$ galaxies Ion2 and J1316+2614
(Figure 2), which are as bright as these $z>6$ galaxies. The LzLCS+ sample
itself also contains only three LCEs with $f_{\rm esc}$ $>0.5$, which limits
the models’ ability to accurately determine predictive relationships in the
high $f_{\rm esc}$ regime (Mascia et al., 2023a). CEERS-1019 may differ from
the LzLCS+ galaxies for additional reasons; it has unusually strong nitrogen
emission and is a possible candidate for supermassive star formation (Marques-
Chaves et al., 2023) or an AGN (Larson et al., 2023). The physical conditions
occurring in extreme galaxies like CEERS-1019 likely are not present within
the LzLCS+ sample. Empirical predictions based on $z=0.3$ galaxies would not
be suitable for such objects.

The predictions for the other strong LCE candidates are likely more reliable,
provided their star-formation and nebular conditions are not significantly
different from the LCEs at $z\sim 3$. Galaxy 9422 has a higher O32 ratio than
the LzLCS+, but the other strong leaker candidates fall within the LzLCS+
parameter space, so that the model is not extrapolating to an unobserved
regime. Hence, our models suggest that several strong emission-line galaxies
identified in the epoch of reionization may be strong LCEs with $f_{\rm esc}$
well above 0.2.

####  4.3.2 JADES-GS-z7-LA

Saxena et al. (2023b) report the detection of a $z=7.3$ galaxy with strong
Ly$\alpha$ emission, JADES-GS-z7-LA, which is presumably located within an
ionized bubble. With $f_{{\rm esc,Ly}\alpha}$$=0.96$, IGM absorption seems to
have had little effect on the galaxy’s Ly$\alpha$ emission. We therefore
choose to incorporate the measured EW(Ly$\alpha$) and $f_{{\rm esc,Ly}\alpha}$
in our model predictions. We can apply most of the models in Table 1 to JADES-
GS-z7-LA, omitting only the TopThree, R50-$\beta$, and $\beta$-Metals models
due to a lack of reported $\Sigma_{\rm SFR}$, radius, and 12+$\log_{\rm
10}$(O/H). We compare the different $f_{\rm esc}$ predictions in Table 9.
JADES-GS-z7-LA does fall outside of the LzLCS+ parameter space, which may
limit the accuracy of the $f_{\rm esc}$ predictions; it is fainter than the
LzLCS+ galaxies by 1.6 mag, and it has a higher EW(Ly$\alpha$) by 144.1 Å.

The different Cox models disagree regarding the $f_{\rm esc}$ of JADES-
GS-z7-LA, with predicted $f_{\rm esc}$ ranging from 0 to 0.40. The only
nonzero predictions come from models that have information about UV dust
attenuation and Ly$\alpha$, but even the extreme Ly$\alpha$ properties of
JADES-GS-z7-LA do not guarantee $f_{\rm esc}$ $>0.1$. The LAE-O32 and
ELG-O32-$\beta$-Ly$\alpha$ models include O32 and stellar mass as well as
Ly$\alpha$ and find $f_{\rm esc}$ $=0.017-0.085$. These two models are also
among the best-performing models for the LzLCS+ sample (Table 3). JADES-
GS-z7-LA is not devoid of dust, with E(B-V)${}_{\rm UV}=0.10$, and its UV
slope $\beta=-2.1$ is not extreme. These properties are consistent with the
properties of $z\sim 0.3$ galaxies with moderate $f_{\rm esc}$; all the LzLCS+
galaxies with higher $f_{\rm esc}$ $>0.1$ have lower E(B-V)${}_{\rm UV}<0.1$,
and a slope of $\beta_{\rm 1550}=-2.1$ matches the median value for the
moderate LCEs with $f_{\rm esc}$ $=0.05-0.1$.

Furthermore, as discussed above, given JADES-GS-z7-LA’s low mass ($\log_{\rm
10}$($M_{*}/$ M☉)=$7.15$), its O32 value of 8.8 may not be extreme. Even its
strong Ly$\alpha$ emission does not necessarily imply a low optical depth
along the line of sight, since Ly$\alpha$ photons can scatter. For instance,
the non-leaking $z\sim 0.3$ galaxy J1248+4259 (Izotov et al., 2018b) has
$f_{\rm esc}$ $\leq 0.013$ and EW(Ly$\alpha$)$=256$ Å. One piece of evidence
that favors a high $f_{\rm esc}$ for JADES-GS-z7-LA is the small offset
between the Ly$\alpha$ velocity and the systemic redshift (Saxena et al.,
2023b). With the low resolution of most of the LzLCS+ UV spectra, we cannot
currently include this parameter in our models, and it may boost the predicted
$f_{\rm esc}$ for JADES-GS-z7-LA. However, this low velocity offset could
still be consistent with little to no LyC escape; its offset of $120\pm 80$
${\rm\thinspace km}{\rm\thinspace s}^{-1}\,$ also resembles the Ly$\alpha$
profile of the non-leaker J1248+4259, whose red peak is offset from the
systemic velocity by $<200$ ${\rm\thinspace km}{\rm\thinspace s}^{-1}\,$.

In agreement with the assessment of Saxena et al. (2023b), we thus find that
the evidence for LyC escape in JADES-GS-z7-LA is ambiguous. The multivariate
models that incorporate the most information suggest that it likely has a
moderate $f_{\rm esc}$ $\sim 0.017-0.085$ and may therefore not be solely
responsible for producing its ionized bubble in the IGM (e.g., Witstok et al.,
2023).

####  4.3.3 Galaxies with High O32

The case of JADES-GS-z7-LA highlights the fact that the Cox multivariate
predictions can differ from single-variable estimates (e.g., based on
Ly$\alpha$ alone). Galaxies with high O32 provide another example of this
result. For instance, Williams et al. (2023) present JWST observations of a
lensed low-mass ($\log_{\rm 10}$($M_{*}$/ M☉) = 7.7) galaxy, 11027, at
$z=9.51$ and hypothesize that the galaxy has a high $f_{\rm esc}$ $>0.1$ based
on its high O32 ratio of 12 ([O iii]$\lambda\lambda$5007+4959/[O ii]=16).
However Lin et al. (2024) suggest that 11027 is not likely to be a strong LCE
considering the combination of O32, UV magnitude, and UV slope. Our ELG-O32
and ELG-O32-$\beta$ Cox models (Tables 6-7) likewise predict $f_{\rm esc}$ of
only 0.012-0.021. O32 alone is not sufficient to constrain $f_{\rm esc}$, and
high O32 values are common among the lowest mass galaxies in the LzLCS+ sample
(Flury et al., 2022b), including weak LyC emitters. In contrast, the
ELG-O32-$\beta$ model predicts that several galaxies with lower O32 than
galaxy 11027 (e.g., 6355, 10000, 10612, 12637 with O32=6.3-10.6, Schaerer et
al. 2022; Mascia et al. 2023a; Saxena et al. 2023a) are actually more likely
to be LCEs because of their higher luminosities and/or bluer UV slopes (Table
7). The Cox predictions emphasize the fact that not all galaxies with high O32
have high $f_{\rm esc}$ (e.g., Izotov et al., 2018b; Flury et al., 2022b) and
additional properties such as mass, UV luminosity, and dust extinction are
important to consider.

Table 9: $f_{\rm esc}$ Predictions for JADES-GS-z7-LA

Model | $f_{\rm esc}$ | $f_{\rm esc,min}$ | $f_{\rm esc,max}$
---|---|---|---
LAE | 0.400 | 0.142 | 0.736
LAE-O32 | 0.085 | 0.028 | 0.325
LAE-O32-nodust | 0 | 0 | 0.028
ELG-EW | 0 | 0 | 0.026
ELG-O32 | 0 | 0 | 0.012
ELG-O32-$\beta$ | 0 | 0 | 0.022
ELG-O32-$\beta$-Ly$\alpha$ | 0.017 | 0 | 0.059

Note. — Predicted $f_{\rm esc}$ values for JADES-GS-z7-LA from different
models. $f_{\rm esc,min}$ and $f_{\rm esc,max}$ represent the 15.9 and 84.1
percentiles of the model $f_{\rm esc}$ predictions.

Multivariate models are an important tool for predicting $f_{\rm esc}$ in the
epoch of reionization. Our results highlight the fact that multivariate models
can give substantially different predictions from single variable estimates,
and $f_{\rm esc}$ predictions should therefore incorporate as much information
as possible.

##  5 Comparison with Other Studies

###  5.1 Single-Variable Predictions

Among the observable properties measured for the LzLCS+ sample, the UV slope
$\beta_{\rm 1550}$ shows one of the strongest correlations with $f_{\rm esc}$.
Consequently, Chisholm et al. (2022) propose that $\beta_{\rm 1550}$ could
predict $f_{\rm esc}$ for high-redshift galaxies and derive a relationship
between $\beta_{\rm 1550}$ and $f_{\rm esc}$ based on the LzLCS+ dataset. In
Figure 5, we compare the predicted $f_{\rm esc}$ from this method with the
observed $f_{\rm esc}$ for the LzLCS+ and for $z\sim 3$ galaxies (Vanzella et
al., 2016; Bassett et al., 2019; Marques-Chaves et al., 2022). We list the
goodness-of-fit statistics for the LzLCS+ with the Chisholm et al. (2022)
model in Table 10. The Chisholm et al. (2022) model’s high $R^{2}=0.45$ and
low RMS$=0.43$ are comparable to the best-performing Cox models in Tables 3
and 10, which demonstrates that $\beta_{\rm 1550}$ alone can indeed predict
$f_{\rm esc}$ for LCEs reasonably well. Consistent with this result, Paper I
finds that, of the variables accessible at high redshift, $\beta_{\rm 1550}$
is the most important variable to include in the Cox models.

Despite its success for LCEs, however, the Chisholm et al. (2022) model has
more difficulty in constraining $f_{\rm esc}$ in the LzLCS+ non-detections, as
indicated by its lower concordance $C=0.76$. This concordance is lower than
all the Cox models in Tables 3 and 10, which have $C=0.77-0.83$. The
difference between the predictions for detections and non-detections is
apparent in Figure 5. For LzLCS+ galaxies with $f_{\rm esc}$ detections or
upper limits between $0.01-0.1$, the one-to-one relation between the Chisholm
et al. (2022) model predictions and the observations runs right between the
LzLCS+ detections. However, nearly all the non-detections are above this line,
indicating that the model is systematically over-predicting their $f_{\rm
esc}$. In contrast, for one of the better performing Cox models, such as the
ELG-O32 model in Figure 2b, both detections and non-detections in this same
$f_{\rm esc}$ range fall on either side of the one-to-one relation, indicating
that predictions for detections and non-detections are comparable.

Figure 5 also applies the Chisholm et al. (2022) model to two $z\sim 3$ LCEs:
Ion2 and J1316+2614. The model correctly identifies both galaxies as LCEs but
under-predicts their $f_{\rm esc}$, although Ion2’s observed $f_{\rm esc}$ is
consistent with the stated uncertainties in the Chisholm et al. (2022)
predictions. Comparing the Ion2 predictions with those from the Cox models
discussed in §3, we find that the Chisholm et al. (2022) single-variable
$\beta_{\rm 1550}$ model is one of the least accurate at predicting Ion2’s
$f_{\rm esc}$, with a predicted $f_{\rm esc}$ $=0.29$, corresponding to
RMS=$0.40$ dex. As discussed in Chisholm et al. (2022), $\beta_{\rm 1550}$
only tracks the loss of LyC photons due to dust. Scatter below and above the
Chisholm et al. (2022) $\beta_{\rm 1550}$-$f_{\rm esc}$ relation results from
variations in the absorbing H i column. Other properties, such as O32,
$\Sigma_{\rm SFR}$, and UV magnitude, may better track the H i component of
LyC absorption for galaxies like Ion2.

Figure 6 illustrates how including these other variables affects predictions
at $z>4$. We plot the difference between the observed LzLCS+ $f_{\rm esc}$ and
the Chisholm et al. (2022) predictions as a function of $\beta_{\rm 1550}$.
Galaxies are color-coded by their O32 ratio (Figure 6a) or $\Sigma_{\rm SFR}$
(Figure 6b). For the reddest $\beta_{\rm 1550}$, all the model predictions and
observations agree that galaxies have little to no LyC escape. However, at
blue UV slopes, the disagreement increases. The Chisholm et al. (2022) model
assigns all galaxies a single value of $f_{\rm esc}$, whereas the LzLCS+
observations show a range of $f_{\rm esc}$ and deviate from the predictions by
as much as $\Delta$$f_{\rm esc}$ $=0.2$ lower or $\Delta$$f_{\rm esc}$ $=0.6$
higher at the bluest slopes. For the LzLCS+ galaxies, this difference is
associated with other differences in the galaxies’ properties; galaxies with
higher O32 or higher $\Sigma_{\rm SFR}$ tend to be offset to higher $f_{\rm
esc}$, and galaxies with lower $f_{\rm esc}$ than the Chisholm et al. (2022)
predictions tend to have lower O32 or lower $\Sigma_{\rm SFR}$.

The multivariate model predictions reflect these trends with other properties
and assign galaxies a higher or lower $f_{\rm esc}$ at fixed $\beta_{\rm
1550}$ accordingly. In Figure 6, we also plot the difference between the
predicted $f_{\rm esc}$ for $z>4$ galaxies from multivariate models vs. the
single-variable Chisholm et al. (2022) model. All these multivariate models
use $\beta_{\rm 1550}$ in addition to two or more other variables. The Mascia
et al. (2023a) model (squares) predicts $f_{\rm esc}$ using a linear fit to
$\beta_{\rm 1550}$, $\log_{\rm 10}$(O32) and $\log_{\rm 10}$($r_{\rm 50,NUV}$)
for their sample of $z>4$ JWST-GLASS galaxies. We show predictions from two
Cox models in crosses and stars: a Cox model using the three top-ranked, high-
redshift accessible variables from Paper I ($\beta_{\rm 1550}$, O32, and
$\Sigma_{\rm SFR}$) for the Mascia et al. (2023a) sample and the
ELG-O32-$\beta$ Model (§§3-4), which uses $\beta_{\rm 1550}$, O32, and $M_{\rm
1500}$ as inputs to predict $f_{\rm esc}$ for $z>4$ galaxies from Williams et
al. (2023), Schaerer et al. (2022), Mascia et al. (2023a), and Saxena et al.
(2023a). Like the LzLCS+ observations, the multivariate $z>4$ predictions can
significantly deviate from the Chisholm et al. (2022) predictions, especially
at blue UV slopes. High-redshift surveys relying only on $\beta_{\rm 1550}$
could therefore miss some of the physical differences among blue galaxies that
might indicate variations in $f_{\rm esc}$. Properties such as O32 and
$\Sigma_{\rm SFR}$ are sensitive to feedback in galaxies and may complement
$\beta_{\rm 1550}$ measurements by tracing a galaxy’s ability to carve low
column density H i channels.

Figure 5: The $f_{\rm esc}$ predictions for the LzLCS+ and $z\sim 3$ galaxies
from the Chisholm et al. (2022) model, which uses a single predictor variable:
$\beta_{\rm 1550}$. Symbols are the same as in Figure 1. The model reproduces
the $f_{\rm esc}$ of LCEs but over-predicts the $f_{\rm esc}$ for the LzLCS+
non-detections and under-predicts $f_{\rm esc}$ for the strongest LCEs.  Table
10: Goodness-of-Fit Statistics for Literature vs. Cox Models

Model | $R^{2}$ | $R^{2}_{\rm adj}$ | RMS | $C$
---|---|---|---|---
Chisholm et al. (2022) | 0.45 | 0.43 | 0.43 | 0.76
Choustikov et al. (2023) | -1.10 | -1.40 | 0.83 | 0.57
Mascia et al. (2023a) | 0.57 | 0.54 | 0.38 | 0.81
JWST Cox Model (Paper I) | 0.29 | 0.14 | 0.47 | 0.83
Cox Model with $\beta_{\rm 1550}$, $\Sigma_{\rm SFR}$, O32 | 0.34 | 0.29 | 0.46 | 0.83
Cox Model with Choustikov et al. (2023) Variables | 0.28 | 0.16 | 0.48 | 0.82
Cox Model with Mascia et al. (2023a) Variables | 0.40 | 0.35 | 0.44 | 0.82

Note. — A comparison of the goodness-of-fit statistics for $f_{\rm esc}$
predictions for the LzLCS+ sample from the Chisholm et al. (2022), Choustikov
et al. (2023), and Mascia et al. (2023a) literature models and various Cox
models. We list statistics for the JWST Model from Paper I and a model limited
to its three top-ranked variables. We also list statistics for Cox models run
using the same sets of variables as in the Choustikov et al. (2023) and Mascia
et al. (2023a) models.

Figure 6: Circles (detections) and triangles (upper limits) show the
difference between LzLCS+ observations and Chisholm et al. (2022) predicted
$f_{\rm esc}$ as a function of observed $\beta_{\rm 1550}$. Other symbols show
the difference between multivariate model predictions for $z>4$ samples and
the Chisholm et al. (2022) single-variable predictions as a function of
$\beta_{\rm 1550}$. Squares compare the $f_{\rm esc}$ predictions from Mascia
et al. (2023a), using $\beta_{\rm 1550}$, O32, and $r_{\rm 50,NUV}$. Crosses
show $f_{\rm esc}$ predictions for the same $z>4$ Mascia et al. (2023a) sample
for a Cox model using $\beta_{\rm 1550}$, O32, and $\Sigma_{\rm SFR}$. Stars
show the ELG-O32-$\beta$ Cox Model predictions for $z>4$ galaxies from
Williams et al. (2023), Schaerer et al. (2022), Mascia et al. (2023a), and
Saxena et al. (2023a), where the input variables are $\beta_{\rm 1550}$, O32,
and $M_{\rm 1500}$. We color-code galaxies by their observed O32 ratio in
panel (a) and by $\Sigma_{\rm SFR}$ in panel (b). At blue UV slopes, the
LzLCS+ observations and high-redshift model predictions can differ
significantly from the Chisholm et al. (2022) model predictions, depending on
the galaxies’ other properties, such as O32 and $\Sigma_{\rm SFR}$.

###  5.2 Multivariate Predictions

Recently, Choustikov et al. (2023) and Mascia et al. (2023a) have developed
multivariate linear regression models to predict $f_{\rm esc}$ from a set of
observable variables. Choustikov et al. (2023) derive their model for $f_{\rm
esc}$ from synthetic spectra of $z>4$ galaxies from the SPHINX cosmological
radiation hydrodynamics simulation. Mascia et al. (2023a) base their
predictions on the LzLCS+ sample but adopt the upper limit in $f_{\rm esc}$ as
the observed value for non-detections. Here, we use galaxies at low and high
redshift to compare the predictions of the Cox models with the predictions
from these literature models.

In Figure 7a, we show the results of applying the Choustikov et al. (2023)
model to the LzLCS+ sample, and we list the corresponding goodness-of-fit
statistics in Table 10. The Choustikov et al. (2023) model fails to predict
the observed $f_{\rm esc}$ in the LzLCS+ sample. This result does not depend
on our method of calculating $f_{\rm esc}$; we find a similarly poor fit if we
substitute the $f_{\rm esc}$ derived from H$\beta$, as used in the Choustikov
et al. (2023) simulations. In Figure 7b, we show $f_{\rm esc}$ predictions
from a Cox model run using the same set of variables ($\beta_{\rm 1550}$,
E(B-V)neb, $L$(H$\beta$), $M_{\rm 1500}$, R23$=$([O iii]
$\lambda\lambda$5007,4959+[O ii]$\lambda$3727)/H$\beta$, and O32) as in the
Choustikov et al. (2023) linear regression mode, and we list the goodness-of-
fit metrics in Table 10.

This model performs comparably to the other Cox models in Table 3, which shows
that the reason for the Choustikov et al. (2023) model’s difficulty is not its
set of variables. Rather, in the Choustikov et al. (2023) model, $f_{\rm esc}$
anti-correlates with O32, whereas the LzLCS+ data and the Cox model using the
Choustikov et al. (2023) variable set indicate a correlation between $f_{\rm
esc}$ and O32. The compact, lower-mass LCEs in the LzLCS do not have
counterparts in the SPHINX simulations; these compact LCEs often have high O32
ratios $\gtrsim 10$ (Flury et al., 2022b), whereas the SPHINX galaxies with
O32 $>10$ typically have lower values of $f_{\rm esc}$. Hence, while we find
that high O32 increases $f_{\rm esc}$ in our multivariate predictions,
Choustikov et al. (2023) find the opposite, that at fixed $\beta$ and E(B-V),
galaxies with older ages and lower O32 have higher $f_{\rm esc}$. Radiative
feedback with a turbulent gas structure may allow $f_{\rm esc}$ from high O32
galaxies at earlier ages than predicted by cosmological simulations (e.g.,
Kakiichi & Gronke, 2021; Kimm et al., 2019; Choustikov et al., 2023). However,
further studies of $z\sim 3$ galaxies are also necessary to test whether the
LzLCS population and the LzLCS-derived Cox models correctly describe the high-
redshift population.

Figure 7: (a) The $f_{\rm esc}$ predictions from the Choustikov et al. (2023)
literature model compared with the LzLCS+ observations. (b) The $f_{\rm esc}$
predictions from a Cox model run using the same variables as the Choustikov et
al. (2023) model: $\beta_{\rm 1550}$, E(B-V)neb, $L$(H$\beta$), $M_{\rm
1500}$, R23, and O32. Symbols are the same as in Figure 1. The Choustikov et
al. (2023) model does not reproduce the $f_{\rm esc}$ observations from the
LzLCS+, but a Cox model using the same input variables does recover the
observed $f_{\rm esc}$.

Like our Cox model predictions in §4, Choustikov et al. (2023) predict a low
$f_{\rm esc}$ $=0.03$ for JADES-GS-z7-LA, due to its moderate UV slope and
dust content. However, because of the different dependence of $f_{\rm esc}$ on
O32, the Choustikov et al. (2023) model predicts a dramatically lower $f_{\rm
esc}$ for CEERS-44, CEERS-698, CEERS-1019, and CEERS-1027 ($f_{\rm esc}$
$=0.006-0.1$) and for JADES 9422 ($f_{\rm esc}$ $=0.01$; Saxena et al. 2023a)
than our model predictions ($f_{\rm esc}$ $>0.4$ to 1). Choustikov et al.
(2023) point out that their predicted $f_{\rm esc}$ values are consistent with
$f_{{\rm esc,Ly}\alpha}$ $>$ $f_{\rm esc}$ as observed locally (e.g., Flury et
al., 2022b), but IGM effects may complicate any comparison with $f_{{\rm
esc,Ly}\alpha}$. As noted above (§4), the Cox model predictions for some of
these galaxies require extrapolating outside the LzLCS+ parameter space, and
their $f_{\rm esc}$ may not truly be as extreme as this model predicts. Still,
these galaxies remain plausible candidates for high $f_{\rm esc}$ given the
observed trends seen at $z\sim 0.3$ and $z\sim 3$ (§3). LyC observations of
galaxies with similar properties are necessary to better constrain their
$f_{\rm esc}$.

In contrast to the results for the Choustikov et al. (2023) model, the Mascia
et al. (2023a) model (Figure 8a and Table 10) reproduces $f_{\rm esc}$ for the
LzLCS+ galaxies reasonably well. This finding is unsurprising, since the
Mascia et al. (2023a) model is in fact derived from the LzLCS+ dataset. A Cox
model run using the same variables ($\beta_{\rm 1550}$, $r_{\rm 50,NUV}$, and
O32) performs comparably to the Mascia et al. (2023a) model (Figure 8b and
Table 10). The largest difference occurs for weak LCEs ($f_{\rm esc}$ $<0.05$)
and non-detections. In these cases, the Mascia et al. (2023a) model tends to
systematically over-predict $f_{\rm esc}$, whereas by incorporating
information from non-detections, Cox model predictions are more evenly
distributed above and below the observed values. However, we note that the
over-prediction of $f_{\rm esc}$ in the Mascia et al. (2023a) model is small,
$\Delta$$f_{\rm esc}$ of only a few percent.

Figure 8: (a) The $f_{\rm esc}$ predictions from the Mascia et al. (2023a)
literature model compared with the LzLCS+ observations. (b) The $f_{\rm esc}$
predictions from a Cox model run using the same variables as the Mascia et al.
(2023a) model: $\beta_{\rm 1550}$, $r_{\rm 50,NUV}$, and O32. Symbols are the
same as in Figure 1. Both models perform similarly, but the Cox model does
slightly better at predicting $f_{\rm esc}$ in weak LCEs and non-detections.

Despite its derivation from the same dataset, the Cox and Mascia et al.
(2023a) models give different predictions at high redshift. The Mascia et al.
(2023a) model predicts $f_{\rm esc}$ $=0.19$ for Ion2, the only $z\sim 3$ LCE
with the required measurements, whereas the Cox models more accurately match
its observed $f_{\rm esc}$ $>0.5$ (§3). This difference suggests that galaxy
luminosity, incorporated in the Cox models in either $\Sigma_{\rm SFR}$ or in
$M_{\rm 1500}$ may be important in reproducing $f_{\rm esc}$ for the strongest
LCEs.

In Figure 9, we compare the model predictions for GLASS-JWST galaxies (Treu et
al., 2022) at $z=4-8$ from Mascia et al. (2023a) with predictions from Cox
models using similar combinations of variables. For a Cox model using the same
variables as in Mascia et al. (2023a), the predictions from Mascia et al.
(2023a) and the corresponding Cox model track each other (Figure 9a), aside
from the aforementioned tendency of the Mascia et al. (2023a) model to give
higher predictions at low $f_{\rm esc}$. However, the Mascia et al. (2023a)
and Cox models disagree more strongly when $\Sigma_{\rm SFR}$ is used as a
variable in the Cox model instead of $r_{\rm 50,NUV}$ (Figure 9b). The
inclusion of $\Sigma_{\rm SFR}$ shifts some galaxies with compact sizes but
weak star formation to much lower $f_{\rm esc}$. We conclude that including
some measure of luminosity or SFR may be important to accurately identify non-
leakers and extreme LCEs like Ion2.

Figure 9: A comparison of $f_{\rm esc}$ predictions for high-redshift galaxies
from Cox models and the Mascia et al. (2023a) model. (a) A Cox model using
$\beta_{\rm 1550}$, $r_{\rm 50,NUV}$, and O32 gives similar predicted $f_{\rm
esc}$ as the Mascia et al. (2023a) model for $z>4$ galaxies from the GLASS-
JWST (black) and CEERS (blue) surveys (Mascia et al., 2023a, b). (b) If
$\Sigma_{\rm SFR}$ is substituted for $r_{\rm 50,NUV}$ in the Cox model, the
Cox model $f_{\rm esc}$ predictions for the GLASS galaxies deviate from those
of the Mascia et al. (2023a) model. (c) For the CEERS galaxies, the $f_{\rm
esc}$ predictions from the ELG-O32 Cox model disagree with the Mascia et al.
(2023b) predictions. (d) A Cox model using $M_{\rm 1500}$ in addition to the
Mascia et al. (2023a) variables begins to agree more closely with the ELG-O32
Cox model predictions. Adding a measurement of galaxy luminosity in the form
of $\Sigma_{\rm SFR}$ or $M_{\rm 1500}$ can significantly affect the $f_{\rm
esc}$ predictions.

Mascia et al. (2023b) also apply their multivariate model to galaxies from
CEERS and find significantly different results than the predictions from our
ELG-O32 Cox model for these galaxies (§4, Table 6). We compare these
predictions in Figure 9c. Again, the disagreement primarily arises from the
inclusion of luminosity in the ELG-O32 Cox model but not the Mascia et al.
(2023a) model. As with the GLASS sample, when we run a Cox model with the same
input variables ($\beta_{\rm 1550}$, $r_{\rm 50,NUV}$, and O32), our
predictions agree closely with the Mascia et al. (2023b) predictions (Figure
9a). Adding one additional variable, $M_{\rm 1500}$, as a measure of
luminosity begins to bring the predicted $f_{\rm esc}$ into agreement with the
ELG-O32 predictions (Figure 9d). The remaining disagreement between the Mascia
et al. (2023b) and ELG-O32 Cox model predictions is due to the use of radius
as a variable in the Mascia et al. (2023a) model and different estimates of
dust content and ionization from different publications (Mascia et al., 2023b;
Fujimoto et al., 2023; Tang et al., 2023).

Observational and theoretical studies agree that $f_{\rm esc}$ depends on
multiple physical parameters. However, studies have not yet reached a
consensus as to which parameters matter and how. Predictions for $f_{\rm esc}$
in $z\gtrsim 6$ galaxies from this work, Choustikov et al. (2023), Mascia et
al. (2023a) and Mascia et al. (2023b) disagree because of different adopted
scalings with O32 and with luminosity. Observationally testing and
distinguishing between these predictions with larger samples at $z<6$ will be
critical to reliably predict $f_{\rm esc}$ in the epoch of reionization.

##  6 Implications for Reionization

As demonstrated in §3, the empirical Cox models derived from the LzLCS+ show
promise as diagnostics of $f_{\rm esc}$ at high redshift. However, with the
limited measurements at $z>6$ available so far, the models do not decisively
show which galaxy populations dominate reionization (§4.2). Because a
combination of properties regulates a galaxy’s optical depth, we need
estimates of many factors to accurately predict $f_{\rm esc}$. Dust
attenuation, galaxy morphology, ionization, and UV luminosity all play a role
in promoting LyC escape. Based primarily on dust attenuation, we would expect
faint galaxies to be strong LCEs (e.g., Chisholm et al., 2022; Atek et al.,
2023). Nevertheless, at a fixed value of $\beta_{\rm 1550}$, $f_{\rm esc}$
still shows considerable spread in the LzLCS+ and in multivariate predictions
(Figure 6), spanning a range of up to $\sim 0.7$ in $f_{\rm esc}$. Without
estimates of other properties (e.g., $\Sigma_{\rm SFR}$, O32), we therefore
cannot determine whether faint galaxies have $f_{\rm esc}$ above or below the
average LzLCS+ value for the same UV color.

At the same time, we cannot yet rule out the contribution of more luminous
galaxies. Several galaxies with $M_{\rm 1500}\leq-20$ (e.g., CEERS 498, 1027,
698, 1019; JADES 18846; GLASS-JWST 100003, 10021) have the low dust content
and high ionization suggestive of extreme $f_{\rm esc}$ $\geq 0.5$. The LzLCS+
data seem to suggest that for the same O32 ratio and $\beta_{\rm 1500}$, a
more luminous galaxy will have higher $f_{\rm esc}$ (see also Lin et al.
2024). This effect could possibly be an observational bias, as the LyC flux
should be easier to detect for more luminous objects. However, we note that at
fixed O32, brighter LzLCS+ galaxies show higher ratios of 900 Å to 1100 Å
flux, which suggests that their high $f_{\rm esc}$ is genuine. Moreover, the
trend of higher $f_{\rm esc}$ with luminosity only appears in galaxies that
have properties associated with LyC escape, such as high O32. Overall, the
brighter galaxies in the LzLCS+ have fewer LyC detections (Flury et al.,
2022b), which suggests that the detections are not biased toward the brightest
galaxies. Furthermore, including UV luminosity in the multivariate models
appears necessary to reproduce the $f_{\rm esc}$ of Ion2 (§3 and §5), which
suggests that the increase of $f_{\rm esc}$ with luminosity may be a real
phenomenon.

In the multivariate models, $M_{\rm 1500}$ might help to break the degeneracy
between ionization parameter and optical depth in galaxies with high O32. For
instance, low-metallicity galaxies tend to have both low luminosities and
inherently high ionization parameters (e.g., Nagao et al., 2006; Dopita et
al., 2006), such that they may be more likely to have high O32 even at high
optical depth. Alternatively, LCEs with higher UV luminosities could be able
to ionize gas over a wider opening angle, so that we are statistically more
likely to observe a high $f_{\rm esc}$ line of sight.

Of course, to evaluate a galaxy population’s influence on reionization, we
must also consider how many ionizing photons are produced and the total LyC
input into the IGM, not just $f_{\rm esc}$. Unfortunately, the Cox models do
not do as well at predicting the ionizing to non-ionizing UV flux ratio or the
total LyC luminosity (Paper I). Galaxies with high O32 and high nebular EWs
may possess elevated $\xi_{\rm ion}$ values (e.g., Schaerer et al., 2016; Tang
et al., 2019; Maseda et al., 2020; Naidu et al., 2022), such that at fixed
$f_{\rm esc}$, these galaxies will provide more LyC photons to the IGM. These
same properties scale with $f_{\rm esc}$ in the multivariate Cox models, and
our identified candidate strong LCEs at $z>6$ generally have O32 $>10$ and/or
EW([O iii]+H$\beta$) $>1500$ Å. High $f_{\rm esc}$ may therefore be coupled
with high $\xi_{\rm ion}$ (e.g., Schaerer et al., 2016; Naidu et al., 2022).

However, ionizing photon production may vary with other parameters as well,
such as galaxy luminosity (e.g., Bouwens et al., 2016; Finkelstein et al.,
2019). Fujimoto et al. (2023) estimate $\xi_{\rm ion}$ for spectroscopically
confirmed galaxies at $z\sim 8-9$ and find that $\xi_{\rm ion}$ is two times
higher at $M_{\rm UV}\sim-19.5$ compared to $M_{\rm UV}\sim-21.5$. Although
this enhanced efficiency is not enough to compensate for the six times fainter
UV luminosity at $M_{\rm UV}\sim-19.5$ vs. -21.5, these fainter galaxies are
also approximately 40 times more numerous (Bouwens et al., 2015). Atek et al.
(2023) infer similarly high $\xi_{\rm ion}$ for even fainter ($M_{\rm
UV}\sim-17$ to -15) and even more numerous galaxies, which only need $f_{\rm
esc}$$<0.05$ to drive reionization. For the faintest galaxies we consider
($M_{\rm 1500}>-17.5$), the median of the Cox model predicted $f_{\rm esc}$
values is near this threshold, with median $f_{\rm esc}$$\sim 0.04$. However,
these $f_{\rm esc}$ estimates require refinement with additional parameters
(e.g., O32, $\Sigma_{\rm SFR}$). If $f_{\rm esc}$ does not vary strongly with
luminosity, as suggested by our preliminary, limited models, the fainter
population would dominate reionization due to their higher $\xi_{\rm ion}$ and
greater numbers.

Assessing the main contributors to reionization requires progress on several
fronts. At $z>6$, we need estimates of nebular and morphological parameters in
addition to $\beta_{\rm 1550}$ and $M_{\rm 1500}$ for galaxy samples spanning
a wide range of luminosities. Parameters such as O32 and $\Sigma_{\rm SFR}$
have some of the greatest effects on the $f_{\rm esc}$ prediction accuracy for
both the LzLCS+ (Paper I) and galaxies at $z\sim 3$ (§3). The limited
magnitude range of the LzLCS+ reference sample ($M_{\rm 1500}\sim-18.5$ to
$\sim-21.5$) also introduces uncertainty. Our $f_{\rm esc}$ estimates for
spectroscopic samples within this magnitude range are likely reasonable and
suggest moderately high median $f_{\rm esc}$ $\sim 0.04-0.05$, with some
galaxies reaching $f_{\rm esc}$ as high as 0.6-0.7 (see Figure 4 and
predictions in Table 6 for the ELG-O32 model and Table 7 for the
ELG-O32-$\beta$ model). However, our $f_{\rm esc}$ estimates for fainter and
brighter galaxy populations rely on extrapolation and are less trustworthy. To
confirm these estimates, we need to explore $f_{\rm esc}$ and its dependence
on galaxy properties across a wider magnitude range at $z<6$. Lastly,
observations of the IGM and galaxy population at $z>6$ will provide a further
test of $f_{\rm esc}$ predictions. An accurate model of $f_{\rm esc}$ should
reproduce both the timing and topology of reionization with the $z>6$ galaxy
population. While much progress is needed to both improve and confirm $f_{\rm
esc}$ estimates, our current $f_{\rm esc}$ predictions imply that plausible
contributors to reionization appear at all magnitude ranges, and star-forming
galaxies with the required levels of LyC escape do indeed exist at $z>6$.

##  7 Summary

Quantifying the LyC escape fraction ($f_{\rm esc}$) of galaxies is critical to
understand the sources of cosmic reionization. We have developed a flexible
tool for predicting $f_{\rm esc}$ empirically using combinations of observable
variables available in the $z\sim 0.3$ LzLCS+ reference sample (Paper I). We
generate multivariate diagnostics for $f_{\rm esc}$ with the Cox proportional
hazards model (Cox, 1972), a survival analysis technique that appropriately
treats data with upper limits. We test Cox models developed from the $z\sim
0.3$ galaxies on $f_{\rm esc}$ observations at $z\sim 3$, and we apply the
models to several samples of $z\gtrsim 6$ galaxies to predict $f_{\rm esc}$
for galaxies in the epoch of reionization. The Appendix provides all the
information necessary to use the multivariate Cox models to predict $f_{\rm
esc}$ for other samples of galaxies, including samples at $z\sim 3$ and
samples at $z\gtrsim 6$, where the LyC is unobservable.333We have also
developed a python script, available at https://github.com/sflury/LyCsurv,
that allows the user to generate new Cox proportional hazards models to
predict $f_{\rm esc}$ using custom combinations of variables for any given
galaxy population. Version 0.1.0 of the code is archived in Zenodo (Flury et
al., 2024).

We summarize our main findings below.

  1. 1.

The models successfully reproduce the observed $f_{\rm esc}$ values for the
high-redshift $z\sim 3$ galaxies and often have a lower RMS for the $z\sim 3$
samples than for the $z\sim 0.3$ galaxies. The success of these models
suggests that low-redshift and high-redshift LCEs may share similar
properties. (§3)

  2. 2.

The best-performing models for the $z\sim 3$ galaxies include the dust
attenuation inferred from the UV SED (E(B-V)UV) or the UV slope $\beta_{\rm
1550}$, plus either O32 or a morphological measurement ($r_{\rm 50,NUV}$ or
$\Sigma_{\rm SFR}$). However, this result is tentative; each model is tested
on a different subset of high-redshift galaxies, which makes comparing the
models difficult. To determine which variable combinations best predict
$f_{\rm esc}$ at high redshift, we require larger $z\sim 3$ samples of LyC-
emitters with the full suite of input variables. (§3)

  3. 3.

We generate new Cox models based on the LzLCS+ observations, which incorporate
variables measured for $z\gtrsim 6$ samples. One model, using $M_{\rm 1500}$,
$M_{*}$, E(B-V)UV and EW([O iii]+H$\beta$) as input variables, applies to 180
galaxies from Endsley et al. (2021, 2023); Bouwens et al. (2023); Tang et al.
(2023); Fujimoto et al. (2023); Saxena et al. (2023b) and predicts a median
$f_{\rm esc}$ $=0.047$ for this combined sample. Of these 180 galaxies, 27%
have $f_{\rm esc}$ $>0.1$ and only 6% have $f_{\rm esc}$ $>0.2$. However, this
set of variables results in predictions that tend to underestimate $f_{\rm
esc}$, and the galaxy samples mostly consist of photometric measurements. A
second model, using $M_{\rm 1500}$, $M_{*}$, $\beta_{\rm 1550}$, and $r_{\rm
50,NUV}$ as input variables, applies to 278 galaxies, mostly with photometric
redshifts, from Morishita et al. (2023) and Mascia et al. (2023a). This model
predicts a higher median $f_{\rm esc}$ $=0.14$, with 56% and 39% of the
galaxies having $f_{\rm esc}$ $>0.1$ and $f_{\rm esc}$ $>0.2$, respectively.
Many of these galaxies are more compact than the LzLCS+ sample, and this model
therefore requires extrapolating beyond the LzLCS+ parameter space. (§4.1)

  4. 4.

Smaller samples of spectroscopically-confirmed $z\gtrsim 6$ galaxies have