# Multivariate Predictors of Lyman Continuum Escape. I. A Survival Analysis of the Low-redshift Lyman Continuum Survey

Anne E. Jaskot, Anneliese C. Silveyra, Anna Plantinga, Sophia R. Flury, Matthew Hayes, John Chisholm, Timothy Heckman, Laura Pentericci, Daniel Schaerer, Maxime Trebitsch, et al.

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

## ► To cite this version:

Anne E. Jaskot, Anneliese C. Silveyra, Anna Plantinga, Sophia R. Flury, Matthew Hayes, et al. Multivariate Predictors of Lyman Continuum Escape. I. A Survival Analysis of the Low-redshift Lyman Continuum Survey. The Astrophysical Journal, 2024, 972, 10.3847/1538-4357/ad58b9. insu-04838850

HAL Id: insu-04838850

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

Submitted on 15 Dec 2024

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

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

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

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

# Multivariate Predictors of Lyman Continuum Escape. I. A Survival Analysis of the Low-redshift Lyman Continuum Survey*

Anne E. Jaskot<fnref n="1" />, Anneliese C. Silveyra<fnref n="1" />,<fnref n="2" />, Anna Plantinga<fnref n="3" />, Sophia R. Flury<fnref n="4" />, Matthew Hayes<fnref n="5" />, John Chisholm<fnref n="6" />, Timothy Heckman<fnref n="7" />, Laura Pentericci<fnref n="8" />, Daniel Schaerer<fnref n="9" />, Maxime Trebitsch<fnref n="10" />, Anne Verhamme<fnref n="9" />,<fnref n="11" />, Cody Carr<fnref n="12" />,<fnref n="13" />, Henry C. Ferguson<fnref n="14" />, Zhiyuan Ji<fnref n="15" />, Mauro Giavalisco<fnref n="4" />, Alaina Henry<fnref n="14" />, Rui Marques-Chaves<fnref n="9" />, Göran Östlin<fnref n="5" />, Alberto Saldana-Lopez<fnref n="5" />, Claudia Scarlata<fnref n="16" />, Gábor Worseck<fnref n="17" />, and Xinfeng Xu<fnref n="18" />

## Abstract

To understand how galaxies reionized the Universe, we must determine how the escape fraction of Lyman continuum (LyC) photons ($f_{esc}$) depends on galaxy properties. Using the $z \sim 0.3$ Low-redshift Lyman Continuum Survey (LzLCS), we develop and analyze new multivariate predictors of $f_{esc}$. These predictions use the Cox proportional hazards model, a survival analysis technique that incorporates both detections and upper limits. Our best model predicts the LzLCS $f_{esc}$ detections with an rms scatter of 0.31 dex, better than single-variable correlations. According to ranking techniques, the most important predictors of $f_{esc}$ are the equivalent width (EW) of Lyman-series absorption lines and the UV dust attenuation, which track line-of-sight absorption due to H I and dust. The H I absorption EW is uniquely crucial for predicting $f_{esc}$ for the strongest LyC emitters, which show properties similar to weaker LyC emitters and whose high $f_{esc}$ may therefore result from favorable orientation. In the absence of H I information, star formation rate surface density ($\Sigma_{SFR}$) and [O III]/[O II] ratio are the most predictive variables and highlight the connection between feedback and $f_{esc}$. We generate a model suitable for $z > 6$, which uses only the UV slope, $\Sigma_{SFR}$, and [O III]/[O II]. We find that $\Sigma_{SFR}$ is more important in predicting $f_{esc}$ at higher stellar masses, whereas [O III]/[O II] plays a greater role at lower masses. We also analyze predictions for other parameters, such as the ionizing-to-nonionizing flux ratio and Ly$\alpha$ escape fraction. These multivariate models represent a promising tool for predicting $f_{esc}$ at high redshift.

Unified Astronomy Thesaurus concepts: Astrostatistics (1882); Reionization (1383); High-redshift galaxies (734); Starburst galaxies (1570); Interstellar medium (847); Ultraviolet astronomy (1736); Radiative transfer (1335)

## 1. Introduction

Star-forming galaxies likely caused one of the most significant transformations in cosmic history: the reionization of hydrogen in the intergalactic medium (IGM) at $z \gtrsim 6$.

Current constraints on the galaxy and quasar luminosity functions at $z > 6$ and on galaxies’ ionizing photon production efficiencies favor stars as the dominant source of ionizing, Lyman continuum (LyC) photons (e.g., Bouwens et al. 2015, 2016; Finkelstein et al. 2015, 2019; Robertson et al. 2015; Ricci et al. 2017; De Barros et al. 2019; Faucher-Giguère 2020; Shen et al. 2020; Endsley et al. 2021). Exactly which star-forming galaxies contribute these LyC photons remains unclear, however. Studies disagree as to whether low-, intermediate-, or high-luminosity galaxies dominate reionization (e.g., Razoumov & Sommer-Larsen 2010; Wyithe & Loeb 2013; Wise et al. 2014; Paardekooper et al. 2015; Finkelstein et al. 2019; Ma et al. 2020; Naidu et al. 2020; Cain et al. 2021; Izotov et al. 2021; Begley et al. 2022; Matthee et al. 2022; Rosdahl et al. 2022; Saldana-Lopez et al. 2023) or whether

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

<footnote n="1">¹Department of Astronomy, Williams College, Williamstown, MA 01267, USA</footnote>

<footnote n="2">²Department of Physics, University of Nevada, Reno, NV 89557, USA</footnote>

<footnote n="3">³Department of Mathematics & Statistics, Williams College, Williamstown, MA 01267, USA</footnote>

<footnote n="4">⁴Department of Astronomy, University of Massachusetts Amherst, Amherst, MA 01002, USA</footnote>

<footnote n="5">⁵Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="6">⁶Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote n="7">⁷Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="8">⁸INAF—Osservatorio Astronomico di Roma, via Frascati 33, 00078, Monteporzio Catone, Italy</footnote>

<footnote n="9">⁹Observatoire de Genève, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="10">¹⁰Astronomy, Kapteyn Astronomical Institute, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="11">¹¹Univ. Lyon, Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="12">¹²Center for Cosmology and Computational Astrophysics, Institute for Advanced Study in Physics, Zhejiang University, Hangzhou 310058, People's Republic of China</footnote>

<footnote n="13">¹³Institute of Astronomy, School of Physics, Zhejiang University, Hangzhou 310058, People's Republic of China</footnote>

<footnote n="14">¹⁴Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="15">¹⁵Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="16">¹⁶Minnesota Institute for Astrophysics, School of Physics and Astronomy, University of Minnesota, 316 Church St. SE, Minneapolis, MN 55455, USA</footnote>

<footnote n="17">¹⁷VDI/VDE Innovation+Technik, Berlin, Germany</footnote>

<footnote n="18">¹⁸Center for Interdisciplinary Exploration and Research in Astrophysics, Northwestern University, Evanston, IL 60201, USA</footnote>

<footnote>Received 2024 January 11; revised 2024 May 30; accepted 2024 June 14; published 2024 August 27</footnote>

<footnote>* Based 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.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
other properties such as concentrated star formation, nebular ionization, or starburst age demarcate the LyC-emitting galaxy population (e.g., Heckman et al. 2001; Clarke & Oey 2002; Jaskot & Oey 2013; Zastrow et al. 2013; Nakajima & Ouchi 2014; Alexandroff et al. 2015; Ma et al. 2015; Sharma et al. 2016; Trebitsch et al. 2017; Izotov et al. 2018b; Marchi et al. 2018; Naidu et al. 2022).

Observations with the James Webb Space Telescope (JWST) are revealing that proposed LyC-emitting galaxy populations exist at high redshift (e.g., Schaerer et al. 2022b; Endsley et al. 2023; Fujimoto et al. 2023; Tang et al. 2023; Williams et al. 2023; Atek et al. 2024; Mascia et al. 2024), although it remains to be seen whether they are sufficiently numerous at $z > 6$. However, identifying the galaxies responsible for reionization also requires knowing $f_{esc}$, the fraction of LyC photons that escape into the IGM. The value of $f_{esc}$ for $z > 6$ galaxies is not known, and its physical connection with galaxy properties may be complex. LyC escape can depend on a galaxy's interstellar gas geometry, dust content, stellar feedback, and gravitational potential and may vary with time (e.g., Heckman et al. 2001, 2011; Wise et al. 2014; Ma et al. 2015; Sharma et al. 2016; Trebitsch et al. 2017; Chisholm et al. 2018; Barrow et al. 2020; Gazagnes et al. 2020; Mauerhofer et al. 2021; Saldana-Lopez et al. 2022). LyC escape is an inherently multiparameter problem.

Because of attenuation in the IGM, detecting the LyC flux from galaxies becomes unlikely above $z \sim 4$ (e.g., Inoue et al. 2014). As a result, the astronomy community has relied on lower-redshift samples in order to investigate $f_{esc}$ and its dependence on galaxy properties. Studies of $z \sim 2-4$ galaxies find that $f_{esc}$ may be enhanced in galaxies with lower UV luminosities, lower dust attenuation, higher $\text{Ly}\alpha$ equivalent widths (EWs), strong [O III] $\lambda5007$ emission, and/or compact sizes (e.g., Marchi et al. 2018; Steidel et al. 2018; Bassett et al. 2019; Fletcher et al. 2019; Nakajima et al. 2020; Begley et al. 2022; Saxena et al. 2022). More local samples, at $z \sim 0.3$, likewise identify strong $\text{Ly}\alpha$ emission, elevated [O III]/[O II] ratios, high star formation rate (SFR) surface densities ($\Sigma_{SFR}$), and low dust attenuation as characteristics of LyC emitters (LCEs; e.g., Borthakur et al. 2014; Izotov et al. 2016b, 2018b; Verhamme et al. 2017; Chisholm et al. 2018). Nevertheless, $f_{esc}$ can have a range of values, even for galaxies with similar properties (e.g., Vanzella et al. 2016; Rutkowski et al. 2017; Izotov et al. 2018b; Fletcher et al. 2019; Bian & Fan 2020; Marques-Chaves et al. 2021).

The recently completed Low-redshift Lyman Continuum Survey (LzLCS. Flury et al. 2022b) offers a new opportunity to investigate all these physical properties simultaneously in a large statistical sample of LCEs and nonemitters. Through the LzLCS and archival programs (Izotov et al. 2016a, 2016b, 2018a, 2018b, 2021; Wang et al. 2019), 89 galaxies at $z \sim 0.3$ have LyC observations from the Hubble Space Telescope (HST) Cosmic Origins Spectrograph (COS), and 50 of these galaxies are detected in the LyC. This combined sample, hereafter the LzLCS+, covers a wide range of luminosities, metallicities, $H\beta$ and $\text{Ly}\alpha$ EWs, [O III]/[O II], and $\Sigma_{SFR}$ (Flury et al. 2022b), which enables it to more clearly reveal the trends and scatter between $f_{esc}$ and galaxy properties. By quantifying the relationship between $f_{esc}$ and observables, the LzLCS+ can generate predictions for $f_{esc}$ in $z > 6$ galaxies, where direct LyC detections are inaccessible.

The results from the LzLCS+ confirm that $f_{esc}$ correlates with a variety of galaxy properties, from line-of-sight measurements such as H I covering fraction, dust attenuation, and $\text{Ly}\alpha$ escape fraction to global properties such as [O III]/[O II] and $\Sigma_{SFR}$ (e.g., Flury et al. 2022a; Chisholm et al. 2022; Saldana-Lopez et al. 2022; Xu et al. 2023). The latter properties hint at the possible role of mechanical or radiative feedback in creating low optical depth sight lines that allow ionizing photons to escape. Nevertheless, all correlations between $f_{esc}$ and observable properties show significant scatter (Wang et al. 2021; Flury et al. 2022a; Chisholm et al. 2022; Saldana-Lopez et al. 2022; Xu et al. 2023).

A combination of properties might more accurately predict $f_{esc}$ than a single variable alone. For instance, several studies have generated multivariate predictions of $\text{Ly}\alpha$ properties for low- and high-redshift galaxies (e.g., Yang et al. 2017; Trainor et al. 2019; Runnholm et al. 2020). Runnholm et al. (2020) present such an analysis for $\text{Ly}\alpha$ luminosity using galaxies in the $\text{Ly}\alpha$ Reference Sample (LARS; Hayes et al. 2013; Östlin et al. 2014). By applying multivariate regression using physical and observable galaxy properties, they predict the $\text{Ly}\alpha$ luminosities of the LARS galaxies with an rms residual scatter of 0.2–0.3 dex.

Similarly, Maji et al. (2022) and Choustikov et al. (2024) perform multivariate linear fits for $f_{esc}$ using simulated galaxies from the SPHINX cosmological simulations (Rosdahl et al. 2018). Maji et al. (2022) find that a combination of four variables (escaping $\text{Ly}\alpha$ luminosity, gas mass, gas metallicity, and recent SFR) can explain 85% of the variance in escaping LyC luminosity, while three variables (escaping $\text{Ly}\alpha$ luminosity, recent SFR, and gas mass) explain 66% of the variance in $f_{esc}$. However, some of the relevant variables (e.g., $\text{Ly}\alpha$ luminosity and gas mass) will be difficult or impossible to obtain for many galaxies at $z > 6$. In contrast, Choustikov et al. (2024) focus on multivariate $f_{esc}$ predictions using observable properties, such as [O III]/[O II], UV magnitude, and UV slope. While promising, the Maji et al. (2022) and Choustikov et al. (2024) simulation results require further testing and confirmation using observational samples.

Other recent studies have taken a more empirical approach by using the $z \sim 0.3$ LyC observations to derive relationships between $f_{esc}$ and observable properties. Lin et al. (2024) fit for the probability of a galaxy having detectable versus undetected $f_{esc}$ as a function of UV magnitude, UV slope, and [O III]/[O II] and use these results to assess the likelihood of LyC escape from galaxies at $z > 6$. Mascia et al. (2023) develop a model that predicts $f_{esc}$ based on some of the strongest correlating variables in the LzLCS: $\beta_{1550}$, UV half-light radius, and [O III]/[O II].

In this paper, we investigate a variety of empirical multivariate predictions for $f_{esc}$ using the LzLCS+, the largest observational sample of LyC measurements at low redshift, which enables both reliable individual $f_{esc}$ measurements and comprehensive rest-frame UV and optical ancillary data. Whereas the $\text{Ly}\alpha$ data used by Runnholm et al. (2020) include only $\text{Ly}\alpha$ emitters and the simulated LyC data used by Maji et al. (2022) and Choustikov et al. (2024) include measurements down to $f_{esc} = 0$, the LzLCS+ sample includes 39 galaxies with nondetected LyC corresponding to $f_{esc}$ upper limits of 0.03%–5.9%. Standard multivariate regression models do not account for upper limits, but a valid statistical analysis of the full LzLCS+ data set should include information from
both detections and nondetections (e.g., Isobe et al. 1986). Hence, we turn to the statistical approach of survival analysis, which can handle censored data like those of the LzLCS. We apply the Cox proportional hazards model (Cox 1972), a semiparametric survival analysis method, to the full LzLCS+ data set to derive multivariate predictions of $f_{esc}$. By testing different sets of variables and assessing their predictive ability, we explore which physical properties combine to set a galaxy's $f_{esc}$. We adopt a cosmology of $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$.

## 2. Methods

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

To derive multivariate predictors of $f_{esc}$, we use the LzLCS (Flury et al. 2022b), the largest sample of LCEs at low redshift. Flury et al. (2022b) describe the survey sample, observations, data processing, and measurements in full, but we review some of the key details here. The LzLCS is a 134-orbit Cycle 26 HST program (GO-15626; PI: A. E. Jaskot) that obtained far-UV (FUV) spectra, including the rest-frame LyC, with the COS G140L grating for 66 star-forming galaxies at $z \sim 0.3$. The targeted galaxies each fulfill one or more hypothesized selection criteria for LCEs: high nebular ionization (O32 = [O III] $\lambda5007/[O II]$ $\lambda3727 \ge 3$), blue UV slopes (power-law index $\beta < -2$), and/or concentrated star formation ($\Sigma_{SFR} > 0.1 \text{ M}_\odot \text{ yr}^{-1} \text{ kpc}^{-2}$). We combine the LzLCS with archival COS observations from Izotov et al. (2016a, 2016b, 2018a, 2018b, 2021) and Wang et al. (2019). We exclude one galaxy, J1333+6246 (Izotov et al. 2016b), whose [O III] $\lambda\lambda5007, 4959, H\alpha$, and $H\beta$ line fluxes may be inaccurate; these lines appear truncated in the Sloan Digital Sky Survey (SDSS) spectrum, and the galaxy's $H\alpha/H\beta$ and $H\beta/H\gamma$ line ratios are unphysical. Excluding this problematic galaxy, the LzLCS+ sample consists of 88 low-redshift galaxies with LyC observations. The sample spans a range of $10^8-10^{10} \text{ M}_\odot$ in stellar mass, $12 + \log_{10}(\text{O/H}) = 7.5-8.5$, and $-21.5$ to $-18.3$ in observed (not corrected for internal reddening) 1500 Å absolute magnitude.

We process the raw COS spectra using the CALCOS pipeline (v3.3.9) and the FAINTCOS software routines (Worseck et al. 2016; Makan et al. 2021) to reduce and calibrate the spectra and account for the backgrounds due to dark current and scattered geocoronal $Ly\alpha$. We correct all spectra for Milky Way attenuation using the Green et al. (2018) dust maps and Fitzpatrick (1999) attenuation law. The LyC flux measurements represent the flux in a 20 Å wide region near rest-frame 900 Å; we exclude wavelengths above an observed wavelength of 1180 Å because of telluric contamination. For consistency, we reprocess and remeasure the archival observations using this same methodology. Following the criteria in Flury et al. (2022b), we define LyC detections as observations that have a probability $<0.02275$ of originating from background counts. By this definition, 49 of the 88 total galaxies have detected LyC. As in Flury et al. (2022b), we adopt the 84th percentile of the background count distribution as the upper limit for nondetections.

Using these LyC fluxes, we then calculate the absolute $f_{esc}$ from different estimates of the intrinsic LyC, as described in Flury et al. (2022b); for this paper, we adopt the $f_{esc}$ values from Saldana-Lopez et al. (2022), derived from STARBURST99 (Leitherer et al. 2011, 2014) spectral energy distribution (SED)

fits to the UV continuum, following the methods outlined in Chisholm et al. (2019). As discussed in Flury et al. (2022b), the UV SED-fitting method is more reliable than $H\beta$-derived $f_{esc}$ estimates for the diverse LzLCS sample because the presence of underlying, older (>10 Myr old) stellar populations that emit optical emission can bias the $H\beta f_{esc}$ results. The $H\beta$ estimates of intrinsic LyC also do not account for LyC photons absorbed by dust and assume isotropic escape. With our preferred UV-based $f_{esc}$ estimate, the LzLCS+ LyC detections cover $f_{esc} = 0.4\%–88.9\%$, and the nondetections have upper limits ranging from 0.03% to 5.9%. We also consider an alternative, purely empirical measure of LyC escape, the $F_{\lambda LyC}/F_{\lambda 1100}$ flux ratio in Section 4.1.

In addition to LyC measurements, the LzLCS+ data set contains a wealth of ancillary measurements, described in full in Flury et al. (2022b) and Saldana-Lopez et al. (2022). Each galaxy in the survey has rest-frame optical photometry and spectroscopy through SDSS (Blanton et al. 2017) and UV photometry via the Galaxy Evolution Explorer (GALEX. Martin et al. 2003). We derive stellar masses ($M_*$) using PROSPECTOR (Leja et al. 2017; Johnson et al. 2019) fits to the SDSS and GALEX photometry (Flury et al. 2022b). We measure nebular line fluxes by fitting the SDSS spectra with multiple Gaussian profiles and then iteratively derive the electron temperature ($T_e$), electron density, stellar absorption, and nebular $E(B-V)$ from the [O III] $\lambda\lambda5007, 4959, 4363$, [S II] $\lambda\lambda6716, 6731$, and Balmer lines. Where the [S II] doublet or [O III] $\lambda4363$ is undetected, we adopt $n_e = 100 \text{ cm}^{-3}$ and the estimated [O III] $\lambda4363$ flux from the Pilyugin et al. (2006) “ff- relation,” respectively. We then derive oxygen abundances via the direct method as implemented in PYNEB (Luridiana et al. 2015). To estimate the SFR, we use the dust-corrected $H\beta$ luminosities, case B $H\alpha/H\beta$ ratio (Storey & Hummer 1995), and Kennicutt & Evans (2012) SFR calibration.

The COS data provide estimates of additional physical parameters. We use the COS near-UV acquisition images to calculate the UV half-light radius. The sources all appear to be very compact within the central $1''$ diameter of the COS aperture, so that vignetting will not strongly affect our radius measurements. From the measured radius, we also derive $\Sigma_{SFR}$ as

$$ \Sigma_{SFR} = \frac{\text{SFR}}{2\pi r_{50,\text{NUV}}^2}. \qquad (1) $$

From the STARBURST99 SED fits to the UV spectra, we fit for the dust excess $E(B-V)$. We denote the $E(B-V)$ derived from the UV spectral fits as $E(B-V)_{UV}$ to distinguish it from $E(B-V)_{\text{neb}}$ derived from the Balmer line ratios. We also extrapolate the UV SED fits to infer the “observed” (non- extinction-corrected) absolute magnitude at 1500 Å ($M_{1500}$) and the power-law index slope at 1550 Å ($\beta_{1550}$; Saldana-Lopez et al. 2022).

The G140L spectra also cover $Ly\alpha$. For our $Ly\alpha$ measurements, we use spectra extracted using a slightly wider spatial aperture (30 pixels vs. 25 pixels) because the scattered $Ly\alpha$ emission may be more spatially extended than the continuum light (Flury et al. 2022b). After masking the Si III $\lambda1206$ and N V $\lambda1240$ regions, we linearly fit the continuum within 100 Å of $Ly\alpha$ using an iterative sigma clipping algorithm and integrate the $Ly\alpha$ fluxes relative to this continuum level. This linear continuum fit typically agrees
within <5% with the continuum estimated from SED fits to the 25-pixel-aperture extracted spectra (Saldana-Lopez et al. 2022). We do not correct for the small contribution of stellar Lyα absorption, but we conservatively adopt a 25% uncertainty on the Lyα continuum estimate (Flury et al. 2022b). From these fluxes, we derive Lyα luminosities, EWs, and escape fractions ($f_{esc, Lyα}$); we use the dust-corrected Hβ flux and the galaxies' measured electron temperatures and densities to infer the intrinsic Lyα flux. The reported Lyα measurements represent the sum of both absorption and emission along the line of sight. Even galaxies with net Lyα emission may have underlying Lyα absorption troughs, which are commonly detected in samples at $z < 0.1$ (e.g., McKinney et al. 2019; Hu et al. 2023). At the higher redshifts of our sample, $z \sim 0.3$, the COS aperture covers a larger physical area, and these troughs may experience more infilling, due to scattered Lyα emission. Because the radii of our Lyα spectral apertures are $\sim 2.6$ times larger than the UV continuum half-light radius (Flury et al. 2022b), we expect our Lyα measurements to capture most of the scattered emission (e.g., Hayes et al. 2013); future Lyα imaging of LzLCS+ galaxies will test this assumption (HST-GO-17069; PI: M. Hayes). Nine galaxies in our sample have significant Lyα absorption that overlaps with the Si III $\lambda$1206 absorption feature. To account for the uncertain strength of Si III $\lambda$1206 contamination in these galaxies, we remeasure the Lyα fluxes and EWs, excluding wavelengths within $500 \text{ km s}^{-1}$ of Si III $\lambda$1206, and we increase our Lyα uncertainties to account for this possible flux difference. We note that this effect is minor, changing EWs by $<3 \AA$ and $f_{esc, Lyα}$ by $\le 0.017$, with a median change of only 0.001 in $f_{esc, Lyα}$ for the nine affected galaxies.

Lastly, we use the G140L spectra to measure the EWs and residual fluxes near line center ($R_l$) for a variety of low-ionization absorption lines, as described in Saldana-Lopez et al. (2022). The UV absorption-line measurements use the same (25-pixel-aperture) spectral extraction as the LyC measurements, and the SED fits are used to estimate the continuum level. For optically thick gas and well-resolved absorption lines, and assuming a uniform dust screen, the residual flux is related to the gas covering fraction ($C_f$) by

$$R(\lambda_i) = 1 - C_f(\lambda_i) \qquad (2)$$

(see, e.g., Gazagnes et al. 2018, for a discussion of geometry effects on $R$). This equation assumes that all lines are saturated; while this assumption appears to hold for most of the H I lines, the observational uncertainties are too high to draw conclusions about saturation for the metal lines (see Saldana-Lopez et al. 2022 for more information). In this work, we consider some of the strongest low-ionization metal absorption lines (Si II $\lambda$1260, C II $\lambda$1334, H Ly$\beta$ + O I). To improve the signal-to-noise ratio (S/N) of the measurements, we also consider the inverse-variance weighted average of the EWs and $R_l$ values for a given galaxy, calculated for the available H I Lyman series lines (Ly$\beta$–Ly6) or for the available low-ionization state (LIS) metal lines (Si II $\lambda\lambda$1190, 1193, Si II $\lambda$1260, [O I] $\lambda$1302, C II $\lambda$1334; see Saldana-Lopez et al. 2022 for details). We refer to these average absorption-line measurements as EW(H I, abs), $R_l(H I, abs)$, EW(LIS), and $R_l(LIS)$. Because the measured residual intensity depends on resolution, we caution that the measured $R_l$ will not perfectly correspond to the true $R_l$ of the system; for the LzLCS+, this systematic error is of order 10%–20% and is

within the reported uncertainty (Saldana-Lopez et al. 2022). While trends between $R_l$ and $f_{esc}$ will still be apparent (e.g., Saldana-Lopez et al. 2022), future investigations should use care when applying models derived from the LzLCS+ $R_l$ data to spectroscopic data with a different resolution.

## 2.2. Survival Analysis

As described above, one of the strengths of the LzLCS+ data set is its large and diverse sample, spanning a wide range of measured $f_{esc}$. The full sample of 88 galaxies includes 39 nondetections with $f_{esc}$ upper limits of only a few percent. These nondetections can convey essential information about which properties do and do not distinguish LCEs from non-LCEs, and omitting nondetections can bias fitted relations (see Isobe et al. 1986). Standard multivariate linear regression does not account for censored data, that is, data with upper limits. Therefore, we instead apply survival analysis techniques to the LzLCS+ sample, as these techniques properly treat censored data.

As implied in the name, survival analysis originated in the field of medicine (for reviews, see Clark et al. 2003; Bradburn et al. 2003; for applications to astronomy, see Feigelson & Nelson 1985; Isobe et al. 1986). In a medical context, the censored data often consist of known and unknown lifetimes for individuals participating in a medical study. People who are alive at the end of the study have an unknown lifetime; all we know is that they will live longer than their current age. Survival analysis techniques ultimately seek to describe the probability of a particular lifetime for a population or to compare how the survival probability changes between different populations, based on data that include some known lifetimes and some limits.

To set up an analogous scenario for predicting the probability of a particular $f_{esc}$, we do the following. Instead of a “death” or known lifetime, we have a detection. Most survival analyses involve lower limits (“right-censored” data). For simplicity in following common techniques, we therefore adopt the absorbed fraction of LyC ($f_{abs} = 1 - f_{esc}$) and its corresponding lower limits rather than $f_{esc}$ and its upper limits. With this setup, increasing $f_{abs}$ (decreasing $f_{esc}$) is equivalent to increasing the “age” of the study participant. A nondetection, i.e., a lower limit to $f_{abs}$ (and upper limit on $f_{esc}$), is analogous to a lifetime greater than some threshold. Instead of starting at time = 0 and proceeding until we record a death or the study ends, we can imagine trying to observe a galaxy with $f_{esc} = 1$ ($f_{abs} = 0$) and proceeding with deeper and deeper observations until we have a detection or cease our efforts. Just as the threshold between recording a measured lifetime and a limit will depend on factors such as the length of the study and age of the participant, the threshold between detected $f_{esc}$ and an upper limit will depend on aspects of our study such as the observation depth and galaxy brightness.

### 2.2.1. The Cox Proportional Hazards Model

The Cox proportional hazards regression model (Cox 1972) describes the probability of an event in an infinitesimally small window of time, given a combination of independent variables and assuming that the event did not already occur. In the Cox model, this event probability (or “hazard function”) at time t for
<table><caption>**Table 1**
Complete List of Independent Variables</caption><thead><tr><th>Category</th><th>Variables</th></tr></thead><tbody><tr><td>Mass and luminosity</td><td>log<sub>10</sub>(M<sub>*</sub>); log<sub>10</sub>(SFR)-Hβ; M<sub>1500</sub></td></tr><tr><td>Nebular properties</td><td>log<sub>10</sub>(EW(Hβ)); E(B-V)<sub>neb</sub>; 12+log<sub>10</sub>(O/H)</td></tr><tr><td>Nebular ionization</td><td>log<sub>10</sub>([Ne III] λ3869/[O II] λ3727); log<sub>10</sub>(O32 = [O III] λ5007/[O II] λ3727); log<sub>10</sub>([O II] λ3727/Hβ)</td></tr><tr><td>Morphology</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1 missing]; log<sub>10</sub>(r<sub>50,NUV</sub>) [1 missing]</td></tr><tr><td>UV color and attenuation</td><td>β<sub>1200</sub>; β<sub>1550</sub>; E(B-V)<sub>UV</sub></td></tr><tr><td>Lyα</td><td>EW(Lyα); L(Lyα); f<sub>esc,Lyα</sub></td></tr><tr><td>UV absorption lines</td><td>EW(Lyβ) [9 missing]; EW(Si II λ1260) [2 missing]; EW(C II λ1334) [11 missing]; EW(H I,abs); EW(LIS)<br>R<sub>f</sub>(Lyβ) [3 missing]; R<sub>f</sub>(Si II λ1260); R<sub>f</sub>(C II λ1334) [4 missing]; R<sub>f</sub>(H I,abs) [5 missing]; R<sub>f</sub>(LIS) [2 missing]</td></tr></tbody></table>

**Note.** Numbers within square brackets denote the number of galaxies within the LzLCS+ sample missing these measurements. Positive values of EW(Hβ) and EW(Lyα) denote net emission; positive values of the UV absorption lines denote net absorption.

a set of variables *x* is modeled as

$$h(t|x) = h_0(t) \exp \left[ \sum_{i=1}^{n} b_i (x_i - \bar{x}_i) \right], \quad (3)$$

where $b_i$ are fitted coefficients for each variable $x_i$ and where $h_0(t)$ represents the baseline hazard, the probability of an event given average values of all variables ($\bar{x}_i$). An increase in variable $x_i$ relative to that variable's average value in the sample results in an exponential increase or decrease in the detection probability, depending on the sign of the coefficient $b_i$. In our case, instead of modeling the probability of an event at time *t*, the hazard function represents the probability of an LyC detection in an infinitesimally small window of $f_{esc}$ for a set of independent variables, assuming that the galaxy was not already detected at a higher value of $f_{esc}$.

The Cox model is semiparametric. One of its strengths is that it does not assume a particular functional form for the baseline hazard function but rather estimates it nonparametrically. In other words, the event times or detection values do not have to obey a known statistical distribution (Bradburn et al. 2003). Because the Cox model predicts the detection probability at each possible $f_{esc}$ value, rather than predicting the value of $f_{esc}$ itself, the model results are identical for $f_{esc}$ and for $log_{10}(f_{esc})$. In other words, the nonparametric estimate of $h_0(f_{abs})$ for an input array of $f_{esc}$ values would be identical to $h_0(log_{10}(f_{abs}))$ for the equivalent input array of $log_{10}(f_{esc})$ values.

However, the parametric nature of the hazard function equation (Equation (3)) with its reported fitted coefficients makes it straightforward to analyze the relative importance of variables and generate predictions for future data sets. On the other hand, this fixed functional form implicitly assumes that the effect of the independent variables is multiplicative and does not depend on *t*, in our case, $f_{abs}$. This functional form would not be appropriate, for example, if one variable tended to introduce a step function in $f_{esc}$, while the others had more of a smooth correlation. However, the adopted functional form has some flexibility in that independent variables are allowed to take any functional form (e.g., $x_i$ can be replaced with $log_{10}(x_i)$, $x_i^2$, etc). Another limitation of the Cox model is that while it handles censored dependent variable values, it cannot handle censored data in the independent variables. Consequently, we have to limit our analysis to input variables, $x_i$, that are available for all (or nearly all) of the sample galaxies, as we explain below. In summary, the Cox proportional hazards model serves as a multivariate regression model for describing

and predicting dependent variables with upper limits. Many other common statistical methods do not treat upper limits (e.g., multivariate linear regression, principal component analysis), are univariate analyses (e.g., Kaplan-Meier analysis), or are fully parametric, less flexible models (e.g., the Weibull model).

In fitting the Cox model to the LzLCS+, we experiment with different combinations of independent variables from the parent list given in Table 1. Where possible, we take the base-10 logarithm of the relevant variable. We cannot take the logarithm for variables that contain both positive and negative values, such as the UV LIS and Lyα lines; these lines range from absorption to emission within the sample. Putting variables on a logarithmic scale serves two purposes. First, it puts variables on a consistent scale of similar order of magnitude. In addition, the coefficients in Equation (3) have a simple interpretation; increasing a variable by an order of magnitude leads to a factor of $e_i^b$ increase in the probability $h(f_{abs}|x)$. Exploring alternative functional forms for the input variables is outside the scope of this paper, but we test the effects of switching between logarithmic and linear scalings in the following subsection (Section 2.2.2). While some statistical methods require putting all variables on an identical scale, such as scaling from 0 to 1, this further rescaling is not necessary in the Cox model. Equation (3) effectively rescales variables; the sample mean is subtracted from each input variable measurement, and the $b_i$ parameter translates a linear change in variable $x_i$ to a change in the probability of measuring a given $f_{esc}$. Any normalization would change the derived values of $b_i$ but not the statistical significance of input variables or their ultimate effect on $f_{esc}$.

In selecting a subset of variables to use in our fits, we include multiple variables from each category in Table 1, but we avoid any variables that are closely related to each other and that may therefore be highly collinear. Collinearity can result in multiple best-fit solutions and can cause the Cox model to fail to converge. Consequently, we typically consider either $\Sigma_{SFR}$ alone or SFR and $r_{50,NUV}$, but not all three variables at once. Similarly, we choose only one variable that traces nebular ionization, one variable to trace UV color, and limited subsets of the UV absorption-line variables. We discuss the effect of selecting different variables in Section 3.2. As mentioned above, the Cox proportional hazard model cannot account for missing or censored independent variable data. Hence, we exclude a galaxy from our fits if it is missing any of the required variables. Numbers in square brackets in Table 1 indicate the number of galaxies missing a given measurement.
Most of the models discussed in this paper exclude only one galaxy, J1046+5827, a nonleaker that lacks a reported $\Sigma_{SFR}$ and $r_{50,NUV}$. We apply the Cox model to $f_{esc}$ as our primary dependent variable, but we also perform and discuss model fits to the $F_{\lambda,LC}/F_{\lambda,1100}$ ratio (Section 4.1), LyC luminosity (Section 4.1), $f_{esc,Lya}$, and Ly$\alpha$ luminosity (Section 4.2).

We apply the Cox proportional hazards model to our data using the Python package lifelines (Davidson-Pilon 2019). The lifelines CoxPHFitter routine returns the best-fit coefficients for our selected independent variable set, their p-values, and various goodness-of-fit statistics. As shown by Equation (3), the Cox proportional hazards model gives the probability of observing a particular $f_{esc}$ value given a set of physical or observable properties. Instead, our goal is to predict the expected value of $f_{esc}$ given a set of properties. As described below, to find the expected $f_{esc}$, we find the median of the probability distribution, where the true $f_{esc}$ has an equal probability of lying above or below this adopted value (e.g., Bradburn et al. 2003; Davidson-Pilon 2019); $f_{esc}$ is predicted to be above the median 50% of the time and below it the other 50% of the time.

The median $f_{esc} = 1 - f_{abs}$ represents the $f_{esc}$ value where the survival function $S(f_{abs})$, the probability that there is no detection at $f_{abs,detect} < f_{abs}$, reaches 0.5. This probability is equivalent to the probability that $f_{esc,detect}$ is detected at $f_{esc, detect} < f_{esc}$. The survival function is calculated as (e.g., Cox 1972; Bradburn et al. 2003; Davidson-Pilon 2019; McLernon et al. 2023)

$$S(f_{abs}) = \exp[-\text{HF}_0(f_{abs}) \cdot \text{ph}(x)], \quad (4)$$

where $\text{HF}_0$ is the baseline cumulative hazard function,

$$\text{HF}_0(f_{\text{abs}}) = \int_0^{f_{\text{abs}}} h_0(f) df, \quad (5)$$

and $\text{ph}(x)$ is the partial hazard function for a set of variables $x$,

$$\text{ph}(x) = \exp\left[\sum_{i=1}^{n} b_i (x_i - \bar{x}_i)\right]. \quad (6)$$

Occasionally, $S(f_{abs})$ for a set of parameters will never reach 0.5, and the predicted median $f_{esc}$ is indeterminate. This situation corresponds to an arbitrarily small predicted $f_{esc} \sim 0$.

The best-fit coefficients for the Cox proportional hazards model are found by maximizing the partial likelihood, which compares $\text{ph}(x)$ for each detection with the sum of $\text{ph}(x)$ for all galaxies with higher $f_{abs}$ (lower $f_{esc}$), including both detections and nondetections. The resulting coefficients $b_i$ are those that best sort $f_{abs}$ in order. Once the coefficients are determined, Breslow's estimator (Breslow 1972) determines the baseline cumulative hazard function, $\text{HF}_0(f_{abs})$, using the number of detections with values lower than each $f_{abs}$ (i.e., detections with higher values of $f_{esc}$) plus the $\text{ph}(x)$ values for all detections and limits higher than $f_{abs}$ (i.e., with lower $f_{esc}$). The model fits $\text{HF}_0(f_{abs})$ nonparametrically, reporting a value of $\text{HF}_0$ for each of the input $f_{abs}$ values of the LzLCS+ sample.

One can also use the survival function to evaluate the values at which $S$ reaches 0.159 and 0.841, the bounds corresponding to the normal-theory $1\sigma$ uncertainty of the predicted $f_{esc}$. These probabilities account for the scatter in the correlations between $f_{esc}$ and the independent variables. The scatter may arise from observational uncertainty in the measurements, as well as inherent variation among the population. In addition to this method, we have also used a Monte Carlo (MC) method to

explore how the $f_{esc}$ predictions change if we vary each variable within its uncertainties. We randomly resample each independent variable measurement and each dependent variable detection using their uncertainties, rerun the Cox fit, and obtain new estimates of the median predicted $f_{esc}$ for each galaxy. Using the distribution of predicted $f_{esc}$, we then calculate the 15.9th and 84.1th percentiles. For nearly all galaxies, the uncertainties estimated using the survival function are greater than or equal to the uncertainties determined from this MC method, which indicates that the scatter in the correlations is the dominant effect. We conclude that the survival function sufficiently represents the uncertainty in predicted $f_{esc}$ in most cases.

### 2.2.2. Interpreting Cox Model $f_{esc}$ Predictions

As shown by Equation (4), the Cox model does not represent an equation that shows how $f_{esc}$ itself depends on particular variables. Rather, it describes how the probability of observing a given $f_{esc}$ changes for different sets of galaxy properties. The baseline cumulative hazard function $\text{HF}_0(f_{abs})$ describes the expected probability distribution for galaxies that have the average properties of the LzLCS+ data set. In this case, $x_i = \bar{x}_i$ for all variables, and the reference probability that the observed $f_{esc,detect}$ is less than a particular value of $f_{esc}$ is simply $\text{S}_{\text{ref}}(f_{esc}) = \exp[-\text{HF}_0(f_{esc})]$. If we increase one variable by an increment of 1, $x_i - \bar{x}_i = 1$, $\text{ph}(x) = \exp(b_i)$, and the new probability $S(f_{esc}) = \text{S}_{\text{ref}}(f_{esc})^{\exp(b_i)}$. Because the probability $\text{S}_{\text{ref}}(f_{esc}) < 1$, raising it to the power of $\exp(b_i)$ decreases the probability, which means that low values of $f_{esc,detect} < f_{esc}$ are less likely. The probability distribution therefore shifts such that the probability corresponding to $S = 0.5$ occurs at higher $f_{esc}$, resulting in a larger predicted median $f_{esc}$. Changing a second variable, $x_j$, by an increment of 1 changes the new probability again by a power of $\exp(b_j)$, such that $S(f_{esc}) = [\text{S}_{\text{ref}}(f_{esc})^{\exp(b_i)}]^{\exp(b_j)}$.

We illustrate an example of the survival function probabilities, the associated median $f_{esc}$, and the dependence of these quantities on input variables in Figure 1. In the left panel, we show a simple Cox model with three input variables: $\beta_{1550}$, $\log_{10}(\Sigma_{SFR})$, and $\log_{10}(\text{O32})$. Circles represent the predicted median $f_{esc}$ where $S(f_{esc}) = 0.5$ for each value of O32. For $\log_{10}(\text{O32}) = 1$, $S(f_{esc}) = 0.5$ at $f_{esc} = 0.005$. Changing only the $\log_{10}(\text{O32})$ value by 1, from O32 = 1, the blue dotted line, to O32 = 10, the magenta long-dashed line, shifts $S(f_{esc} = 0.005)$ from 0.5 down to $0.5^{\exp(0.996)} = 0.153$, where 0.996 is the best-fit coefficient $b_i$ for the $\log_{10}(\text{O32})$ variable. Thus, according to this model, high $f_{esc}$ values are more common among galaxies with high O32. In this example, we changed only one parameter, but changing the other variables at the same time could either further shift the probability distribution to higher $f_{esc}$ or counteract the change caused by increasing O32. We note that, unlike an equation for $f_{esc}$ as a function of these input variables, which could conceivably reach unphysical values of $f_{esc} > 1$ for extreme parameters, the probability distribution can shift toward higher values of $f_{esc}$, but it never extends beyond $f_{esc} = 1$.

Figure 1 also shows how the predicted probability responds to linear versus logarithmic variables. As outlined above, the probabilities shift by a power of $\exp(b_i)$ for a step size of 1 in $x_i$, which corresponds to an order-of-magnitude increase if $x_i$ represents the base-10 logarithm of a variable. In the left panel the model uses $\log_{10}(\text{O32})$ as the input variable, and in the right
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Examples of the survival function probabilities predicted from the Cox model. The left panel uses $\beta_{1550}$, $\log_{10}(\Sigma_{SFR})$, and $\log_{10}(\text{O32})$ as input variables, and the right panel uses $\beta_{1550}$, $\log_{10}(\Sigma_{SFR})$, and O32 instead of $\log_{10}(\text{O32})$. Lines show the probability that the detected $f_{\text{esc,detect}}$ is less than a particular value of $f_{\text{esc}}$ or, equivalently, the probability that $f_{\text{abs,detect}}$ is not less than a given $f_{\text{abs}}$. The lines are plotted at $f_{\text{esc}}$ values corresponding to the measured $f_{\text{esc}}$ for the LzLCS+ data set, from which the baseline cumulative hazard function is estimated nonparametrically. Data points show the median predicted $f_{\text{esc}}$, where $S(f_{\text{esc}}) = 0.5$. The different lines show predictions for different input values of O32. As O32 increases, the probability distribution shifts to higher values of $f_{\text{esc}}$, and the extent of the shift depends on whether the variable is linear or logarithmic. The black solid line in the left panel ($\log_{10}(\text{O32}) = -1$) shows an example where $S$ never reaches 0.5, indicating an inferred median $f_{\text{esc}} \sim 0$.</figcaption></figure>

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

panel we show the effect of changing the form of the input variable from $\log_{10}(\text{O32})$ to linear O32. The differences are most apparent where a large linear spacing (e.g., changing O32 from 10 to 25) is not equivalent to a large logarithmic spacing and vice versa (e.g., O32 from 0.1 to 1). However, the change in the predicted median $f_{\text{esc}}$ values is largely minor. In changing to a linear scale, $f_{\text{esc}}$ changes from 0 to 0.006 for O32 = 0.1 and from 0.027 to 0.059 for O32 = 25. For the other plotted values, the change in $f_{\text{esc}}$ is <0.002. The best-fit coefficients $b_i$ and $\text{HF}_0(f_{\text{esc}})$ values also differ between the two models and are optimized to match the observed distribution of $f_{\text{esc}}$ in the LzLCS+. Consequently, although the choice of variable form does affect predictions for galaxies at the extremes of the input parameter space, our main results, including the overall quality of our model predictions, are not highly sensitive to the functional form of the input variables.

We do not thoroughly explore functional forms here, but we experiment with changing the interstellar medium (ISM) absorption measurements (EW(H I,abs) and $R_I(\text{H I,abs})$) and Ly$\alpha$ variables (EW(Ly$\alpha$), $L(\text{Ly}\alpha)$, $f_{\text{esc,Ly}\alpha}$) to logarithmic forms. For the logarithmic Ly$\alpha$ variables, we exclude all galaxies with net Ly$\alpha$ absorption from the data set. We find that the same variables are statistically significant in our models, although their exact coefficients and *p*-values necessarily change. The goodness-of-fit metrics (Section 2.3) are also similar, typically changing by $\le 0.02$. The change in variable form also does not affect our main conclusions regarding the most important ranked variables (Section 2.4), although it can result in minor changes in ranked variable order.

Observationally, for most variables, $\log_{10}(f_{\text{esc}})$ does seem to change significantly with the logarithm of a variable (Flury et al. 2022a). These trends could indicate that $f_{\text{esc}}$ has a power-law dependence on many variables, although, as noted above, such relationships could only apply over the range where the predicted $f_{\text{esc}} \le 1$. One exception to this logarithmic

dependence is the UV absorption-line measurements, which do show a dependence of $\log_{10}(f_{\text{esc}})$ on the linear form of the variables (Saldana-Lopez et al. 2022). This functional form may indicate a more complex dependence of $f_{\text{esc}}$ on the gas geometry. We find that our predictions using the linear EW(H I, abs) better match the observed $f_{\text{esc}}$ by every metric compared to models using $\log_{10}(\text{EW(H I,abs)})$, and we therefore choose to keep EW(H I,abs) in the linear form. Models using the linear or logarithmic form of $R_I(\text{H I,abs})$ perform comparably well, but with the linear form better matching the observed $f_{\text{esc}}$ for the strongest LCEs. Given their comparable or improved performance and for consistency with the LIS line measurements, which can reach negative values, we therefore keep all absorption-line variables in their linear form. We also choose to keep the Ly$\alpha$ measurements in a linear form so as to not bias our data set by excluding non-Ly$\alpha$ emitters. We emphasize that our fit quality and main results are not sensitive to these choices regarding variable form.

In summary, the Cox model predicts the probability of observing particular $f_{\text{esc}}$ values for different galaxy populations and can identify which variables most affect that probability. Statistically significant variables signify that the $f_{\text{esc}}$ probability is highly responsive to an incremental change in that variable: either a linear or logarithmic increase, depending on the variable form. While this work identifies these significant variables, future theoretical or observational programs could endeavor to derive the exact functional dependence of $f_{\text{esc}}$ on these parameters.

## 2.3. Model Assessment

We assess the goodness of fit for our models in a few different ways. Our primary method is the concordance index (Harrell et al. 1982), which applies to censored data. The concordance index considers all possible pairs of data points and how the observed rank order of $f_{\text{esc}}$ compares to the order predicted by the model. The model classifies each pair as
concordant if the data point with a higher observed $f_{esc}$ also has a higher predicted $f_{esc}$ and discordant if it does not. Pairs can also be tied, if their predicted $f_{esc}$ values are identical. For some pairs with upper limits, their rank is ambiguous. The concordance index is

$$C = \frac{n_c + 0.5n_t}{n_c + n_d + n_t}, \qquad (7)$$

where $n_c$ is the number of concordant pairs, $n_t$ is the number of tied pairs, and $n_d$ is the number of discordant pairs. The concordance index ranges from 0 (perfect disagreement) to 0.5 (perfectly random) to 1.0 (perfect concordance; e.g., Davidson-Pilon 2019).<fnref n="19" /> We also calculate the $R^2$ statistic for the galaxies with $f_{esc}$ detections as

$$R^2 = 1 - \frac{\sum_i (y_i - f_i)^2}{\sum_i (y_i - \bar{y})^2}, \qquad (8)$$

where $y_i$ are the observed values of $\log_{10}(f_{esc})$, $\bar{y}$ is their mean value, and $f_i$ are the model-predicted $\log_{10}(f_{esc})$ values. Following Maji et al. (2022), we also calculate the adjusted $R^2$:

$$R_{\text{adj}}^2 = 1 - (1 - R^2) \frac{n-1}{n-p-1}, \qquad (9)$$

where $n$ is the number of data points and $p$ is the number of independent variables in the model. The $R_{\text{adj}}^2$ parameter measures whether additional variables improve the model fit more than would be expected by random chance. Lastly, we also calculate the rms dispersion about the model predictions:

$$\text{rms} = \sqrt{\frac{\sum_i (y_i - f_i)^2}{n}}. \qquad (10)$$

We can only calculate these three quantities, $R^2$, $R_{\text{adj}}^2$, and rms, for the detections in our sample. Although they cannot provide a complete picture of model performance, they can tell us how well our predictions work for the galaxies with known $f_{esc}$ and can quantify our model performance at the high-$f_{esc}$ end. Each of these metrics correlates strongly with the others and with the concordance index, and our main results are insensitive to the exact metric used.

We choose to evaluate both $R^2$ and rms using $\log_{10}(f_{esc})$ rather than the linear $f_{esc}$, because the scatter in the predicted $f_{esc}$ values is approximately constant in logarithmic space over the full observed $f_{esc}$ range. A single rms value is therefore representative of both low- and high-$f_{esc}$ galaxies. In contrast, in linear space, the rms changes systematically from low to high observed $f_{esc}$. For instance, for our fiducial model (Section 3.1), the rms is 0.02 for galaxies with $f_{esc}$ in the lowest third of the sample, 0.16 for the next third, and 0.22 for the final third, whereas the logarithmic rms is ~0.3–0.4 dex across the full sample range. Because $R^2$ and rms use $\log_{10}(f_{esc})$, these statistics only include galaxies with both detected $f_{esc}$ and nonzero predicted $f_{esc}$. In contrast, the concordance index $C$ incorporates the full data set, including

galaxies with observed upper limits and those with predicted $f_{esc} = 0$.

We also assess the predictive ability of our models using cross-validation, applying the model to data that were not used in generating the model itself. We perform a k-fold cross-validation, in which we randomly split the data set into k groups. We combine k-1 of the groups to become our training set, which we use to fit the model. We then test the model on the remaining group using the $C$, $R^2$, and rms metrics described above. ($R_{\text{adj}}^2$ is often undefined, due to runs where $n$ for the small test set is equal to $p$). We then repeat this process until each group has been used as the test group with the remaining groups used as the training set. Following Runnholm et al. (2020), we adopt $k=3$, repeat the group selection 100 times, and average the final results. A value of $k=3$ gives us sufficient galaxies in the test set to calculate statistics, and the repetition of the process ensures that our results are not sensitive to the exact group selected (Runnholm et al. 2020). We provide the results of these tests in our evaluation of the models and distinguish the metrics derived from this cross-validation analysis by the subscript “CV.”

## 2.4. Variable Selection

We explore different subsets of independent variables and discuss the resulting model fits in Section 3.2. With these different selections, we test the effect of switching the absorption line probed, of including or excluding $\text{Ly}\alpha$, and of limiting the variable set to observables accessible to JWST (Section 3.4). However, we also wish to analyze the relative importance of each independent variable and identify the variable subsets that most accurately predict $f_{esc}$. To evaluate the importance of each variable, we turn to the tools of forward and backward selection (e.g., Runnholm et al. 2020; Maji et al. 2022).

In forward selection, we run the Cox model with each independent variable individually in turn and determine which variable provides a fit with the highest concordance. We select this variable and then combine it with each remaining variable in turn to determine which combination of two variables provides the highest concordance. We proceed in this manner until we have a ranked order of the most significant variables. Conversely, with backward selection, we start with the full set of independent variables and remove each one in turn. The combination of variables that gives the highest concordance identifies the least important excluded variable. We exclude it and proceed with the remaining variables, identifying and removing the next least important variable each time to obtain a rank-ordered list. Importantly, each round of forward or backward selection only ranks variables with respect to the current best-performing model. Hence, some variables may be poorly ranked if they provide similar information to an already-selected variable. Variable rankings can reveal the most important predictors of $f_{esc}$, but a poor ranking does not necessarily imply that a variable is uncorrelated with $f_{esc}$.

When we perform our forward and backward selections, we do not use the full list of variables in Table 1 but limit it in the following ways. First, we include only the EW(H I,abs), EW(LIS), $R_l(\text{H I,abs})$, and $R_l(\text{LIS})$ measurements, which are averages of multiple lines, rather than the measurements of individual absorption lines. This choice avoids overly limiting our sample sizes, since 22 galaxies are missing one or more individual absorption-line measurements. Second, we avoid

<footnote n="19">¹⁹ For the concordance index calculation in the lifelines package, pairs are ranked by their partial hazards function $\text{ph}(x)$, rather than by median $f_{esc}$ or $f_{abs}$. The $\text{ph}(x)$ will scale with $f_{esc}$, except for the cases where $S(f_{esc})$ is indeterminate and median $f_{esc} \sim 0$. Two galaxies that both have $f_{esc} \sim 0$ could still have slightly different $\text{ph}(x)$, corresponding to different probability distributions for $f_{esc} > 0$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 2.** The $f_{esc}$ predictions from the fiducial Cox model compared with the observed $f_{esc}$ values for the LzLCS+ sample. Red circles indicate detections, and blue triangles represent the 1σ upper limit on nondetections. The dashed line shows a one-to-one correspondence. The black cross in the upper left corner indicates the median size of the uncertainties in the observed and predicted $f_{esc}$. The uncertainties in the predicted $f_{esc}$ represent the 15.9th and 84.1th percentiles of the $f_{esc}$ probability distribution. We plot galaxies with predicted $f_{esc} \sim 0$ at the bottom of the figure with their predicted $f_{esc}$ set to $10^{-3}$. The bottom panel shows a plot of the residuals on a linear scale. This model predicts $f_{esc}$ for the LyC detections with an rms scatter of 0.36 dex.</figcaption></figure>

<table><caption>**Table 2**
Fiducial Model Coefficients</caption><thead><tr><th>Variable</th><th>Coefficient <i>b</i></th><th><i>p</i>-value</th></tr></thead><tbody><tr><td><i>f</i><sub>esc,Lyα</sub></td><td>7.00</td><td>6.0E-9</td></tr><tr><td><i>E</i>(<i>B</i> − <i>V</i>)<sub>UV</sub></td><td>−14.90</td><td>6.3E-5</td></tr><tr><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>1.83</td><td>1.1E-3</td></tr><tr><td><i>E</i>(<i>B</i> − <i>V</i>)<sub>neb</sub></td><td>6.82</td><td>2.8E-3</td></tr><tr><td>log<sub>10</sub>(O32)</td><td>2.65</td><td>1.6E-2</td></tr><tr><td>log<sub>10</sub>(<i>M</i><sub>*</sub>)</td><td>0.63</td><td>0.18</td></tr><tr><td>12+log<sub>10</sub>(O/H)</td><td>1.46</td><td>0.22</td></tr><tr><td><i>M</i><sub>1500</sub></td><td>−0.46</td><td>0.21</td></tr><tr><td>EW(LIS)</td><td>−0.26</td><td>0.32</td></tr><tr><td>log<sub>10</sub>(EW(Hβ))</td><td>−0.04</td><td>0.98</td></tr></tbody></table>

**Note.** Variables are listed in order of their *p*-value significance. Positive values of EW(LIS) represent net absorption; positive values of EW(Hβ) represent net emission.

highly collinear variables and variables that measure very similar properties. We include only $\Sigma_{SFR}$ rather than SFR and $r_{50,NUV}$, we use only one ionization-sensitive line ratio (O32 = [O III]/[O II]), and we use only one measure of UV dust attenuation ($\beta_{1550}$ or $E(B-V)_{UV}$). When included separately, alternative measures of these properties end up with the same or nearly the same ranks, which suggests that they are indeed largely interchangeable. We explore three different sets of variables in our rankings: one set includes $\beta_{1550}$ and all remaining variables, one set corresponds to the variables in our best-performing Cox model, and one set is limited to variables accessible at high redshift. We also test our final variable rankings using an MC method. We rerun the

forward and backward selection processes after randomly sampling the independent and dependent variables 100 times using their uncertainties. We present the results of the forward and backward selection methods in Section 3.3.

### 3. Predicting $f_{esc}$

#### 3.1. Fiducial Model

For our fiducial model, we choose the following variables: $\log_{10}(M_*)$, $M_{1500}$, $\log_{10}(\text{EW}(H\beta))$, $E(B-V)_{\text{neb}}$, $12+\log_{10}(O/H)$, $\log_{10}(O32)$, $\log_{10}(\Sigma_{SFR})$, $E(B-V)_{\text{UV}}$, $f_{\text{esc,Lyα}}$, and EW(LIS). We present the fitted coefficients and significance for these variables in Table 2 and compare the model-predicted median $f_{esc}$ with the observations in Figure 2. We list the goodness-of-fit metrics for this model in Table 3. Statistically significant variables, with *p*-values $\le 0.05$, are $f_{\text{esc,Lyα}}$, $E(B-V)_{\text{UV}}$, $\log_{10}(\Sigma_{SFR})$, $E(B-V)_{\text{neb}}$, and $\log_{10}(O32)$. Not surprisingly, these same variables individually correlate or anticorrelate with $f_{\text{esc}}$, as discussed in Flury et al. (2022a) and Saldana-Lopez et al. (2022). We choose to include both O32 and EW(Hβ) in the fiducial model, even though they correlate with each other (e.g., Flury et al. 2022a), because they may have subtly different relationships with $f_{\text{esc}}$; a young starburst may have both high EW(Hβ) and high O32, but high global $f_{\text{esc}}$ will tend to increase O32, while decreasing EW(Hβ) (e.g., Zackrisson et al. 2013; Nakajima & Ouchi 2014; Jaskot & Ravindranath 2016). However, we find that EW(Hβ) is not significant and its coefficient is near zero, which indicates that EW(Hβ) does not contribute meaningful information to the fit beyond what O32 already provides. For the fiducial model, we choose EW(LIS) instead of EW(Hβ), because EW(LIS) is less affected by IGM and circumgalactic medium opacity and is
<table><caption>**Table 3**
Goodness of Fit for Cox $f_{esc}$ Models</caption><thead><tr><th>Model</th><th>Sample Size</th><th>R<sup>2</sup></th><th>R<sup>2</sup><sub>adj</sub></th><th>rms</th><th>C</th><th>R<sup>2</sup><sub>CV</sub></th><th>rms<sub>CV</sub></th><th>C<sub>CV</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>87</td><td>0.60</td><td>0.49</td><td>0.36</td><td>0.89</td><td>0.55</td><td>0.36</td><td>0.85</td></tr><tr><td colspan="9">Fiducial Model with Modifications<sup>a</sup></td></tr><tr><td>EW(H I,abs) instead of EW(LIS)</td><td>87</td><td>0.69</td><td>0.60</td><td>0.31</td><td>0.91</td><td>0.64</td><td>0.31</td><td>0.88</td></tr><tr><td>With both EW(Lyα) and f<sub>esc,Lyα</sub></td><td>87</td><td>0.61</td><td>0.48</td><td>0.36</td><td>0.89</td><td>0.56</td><td>0.35</td><td>0.86</td></tr><tr><td>EW(Lyβ) instead of EW(LIS)</td><td>78<sup>b</sup></td><td>0.64</td><td>0.52</td><td>0.34</td><td>0.89</td><td>0.58</td><td>0.34</td><td>0.85</td></tr><tr><td>EW(C II) instead of EW(LIS)</td><td>77<sup>b</sup></td><td>0.59</td><td>0.46</td><td>0.37</td><td>0.89</td><td>0.54</td><td>0.36</td><td>0.85</td></tr><tr><td>EW(Si II) instead of EW(LIS)</td><td>85<sup>b</sup></td><td>0.63</td><td>0.52</td><td>0.33</td><td>0.89</td><td>0.56</td><td>0.33</td><td>0.85</td></tr><tr><td>[O II]/Hβ instead of O32</td><td>87</td><td>0.62</td><td>0.51</td><td>0.35</td><td>0.88</td><td>0.56</td><td>0.34</td><td>0.85</td></tr><tr><td>R<sub>f(H I,abs)</sub> instead of EW(LIS)</td><td>82<sup>b</sup></td><td>0.61</td><td>0.50</td><td>0.34</td><td>0.88</td><td>0.56</td><td>0.33</td><td>0.85</td></tr><tr><td>β<sub>1550</sub> instead of E(B-V)<sub>UV</sub></td><td>87</td><td>0.58</td><td>0.46</td><td>0.37</td><td>0.88</td><td>0.53</td><td>0.37</td><td>0.85</td></tr><tr><td>R<sub>f(Si II)</sub> instead of EW(LIS)</td><td>87</td><td>0.58</td><td>0.47</td><td>0.37</td><td>0.88</td><td>0.53</td><td>0.37</td><td>0.85</td></tr><tr><td>R<sub>f(C II)</sub> instead of EW(LIS)</td><td>84<sup>b</sup></td><td>0.55</td><td>0.41</td><td>0.38</td><td>0.88</td><td>0.49</td><td>0.37</td><td>0.85</td></tr><tr><td>[Ne III]/[O II] instead of O32</td><td>87</td><td>0.57</td><td>0.44</td><td>0.38</td><td>0.88</td><td>0.51</td><td>0.37</td><td>0.85</td></tr><tr><td>R<sub>f(Lyβ)</sub> instead of EW(LIS)</td><td>84<sup>b</sup></td><td>0.52</td><td>0.37</td><td>0.39</td><td>0.88</td><td>0.45</td><td>0.39</td><td>0.84</td></tr><tr><td>R<sub>f(LIS)</sub> instead of EW(LIS)</td><td>85<sup>b</sup></td><td>0.60</td><td>0.49</td><td>0.36</td><td>0.88</td><td>0.54</td><td>0.36</td><td>0.84</td></tr><tr><td>Without O32</td><td>87</td><td>0.57</td><td>0.46</td><td>0.38</td><td>0.88</td><td>0.51</td><td>0.37</td><td>0.85</td></tr><tr><td>Without E(B-V)<sub>neb</sub></td><td>87</td><td>0.54</td><td>0.42</td><td>0.38</td><td>0.88</td><td>0.48</td><td>0.38</td><td>0.84</td></tr><tr><td>EW(Lyα) instead of f<sub>esc,Lyα</sub></td><td>87</td><td>0.54</td><td>0.41</td><td>0.39</td><td>0.88</td><td>0.49</td><td>0.39</td><td>0.83</td></tr><tr><td>β<sub>1200</sub> instead of E(B-V)<sub>UV</sub></td><td>87</td><td>0.51</td><td>0.36</td><td>0.40</td><td>0.88</td><td>0.43</td><td>0.39</td><td>0.85</td></tr><tr><td>Without Σ<sub>SFR</sub></td><td>88<sup>b</sup></td><td>0.49</td><td>0.37</td><td>0.41</td><td>0.87</td><td>0.41</td><td>0.40</td><td>0.84</td></tr><tr><td>Without E(B-V)<sub>UV</sub></td><td>87</td><td>0.30</td><td>0.13</td><td>0.48</td><td>0.86</td><td>0.23</td><td>0.47</td><td>0.83</td></tr><tr><td>L(Lyα) instead of f<sub>esc,Lyα</sub></td><td>87</td><td>0.41</td><td>0.24</td><td>0.44</td><td>0.85</td><td>0.34</td><td>0.44</td><td>0.82</td></tr><tr><td>Without f<sub>esc,Lyα</sub></td><td>87</td><td>0.35</td><td>0.19</td><td>0.46</td><td>0.84</td><td>0.27</td><td>0.46</td><td>0.79</td></tr><tr><td colspan="9">JWST Models</td></tr><tr><td>Full JWST model</td><td>87</td><td>0.29</td><td>0.14</td><td>0.47</td><td>0.83</td><td>0.21</td><td>0.47</td><td>0.79</td></tr><tr><td>β<sub>1550</sub>+log<sub>10</sub>(Σ<sub>SFR</sub>+log<sub>10</sub>(O32))</td><td>87</td><td>0.34</td><td>0.29</td><td>0.46</td><td>0.83</td><td>0.27</td><td>0.45</td><td>0.81</td></tr><tr><td>β<sub>1550</sub>+log<sub>10</sub>(Σ<sub>SFR</sub>+log<sub>10</sub>([Ne III]/[O II]))</td><td>87</td><td>0.40</td><td>0.35</td><td>0.44</td><td>0.83</td><td>0.32</td><td>0.43</td><td>0.82</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Models are listed in order of decreasing concordance, *C*. We substitute different measures of dust, UV absorption lines, Lyα, and nebular ionization as indicated. We also experiment with removing each of the statistically significant variables in the fiducial model.

<sup>b</sup> These models apply to fewer or more galaxies than the fiducial model, and some differences in goodness-of-fit metrics relative to the fiducial model may result from the changed sample.

therefore a viable indirect indicator across a wider redshift range. We consider the effect of substituting EW(H I) or other measurements of absorption-line strength in Section 3.2.

Although most variables show the behavior expected from their individual relationship with $f_{esc}$, one exception is the nebular attenuation, $E(B-V)_{neb}$, which correlates with $f_{esc}$ in the Cox model instead of anticorrelating as might be expected. However, the model already contains the UV $E(B-V)$ and $f_{esc,Lyα}$, two parameters that show strong trends with $E(B-V)_{neb}$. The nebular $E(B-V)$ therefore essentially operates as a second-order effect, where the nebular attenuation would increase $f_{esc}$ at a fixed value of line-of-sight UV attenuation or $f_{esc,Lyα}$. This extra nebular $E(B-V)$ correlation could indicate a trend with some other physical property not included in the model, such as starburst age, global gas content, or gas clumpiness. In the case of $f_{esc,Lyα}$, $E(B-V)_{neb}$ may help quantify how much of the Lyα escape is due to a low H I optical depth. At fixed $f_{esc,Lyα}$, a higher $E(B-V)_{neb}$ implies a greater contribution of dust and weaker contribution of H I to the Lyα absorption. As expected, if we remove both Lyα and E (B-V)UV from the model, the correlation with $E(B-V)_{neb}$ disappears.

As discussed in Section 2.2.2, the best-fit coefficients in Table 2 quantify how the probability of observing a given $f_{esc}$ responds to an incremental change in each variable. The

probability *S* of observing $f_{esc,detect} < f_{esc}$ changes from the original probability $S_{ref}$ by $S = S_{ref}^{exp(b_i \Delta x_i)}$ for a change in variable $x_i$. The coefficients in Table 2 indicate that changing $f_{esc,Lyα}$ by 0.1 results in raising the probability to the power of $exp(0.7) = 2.0$, leading to a lower probability of low $f_{esc}$ values (see Section 2.2.2). Conversely, increasing $E(B-V)_{UV}$ by 0.1 increases the probability of low $f_{esc}$ values by an even greater factor; the probability is raised to $exp(-1.49) \sim 1/5$. For the other statistically significant variables, a 0.1 change in $E(B-V)_{neb}$ also raises the probability to the power of $\sim 2$, while 0.1 changes in $log_{10}(\Sigma_{SFR})$ and $log_{10}(O32)$ raise the probability to powers of 1.2 and 1.3, respectively.

As illustrated in Figure 2 and Table 3, the fiducial model predicts the observed $f_{esc}$ values for detections with an rms scatter of 0.36 dex for both the full sample and when the data set is split into test and training sets ($rms_{CV} = 0.36$; see Section 2.3). The model’s high concordance, $C = 0.89$ ($C_{CV} = 0.85$), shows that it successfully ranks galaxies on the basis of $f_{esc}$. The $R^2$, $R^2_{adj}$, and $R^2_{CV}$ metrics are 0.60, 0.49, and 0.55, respectively. This successful fit shows that including multiple variables results in a substantial improvement over predictions using a single variable. Even for some of the best single-variable correlations, $f_{esc}$ spans 2 dex at a given galaxy property (e.g., Flury et al. 2022a; Chisholm et al. 2022).
<figure><img src="image_8.png" /><figcaption>**Figure 3.** The $f_{\text{esc}}$ predictions from the Cox model improve when we include EW(H I,abs) instead of EW(LIS). Symbols are the same as in Figure 2.</figcaption></figure>

The improved accuracy of the fiducial multivariate model is driven by only a handful of variables. Table 2 shows that only four variables are statistically significant in the fiducial model, and limiting the model to only these four variables achieves comparable accuracy, with the $R^2$, rms, and C metrics changing by 0.01–0.02. As can be seen in Table 2, highly insignificant variables tend to have best-fit coefficients near 0, which means that the variable plays almost no role in the prediction. Excluding these variables therefore has little effect on the model, and because these extra parameters have negligible effects, overfitting is not a major concern.

Despite some overall improvement, the scatter in the fiducial Cox model is still significant, which demonstrates the difficulty in predicting the line-of-sight $f_{\text{esc}}$ even given a large amount of information. The residuals in the bottom panel of Figure 2 also show that the fiducial model systematically underpredicts the value of $f_{\text{esc}}$ for the strongest LCEs, even though it does typically identify them as having $f_{\text{esc}} > 0.1$. Mascia et al. (2023) find a similar result and suggest that it may arise from the limited number of strong LCEs in the LzLCS+. We suggest that this underprediction arises from a different effect, where the strongest LCEs represent a subset of galaxies with favorably oriented low column density channels but whose global $f_{\text{esc}}$ may actually be substantially lower. We elaborate on this possibility in Sections 3.2 and 3.6. Because LyC may escape from narrow channels, the chance orientation of a galaxy can introduce scatter in relationships between global galaxy properties and the measured line-of-sight $f_{\text{esc}}$ (e.g., Cen & Kimm 2015; Flury et al. 2022a; Seive et al. 2022). Although the fiducial model explicitly includes variables that trace line-of-sight properties, such as EW(LIS) and $E(B-V)_{\text{UV}}$, the LIS metals may imperfectly trace the LyC-absorbing H I gas.

Other effects may also contribute to the observed scatter. Simulations show that $f_{\text{esc}}$ may fluctuate in time, which may introduce scatter between $f_{\text{esc}}$ and observables that are sensitive to a different timescale (e.g., Trebitsch et al. 2017;

Barrow et al. 2020). Other possible causes of the model scatter and underprediction of $f_{\text{esc}}$ are uncertainties in the observed $f_{\text{esc}}$ values or a disconnect between the global properties we measure and the local properties of the LyC-emitting region (e.g., Martin et al. 2015; Kim et al. 2023). Rivera-Thorsen et al. (2022) point out that a small number of massive stars can potentially dominate the escaping LyC emission; the properties of this LyC source region may differ from both the global galaxy properties and the properties inferred from the nonionizing UV spectrum. Finally, because different properties may regulate LyC escape in different types of galaxies (e.g., Flury et al. 2022a), one relation may not suffice to predict $f_{\text{esc}}$ in all types of galaxies. We explore this possibility further in Section 3.5.

## 3.2. Modifications to the Fiducial Model

Modifications to the fiducial model demonstrate that properties sensitive to dust and H I column density are the most useful predictors of $f_{\text{esc}}$. We summarize the performance of the fiducial model and compare it with different modifications in Table 3. We do not test all possible combinations of the variables in Table 1, but we experiment with modifying the fiducial model by swapping alternative measures of some properties and by dropping each statistically significant variable.

We first substitute different absorption-line tracers (Table 1) for EW(LIS). Of these different tracers, only EW(Ly$\beta$), EW(H I,abs), and $R_I(\text{H I,abs})$ give statistically significant coefficients with $p < 0.05$. Because the sample size changes slightly for the other low-ionization metal absorption-line selections (see Table 1), we cannot compare the metrics for the other models in detail. However, the fit quality appears comparable to the fiducial model ($\text{rms} = 0.33–0.39$ and $C = 0.88–0.89$). Only EW(H I,abs) clearly improves the model fit, with rms = 0.31 and C = 0.91.
The fit with EW(H I,abs) (Figure 3) is also our best-
performing model, and it scores better than the fiducial model
on all metrics (Table 3). Moreover, this model more accurately
predicts *f*<sub>esc</sub> for strong LCEs (Figure 3). Hence, the *f*<sub>esc</sub>
underprediction in the fiducial model arises, at least in part,
from the fact that EW(LIS) and the other variables imperfectly
trace the line-of-sight H I gas. The Cox model is predicting the
median *f*<sub>esc</sub> for a given set of parameters, and the strongest
LCEs may not be noticeably distinct from more moderate LCEs
in the fiducial model's set of parameters. In contrast, the low
EW(H I,abs) values of the strongest leakers do distinguish them
from almost all of the moderate LCEs (e.g., Saldana-Lopez
et al. 2022). The exact conversion from EW(LIS) to H I
absorption and *f*<sub>esc</sub> depends on metallicity and gas geometry,
since the trace metals will only show up in absorption for a
sufficiently large H I column. While a range of physical
conditions could cause weak LIS absorption, weak H I
absorption more directly indicates LyC escape. Lastly, the
observational uncertainty in EW(LIS) is also higher than in
EW(H I,abs), and these uncertainties will be the most
problematic for the strongest LCEs with the weakest absorption
lines. In conclusion, while LIS lines can identify LCEs, H I
absorption lines appear necessary to accurately determine their
line-of-sight *f*<sub>esc</sub>.

Like EW(H I,abs), Lyα measurements can contribute key
information about H I column density. Indeed, Maji et al.
(2022) find that Lyα luminosity, *L*(Lyα), is the most important
variable in predicting the LyC *f*<sub>esc</sub> from simulations. In our
fiducial model, *f*<sub>esc,Lyα</sub> was the most statistically significant
coefficient, consistent with observational and theoretical
studies that find a connection between LyC and Lyα (e.g.,
Dijkstra et al. 2016; Verhamme et al. 2017; Steidel et al. 2018;
Izotov et al. 2020). In the fiducial model, *f*<sub>esc,Lyα</sub> serves as the
best estimate of the H I optical depth. When EW(H I,abs) is
substituted in the model, *f*<sub>esc,Lyα</sub> remains statistically sig-
nificant, but its *p*-value drops from 6E-9 to 3E-3 and its best-fit
coefficient declines from 7 to 4. Because it is less sensitive to
scattered emission, EW(H I,abs) more directly traces the line-
of-sight H I optical depth, which reduces the model's reliance
on *f*<sub>esc,Lyα</sub>. Without EW(H I,abs), however, *f*<sub>esc,Lyα</sub> provides
crucial information in the fiducial model. If we exclude it, the
fit quality worsens by all metrics; the rms rises from 0.36 to
0.46, *C* drops from 0.89 to 0.84, and *R*² drops from 0.60 to
0.35. Substituting *L*(Lyα) for *f*<sub>esc,Lyα</sub> also decreases the quality
of the fit (rms = 0.44, *C* = 0.85, *R*² = 0.41), whereas EW(Lyα)
performs almost as well as *f*<sub>esc,Lyα</sub> (rms = 0.39, *C* = 0.88,
*R*² = 0.54). Interestingly, when included together, both
EW(Lyα) and *f*<sub>esc,Lyα</sub> show up as significant coefficients. The
two variables may provide slightly different information
relevant to LyC escape. While *f*<sub>esc,Lyα</sub> may be more directly
linked to the fraction of escaping LyC, the EW provides
additional information about starburst properties such as the
intrinsic Lyα production and underlying continuum. However,
including both EW(Lyα) and *f*<sub>esc,Lyα</sub> results in marginal, if any,
improvement in the fit. All metrics change by only 0–0.01.

Of all the variables in the fiducial model, excluding *E* (*B* − *V*)<sub>UV</sub> or *f*<sub>esc,Lyα</sub> has the most detrimental effect on the model fit. Without information on the line-of-sight dust attenuation from *E*(*B* − *V*)<sub>UV</sub>, *R*<sup>2</sup> drops to 0.30 and the rms scatter rises by ~0.1 dex to 0.48. The concordance also drops slightly to *C* = 0.86. Excluding *f*<sub>esc,Lyα</sub> has the strongest negative impact on the concordance (*C* = 0.84) and the

second-largest impact on $R^2$ and rms. The concordance metric
includes nondetections, whereas the $R^2$ and rms parameters do
not. Hence, the larger effect of $f_{esc, Ly\alpha}$ on concordance may
illustrate that, even without information on dust, a lack of $Ly\alpha$
emission is generally sufficient to identify nonleakers in the
LzLCS+, even if the presence of $Ly\alpha$ does not necessarily
imply LyC escape. Conversely, predicting accurate $f_{esc}$ for
LCEs requires measurements of the line-of-sight dust
attenuation.

Substituting alternate measures of dust attenuation or ionization has little effect on the fiducial model. Using β<sub>1550</sub> instead of E(B − V)<sub>UV</sub> changes the goodness-of-fit metrics by only 0–0.03, consistent with the observed tight correlation between β<sub>1550</sub> and E(B − V)<sub>UV</sub> in the LzLCS+ (Chisholm et al. 2022). However, the β<sub>1200</sub> parameter is not as successful at tracing dust content (e.g., Chisholm et al. 2022) and worsens the fit quality by all metrics. Likewise, alternate measures of ionization such as [Ne III]/[O II] or [O II]/Hβ work as well as O32, with no change in C and minor (≤ 0.03) changes in R<sup>2</sup> and rms. Of the three measures of ionization, [O II]/Hβ performs best in all metrics, but this improvement is marginal.

3.3. Most Important Variables

As discussed in Section 2.4, we perform forward and
backward selection to determine which variables have the
greatest effect on the model fit quality. In Table 4, we show a
ranked ordering of variables based on forward and backward
selection for three representative models. The “Full Model”
includes β₁₅₅₀ instead of E(B − V)ₐᵥ and includes some
variables that measure similar but not highly collinear
properties, in order to test which ones perform best. For
instance, we include *f*<sub>esc,Lyα</sub>, EW(Lyα), and *L*(Lyα), which are
related but may have different relationships with *f*<sub>esc</sub>. Our
second set of variables is the more limited list from the Fiducial
+HI model, our best-performing Cox model. This model also
differs from the Full Model by using E(B − V)ₐᵥ. The third
model, the “JWST Model,” excludes absorption-line measure-
ments, which are difficult or impossible to measure for most
galaxies at *z* > 6, and Lyα measurements, which are heavily
affected by the IGM at *z* > 6; we discuss this model in
Section 3.4. We also list the mean MC ranks for each variable
obtained by sampling the observational uncertainties and
rerunning the ranking process 100 times. We plot the
distribution of these ranks in Figure 4.

We find that the most important factors for predicting *f*<sub>esc</sub> are almost always the strength of H I absorption (EW(H I,abs)) and dust attenuation (E(B − V)<sub>UV</sub> or β<sub>1550</sub>). These factors are naturally connected to the line-of-sight *f*<sub>esc</sub>, since the two sources of LyC absorption are H I atoms and dust. Previous studies have used the combination of H I absorption lines and dust to derive predictions for *f*<sub>esc* (e.g., Reddy et al. 2016; Chisholm et al. 2018; Saldana-Lopez et al. 2022), and the β<sub>1550</sub> parameter alone predicts much of the variation in *f*<sub>esc* (Chisholm et al. 2022). Importantly, unlike other tracers of H I content and dust such as Lyα and the nebular E(B − V), the H I absorption lines and UV measures of dust attenuation are more direct tracers of the H I and dust absorption along the line of sight to the UV-emitting stars and can thus substantially aid the prediction of the line-of-sight *f*<sub>esc</sub>. The combination of EW(H I,abs) and one of the UV dust measures alone achieves a concordance *C* = 0.84.
<table><caption>**Table 4**
Variables by Ranked Order</caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td>1</td><td>EW(H I,abs) [1.52]</td><td>EW(H I,abs) [3.91]</td><td>EW(H I,abs) [1.21]</td><td>f<sub>esc,Lyα</sub> [3.96]</td><td>β<sub>1550</sub> [2.55]</td><td>β<sub>1550</sub> [2.31]</td></tr><tr><td>2</td><td>β<sub>1550</sub> [3.42]</td><td>β<sub>1550</sub> [4.14]</td><td>E(B-V)<sub>UV</sub> [2.58]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [5.40]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)<sup>a</sup> [2.23]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)<sup>a</sup> [2.33]</td></tr><tr><td>3</td><td>M<sub>1500</sub> [9.41]</td><td>M<sub>1500</sub> [8.19]</td><td>log<sub>10</sub>(M<sub>*</sub>)<sup>a</sup> [7.32]</td><td>E(B-V)<sub>UV</sub> [3.79]</td><td>log<sub>10</sub>(O32) [2.19]</td><td>log<sub>10</sub>(O32) [2.05]</td></tr><tr><td>4</td><td>R<sub>f</sub>(H I,abs) [7.01]</td><td>R<sub>f</sub>(H I,abs) [7.86]</td><td>log<sub>10</sub>(EW(Hβ)) [7.39]</td><td>EW(H I,abs) [2.63]</td><td>log<sub>10</sub>(EW(Hβ)) [4.29]</td><td>log<sub>10</sub>(EW(Hβ)) [4.49]</td></tr><tr><td>5</td><td>EW(Lyα)<sup>a</sup> [8.82]</td><td>EW(Lyα)<sup>a</sup> [8.56]</td><td>f<sub>esc,Lyα</sub> [4.32]</td><td>log<sub>10</sub>(M<sub>*</sub>)<sup>a</sup> [7.65]</td><td>log<sub>10</sub>(M<sub>*</sub>) [6.04]</td><td>log<sub>10</sub>(M<sub>*</sub>) [5.50]</td></tr><tr><td>6</td><td>log<sub>10</sub>(EW(Hβ)) [9.18]</td><td>log<sub>10</sub>(M<sub>*</sub>) [9.72]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.43]</td><td>log<sub>10</sub>(EW(Hβ)) [6.90]</td><td>E(B-V)<sub>neb</sub> [6.04]</td><td>E(B-V)<sub>neb</sub> [6.11]</td></tr><tr><td>7</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [7.80]</td><td>R<sub>f</sub>(LIS) [10.21]</td><td>log<sub>10</sub>(O32) [6.52]</td><td>log<sub>10</sub>(O32) [4.99]</td><td>12+log<sub>10</sub>(O/H) [5.75]</td><td>12+log<sub>10</sub>(O/H) [6.16]</td></tr><tr><td>8</td><td>L(Lyα) [6.13]</td><td>EW(LIS) [9.72]</td><td>12+log<sub>10</sub>(O/H) [7.92]</td><td>12+log<sub>10</sub>(O/H) [7.34]</td><td>M<sub>1500</sub> [6.91]</td><td>M<sub>1500</sub> [7.05]</td></tr><tr><td>9</td><td>R<sub>f</sub>(LIS) [10.26]</td><td>f<sub>esc,Lyα</sub> [7.32]</td><td>E(B-V)<sub>neb</sub> [7.02]</td><td>E(B-V)<sub>neb</sub> [6.23]</td><td>...</td><td>...</td></tr><tr><td>10</td><td>log<sub>10</sub>(M<sub>*</sub>) [9.62]</td><td>12+log<sub>10</sub>(O/H) [9.66]</td><td>M<sub>1500</sub> [6.29]</td><td>M<sub>1500</sub> [6.11]</td><td>...</td><td>...</td></tr><tr><td>11</td><td>EW(LIS) [10.12]</td><td>E(B-V)<sub>neb</sub> [8.10]</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>12</td><td>f<sub>esc,Lyα</sub> [8.22]</td><td>log<sub>10</sub>(O32) [6.51]</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>13</td><td>12+log<sub>10</sub>(O/H) [10.13]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [9.14]</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>14</td><td>E(B-V)<sub>neb</sub> [9.07]</td><td>L(Lyα) [8.18]</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>15</td><td>log<sub>10</sub>(O32) [9.29]</td><td>log<sub>10</sub>(EW(Hβ)) [8.76]</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Notes. Variables in rank order from 1 (most important) to 15 (least important) by forward or backward selection. Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs. We plot the full distribution of ranks in the MC runs in Figure 4. The Full Model includes most variables in Table 1 aside from those that are highly collinear. The Fiducial+HI Model includes the variables in the best-performing Cox Model. The JWST Model includes only variables accessible at $z > 6$.

a After this variable, adding variables improves C by <0.01.

The MC method also finds that these two parameters, HI absorption and UV dust attenuation, have the top mean ranking. The MC mean ranks for the other variables are not necessarily consistent with their rank order in Table 4, which reflects the fact that the ranked order varies considerably among different MC runs (Figure 4). Since C barely changes after adding the first three to five variables, the order of most variables may be largely random, consistent with the broad ranking distribution in Figure 4. As also noted in Section 3.1, only the most significant variables dominate the model predictions, and limiting the model to this smaller set of variables gives comparable results to a model generated using a larger variable set. Consequently, we recommend using only the most statistically significant or top-ranked variables when using Cox models to predict $f_{esc}$.

Table 4 and Figure 4 show that Lyα measurements are often some of the more important variables, although generally not as important as EW(H I,abs). For all rankings that involve Lyα, including the MC rankings, an Lyα measurement, either EW(Lyα), L(Lyα), or $f_{esc,Lyα}$, achieves a rank between 1 and 5. In the MC runs, for the Fiducial+HI model, $f_{esc,Lyα}$ is the third most important variable after EW(H I,abs) and $E(B-V)_{UV}$. In Section 3.2, we found that excluding $f_{esc,Lyα}$ had one of the most detrimental effects on the fiducial model. Crucially, the fiducial model lacks EW(H I,abs), so $f_{esc,Lyα}$ provided essential information on the H I optical depth in its place. With EW(H I, abs) included, $f_{esc,Lyα}$ plays a lesser role but is still ranked among the top variables.

The rankings for the Full Model in Table 4 show that EW(H I,abs) is a better predictor of $f_{esc}$ in the LzLCS+ sample than the residual intensity $R_f(H I,abs)$. This ranking may result from the low resolution of the FUV spectra. Low spectral resolution, which depends on the spatial size of the source, will artificially increase the observed residual intensity, $R_f(H I,abs)$. In contrast, the EW(H I,abs) is less sensitive to spectral resolution and may serve as a more accurate measure of the

line-of-sight gas. Indeed, Saldana-Lopez et al. (2022) find that at a fixed $R_f$ stronger LCEs have lower EWs. At the same time, the observational uncertainties are also higher in $R_f$, such that $f_{esc}$ trends may appear more readily with EW. To investigate the effect of observational uncertainty on the rankings, we compare each variable's median uncertainty with the standard deviation of that variable in the LzLCS+ sample. If the measurement uncertainty is comparable to the standard deviation, we may not be able to discern trends with that variable across the LzLCS+. For most variables, the median uncertainty is much lower, ≲20% of the standard deviation. These variables span the full range of possible ranks and MC ranks, which suggests that we can distinguish variables that do and do not affect the $f_{esc}$ predictions. However, a few variables ($R_f(H I,abs)$, $log_{10}(M_*)$, EW(LIS), and $R_f(LIS)$) have higher ratios of median uncertainty to standard deviation (0.50–1.03). These variables all have high MC ranks in the Full Model (6.73–10.81), and any genuine trends with $f_{esc}$ may be hidden by the uncertainties in their measurements. As a test, we insert $log_{10}(f_{esc})$ as a dummy variable, as it should correlate perfectly with itself. Its MC rank begins to deviate from 1.00 when we give it an uncertainty ≳0.6 times the standard deviation. To further test the effect of uncertainty, we rerun the MC rankings after doubling the uncertainty in the EW(H I,abs) measurements, such that the EW(H I,abs) and $R_f(H I,abs)$ variables have similar ratios of uncertainty to standard deviation. The EW(H I, abs) MC ranks of 1.21–4.90 increase to 2.6–7.2 when the uncertainty is doubled. We conclude that the higher uncertainty in $R_f(H I,abs)$, $log_{10}(M_*)$, EW(LIS), and $R_f(LIS)$ may prevent us from determining the importance of these variables. Higher-resolution data are necessary to test whether EW(H I,abs) or $R_f(H I,abs)$ better predicts $f_{esc}$.

In their multivariate analysis of $f_{esc}$ from cosmological simulations, Maji et al. (2022) also use forward and backward selection to rank the importance of variables. They find that the three most important predictors of $f_{esc}$ are L(Lyα), the SFR, and
<figure><img src="image_9.png" /></figure>

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

Forward MC Rank Fiducial+HI

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

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

Forward MC Rank JWST Model

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

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

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

Backward MC Rank Fiducial+HI

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

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

Backward MC Rank JWST Model

<figure><img src="image_18.png" /><figcaption>**Figure 4.** The distribution of variable rankings from 100 forward selection (left) and backward selection (right) MC runs, after resampling each variable using its observational uncertainties. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.</figcaption></figure>

galaxy gas mass. Our rankings share some broad similarities with these simulation results. We also find that gas content (as measured by EW(HI,abs)), SFR (in the form of the UV

luminosity), and Lyα emerge as important parameters. However, our sample and the Maji et al. (2022) simulated galaxies have some crucial differences. First, we measure the
<figure><img src="image_19.png" /><figcaption>**Figure 5.** The $f_{\text{esc}}$ predictions from JWST models with (a) the full set of variables and (b) top-ranked variables only. Symbols are the same as in Figure 2. The scatter in these models is higher than in the fiducial model, but they still reproduce the observed $f_{\text{esc}}$ with an rms scatter of 0.46–0.47 dex.</figcaption></figure>

line-of-sight $f_{\text{esc}}$, whereas Maji et al. (2022) measure the total global $f_{\text{esc}}$ through all sight lines; consequently, the relevant parameter for our analysis is the line-of-sight H I rather than the global H I. Second, the galaxies in the Maji et al. (2022) sample are much less massive than the LzLCS+ sample, with a median stellar mass of $\log_{10}(M_*/M_\odot) = 6.41$ versus 8.8. Dust may play a greater role in determining $f_{\text{esc}}$ in the higher-mass, more enriched galaxies in our sample, which would account for its higher importance in our predictions.

### 3.4. Predictions at $z > 6$ with JWST

Although the fiducial model and the models using the top-ranked variables predict $f_{\text{esc}}$ reasonably well, we cannot apply these models to JWST observations of galaxies in the epoch of reionization. The Gunn–Peterson trough and Lyman series IGM absorption at lower redshifts will prevent H I absorption-line measurements, and the partially neutral IGM at $z > 6$ can suppress Ly$\alpha$ emission (e.g., Stark et al. 2011; Schenker et al. 2014). Measuring LIS absorption lines instead requires high-S/N observations of the rest-UV continuum, which will be difficult for faint galaxies. Consequently, we explore alternative models that use parameters that can be derived from $z > 6$ JWST observations.

We create a JWST model by modifying the fiducial model to exclude EW(LIS) and $f_{\text{esc},Ly\alpha}$. We also choose $\beta_{1550}$ instead of $E(B-V)_{\text{UV}}$, as it is more easily inferred from observations without any required stellar population modeling (e.g., Chisholm et al. 2022) or assumptions about the dust attenuation law. The JWST model thus includes the following variables: $\log_{10}(M_*)$, $M_{1500}$, $\log_{10}(\text{EW}(H\beta))$, $E(B-V)_{\text{neb}}$, $12+\log_{10}(\text{O/H})$, $\log_{10}(\text{O32})$, $\log_{10}(\Sigma_{\text{SFR}})$, and $\beta_{1550}$. We show the resulting Cox model fit in Figure 5(a) and list the goodness-of-fit metrics in Table 3. The scatter is noticeably higher by 0.1 dex than in the fiducial model, with rms = 0.47 dex. Like the fiducial model, the JWST model also tends to underpredict $f_{\text{esc}}$ in several of the strongest LCEs with $f_{\text{esc}} > 0.2$. This reduced fit quality shows that information about the line-of-sight H I content from absorption lines and Ly$\alpha$ is essential to precisely predict $f_{\text{esc}}$. Nevertheless, other observable properties can still provide rough $f_{\text{esc}}$ estimates and

distinguish strong LCEs from nonleakers. Since the JWST model relies more on global properties than line-of-sight H I to predict $f_{\text{esc}}$, the underprediction for the most extreme LCEs suggests that their global properties are not distinct from more moderate LCEs. The model-predicted $f_{\text{esc}}$ may indicate the typical $f_{\text{esc}}$ value for this combination of parameters. Moreover, if the most extreme $f_{\text{esc}}$ occurs only along favorable, nearly transparent sight lines, the $f_{\text{esc}}$ predicted from global parameters could be a better estimate of these galaxies’ global average $f_{\text{esc}}$. The strongest LCEs in the LzLCS+ sample do tend to have high nebular EWs (EW(H$\beta \gtrsim 200$ Å; Flury et al. 2022a), which shows that they cannot be devoid of absorbing gas in all directions.

Only three variables are statistically significant coefficients in the JWST model: $\beta_{1550}$, $\log_{10}(\Sigma_{\text{SFR}})$, and $\log_{10}(\text{O32})$. These same three parameters are ranked as the most important in forward and backward selection and have the lowest mean ranks from forward and backward selection with MC sampling (Table 4, Figure 4). When we limit the model to just these three variables, the fit quality is similar to or better than the JWST model with the full set of variables (Table 3). We also find that the $C_{\text{CV}}$ for the cross-validation analysis is closer to the C derived from the full sample; with fewer variables in the fit, even a smaller training sample can generate a reliable model. We display the predicted $f_{\text{esc}}$ from this model in Figure 5(b). Using [Ne III]/[O II] instead of O32 results in a similar quality fit (Table 3) and may be more suitable for observations covering a limited spectral range, observations at $z \gtrsim 11$ where [O III] is redshifted out of the NIRSpec wavelength range, or observations with uncertain nebular dust attenuation (Levesque & Richardson 2014). These fits demonstrate that the multivariate Cox model can predict $f_{\text{esc}}$ for high-redshift galaxies using JWST observables. We apply Cox models to high-redshift galaxy samples in a forthcoming paper (Jaskot et al. 2024), where we also provide all parameters needed to apply these models to future samples.

Our JWST model includes similar variables to the multi-variate $f_{\text{esc}}$ predictions derived by Choustikov et al. (2024) from the SPHINX simulations, but we find a different dependence on some of these variables. Choustikov et al. (2024) use the following observables: UV slope $\beta$, $E(B-V)_{\text{neb}}$, H$\beta$
luminosity, EW(Hβ), UV magnitude, R23 = ([O III] λλ5007, 4959 + [O II] λ3727)/Hβ, O32, and half-light radius. Like the Choustikov et al. (2024) model, we find that β1550 is statistically significant and that high $f_{esc}$ is associated with blue UV slopes. We likewise find that O32 is statistically significant. However, in the LzLCS+ sample, O32 correlates with $f_{esc}$, whereas the SPHINX simulations find an anticorrelation. Strong LCEs with high O32 do not appear within the SPHINX galaxy population (Choustikov et al. 2024). This disagreement may result from the different properties of the observed versus simulated galaxy populations, such as the stellar mass range. Alternatively, it may indicate the need for more efficient radiative feedback and/or resolving smaller-scale turbulent gas structure, to allow LyC escape from younger clusters (e.g., Kimm et al. 2019; Kakiichi & Gronke 2021; Choustikov et al. 2024).

Our model also disagrees with Choustikov et al. (2024) regarding the role of Σ<sub>SFR</sub>. Whereas we find that Σ<sub>SFR</sub> is one of the most important predictors of $f_{esc}$ in the JWST model and is statistically significant, Choustikov et al. (2024) find that compactness is insignificant and shows little relationship with $f_{esc}$. One of the main differences between the LzLCS+ and SPHINX galaxy data sets is the mass range probed; the LzLCS+ sample has a median mass of $log_{10}(M_*/M_\odot) = 8.8$ and extends up to $log_{10}(M_*/M_\odot) = 10.8$, whereas the majority of the SPHINX galaxies have masses $log_{10}(M_*/M_\odot) < 9$. At even lower masses, $log_{10}(M_*/M_\odot) < 7$, simulations suggest that LyC escape may occur in a more extended mode, driven by star clusters in low column density regions in the outskirts of galaxies (Kostyuk et al. 2023). The LzLCS+ likewise hints that the link between Σ<sub>SFR</sub> and $f_{esc}$ may be mass dependent. Several of the LCEs with high Σ<sub>SFR</sub> in the LzLCS+ sample have high masses ($log_{10}(M_*/M_\odot) = 9.7–10.5$), and as we discuss in the following section, Σ<sub>SFR</sub> may be a particularly significant predictor of $f_{esc}$ in high-mass galaxies.

### 3.5. Variations with Galaxy Properties

Flury et al. (2022a) discuss the possibility that the LyC escape process may vary among different types of galaxies. For example, the dominant feedback mechanism may shift from radiative feedback to mechanical supernova feedback as a starburst ages or in galaxies with different masses or metallicities (e.g., Jaskot et al. 2019; Kimm et al. 2019; Jecmen & Oey 2023). In addition, the strongest LCEs could have a nearly density-bounded gas geometry, whereas LyC photons may escape along narrow channels in weaker LCEs (e.g., Gazagnes et al. 2020; Flury et al. 2022a). We therefore explore whether the Cox model predictions change if we limit our analysis to different subpopulations of the larger LzLCS+ sample. We first divide our galaxy sample into two bins, above and below the median value of stellar mass, $log_{10}(M_*/M_\odot) = 8.8$.

Different variables appear important for the low- versus high-mass sample fits. For the fiducial model, Σ<sub>SFR</sub> is statistically significant in the high-mass fit only, with a coefficient value ~5× higher than that found in the low-mass fit. Conversely, M<sub>1500</sub> and O32 are only statistically significant in the low-mass fit, with coefficients ~8× and ~5× stronger than in the high-mass fit. The same pattern occurs with the JWST models, with Σ<sub>SFR</sub> significant in the high-mass fits and M<sub>1500</sub> and O32 significant in the low-mass fits.

The variable ranks by forward and backward selection (Tables 5–6, Figures 6–7) likewise suggest that different

variables are important for the two subsamples. As with the full sample (Table 4), variables tracing dust attenuation, H I absorption, and Lyα emission are some of the most important variables for both the high- and low-mass subsamples. However, for the high-mass subsample (Table 5), $log_{10}(\Sigma_{SFR})$ consistently ranks as the most important variable in all models and in its mean MC rank values. In the low-mass subsample (Table 6), $log_{10}(O32)$ actually appears at the bottom of the ranked list for the full set of variables (rank 15), but by mean MC rank, $log_{10}(O32)$ is one of the four most important variables in all models and is the most important variable in the JWST model runs. This seemingly discrepant rank may suggest that O32 provides similar information to some of the other top-ranked variables. Figure 7 shows that O32 replaces EW(H I, abs) as the top-ranked variable ~10%–25% of the time, and high O32 may generally imply weak H I absorption. In the JWST model, when EW(H I, abs) is not included, O32 takes over as the most important variable in a majority of the MC runs. To test whether the smaller sample size of the high- and low-mass subsets affects these results, we rerun the MC variable ranking using bootstrap resampling to change the selected sample. The resulting ranking distributions are similar to those in Figures 6–7, with only minor changes to the top four ranked variables. Our main conclusions remain unchanged. Σ<sub>SFR</sub> consistently holds the top MC rank for every model for the high-mass subset. Dust, H I absorption, Lyα, and O32 are important for the low-mass subset, and O32 is the top-ranked variable for the low-mass JWST model.

The separate high- and low-mass models do a slightly better job at reproducing the observed $f_{esc}$ values in LCEs than the full fiducial model with all galaxies included. The separate models show a slightly higher $R^2 = 0.64–0.70$ and lower rms = 0.28–0.32. However, the $R^2_{adj}$ values of 0.31–0.43 are lower than the full model $R^2_{adj} = 0.49$, which suggests that the improvement mainly comes from fitting a smaller sample. One exception is the JWST Model, where the high-mass subsample predictions are substantially improved compared to the full model; $R^2$ rises to 0.54 from 0.29, $R^2_{adj}$ is 0.23 instead of 0.14, the rms drops from 0.47 to 0.34, and C increases from 0.83 to 0.90. The Σ<sub>SFR</sub> variable takes on greater importance in the JWST model, since it lacks H I and Lyα information. The high-mass subsample $f_{esc}$ depends strongly on Σ<sub>SFR</sub>, and the best-fit coefficient for Σ<sub>SFR</sub> increases by a factor of four compared to the full sample JWST model. Evidently Σ<sub>SFR</sub> is key for accurately predicting $f_{esc}$ in the higher-mass galaxies within the LzLCS+.

Interestingly, the importance of Σ<sub>SFR</sub> versus O32 for different subsamples only emerges when we split the sample by stellar mass, not by UV luminosity. Almost all the highly ionized galaxies, with O32 ≳ 5, are in the lower stellar mass bin. Some of these galaxies can be quite luminous for their mass, so they are not as distinguishable by luminosity. The galaxies with the highest O32 ratios also have high nebular emission line EWs, suggestive of young ages. Lower-mass galaxies may have more bursty star formation histories in general (e.g., Lee et al. 2007). Among the low-mass galaxies in the LzLCS+, the high-O32, high-EW galaxies may represent currently bursting galaxies with extremely young average ages for their UV-emitting stellar populations (e.g., Izotov et al. 2011).

The noticeably different roles of the Σ<sub>SFR</sub> and O32 variables in the high- versus low-mass samples suggest that the causes of LyC escape may indeed differ between these galaxy
<table><caption>**Table 5**
Most Important Variables for High-mass Sample ($\log_{10}(M_*) \ge 8.8$)</caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top-ranked Variables</td></tr><tr><td>1</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td></tr><tr><td>2</td><td>f<sub>esc,Lyα</sub></td><td>f<sub>esc,Lyα</sub></td><td>f<sub>esc,Lyα</sub></td><td>f<sub>esc,Lyα</sub></td><td>β<sub>1550</sub></td><td>β<sub>1550</sub></td></tr><tr><td>3</td><td>β<sub>1550</sub></td><td>β<sub>1550</sub></td><td>E(B-V)<sub>UV</sub></td><td>E(B-V)<sub>UV</sub></td><td>log<sub>10</sub>(M<sub>*</sub>)</td><td>log<sub>10</sub>(M<sub>*</sub>)</td></tr><tr><td>4</td><td>R<sub>i</sub>(H I,abs)</td><td>R<sub>i</sub>(H I,abs)</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>M<sub>1500</sub></td><td>M<sub>1500</sub></td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top Variables by MC Rank Order</td></tr><tr><td></td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [3.24]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.35]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.79]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [2.98]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.10]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.17]</td></tr><tr><td></td><td>β<sub>1550</sub> [5.34]</td><td>log<sub>10</sub>(O32) [6.75]</td><td>f<sub>esc,Lyα</sub> [3.53]</td><td>EW(H I,abs) [4.33]</td><td>β<sub>1550</sub> [2.64]</td><td>β<sub>1550</sub> [2.93]</td></tr><tr><td></td><td>EW(H I,abs) [5.38]</td><td>L(Lyα) [7.17]</td><td>E(B-V)<sub>UV</sub> [3.68]</td><td>f<sub>esc,Lyα</sub> [4.37]</td><td>log<sub>10</sub>(M<sub>*</sub>) [5.06]</td><td>log<sub>10</sub>(O32) [4.61]</td></tr><tr><td></td><td>f<sub>esc,Lyα</sub> [6.18]</td><td>β<sub>1550</sub> [7.41]</td><td>EW(H I,abs) [3.99]</td><td>E(B-V)<sub>UV</sub> [4.76]</td><td>log<sub>10</sub>(EW(Hβ)) [5.36]</td><td>E(B-V)<sub>neb</sub> [5.14]</td></tr></tbody></table>

**Note.** Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

<table><caption>**Table 6**
Most Important Variables for Low-mass Sample ($\log_{10}(M_*) < 8.8$)</caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top-ranked Variables</td></tr><tr><td>1</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>β<sub>1550</sub></td><td>β<sub>1550</sub></td></tr><tr><td>2</td><td>β<sub>1550</sub></td><td>β<sub>1550</sub></td><td>E(B-V)<sub>UV</sub></td><td>E(B-V)<sub>UV</sub></td><td>log<sub>10</sub>(O32)</td><td>log<sub>10</sub>(O32)</td></tr><tr><td>3</td><td>f<sub>esc,Lyα</sub></td><td>M<sub>1500</sub></td><td>f<sub>esc,Lyα</sub></td><td>M<sub>1500</sub></td><td>M<sub>1500</sub></td><td>M<sub>1500</sub></td></tr><tr><td>4</td><td>log<sub>10</sub>(M<sub>*</sub>)</td><td>EW(Lyα)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(O32)</td><td>log<sub>10</sub>(EW(Hβ))</td><td>E(B-V)<sub>neb</sub></td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top Variables by MC Rank Order</td></tr><tr><td></td><td>EW(H I,abs) [3.08]</td><td>β<sub>1550</sub> [5.55]</td><td>EW(H I,abs) [2.09]</td><td>EW(H I,abs) [3.86]</td><td>log<sub>10</sub>(O32) [1.50]</td><td>log<sub>10</sub>(O32) [1.88]</td></tr><tr><td></td><td>β<sub>1550</sub> [4.70]</td><td>EW(H I,abs) [6.32]</td><td>E(B-V)<sub>UV</sub> [3.07]</td><td>E(B-V)<sub>UV</sub> [3.91]</td><td>β<sub>1550</sub> [3.02]</td><td>β<sub>1550</sub> [2.86]</td></tr><tr><td></td><td>log<sub>10</sub>(O32) [7.50]</td><td>log<sub>10</sub>(O32) [6.39]</td><td>f<sub>esc,Lyα</sub> [4.18]</td><td>log<sub>10</sub>(O32) [4.15]</td><td>M<sub>1500</sub> [4.38]</td><td>M<sub>1500</sub> [3.52]</td></tr><tr><td></td><td>R<sub>i</sub>(H I,abs) [7.89]</td><td>f<sub>esc,Lyα</sub> [7.00]</td><td>log<sub>10</sub>(O32) [4.89]</td><td>f<sub>esc,Lyα</sub> [4.52]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.85]</td><td>E(B-V)<sub>neb</sub> [5.06]</td></tr></tbody></table>

**Note.** Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

populations. In the low-mass sample, radiative feedback in young starbursts may drive high LyC escape over wide opening angles (e.g., Flury et al. 2022a). Galaxies with high luminosities and highly ionized gas would have higher $f_{esc}$. At the same time, if the youngest stellar population dominates the UV light output, the measured line-of-sight dust attenuation and HI absorption would also be representative of the conditions near the LyC-emitting stars and would be sufficient to constrain $f_{esc}$. In contrast, escape in higher-mass galaxies, with potentially longer-duration starburst episodes, may rely on the cumulative effect of supernova feedback punching holes in the ISM (e.g., Flury et al. 2022a). The ability to drive these required outflows may be linked to the galaxy's $\Sigma_{SFR}$ (e.g., Heckman et al. 2001; Kim et al. 2020). If the size of these holes is sufficiently small, the measured dust attenuation and HI EW across the full starburst may not reflect the conditions at the LyC escape site, especially if the youngest LyC-emitting populations do not dominate the total UV light. This disconnect could account for the reduced importance of $\beta_{1550}$ and EW(H I, abs) and enhanced importance of $\Sigma_{SFR}$ in the high-mass models.

To explore variations related to feedback, we also split the sample based on the median values of $\log_{10}(\Sigma_{SFR}) = 0.77$ and

$O32 = 3.4$. While splitting the LzLCS+ sample on mass resulted in slightly different variable selections, the high-O32 and high-$\Sigma_{SFR}$ subsample fiducial fits more closely resemble the full model, likely because they contain nearly the same population of strong LCEs; the high-$\Sigma_{SFR}$ and high-O32 subsets each include all but two of the strong LCEs with $f_{esc} > 0.1$. The two most important predictors of $f_{esc}$ remain HI absorption (either $R_i(H I, abs)$ or EW(H I, abs)) and dust attenuation. Interestingly, $\Sigma_{SFR}$ is no longer a statistically significant variable in the fiducial and JWST models for the high-$\Sigma_{SFR}$ and high-O32 subsamples. $\Sigma_{SFR}$ may be more useful for identifying nonleakers than in actually quantifying $f_{esc}$. Indeed, LCEs in the LzLCS+ sample nearly all lie above a threshold value of $\Sigma_{SFR} = 10 M_{\odot} \text{ yr}^{-1} \text{ kpc}^{-2}$, but above this value no discernible relationship exists between $f_{esc}$ and $\Sigma_{SFR}$ (Flury et al. 2022a).

For the low-$\Sigma_{SFR}$ subsample, only $f_{esc,Lyα}$ and $E(B-V)_{UV}$ are statistically significant in the fiducial model, while these same variables plus $E(B-V)_{neb}$ are significant for the low-O32 subset. As usual, in the variable rankings, EW(H I, abs) and dust measurements are typically important variables. However, Lyα measurements also appear to be relevant predictors. For the Full Model variable set, $L(Lyα)$ is the top-ranked variable for
<figure><img src="image_20.png" /><figcaption>**Figure 6.** The distribution of MC variable rankings from forward selection (left) and backward selection (right) for galaxies with $\log_{10}(M_*/M_\odot) \ge 8.8$. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.</figcaption></figure>

the low-O32 subsample. Similarly, for the low-$\Sigma_{SFR}$ sub-sample, $f_{esc,LY\alpha}$ is the top-ranked variable for the Fiducial+HI model and the top-ranked variable in forward modeling for the

Full Model. Because these subsamples contain few strong LCEs, this result emphasizes that Ly$\alpha$ emission may be useful in distinguishing weak LCEs from nonleakers.
<figure><img src="image_21.png" /><figcaption>**Figure 7.** The distribution of MC variable rankings from forward selection (left) and backward selection (right) for galaxies with $\log_{10}(M_*/M_\odot) < 8.8$. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.</figcaption></figure>

### 3.6. Outliers in $f_{esc}$ Trends

Many studies have identified relationships between $f_{esc}$ and physical or observable parameters (e.g., Verhamme et al. 2017;

Izotov et al. 2018b; Flury et al. 2022a; Chisholm et al. 2022; Saldana-Lopez et al. 2022), yet even the best trends show significant scatter, with individual galaxies showing $f_{esc}$
<figure><img src="image_22.png" /><figcaption>**Figure 8.** The $f_{\text{esc}}$ predictions from the fiducial model (a) modified to use EW(H I,abs) and (b) modified to exclude information about the line of sight from $Ly\alpha$, absorption lines, dust attenuation, or the observed UV luminosity. Black crosses show strong LCEs that are outliers by one or more selection criteria; gray circles and triangles show weak and non-LCEs that are outliers by one or more selection criteria. Small red circles and blue triangles represent the other detections and nondetections in the LzLCS+ sample, respectively. The fiducial+EW(H I,abs) model reproduces the $f_{\text{esc}}$ of the outliers, whereas the model without line-of-sight information does not.</figcaption></figure>

significantly above or below the expected value. Schaerer et al. (2022a) highlight one outlier galaxy in particular, J1248+4259. J1248+4259's high nebular ionization and strong $Ly\alpha$ emission resemble those of other strong low-redshift LCEs, yet its $f_{\text{esc}}$ is $\le 1.3\%$. Here we investigate whether the Cox model improves $f_{\text{esc}}$ predictions for J1248+4259 and other outlier galaxies.

We select outliers based on the variables that show the strongest trends with $f_{\text{esc}}$ in the Flury et al. (2022a) analysis of the LzLCS+ sample: $r_{50,\text{NUV}}$, O32, $\Sigma_{\text{SFR}}/M_*$, EW($Ly\alpha$), $\Sigma_{\text{SFR}}$, and $f_{\text{esc},Ly\alpha}$. From the strong LCEs ($f_{\text{esc}} \ge 0.05$) in the LzLCS+ sample, we choose the three with the highest $r_{50,\text{NUV}}$, the three lowest in O32, the three lowest in $\Sigma_{\text{SFR}}/M_*$, the three lowest in EW($Ly\alpha$), the three lowest in $\Sigma_{\text{SFR}}$, and the three lowest in $f_{\text{esc},Ly\alpha}$. Conversely, we choose three weak or non-LCEs ($f_{\text{esc}} < 0.05$) with the lowest $r_{50,\text{NUV}}$, highest O32, highest $\Sigma_{\text{SFR}}/M_*$, highest EW($Ly\alpha$), highest $\Sigma_{\text{SFR}}$, or highest $f_{\text{esc},Ly\alpha}$. This selection gives us a list of 10 strong LCEs and 13 weak or non-LCEs that are outliers by one or more selection. In other words, applying a simple, single-variable trend as in Flury et al. (2022a) would incorrectly identify these galaxies as strong LCEs or as nonleakers.

Figure 8(a) shows that the fiducial model modified to include EW(H I,abs) predicts $f_{\text{esc}}$ for these “outliers” quite well. In fact, the outlier predictions are generally no more or less accurate than the predictions for other galaxies in the LzLCS+ sample. For example, the outlier galaxy J1248+4259, which “should” be a strong LCE based on its high O32 ratio and high EW($Ly\alpha$), has a predicted $f_{\text{esc}}$ of 0.006 according to the model, consistent with its observed limit of $f_{\text{esc}} \le 0.013$.

Crucially, the fiducial+EW(H I,abs) model includes information about the optical depth along the line of sight, which is necessary to predict the $f_{\text{esc}}$ observed along the same line of sight. As discussed in Sections 3.2 and 3.3, EW(H I,abs), $f_{\text{esc}}$, $Ly\alpha$, and $E(B-V)_{\text{UV}}$, all of which are sight line dependent, have some of the largest effects on the accuracy of the model predictions. If we remove these and other sight-line-dependent parameters ($E(B-V)_{\text{neb}}$, $M_{1500}$) from the fit (Figure 8(b)), we

can no longer reproduce the line-of-sight $f_{\text{esc}}$ of many galaxies. Figure 8(b) indicates that we would expect a different $f_{\text{esc}}$ for these galaxies based on their global properties. Hence, if the only cause of the scatter in Figure 8(b) is chance galaxy orientation, the predicted $f_{\text{esc}}$ should represent the global average across all sight lines. By this model, J1248+4259 should indeed be an LCE, although not an exceptionally strong one, with $f_{\text{esc}} = 0.051$.

At high $f_{\text{esc}}$, the model without line-of-sight information seems to systematically underpredict the observed $f_{\text{esc}}$, which may suggest that the global properties of the strongest LCEs are not distinct from those of weaker LCEs. The observed properties of the strongest LCEs support this interpretation. For example, the strongest LCEs all tend to have high O32 (Flury et al. 2022a), yet the median observed $f_{\text{esc}}$ for the 11 galaxies with O32 $\ge 10$ is only 0.05. For these same 11 galaxies, the model without line-of-sight information finds the median of their predicted $f_{\text{esc}}$ values to be 0.04, comparable to the observations. The only model that successfully reproduces the $f_{\text{esc}}$ of the strongest LCEs is the fiducial+EW(H I,abs) model. Hence, of the variables we consider, only EW(H I,abs) distinguishes the strongest LCEs from other galaxies. Either the strongest LCEs differ in some other unknown global property, or their high $f_{\text{esc}}$ results from a favorable orientation compared to weaker LCEs. The high-O32 galaxies may have $f_{\text{esc}} \sim 0.04-0.05$ along most sight lines but may also contain several extremely optically thin channels, which lead to much higher observed $f_{\text{esc}}$, $\gtrsim 0.5$ in some cases.

Although orientation bias likely causes some of the scatter in $f_{\text{esc}}$ predictions, the Cox modeling results also suggest that some of the “outlier” galaxies are not outliers at all, at least not when multiple global variables are considered. For example, J0919+4906 and J0232-0426 have high O32 ratios of 12.7 and 9.9 but $f_{\text{esc}}$ of 0.049 and $<0.047$, respectively. However, both galaxies also have low stellar masses ($\log_{10}(M_*) = 7.5$ and 7.4), such that their O32 values are not extreme compared to other galaxies of similar mass (see Figure 23 in Flury et al. 2022a). Furthermore, other criteria, such as their moderate or low $\Sigma_{\text{SFR}}$
values ($\log_{10}(\Sigma_{SFR}) = 0.66$ and $-0.39$), would argue against a high $f_{esc}$. The combination of multiple properties leads to a moderate predicted $f_{esc}$ for J0919+4906 ($f_{esc} = 0.048$) and little escaping LyC ($f_{esc} = 0.005$) for J0232-0426, even when the model includes only global, and not sight-line-dependent, properties. Similarly, several galaxies have concentrated star formation (J0910+6105, J1209+3053, J1248+1234, J1349+5631, J1503+3644, J1517+3705) or strong Ly$\alpha$ (J0122+0520, J0811+4141, J0901+2119, J1648+4957) but $f_{esc} < 0.05$; these galaxies' other properties yield consistent predictions of $f_{esc} < 0.05$, even in models that neglect line-of-sight information.

The Cox model's success at reproducing $f_{esc}$ even for these outliers demonstrates that accurately predicting $f_{esc}$ requires multiple input variables. Moreover, because the observed $f_{esc}$ depends on our line of sight, the best models also require sight-line-specific information about optical depth from absorption lines, UV dust attenuation, and Ly$\alpha$. Although Ly$\alpha$ and H I absorption lines are inaccessible for most $z > 6$ galaxies, high-redshift $f_{esc}$ predictions can include information about UV dust attenuation as in the JWST models discussed in Section 3.4. Information about dust attenuation is critical to predict $f_{esc}$, and large sample sizes can probe multiple random orientations to estimate the average $f_{esc}$ for a galaxy population (e.g., Cen & Kimm 2015; Saldana-Lopez et al. 2023).

# 4. Predicting Alternative Quantities

## 4.1. Alternative Measures of LyC: Predicting $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ and $L(LyC)$

As shown above, multivariate Cox models can successfully predict $f_{esc}$; we now consider whether they can predict alternative measurements of the escaping LyC. To model how galaxies reionized the Universe, we need to know the input rate of LyC photons into the IGM from different galaxy populations. Deriving this quantity from $f_{esc}$ requires knowing the production rate of LyC photons in each galaxy, which is typically derived from SED fitting or Balmer line observations. To avoid this process and its associated systematics, we could instead predict a more directly observable quantity, such as the $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ ratio or the total escaping LyC luminosity $\log_{10}(L(\text{LyC}))$. The $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ flux ratio measures the ratio of observed LyC flux to the observed flux at 1100 Å. As such, it depends on observed quantities and does not require an assumed stellar population model or dust attenuation law. To predict $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ using the Cox proportional hazards model, we follow the same procedure as for $f_{esc}$ and convert the $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ measurements to right-censored data (i.e., detections and lower limits) by predicting $1 - F_{\lambda_{LyC}}/F_{\lambda_{1100}}$. Unlike $f_{esc}$ and $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$, which cannot be higher than 1, $L(\text{LyC})$ has no obvious maximum allowed value. To convert it to a right-censored format, we therefore adopt a maximum of $\log_{10}(L(\text{LyC})) = 41$, greater than the maximum luminosity in the sample of $\log_{10}(L(\text{LyC})) = 40.6$, and predict $41 - \log_{10}(L(\text{LyC})) = 0.4$.

The $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ models perform worse than the $f_{esc}$ models for every variable combination we consider (see Figure 9) and by all three metrics ($R^2$, rms, C). Although $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ is simpler to measure, $f_{esc}$ appears simpler to predict. It may be more directly linked to galaxy physical

<table><caption>**Table 7**
Goodness-of-fit Metrics for Different Dependent Variables</caption><thead><tr><th>Dependent Variable</th><th>R<sup>2</sup></th><th>R<sup>2</sup><sub>adj</sub></th><th>rms</th><th>C</th><th>R<sup>2</sup><sub>CV</sub></th><th>rms<sub>CV</sub></th><th>C<sub>CV</sub></th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Fiducial Model<sup>a</sup></td></tr><tr><td>f<sub>esc</sub></td><td>0.60</td><td>0.49</td><td>0.36</td><td>0.89</td><td>0.55</td><td>0.36</td><td>0.85</td></tr><tr><td>F<sub>λ<sub>LyC</sub></sub>/F<sub>λ<sub>1100</sub></sub></td><td>0.07</td><td>-0.22</td><td>0.40</td><td>0.85</td><td>-0.04</td><td>0.40</td><td>0.81</td></tr><tr><td>L(LyC)</td><td>-0.20</td><td>-0.55</td><td>0.38</td><td>0.81</td><td>-0.33</td><td>0.37</td><td>0.76</td></tr><tr><td>f<sub>esc,Lyα</sub></td><td>0.33</td><td>0.24</td><td>0.32</td><td>0.75</td><td>0.29</td><td>0.32</td><td>0.71</td></tr><tr><td>L(Lyα)<sup>b</sup></td><td>0.58</td><td>0.53</td><td>0.23</td><td>0.76</td><td>0.55</td><td>0.23</td><td>0.71</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">JWST Model</td></tr><tr><td>f<sub>esc</sub></td><td>0.29</td><td>0.14</td><td>0.47</td><td>0.83</td><td>0.21</td><td>0.47</td><td>0.79</td></tr><tr><td>F<sub>λ<sub>LyC</sub></sub>/F<sub>λ<sub>1100</sub></sub></td><td>-0.18</td><td>-0.48</td><td>0.49</td><td>0.77</td><td>-0.38</td><td>0.48</td><td>0.73</td></tr><tr><td>L(LyC)</td><td>-0.25</td><td>-0.53</td><td>0.39</td><td>0.78</td><td>-0.34</td><td>0.38</td><td>0.74</td></tr><tr><td>f<sub>esc,Lyα</sub></td><td>0.07</td><td>-0.04</td><td>0.38</td><td>0.75</td><td>-0.01</td><td>0.37</td><td>0.72</td></tr><tr><td>L(Lyα)<sup>b</sup></td><td>0.53</td><td>0.48</td><td>0.25</td><td>0.77</td><td>0.50</td><td>0.24</td><td>0.73</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The $f_{esc,Lyα}$ variable is not included in the fiducial model for the $f_{esc,Lyα}$ and $L(Lyα)$ predictions.

<sup>b</sup> The $L(Lyα)$ model sample excludes six galaxies with Ly$\alpha$ in absorption.

properties such as optical depth, porosity, and feedback. Furthermore, while measuring $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ does not depend on assumptions about the galaxy SED, predicting it does. The galaxy SED shape changes with age and metallicity, which affects the ratio of $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ even for $f_{esc} = 1$ (see, e.g., Chisholm et al. 2019).

Although we derive both $E(B-V)_{UV}$ and $f_{esc}$ from the COS UV spectra, the greater accuracy of the $f_{esc}$ Cox models does not appear to result from this interdependency. First, we see the same improvement in accuracy for $f_{esc}$ versus $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ when we use $\beta_{1550}$ instead of $E(B-V)_{UV}$; $\beta_{1550}$ is directly related to the observed spectral slope and is not sensitive to the choice of stellar population model or star formation history (Chisholm et al. 2022). Second, the $f_{esc}$ models have higher $R^2$ and C than the $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ predictions, even for the models that contain no UV spectral information at all. We conclude that although $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ is a useful quantity, accurately predicting it is difficult given its dependence on the intrinsic SED shape. We list and discuss the top-ranked variables for the $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ predictions in the Appendix.

The $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ predictions rely on much the same variables as the $f_{esc}$ predictions, except for a weaker dependence on the UV dust attenuation, as expected. For the fiducial model, the same variables are statistically significant except for $E(B-V)_{UV}$, and excluding $E(B-V)_{UV}$ from the fiducial model changes $R^2$, rms, and C by only $\le 0.03$. However, when we substitute EW(H I,abs) in the fiducial model, $E(B-V)_{UV}$ does appear statistically significant, indicating that its influence is not always negligible. Table 8 in the Appendix lists the top-ranked variables and likewise indicates a reduced role for dust attenuation in predicting $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ compared to its role in predicting $f_{esc}$.

Like $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$, the predictions for $L(LyC)$ (Figure 10, Table 7) are also less accurate than for the $f_{esc}$ predictions. The concordance drops for every model, and $R^2$ is negative in all cases, meaning that simply adopting the mean $L(LyC)$ for the detections would result in a better fit than using the model predictions. The negative $R^2$ arises from the fact that the $R^2$
<figure><img src="image_23.png" /><figcaption>**Figure 9.** Examples of $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ model predictions using the variable sets for (a) the fiducial model and (b) the JWST model. Symbols are the same as in Figure 2. Predictions of $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ are not as accurate as the model predictions for $f_{esc}$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 10.** Examples of $L(\text{LyC})$ model predictions using the variable sets for (a) the fiducial model and (b) the JWST model. We plot galaxies with predicted $L(\text{LyC}) \sim 0$ at the bottom of each panel. Symbols are the same as in Figure 2. The models fail to accurately predict $L(\text{LyC})$.</figcaption></figure>

calculation (Equation (8)) uses only the detections and their mean, whereas the Cox model fits to the full set of galaxies; $R^2$ becomes positive if we use the median $L(\text{LyC})$ from the full sample of limits and detections, instead of the mean $L(\text{LyC})$ of the detections alone. Not surprisingly, other measures of luminosity show up as important variables for the $L(\text{LyC})$ models. $M_{1500}$ is statistically significant in the fiducial model and is one of the top three ranked parameters by forward and backward selection, while $L(\text{Ly}\alpha)$ holds the top place in forward selection for the full variable list (see Table 9 in the Appendix for a list of top-ranked variables). However, the $L(\text{LyC})$ models’ poor performance indicates that key information is missing. Specifically, a lack of detailed information about the SED likely prevents an estimate of the intrinsic $L(\text{LyC})$ from $M_{1500}$. The escaping $L(\text{Ly}\alpha)$ also imperfectly traces the escaping $L(\text{LyC})$ because of scattering. In addition, the production of the Ly$\alpha$ recombination line requires LyC

absorption, which could cause $L(\text{Ly}\alpha)$ to drop for the few extreme LCEs in the LzLCS+ sample (e.g., Nakajima & Ouchi 2014). In conclusion, neither the $L(\text{LyC})$ nor the $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ models provide satisfactory alternatives to the $f_{esc}$ predictions.

## 4.2. Predicting $f_{esc, Ly\alpha}$

Although closely related to the H I optical depth, the escape of Ly$\alpha$ may depend on different parameters than the LyC escape fraction does. Here we apply the Cox model to the $f_{esc, Ly\alpha}$ data for the LzLCS+ sample and assess the differences between predictions for $f_{esc, Ly\alpha}$ versus $f_{esc}$. We use the same input variables, with one exception: we exclude any variables containing Ly$\alpha$ measurements ($f_{esc, Ly\alpha}$, EW(Ly$\alpha$), $L(\text{Ly}\alpha)$) from the models. In addition, unlike our $f_{esc}$ data, the $f_{esc, Ly\alpha}$ measurements include no upper limits, since Ly$\alpha$ can appear in absorption or emission. For the negative values of $f_{esc, Ly\alpha}$, representing net absorption, we cannot calculate $\log_{10}(f_{esc, Ly\alpha})$
<figure><img src="image_25.png" /><figcaption>**Figure 11.** Examples of $f_{esc, Ly\alpha}$ model predictions using the variable sets for (a) the fiducial model, (b) the fiducial model with EW(H I, abs), and (c) the JWST model. Symbols are the same as in Figure 2. We plot galaxies with negative $f_{esc, Ly\alpha}$ at $f_{esc, Ly\alpha} = 2\text{E}-3$. Of the three models, the fiducial model with EW(H I, abs) predicts $f_{esc, Ly\alpha}$ with the lowest rms scatter.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 12.** Examples of $L(Ly\alpha)$ model predictions using the variable sets for (a) the fiducial model and (b) the JWST model. Symbols are the same as in Figure 2. We exclude galaxies with negative $L(Ly\alpha)$. The models do a better job at predicting $L(Ly\alpha)$ compared to $f_{esc, Ly\alpha}$.</figcaption></figure>

and do not include these galaxies in the reported $R^2$ or rms; they are included in the concordance.

The fiducial model is not very successful at predicting $f_{esc, Ly\alpha}$ (Figure 11(a), Table 7), with $R^2 = 0.33$, $R_{adj}^2 = 0.24$, rms = 0.32, and $C = 0.75$. The lack of upper limits and need to precisely rank galaxies with even very low $f_{esc, Ly\alpha}$ may explain the reduced $C$ values for the $f_{esc, Ly\alpha}$ model. Including EW(H I, abs) improves the $f_{esc, Ly\alpha}$ predictions ($R^2 = 0.44$, rms = 0.30, $C = 0.78$; Figure 11(b)). As with the $f_{esc}$ predictions, the JWST $f_{esc, Ly\alpha}$ model performs worse than the fiducial model (Figure 11(c)), which demonstrates that the observed $f_{esc, Ly\alpha}$ does depend on line-of-sight H I properties such as the neutral gas covering fraction or H I column density. Unlike the $f_{esc}$ models, for $f_{esc, Ly\alpha}$ the fiducial model without $E(B-V)_{UV}$ performs comparably to the fiducial model ($R^2 = 0.37$, rms = 0.29, $C = 0.75$), because $Ly\alpha$ is more sensitive to the nebular attenuation than to the stellar attenuation.

We also investigate $L(Ly\alpha)$, using a similar method to that for $L(LyC)$ and predicting $43.5 - \log_{10}(L(Ly\alpha))$. We exclude all galaxies with net $Ly\alpha$ absorption and predict only the luminosities of $Ly\alpha$ emitters. Whereas we found that $L(LyC)$

was harder to predict than $f_{esc}$, for $L(Ly\alpha)$ the situation is reversed: the $L(Ly\alpha)$ models are actually slightly better than the $f_{esc, Ly\alpha}$ models (Table 7 and Figure 12). The $L(LyC)$ predictions suffered from the large uncertainty in estimating the intrinsic LyC production from the FUV continuum. The available Balmer line measurements make estimating the intrinsic $Ly\alpha$ luminosity less of a concern. Another explanation for the better performance of the $L(Ly\alpha)$ predictions may be mundane; since it depends on predicting the intrinsic $Ly\alpha$, $f_{esc, Ly\alpha}$ is a more model-dependent measurement than $L(Ly\alpha)$ and consequently has a greater uncertainty (26% on average for $f_{esc, Ly\alpha}$ vs. 21% for the $Ly\alpha$ flux alone). Again, the $L(Ly\alpha)$ JWST model (Table 7 and Figure 12(b)) has slightly worse $R^2$ and rms metrics than the fiducial model because of the lack of information about line-of-sight gas.

Tables 10 and 11 in the Appendix list the top-ranked variables for predicting $f_{esc, Ly\alpha}$ and $L(Ly\alpha)$. For $f_{esc, Ly\alpha}$, we find that H I absorption strength and nebular dust attenuation are two of the most important predictive variables. In contrast, the UV dust attenuation is less important, always ranking after $E(B-V)_{neb}$ in forward and backward selection. As a nebular
emission line, Ly$\alpha$ escape is more closely linked to the nebular attenuation rather than the attenuation experienced by the stellar light.

For $L(\text{Ly}\alpha)$, the most important variables include $\Sigma_{\text{SFR}}$, $M_{1500}$, $E(B-V)_{\text{neb}}$, EW(H$\beta$), and EW(HI,abs). These five variables are also statistically significant in the fiducial model ($\Sigma_{\text{SFR}}$, the fiducial+EW(HI,abs) model (EW(H$\beta$), EW(HI, abs)), or both ($E(B-V)_{\text{neb}}$ and $M_{1500}$). The dependence on UV luminosity reflects the fact that higher SFRs enhance Ly$\alpha$ production, while the dependence on $E(B-V)_{\text{neb}}$ and EW(HI, abs) arises because lower dust attenuation and lower HI column densities enhance Ly$\alpha$ escape. As a recombination line, the H$\beta$ EW is also directly linked to nebular Ly$\alpha$ production. While the dependence on $\Sigma_{\text{SFR}}$ may reflect the role of feedback in Ly$\alpha$ escape (Kim et al. 2021), it may also serve to quantify aperture losses, as more extended galaxies may have more of their Ly$\alpha$ emission scatter outside the COS aperture.

Runnholm et al. (2020) also assess the most important variables in predicting $L(\text{Ly}\alpha)$ using the LARS. They find that the top five ranked variables by both forward and backward selection are SFR, $E(B-V)_{\text{neb}}$, $M_*$, $r_{50}$, and the gas covering fraction. If using observables, the top variables are the FUV luminosity and the UV size. Despite the different methodologies between the LARS Ly$\alpha$ measurements and the LzLCS+ Ly$\alpha$ measurements, our results generally agree with their findings, as $\Sigma_{\text{SFR}}$, $E(B-V)_{\text{neb}}$, and the galaxy FUV luminosity also appear among our most important predictors.

## 5. Discussion

Both theoretical and observational studies have long recognized that LyC escape is likely anisotropic owing to the inherently inhomogeneous structure of the ISM (e.g., Gnedin et al. 2008; Wise & Cen 2009; Kim et al. 2013, 2023; Zastrow et al. 2013; Cen & Kimm 2015; Paardekooper et al. 2015; Keenan et al. 2017; Rivera-Thorsen et al. 2017, 2019; Trebitsch et al. 2017). LyC may escape in a “picket-fence” geometry, emerging only through holes in the ISM (e.g., Heckman et al. 2001), and simulations find that supernova feedback may carve these low-density paths in the ISM (e.g., Wise & Cen 2009; Cen & Kimm 2015; Paardekooper et al. 2015; Trebitsch et al. 2017). Observationally, the detection of residual flux from saturated low-ionization absorption lines (e.g., Heckman et al. 2001, 2011; Gazagnes et al. 2018) and the coexistence of narrow Ly$\alpha$ emission with underlying broad absorption (e.g., McKinney et al. 2019) provide evidence for optical depth variations within galaxies. Highly ionized channels observed in some low-redshift galaxies (e.g., Zastrow et al. 2013; Bik et al. 2018) and spatially resolved LyC observations in the lensed Sunburst Arc at $z=2.37$ (Rivera-Thorsen et al. 2019; Kim et al. 2023) likewise imply that LyC may escape through narrow opening angles, perhaps driven by outflows (e.g., Amorín et al. 2024). Even when $f_{\text{esc}}$ is large, with escape occurring in most or all directions, HI column density variations will still lead to variations in the observed $f_{\text{esc}}$ along different lines of sight (e.g., Gazagnes et al. 2020).

Because of the anisotropy of LyC escape, estimating global $f_{\text{esc}}$ values from line-of-sight LyC detections is difficult, particularly since the observed $f_{\text{esc}}$ depends both on galaxy properties and on orientation. The multivariate Cox models present new evidence for the anisotropy of LyC escape by showing that the global properties of the strongest LCEs ($f_{\text{esc}} \gtrsim$

0.4) do not distinguish them from weaker LCEs ($f_{\text{esc}} \sim 0.04$; see Section 3.6). At $z \sim 3$, Nakajima et al. (2020) reach a similar conclusion, finding that high-O32 galaxies can have wide variations in $f_{\text{esc}}$, despite sharing comparable spectral properties. We find that the only variable that successfully predicts the high observed $f_{\text{esc}}$ of the strongest LCEs is EW(HI, abs), which traces the line-of-sight HI absorption. This result supports the picture described by Gazagnes et al. (2020), where galaxies with high O32 may have sufficient ionizing radiation to ionize most lines of sight, but where the highest $f_{\text{esc}}$ occurs in only the most diffuse channels.

Multivariate survival analysis models such as those presented here may offer a tool for quantifying the line-of-sight variation in $f_{\text{esc}}$ among galaxy populations. Given a large reference sample of galaxies with measurements of relevant global properties, the Cox models generate the probability distribution of $f_{\text{esc}}$ for a given set of variables. If the variation in the observed $f_{\text{esc}}$ for similar galaxies results only from orientation effects, the derived probability distributions should reflect the distribution of $f_{\text{esc}}$ across all sight lines.

Although the outsized importance of EW(HI,abs) in accurately predicting high observed $f_{\text{esc}}$ suggests that favorable orientation accounts for the strongest LCEs in the LzLCS+, orientation alone may not fully explain the scatter between predicted and observed $f_{\text{esc}}$ across the entire LzLCS+ sample. LyC escape will depend on local properties, which global variables may not accurately capture (e.g., Rivera-Thorsen et al. 2019). Observational studies at higher spatial resolution and simulations that track the evolution of $f_{\text{esc}}$ with high spatial resolution will help reveal the origin of $f_{\text{esc}}$ variations in otherwise similar galaxies.

The variable rankings in Section 3.3 show that the two most essential variables in predicting the line-of-sight $f_{\text{esc}}$ are variables sensitive to the line-of-sight HI gas (EW(HI,abs)) and the line-of-sight dust ($E(B-V)_{\text{UV}}$ or $\beta_{1550}$). Although $z>6$ observations will not be able to measure the former, $\beta_{1550}$ is easily accessible. JWST observations of the epoch of reionization are rapidly expanding the available measurements of $\beta_{1550}$ at high redshift (e.g., Schaerer et al. 2022b; Mascia et al. 2023; Morishita et al. 2024; Saxena et al. 2024), and this variable alone can provide useful constraints on $f_{\text{esc}}$ (Chisholm et al. 2022). By examining the spread of $\beta_{1550}$ values for a fixed intrinsic $M_{1500}$ or set of galaxy properties, high-redshift studies could gain insight into the possible variation in the line-of-sight $f_{\text{esc}}$ among a population, at least the variation due to changing levels of dust attenuation. To understand the HI component of the LyC absorption and its variation, high-redshift studies will need to rely on simulations and lower-redshift observations.

Our models provide some clues as to which physical properties most affect a galaxy's HI distribution and resulting $f_{\text{esc}}$. In the absence of EW(HI,abs) and Ly$\alpha$ measurements, which more directly trace HI, the models find that $\Sigma_{\text{SFR}}$ and O32 are the most important variables in predicting $f_{\text{esc}}$, given constraints on the dust attenuation from $\beta_{1550}$ (Section 3.4). These variables suggest that feedback is key in determining a galaxy's average HI column density. Interestingly, we find in Section 3.5 that O32 may be more important in lower-mass galaxies ($\log_{10}(M_*/M_{\odot}) < 8.8$), whereas $\Sigma_{\text{SFR}}$ takes on greater importance in the higher-mass LzLCS+ galaxies. Previous studies have suggested that the mechanism responsible for LyC escape may be different in different galaxy populations (e.g., Flury et al. 2022a; Katz et al. 2023). Here we find that
differences in the multivariate $f_{esc}$ models are specifically associated with galaxy mass, rather than luminosity. We hypothesize that the typical star formation histories of the low- and high-mass LzLCS+ galaxies may differ, with the lower-mass starbursts exhibiting burstier star formation (e.g., Lee et al. 2007). Radiative feedback in extremely young, strong starbursts may drive LyC escape over wide opening angles in these low-mass galaxies, whereas the cumulative effect of supernovae in longer-lived concentrated starbursts could regulate LyC escape at higher masses (e.g., Bremer & Dayal 2023).

While the same physical processes, such as bursty star formation with radiative feedback or an extended star formation episode with supernova feedback, may occur at high redshift, they may not occur in exactly the same type of galaxies as the LzLCS+. In particular, at a given stellar mass, a high-redshift galaxy could have a different gas mass, total mass, or star formation history compared to the same stellar mass galaxy at low redshift. Nevertheless, the distinct models for the low-versus high-mass subsets in the LzLCS+ demonstrate that different parameters and physical mechanisms may regulate LyC escape in different types of galaxies. These differences could explain why some studies disagree regarding the properties of LCEs (e.g., Alexandroff et al. 2015; Rutkowski et al. 2017; Izotov et al. 2018b; Naidu et al. 2018). Depending on the sample probed, high O32 or high $\Sigma_{SFR}$ may or may not be a relevant predictor of $f_{esc}$. Understanding the origin of the relationships between O32, $\Sigma_{SFR}$, and $f_{esc}$ in low- versus high-mass LzLCS+ galaxies will require more observational and theoretical studies examining the interplay between star formation history, geometry, and feedback in these types of galaxies.

## 6. Summary

LyC escape depends on a range of galaxy properties, from gas and dust distributions in the ISM to galaxy star formation histories and the resulting stellar feedback. Because of this complexity, $f_{esc}$ correlates observationally with a variety of parameters, and correlations with individual variables show high scatter. Using the results of the LzLCS and other $z \sim 0.3$ LyC observations from the literature, we have generated new multivariate models to predict $f_{esc}$. We use the statistical technique of survival analysis, specifically the Cox proportional hazards model (Cox 1972), which accounts for data from both LyC detections and $f_{esc}$ upper limits from nondetections. With these models, we explore which combinations of variables best predict $f_{esc}$ and discuss the possible physical basis for these trends.

We summarize our main findings below.

1. The $f_{esc}$ predictions from our fiducial model match the observed $f_{esc}$ of $z \sim 0.3$ LCEs with an rms of 0.36 dex, and our best-performing model, which includes H I absorption lines, reaches an rms of 0.31 dex. These multivariate models offer a substantial improvement over correlations with single variables, which can span $\gtrsim 2$ dex in $f_{esc}$ for a given galaxy property (e.g., Wang et al. 2021; Flury et al. 2022a; Chisholm et al. 2022; Saldana-Lopez et al. 2022). (See Sections 3.1–3.2.)

2. Our best-performing model accurately predicts $f_{esc}$ even for galaxies that are outliers in a single-variable selection. This model, which includes H I absorption-line

measurements, is the only model that reproduces the high observed $f_{esc}$ values of the strongest LCEs in the LzLCS+. Among the variables we consider, weak H I absorption is the only distinctive property of these strong LCEs, which suggests that these high $f_{esc}$ values may occur only along fortuitous sight lines. (See Section 3.6.)

3. The most important variables in the models, as determined by forward and backward selection, are the EW of the H I Lyman series absorption lines and either $E(B-V)_{UV}$ or the UV slope $\beta_{1550}$. These properties trace the line-of-sight H I absorption and UV dust attenuation, which determine the line-of-sight $f_{esc}$ that we measure (e.g., Saldana-Lopez et al. 2022). In our fiducial model, $E(B-V)_{UV}$ is a statistically significant variable, along with $f_{esc, Ly\alpha}$, $\Sigma_{SFR}$, $E(B-V)_{neb}$, and O32. Excluding $E(B-V)_{UV}$ or $f_{esc, Ly\alpha}$ has the greatest adverse effects on model performance, and including $f_{esc, Ly\alpha}$ appears to be particularly important for identifying non-LyC emitters. These analyses show which variables are the most important to include for accurate predictions; models using only the statistically significant variables or the top three to five ranked variables achieve comparable accuracy to models with a longer list of inputs (See Sections 3.2–3.3.)

4. We generate a “JWST model” for predicting $f_{esc}$, which can apply to $z > 6$ observations. The model excludes $Ly\alpha$ and absorption-line measurements, which will be inaccessible for most $z > 6$ galaxies. When compared with the $z \sim 0.3$ LCEs, this model has a higher scatter (rms = 0.47 dex) than the fiducial model, but it can distinguish LCEs from nonemitters. Three variables in this model are statistically significant: $\beta_{1550}$, $\Sigma_{SFR}$, and O32. These three variables alone can predict $f_{esc}$ with an rms of 0.46 dex in the $z \sim 0.3$ sample. (See Section 3.4.)

5. We find that $f_{esc}$ may depend on different variables for high- versus low-mass galaxies. $\Sigma_{SFR}$ is statistically significant for the high-mass subsample but not the low-mass subsample and is the top-ranked variable for the high-mass subsample by forward and backward selection. In contrast, $M_{1500}$ and O32 are statistically significant for the low-mass subsample but not for the high-mass subsample. These differences suggest that the physical causes of LyC escape and the feedback mechanisms involved could differ across the galaxy population, with ionization and radiative feedback potentially playing a greater role in low-mass galaxies. (See Section 3.5.)

6. The multivariate Cox models are less successful at accurately predicting alternative measures of LyC escape, such as $F_{\lambda_{LyC}}/F_{\lambda_{1100}}$ and the LyC luminosity. We hypothesize that accurately predicting these properties requires information about the unattenuated SED shape, which is not included in our set of variables. (See Section 4.1.)

7. We also create Cox models to predict $Ly\alpha$, specifically the $f_{esc, Ly\alpha}$ and $Ly\alpha$ luminosity. Unlike the LyC $f_{esc}$ predictions, H I absorption EW and $E(B-V)_{neb}$, rather than $E(B-V)_{UV}$, are two of the most important variables for predicting $f_{esc, Ly\alpha}$. For $L(Ly\alpha)$, top-ranked and statistically significant variables are the H I absorption EW, $E(B-V)_{neb}$, $M_{1500}$, EW(Hβ), and $\Sigma_{SFR}$, which are similar to results found for the LARS (Runnholm et al. 2020). (See Section 4.2.)
Our models highlight the complexity of LyC escape and the dependence of $f_{esc}$ on multiple variables. Multivariate Cox models represent a new tool to investigate the relationships between $f_{esc}$ and galaxy properties and to predict $f_{esc}$ in the epoch of reionization. These models offer a substantial improvement over $f_{esc}$ predictions based on a single variable, and the models presented here can be refined and expanded to incorporate future observations. Our results show which variables have the greatest effect on $f_{esc}$ predictions for the LzLCS+ sample and for galaxies in different mass ranges. By exploring these same variables, future simulations could try to reproduce these scalings and determine their physical origin. These flexible multivariate models also offer an empirical approach to predicting $f_{esc}$ at high redshift. In a follow-up paper (Jaskot et al. 2024), we apply Cox models derived from the LzLCS+ to $z \sim 3$ and $z > 6$ samples and provide the best-fit models to the community for application to future observations.<fnref n="20" />

Acknowledgments

We thank Kate Follette for interesting discussions that led to this project, and we thank the anonymous referee for comments that improved the clarity of this work. A.E.J. and S.R.F. acknowledge support from NASA through grant HST-GO-15626. STScI is operated by AURA under NASA contract NAS-5-26555. A.S.-L. acknowledges support from the Knut and Alice Wallenberg Foundation. M.T. acknowledges support from the NWO grant 0.16.VIDI.189.162 ("ODIN").

Appendix

Tables 8–11 list the top four ranked variables by forward selection, backward selection, and an MC rank ordering for the $F_{\lambda\text{LyC}}/F_{\lambda\text{1100}}$, $L(\text{LyC})$, $f_{\text{esc,Ly}\alpha}$, and $L(\text{Ly}\alpha)$ predictions. We describe these ranking methods in Section 2.4 and discuss the results in Sections 4.1 and 4.2. The “Full Model” uses $\beta_{1550}$ as

<table><caption>**Table 8**
Most Important Variables for $F_{\lambda\text{LyC}}/F_{\lambda\text{1100}}$ Predictions</caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top-ranked Variables</td></tr><tr><td>1</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td></tr><tr><td>2</td><td>R<sub>λ</sub>(H I,abs)</td><td>M<sub>1500</sub></td><td>M<sub>1500</sub></td><td>M<sub>1500</sub></td><td>β<sub>1550</sub></td><td>log<sub>10</sub>(O<sub>32</sub>)</td></tr><tr><td>3</td><td>M<sub>1500</sub></td><td>EW(Lyα)</td><td>f<sub>esc,Lyα</sub></td><td>f<sub>esc,Lyα</sub></td><td>12+log<sub>10</sub>(O/H)</td><td>log<sub>10</sub>(EW(Hβ))</td></tr><tr><td>4</td><td>β<sub>1550</sub></td><td>β<sub>1550</sub></td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>E(B-V)<sub>neb</sub></td><td>log<sub>10</sub>(EW(Hβ))</td><td>β<sub>1550</sub></td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top Variables by MC Rank Order</td></tr><tr><td></td><td>EW(H I,abs) [1.01]</td><td>EW(H I,abs) [2.87]</td><td>EW(H I,abs) [1.00]</td><td>EW(H I,abs) [1.94]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.13]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.09]</td></tr><tr><td></td><td>R<sub>λ</sub>(H I,abs) [2.86]</td><td>R<sub>λ</sub>(H I,abs) [5.83]</td><td>f<sub>esc,Lyα</sub> [3.57]</td><td>f<sub>esc,Lyα</sub> [3.64]</td><td>log<sub>10</sub>(O<sub>32</sub>) [3.22]</td><td>log<sub>10</sub>(O<sub>32</sub>) [2.14]</td></tr><tr><td></td><td>f<sub>esc,Lyα</sub> [7.20]</td><td>log<sub>10</sub>(EW(Hβ)) [5.95]</td><td>M<sub>1500</sub> [3.79]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.58]</td><td>log<sub>10</sub>(EW(Hβ)) [5.45]</td><td>log<sub>10</sub>(EW(Hβ)) [3.11]</td></tr><tr><td></td><td>EW(Lyα) or log<sub>10</sub>(Σ<sub>SFR</sub>) [7.50]</td><td>EW(Lyα) [6.68]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.21]</td><td>log<sub>10</sub>(O<sub>32</sub>) [4.85]</td><td>β<sub>1550</sub> [3.97]</td><td>β<sub>1550</sub> [4.73]</td></tr></tbody></table>

Note. Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

<table><caption>**Table 9**
Most Important Variables for $L(\text{LyC})$ Predictions</caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top-ranked Variables</td></tr><tr><td>1</td><td>L(Lyα)</td><td>EW(Lyα)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>EW(H I,abs)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(O<sub>32</sub>)</td></tr><tr><td>2</td><td>E(B-V)<sub>neb</sub></td><td>M<sub>1500</sub></td><td>f<sub>esc,Lyα</sub></td><td>M<sub>1500</sub></td><td>β<sub>1550</sub></td><td>M<sub>1500</sub></td></tr><tr><td>3</td><td>M<sub>1500</sub></td><td>β<sub>1550</sub></td><td>M<sub>1500</sub></td><td>log<sub>10</sub>(O<sub>32</sub>)</td><td>M<sub>1500</sub></td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td></tr><tr><td>4</td><td>EW(Lyα)</td><td>EW(H I,abs)</td><td>log<sub>10</sub>(O<sub>32</sub>)</td><td>f<sub>esc,Lyα</sub></td><td>log<sub>10</sub>(O<sub>32</sub>)</td><td>log<sub>10</sub>(EW(Hβ))</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Top Variables by MC Rank Order</td></tr><tr><td></td><td>L(Lyα) [1.00]</td><td>M<sub>1500</sub> [3.05]</td><td>M<sub>1500</sub> [2.61]</td><td>M<sub>1500</sub> [2.95]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.17]</td><td>log<sub>10</sub>(O<sub>32</sub>) [1.71]</td></tr><tr><td></td><td>M<sub>1500</sub> [4.43]</td><td>EW(Lyα) [4.78]</td><td>EW(H I,abs) [3.28]</td><td>EW(H I,abs) [3.39]</td><td>β<sub>1550</sub> [3.40]</td><td>M<sub>1500</sub> [2.55]</td></tr><tr><td></td><td>E(B-V)<sub>neb</sub> [6.21]</td><td>EW(H I,abs) [6.08]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [3.30]</td><td>f<sub>esc,Lyα</sub> [4.27]</td><td>M<sub>1500</sub> [3.53]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [2.87]</td></tr><tr><td></td><td>EW(Lyα) [6.96]</td><td>log<sub>10</sub>(O<sub>32</sub>) [6.52]</td><td>f<sub>esc,Lyα</sub> [3.32]</td><td>log<sub>10</sub>(O<sub>32</sub>) [4.54]</td><td>log<sub>10</sub>(O<sub>32</sub>) [4.38]</td><td>log<sub>10</sub>(EW(Hβ)) [4.58]</td></tr></tbody></table>

Note. Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

<footnote n="20">²⁰ We have also created a Python script, available at https://github.com/sfury/LyCsurf, that generates Cox model $f_{esc}$ predictions using any desired combination of variables from the LzLCS+. The LzLCS+ data set used in this paper is available at the same page, and version 0.1.0 of the code is archived in Zenodo (Flury et al. 2024).</footnote>
<table><caption>**Table 10**
Most Important Variables for $f_{esc,Ly\alpha}$ Predictions<sup>a</sup></caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7">Top-ranked Variables</td></tr><tr><td>1</td><td>log<sub>10</sub>(M<sub>*</sub>)</td><td>E(B-V)<sub>neb</sub></td><td>log<sub>10</sub>(M<sub>*</sub>)</td><td>EW(H I,abs)</td><td>log<sub>10</sub>(M<sub>*</sub>)</td><td>β<sub>1550</sub></td></tr><tr><td>2</td><td>E(B-V)<sub>neb</sub></td><td>EW(LIS)</td><td>EW(H I,abs)</td><td>E(B-V)<sub>neb</sub></td><td>E(B-V)<sub>neb</sub></td><td>E(B-V)<sub>neb</sub></td></tr><tr><td>3</td><td>EW(H I,abs)</td><td>EW(H I,abs)</td><td>E(B-V)<sub>neb</sub></td><td>E(B-V)<sub>UV</sub></td><td>12+log<sub>10</sub>(O/H)</td><td>12+log<sub>10</sub>(O/H)</td></tr><tr><td>4</td><td>EW(LIS)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>E(B-V)<sub>UV</sub></td><td>12+log<sub>10</sub>(O/H)</td><td>β<sub>1550</sub></td><td>log<sub>10</sub>(M<sub>*</sub>)</td></tr><tr><td colspan="7">Top Variables by MC Rank Order</td></tr><tr><td></td><td>E(B-V)<sub>neb</sub> [1.61]</td><td>E(B-V)<sub>neb</sub> [2.01]</td><td>E(B-V)<sub>neb</sub> [2.36]</td><td>E(B-V)<sub>neb</sub> [2.43]</td><td>E(B-V)<sub>neb</sub> [2.18]</td><td>E(B-V)<sub>neb</sub> [1.93]</td></tr><tr><td></td><td>EW(H I,abs) [3.19]</td><td>EW(H I,abs) [2.85]</td><td>EW(H I,abs) [2.89]</td><td>EW(H I,abs) [2.51]</td><td>12+log<sub>10</sub>(O/H) [3.34]</td><td>12+log<sub>10</sub>(O/H) [2.91]</td></tr><tr><td></td><td>EW(LIS) [5.65]</td><td>12+log<sub>10</sub>(O/H) [5.87]</td><td>E(B-V)<sub>UV</sub> [3.92]</td><td>12+log<sub>10</sub>(O/H) [4.21]</td><td>log<sub>10</sub>(O32) [3.40]</td><td>β<sub>1550</sub> [3.51]</td></tr><tr><td></td><td>12+log<sub>10</sub>(O/H) [6.24]</td><td>EW(LIS) [6.39]</td><td>log<sub>10</sub>(O32) [3.97]</td><td>E(B-V)<sub>UV</sub> [4.29]</td><td>β<sub>1550</sub> [3.45]</td><td>log<sub>10</sub>(M<sub>*</sub>) [4.61]</td></tr></tbody></table>

**Notes.** Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

<sup>a</sup> Variable lists exclude any measures of Lyα.

<table><caption>**Table 11**
Most Important Variables for $L(Ly\alpha)$ Predictions<sup>a</sup></caption><thead><tr><th rowspan="2">Rank</th><th colspan="2">Full Model</th><th colspan="2">Fiducial+HI Model</th><th colspan="2">JWST Model</th></tr><tr><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th><th>Forward</th><th>Backward</th></tr></thead><tbody><tr><td colspan="7">Top-ranked Variables</td></tr><tr><td>1</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>M<sub>1500</sub></td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(EW(Hβ))</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td><td>log<sub>10</sub>(EW(Hβ))</td></tr><tr><td>2</td><td>E(B-V)<sub>neb</sub></td><td>log<sub>10</sub>(EW(Hβ))</td><td>E(B-V)<sub>neb</sub></td><td>M<sub>1500</sub></td><td>E(B-V)<sub>neb</sub></td><td>M<sub>1500</sub></td></tr><tr><td>3</td><td>M<sub>1500</sub></td><td>EW(H I,abs)</td><td>M<sub>1500</sub></td><td>EW(H I,abs)</td><td>M<sub>1500</sub></td><td>12+log<sub>10</sub>(O/H)</td></tr><tr><td>4</td><td>12+log<sub>10</sub>(O/H)</td><td>E(B-V)<sub>neb</sub></td><td>12+log<sub>10</sub>(O/H)</td><td>E(B-V)<sub>UV</sub></td><td>12+log<sub>10</sub>(O/H)</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>)</td></tr><tr><td colspan="7">Top Variables by MC Rank Order</td></tr><tr><td></td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.73]</td><td>M<sub>1500</sub> [1.93]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [2.10]</td><td>M<sub>1500</sub> [1.98]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [1.00]</td><td>M<sub>1500</sub> [1.94]</td></tr><tr><td></td><td>E(B-V)<sub>neb</sub> [3.24]</td><td>log<sub>10</sub>(EW(Hβ)) [4.25]</td><td>E(B-V)<sub>neb</sub> [2.71]</td><td>EW(H I,abs) [4.15]</td><td>E(B-V)<sub>neb</sub> [2.07]</td><td>log<sub>10</sub>(EW(Hβ)) [3.71]</td></tr><tr><td></td><td>M<sub>1500</sub> [3.50]</td><td>EW(H I,abs) [4.88]</td><td>M<sub>1500</sub> [3.24]</td><td>log<sub>10</sub>(EW(Hβ)) [4.21]</td><td>M<sub>1500</sub> [3.61]</td><td>log<sub>10</sub>(Σ<sub>SFR</sub>) [4.36]</td></tr><tr><td></td><td>12+log<sub>10</sub>(O/H) [5.08]</td><td>12+log<sub>10</sub>(O/H) [5.63]</td><td>EW(H I,abs) [4.78]</td><td>12+log<sub>10</sub>(O/H) [5.18]</td><td>12+log<sub>10</sub>(O/H) [3.97]</td><td>12+log<sub>10</sub>(O/H) [4.58]</td></tr></tbody></table>

**Notes.** Numbers in square brackets indicate the mean of the variable's ranks in the 100 MC runs.

<sup>a</sup> Variable lists exclude any measures of Lyα.

a measure of UV dust attenuation and includes multiple variables that trace Lyα or UV absorption-line strength. The “Fiducial+HI Model” uses $E(B-V)_{UV}$ to quantify dust attenuation, $f_{esc,Ly\alpha}$ as a measure of Lyα, and EW(H I,abs) as the only absorption-line variable. The “JWST Model” excludes absorption-line and Lyα variables.

## ORCID iDs

Anne E. Jaskot https://orcid.org/0000-0002-6790-5125

Anna Plantinga https://orcid.org/0000-0001-6281-7727

Sophia R. Flury https://orcid.org/0000-0002-0159-2613

Matthew Hayes https://orcid.org/0000-0001-8587-218X

John Chisholm https://orcid.org/0000-0002-0302-2577

Timothy Heckman https://orcid.org/0000-0001-6670-6370

Laura Pentericci https://orcid.org/0000-0001-8940-6768

Daniel Schaerer https://orcid.org/0000-0001-7144-7182

Maxime Trebitsch https://orcid.org/0000-0002-6849-5375

Anne Verhamme https://orcid.org/0000-0002-2201-1865

Cody Carr https://orcid.org/0000-0003-4166-2855

Henry C. Ferguson https://orcid.org/0000-0001-7113-2738

Zhiyuan Ji https://orcid.org/0000-0001-7673-2257

Mauro Giavalisco https://orcid.org/0000-0002-7831-8751

Alaina Henry https://orcid.org/0000-0002-6586-4446

Rui Marques-Chaves https://orcid.org/0000-0001-8442-1846

Göran Östlin https://orcid.org/0000-0002-3005-1349

Alberto Saldana-Lopez https://orcid.org/0000-0001-8419-3062

Claudia Scarlata https://orcid.org/0000-0002-9136-8876

Gábor Worsack https://orcid.org/0000-0003-0960-3580

Xinfeng Xu https://orcid.org/0000-0002-9217-7051

## References

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

Amorín, R. O., Rodríguez-Henríquez, M., Fernández, V., et al. 2024, A&A, 682, L25

Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Natur, 626, 975

Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJL, 902, L39

Bassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2019, MNRAS, 483, 5223

Begley, R., Cullen, F., McLure, R. J., et al. 2022, MNRAS, 513, 3510
Bian, F., & Fan, X. 2020, MNRAS, 493, L65

Bik, A., Östlin, G., Menacho, V., et al. 2018, A&A, 619, A131

Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28

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

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

Bouwens, R. J., Smit, R., Labbé, I., et al. 2016, ApJ, 831, 176

Bradburn, M. J., Clark, T. G., Love, S. B., & Altman, D. G. 2003, Br. J. Cancer, 89, 431

Bremer, J., & Dayal, P. 2023, MNRAS, 524, 118

Breslow, N. E. 1972, J. R. Stat. Soc. B Stat. Method., 34, 216

Cain, C., D'Aloisio, A., Gangolli, N., & Becker, G. D. 2021, ApJL, 917, L37

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

Chisholm, J., Gazagnes, S., Schaerer, D., et al. 2018, A&A, 616, A30

Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182

Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104

Choustikov, N., Katz, H., Saxena, A., et al. 2024, MNRAS, 529, 3751

Clark, T. G., Bradburn, M. J., Love, S. B., & Altman, D. G. 2003, Br. J. Cancer, 89, 232

Clarke, C., & Oey, M. S. 2002, MNRAS, 337, 1299

Cox, D. R. 1972, J. R. Stat. Soc. B Stat. Method., 34, 187

Davidson-Pilon, C. 2019, JOSS, 4, 1317

De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355

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

Endsley, R., Stark, D. P., Chevillard, J., & Charlot, S. 2021, MNRAS, 500, 5229

Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312

Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614

Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192

Finkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36

Finkelstein, S. L., Ryan, R. E., Jr+, Papovich, C., et al. 2015, ApJ, 810, 71

Fitzpatrick, E. L. 1999, PASP, 111, 63

Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87

Flury, S., Jaskot, A., & Silveyra, A. 2024, LyCsurv, v0.1.0, Zenodo, doi:10.5281/zenodo.11392442

Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJ, 930, 126

Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022b, ApJS, 260, 1

Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJL, 949, L25

Gazagnes, S., Chisholm, J., Schaerer, D., Verhamme, A., & Izotov, Y. 2020, A&A, 639, A85

Gazagnes, S., Chisholm, J., Schaerer, D., et al. 2018, A&A, 616, A29

Gnedin, N. Y., Kravtsov, A. V., & Chen, H.-W. 2008, ApJ, 672, 765

Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651

Harrell, F. E., Jr+, Calififf, R. M., Pryor, D. B., Lee, K. L., & Rosati, R. A. 1982, JAMA, 247, 2543

Hayes, M., Östlin, G., Schaerer, D., et al. 2013, ApJL, 765, L27

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

Heckman, T. M., Sembach, K. R., Meurer, G. R., et al. 2001, ApJ, 558, 56

Hu, W., Martin, C. L., Gronke, M., et al. 2023, ApJ, 956, 39

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

Isobe, T., Feigelson, E. D., & Nelson, P. I. 1986, ApJ, 306, 490

Izotov, Y. I., Guseva, N. G., & Thuan, T. X. 2011, ApJ, 728, 161

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

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

Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2018a, MNRAS, 474, 4514

Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2020, MNRAS, 491, 468

Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018b, MNRAS, 478, 4851

Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2021, MNRAS, 503, 1734

Jaskot, A. E., Dowd, T., Oey, M. S., Scarlata, C., & McKinney, J. 2019, ApJ, 885, 96

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

Jaskot, A. E., & Ravindranath, S. 2016, ApJ, 833, 136

Jaskot, A. E., Silveyra, A. C., Plantinga, A., et al. 2024, arXiv:2406.10179

Jecmen, M. C., & Oey, M. S. 2023, ApJ, 958, 149

Johnson, B. D., Leja, J. L., Conroy, C., & Speagle, J. S. 2019, Prospector: Stellar population inference from spectra and SEDs, Astrophysics Source Code Library, ascl:1905.025

Kakiuchi, K., & Gronke, M. 2021, ApJ, 908, 30

Katz, H., Saxena, A., Rosdahl, J., et al. 2023, MNRAS, 518, 270

Keenan, R. P., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 848, 12

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kim, J.-h., Krumholz, M. R., Wise, J. H., et al. 2013, ApJ, 775, 109

Kim, K., Malhotra, S., Rhoads, J. E., Butler, N. R., & Yang, H. 2020, ApJ, 893, 134

Kim, K. J., Bayliss, M. B., Rigby, J. R., et al. 2023, ApJL, 955, L17

Kim, K. J., Malhotra, S., Rhoads, J. E., & Yang, H. 2021, ApJ, 914, 2

Kimm, T., Blaizot, J., Garel, T., et al. 2019, MNRAS, 486, 2215

Kostyuk, I., Ciardi, B., & Ferrara, A. 2023, arXiv:2308.01476

Lee, J. C., Kennicutt, R. C. F., José, G. S. J., Sakai, S., & Akiyama, S. 2007, ApJL, 671, L113

Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leitherer, C., Schaerer, D., Goldader, J., et al. 2011, Starburst99: Synthesis Models for Galaxies with Active Star Formation, Astrophysics Source Code Library, ascl:1104.003

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Levesque, E. M., & Richardson, M. L. A. 2014, ApJ, 780, 100

Lin, Y.-H., Scarlata, C., Williams, H., et al. 2024, MNRAS, 527, 4173

Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42

Ma, X., Kasen, D., Hopkins, P. F., et al. 2015, MNRAS, 453, 960

Ma, X., Quataert, E., Wetzel, A., et al. 2020, MNRAS, 498, 2001

Maji, M., Verhamme, A., Rosdahl, J., et al. 2022, A&A, 663, A66

Makan, K., Worseck, G., Davies, F. B., et al. 2021, ApJ, 912, 38

Marchi, F., Pentericci, L., Guaita, L., et al. 2018, A&A, 614, A11

Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2021, MNRAS, 507, 524

Martin, C., Barlow, T., Barnhart, W., et al. 2003, Proc. SPIE, 4854, 336

Martin, C. L., Dijkstra, M., Henry, A., et al. 2015, ApJ, 803, 6

Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155

Mascia, S., Pentericci, L., Calabrò, A., et al. 2024, A&A, 685, A3

Matthee, J., Naidu, R. P., Pezzulli, G., et al. 2022, MNRAS, 512, 5960

Mauerhofer, V., Verhamme, A., Blaizot, J., et al. 2021, A&A, 646, A80

McKinney, J., Jaskot, A., Oey, M. S., et al. 2019, ApJ, 874, 52

McLernon, D. J., Giardiello, D., Van Calster, B., et al. 2023, Ann. Intern. Med., 176, 105

Morishita, T., Stiavelli, M., Chary, R.-R., et al. 2024, ApJ, 963, 9

Naidu, R. P., Forrest, B., Oesch, P. A., Tran, K.-V. H., & Holden, B. P. 2018, MNRAS, 478, 791

Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2022, MNRAS, 510, 4582

Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109

Nakajima, K., Ellis, R. S., Robertson, B. E., Tang, M., & Stark, D. P. 2020, ApJ, 889, 161

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

Östlin, G., Hayes, M., Duval, F., et al. 2014, ApJ, 797, 11

Paardekooper, J.-P., Khochfar, S., & Dalla Vecchia, C. 2015, MNRAS, 451, 2544

Pilyugin, L. S., Thuan, T. X., & Vílchez, J. M. 2006, MNRAS, 367, 1139

Razoumov, A. O., & Sommer-Larsen, J. 2010, ApJ, 710, 1239

Reddy, N. A., Steidel, C. C., Pettini, M., Bogosavljević, M., & Shapley, A. E. 2016, ApJ, 828, 108

Ricci, F., Marchesi, S., Shankar, F., La Franca, F., & Civano, F. 2017, MNRAS, 465, 1915

Rivera-Thorsen, T. E., Dahle, H., Chisholm, J., et al. 2019, Sci, 366, 738

Rivera-Thorsen, T. E., Hayes, M., & Melinder, J. 2022, A&A, 666, A145

Rivera-Thorsen, T. E., Östlin, G., Hayes, M., & Puschnig, J. 2017, ApJ, 837, 29

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

Rósdahl, J., Blaizot, J., Katz, H., et al. 2022, MNRAS, 515, 2386

Rósdahl, J., Katz, H., Blaizot, J., et al. 2018, MNRAS, 479, 994

Runnholm, A., Hayes, M., Melinder, J., et al. 2020, ApJ, 892, 48

Rutkowski, M. J., Scarlata, C., Henry, A., et al. 2017, ApJL, 841, L27

Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2022, A&A, 663, A59

Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2023, MNRAS, 522, 6295

Saxena, A., Bunker, A. J., Jones, G. C., et al. 2024, A&A, 684, A84

Saxena, A., Pentericci, L., Ellis, R. S., et al. 2022, MNRAS, 511, 120

Schaerer, D., Izotov, Y. I., Worseck, G., et al. 2022a, A&A, 658, L11

Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022b, A&A, 665, L4

Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, ApJ, 795, 20

Seive, T., Chisholm, J., Leclercq, F., & Zeimann, G. 2022, MNRAS, 515, 5556

Sharma, M., Theuns, T., Frenk, C., et al. 2016, MNRAS, 458, L94

Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252

Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, ApJL, 728, L2

Stelidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41

Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657

Trainor, R. F., Strom, A. L., Steidel, C. C., et al. 2019, ApJ, 887, 85

Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224

Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41
Verhamme, A., Orlitová, I., Schäerer, D., et al. 2017, A&A, 597, A13

Wang, B., Heckman, T. M., Amorín, R., et al. 2021, ApJ, 916, 3

Wang, B., Heckman, T. M., Leitherer, C., et al. 2019, ApJ, 885, 57

Williams, H., Kelly, P. L., Chen, W., et al. 2023, Sci, 380, 416

Wise, J. H., & Cen, R. 2009, ApJ, 693, 984

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

Worseck, G., Prochaska, J. X., Hennawi, J. F., & McQuinn, M. 2016, ApJ, 825, 144

Wyithe, J. S. B., & Loeb, A. 2013, MNRAS, 428, 2741

Xu, X., Henry, A., Heckman, T., et al. 2023, ApJ, 943, 94

Yang, H., Malhotra, S., Gronke, M., et al. 2017, ApJ, 844, 171

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

Zastrow, J., Oey, M. S., Veilleux, S., & McDonald, M. 2013, ApJ, 779, 76