JWST is uncovering the properties of everincreasing numbers of galaxies at z gt 6 during the epoch of reionization. Connecting these observed populations to the process of reionization requires understanding how efficiently they produce Lyman continuum LyC photons and what fraction f SUBescSUB of these photons escape into the intergalactic medium. By applying the Cox proportional hazards model a survival analysis technique to the Lowredshift Lyman Continuum Survey LzLCS we develop new empirical multivariate predictions for f SUBescSUB. The models developed from the LzLCS reproduce the observed f SUBescSUB for z  3 samples which suggests that LyC emitters may share similar properties at low and high redshift. Our bestperforming models for the z  3 galaxies include information about dust attenuation ionization andor morphology. We then apply these models to z  6 galaxies. For large photometric samples we find a median predicted f SUBescSUB  0.0470.14. For smaller spectroscopic samples which may include stronger emissionline galaxies we find that 33 of the galaxies have f SUBescSUB gt 0.2 and we identify several candidate extreme leakers with f SUBescSUB  0.5. The current samples show no strong trend between predicted f SUBescSUB and UV magnitude but limited spectroscopic information makes this result uncertain. Multivariate predictions can give significantly different results from singlevariable predictions and the predicted f SUBescSUB for highredshift galaxies can differ significantly depending on whether star formation rate surface density or radius is used as a measure of galaxy morphology. We provide all parameters necessary to predict f SUBescSUB for additional samples of highredshift galaxies using these models. SUPSUP Based on observations made with the NASAESA 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 526555. These observations are associated with programs GO15626 GO13744 GO14635 GO15341 and GO15639.
['Astrostatistics', 'Reionization', 'High-redshift galaxies', 'Starburst galaxies', 'Interstellar medium', 'Ultraviolet astronomy', 'Radiative transfer', '1882', '1383', '734', '1570', '847', '1736', '1335', 'Astrophysics - Astrophysics of Galaxies']
Multivariate Predictors of Lyman Continuum Escape. II. Predicting Lyman Continuum Escape Fractions for Highredshift Galaxies
{'Multivariate Predictors of LyC Escape II: Predicting LyC Escape Fractions for High-Redshift Galaxies ∗': "Anne E. Jaskot, 1 Anneliese C. Silveyra, 1, 2 Anna Plantinga, 3 Sophia R. Flury, 4 Matthew Hayes, 5 John Chisholm, 6 Timothy Heckman, 7 Laura Pentericci, 8 Daniel Schaerer, 9 Maxime Trebitsch, 10 Anne Verhamme, 9, 11 Cody Carr, 12, 13 Henry C. Ferguson, 14 Zhiyuan Ji, 15 Mauro Giavalisco, 4 Alaina Henry, 14 Rui Marques-Chaves, 9 Goran Ostlin, 5 Alberto Saldana-Lopez, 5 Claudia Scarlata, 16 G'abor Worseck, 17 and Xinfeng Xu 18  \n1 Department of Astronomy, Williams College, Williamstown, MA 01267, USA  \n2 Department of Physics, University of Nevada, Reno, NV 89557, USA  \n3 Department of Mathematics & Statistics, Williams College, Williamstown, MA 01267, USA  \n4 Department of Astronomy, University of Massachusetts Amherst, Amherst, MA 01002, USA  \n5 Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden  \n6 Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA  \n7 Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  \n8 INAF - Osservatorio Astronomico di Roma, via Frascati 33, 00078, Monteporzio Catone, Italy  \n9 Observatoire de Gen'eve, Universit'e de Gen'eve, Chemin Pegasi 51, 1290 Versoix, Switzerland  \n10 Astronomy, Kapteyn Astronomical Institute, Landleven 12, 9747 AD Groningen, The Netherlands  \n11 Univ. Lyon, Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France  \n12 Center for Cosmology and Computational Astrophysics, Institute for Advanced Study in Physics, Zhejiang University, Hangzhou 310058, China  \n13 Institute of Astronomy, School of Physics, Zhejiang University, Hangzhou 310058, China  \n14 Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  \n15 Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  \n16 Minnesota Institute for Astrophysics, School of Physics and Astronomy, University of Minnesota, 316 Church St. SE, Minneapolis, MN 55455, USA  \n17 VDI/VDE Innovation+Technik, Berlin, Germany  \n18 Center for Interdisciplinary Exploration and Research in Astrophysics, Northwestern University, Evanston, IL 60201, USA", 'ABSTRACT': 'JWST is uncovering the properties of ever increasing numbers of galaxies at z > 6, during the epoch of reionization. Connecting these observed populations to the process of reionization requires understanding how efficiently they produce Lyman continuum (LyC) photons and what fraction ( f esc ) of these photons escape into the intergalactic medium. By applying the Cox proportional hazards model, a survival analysis technique, to the Low-redshift Lyman Continuum Survey (LzLCS), we develop new, empirical, multivariate predictions for f esc . The models developed from the LzLCS reproduce the observed f esc for z ∼ 3 samples, which suggests that LyC emitters may share similar properties at low and high redshift. Our best-performing models for the z ∼ 3 galaxies include information about dust attenuation, ionization, and/or morphology. We then apply these models to z ≳ 6 galaxies. For large photometric samples, we find a median predicted f esc =0.047-0.14. For smaller spectroscopic samples, which may include stronger emission line galaxies, we find that ≥ 33% of the galaxies have f esc > 0 . 2, and we identify several candidate extreme leakers with f esc ≥ 0 . 5. The current samples show no strong trend between predicted f esc and UV magnitude, but limited spectroscopic information makes this result uncertain. Multivariate predictions can give significantly different results from single variable predictions, and the predicted f esc for high-redshift galaxies can differ significantly depending on whether star formation rate surface density or radius is used as a measure of galaxy morphology. We provide all parameters necessary to predict f esc for additional samples of high-redshift galaxies using these models.  \n- ∗ 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.', '1. INTRODUCTION': "Reionization represents a fundamental transformation of the universe's hydrogen gas and dramatically illustrates the effect galaxies can have on their surroundings. The presence of Ly α absorption in quasar spectra suggests that reionization ended sometime near z ∼ 6 (e.g., Fan et al. 2006; McGreer et al. 2015; Robertson 2022). The electron scattering optical depth of the cosmic microwave background is consistent with this general picture, constraining the midpoint of reionization to z ∼ 7 . 8 (Planck Collaboration et al. 2020). However, the exact timeline of reionization remains uncertain, in part because of the inhomogeneous nature of the reionization process (e.g., Eilers et al. 2018; Jung et al. 2020; Becker et al. 2021).  \nThe recently launched James Webb Space Telescope (JWST) is rapidly expanding our knowledge of galaxy properties during the epoch of reionization. Surveys are detecting more bright galaxies than anticipated, challenging our understanding of early galaxy evolution (e.g., Castellano et al. 2022; Naidu et al. 2022; Finkelstein et al. 2023; Donnan et al. 2023; Harikane et al. 2023; McLeod et al. 2024). JWST images suggest that galaxies in the reionization era are morphologically compact (e.g., Robertson et al. 2023; Ormerod et al. 2024; Morishita et al. 2024), and spectroscopic observations are tracing the evolution of nebular ionization and metallicity to redshifts well above 7 (e.g., Schaerer et al. 2022; Sanders et al. 2023; Tang et al. 2023; Fujimoto et al. 2023; Curti et al. 2023; Backhaus et al. 2024). Measurements of nebular emission lines in early galaxies are also constraining their ionizing, Lyman continuum (LyC) photon production rate, ξ ion , a key quantity needed to understand reionization. Such observations suggest that faint galaxies with bursty star formation histories may have elevated ξ ion (e.g., Atek et al. 2024; Simmonds et al. 2024; Saxena et al. 2024).  \nSeveral complementary approaches can give insights as to which galaxy populations dominated the reionization process. A galaxy's contribution to reionization depends on its ξ ion and the fraction f esc of these LyC photons that escape into the intergalactic medium (IGM). Using constraints on galaxy number densities at z > 6, some studies explore different distributions of ξ ion and f esc among the galaxy population and seek to match the observed constraints on the reionization timeline. Depending on the model assumptions, either moderately bright or very faint galaxies may provide the majority of the ionizing photons (e.g., Finkelstein et al. 2019; Naidu et al. 2020).  \nCosmological hydrodynamical simulations can also model the progress of reionization over time and can ex-  \nplore how and why LyC photons escape simulated galaxies. Simulations consistently find that feedback plays an essential role in clearing out obscuring material (e.g., Wise & Cen 2009; Cen & Kimm 2015; Paardekooper et al. 2015; Trebitsch et al. 2017). LyC input may peak shortly following a burst of star formation after supernovae explode (e.g., Ma et al. 2015; Trebitsch et al. 2017), although other factors such as binary star evolution or an extended star formation history may affect this timeline (e.g., Ma et al. 2015; Barrow et al. 2020; Katz et al. 2023). The geometry and mechanism of LyC escape may also vary across the galaxy population (e.g., Katz et al. 2023; Bremer & Dayal 2023). Despite general agreement on the importance of feedback, however, simulation predictions have yet to be confirmed observationally and the galaxies that are the primary drivers of reionization have not yet been conclusively identified.  \nDirectly detecting LyC and constraining f esc observationally becomes difficult above z > 4 (e.g., Inoue et al. 2014) due to high IGM attenuation. Fortunately, LyC observations at lower redshifts can help test simulation predictions regarding which galaxy properties regulate LyC escape. Over the past decade, the number of LyC detections has grown rapidly, with dozens of LyC Emitters (LCEs) now known at z ∼ 2 -3 (e.g., Mostardi et al. 2015; Shapley et al. 2016; Vanzella et al. 2016, 2018; Bassett et al. 2019; Fletcher et al. 2019; RiveraThorsen et al. 2019; Ji et al. 2020; Saxena et al. 2022), multiple detections of LyC emission in stacked samples at z ∼ 2 -4 (e.g., Marchi et al. 2018; Steidel et al. 2018; Bian & Fan 2020; Nakajima et al. 2020), and more than 50 known LCEs at z < 0 . 5 (e.g., Leitet et al. 2011, 2013; Borthakur et al. 2014; Izotov et al. 2016b, 2018b, 2021; Leitherer et al. 2016; Wang et al. 2019; Flury et al. 2022a). At both low and intermediate redshift, LCEs appear deficient in absorbing gas and dust (e.g., Gazagnes et al. 2018; Chisholm et al. 2018; SaldanaLopez et al. 2022; Steidel et al. 2018; Ji et al. 2020), morphologically compact (e.g., Borthakur et al. 2014; Izotov et al. 2018b; Flury et al. 2022b; Vanzella et al. 2016; Marchi et al. 2018; Rivera-Thorsen et al. 2019), and bright in higher ionization emission lines such as [O iii ] λ 5007 (e.g., Izotov et al. 2018b; Flury et al. 2022b; Vanzella et al. 2016; Fletcher et al. 2019; Nakajima et al. 2020).  \nBuilding on these observational and theoretical efforts, several studies have proposed diagnostics to predict f esc at z > 6 based on observable properties. Using simulated galaxies from the SPHINX cosmological radiation hydrodynamics simulation, Choustikov et al. (2024) develop a method to predict f esc from a linear combination of observables, including the UV slope, E(B-V),  \nH β luminosity, UV magnitude, and nebular line ratios. Most other studies have taken an empirical approach, constructing f esc diagnostics based on the large observational sample of LCEs at z ∼ 0 . 3 (e.g., Verhamme et al. 2017; Wang et al. 2021; Flury et al. 2022b; Xu et al. 2023). The largest of these z ∼ 0 . 3 samples is the Low-redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022a), a set of 66 galaxies with LyC measurements from the Hubble Space Telescope (HST) and ancillary ultraviolet and optical data from HST and the Sloan Digital Sky Survey (SDSS; Blanton et al. 2017). Combined with archival datasets (Izotov et al. 2016a,b, 2018a,b, 2021; Wang et al. 2019), this z ∼ 0 . 3 sample, hereafter the LzLCS+, consists of 89 galaxies with measured LyC or stringent upper limits. Based on an analysis of the LzLCS+, Chisholm et al. (2022) find that the UV slope β 1550 shows one of the strongest correlations with f esc and propose that this single variable can serve as a predictor of f esc at high redshift. Other recent studies of the LzLCS+ consider f esc diagnostics that incorporate information from multiple variables. SaldanaLopez et al. (2022) generate an equation to predict f esc from E(B-V) and low-ionization UV absorption lines. By tracing the gas and dust that destroy LyC photons, these parameters closely track f esc . Employing this method at high redshift may be a challenge, however, as measuring weak absorption lines requires high signal-to-noise observations of galaxy continua (but see Saldana-Lopez et al. 2023). Emission lines and photometry offer a simpler, if less direct, means of predicting f esc . Mascia et al. (2023) propose a new f esc diagnostic using β 1550 , [O iii ] λ 5007/[O ii ] λ 3727=O32, and halflight radius, three of the variables that correlate strongly with f esc in the LzLCS+. By applying this diagnostic at high redshift using JWST observations from the GLASS and CEERS surveys, they find predicted f esc of ∼ 0 . 1 for galaxies at z > 6 (Mascia et al. 2023, 2024). With a similar combination of variables ( β 1550 , O32, and UV magnitude), Lin et al. (2024) develop a regression model for the probability of LyC escape. Their model suggests that f esc may be high ( ∼ 0 . 2) in brighter galaxies in the epoch of reionization (Lin et al. 2024). These studies all agree that plausible LCE candidates exist at high redshift, although they differ in their f esc prediction methods.  \nIn Jaskot et al. (2024), hereafter Paper I, we use the LzLCS+ to develop new empirical multivariate models for predicting f esc . Because the LzLCS+ contains both LyC detections and upper limits, we adopt the statistical techniques of survival analysis, which are suitable for such censored data. Specifically, we employ the semiparametric Cox proportional hazards model to generate  \nf esc predictions based on a desired set of input observables. We show that a model limited to observables accessible at z > 6 can reproduce the observed f esc in the LzLCS+ with a root-mean-square scatter of 0.46 dex. Of these observables, three variables (O32, β 1550 , and the star formation rate surface density) are statistically significant in the fit, and a model limited to these three input observables predicts f esc as well as the full model. The Cox model technique can be customized to include any combination of variables that are available for most of the LzLCS+ galaxies and hence offers a flexible tool for predicting f esc at high redshift.  \nIn this paper, we apply Cox models for f esc to samples of high-redshift galaxies. Following the techniques in Paper I, which we summarize in § 2, we generate new models for the sets of variables in published high-redshift samples. In § 3, we first test the models' performance at z ∼ 3 using samples with published LyC measurements. With published samples at z ≳ 6, we then generate f esc predictions for galaxies in the epoch of reionization in § 4. In § 5, we compare our models with alternative proposed f esc prediction methods from the literature, and we discuss the implications of our results for studies of reionization in § 6. We summarize our conclusions in § 7. In the Appendix, we provide parameters for all models as well as examples on how to apply these models to future samples. We adopt a cosmology of H 0 = 70 km s -1 Mpc -1 , Ω m = 0 . 3, and Ω Λ = 0 . 7.", '2.1. Sample: The Low-redshift Lyman Continuum Survey': "We derive our empirical predictions of f esc from the LzLCS+, a combined, homogeneously processed dataset consisting of the Low-redshift Lyman Continuum Survey (Flury et al. 2022a) and archival samples with HST Cosmic Origins Spectrograph (COS) LyC observations (Izotov et al. 2016a,b, 2018a,b, 2021; Wang et al. 2019). A full description of the data processing and UV and optical measurements appears in Flury et al. (2022a), Saldana-Lopez et al. (2022), and Paper I, but we summarize key points here.  \nThe LzLCS+ contains 89 galaxies at z ∼ 0 . 3, a distance where the LyC redshifts into a sensitive wavelength range for the COS detector. The LzLCS targets were selected on properties proposed to indicate LyC escape: high O32 ratios (O32 ≥ 3), high star formation rate surface densities (Σ SFR > 0 . 1 M ⊙ yr -1 kpc -2 ), and/or blue UV slopes ( β < -2). The LzLCS+ covers a wide range of parameter space, spanning ∼ 2 dex in O32, ∼ 2 dex in Σ SFR , and an observed (not dust-corrected) UV absolute magnitude range of M 1500 = -18 . 3 to -  \n21.5. In Paper I and this study, we exclude one galaxy from the LzLCS+ sample: J1333+6246 (Izotov et al. 2016b). This galaxy has visibly truncated emission lines in its SDSS spectrum and unphysical Balmer line ratios, which together suggest that its nebular line flux measurements may be inaccurate. For this work, our total sample therefore includes 88 z ∼ 0 . 3 galaxies, 49 of which have detected LyC.  \nTo measure the LyC, we use COS G140L observations, processed using the calcos pipeline (v3.3.9) and the FaintCOS software routines (Worseck et al. 2016; Makan et al. 2021). We measure the LyC in a 20 ˚ A-wide wavelength bin near rest-frame 900 ˚ A, while excluding any wavelengths above 1180 ˚ A in the observed frame because of telluric emission. We follow the definitions of LyC detections and upper limits from Flury et al. (2022a), where detections are observations with a probability < 0 . 02275 of originating from background counts and the upper limit for non-detections represents the the 84th percentile of the background count distribution. To correct for Milky Way attenuation, we adopt the Green et al. (2018) dust maps and Fitzpatrick (1999) attenuation law.  \nFlury et al. (2022a) investigate several different measures of f esc . Here, as in Paper I, we adopt the absolute f esc , with Starburst99 (Leitherer et al. 2011, 2014) spectral energy distribution (SED) fits to the UV continuum providing the estimate of the initial intrinsic LyC (Chisholm et al. 2019; Saldana-Lopez et al. 2022). Alternative estimates of the intrinsic LyC from H β neglect LyC photons absorbed by dust, assume an isotropic geometry, and require an assumed simplistic star formation history (e.g., Flury et al. 2022a). Another alternative measure of LyC escape is the F λ LyC /F λ 1100 flux ratio. As a ratio of two observed fluxes, this method requires no model assumptions. However, we found in Paper I that this quantity was more difficult to predict from easily accessible observables with the survival analysis models. Hence, for the f esc models in this paper, we proceed with the absolute f esc estimates from the UV SED fits for the LzLCS+.  \nOptical and UV observations supply information about numerous other properties of the LzLCS+ galaxies (see Flury et al. 2022a and Saldana-Lopez et al. (2022) for details). We obtain stellar mass ( M ∗ ) estimates from Prospector (Leja et al. 2017; Johnson et al. 2019) fits to the SDSS and Galaxy Evolution Explorer ( GALEX ; Martin et al. 2003) photometry (Flury et al. 2022a, Ji et al. in prep.). With multi-Gaussian fits to nebular lines, we measure nebular line fluxes and equivalent widths (EWs), and we estimate the nebular dust attenuation, E(B-V) neb , from Balmer line ratios.  \nWe calculate the oxygen abundance using the direct method and the pyneb package (Luridiana et al. 2015), adopting n e = 100 cm -3 and the estimated [O iii ] λ 4363 flux from the Pilyugin et al. (2006) 'ff-relation' in cases where the [S ii ] λλ 6716,6731 doublet or [O iii ] λ 4363 auroral line are undetected. We estimate star formation rates (SFRs) from the Kennicutt & Evans (2012) SFR calibration using the dust-corrected H β luminosities and Case B H α /H β ratio (Storey & Hummer 1995). The COS near-UV acquisition images allow us to measure the UV half-light radius r 50 , NUV , which we then use to calculate Σ SFR as  \nΣ SFR = SFR 2 πr 2 50 , NUV . (1)  \nIn addition to f esc , we derive several other parameters from the HST COS UV spectra. We obtain estimates of the dust attenuation E(B-V) from the UV spectrum Starburst99 SED fits; we label this parameter as E(BV) UV to distinguish it from the nebular dust attenuation E(B-V) neb derived from the Balmer lines. Although the G140L spectra do not extend to rest-frame 1500 ˚ A, we can estimate the 'observed' (non-extinction corrected) absolute magnitude at 1500 ˚ A ( M 1500 ) and the power law index slope at 1550 ˚ A ( β 1550 ) by extrapolating the SED fits to longer wavelengths (Saldana-Lopez et al. 2022). The inferred β 1550 values do match the observed values for the few galaxies with existing longer wavelength UV spectra (Chisholm et al. 2022). We also measure Ly α properties from the COS spectra. We derive Ly α fluxes and EWs by linearly fitting the continuum within 100 ˚ A of Ly α , excluding regions affected by nebular or stellar features, and integrating all emission above the continuum. We then derive the Ly α escape fraction ( f esc , Ly α ) using the dust-corrected H β flux and intrinsic Case B Ly α /H β ratio (Storey & Hummer 1995) appropriate for the galaxies' measured electron temperatures and densities. We note that the Ly α measurements represent the net sum of both underlying absorption along the line of sight and scattered Ly α emission within the aperture. Nine galaxies have detectable Ly α absorption troughs that overlap with the Si iii λ 1206 absorption feature. We increase our uncertainties to account for the change in flux from omitting wavelengths within 500 km s -1 of the Si iii line. The inclusion or exclusion of this region has only minor effects, changing the Ly α EW by < 3 ˚ A and with typical changes to f esc , Ly α of only 0.001.", '2.2. Multivariate Survival Analysis: The Cox Proportional Hazards Model': "In order to generate multivariate diagnostics for f esc , we need to incorporate information from both the f esc  \ndetections and the upper limits. Within the field of statistics, survival analysis techniques are appropriate for censored datasets that contain limits. One such survival analysis method is the Cox proportional hazards model (Cox 1972; see Clark et al. 2003 and Bradburn et al. 2003 for reviews and Feigelson & Nelson 1985 and Isobe et al. 1986 for examples in astronomy). Here we describe the basic form of the Cox model and its assumptions. We refer the reader to Paper I for a more thorough discussion of this model and its application to the LzLCS+. To implement the Cox model, we use the CoxPHFitter routine in the lifelines python package (Davidson-Pilon 2019).  \nThe Cox proportional hazards model predicts the probability of a particular f esc value given a set of input variables. Like many other survival analysis techniques, implementations of the Cox model typically assume the dataset contains measurements and lower limits (so-called 'right-censored' data). In contrast, our dependent variable data consist of f esc measurements and associated upper limits. Consequently, we transform our f esc values to the absorbed fraction of LyC, f abs = 1 f esc , for use in the Cox model. However, for ease of interpretation, we put the f abs results back into the form of f esc in all figures and in the tabulated results in the Appendix.  \nAs applied to the LzLCS+, the Cox proportional hazards regression model fits for the probability of a LyC detection in an infinitesimally small increment of f abs , for a given set of independent variables and assuming no detection at a lower value of f abs (higher value of f esc ). The Cox model assumes a particular functional form for this f esc probability, which is known as the 'hazard function':  \nh ( f abs | x ) = h 0 ( f abs ) exp[ n ∑ i =1 b i ( x i -¯ x i )] . (2)  \nIn this equation, b i are the best-fit coefficients for each input variable x i , and ¯ x i is the mean value of each input variable in the reference LzLCS+ dataset. The term h 0 ( f abs ) is the baseline hazard function, the probability of having f abs in the case where all input variables match their average values within the LZLCS+. The Cox model is semi-parametric, in that the dependence on the input variables has a fixed exponential functional form, but the baseline hazard function, h 0 , is estimated non-parametrically. Despite its fixed functional form for the input variable dependence, the Cox model allows the input variables to take any form. For instance, we could define an input variable x i as a measured value j , the logarithm of that measurement log 10 ( j ), the square of that measurement j 2 , or any other functional form of our  \nchoosing. In our case, we opt to use logarithmic forms of our input variables where possible (e.g, log 10 (O32), M 1500 , log 10 ( M ∗ )). In this logarithmic form, most variables have a similar order of magnitude and scale in the same manner, and an order of magnitude increase in a particular input variable simply translates into increasing the probability h ( f abs | x ) by a factor of e b i . In this paper, the only variables that we do not use in a logarithmic fashion are EW(Ly α ) and f esc , Ly α , because they range from negative (Ly α absorption) to positive (Ly α emission) within the LzLCS+.  \nAlthough Equation 2 gives the probability of f abs , we would like to predict the expected values of f abs and f esc and their associated uncertainty. We adopt the median of the probability distribution as this expected value (e.g., Bradburn et al. 2003; Davidson-Pilon 2019); the model predicts that f esc will be above this value 50% of the time and below it 50% of the time. To determine this median f abs mathematically, we first calculate the survival function, S ( f abs ), the probability that we do not detect LyC at f abs , detect < f abs (i.e., at f esc , detect > f esc ). In the Cox model, the survival function is  \nS ( f abs ) = exp[ -HF 0 ( f abs ) · ph( x )] , (3)  \nwhere HF 0 is the baseline cumulative hazard function  \nHF 0 ( f abs ) = ∫ f abs 0 h 0 ( f ) df (4)  \nand ph( x ) is the partial hazards function, which describes how the probability scales with the set of input variables x :  \nph( x ) = exp[ n ∑ i =1 b i ( x i -¯ x i )] . (5)  \n(e.g., Cox 1972; Bradburn et al. 2003; Davidson-Pilon 2019; McLernon et al. 2023). The median f abs , which corresponds to our predicted f esc value, is the f abs value where S ( f abs )= 0 . 5. In some cases, S ( f abs ) > 0 . 5, even for our largest tabulated value of f abs , which implies that f abs is closer to 1 than we can determine and f esc ∼ 0. In this circumstance, we infer that f esc is arbitrarily small and report a predicted f esc = 0.  \nThe lifelines CoxPhFitter returns the best-fit coefficients b i and the cumulative baseline hazard HF 0 ( f abs ) for each value of f abs corresponding to a LyC detection in the LzLCS+ (see Paper I for more details about the lifelines methodology). In the Appendix, we provide these best-fit parameters for the models in Paper I and the models in this paper, and we give examples of how to use these models to calculate the predicted f esc for a set of observed input variables. The coefficients and  \nhazards can be used to predict the expected f esc for any galaxy, as long as it has estimated values for each of the independent variables used in the model.", '2.2.1. Uncertainty in the Predicted f esc': 'Because the survival function represents a probability distribution, we can also use it to calculate the expected uncertainty in our predicted f esc estimate. The f abs values where S ( f abs ) = 0.159 and S ( f abs ) = 0.841 represent the range in the predicted f esc corresponding to the Normal theory 1σ uncertainty. This uncertainty range reflects the inherent scatter of the relationships between f esc and the input variables, where this scatter can come both from measurement uncertainties and from genuine variation among the galaxy population. In Paper I, we tested the effect of measurement uncertainty by performing a Monte Carlo (MC) resampling of each independent and dependent variable according to its observational uncertainty. We found that the distribution of predicted f esc from the resampled inputs was nearly always smaller than the 1σ bounds inferred from the survival function. In other words, the inherent scatter in the correlations, not the measurement error, is the dominant source of uncertainty in the predicted f esc , and the survival function bounds serve as a reasonable estimate of the uncertainty in the f esc predictions.', '2.2.2. Goodness-of-Fit Metrics': "As in Paper I, we evaluate the Cox models' performance using several complementary metrics. The concordance index, C , is particularly useful, as it is appropriate for censored data. The concordance index assesses whether the model f esc predictions correctly sort the dataset in the order of its observed f esc . To evaluate this sorting, the concordance index calculation compares each possible pair of data points. Concordant pairs are those where the galaxy with higher observed f esc also has higher predicted f esc , discordant pairs are the opposite, and tied pairs have identical predicted f esc . Some pairs with upper limits in f esc lead to ambiguous rankings and do not appear in the concordance index calculation. With all pairs evaluated, the concordance index is calculated as  \nC = n c +0 . 5 n t n c + n d + n t , (6)  \nwhere n c , n d , and n t are the number of concordant, discordant, and tied pairs, respectively. A value of C = 1 . 0 indicates a perfect rank ordering, 0.5 is perfectly random ordering, and 0 is perfect disagreement. The major advantage of C is that it includes both LyC detections and non-detections. However, it only assesses the relative order of the predicted f esc values, not their accuracy. To  \nmeasure the latter, we turn to alternative quantities. One such metric is the R 2 statistic  \nR 2 = 1 -∑ i ( y i -f i ) 2 ∑ i ( y i -¯ y ) 2 , (7)  \nwhere y i are the observed values of log 10 ( f esc ), ¯ y is their mean value, and f i are the predicted f esc values from the Cox model. As in Paper I and Maji et al. (2022), we also calculate a variant of R 2 , the adjusted R 2 , which accounts for the number of free parameters p in the model and number of data points n :  \nR 2 adj = 1 -(1 -R 2 ) n -1 n -p -1 . (8)  \nWhen increasing the number of variables, R 2 adj increases only if a variable improves the predictions more than expected by chance. Finally, we report the root-meansquare (RMS) dispersion  \nRMS = √ ∑ i ( y i -f i ) 2 n . (9)  \nAs explained in Paper I, we evaluate these three quantities ( R 2 , R 2 adj , and RMS) using log 10 ( f esc ). The scatter in the log 10 ( f esc ) predictions is relatively consistent across the full range of observed f esc , whereas the scatter in the linear f esc changes systematically across this range. We can only calculate the R 2 , R 2 adj , and RMS metrics for the LzLCS+ galaxies with both detected LyC and non-zero predicted f esc , since these metrics require observed values of log 10 ( f esc ) and finite predicted log 10 ( f esc ). Hence, these metrics only indicate the accuracy of the model predictions for LCEs. However, the most successful models according to the R 2 , R 2 adj , and/or RMS metrics also tend to have high C and vice versa (see Paper I).", '2.2.3. Choice of Input Variables': 'In generating Cox models from the LzLCS+, we consider a limited set of input variables. Although the dependent variable in the Cox model can contain measurement limits, the independent variables cannot. We therefore only use independent variables that have measurements for nearly all the LzLCS+ galaxies, which requires us to exclude fainter emission line measurements (e.g., [O i ] λ 6300 and [S ii ] λλ 6716,6731) or measurements not widely available for the full sample (e.g., Ly α velocity peak separation). In this paper, all 88 galaxies in the LzLCS+ have measurements for our chosen input variables with one exception. The non-leaker J1046+5827 does not have a reported Σ SFR or r 50 , NUV , and we do not include it in deriving Cox models that use these variables. The Cox model may also fail to  \nconverge if we include highly collinear variables, which trace nearly identical properties. Given this limitation, we choose to use only one measure of UV dust attenuation (E(B-V) UV or β 1550 ) per model, and we do not include both Σ SFR and r 50 , NUV in a single model.  \nWe explored a variety of variable combinations in Paper I, and our best model attained R 2 = 0 . 69, R 2 adj = 0 . 60, RMS = 0 . 31 dex, and C = 0 . 91. We provide the parameters for this model in the Appendix. However, this model included the average EW of Lyman-series absorption lines (EW(H i ,abs)) and f esc , Ly α , both of which will be affected by the IGM at z > 6. Models in Paper I without UV absorption lines or Ly α showed higher scatter but an overall ability to identify LCEs, with R 2 = 0 . 29 -0 . 40, R 2 adj = 0 . 14 -0 . 35, RMS = 0 . 44 -0 . 47 dex, and C = 0 . 83 (Paper I; see Appendix for these model parameters). In this paper, we are specifically concerned with variables that are easily observable at z > 6 or measured in large samples at z ∼ 3. Hence, we likewise omit UV absorption line measurements, and we avoid Ly α measurements for most of the models tailored to z > 6 galaxies.', '3. TESTING THE COX MODELS AT Z ∼ 3': 'We apply the Cox models developed on the LzLCS+ sample to galaxies at z ∼ 3 and z ≳ 6. We use the z ∼ 3 galaxies to test whether the models based on the lowredshift LzLCS+ galaxies can successfully predict f esc for galaxies at high redshift. We then apply the Cox models to z ≳ 6 samples to predict f esc for galaxies in the epoch of reionization.', '3.1. High-redshift Datasets and Models': "First, we compile samples of z ∼ 3 galaxies with reported global absolute LyC f esc and at least three of the input variables from Paper I, which include stellar mass M ∗ , M 1500 , nebular EWs, metallicity, optical nebular line ratios, Σ SFR , r 50 , NUV , E(B-V) UV , β 1550 , and Ly α measurements. The z ∼ 3 LyC measurements include individual detections, reported LyC upper limits, and stacked samples that average over variations in IGM attenuation (Hainline et al. 2009; Vasei et al. 2016; James et al. 2018; de Barros et al. 2016; Vanzella et al. 2016; Steidel et al. 2018; Pahl et al. 2021; Bassett et al. 2019; Fletcher et al. 2019; Nakajima et al. 2020; Bian & Fan 2020; Marques-Chaves et al. 2021, 2022; Liu et al. 2023; Kerutt et al. 2024). We exclude one AGN from the Fletcher et al. (2019) sample and one shock-dominated galaxy from the Bassett et al. (2019) sample. Some wellknown z ∼ 3 LCEs do not appear in our sample, because they lack published global absolute f esc measurements (e.g., Ion3 and the Sunburst Arc). We also note that  \nsome of the high-redshift samples fall outside of the parameter space probed by the LzLCS+, such that our model predictions will extrapolate for these galaxies. The z ∼ 3 data has limitations as well; the compiled z ∼ 3 samples differ in the methods used for f esc and input variable measurements, and individual z ∼ 3 f esc measurements have high uncertainty due to unknown variations in IGM attenuation. We discuss the limitations of the z ∼ 3 comparison further in the following sections.  \nIn Paper I, we introduced a JWST Cox Model, which could be applied to an ideal z > 6 sample, with eight relevant variables from the LzLCS+ included. We also found that a 'limited JWST Model', fit using only the three top-ranked variables ( β 1550 , log 10 (O32), and log 10 (Σ SFR )) performed equally well for predicting f esc in the LzLCS+ sample. Unfortunately, many of the required variables for both the full and limited JWST models have not been measured for large samples of z ∼ 3 LCEs or z > 6 galaxies. Consequently, we run new models limited to the variables available for the z ∼ 3 and for z ≳ 6 samples. In Table 1 and Table 2, we list these models, the variables they include, and the samples to which they apply.  \nFor the z ≳ 6 samples, we prioritize models that can apply to large samples (e.g., Endsley et al. 2021, 2023; Morishita et al. 2024), models that can apply to faint galaxies (Atek et al. 2024), and models that have at least two of the most statistically significant variables in Paper I (a measure of dust attenuation plus a measurement of ionization or morphology). The 'TopThree' model includes only the three top-ranked variables from Paper I: dust, O32, and Σ SFR . In the TopThree model, we use E(B-V) UV rather than β 1550 , since it enables us to compare with a larger z ∼ 3 sample. As noted in Chisholm et al. (2022), the β 1550 and E(B-V) UV measurements for the LzLCS+ track each other almost perfectly and provide equivalent information. Most models in Table 1 use E(B-V) UV ; models with ' β ' in their name use β 1550 instead, and one model, LAE-O32-nodust, has no dust attenuation measurement. All models except the TopThree model include M 1500 , and most models include M ∗ as well. The only models without M ∗ are the LAE, LAE-O32-nodust, and ELG-O32β models.  \nWe provide the best-fit coefficients and cumulative baseline hazards for each model in the Appendix, which can be used to derive the median predicted f esc for a given galaxy and the uncertainty from the 16-84th percentiles of the f esc probability distribution. As with the LzLCS+ galaxies, we use an MC method to sample the variable uncertainties for the z ∼ 3 and z ≳ 6 observations and re-generate the predicted f esc . We again  \nfind that the uncertainty estimated from the Cox model survival function dominates over the uncertainty from sampling the input variables in nearly all cases.  \nTable 1. Cox Models for High-Redshift Predictions", 'T able 2 . Co x Mo del High-Redshift Samples': 'Note -(a) F rom de Barros et al. ( 2016 ); V anzella et al. ( 2016 , 2020 ). (b) F rom Hainline et al. ( 2009 ); V asei et al. ( 2016 ); James et al. ( 2018 ). W e adopt the a v erage dust-corrected O32 ratio o f the t w o similar-flux regio ns in Hainline et al. ( 2009 ). The E(B-V) UV and Σ SFR v alues are flux-w eigh ted a v erages of the regions in Jame s et al. ( 2018 ). (c) F rom Marques-Cha v es et al. ( 2022 ). (d) F rom Fletc her et al. ( 2019 ); Nak a jima et al. ( 2020 ). (e) W e adopt the f esc determined from the SED-fitting metho d in Liu et al. ( 2023 ). (f ) F rom Marques-Cha v es et al. ( 2021 ). (g) F rom Steidel et al. ( 2018 ); P ahl et al. ( 2021 ). (h) W e adopt M ∗ calculated using the burst y non-parametric star formation history , but include the alternativ e v alues in the unce rtain ties. W e use the direct metho d meta llicities for the t w o ga lax i es with secure [O i i i ] λ 4363 detections and the mo deled metallicities for the others. (i) W e correct M ∗ for magnification, usin g the v alues in Sc haerer et al. ( 2022 ). F or the source without a direct metho d metallicit y , w e adopt the strong -line metho d 12+log 10 (O/H). (j) W e exclude so urces flagged as ha ving significan t residuals in the morp hological fits. (k) Metallicities calculated using strong-line metho ds.  \nWe first test the Cox model predictions on z ∼ 3 samples to see whether the high-redshift LCEs behave similarly to the LzLCS+ sample. We derive Cox models for the variable sets in Table 1 using only the LzLCS+ data; the high-redshift samples are not included in the fitting process. We then use the models to predict f esc for each set of z ∼ 3 measurements, and we calculate goodness-of-fit metrics for the LzLCS+ sample alone, the high-redshift sample alone, and for the combined sample of low- and high-redshift galaxies. We list these goodness-of-fit metrics in Table 3.', '3.2. Model Performance': 'Figures 1 and 2 show the predicted vs. observed f esc for the LzLCS+ and z ∼ 3 samples. As seen from the plots and the goodness-of-fit statistics in Table 3, these models do not perform as well as the fiducial model from Paper I, but several (TopThree, LAE-O32, ELG-O32, ELG-O32β , ELG-O32β -Ly α , R50β ) are comparable to or better than the Paper I JWST model (see Tables 11-13 for the list of input variables used in the Paper I models). The metrics for the LzLCS+ sample are R 2 = 0 . 02 -0 . 42 (vs. 0.59 for the fiducial model and 0.29 for JWST), R 2 adj = -0 . 06 to 0.36 (vs. 0.48 for the fiducial model and 0.14 for JWST), RMS= 0 . 44 -0 . 59 (vs. 0.37 for the fiducial model and 0.47 for JWST), and C = 0 . 77 to 0.84 (vs. 0.88 for the fiducial model and 0.83 for JWST).  \nThe models developed on the LzLCS+ sample generally reproduce the observed f esc values of the z ∼ 3 LCEs. In fact, the high-redshift samples often have a lower RMS scatter than the predictions for the LzLCS+, and in several models, the R 2 and C values for the combined high- and low-redshift sample are comparable to or higher than the R 2 and C values for the LzLCS+ sample alone. Notably, predictions for the two strongest high-redshift LCEs, Ion2 and J1316+2614, match the observations nearly perfectly in most cases (Figures 1 and 2), demonstrating the success of the model in identifying extreme LCEs.  \nThe goodness-of-fit metrics give clues as to which variables are most important in predicting f esc . As seen in Table 3, the models with the highest R 2 and lowest RMS for the LzLCS+ sample are the TopThree, LAE-O32, ELG-O32 and variants, and the R50β models. These same models are the only ones that include both a measure of UV dust attenuation and either O32 or radius as variables. Turning to the C metric, which includes nondetections, the best-performing models for the LzLCS+ sample are the R50β , TopThree, LAE-O32, ELG-O32-  \nβ -Ly α , and LAE models. To properly estimate f esc for weak and non-LCEs, UV dust attenuation again appears important, but each of the best-performing models also includes either morphological information or Ly α measurements. These additional variables may help distinguish non-LCEs from LCEs.  \nBy all metrics, the LAE-O32-nodust and β -Metals models perform the worst for the LzLCS+ galaxies. The LAE-O32-nodust model is the only one that lacks E(BV) UV or β 1550 , which highlights the crucial role of dust extinction in LyC escape. Conversely, the β -Metals model relies almost exclusively on dust attenuation to infer f esc and lacks information on O32, Ly α , or morphology, which may better constrain the LyC absorption due to H i .', '3.2.1. Model Performance for z ∼ 3 Samples': "Comparing the models' performance for the highredshift and combined low- and high-redshift samples is difficult, because each high-redshift model applies to a different set of galaxies. For example, the ELG-O32β model has the lowest RMS for the high-redshift galaxies, but this RMS is based on a single galaxy: Ion2. Ion2 is the only galaxy that is included in most of the models (all but ELG-O32β -Ly α , where it has an input limit). For Ion2, we adopt an observed absolute f esc =0.75, in the middle of the 0.5-1 bounds reported by Vanzella et al. (2016). Using this f esc , Ion2's LyC escape is predicted most accurately by the ELG-O32β (RMS=0.02 dex) and ELG-O32 (RMS=0.04 dex) models. For Ion2, dust attenuation, luminosity, and ionization appear to be key factors in predicting its f esc , although most other models are still consistent with its reported range of f esc (RMS=0.09-0.15 dex for the next six models). The β -Metals model, which lacks information about O32 and morphology, is the only inaccurate model (RMS=0.43 dex). More surprisingly, however, the TopThree model gives the lowest predicted f esc =0.53 (RMS=0.15 dex), which highlights the role of UV magnitude in the f esc predictions, a point we discuss further in § 4.  \nNine galaxies, including four LCEs, appear in each of the LAE, LAE-O32, LAE-O32-nodust, and ELG-O32 model test samples. Although the LAE-O32-nodust model performed worst for the LzLCS+, for this common subsample of nine high-redshift galaxies, the LAEO32-nodust model has the lowest RMS (0.09 dex) and highest C (0.76), compared to RMS values of 0.27-0.32 dex and C values of 0.65-0.71 for these galaxies in the other three models. Given the generally successful performance of all four of these models, models with O32 and/or EW(Ly α ) appear able to predict f esc for highredshift galaxies equally well. Hence, while models with  \nTable 3. Metrics for High-Redshift Cox Modelsz ∼ 3 Sample", 'Combined Sample': 'Note -N Gal is the number of galaxies assessed in the fit and does not include any galaxies with limits for input variables. N Detect is the number of galaxies with LyC detections, finite f esc predictions, and no upper or lower limits for input variables; the R 2 and RMS metrics use only these galaxies. The R 2 statistic measures how well the predictive model explains the observed variance in the log 10 ( f esc ) data, with higher R 2 values corresponding to more accurate models. The R 2 adj metric accounts for the number of parameters used in the model and increases only if a variable improves the fit more than expected by chance. RMS is the rootmean-square dispersion of the predicted vs. observed log 10 ( f esc ) for the LyC detections. Higher C values indicate that the model more accurately sorts the observations in the correct order of increasing f esc . C includes both detections and galaxies with f esc upper limits. See § 2.2.2 for a full description of these statistics.  \nFigure 1. The f esc predictions from the TopThree (a), LAE (b), LAE-O32 (c), and LAE-O32-nodust (d) models. See Table 1 for model descriptions. Red circles represent LzLCS+ LyC detections, and blue triangles represent upper limits. We plot z ∼ 3 galaxies in green, with stars denoting LyC detections and triangles denoting upper limits. We identify the two strongest highredshift LCEs, Ion2 and J1316+2614, by a teal and a green diamond, respectively. Light green symbols indicate high-redshift galaxies that have a limit for one or more input variables. We draw circles around data points representing high-redshift galaxy stacks, which are less subject to uncertainty in IGM attenuation. The error bar in the upper left corner indicates the median size of the uncertainties in the observed and predicted f esc for the combined sample of low- and high-redshift galaxies. The dashed line shows a one-to-one correspondence. Several of the Cox models predict f esc for both low- and high-redshift galaxies with comparable accuracy.  \n<!-- image -->  \nELG-EW Model  \n<!-- image -->  \nELG-O32-  Model  \n<!-- image -->  \nR50-  Model  \n2  \n10  \nObserved  \nf  \nesc(LyC)  \n(e)  \n<!-- image -->  \nELG-O32 Model  \n<!-- image -->  \nFigure 2. The f esc predictions from the ELG-EW (a), ELG-O32 (b), ELG-O32β (c), ELG-O32β -Ly α (d), R50β (e), and β -Metals (f) models. See Table 1 for model descriptions. Ion2 is shown by an open symbol in panel d, because its prediction is based on a limit in f esc , Ly α . Other symbols are the same as in Figure 1.  \n<!-- image -->  \n10  \n)  \nC  \ny  \nL  \n(  \nc  \ns  \ne  \nf  \nd  \ne  \nt  \nc  \ni  \nd  \ne  \nr  \nP  \n10  \n10  \n10  \n10  \nesiduals  \nR  \n0  \n1  \n2  \n3  \n1  \n0  \n1  \n10  \n3  \n1  \n10  \n0  \nE(B-V) UV and O32 or morphology work best for the LzLCS+ sample, we cannot rule out the possibility that EW(Ly α ) may also be important in predicting f esc at z ∼ 3.  \nNevertheless, using EW(Ly α ) as the sole nebular measurement is not sufficient to accurately predict f esc . The LAE model (Figure 1b), which applies to the largest set of z ∼ 3 LyC measurements, is one of the least successful models at predicting f esc in the high-redshift galaxies. As seen in Table 3, this model has the largest RMS scatter for the z ∼ 3 sample (0.9 dex), lowest C (0.49), and the second worst R 2 (-4.2). The other model with comparably poor R 2 and C for z ∼ 3 galaxies is the ELG-EW model, which is the only other model that exclusively relies on an emission line EW for its input nebular information. Both models also have some of the lowest R 2 values for the LzLCS+ sample, which reflects the fact that these models substantially under-predict the f esc of the strongest LCEs in the LzLCS+. The LAE model shows this same tendency towards underprediction for the z ∼ 3 LCEs as well; 13 z ∼ 3 galaxies with observed f esc ≥ 0 . 2 are predicted to have f esc of only ≤ 0 . 03 by the model.  \nThe poor performance of the LAE model may result from limitations in both the model itself and in the z ∼ 3 data. The fact that both the LAE and ELG-EW models fare poorly for the LzLCS+ and for the high-redshift galaxies, especially for strong LCEs, suggests that nebular EW is a flawed tracer of f esc . Indeed, [O iii ] EW and Ly α EW should be low both for the weakest and for the strongest LCEs. In non-leakers, a high optical depth and corresponding lack of Ly α escape should lead to a low Ly α EW. Low [O iii ] and Ly α EWs could also indicate a weak current burst of star formation, without significant feedback or LyC production, which might result in undetected LyC. At lower optical depths and moderate f esc , Ly α EW may increase due to enhanced escape, and [O iii ] EWs may likewise increase if these lower optical depths are preferentially associated with stronger starbursts. However, the [O iii ] and Ly α EWs will decrease again in the strongest LCEs, as a general lack of nebular gas results in the limited production of nebular emission lines (e.g., Zackrisson et al. 2013; Nakajima & Ouchi 2014). Thus, EW may be an ambiguous indicator of LyC escape.  \nThe z ∼ 3 data also have limitations. The strong LCEs that are most severely under-predicted by the LAE model all come from the Liu et al. (2023) and Kerutt et al. (2024) samples of individual LAEs with LyC detections. Both papers note that their LyC detections represent the extreme of the population and may not be representative of the average galaxy population  \nwith these parameters, whose f esc may be much lower. In Paper I, we found a similar result for the strongest LCEs in the LzLCS+, whose only distinguishing feature compared to more moderate LCEs is their low line-ofsight H i content, suggestive of a favorable orientation. Another limitation of the high-redshift data is the uncertainty in the IGM attenuation. At z ∼ 3, this attenuation is significant and varies along different lines of sight (e.g., Rudie et al. 2013; Inoue et al. 2014; Vanzella et al. 2016; Steidel et al. 2018). This unknown sightline dependence leads to additional uncertainty in the reported f esc for individual galaxies and could cause some genuine LCEs to appear as non-detections or to have overestimated f esc .  \nThe stacked samples (Steidel et al. 2018; Nakajima et al. 2020; Bian & Fan 2020) shown circled in Figure 1, average over these variations and are less subject to this uncertainty. In the LAE model in Figure 1b, we see that the stacked samples do indeed show less scatter than the individual high-redshift detections. However, the model does not perfectly predict f esc for the stacked samples, and their scatter is still comparable to the scatter in the individual LzLCS+ galaxies. In addition to averaging over IGM attenuation variations, the stacks also average over any galaxy-to-galaxy variation in physical properties. Hence, the input variables may not represent the true set of properties of an individual system, leading to some uncertainty in the predicted f esc .  \nUnfortunately, the small sample sizes in most models and the lack of a common high-redshift sample across the models make it difficult to discern which model parameters are most important for predicting f esc at high redshift. Our results suggest that if high-redshift galaxies behave like their lower redshift counterparts in the LzLCS+ sample, E(B-V) UV or β 1550 , O32, and Σ SFR or half-light radius are essential variables to include. For the existing z ∼ 3 samples, models with at least two of these parameters do successfully predict the f esc of the available high-redshift galaxies to within ∼ 0.3 dex.  \nThis performance also shows that the LzLCS+ galaxies may indeed be reasonable analogs for high-redshift galaxies and that the same observable and physical properties may correlate with LyC escape at both low and high redshift (e.g., Saldana-Lopez et al. 2023; Schaerer et al. 2022; Mascia et al. 2023). This agreement is not trivial, as in principle, the z ∼ 0 . 3 and z ∼ 3 galaxies could have different star formation histories, dust properties, morphologies, or other properties, any of which could affect f esc . The LzLCS+ parameter space does cover the properties of the z ∼ 3 samples, with some exceptions. Seventeen z ∼ 3 targets (8 LCEs and 9 non-LCEs) are brighter than the LzLCS+ UV magni-  \nge by 0.1-3.18 magnitudes, including the strong LCEs Ion2 and J1316+2614. An additional two LCEs in the LAE model are fainter than the LzLCS+ by 0.50.7 magnitudes. The predictions for most models, all those with M 1500 as a variable, therefore extrapolate to an unobserved part of parameter space, yet still perform well.  \nThe data in Figures 1 and 2 and metrics in Table 3 include z ∼ 3 galaxies that have measured values for all of the input independent variables in the models 1 . However, additional galaxies in these samples have upper or lower limits for some independent variables, which can still provide potentially useful constraints on f esc . For these galaxies, we adopt the limit as the value of the independent variable and predict their f esc using the LzLCS+ models. We plot these approximate estimates in Figures 1 and 2 as light green symbols and an open symbol for Ion2, but we do not include these galaxies in the goodness-of-fit metrics.  \nThe galaxies with limits have one or more of the following: upper limits in M ∗ , lower limits in Σ SFR , lower limits in EW([O iii ]+H β ), lower limits in O32, lower limits in EW(Ly α ), lower limits in M 1500 , and lower limits in f esc , Ly α . One galaxy has an upper limit in EW([O iii ]+H β ) rather than a lower limit. Given the coefficients for these variables in the models, the mass and magnitude limits would cause the model to overestimate f esc , while the lower limits in the nebular emission lines should generally cause the model to underestimate f esc . For example, a galaxy with a lower limit of EW(Ly α ) > 100 ˚ Ais an even stronger Ly α emitter than we assume, and its predicted f esc will be an underestimate.  \nMost of the z ∼ 3 LCEs with input limits in the LAEO32 and LAE-O32-nodust models only have lower limits in O32, such that their predicted f esc should be an underestimate. For galaxies with lower limits in both a nebular property and M 1500 , we cannot easily interpret the predicted f esc as a lower or upper limit. However, the best-fit coefficients in the Appendix show that the models have a steeper dependence on the nebular lines than they do on M 1500 . Figures 1 and 2 show that for most of the galaxies with input limits, the measured limits still constrain f esc as accurately as for the rest of the sample. However, a few galaxies (90675, 101846, 105937 from Fletcher et al. 2019) have strongly underpredicted f esc , which suggests that their nebular lines may be much stronger than the reported limit.  \nFinally, we note that our model predictions do not account for systematic uncertainties, including differences in methodology. We adopt the published values for all data. However, each paper makes different assumptions, which could affect the tabulated values. For instance, papers differ in their adopted IGM transmission models, the models used in SED fits, and adopted dust attenuation laws. These assumptions could lead to systematic differences in properties like f esc , M ∗ , and E(B-V) UV among the different samples. In addition, LyC measurements may be photometric or spectroscopic and may cover different wavelength ranges. Despite these systematics, the models work well in predicting f esc . Adopting consistent methodologies could potentially result in better predictions; however, as shown by the scatter in the LzLCS+ sample, which has a consistent methodology, the inherent uncertainty in the model itself also limits the possible accuracy of predictions.', '3.2.2. Summary of z ∼ 3 Results': 'In conclusion, Cox models derived using the z ∼ 0 . 3 LzLCS+ sample can successfully predict f esc in z ∼ 3 LCEs. This agreement suggests that LCEs may have similar physical properties at both low and high redshift (e.g. Saxena et al. 2023). We find that the most accurate models include E(B-V) UV and O32 or morphology measurements as variables. However, larger high-redshift samples with a full suite of measurements are required to test this result. Future observations of z ∼ 3 LCEs with JWST will further clarify whether the relationship between f esc and physical properties evolves with redshift or remains constant. Preliminary JWST observations suggest that z ≳ 6 galaxies may indeed share numerous physical properties with low-redshift analog samples like the LzLCS (e.g., Schaerer et al. 2022; Endsley et al. 2023; Mascia et al. 2023; Lin et al. 2024). The models developed on the LzLCS+ sample and tested on z ∼ 3 LCEs may therefore apply equally well at the reionization epoch.', '4.1. Model Predictions for z ≳ 6 Samples': "Given the success of our models in predicting f esc at z ∼ 3, we now use the Cox models to estimate f esc for galaxies at z ≳ 6, where LyC is not detectable because of the IGM opacity. We consider several models: the 'ELG' models, R50β , and β -Metals. The β -Metals model is the most limited, as it only includes measurements of luminosity, mass, UV slope, and metallicity and has no information on nebular line strength or morphology. Based on the LzLCS+ galaxies, these models have an RMS scatter of 0.44-0.55 dex in f esc (see Sec-  \nand Table 3), with the ELG-O32 or R50β models giving the most accurate predictions, depending on the metric considered.  \nWe can apply these models to several samples of reionization-era galaxies, which span z ∼ 6 -14: photometric samples from Endsley et al. (2021), Endsley et al. (2023), Bouwens et al. (2023), and Morishita et al. (2024) and spectroscopic samples from Williams et al. (2023), Schaerer et al. (2022), Tang et al. (2023), Fujimoto et al. (2023), Saxena et al. (2023, 2024), Mascia et al. (2023), and Atek et al. (2024). The Endsley et al. (2021) sample is based on Spitzer , HST , and groundbased observations; all other samples incorporate JWST NIRCam photometry and several (Williams et al. 2023; Schaerer et al. 2022; Tang et al. 2023; Fujimoto et al. 2023; Saxena et al. 2023, 2024; Mascia et al. 2023; Atek et al. 2024) use JWST NIRSpec spectroscopy as well. We select all galaxies at z > 5 . 9 from these samples to investigate f esc in the epoch of reionization.", '4.1.1. Photometric Samples': 'We list the predicted f esc values and their associated uncertainties for the z ≳ 6 galaxies in Tables 4-8. The ELG-EW and R50β models (Tables 4 and 5) apply to the largest compilations of z ≳ 6 galaxies, with 183 measurements corresponding to 180 unique galaxies for the ELG-EW samples and 278 measurements at z > 5 . 9 for the R50β samples, although most of the galaxies for both models only have photometric redshifts. Most of the galaxies in the ELG-EW samples are from Endsley et al. (2021, 2023), and most of the R50β sample galaxies come from Morishita et al. (2024).  \nWeshow histograms of the ELG-EW and R50β model f esc predictions in Figure 3. The median of the f esc predictions for the ELG-EW model samples is fairly low, f esc = 0 . 047, which is lower than the average f esc value of 0.1-0.2 required to reionize the universe assuming a canonical ionizing photon production efficiency log 10 ( ξ ion ) ∼ 25 . 3(e.g., Finkelstein et al. 2015, 2019; Robertson et al. 2015; Naidu et al. 2020). Of these z ≳ 6 galaxies, 27% have f esc ≥ 0 . 1 and only 6% have f esc ≥ 0 . 2. However, as shown in Figure 2a, the ELGEW model is one of the less accurate models, with a greater tendency to underpredict the true f esc of the LzLCS+ sample. This underprediction implies that the ELG-EW variable set, which includes EW([O iii ]+H β ) but not O32, does not distinguish strong LCEs from weaker LCEs. The more accurate R50β model (see Table 3 and Figure 2) implies a higher fraction of LCEs, with a median value of f esc = 0 . 14, and 56% and 39% of the galaxies with f esc > 0 . 1 and f esc > 0 . 2, respectively. Both models find a substantial fraction of weak  \nor non-leakers, and 25% of the galaxies in both models have f esc < 0 . 03 -0 . 04. Taken at face value, our preliminary results would suggest that f esc values > 0 . 1 are common but not ubiquitous among moderate to bright galaxies ( M 1500 < -18) in the epoch of reionization (see also Mascia et al. 2023, 2024).  \nWe caution that some of the z ≳ 6 galaxies in the ELG-EW sample and many in the R50β sample fall outside of the parameter space probed by the LzLCS+, although a similar extrapolation did not seem to adversely affect the f esc predictions for the z ∼ 3 galaxies (Section 3). Of the 183 measurements in the ELG-EW model sample, 14 galaxies have brighter UV magnitudes than the LzLCS+ galaxies, one is fainter, and 17 have higher EW([O iii ]+H β ). Likewise, the R50β sample contains 8 galaxies brighter than and 30 galaxies fainter than the LzLCS+ sample. More concerningly, 114 of the 278 galaxies in the R50β compilation are more compact than the LzLCS+ galaxies, which could mean that their high inferred f esc values result from an incorrect extrapolation into this compact regime.', '4.1.2. Spectroscopic Samples': 'The ELG-O32 and ELG-O32β model predictions (Tables 6 and 7) apply to fewer z ≳ 6 galaxies, only 17 and 27 galaxies, respectively, but all of these galaxies have spectroscopic redshifts. Half the galaxies in the ELGO32 model have predicted f esc < 0 . 05, although most of these are only lower limits in f esc , and 25% of the ELGO32β sample galaxies have f esc ≤ 0 . 01. Like the R50β model predictions, low f esc is common, but the overall distribution also extends to very high f esc . Both models include some galaxies with lower limits in O32, whose predicted f esc values therefore also correspond to lower limits. Depending on how high the true f esc for these galaxies are, 35-82% of the ELG-O32 model galaxies and 33-41% of the ELG-O32β galaxies have f esc ≥ 0 . 2.  \nBecause these galaxies are all spectroscopically confirmed, these samples could be biased toward galaxies with stronger emission lines and hence higher f esc ; the median EW([O iii ]+H β ) of the spectroscopic ELG-O32 sample is 1790 ˚ A (Fujimoto et al. 2023; Tang et al. 2023; Saxena et al. 2023), compared with 690 ˚ A for the photometric samples in the ELG-EW model, for example (Endsley et al. 2021, 2023; Bouwens et al. 2023). In addition, these models may not underpredict f esc for strong LCEs to the same extent as the ELG-EW model (see Figure 2 and Table 3). For the six galaxies with pre-  \nTable 4. Predictions for z ≳ 6 Galaxies from the ELG-EW Model  \nNote -Predicted f esc values from the ELG-EW model for the z ≳ 6 galaxies from Endsley et al. (2021, 2023), Bouwens et al. (2023), Tang et al. (2023), Fujimoto et al. (2023), and Saxena et al. (2023). z phot and z spec are the photometric and spectroscopic redshifts. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions. For galaxies with upper limits on EW([O iii ]+H β ), the f esc predictions are also upper limits and are marked accordingly. The Reference column lists the publication used for the model input variables. The full, machine-readable version of this table is available online. We show the first five rows as an example here.  \nTable 5. Predictions for z ≳ 6 Galaxies from the R50β Model  \nNote -Predicted f esc values from the R50β model. z phot and z spec are the photometric and spectroscopic redshifts. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions. The Reference column lists the publication used for the model input variables. The full, machine-readable version of this table is available online. We show the first five rows as an example here.  \ndicted f esc > 0 . 2 in the ELG-O32 model, the ELG-EW predictions are indeed lower by ∆ f esc =0.26 on average, which suggests that model differences also account for some of the higher f esc values compared to the ELG-EW distribution.  \nWhile we have generated initial predictions for z ≳ 6 galaxies using our models, predictions for the largest samples are limited by a lack of spectroscopic information. We do find evidence of high f esc > 0 . 2 among smaller spectroscopic samples at z ≳ 6, but these samples may be biased toward stronger emission-line galaxies. To obtain more accurate predictions, high-redshift  \nstudies should prioritize observations of nebular line ratios and galaxy sizes for larger samples. As JWST surveys continue, larger, more representative spectroscopic samples will improve estimates of the average f esc at z > 6.', '4.2. Trends with UV Magnitude': "We examine the predicted f esc values as a function of magnitude in Figure 4. We plot predictions from the three models that apply to the largest sample sizes (ELG-EW, ELG-O32β , and R50β ) and from the model that applies to the faintest galaxies ( β -Metals; Table 8). We see no strong trend with magnitude, although the few brightest sources with M 1500 ≲ -21 . 5 do appear to have higher f esc on average, with a median f esc = 0 . 18 vs. 0.07 for fainter sources. Aside from these  \nTable 6. Predictions for z ≳ 6 Galaxies from the ELG-O32 Model  \nNote -Predicted f esc values from the ELG-O32 model. z spec is the spectroscopic redshift. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions. For galaxies with lower limits on O32, the f esc predictions are also lower limits and are marked accordingly. The Reference column lists the publication used for the model input variables.  \nrare sources, the ELG-EW, ELG-O32β , and β -Metals samples tend to show f esc within ∼ 0 . 02 -0 . 3 across a wide range of magnitudes, from -16 to -22. The R50β model (Figure 4b) also shows a flat distribution of f esc with magnitude, but it predicts many more galaxies at higher f esc . However, as previously noted, 41% of the galaxies in this model require extrapolating predictions to galaxies with extremely compact morphologies, which makes these predictions uncertain.  \nAny trends with magnitude or lack thereof in the other models should also be taken with caution. All the models have a tendency to under-predict f esc , but this tendency is most pronounced in the ELG-EW and β -Metals models (Figure 2), which trace the extremes of the magnitude range in Figure 4. The more reliable ELG-O32β model, which better reproduces the LzLCS+ f esc values, suggests that some of the galaxies in the M 1500 = -19 to -21 magnitude range may have significantly high f esc . The ELG-O32β sample does not include enough faint  \ngalaxies to reveal whether or not a trend between f esc and magnitude exists, however.  \nAt the faintest magnitudes, the Atek et al. (2024) sample, included in the β -Metals model includes several extremely low-luminosity lensed galaxies observed as part of the JWST UNCOVERsurvey (Bezanson et al. 2022; Weaver et al. 2024). For these galaxies, the β -Metals model therefore extrapolates the LzLCS+ sample to an unobserved parameter space. The only statistically significant predictor variable in the β -Metals model is β 1550 , and Figure 4 illustrates that these galaxies' blue slopes suggest at least moderate f esc (0.02-0.18; Table 8). Crucially, even these moderate f esc values, which may be under-predicted, are still greater than what is required for such faint galaxies to dominate reionization given their high LyC photon production rates (Atek et al. 2024).  \nWhile moderate or high f esc for faint, blue galaxies matches past predictions (e.g., Chisholm et al. 2022), the high f esc in the brighter sources in Figure 4 may seem surprising at first glance. In all four models, the fit-  \nTable 7. Predictions for z ≳ 6 Galaxies from the ELG-O32β Model  \nNote -Predicted f esc values from the ELG-O32β model. z spec is the spectroscopic redshift. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions. For galaxies with lower limits on O32, the f esc predictions are also lower limits and are marked accordingly. The Reference column lists the publication used for the model input variables.  \nted coefficients (see the Appendix) show that a brighter M 1500 results in a higher f esc at fixed values of the other inputs (fixed E(B-V), fixed O32, etc). This effect does not come from the association of higher observed UV luminosity with lower dust content. The data do not show such a correlation, and, because the models already include dust measurements, the model implies that brighter UV luminosities increase f esc at fixed dust attenuation. In addition, we see a similar relationship with M ∗ , where higher M ∗ is also associated with higher f esc at fixed values of the other parameters. Lin et al. (2024) find a similar result using a different technique. Using the LzLCS+ data, they fit a logistic regression  \nmodel for the probability of having LyC escape as a function of M 1500 , β 1550 , and O32 and likewise find a brighter M 1500 increases the probability of LyC escape for fixed values of the other variables. This scaling does not imply that higher UV luminosities correlate with f esc in the LzLCS+ sample; if anything, the LzLCS+ indicates a slight trend in the opposite direction (Flury et al. 2022b). Instead, this scaling implies that a bright galaxy with LCE-like properties, such as high O32 and blue UV slope, is a more extreme object than a faint galaxy with identical properties. For instance, at fixed O32, the LzLCS+ galaxies with higher luminosities tend to have higher f esc , but high O32 values are also uncom-  \nTable 8. Predictions for z ≳ 6 Galaxies from the β -Metals Model  \nNote -Predicted f esc values from the β -Metals model. z spec is the spectroscopic redshift. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions. The Reference column lists the publication used for the model input variables.  \nmon among the bright galaxies as a whole, leading to an overall trend of decreasing f esc with luminosity (e.g., Flury et al. 2022b).  \nThe enhanced f esc for bright LCE candidates affects the predictions in Figure 4 by boosting the f esc of a bright, dust-poor emission-line galaxy or compact galaxy relative to an otherwise similar faint galaxy. For instance, for the ELG-EW model, the brightest galaxies should have higher predicted f esc only if they have equally low E(B-V) and equally high EW as the fainter galaxies. Indeed, the ELG-EW samples, mostly galaxies from Endsley et al. (2021, 2023), show no trend between E(B-V) UV and M 1500 or between EW([O iii ]+H β ) and M 1500 ; with comparable dust and emission line strengths, the brighter galaxies therefore end up with higher predicted f esc . Endsley et al. (2023) note that lower metallicities or higher f esc among the fainter galaxies could suppress their EWs and account for the lack of an observed trend with M 1500 . Hence, other properties, such as O32, might be necessary to distinguish faint galaxies with high f esc . The EWs in the Endsley et al. (2021, 2023) samples also come from photometry, which may be more uncertain than spectroscopic measurements (e.g., Duan et al. 2024).  \nThe ELG-O32β samples have O32 measurements, but they similarly show no strong trend between β 1550 and M 1500 or between O32 and M 1500 . Likewise, the Morishita et al. (2024) sample, which constitutes the bulk of the R50β model sample, exhibits a relatively flat trend between radius and M 1500 . Only the β -Metals samples,  \nmostly from Atek et al. (2024), show any trend between physical properties and M 1500 , with the lowest metallicities and bluest slopes appearing among the faintest galaxies.  \nThe limitations of the models and the available z ≳ 6 datasets make it difficult to discern any trends between M 1500 and f esc . Limited measurements for galaxies at the brightest and faintest magnitudes restrict us to some of the less accurate Cox models. The models also require extrapolating to make predictions for galaxies outside the magnitude or radius range of the LzLCS+ sample. Better constraints on f esc within the z > 6 population will require larger spectroscopic samples across a wide range of galaxy magnitudes as well as better observational constraints on f esc in a more diverse set of galaxies at lower redshifts.", '4.3. Notes on Individual Sources': "4.3.1. Strong LCES  \nSeveral galaxies in the ELG-O32 and ELG-O32β models appear to be extremely strong LCEs (Tables 6 and 7). One unusual galaxy, CEERS-1019, has a predicted f esc = 1. Such a high f esc would seemingly conflict with the presence of strong emission lines, since it should correspond to a complete absence of absorbing gas. However, given the uncertainty in the model predictions, the f esc of CEERS-1019 could be as low as 0.6. Two other galaxies in the ELG-O32 model sample, CEERS-698 and CEERS-1027, have f esc ≥ 0 . 5, and f esc as high as 1 is consistent with their uncertainties. In-  \nFigure 3. The distribution of predicted f esc from the ELG-EW model (a) for z ≳ 6 galaxies from Endsley et al. (2021, 2023), Bouwens et al. (2023), Tang et al. (2023), Fujimoto et al. (2023), and Saxena et al. (2023) and from the R50β model (b) for z ≳ 6 galaxies from Morishita et al. (2024) and Mascia et al. (2023). The vertical dotted lines show the median and quartiles of the distributions ( f esc = 0 . 027 , 0 . 047 , 0 . 105 for ELG-EW and f esc = 0 . 039 , 0 . 141 , 0 . 467 for R50β ). Both models have some limitations but suggest that many, but not all, high-redshift galaxies may have f esc > 0 . 1.  \n<!-- image -->  \nFigure 4. Predicted f esc as a function of M 1500 for z ≳ 6 galaxies. Panel a shows the predicted f esc from the ELG-EW model (circles) for galaxies with photometric redshifts (gray) and spectroscopic redshifts (black), from the ELG-O32β model (blue stars), and from the β -Metals model (green diamonds). Panel b shows the predicted f esc from the R50β model (squares) with photometric (gray) and spectroscopic (black) redshifts. The crosses in the upper corners show the median uncertainties for each panel. Dark gray dotted lines connect the predictions for the same galaxy in different models. We see no strong trend between f esc and M 1500 , but many of the predictions suffer from a limited set of variables (ELG-EW and β -Metals models) or require extrapolation outside the LzLCS+ parameter space (R50β model).  \n<!-- image -->  \nPredicted  \nf  \nesc vs.  \nM  \n1500 for  \nz  \n6 Galaxies  \n18  \n19  \nM  \n(b)  \n20  \n1500  \nR50-  Model  \n21  \n22  \nc  \ns  \ne  \nf  \nd  \ne  \nt  \nc  \ni  \nd  \ne  \nr  \nP  \n10  \n10  \n10  \n10  \n0  \n1  \n2  \n3  \n17  \n23  \nterestingly, all three of these strong LCE candidates are strong Ly α emitters, which suggests that they reside within ionized bubbles (Tang et al. 2023). Similarly, two galaxies in the ELG-O32β samples have predicted f esc > 0 . 5: 18846 and 9422 from JADES (Saxena et al. 2024). With O32= 70 . 6, the ELG-O32β model infers f esc = 0 . 85 for 9422, much higher than the f esc = 0 . 01 value derived by Saxena et al. (2024) using the multivariate prediction method of Choustikov et al. (2024); we discuss the differences between our predictions and those of other models in § 5.  \nThe true f esc of some of these extreme galaxies may not be quite as high as the models predict, however. Because two of these strong LCEs, CEERS-1019 ( f esc = 1) and CEERS-698 ( f esc ∼ 0 . 7), have brighter UV magnitudes than the LzLCS+ sample by 0.4-0.9, their high predicted f esc could be unreliable. However, these same models successfully predicted the high f esc of the z ∼ 3 galaxies Ion2 and J1316+2614 (Figure 2), which are as bright as these z > 6 galaxies. The LzLCS+ sample itself also contains only three LCEs with f esc > 0 . 5, which limits the models' ability to accurately determine predictive relationships in the high f esc regime (Mascia et al. 2023). CEERS-1019 may differ from the LzLCS+ galaxies for additional reasons; it has unusually strong nitrogen emission and is a possible candidate for supermassive star formation (Marques-Chaves et al. 2024) or an AGN (Larson et al. 2023). The physical conditions occurring in extreme galaxies like CEERS-1019 likely are not present within the LzLCS+ sample. Empirical predictions based on z = 0 . 3 galaxies would not be suitable for such objects.  \nThe predictions for the other strong LCE candidates are likely more reliable, provided their star-formation and nebular conditions are not significantly different from the LCEs at z ∼ 3. Galaxy 9422 has a higher O32 ratio than the LzLCS+, but the other strong leaker candidates fall within the LzLCS+ parameter space, so that the model is not extrapolating to an unobserved regime. Hence, our models suggest that several strong emission-line galaxies identified in the epoch of reionization may be strong LCEs with f esc well above 0.2.", '4.3.2. JADES-GS-z7-LA': "Saxena et al. (2023) report the detection of a z = 7 . 3 galaxy with strong Ly α emission, JADES-GS-z7-LA, which is presumably located within an ionized bubble. With f esc , Ly α = 0 . 96, IGM absorption seems to have had little effect on the galaxy's Ly α emission. We therefore choose to incorporate the measured EW(Ly α ) and f esc , Ly α in our model predictions. We can apply most of the models in Table 1 to JADES-GS-z7-LA, omitting  \nonly the TopThree, R50β , and β -Metals models due to a lack of reported Σ SFR , radius, and 12+log 10 (O/H). We compare the different f esc predictions in Table 9. JADES-GS-z7-LA does fall outside of the LzLCS+ parameter space, which may limit the accuracy of the f esc predictions; it is fainter than the LzLCS+ galaxies by 1.6 mag, and it has a higher EW(Ly α ) by 144.1 ˚ A.  \nThe different Cox models disagree regarding the f esc of JADES-GS-z7-LA, with predicted f esc ranging from 0 to 0.40. The only nonzero predictions come from models that have information about UV dust attenuation and Ly α , but even the extreme Ly α properties of JADES-GS-z7-LA do not guarantee f esc > 0 . 1. The LAE-O32 and ELG-O32β -Ly α models include O32 and stellar mass as well as Ly α and find f esc = 0 . 017 -0 . 085. These two models are also among the best-performing models for the LzLCS+ sample (Table 3). JADES-GSz7-LA is not devoid of dust, with E(B-V) UV = 0 . 10, and its UV slope β = -2 . 1 is not extreme. These properties are consistent with the properties of z ∼ 0 . 3 galaxies with moderate f esc ; all the LzLCS+ galaxies with higher f esc > 0 . 1 have lower E(B-V) UV < 0 . 1, and a slope of β 1550 = -2 . 1 matches the median value for the moderate LCEs with f esc = 0 . 05 -0 . 1.  \nFurthermore, as discussed above, given JADES-GSz7-LA's low mass (log 10 ( M ∗ / M ⊙ )=7 . 15), its O32 value of 8.8 may not be extreme. Even its strong Ly α emission does not necessarily imply a low optical depth along the line of sight, since Ly α photons can scatter. For instance, the non-leaking z ∼ 0 . 3 galaxy J1248+4259 (Izotov et al. 2018b) has f esc ≤ 0 . 013 and EW(Ly α )= 256 ˚ A. One piece of evidence that favors a high f esc for JADESGS-z7-LA is the small offset between the Ly α velocity and the systemic redshift (Saxena et al. 2023). With the low resolution of most of the LzLCS+ UV spectra, we cannot currently include this parameter in our models, and it may boost the predicted f esc for JADES-GS-z7LA. However, this low velocity offset could still be consistent with little to no LyC escape; its offset of 120 ± 80 kms -1 also resembles the Ly α profile of the non-leaker J1248+4259, whose red peak is offset from the systemic velocity by < 200 kms -1 .  \nIn agreement with the assessment of Saxena et al. (2023), we thus find that the evidence for LyC escape in JADES-GS-z7-LA is ambiguous. The multivariate models that incorporate the most information suggest that it likely has a moderate f esc ∼ 0 . 017 -0 . 085 and may therefore not be solely responsible for producing its ionized bubble in the IGM (e.g., Witstok et al. 2024).  \nThe case of JADES-GS-z7-LA highlights the fact that the Cox multivariate predictions can differ from singlevariable estimates (e.g., based on Ly α alone). Galaxies with high O32 provide another example of this result. For instance, Williams et al. (2023) present JWST observations of a lensed low-mass (log 10 ( M ∗ / M ⊙ ) = 7.7) galaxy, 11027, at z = 9 . 51 and hypothesize that the galaxy has a high f esc > 0 . 1 based on its high O32 ratio of 12 ([O iii ] λλ 5007+4959/[O ii ]=16). However Lin et al. (2024) suggest that 11027 is not likely to be a strong LCE considering the combination of O32, UV magnitude, and UV slope. Our ELG-O32 and ELGO32β Cox models (Tables 6-7) likewise predict f esc of only 0.012-0.021. O32 alone is not sufficient to constrain f esc , and high O32 values are common among the lowest mass galaxies in the LzLCS+ sample (Flury et al. 2022b), including weak LyC emitters. In contrast, the ELG-O32β model predicts that several galaxies with lower O32 than galaxy 11027 (e.g., 6355, 10000, 10612, 12637 with O32=6.3-10.6, Schaerer et al. 2022; Mascia et al. 2023; Saxena et al. 2024) are actually more likely to be LCEs because of their higher luminosities and/or bluer UV slopes (Table 7). The Cox predictions emphasize the fact that not all galaxies with high O32 have high f esc (e.g., Izotov et al. 2018b; Flury et al. 2022b) and additional properties such as mass, UV luminosity, and dust extinction are important to consider.  \nMultivariate models are an important tool for predicting f esc in the epoch of reionization. Our results highlight the fact that multivariate models can give substantially different predictions from single variable estimates, and f esc predictions should therefore incorporate as much information as possible.", '5.1. Single-Variable Predictions': "Among the observable properties measured for the LzLCS+ sample, the UV slope β 1550 shows one of the strongest correlations with f esc . Consequently, Chisholm et al. (2022) propose that β 1550 could predict f esc for high-redshift galaxies and derive a relationship between β 1550 and f esc based on the LzLCS+ dataset. In Figure 5, we compare the predicted f esc from this method with the observed f esc for the LzLCS+ and for z ∼ 3 galaxies (Vanzella et al. 2016; Bassett et al. 2019; Marques-Chaves et al. 2022). We list the goodness-offit statistics for the LzLCS+ with the Chisholm et al. (2022) model in Table 10. The Chisholm et al. (2022) model's high R 2 = 0 . 45 and low RMS= 0 . 43 are comparable to the best-performing Cox models in Tables 3  \nand 10, which demonstrates that β 1550 alone can indeed predict f esc for LCEs reasonably well. Consistent with this result, Paper I finds that, of the variables accessible at high redshift, β 1550 is the most important variable to include in the Cox models.  \nDespite its success for LCEs, however, the Chisholm et al. (2022) model has more difficulty in constraining f esc in the LzLCS+ non-detections, as indicated by its lower concordance C = 0 . 76. This concordance is lower than all the Cox models in Tables 3 and 10, which have C = 0 . 77 -0 . 83. The difference between the predictions for detections and non-detections is apparent in Figure 5. For LzLCS+ galaxies with f esc detections or upper limits between 0 . 01 -0 . 1, the one-to-one relation between the Chisholm et al. (2022) model predictions and the observations runs right between the LzLCS+ detections. However, nearly all the non-detections are above this line, indicating that the model is systematically over-predicting their f esc . In contrast, for one of the better performing Cox models, such as the ELG-O32 model in Figure 2b, both detections and non-detections in this same f esc range fall on either side of the one-toone relation, indicating that predictions for detections and non-detections are comparable.  \nFigure 5 also applies the Chisholm et al. (2022) model to two z ∼ 3 LCEs: Ion2 and J1316+2614. The model correctly identifies both galaxies as LCEs but underpredicts their f esc , although Ion2's observed f esc is consistent with the stated uncertainties in the Chisholm et al. (2022) predictions. Comparing the Ion2 predictions with those from the Cox models discussed in § 3, we find that the Chisholm et al. (2022) single-variable β 1550 model is one of the least accurate at predicting Ion2's f esc , with a predicted f esc = 0 . 29, corresponding to RMS=0 . 40 dex. As discussed in Chisholm et al. (2022), β 1550 only tracks the loss of LyC photons due to dust. Scatter below and above the Chisholm et al. (2022) β 1550 -f esc relation results from variations in the absorbing H i column. Other properties, such as O32, Σ SFR , and UV magnitude, may better track the H i component of LyC absorption for galaxies like Ion2.  \nFigure 6 illustrates how including these other variables affects predictions at z > 4. We plot the difference between the observed LzLCS+ f esc and the Chisholm et al. (2022) predictions as a function of β 1550 . Galaxies are color-coded by their O32 ratio (Figure 6a) or Σ SFR (Figure 6b). For the reddest β 1550 , all the model predictions and observations agree that galaxies have little to no LyC escape. However, at blue UV slopes, the disagreement increases. The Chisholm et al. (2022) model assigns all galaxies a single value of f esc , whereas the LzLCS+ observations show a range of f esc and deviate  \nTable 9. f esc Predictions for JADES-GS-z7-LA  \nNote -Predicted f esc values for JADES-GS-z7LA from different models. f esc , min and f esc , max represent the 15.9 and 84.1 percentiles of the model f esc predictions.  \nfrom the predictions by as much as ∆ f esc = 0 . 2 lower or ∆ f esc = 0 . 6 higher at the bluest slopes. For the LzLCS+ galaxies, this difference is associated with other differences in the galaxies' properties; galaxies with higher O32 or higher Σ SFR tend to be offset to higher f esc , and galaxies with lower f esc than the Chisholm et al. (2022) predictions tend to have lower O32 or lower Σ SFR .  \nThe multivariate model predictions reflect these trends with other properties and assign galaxies a higher or lower f esc at fixed β 1550 accordingly. In Figure 6, we also plot the difference between the predicted f esc for z > 4 galaxies from multivariate models vs. the singlevariable Chisholm et al. (2022) model. All these multivariate models use β 1550 in addition to two or more other variables. The Mascia et al. (2023) model (squares) predicts f esc using a linear fit to β 1550 , log 10 (O32) and log 10 ( r 50 , NUV ) for their sample of z > 4 JWST-GLASS galaxies. We show predictions from two Cox models in crosses and stars: a Cox model using the three topranked, high-redshift accessible variables from Paper I ( β 1550 , O32, and Σ SFR ) for the Mascia et al. (2023) sample and the ELG-O32β Model ( §§ 3-4), which uses β 1550 , O32, and M 1500 as inputs to predict f esc for z > 4 galaxies from Williams et al. (2023), Schaerer et al. (2022), Mascia et al. (2023), and Saxena et al. (2024). Like the LzLCS+ observations, the multivariate z > 4 predictions can significantly deviate from the Chisholm et al. (2022) predictions, especially at blue UV slopes. Highredshift surveys relying only on β 1550 could therefore miss some of the physical differences among blue galaxies that might indicate variations in f esc . Properties such as O32 and Σ SFR are sensitive to feedback in galaxies and may complement β 1550 measurements by tracing a galaxy's ability to carve low column density H i channels.", '5.2. Multivariate Predictions': "Recently, Choustikov et al. (2024) and Mascia et al. (2023) have developed multivariate linear regression models to predict f esc from a set of observable variables. Choustikov et al. (2024) derive their model for f esc from synthetic spectra of z > 4 galaxies from the SPHINX cosmological radiation hydrodynamics simulation. Mascia et al. (2023) base their predictions on the LzLCS+ sample but adopt the upper limit in f esc as the observed value for non-detections. Here, we use galaxies at low and high redshift to compare the predictions of the Cox models with the predictions from these literature models.  \nIn Figure 7a, we show the results of applying the Choustikov et al. (2024) model to the LzLCS+ sample, and we list the corresponding goodness-of-fit statistics in Table 10. The Choustikov et al. (2024) model fails to predict the observed f esc in the LzLCS+ sample. This result does not depend on our method of calculating f esc ; we find a similarly poor fit if we substitute the f esc derived from H β , as used in the Choustikov et al. (2024) simulations. In Figure 7b, we show f esc predictions from a Cox model run using the same set of variables ( β 1550 , E(B-V) neb , L (H β ), M 1500 , R23=([O iii ] λλ 5007,4959+[O ii ] λ 3727)/H β , and O32) as in the Choustikov et al. (2024) linear regression mode, and we list the goodness-of-fit metrics in Table 10.  \nThis model performs comparably to the other Cox models in Table 3, which shows that the reason for the Choustikov et al. (2024) model's difficulty is not its set of variables. Rather, in the Choustikov et al. (2024) model, f esc anti-correlates with O32, whereas the LzLCS+ data  \nFigure 5. The f esc predictions for the LzLCS+ and z ∼ 3 galaxies from the Chisholm et al. (2022) model, which uses a single predictor variable: β 1550 . Symbols are the same as in Figure 1. The model reproduces the f esc of LCEs but over-predicts the f esc for the LzLCS+ non-detections and under-predicts f esc for the strongest LCEs.  \n<!-- image -->  \nTable 10. Goodness-of-Fit Statistics for Literature vs. Cox Models  \nNote -A comparison of the goodness-of-fit statistics for f esc predictions for the LzLCS+ sample from the Chisholm et al. (2022), Choustikov et al. (2024), and Mascia et al. (2023) literature models and various Cox models. We list statistics for the JWST Model from Paper I and a model limited to its three top-ranked variables. We also list statistics for Cox models run using the same sets of variables as in the Choustikov et al. (2024) and Mascia et al. (2023) models.  \n<!-- image -->  \nFigure 6. Circles (detections) and triangles (upper limits) show the difference between LzLCS+ observations and Chisholm et al. (2022) predicted f esc as a function of observed β 1550 . Other symbols show the difference between multivariate model predictions for z > 4 samples and the Chisholm et al. (2022) single-variable predictions as a function of β 1550 . Squares compare the f esc predictions from Mascia et al. (2023), using β 1550 , O32, and r 50 , NUV . Crosses show f esc predictions for the same z > 4 Mascia et al. (2023) sample for a Cox model using β 1550 , O32, and Σ SFR . Stars show the ELG-O32β Cox Model predictions for z > 4 galaxies from Williams et al. (2023), Schaerer et al. (2022), Mascia et al. (2023), and Saxena et al. (2024), where the input variables are β 1550 , O32, and M 1500 . We color-code galaxies by their observed O32 ratio in panel (a) and by Σ SFR in panel (b). At blue UV slopes, the LzLCS+ observations and high-redshift model predictions can differ significantly from the Chisholm et al. (2022) model predictions, depending on the galaxies' other properties, such as O32 and Σ SFR .  \n<!-- image -->  \nand the Cox model using the Choustikov et al. (2024) variable set indicate a correlation between f esc and O32. The compact, lower-mass LCEs in the LzLCS do not have counterparts in the SPHINX simulations; these compact LCEs often have high O32 ratios ≳ 10 (Flury et al. 2022b), whereas the SPHINX galaxies with O32 > 10 typically have lower values of f esc . Hence, while we find that high O32 increases f esc in our multivariate predictions, Choustikov et al. (2024) find the opposite, that at fixed β and E(B-V), galaxies with older ages and lower O32 have higher f esc . Radiative feedback with a turbulent gas structure may allow f esc from high O32 galaxies at earlier ages than predicted by cosmological simulations (e.g., Kakiichi & Gronke 2021; Kimm et al. 2019; Choustikov et al. 2024). However, further studies of z ∼ 3 galaxies are also necessary to test whether the LzLCS population and the LzLCS-derived Cox models correctly describe the high-redshift population.  \nLike our Cox model predictions in § 4, Choustikov et al. (2024) predict a low f esc = 0 . 03 for JADES-GS-z7-LA, due to its moderate UV slope and dust content. However, because of the different dependence of f esc on O32, the Choustikov et al. (2024) model predicts a dramatically lower f esc for CEERS-44, CEERS-698, CEERS1019, and CEERS-1027 ( f esc = 0 . 006 -0 . 1) and for JADES 9422 ( f esc = 0 . 01; Saxena et al. 2024) than our model predictions ( f esc > 0 . 4 to 1). Choustikov et al. (2024) point out that their predicted f esc values are consistent with f esc , Ly α > f esc as observed locally (e.g., Flury et al. 2022b), but IGM effects may complicate any comparison with f esc , Ly α . As noted above ( § 4), the Cox model predictions for some of these galaxies require extrapolating outside the LzLCS+ parameter space, and their f esc may not truly be as extreme as this model predicts. Still, these galaxies remain plausible candidates for high f esc given the observed trends seen at z ∼ 0 . 3 and z ∼ 3 ( § 3). LyC observations of galaxies with similar properties are necessary to better constrain their f esc .  \nIn contrast to the results for the Choustikov et al. (2024) model, the Mascia et al. (2023) model (Figure 8a and Table 10) reproduces f esc for the LzLCS+ galaxies reasonably well. This finding is unsurprising, since the Mascia et al. (2023) model is in fact derived from the LzLCS+ dataset. A Cox model run using the same variables ( β 1550 , r 50 , NUV , and O32) performs comparably to the Mascia et al. (2023) model (Figure 8b and Table 10). The largest difference occurs for weak LCEs ( f esc < 0 . 05) and non-detections. In these cases, the Mascia et al. (2023) model tends to systematically overpredict f esc , whereas by incorporating information from non-detections, Cox model predictions are more evenly  \ndistributed above and below the observed values. However, we note that the over-prediction of f esc in the Mascia et al. (2023) model is small, ∆ f esc of only a few percent.  \nDespite its derivation from the same dataset, the Cox and Mascia et al. (2023) models give different predictions at high redshift. The Mascia et al. (2023) model predicts f esc = 0 . 19 for Ion2, the only z ∼ 3 LCE with the required measurements, whereas the Cox models more accurately match its observed f esc > 0 . 5 ( § 3). This difference suggests that galaxy luminosity, incorporated in the Cox models in either Σ SFR or in M 1500 may be important in reproducing f esc for the strongest LCEs.  \nIn Figure 9, we compare the model predictions for GLASS-JWST galaxies (Treu et al. 2022) at z = 4 -8 from Mascia et al. (2023) with predictions from Cox models using similar combinations of variables. For a Cox model using the same variables as in Mascia et al. (2023), the predictions from Mascia et al. (2023) and the corresponding Cox model track each other (Figure 9a), aside from the aforementioned tendency of the Mascia et al. (2023) model to give higher predictions at low f esc . However, the Mascia et al. (2023) and Cox models disagree more strongly when Σ SFR is used as a variable in the Cox model instead of r 50 , NUV (Figure 9b). The inclusion of Σ SFR shifts some galaxies with compact sizes but weak star formation to much lower f esc . We conclude that including some measure of luminosity or SFR may be important to accurately identify non-leakers and extreme LCEs like Ion2.  \nMascia et al. (2024) also apply their multivariate model to galaxies from CEERS and find significantly different results than the predictions from our ELG-O32 Cox model for these galaxies ( § 4, Table 6). We compare these predictions in Figure 9c. Again, the disagreement primarily arises from the inclusion of luminosity in the ELG-O32 Cox model but not the Mascia et al. (2023) model. As with the GLASS sample, when we run a Cox model with the same input variables ( β 1550 , r 50 , NUV , and O32), our predictions agree closely with the Mascia et al. (2024) predictions (Figure 9a). Adding one additional variable, M 1500 , as a measure of luminosity begins to bring the predicted f esc into agreement with the ELGO32 predictions (Figure 9d). The remaining disagreement between the Mascia et al. (2024) and ELG-O32 Cox model predictions is due to the use of radius as a variable in the Mascia et al. (2023) model and different estimates of dust content and ionization from different publications (Mascia et al. 2024; Fujimoto et al. 2023; Tang et al. 2023).  \nObservational and theoretical studies agree that f esc depends on multiple physical parameters. However,  \n<!-- image -->  \nFigure 7. (a) The f esc predictions from the Choustikov et al. (2024) literature model compared with the LzLCS+ observations. (b) The f esc predictions from a Cox model run using the same variables as the Choustikov et al. (2024) model: β 1550 , E(B-V) neb , L (H β ), M 1500 , R23, and O32. Symbols are the same as in Figure 1. The Choustikov et al. (2024) model does not reproduce the f esc observations from the LzLCS+, but a Cox model using the same input variables does recover the observed f esc .  \n<!-- image -->  \n<!-- image -->  \nFigure 8. (a) The f esc predictions from the Mascia et al. (2023) literature model compared with the LzLCS+ observations. (b) The f esc predictions from a Cox model run using the same variables as the Mascia et al. (2023) model: β 1550 , r 50 , NUV , and O32. Symbols are the same as in Figure 1. Both models perform similarly, but the Cox model does slightly better at predicting f esc in weak LCEs and non-detections.  \n<!-- image -->  \nstudies have not yet reached a consensus as to which parameters matter and how. Predictions for f esc in z ≳ 6 galaxies from this work, Choustikov et al. (2024), Mascia et al. (2023) and Mascia et al. (2024) disagree because of different adopted scalings with O32 and with luminosity. Observationally testing and distinguishing between these predictions with larger samples at z < 6  \nwill be critical to reliably predict f esc in the epoch of reionization.", '6. IMPLICATIONS FOR REIONIZATION': 'As demonstrated in § 3, the empirical Cox models derived from the LzLCS+ show promise as diagnostics of f esc at high redshift. However, with the limited measurements at z > 6 available so far, the models do', 'Cox Model Variables: O32, r50,  1550': 'Figure 9. A comparison of f esc predictions for high-redshift galaxies from Cox models and the Mascia et al. (2023) model. (a) A Cox model using β 1550 , r 50 , NUV , and O32 gives similar predicted f esc as the Mascia et al. (2023) model for z > 4 galaxies from the GLASS-JWST (black) and CEERS (blue) surveys (Mascia et al. 2023, 2024). (b) If Σ SFR is substituted for r 50 , NUV in the Cox model, the Cox model f esc predictions for the GLASS galaxies deviate from those of the Mascia et al. (2023) model. (c) For the CEERS galaxies, the f esc predictions from the ELG-O32 Cox model disagree with the Mascia et al. (2024) predictions. (d) A Cox model using M 1500 in addition to the Mascia et al. (2023) variables begins to agree more closely with the ELG-O32 Cox model predictions. Adding a measurement of galaxy luminosity in the form of Σ SFR or M 1500 can significantly affect the f esc predictions.  \n<!-- image -->  \n(a)  \nComparison of  \nf  \nesc Predictions for  \nz  \n>6 CEERS Galaxies  \n<!-- image -->', 'Cox Model Variables: O32,  SFR,  1550': "3  \nPredicted  \nf  \nesc from Mascia et al. (2023)  \nComparison of  \nf  \nesc Predictions for  \nz  \n>6 CEERS Galaxies  \n10  \n3  \nPredicted  \nf  \nesc from  1550,  \nr  \n50, O32,  \nM  \n1500 Cox Model  \n(d)  \nnot decisively show which galaxy populations dominate reionization ( § 4.2). Because a combination of properties regulates a galaxy's optical depth, we need estimates of many factors to accurately predict f esc . Dust attenuation, galaxy morphology, ionization, and UV luminosity all play a role in promoting LyC escape. Based primarily on dust attenuation, we would expect faint galaxies to be strong LCEs (e.g., Chisholm et al. 2022; Atek et al. 2024). Nevertheless, at a fixed value of β 1550 , f esc still  \nshows considerable spread in the LzLCS+ and in multivariate predictions (Figure 6), spanning a range of up to ∼ 0 . 7 in f esc . Without estimates of other properties (e.g., Σ SFR , O32), we therefore cannot determine whether faint galaxies have f esc above or below the average LzLCS+ value for the same UV color.  \nAt the same time, we cannot yet rule out the contribution of more luminous galaxies. Several galaxies with M 1500 ≤ -20 (e.g., CEERS 498, 1027, 698, 1019;  \n10  \n2  \n10  \n1  \n10  \n0  \n(b)  \n10  \n2  \n10  \n1  \nl  \ne  \nd  \no  \nM  \nx  \no  \nC  \nm  \no  \nr  \nf  \nc  \ns  \ne  \nf  \nd  \ne  \nt  \nc  \ni  \nd  \ne  \nr  \nP  \nl  \ne  \nd  \no  \nM  \nx  \no  \nC  \n2  \n3  \nO  \n-  \nG  \nL  \nE  \nm  \no  \nr  \nf  \nc  \ns  \ne  \nf  \nd  \ne  \nt  \nc  \ni  \nd  \ne  \nr  \nP  \n10  \n10  \n10  \n10  \n10  \n10  \n10  \n10  \n0  \n1  \n2  \n3  \n0  \n1  \n2  \n3  \n10  \n10  \n0  \nJADES 18846; GLASS-JWST 100003, 10021) have the low dust content and high ionization suggestive of extreme f esc ≥ 0 . 5. The LzLCS+ data seem to suggest that for the same O32 ratio and β 1500 , a more luminous galaxy will have higher f esc (see also Lin et al. 2024). This effect could possibly be an observational bias, as the LyC flux should be easier to detect for more luminous objects. However, we note that at fixed O32, brighter LzLCS+ galaxies show higher ratios of 900 ˚ A to 1100 ˚ A flux, which suggests that their high f esc is genuine. Moreover, the trend of higher f esc with luminosity only appears in galaxies that have properties associated with LyC escape, such as high O32. Overall, the brighter galaxies in the LzLCS+ have fewer LyC detections (Flury et al. 2022b), which suggests that the detections are not biased toward the brightest galaxies. Furthermore, including UV luminosity in the multivariate models appears necessary to reproduce the f esc of Ion2 ( § 3 and § 5), which suggests that the increase of f esc with luminosity may be a real phenomenon.  \nIn the multivariate models, M 1500 might help to break the degeneracy between ionization parameter and optical depth in galaxies with high O32. For instance, lowmetallicity galaxies tend to have both low luminosities and inherently high ionization parameters (e.g., Nagao et al. 2006; Dopita et al. 2006), such that they may be more likely to have high O32 even at high optical depth. Alternatively, LCEs with higher UV luminosities could be able to ionize gas over a wider opening angle, so that we are statistically more likely to observe a high f esc line of sight.  \nOf course, to evaluate a galaxy population's influence on reionization, we must also consider how many ionizing photons are produced and the total LyC input into the IGM, not just f esc . Unfortunately, the Cox models do not do as well at predicting the ionizing to nonionizing UV flux ratio or the total LyC luminosity (Paper I). Galaxies with high O32 and high nebular EWs may possess elevated ξ ion values (e.g., Schaerer et al. 2016; Tang et al. 2019; Maseda et al. 2020; Naidu et al. 2022), such that at fixed f esc , these galaxies will provide more LyC photons to the IGM. These same properties scale with f esc in the multivariate Cox models, and our identified candidate strong LCEs at z > 6 generally have O32 > 10 and/or EW([O iii ]+H β ) > 1500 ˚ A. High f esc may therefore be coupled with high ξ ion (e.g., Schaerer et al. 2016; Naidu et al. 2022).  \nHowever, ionizing photon production may vary with other parameters as well, such as galaxy luminosity (e.g., Bouwens et al. 2016; Finkelstein et al. 2019). Fujimoto et al. (2023) estimate ξ ion for spectroscopically confirmed galaxies at z ∼ 8 -9 and find that  \nξ ion is two times higher at M UV ∼ -19 . 5 compared to M UV ∼ -21 . 5. Although this enhanced efficiency is not enough to compensate for the six times fainter UV luminosity at M UV ∼ -19 . 5 vs. -21.5, these fainter galaxies are also approximately 40 times more numerous (Bouwens et al. 2015). Atek et al. (2024) infer similarly high ξ ion for even fainter ( M UV ∼ -17 to -15) and even more numerous galaxies, which only need f esc < 0 . 05 to drive reionization. For the faintest galaxies we consider ( M 1500 > -17 . 5), the median of the Cox model predicted f esc values is near this threshold, with median f esc ∼ 0 . 04. However, these f esc estimates require refinement with additional parameters (e.g., O32, Σ SFR ). If f esc does not vary strongly with luminosity, as suggested by our preliminary, limited models, the fainter population would dominate reionization due to their higher ξ ion and greater numbers.  \nAssessing the main contributors to reionization requires progress on several fronts. At z > 6, we need estimates of nebular and morphological parameters in addition to β 1550 and M 1500 for galaxy samples spanning a wide range of luminosities. Parameters such as O32 and Σ SFR have some of the greatest effects on the f esc prediction accuracy for both the LzLCS+ (Paper I) and galaxies at z ∼ 3 ( § 3). The limited magnitude range of the LzLCS+ reference sample ( M 1500 ∼ -18 . 5 to ∼ -21 . 5) also introduces uncertainty. Our f esc estimates for spectroscopic samples within this magnitude range are likely reasonable and suggest moderately high median f esc ∼ 0 . 04 -0 . 05, with some galaxies reaching f esc as high as 0.6-0.7 (see Figure 4 and predictions in Table 6 for the ELG-O32 model and Table 7 for the ELG-O32β model). However, our f esc estimates for fainter and brighter galaxy populations rely on extrapolation and are less trustworthy. To confirm these estimates, we need to explore f esc and its dependence on galaxy properties across a wider magnitude range at z < 6. Lastly, observations of the IGM and galaxy population at z > 6 will provide a further test of f esc predictions. An accurate model of f esc should reproduce both the timing and topology of reionization with the z > 6 galaxy population. While much progress is needed to both improve and confirm f esc estimates, our current f esc predictions imply that plausible contributors to reionization appear at all magnitude ranges, and star-forming galaxies with the required levels of LyC escape do indeed exist at z > 6.", '7. SUMMARY': "Quantifying the LyC escape fraction ( f esc ) of galaxies is critical to understand the sources of cosmic reionization. We have developed a flexible tool for predicting  \nf esc empirically using combinations of observable variables available in the z ∼ 0 . 3 LzLCS+ reference sample (Paper I). We generate multivariate diagnostics for f esc with the Cox proportional hazards model (Cox 1972), a survival analysis technique that appropriately treats data with upper limits. We test Cox models developed from the z ∼ 0 . 3 galaxies on f esc observations at z ∼ 3, and we apply the models to several samples of z ≳ 6 galaxies to predict f esc for galaxies in the epoch of reionization. The Appendix provides all the information necessary to use the multivariate Cox models to predict f esc for other samples of galaxies, including samples at z ∼ 3 and samples at z ≳ 6, where the LyC is unobservable. 2  \nWe summarize our main findings below.  \n- 1. The models successfully reproduce the observed f esc values for the high-redshift z ∼ 3 galaxies and often have a lower RMS for the z ∼ 3 samples than for the z ∼ 0 . 3 galaxies. The success of these models suggests that low-redshift and highredshift LCEs may share similar properties. ( § 3)\n- 2. The best-performing models for the z ∼ 3 galaxies include the dust attenuation inferred from the UV SED (E(B-V) UV ) or the UV slope β 1550 , plus either O32 or a morphological measurement ( r 50 , NUV or Σ SFR ). However, this result is tentative; each model is tested on a different subset of high-redshift galaxies, which makes comparing the models difficult. To determine which variable combinations best predict f esc at high redshift, we require larger z ∼ 3 samples of LyC-emitters with the full suite of input variables. ( § 3)\n- 3. We generate new Cox models based on the LzLCS+ observations, which incorporate variables measured for z ≳ 6 samples. One model, using M 1500 , M ∗ , E(B-V) UV and EW([O iii ]+H β ) as input variables, applies to 180 galaxies from Endsley et al. (2021, 2023); Bouwens et al. (2023); Tang et al. (2023); Fujimoto et al. (2023); Saxena et al. (2023) and predicts a median f esc = 0 . 047 for this combined sample. Of these 180 galaxies, 27% have f esc > 0 . 1 and only 6% have f esc > 0 . 2. However, this set of variables results in predictions that tend to underestimate f esc , and the galaxy samples mostly consist of photometric measurements. A second model, using M 1500 , M ∗ , β 1550 , and  \nr 50 , NUV as input variables, applies to 278 galaxies, mostly with photometric redshifts, from Morishita et al. (2024) and Mascia et al. (2023). This model predicts a higher median f esc = 0 . 14, with 56% and 39% of the galaxies having f esc > 0 . 1 and f esc > 0 . 2, respectively. Many of these galaxies are more compact than the LzLCS+ sample, and this model therefore requires extrapolating beyond the LzLCS+ parameter space. ( § 4.1)  \n- 4. Smaller samples of spectroscopically-confirmed z ≳ 6 galaxies have higher predicted f esc , likely because they tend to include stronger emissionline galaxies. We use a model with β 1550 , M 1500 , and O32 to predict f esc for 27 spectroscopicallyconfirmed galaxies (Williams et al. 2023; Schaerer et al. 2022; Mascia et al. 2023; Saxena et al. 2023, 2024) and find that 33-41% of these galaxies have f esc > 0 . 2. For a smaller sample of 17 galaxies (Williams et al. 2023; Tang et al. 2023; Fujimoto et al. 2023; Saxena et al. 2023), using M ∗ as an additional input variable and E(B-V) UV instead of β 1550 , we find that 46-85% of the galaxies have f esc > 0 . 2. These models identify five galaxies at z > 5 . 9 that may have a line-of-sight f esc as high as 0.5-1: CEERS-698, CEERS-1019, CEERS-1027 from Tang et al. (2023) and JADES 18846 and JADES 9422 from Saxena et al. (2024). ( §§ 4.1, 4.3).\n- 5. Currently, the predicted f esc values for different galaxy samples show no strong trend with UV magnitude across the range of M 1500 = -16 to -22. However, many of these models are limited by a lack of spectroscopic information. The models suggest that low-luminosity galaxies ( M 1500 > -18) have at least moderate f esc ∼ 0 . 05. If this f esc is coupled with a high ionizing photon production efficiency, such faint galaxies could substantially contribute to reionization (e.g., Atek et al. 2024; Simmonds et al. 2024). Additional measurements of variables such as O32 and Σ SFR for more z ≳ 6 galaxies could enable more accurate f esc predictions for both bright and faint galaxies and might increase the numbers of suspected strong LCEs. ( § 4.2)\n- 6. The multivariate predictions for z ≳ 6 galaxies can differ strongly from f esc predictions based on a single variable. We predict f esc for the strong Ly α Emitter, JADES-GS-z7-LA, discovered at z = 7 . 278 (Saxena et al. 2023). Despite its exceedingly high f esc , Ly α , which would seem to imply a high f esc , the two Cox models that include  \nthe most available information (dust, Ly α , O32, and luminosity) predict that JADES-GS-z7-LA is a moderate LCE with f esc = 0 . 017 -0 . 085. Another example is the galaxy 11027 from Williams et al. (2023), which has O32 = 12, similar to the z ≳ 6 galaxies with predicted f esc ≥ 0 . 6. Yet, models incorporating other information about the galaxy, such as its mass (log 10 ( M ∗ / M ⊙ )=7.7), UV magnitude ( M 1500 = -17 . 4), and dust attenuation ( β 1550 = -2) predict f esc ≤ 0 . 02. However, we caution that some z ≳ 6 galaxies, including these ones, fall outside the parameter space covered by the z ∼ 0 . 3 galaxy sample, which make make our f esc predictions less reliable. ( § 4.3)  \n- 7. We compare our model predictions to single variable predictions (Chisholm et al. 2022) and to other multivariate f esc predictions derived from simulations (Choustikov et al. 2024) or observations (Mascia et al. 2023, 2024). The Cox models more accurately predict f esc in non-detections than the β 1550 model from Chisholm et al. (2022), and the inclusion of feedback-sensitive variables such as O32 and Σ SFR may better trace the effect of absorption from H i in LCEs with low dust content. Turning to multivariate models, the LzLCS+ sample implies a correlation between O32 and f esc , rather than the anti-correlation adopted in the Choustikov et al. (2024) model. As a result, we find different predicted f esc for high-redshift galaxies, and the Choustikov et al. (2024) model fails to reproduce the f esc of the LzLCS+ galaxies. In contrast, our predictions agree with the Mascia et al. (2023) model, when we use their same set of input variables. However, our high-redshift f esc predic-  \nons differ strongly when we add UV luminosity or Σ SFR as input variables.  \nLyC-emitting galaxies show similar physical characteristics at both low and high redshift, which suggests that the same processes may govern LyC escape across cosmic time. Because a variety of factors can influence LyC escape, accurately predicting f esc requires information about galaxy luminosities, dust, nebular properties, and ionization. Direct measurements of LyC at low redshift, combined with multivariate statistical models for f esc , can give insight into the possible f esc of galaxies in the epoch of reionization. Nevertheless, additional observations are necessary before these techniques can reach their full potential. At low redshift, LyC observations that push to new parameter space, especially fainter and brighter UV magnitudes, will aid the application of low-redshift results to high-redshift samples. Measurements of relevant galaxy properties, especially in the rest-frame optical, for larger samples of z ∼ 3 LCEs will better test the relative performance of these models at high redshift. Based on our current results, models using O32 and Σ SFR appear promising but are not yet applicable to large samples at z > 6. By measuring these properties for larger, representative galaxy samples at z > 6, high-redshift surveys can connect galaxies in the epoch of reionization with their LyCemitting counterparts at low redshift.  \nWe thank Aayush Saxena for providing information about the dust attenuation in JADES-GS-z7-LA, and we thank the anonymous referee for feedback that improved the paper's readability. AEJ and SRF acknowledge support from NASA through grant HST-GO-15626. STScI is operated by AURA under NASA contract NAS-526555. ASL acknowledges support from Knut and Alice Wallenberg Foundation. MT acknowledges support from the NWO grant 0.16.VIDI.189.162 ('ODIN').", 'APPENDIX': "In Tables 11-23, we provide the information necessary to predict f esc using the models discussed in this paper, and we give an example of how to use these models below. We have also developed a python script, available on GitHub 3 , that allows the user to generate and apply new Cox proportional hazards models to predict f esc using a desired set of observed variables (version 0.1.0 of the code is archived in Zenodo; Flury et al. 2024).  \nHere, we provide the parameters for the main models presented in Paper I and in this work. The data for the fiducial model from Paper I appears in Table 11. The best-performing model for the LzLCS+ sample is the fiducial model modified to use EW(H i ,abs), and these model parameters appear in Table 12. In Table 13, we list the parameters for  \nthe Paper I JWST model, which only uses variables accessible at z > 6. Two alternative JWST models use only the top-ranked accessible variables: β 1550 , Σ SFR , and either O32 or [Ne iii ]/[O ii ]. We provide the data for these models in Tables 14 and 15. Finally, Tables 16-25 present the parameters for the models used for the z ∼ 3 and z ≳ 6 f esc predictions discussed in Sections §§ 3-4.  \nEach table first lists the goodness-of-fit statistics for the LzLCS+ galaxies as a measure of the model performance. It then provides the fitted coefficients ( b i ) for each included variable and the reference values ¯ x i , which are the mean of the LzLCS+ x i values, where x i is the input variable ( β 1550 , log 10 (O32), etc). Finally, the table lists the baseline cumulative hazard function, HF 0 , calculated by the lifelines Cox Model fitting routine for each of the observed f esc values for LCEs in the LzLCS+. Together, these parameters predict f esc for a given galaxy with a set of input variables x as follows (see Section § 2.2 for more detail).  \nTo predict f esc for a galaxy, we first calculate the partial hazards function ph( x ) for its set of variables:  \nph( x ) = exp[ n ∑ i =1 b i ( x i -¯ x i )] . (1)  \nWe use the coefficients b i given in the model table, the galaxy's observed values x i for each variable, and the mean values ¯ x i also given in the model table. For example, for the galaxy Ion2, we have observed values of E(B-V) UV = 0, log 10 (O32)=0.851, and log 10 (Σ SFR /( M ⊙ yr -1 kpc -2 ))=2.03. Using b i (-10.083, 1.240, 1.369) and ¯ x i (0.158, 0.521, 0.705) from Table 16, we calculate  \nph( x ) = exp[ -10 . 083 ∗ (0 -0 . 158) + 1 . 240 ∗ (0 . 851 -0 . 521) + 1 . 369 ∗ (2 . 03 -0 . 705)] = 45 . 43 . (2)  \nWe then use this value of ph( x ) to scale the baseline cumulative hazard function HF 0 ( f esc ) and calculate the Survival Function S ( f esc ), which represents the probability that the true escape fraction, f esc , true , is lower than each tabulated f esc :  \nS ( f esc ) = exp[ -HF 0 ( f esc ) ∗ ph( x )] . (3)  \nFrom Table 16, for f esc = 0 . 8890, HF 0 is 0.003949 and S (0.8890) is therefore exp[ -0 . 003949 ∗ 45 . 43] = 0 . 8358. For f esc = 0 . 6247, HF 0 = 0 . 008253 and S (0 . 6247) = exp[ -0 . 008253 ∗ 45 . 43] = 0 . 6873, and so on, for each value of f esc . Thus, according to the 'TopThree model, there is a 68.73% probability that Ion2 has a measured f esc < 0 . 6247 and an 83.58% probability that Ion2's f esc < 0 . 8890. To find the median predicted f esc , we calculate S for each value of f esc in the table and find the value where S reaches 0.5; f esc is predicted to be above this value 50% of the time and below this value 50% of the time. For Ion2, we see that S (0 . 5838) = 0 . 5629 and S (0 . 4911) = 0 . 4577. We linearly interpolate to find where S ( f esc )= 0 . 5, finding that S (0 . 53) ∼ 0 . 5. The predicted f esc of Ion2 is then f esc = 0 . 53. Similarly, to determine the uncertainties in the predicted f esc , we find where S ( f esc ) reaches 0.159 and 0.841. According to the model, f esc will be between these f esc values 68% of the time. For some galaxies, S ( f esc ) is always > 0 . 5 for the tabulated f esc values, indicating a > 50% probability that f esc is smaller than the smallest tabulated f esc . This situation corresponds to an arbitrarily small predicted f esc , f esc ∼ 0. Conversely, if S ( f esc ) is always < 0 . 5, f esc is arbitrarily large and f esc ∼ 1.  \nTable 11. Fiducial Model (Paper I)Baseline Cumulative Hazard  \nTable 12. Fiducial Model with EW(H i ,abs) (Paper I)Baseline Cumulative Hazard  \nTable 13. Full JWST Model (Paper I)Baseline Cumulative Hazard  \nStatistics  \nR  \n2  \n0.34  \n0.29  \n0.46  \n0.83  \nVariable  \nb  \nlog  \nlog  \nβ  \n10  \n10  \n([O  \n(Σ  \n1550  \n]  \nλ  \n3727)  \n0.996  \n0.521  \n-  \n2  \nii  \nkpc  \n)  \n1.404  \n0.705  \n-2.274  \n-1.810", 'Baseline Cumulative Hazard': "Table 18. LAE-O32 Model  \na Positive EW(Ly α ) denotes net emission.  \nClark, T. G., Bradburn, M. J., Love, S. B., & Altman, D. G. 2003, British Journal of Cancer, 89, 232, doi: 10.1038/sj.bjc.6601118  \nCox, D. R. 1972, Journal of the Royal Statistical Society: Series B (Methodological), 34, 187,  \ndoi: https://doi.org/10.1111/j.2517-6161.1972.tb00899.x  \nCurti, M., D'Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425, doi: 10.1093/mnras/stac2737  \nDavidson-Pilon, C. 2019, Journal of Open Source Software, 4, 1317, doi: 10.21105/joss.01317  \nde Barros, S., Vanzella, E., Amor'ın, R., et al. 2016, A&A, 585, A51, doi: 10.1051/0004-6361/201527046  \nModel Parameters  \ni  \n-0.935  \n-19.84  \n¯ x  \ni  \n⊙  \nR  \n2  \nadj  \nStatistics  \nRMS  \nC  \n)  \n0.961  \n8.920  \n-15.987  \n0.157  \nii  \nTable 19. LAE-O32-nodust Model  \nDonnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023,  \nMNRAS, 518, 6011, doi: 10.1093/mnras/stac3472 Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJS, 167, 177, doi: 10.1086/508261 Duan, Q., Conselice, C. J., Li, Q., et al. 2024, MNRAS, 529, 4728, doi: 10.1093/mnras/stae872  \nEilers, A.-C., Davies, F. B., & Hennawi, J. F. 2018, ApJ, 864, 53, doi: 10.3847/1538-4357/aad4fd  \nEndsley, R., Stark, D. P., Chevallard, J., & Charlot, S. 2021, MNRAS, 500, 5229, doi: 10.1093/mnras/staa3370 Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312, doi: 10.1093/mnras/stad1919 Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415, doi: 10.1146/annurev.astro.44.051905.092514 Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192,  \ndoi: 10.1086/163225  \nTable 20. ELG-EW ModelBaseline Cumulative Hazard  \nFinkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015, ApJ, 810, 71, doi: 10.1088/0004-637X/810/1/71 Finkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36, doi: 10.3847/1538-4357/ab1ea8 Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13, doi: 10.3847/2041-8213/acade4 Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293 Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87, doi: 10.3847/1538-4357/ab2045  \nFlury, S., Jaskot, A., & Silveyra, A. 2024, LyCsurv, 0.1.0, Zenodo, doi: 10.5281/zenodo.11392442  \nFlury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJS, 260, 1, doi: 10.3847/1538-4365/ac5331 -. 2022b, ApJ, 930, 126, doi: 10.3847/1538-4357/ac61e4 Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJL, 949, L25, doi: 10.3847/2041-8213/acd2d9 Gazagnes, S., Chisholm, J., Schaerer, D., et al. 2018, A&A, 616, A29, doi: 10.1051/0004-6361/201832759  \nTable 21. ELG-O32 ModelBaseline Cumulative Hazard  \nGreen, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651, doi: 10.1093/mnras/sty1008 Hainline, K. N., Shapley, A. E., Kornei, K. A., et al. 2009, ApJ, 701, 52, doi: 10.1088/0004-637X/701/1/52 Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9 Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805, doi: 10.1093/mnras/stu936  \n- Isobe, T., Feigelson, E. D., & Nelson, P. I. 1986, ApJ, 306, 490, doi: 10.1086/164359\n- Izotov, Y. I., Orlitov'a, I., Schaerer, D., et al. 2016a, Nature, 529, 178, doi: 10.1038/nature16456\n- Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683, doi: 10.1093/mnras/stw1205  \nIzotov, Y. I., Schaerer, D., Worseck, G., et al. 2018a, MNRAS, 474, 4514, doi: 10.1093/mnras/stx3115 Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2021, MNRAS, 503, 1734, doi: 10.1093/mnras/stab612 -. 2018b, MNRAS, 478, 4851, doi: 10.1093/mnras/sty1378  \nTable 22. ELG-O32β ModelBaseline Cumulative Hazard  \nJames, B. L., Auger, M., Pettini, M., et al. 2018, MNRAS, 476, 1726, doi: 10.1093/mnras/sty315  \n- Ji, Z., Giavalisco, M., Vanzella, E., et al. 2020, ApJ, 888, 109, doi: 10.3847/1538-4357/ab5fdc\n- Johnson, B. D., Leja, J. L., Conroy, C., & Speagle, J. S. 2019, Prospector: Stellar population inference from spectra and SEDs, Astrophysics Source Code Library, record ascl:1905.025. http://ascl.net/1905.025\n- Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144, doi: 10.3847/1538-4357/abbd44  \nKakiichi, K., & Gronke, M. 2021, ApJ, 908, 30, doi: 10.3847/1538-4357/abc2d9 Katz, H., Saxena, A., Rosdahl, J., et al. 2023, MNRAS, 518, 270, doi: 10.1093/mnras/stac3019 Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, doi: 10.1146/annurev-astro-081811-125610 Kerutt, J., Oesch, P. A., Wisotzki, L., et al. 2024, A&A, 684, A42, doi: 10.1051/0004-6361/202346656  \n- Kimm, T., Blaizot, J., Garel, T., et al. 2019, MNRAS, 486, 2215, doi: 10.1093/mnras/stz989  \nTable 23. ELG-O32β -Ly α ModelBaseline Cumulative Hazard  \nLarson, R. L., Finkelstein, S. L., Kocevski, D. D., et al.  \n2023, ApJL, 953, L29, doi: 10.3847/2041-8213/ace619 Leitet, E., Bergvall, N., Hayes, M., Linn'e, S., & Zackrisson, E. 2013, A&A, 553, A106, doi: 10.1051/0004-6361/201118370 Leitet, E., Bergvall, N., Piskunov, N., & Andersson, B. G. 2011, A&A, 532, A107, doi: 10.1051/0004-6361/201015654  \nLeitherer, C., Ekstrom, S., Meynet, G., et al. 2014, ApJS, 212, 14, doi: 10.1088/0067-0049/212/1/14  \nLeitherer, C., Hernandez, S., Lee, J. C., & Oey, M. S. 2016, ApJ, 823, 64, doi: 10.3847/0004-637X/823/1/64  \nLeitherer, C., Schaerer, D., Goldader, J., et al. 2011, Starburst99: Synthesis Models for Galaxies with Active Star Formation, Astrophysics Source Code Library, record ascl:1104.003. http://ascl.net/1104.003  \nTable 24. R50β ModelBaseline Cumulative Hazard  \nLeja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170,  \ndoi: 10.3847/1538-4357/aa5ffe  \nLin, Y.-H., Scarlata, C., Williams, H., et al. 2024, MNRAS, 527, 4173, doi: 10.1093/mnras/stad3483  \nLiu, Y., Jiang, L., Windhorst, R. A., Guo, Y., & Zheng, Z.-Y. 2023, ApJ, 958, 22, doi: 10.3847/1538-4357/acf9fa Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42, doi: 10.1051/0004-6361/201323152  \nMa, X., Kasen, D., Hopkins, P. F., et al. 2015, MNRAS, 453, 960, doi: 10.1093/mnras/stv1679  \nMaji, M., Verhamme, A., Rosdahl, J., et al. 2022, A&A, 663, A66, doi: 10.1051/0004-6361/202142740  \nMakan, K., Worseck, G., Davies, F. B., et al. 2021, ApJ, 912, 38, doi: 10.3847/1538-4357/abee17  \nMarchi, F., Pentericci, L., Guaita, L., et al. 2018, A&A, 614, A11, doi: 10.1051/0004-6361/201732133  \nTable 25. β -Metals ModelBaseline Cumulative Hazard  \nAlvarez-M'arquez, J.,  \nMarques-Chaves, R., Schaerer, D., '  \net al. 2021, MNRAS, 507, 524, doi: 10.1093/mnras/stab2187  \n- -. 2022, MNRAS, 517, 2972, doi: 10.1093/mnras/stac2893  \nMarques-Chaves, R., Schaerer, D., Kuruvanthodi, A., et al. 2024, A&A, 681, A30, doi: 10.1051/0004-6361/202347411  \nMartin, C., Barlow, T., Barnhart, W., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4854, Future EUV/UV and Visible Space Astrophysics Missions and Instrumentation., ed. J. C. Blades & O. H. W. Siegmund, 336-350, doi: 10.1117/12.460034  \nMascia, S., Pentericci, L., Calabr'o, A., et al. 2023, A&A, 672, A155, doi: 10.1051/0004-6361/202345866  \n- -. 2024, A&A, 685, A3, doi: 10.1051/0004-6361/202347884  \nWitstok, J., Smit, R., Saxena, A., et al. 2024, A&A, 682, A40, doi: 10.1051/0004-6361/202347176  \nWorseck, G., Prochaska, J. X., Hennawi, J. F., & McQuinn,  \nM. 2016, ApJ, 825, 144,  \ndoi: 10.3847/0004-637X/825/2/144  \nXu, X., Henry, A., Heckman, T., et al. 2023, ApJ, 943, 94, doi: 10.3847/1538-4357/aca89a  \nZackrisson, E., Inoue, A. K., & Jensen, H. 2013, ApJ, 777, 39, doi: 10.1088/0004-637X/777/1/39", 'Model Parameters': 'i  \n¯  \nx  \ni  \nR  \n2  \nadj  \nRMS  \nC  \n⊙  \nStatistics  \nR  \n2  \n0.40  \n0.35  \n0.44  \n0.83  \nVariable  \nb  \nlog  \nlog  \nβ  \n10  \n10  \n([Ne  \n(Σ  \n1550  \n]  \nλ  \n3727)  \n1.315  \n-0.532  \n-  \n2  \n]  \nλ  \n3869/[O  \niii  \nSFR  \n/( M  \n-  \n1  \nyr  \nii  \nkpc  \n)  \n1.343  \n0.705  \n-2.192  \n-1.810  \nBaseline Cumulative Hazard  \nTable 15. Limited JWST Model (Paper I): β 1550 , log 10 (Σ SFR ), log 10 ([Ne iii ]/[O ii ])', 'REFERENCES': "Atek, H., Labb'e, I., Furtak, L. J., et al. 2024, Nature, 626, 975, doi: 10.1038/s41586-024-07043-6  \nBackhaus, B. E., Trump, J. R., Pirzkal, N., et al. 2024,  \nApJ, 962, 195, doi: 10.3847/1538-4357/ad1520  \nBarrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJL, 902, L39, doi: 10.3847/2041-8213/abbd8e  \nBassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2019, MNRAS, 483, 5223, doi: 10.1093/mnras/sty3320  \nModel Parameters  \ni  \n¯  \nx  \ni  \nR  \n2  \nadj  \nRMS  \nC  \n⊙  \nTable 16. TopThree ModelBaseline Cumulative Hazard  \nBecker, G. D., D'Aloisio, A., Christenson, H. M., et al. 2021, MNRAS, 508, 1853, doi: 10.1093/mnras/stab2696  \nBezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv e-prints, arXiv:2212.04026, doi: 10.48550/arXiv.2212.04026 Bian, F., & Fan, X. 2020, MNRAS, 493, L65, doi: 10.1093/mnrasl/slaa007  \nBlanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/aa7567  \nBorthakur, S., Heckman, T. M., Leitherer, C., & Overzier, R. A. 2014, Science, 346, 216, doi: 10.1126/science.1254214 Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140, doi: 10.1088/0004-637X/811/2/140 Bouwens, R. J., Smit, R., Labb´e, I., et al. 2016, ApJ, 831, 176, doi: 10.3847/0004-637X/831/2/176  \nBouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145  \nTable 17. LAE ModelBaseline Cumulative Hazard  \nBradburn, M. J., Clark, T. G., Love, S. B., & Altman, D. G. 2003, British Journal of Cancer, 89, 431, doi: 10.1038/sj.bjc.6601119 Bremer, J., & Dayal, P. 2023, MNRAS, 524, 118, doi: 10.1093/mnras/stad1844 Castellano, M., Pentericci, L., Cupani, G., et al. 2022, A&A, 662, A115, doi: 10.1051/0004-6361/202243348 Cen, R., & Kimm, T. 2015, ApJL, 801, L25,  \ndoi: 10.1088/2041-8205/801/2/L25  \nChisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182, doi: 10.3847/1538-4357/ab3104  \nChisholm, J., Gazagnes, S., Schaerer, D., et al. 2018, A&A, 616, A30, doi: 10.1051/0004-6361/201832758  \nChisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104, doi: 10.1093/mnras/stac2874  \nChoustikov, N., Katz, H., Saxena, A., et al. 2024, MNRAS, 529, 3751, doi: 10.1093/mnras/stae776  \nR  \n2  \n0.37  \n0.29  \n0.45  \n0.82  \nVariable  \nb  \nM  \n1500  \nlog  \n10  \n(  \nM  \nE(B-V)  \nlog  \n10  \n∗  \n/ M  \nUV  \n([O  \nEW(Ly  \nα  \n)(  \niii  \n]  \nλ  \n5007/[O  \n˚  \nA)  \na  \n]  \nλ  \n3727)  \n1.863  \n0.523  \n1.506e-02  \n68.16"}