# Inferring the HII region escape fraction of ionizing photons from infrared emission lines in metal-poor star-forming dwarf galaxies

L. Ramambason, Vianney Lebouteiller, A. Bik, C. T. Richardson, F. Galliano, D. Schaerer, C. Morisset, F. L. Polles, S. C. Madden, M. Chevance, et al.

► To cite this version:

L. Ramambason, Vianney Lebouteiller, A. Bik, C. T. Richardson, F. Galliano, et al. Inferring the HII region escape fraction of ionizing photons from infrared emission lines in metal-poor star-forming dwarf galaxies. Astronomy & Astrophysics - A&A, 2022, 10.1051/0004-6361/202243866. hal-03747400

HAL Id: hal-03747400

https://hal.science/hal-03747400v1

Submitted on 10 Feb 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Inferring the HII region escape fraction of ionizing photons from infrared emission lines in metal-poor star-forming dwarf galaxies

L. Ramambason<fnref n="1" />, V. Lebouteiller<fnref n="1" />, A. Bik<fnref n="2" />, C. T. Richardson<fnref n="3" />, F. Galliano<fnref n="1" />, D. Schaerer<fnref n="4" />, C. Morisset<fnref n="5" />, F. L. Polles<fnref n="6" />, S. C. Madden<fnref n="1" />, M. Chevance<fnref n="7" />,<fnref n="8" />, and I. De Looze<fnref n="9" />,<fnref n="10" />

Received 26 April 2022 / Accepted 11 July 2022

## ABSTRACT

**Context.** Local metal-poor galaxies stand as ideal laboratories for probing the properties of the interstellar medium (ISM) in chemically unevolved conditions. Detailed studies of this primitive ISM can help gain insights into the physics of the first primordial galaxies that may be responsible for the reionization. Quantifying the ISM porosity to ionizing photons in nearby galaxies may improve our understanding of the mechanisms leading to Lyman continuum photon leakage from galaxies. The wealth of infrared (IR) tracers available in local galaxies and arising from different ISM phases allows us to constrain complex models in order to estimate physical quantities.

**Aims.** Primitive galaxies with low metal and dust content have been shown to host a patchier and more porous ISM than their high-metallicity counterparts, with numerous density-bounded regions from where ionizing photons might leak out. To what extent this peculiar structure contributes to the leakage of ionizing photons remains to be quantitatively studied. Such effects can only be investigated by accounting for the complexity and inhomogeneity of the ISM. We aim to provide a new statistical framework to quantify various galactic observables by constraining a representative multiphase and multisector topology using a combination of 1D models.

**Methods.** To address these questions, we built a refined grid of models that include density-bounded regions and a possible contribution of an X-ray source. Using MULTIGRIS, a new Bayesian code based on Monte Carlo sampling, we combined the models as sectors under various assumptions to extract the probability density functions of the parameters and infer the corresponding escape fractions from H II regions ($f_{esc,HI}$). We applied this new code to a sample of 39 well-known local starbursting dwarf galaxies from the Dwarf Galaxy Survey.

**Results.** We confirm previous results that hinted at an increased porosity to ionizing photons of the ISM in low-metallicity galaxies and provide, for the first time, quantitative predictions for $f_{esc,HI}$. The predicted $f_{esc,HI}$ for low-metallicity objects span a large range of values, up to ~60%, while the values derived for more metal-rich galaxies are globally lower. We also examine the influence of other parameters on the escape fractions, and find that the specific star-formation rate correlates best with $f_{esc,HI}$. Finally, we provide observational line ratios that could be used as tracers of the photons escaping from density-bounded regions. Among others, we discuss the possible caveats of diagnostics based on [C II]158 µm in low-metallicity environments as we find a strong metallicity dependence of the fraction of [C II]158 µm emitted in the different ISM phases.

**Conclusions.** The new framework presented in this paper allows us to use suites of unresolved IR emission lines to constrain various galactic parameters, including the escape fraction of ionizing photons from H II regions. Although this multisector modeling remains too simple to fully capture the ISM complexity, it can be used to preselect galaxy samples with potential leakage of ionizing photons based on current and forthcoming spectral data in unresolved surveys of local and high-redshift galaxies.

**Key words.** galaxies: starburst – galaxies: dwarf – ISM: structure – radiative transfer – infrared: ISM – methods: numerical

## 1. Introduction

Young stars that have just been formed irradiate the surrounding gas and ionize the interstellar medium (ISM). Those pockets of ionized gas (i.e., H II regions) often dominate the emission at galactic scales in star-forming galaxies. At equilibrium, H II

regions are surrounded by atomic neutral hydrogen, the photodissociation regions (PDRs), which may recombine further from the stars to form H$_2$ in dense molecular clouds. In some cases, part of the ionizing radiation (Lyman continuum below 912 Å; LyC) can leak out of H II regions and irradiate a diffuse ionized gas (DIG) reservoir (Zurita et al. 2002; Weilbacher et al. 2018;

<footnote n="1">¹ Université Paris-Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: lise. ramambason@cea.fr</footnote>

<footnote n="2">² Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova University Centre, 106 91 Stockholm, Sweden</footnote>

<footnote n="3">³ Physics Department, Elon University, 100 Campus Drive CB 2625, Elon, NC 27244, USA</footnote>

<footnote n="4">⁴ Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland</footnote>

<footnote n="5">⁵ Instituto de Astronomía, Universidad Nacional Autónoma de México, AP 106, 22800 Ensenada, BC, Mexico</footnote>

<footnote n="6">⁶ SOFIA Science Center, USRA, NASA Ames Research Center, M.S. N232-12, Moffett Field, CA 94035, USA</footnote>

<footnote n="7">⁷ Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstrasse 12-14, 69120 Heidelberg, Germany</footnote>

<footnote n="8">⁸ Institut fur Theoretische Astrophysik, Zentrum für Astronomie, Universität Heidelberg, 69120 Heidelberg, Germany</footnote>

<footnote n="9">⁹ Sterrenkundig Observatorium, Ghent University, Krijgslaan 281 – S9, 9000 Ghent, Belgium</footnote>

<footnote n="10">¹⁰ Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>
Herenz et al. 2017; Bik et al. 2018; Menacho et al. 2019, 2021). Depending on the exact morphology of the gas distribution, the ionizing photons can freely travel on large scales to escape in the surrounding circum-galactic medium (CGM) and potentially in the inter-galactic medium (IGM).

The contribution of LyC-leaking galaxies to the total ion-
izing budget in the epoch of reionization (EoR, *z* ~ 6–9) is a
key element in our understanding of the reionization process.
Recent simulations indicate that such populations of numerous,
low-mass, LyC-leaking galaxies with average escape fractions
(*f*<sub>esc</sub>(LyC)) of 10–20% would be sufficient to reionize the whole
universe without invoking any other contribution from ionizing
sources such as active galactic nuclei (AGN; Robertson et al.
2013, 2015). Under favorable assumptions on ionizing photons
production and accounting for a subdominant contribution from
AGN, Finkelstein et al. (2019) find that even an average escape
fraction below 5% throughout the bulk of the EoR would be
enough to match observational constraints. Naidu et al. (2020)
propose an alternative model where reionization is not driven by
the lower mass galaxies but by a few (<5%) highly star-forming
galaxies with stellar masses above 10<sup>8</sup> M<sub>☉</sub> and extreme escape
fractions (the “oligarchs”).

Regardless of which galaxy population drives the reion-
ization process, the inclusion of binary stars, which provide
ionizing photons at later stellar evolution stages than sin-
gle stars, might also play a crucial role in providing ener-
getic photons over large timescales. Simulations of Ma et al.
(2016) and Rosdahl et al. (2018) find that this leads to sig-
nificantly higher time-averaged escape fractions of ionizing
photons. In low-metallicity environments hosting very mas-
sive stars, stellar-mass black holes might also contribute sig-
nificantly to the ionizing photon production (e.g., Mirabel et al.
2011). While escaping LyC photons have been directly observed
in the UV domain at redshifts below 0.5 (e.g., Bergvall et al.
2006; Leitet et al. 2013; Borthakur et al. 2014; Leitherer et al.
2016; Izotov et al. 2016a,b, 2018a,b; Wang et al. 2019, 2021;
Izotov et al. 2021; Flury et al. 2022a,b), and at *z* ~ 2–3 (e.g.,
Vanzella et al. 2015, 2016, 2018, 2020; Shapley et al. 2016;
de Barros et al. 2016; Steidel et al. 2018; Bian et al. 2017;
Fletcher et al. 2019; Rivera-Thorsen et al. 2019; Pahl et al.
2021, 2022) with observed *f*<sub>esc</sub>(LyC) ranging from 2 to 72%,
direct detections of leaking LyC radiation is not possible – or
extremely unlikely – above *z* ~ 4 due to the absorption by
neutral hydrogen in the IGM, preventing any direct observa-
tion of potentially LyC-leaking galaxies directly within the EoR.
This observational barrier makes it difficult to perform quan-
titative studies of primordial LyC-leaking galaxies, which are
much needed to understand their role in the reionization process.
Moreover, since LyC detections only probe a single line of sight,
the measured values are sensitive to viewing angle dependences.

To overcome these constraints, several indirect methods to
trace the escape fraction have been explored. They rely on the
fact that photons escaping from H II regions are quite sensi-
tive to the structure and properties of the surrounding gas. This
view is supported by hydrodynamical simulations of the ISM,
which account for an inhomogeneous gas distribution produced
by turbulence and/or stellar feedback (e.g., Fujita et al. 2003;
Trebitsch et al. 2017; Kimm et al. 2017, 2019; Kim et al. 2018;
Kakiichi & Gronke 2021; Yoo et al. 2020). First, far-UV absorp-
tion lines can serve as a promising proxy to infer the covering
fraction of H I gas, which places an upper limit on the amount
of escaping photons (e.g., Reddy et al. 2016; Gazagnes et al.
2018, 2020; Chisholm et al. 2018; Saldana-Lopez et al. 2022).
Another approach consisting in looking at far-UV colors selec-

tion diagrams has been proposed in, for example Vanzella et al. (2015) and Naidu et al. (2017). Finally, emission lines arising from different phases of the ISM can also be very useful probes of the global porosity to ionizing photons. In the UV domain, the Lyman alpha (Lyα) line profile has proven to be an interesting proxy with the presence of double-peaked or triple-peaked profiles being associated with LyC and Lyα leakage (e.g., Verhamme et al. 2015, 2017; Henry et al. 2015; Izotov et al. 2020; Maji et al. 2022). In particular, the velocity separation of the blue peak and red peak has been shown to strongly anti-correlate with the measured escape fraction of Lyα photons. Recent studies have also highlighted the interest of using proxies such as weak helium lines (Izotov et al. 2017) or the Mg II λλ2796,2803 Å doublet (e.g., Henry et al. 2018; Chisholm et al. 2020; Xu et al. 2022; Katz et al. 2022a) for probing the low-density lines of sight necessary to allow LyC-photons to escape.

In the optical range, several proxies can serve as indi-
cators of LyC leakage. In particular, line ratios involving
ions with different ionization potentials, produced at differ-
ent depths, have been proposed as indicators to discriminate
between radiation-bounded and density-bounded H II regions.
Radiation-bounded regions correspond to ionized spheres delim-
ited by their Strömgren radii set by the equilibrium between
production of photons by stars and ionization of the surround-
ing gas. Density-bounded regions are instead delimited by the
lack of matter, which sets their outer radius before the Ström-
gren radius. Hence, they allow part of the LyC-photons pro-
duced by stars to escape from H II regions. The oxygen line
ratio [O III]λ5007 Å/[O II] λλ3726, 3728 Å (O32) proposed by
Jaskot & Oey (2013) and Nakajima & Ouchi (2014) was suc-
cessfully used to select LyC-leaking candidates but no strong
correlation was found with the measured values of escape frac-
tion (see Izotov et al. 2018b; Naidu et al. 2018; Bassett et al.
2019; Nakajima et al. 2020, and discussions therein). Based on
a similar idea, the lack of emission from ions with low ioniza-
tion potentials like [S II] λλ6716,6731 Å has also been proposed
to target leaking candidates (Wang et al. 2019, 2021; Katz et al.
2020). This lack of emission of some low ionization species
was first interpreted as the signature of a density-bounded
galaxy where the outer part of H II regions were completely
stripped out. However, using simple photoionization models,
Stasińska et al. (2015) have shown that on average galaxies with
high O32 cannot have massive escapes of ionizing photons, since
low ionization lines like [O I]6300 Å are often also detected
in these galaxies, implying the presence of radiation-bounded
regions. Subsequently Plat et al. (2019) and Ramambason et al.
(2020) noted that several strong LyC-emitters show surprisingly
strong [O I]6300 Å emission, and proposed several explanations.
While Stasińska et al. (2015) and Plat et al. (2019) suggested
that such emission could be powered by the presence of AGN
or radiative-shocks, we proposed in Ramambason et al. (2020)
a 2-component model combining both density- and ionization-
bounded regions.

A complementary picture has been provided by studies
of galaxies in the local universe with resolved H II regions.
Recently, Della Bruna et al. (2021) estimated the average escape
fraction of ionizing photons from H II regions in NGC 7793
to be 67% from MUSE observations with a ~10 pc resolu-
tion. A large fraction of those escaping photons is, however,
likely reabsorbed within galactic scales and contributes to cre-
ate a DIG reservoir also seen with MUSE (Della Bruna et al.
2020). This picture is in line with recent PHANGS-MUSE
observations of resolved H II regions in nearby spiral galaxies (Belfiore et al. 2022; Chevance et al. 2022). Other indirect methods applied to local objects based on the mapping of the ionization parameter (e.g., Zastrow et al. 2011, 2013), on the estimation of the intrinsic ionizing photon production rate from resolved stars (Choi et al. 2020), and on the ionized gas kinematics (Eggen et al. 2021) are also suggestive of large escape fractions of ionizing photons from H II regions. Additionally, results from Polles et al. (2019) on the local starbursting galaxy IC10 indicate that the derived porosity depends on the spatial scale, with most clouds being matter-bounded at small scales (~60 to 200 pc), while larger regions become more and more radiation-bounded at galactic scales. This result highlights the complexity of the ISM in which the energy produced by feedback is deposited at various spatial scales and over different dynamical timescales, hence producing a highly inhomogeneous internal structure.

While the indirect methods relying on spatially resolved information can be applied in the local universe, they are not applicable in more distant unresolved galaxies. Disentangling the contribution from the ionized and neutral gas phase is a crucial step to better understand the interplay between ionizing photons and the surrounding ISM. To do so, the infrared (IR) domain offers interesting tracers, not only arising from H II regions but also from the PDR and molecular phases. Unfortunately, this part of the spectra is often inaccessible in samples of known LyC-leaking galaxies and only a few objects with LyC detections were also observed in the IR. As an alternative, recent studies (e.g., Cormier et al. 2012, 2015, 2019) have tried to constrain the covering factor of neutral gas (PDR covering factor) in local galaxies to quantify the porosity to ionizing radiation by estimating the fraction of gas residing in neutral atomic and molecular phases. This complex approach was first introduced in Péquignot (2008) that provided an unprecedentedly detailed analysis of the proto-typical, low-metallicity galaxy IZw 18 by developing refined multisector topological models representing the contribution of each phase to the total emission. Similar models were adapted and successfully applied to local objects (e.g., Haro 11: Cormier et al. 2012, IZw 18: Lebouteiller et al. 2017, IC 10: Polles et al. 2019), to a sample of local dwarf galaxies (Cormier et al. 2019) drawn from the Dwarf Galaxy Survey (DGS; Madden et al. 2013) and to resolved regions in the Small and Large Magellanic Clouds (SMC/LMC; Lambert-Huyghe et al. 2022). The results from those studies indicate that nonunity PDR covering factors of neutral gas are necessary to reproduce the emission lines of most local objects. Such findings appear at odds with the few UV observations that detect little to no LyC-leakage in the local universe (e.g., Bergvall et al. 2006; Leitet et al. 2013; Borthakur et al. 2014; Leitherer et al. 2016).

In this context, it becomes crucial to understand what properties of the ISM are responsible for its porosity to ionizing radiation and determine if and how well integrated emission lines of unresolved galaxies can be used to constrain their escape fractions of ionizing photons. The question is especially interesting in the context of high-redshift studies, as more and more galaxies are detected with facilities like ALMA and NOEMA at z ~ 4–9 (e.g., Inoue et al. 2016; Carniani et al. 2017; Walter et al. 2018; De Breuck et al. 2019; Hashimoto et al. 2019; Harikane et al. 2020; Falkendal et al. 2021; Bakx et al. 2020; Meyer et al. 2022). Such observatories and the advance of future ones such as the James Webb Space Telescope (JWST) are opening a new window to observe galaxies close to or within the EoR.

In this paper, we present a first application of MULTIGRIS (Lebouteiller & Ramambason 2022, hereafter LR22), a new Bayesian code designed to constrain multisector models using spectra of unresolved galaxies. Our work builds on previous studies (e.g., Cormier et al. 2012, 2019; Lebouteiller et al. 2017; Polles et al. 2019; Lambert-Huyghe et al. 2022) in which multisector models were constrained using frequentist methods. In particular, this paper is a direct continuation of Cormier et al. (2019) that used a $\chi^2$ minimization method to select the best-fitting configurations between 1- and 2-sector models and derived PDR covering factors. We revisit those results using the same sample of galaxies but adopting a new method. Using a Bayesian framework allows us to overcome some major issues of the $\chi^2$ method: difficulty to derive errorbars, sensitivity to outliers, impossibility to include complex priors etc. Most importantly, it allows us to infer probability density functions (PDF) of various parameters, including for the first time the escape fraction of ionizing photons from H II regions ($f_{esc,HI}$).

The paper is organized as follows. In Sect. 2 we present our sample and the tracers used in the analysis. The grid of models and the Bayesian code are presented in Sects. 3 and 4 and our results in Sects. 5 and 6. We discuss the limits and possible improvements of this new framework in Sect. 7. Our main conclusions are summarized in Sect. 8.

## 2. Sample

### 2.1. Overview

Our sample is drawn from the Herschel Dwarf Galaxy Survey (DGS; Madden et al. 2013), which gathers photometric and spectroscopic observations of 50 nearby (0.5–191 Mpc) galaxies in the far-infrared (FIR) and submillimeter domains performed with the Herschel Space Telescope. All of these galaxies were also observed in the mid-infrared (MIR) domain with the Spitzer observatory and spectroscopic measurements are available for all but 5 objects. We focus on a sub-sample used in Cormier et al. (2019) that selected 38 compact galaxies with at least three spectral lines detected in the MIR and FIR domain (~5 to 120 µm) among the 50 observed galaxies. The DGS sample exhibits a wide range of physical properties, which makes it an ideal laboratory to study the variation of ISM properties over a range of physical and chemical conditions. In particular, this sample is ideally suited to study how the escape fraction evolves in local, low-metallicity galaxies whose ISM may resemble primordial galaxies from the early universe. More specifically, the intense sources of radiation (see Sect. 2.2) as well as the low masses, compact sizes and metal-poor gas reservoirs (see Sect. 2.3) of these galaxies may, to some extent, resemble the physical and chemical conditions of the primordial dwarf galaxies. We note, however, that our selection criterion favors IR-bright galaxies hosting actively star-forming regions that have formed in a previously enriched ISM. The possible analogy with primordial galaxies should hence be taken with caution.

We use the line fluxes provided in Cormier et al. (2019) which combine Herschel/PACS (60–210 µm) data with Spitzer/IRS (5–35 µm) data. The latter are available for all galaxies in our sample except three (HS 0017+1055, UGC 4483 and UM 133). The fluxes correspond to integrated measurements for galaxies that were fully covered by the instrumental apertures, except for one galaxy. The only exception is NGC 4214, which was observed in two pointings (central and southern star-forming regions) that are studied separately in this study. The extraction procedures and corrections applied to extended
sources are detailed in Cormier et al. (2015). We corrected one line flux from Cormier et al. (2019) ([Ne V]14 µm for NGC 5253) that corresponded to a false detection due to an error in the fitting process. Finally, we instead use a detection upper limit at $2-\sigma$ of $0.112 \times 10^{16} \text{W m}^{-2}$. We also include the total IR luminosities ($L_{\text{TIR}}$) derived by modeling the dust spectral energy distributions (SEDs) in Rémy-Ruyer et al. (2015). All the available upper limits were used in the analysis, contrary to Cormier et al. (2019) that manually selected a suite of classical emission lines arising from H II regions and PDR. Our study aims to extend the multiphase picture of the ISM that was provided in Cormier et al. (2019) by including more lines arising from different phases of the ISM and tracing different physical processes. The list of observables used as constraints is summarized in Table 1.

Although it requires several lines to derive the various parameter values with reasonable uncertainties, the code can run properly as long as it is provided with at least one line and one upper limit (which are necessary to set the prior distributions). In this case, the given solution corresponds to a relatively wide PDF (defined in Sect. 4.3) and large errorbars. While little can be said about the parameter values of such individual objects, they do not, however, bias the analysis of the global trends in our sample. We thus decided to include all objects with at least one detection and one upper limit, which adds one more galaxy to the sample of Cormier et al. (2019) (HS 2352+2733, detected in [O III]88 µm with two upper limits on [C II]158 µm and $L_{\text{TIR}}$), hence leading to a total of 39 galaxies.

In addition to including upper limits for the first time, we can also enlarge the selection of lines since our method is more robust to outliers than the previously used $\chi^2$ method. The number of constraints available for each galaxy varies from 1 single detection (in addition to two upper limits) up to 22 detected emission lines. The exact number of detections and instrumental upper limits available for each object is provided in Fig. A.4. The choice of the suite of lines used in the analysis has an impact on the best solution that is selected by our code and, in turn, on the parameter values that are derived. Choosing the optimal suite of lines to be used and the minimal number needed to constrain a given parameter is a complex problem and we refer to LR22 for a more detailed discussion.

In this study we focus only on IR lines although most galaxies in our sample are also detected in the optical domain. We check a posteriori that our models predict values consistent with the available $H\alpha$ measurements (see Sect. 5.2). Combining optical and IR lines is possible with MULTIGRIS but would require an additional treatment of the dust attenuation and of systematic uncertainties due to instruments. Increasing the number of lines used as constraints can also raise specific issues regarding redundancies and the risk to over-constrain some parameters. We postpone this study to a future work and refer to LR22 for an example combining both IR and optical lines.

## 2.2. Radiation sources and feedback

The DGS galaxies are starbursting galaxies with prominent MIR and FIR emission lines that hint at the presence of a population of young UV emitting stars, which strongly irradiate the ISM (Madden et al. 2013). Some galaxies in our sample (15/39) have also been reported to host a population of Wolf-Rayet stars and signatures associated with massive stars have been reported in Schaerer et al. (1999).

Additionally, evidence suggests that several galaxies in our sample may host X-ray sources with the claimed detections

<table><caption>**Table 1.** IR tracers used as constraints and corresponding ionization potentials for ionic lines.</caption><thead><tr><th></th><th>Tracers</th></tr></thead><tbody><tr><td>Molecular hydrogen</td><td>H<sub>2</sub> S(0), H<sub>2</sub> S(1), H<sub>2</sub> S(2), H<sub>2</sub> S(3)</td></tr><tr><td>Neutral and ionized gas tracers</td><td>[O I]λ63,145 µm<br>[Fe II]λ17,25 µm (7.9 eV),<br>[Si II]λ34 µm (8.2 eV),<br>[C II]λ158 µm (11.3 eV),<br>Hαλ12 µm (13.6 eV),<br>[N II]λ122,205 µm (14.5 eV),<br>[Ar II]λ7 µm (15.7 eV),<br>[Fe III]λ23 µm (16.2 eV),<br>[Ne II]λ12 µm (21.6 eV),<br>[S III]λ18,33 µm (23.3 eV),<br>[Ar III]λ9,21 µm (27.6 eV),<br>[N III]λ57 µm (29.6 eV),<br>[S IV]λ10 µm (34.7 eV),<br>[O III]λ88 µm (35.1 eV),<br>[Ne III]λ15 µm (40.9 eV)<br>[O IV]λ26 µm (54.9 eV),<br>[Ne V]λ14,24 µm (97.1 eV),</td></tr><tr><td>Total IR luminosity</td><td>L<sub>TIR</sub> (1 µm–1000 µm)</td></tr></tbody></table>

**Notes.** We report the ionization potentials corresponding to the energy thresholds required to create the ion producing a given emission line, either by de-excitation or by recombination.

reported in the following papers: HS 1442+4250, VII Zw 403 (Papaderos et al. 1994; Kaaret et al. 2011; Brorby et al. 2014), NGC 1569, NGC 5253, NGC 4214 (Ott et al. 2005; Binder et al. 2015; McQuinn et al. 2018) and NGC 625 (McQuinn et al. 2018). Some of them even host Ultra-Luminous X-ray sources (ULX) with measured X-ray luminosities above $10^{39} \text{ erg s}^{-1}$: Haro 2 (Otíf-Floranes et al. 2012), Haro 11 (Prestwich et al. 2015; Gross et al. 2021), He 2–10 (Ott et al. 2005; Reines et al. 2011), I Zw 18 (Thuan et al. 2004; Ott et al. 2005; Kaaret et al. 2011; Kaaret & Feng 2013; Brorby et al. 2014) and SBS 0335-052 (Thuan et al. 2004; Prestwich et al. 2013). We note that in some objects, other physical mechanisms are considered to explain the X-ray emission and the claimed detections of compact objects have been actively debated (e.g., He 2–10; Cresci et al. 2017, HS 1442+4250; Senchyna et al. 2020, and NGC 5253; Zastrow et al. 2011, and references therein). Regardless of their exact nature (e.g., high-mass X-ray binaries, intermediate mass black hole or AGN), the contribution of such sources may have an important impact in low-metallicity, transparent environments in which high energy photons can travel freely over large scales. In I Zw 18, the galaxy with the lowest metallicity in our sample, Lebouteiller et al. (2017) have shown that the X-ray radiation emitted by a single point source ULX dominates the energy balance in the neutral gas over galactic scales.

We also expect X-ray heating of the ISM to produce specific spectroscopic signatures such as the presence of emission lines associated with very high ionization potentials (e.g., above 40eV) ions. Among the X-ray-sensitive tracers that we consider in this study (see Table 1), [Ne III] is detected in 33/39 sources and [O IV] detected in 17/39 sources. While we only have upper limits on [Ne V]λ24,14 µm, other lines produced by ions with high ionization potentials have been detected in the optical domain (e.g., [Fe IV], [Fe V], [Ne V] and even [Fe VI],
[Fe VII]; Izotov et al. 2001, 2004a,b). We note, however, that the origin of such lines is debated with two main hypothesis being either the presence of X-ray sources (e.g., Lebouteiller et al. 2017; Schaerer et al. 2019; Simmonds et al. 2021) or the presence of fast radiative shocks (e.g., Allen et al. 2008; Izotov et al. 2012). In this study we do not take into account shocks since most of the lines that we consider arise from the H II region and PDR where the dominant source of feedback is likely radiation pressure (Lee et al. 2016, 2019). In any case, although shocks may contribute to the emission of some of the lines we consider (e.g., lines associated with very high ionization potentials and H₂ lines), we do not have enough constraints to disentangle the contribution of shocks to the emission in the current study. For a similar reason, we do not study the impact of varying the cosmic rays (CR) rate in our models although CR might significantly contribute to the heating of neutral gas, especially at low-metallicity (Lebouteiller et al. 2017).

Even when nonstellar sources do not significantly contribute to the energy balance, their presence can be linked to mechanical feedback mechanisms that may facilitate the escape of LyC photons produced by stars. Any additional kind of feedback that modifies the gas reservoirs surrounding the stars (e.g., through fragmentation or creation of low-density channels) might affect the resulting global porosity to ionizing photons. The DGS sample exhibits a wide variety of feedback mechanisms, which impact the gas structure and kinematics. In particular, several indications of the ISM being disrupted have been discussed in previous studies that have reported signatures of outflows (in Haro 2: Otí-Floranes et al. 2012, Haro 11: Menacho et al. 2019), ionization cones (NGC 5253: Zastrow et al. 2011), superbubbles possibly associated with galactic winds (NGC 1569: Westmoquette et al. 2008; Sánchez-Cruces et al. 2015), and other signatures of strong feedback phenomena (NGC 1705: Zastrow et al. 2013, NGC 625: Cannon et al. 2004; McQuinn et al. 2018, UM 461: Carvalho & Plana 2018, NGC 4214: Martin 1998; McQuinn et al. 2018 and Pox 186: Eggen et al. 2021). Such dynamical effects strongly affect the chemical and physical conditions of the surrounding gas. They might even lower the metallicity of galaxies since metal-enriched outflows can remove newly formed metals from H II regions (e.g., Amorín et al. 2010; Hogarth et al. 2020). Another direct effect might also be the creation of low-density channels formed by gas ejection that may favor escaping photons.

## 2.3. Gas and dust properties

Our sample spans a large range of sub-solar<fnref n="1" /> metallicities ranging from $12+\log(O/H) = 7.14$ (~1/35 $Z_\odot$) up to 8.43 (~1/2 $Z_\odot$), derived from empirical strong line methods (Madden et al. 2013; Rémy-Ruyer 2013). Their dust-to-gas mass ratios also span a large range of values (0.07–0.33) and have been carefully investigated in Rémy-Ruyer et al. (2014, 2015) and Galliano et al. (2021) using continuum measurements from Spitzer and Herschel.

Madden et al. (2020) have shown that in such environments, the UV photons can penetrate deeper in the clouds and photodissociate CO, hence creating a layer of CO-dark gas. While H₂ and CO are only detected in a few galaxies in our sample, evidence suggests that their molecular gas reservoir could be largely underestimated when using CO lines as a tracer. Using in particular the [C II]158um line, Madden et al. (2020) developed a

method to estimate masses of CO-dark H₂ gas residing in the envelop where CO is photodissociated by strong radiation field and find that the DGS galaxies most likely host unseen molecular gas reservoir (>70% of the total H₂ in all the galaxies in our sample) that can explain their high star-formation rates (SFR; -2.2 < log SFR < 1.4) compared to the CO-based estimates of molecular gas content.

Additionally, some of the DGS galaxies are associated with large H I reservoirs with masses ranging from $10^7$ to ~$10^{11} M_\odot$ (Rémy-Ruyer et al. 2014) and specific gas masses ($M_{HI}/M_\star$) ranging from 0.03 to 17.3. Although this neutral component is important, little is known on the actual distribution of this gas in our sample. Some of the objects observed in absorption have been associated with large column densities reaching values greater than $10^{21} \text{ cm}^{-2}$ (see Table A.3 and associated references). Such high values completely rule out photons escaping along these lines of sight. However, previous studies have also shown the inhomogeneous distribution of the neutral component and line-of-sight effects (e.g., Gazagnes et al. 2020).

# 3. Models

## 3.1. Modeling strategy

One of the main challenges driving our modeling strategy is the necessity to deal with spatially unresolved observations in which the structure of the ISM is not directly accessible. In such cases, emission arising from the different phases of the ISM (e.g., H II regions, PDRs, molecular gas, DIG) is blended into one single beam.

To overcome this problem, we recover the underlying gas distribution in the different ISM phases from unresolved spectra. Using a multisector topology as a representative view of the galaxy, such as proposed in Péquignot (2008), it becomes possible to disentangle the relative contribution of each phase. Each “sector” corresponds to a fraction of a sphere of gas where radiative transfers are performed in a 1D continuous way throughout the H II region, PDR, and molecular zone. While this model simplifies the actual geometry of a galaxy, there is much to gain from combining sectors (even independent ones) as compared to a single model, especially when dealing with tracers coming from different phases or depths. The need to automate and generalize this multisector topological approach led to the development of a new code, MULTIGRIS, which allows a flexible combination of models within a grid. We briefly describe the code in Sect. 4 and refer to LR22 for a detailed description of the general strategy used in combining sectors.

The models used in the combination must account for the main physical and chemical mechanisms producing the line emission. Photoionization and photodissociation codes with complex chemical networks and refined prescriptions for radiative transfer are well suited to that purpose. They have been extensively used to study H II regions (e.g., with Cloudy; Ferland et al. 2017, MAPPINGS V; Sutherland et al. 2018) or PDR (e.g., with MeudonPDR; Le Petit et al. 2006). We chose to use Cloudy models that allow a consistent treatment of the emission line physics throughout the ionized phase, PDR, and molecular phase. Combining such models into a representative topology somewhat compensates for their simplistic geometry, which, by itself, cannot capture the full complexity of the multiphase ISM. In particular, several studies pointed out the necessity to combine several regions having different physical depths to successfully reproduce the emission lines of local objects (e.g., Binette et al. 1996; Péquignot 2008;

<footnote n="1">¹ We use the solar value from Asplund et al. (2009) of $12+\log(O/H)=8.69$.</footnote>
Lebouteiller et al. 2017; Cormier et al. 2012, 2019; Polles et al. 2019; Ramambason et al. 2020; Lambert-Huyghe et al. 2022). The present study is a direct continuation of the previous analysis from Cormier et al. (2019) of the DGS sample (presented in Sect. 2.1) using topological models. Based on 1-sector and 2-sector models, with varying PDR covering factors, they selected the best fitting models using a $\chi^2$ minimization. In each model, at least one sector was computed until $A_V = 5$ and the PDR scaling factor was defined a posteriori as a linear scaling, between 0 and 1, applied to the lines arising from the PDR. For 2-sector models, an additional sector stopping at the ionization front was added. Based on this definition, they find nonunity covering factors for most of the galaxies in our sample, and a strong correlation between metallicity and covering factor.

Similarly to Cormier et al. (2019), we model each galaxy as a weighted combination of sectors where the mixing weights represent the contribution of each sector (see Sect. 4.1). The models are computed with a fixed cluster luminosity and scaled to match the observed line fluxes used as constraints with a free scaling factor defined within the Bayesian model (see Sect. 4.2). The prescriptions used for this new grid are inspired by Cormier et al. (2019) with a few modifications. One major change is that the line luminosities are saved for each tracer in a cumulative mode, for each depth computed in Cloudy, allowing us to include density-bounded sectors (i.e., models stopping before the ionization front), which was not the case in Cormier et al. (2019). This setting also allows a refined sampling of the parameter controlling the depth of each model (see Sect. 3.2.3). We discuss the changes in the model grid and their implications in the following section.

## 3.2. Cloudy models

The grid used in this article was built using the photoionization and photodissociation code Cloudy v17.02 (Ferland et al. 2017). Each model consists of a spherical shell of gas placed at a fixed inner radius of the incident radiation source. We use a closed, spherical geometry taking into account the transmitted and reflected radiation, assuming a unity covering factor of the gas. The radiative transfer is computed along each line of sight (1D) in a continuous way throughout the H II region, PDR, and molecular zone. We summarize the main input parameters and the range of values spanned for each parameter in Table 2. The grid contains 28 800 models with an X-ray source and 3200 models without X-ray source, which results in a total of 32 000 Cloudy models. Each model is then truncated a posteriori with 17 cuts to create 544 000 sub-models. The current grid includes predictions for 516 emission lines from UV to IR.

### 3.2.1. Radiation field

In our grid we consider a stellar component and an X-ray component as represented in Fig. 1. The stellar population consists of a single burst population from BPASSv2.1 (Eldridge et al. 2017) that includes the contribution from binary stars. We note that the stellar population in Cormier et al. (2019) was instead computed using a continuous star-formation history (SFH) for a 10 Myr old cluster simulated with Starburst99 (Leitherer et al. 1999). Switching to a continuous SFH instead of single-bursts would result in a change in our ionizing spectrum as young O/B stars would provide ionizing photons over longer timescales. Since we consider only one single burst, our grid spans ages below 10 Myr, after which the H$\alpha$ emission drops. We stress, however, that for the lowest metallicity models, the contribution

<table><caption>**Table 2.** Input parameters of Cloudy models.</caption><thead><tr><th></th><th>Fixed parameters</th></tr></thead><tbody><tr><td>Stellar population</td><td>BPASS v2.1<sup>(1)</sup> (Eldridge et al. 2017)<br>Broken power-law IMF<br>α<sub>1</sub> (0.1 M<sub>☉</sub> − 0.5 M<sub>☉</sub>) = −1.30<br>α<sub>2</sub> (0.5 M<sub>☉</sub> − 100 M<sub>☉</sub>) = −2.35<br>Single burst, including binary stars</td></tr><tr><td>X-ray component</td><td>Multicolor blackbody<br>(Mitsuda et al. 1984),<br>T<sub>out</sub> = 10<sup>3</sup>K, variable T<sub>in</sub></td></tr><tr><td>Other fields</td><td>Cosmic microwave background<br>Cosmic ray background 0.5log</td></tr><tr><td>Grains</td><td>SMC graphite/silicate model<br>(Weingartner & Draine 2001)</td></tr><tr><td>Z<sub>dust</sub>, q<sub>PAH</sub></td><td>Values scaled with metallicity<br>based on the analytical fits from<br>Galliano et al. (2021)</td></tr><tr><td>X/O profiles <sup>(2)</sup></td><td>N<sup>(1)</sup>, C<sup>(1)</sup>, Ne<sup>(2)</sup>, S<sup>(2)</sup>, Ar<sup>(2)</sup>, Fe<sup>(2)</sup>,<br>Si<sup>(3)</sup>, Cl<sup>(4)</sup></td></tr><tr><td>Other elements</td><td>Table ISM from Cloudy</td></tr><tr><td>Density law</td><td>n<sub>H</sub> = n<sub>H,III</sub> × (1 + N(H)/[10<sup>21</sup> cm<sup>−2</sup>])</td></tr><tr><td>Stopping criteria</td><td>A<sub>V</sub>=10 or T<sub>e</sub> = 10K<br>8 varied parameters</td></tr><tr><td>log(L<sub>*</sub>/L<sub>☉</sub>)</td><td>[7, 9]</td></tr><tr><td>L<sub>X</sub>/L<sub>*</sub></td><td>[0, 0.001, 0.01, 0.1]</td></tr><tr><td>log(T<sub>X,in</sub>)</td><td>[5, 6, 7] K</td></tr><tr><td>log([O/H])</td><td>[-2.190, −1.889, −1.190,<br>−0.889, −0.667, −0.491,<br>−0.190, 0.111]</td></tr><tr><td>Age of burst</td><td>[1, 2, 3, 4, 5, 6, 8, 10] Myr</td></tr><tr><td>log(n<sub>H</sub>)</td><td>[0, 1, 2, 3, 4] cm<sup>−3</sup></td></tr><tr><td>log(U<sub>in</sub>) <sup>(3)</sup></td><td>[-4, −3, −2, −1, 0]</td></tr><tr><td>Cut <sup>(4)</sup></td><td>[0, 4], step=0.25</td></tr></tbody></table>

**Notes.** <sup>(1)</sup>Metallicity bins in BPASS are $Z = 10^{-5}$, $10^{-4}$, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.014, 0.020, 0.030 and 0.040 (i.e. 0.05% of solar to twice solar). We note that $Z_{☉} = 0.02$ while the default value we use here is $Z_{☉} = 0.014$. <sup>(2)</sup>The input ionization parameter, $U_{in}$, is set by varying the inner radius for a given luminosity. It is described in Sect. 3.2.3. <sup>(3)</sup>The cut parameter controls the stopping depth of a given model. It is described in Sect. 3.2.3.

**References.** <sup>(2)</sup>References for elemental abundances. (1): Nicholls et al. 2017, (2): Izotov et al. 2006, (3): Izotov & Thuan 1999, (4): fixed at −3.4.

from binary stars delays this drop in H$\alpha$ and that ages up to ~30 Myr (Xiao et al. 2018) should be explored. Considering that our solutions favor ages below 6 Myr and that the addition of a new bin in stellar age would significantly increase the size of our grid of models, we postpone this potential improvement to a future work. In practice, an older population of stars is present in many of the galaxies in our sample but single bursts of later ages alone would not match the spectral signatures that we study here. We note that considering a combination of bursts instead of single-burst model could have a substantial impact on the ionizing continuum, since mixed-age models, which allows contributions from extremely young stars, typically produce more ionizing photons and harder ionizing spectra than single-burst
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Stellar (**left**) and X-ray (**right**) incident SEDs of Cloudy models for a solar metallicity and bolometric luminosity of $10^9 L_\odot$. $Q_E$ shows the number of photons produced above a given energy $E_\nu$. The vertical dashed lines represent the ionization potentials of $H^+$ (13.6 eV), $O^{2+}$ (35.1 eV), $O^{3+}$ (54.9 eV) and $Ne^{4+}$ (97.1 eV). **Left-hand side:** the dashed lines represent BPASS models with binary stars while the solid lines of the same color correspond to BPASS models without binary stars for a single stellar population of the same age. The change in the spectra due to the inclusion of binary stars is represented by the shaded area in between both lines. The insert shows a zoom around the Lyman edge at 912 Å where the additional contribution to ionizing photons is visible for ages above 3 Myr and increases with the age of the burst. **Right-hand side:** the colors represent different inner temperature of the multicolor blackbody and the different linestyles to different percentage of stellar luminosity: 0.001%, 0.01% and 0.1%.

and constant star-formation rate models corresponding to similar ages (Chisholm et al. 2019).

The stellar evolutionary tracks incorporate the effect of mass transfers between members of binary systems and the stellar initial mass function (IMF) includes a distribution of binaries tuned to reproduce the binary fractions observed in the local Universe. We use a broken power-law IMF with two indices of -1.3 and -2.35 with a change of slope at 0.5 $M_\odot$, which is the default IMF in Eldridge et al. (2017). Although very massive stars with masses above 100 $M_\odot$ may exist in local, low-metallicity galaxies (Crowther et al. 2010; Wofford et al. 2021), such objects remain largely unconstrained in models. Since we do not need to invoke very massive stars to reproduce the IR lines considered in this study, we choose to use the BPASS default mass cut-off at 100 $M_\odot$.

In Fig. 1, we illustrate the effects of including binary stars by comparing the single-star and binary-stars BPASS (Eldridge et al. 2017) SEDs. A clear difference is visible at ages above 3–4 Myr where models that include binary stars produce more ionizing photons. This feature has also been pointed out in Xiao et al. (2018): while single-star and binary-star populations produce fairly similar hydrogen- and helium-ionizing spectra at young ages, the inclusion of binaries produces a shallower drop in ionizing flux than for their single-star counterparts at later ages. The inclusion of binaries with ages between 1–10 Myr has a most profound effect on lines with ionization potentials between 13.6 up to ~54 eV (see Table 2). Indeed, X-ray binaries provide high energy photons with an energy sufficient to power the emission of several species with very high ionization potentials.

However, it is well-known that these high ionization lines cannot be reproduced by classical photoionization models (e.g.,

Stasińska et al. 2015), even when including effects of binary stars (cf. Stanway & Eldridge 2019). Other sources producing a harder ionizing continuum need to be invoked to reproduce the high ionization potential lines observed in local and high-redshift galaxies (e.g., He II: Kehrig et al. 2018; Schaerer et al. 2019; Stanway & Eldridge 2019; Senchyna et al. 2020, 2021; or [C IV]: Stark et al. 2015; Senchyna et al. 2021, 2022). For example, high-mass X-ray binaries have been proposed in Schaerer et al. (2019) to explain the presence of ubiquitous He II emission. While Senchyna et al. (2020) find that the addition of an accretion disk associated with high-mass X-ray binaries is insufficient to fully explain the He II emission, several recent studies have argued that the addition of an X-ray source is still needed to simultaneously reproduce the emission lines arising from different ISM phases, including the emission from ions with high ionization potential (e.g., Simmonds et al. 2021; Olivier et al. 2021; Umeda et al. 2022).

Similarly, we find that in order to produce both [O IV] and [Ne V] emission, we need to add an X-ray source that produces harder ionizing photons than BPASS models alone. As discussed in Sect. 2.1, the exact nature of such compact objects is unknown (e.g., high-mass X-ray binaries, intermediate mass black hole or AGN). Moreover, the X-ray spectra themselves are widely unconstrained and rely on strong modeling assumptions (Simmonds et al. 2021). Hence, we choose to use a general prescription to model a compact object surrounded by an accretion disk. To do so, we used a multicolor blackbody spectrum as defined in Mitsuda et al. (1984). This spectrum is defined by an outer temperature fixed at $10^3$ K and an inner temperature varying from $10^5$ to $10^7$ K. The luminosity is set with respect to the stellar luminosity and varies from 0% to 10% of the stellar cluster luminosity. The right-hand side panel of Fig. 1 shows
the X-ray component when varying the inner temperature of the disk and the relative luminosity. Although simplistic, our prescription for the X-ray source is more general than that used in previous studies that considered single blackbody spectra (e.g., Lebouteiller et al. 2017; Cormier et al. 2019).

Finally, we include the cosmic microwave background (CMB) and CR background. Following Cormier et al. (2019) we use a CR rate ~3 times higher than the standard CR rate (2×10⁻¹⁶ s⁻¹; Indriolo et al. 2007) to account for the recent star-formation history in the DGS and match the observed [O I] IR line ratio. We note that CRs and X-rays both impact the ionization and heating in the PDR but we do not have means to disentangle both effects. This somewhat arbitrary choice for the CR rate has a strong effect at low metallicities (below 1/10 Z☉) where the emission of low ionization potential species and recombination lines is boosted even at large A_V, deep inside the molecular zone. However, in the range of A_V, density and temperature favored in the results presented here, the CR effects should not be significant. Because we have no means to discriminate between CR and X-ray effect, we do not further discuss their potential contribution in this paper.

### 3.2.2. Metal and dust abundances

Our grid is designed to match the abundance patterns of low-metallicity dwarf galaxies that belong to the group of blue compact dwarf (BCD) galaxies. The prescriptions used in our models are summarized in Table 2. The abundances for nitrogen and carbon are based on Nicholls et al. (2017) who derive analytical curves accounting for primary and secondary production. Their fit is derived from a large sample of stellar measurements in the Milky Way spanning a wide range of metallicities (6 ≤ 12+log(O/H) ≤ 9). Their analytical curves are compatible with the gas-phase measurements for carbon and nitrogen in the BCDs from Izotov & Thuan (1999), which have little depletion due to their poor dust content. Neon, sulfur, argon, and iron abundance profiles follow the regressions of Izotov et al. (2006), based on a sample of low-metallicity BCDs. To avoid extrapolating at high metallicities, we used flat profiles for metallicities above 8.2 for those four elements. Because of the low statistics of Izotov et al. (2006) for chlorine, we fixed the [Cl/O]<fnref n="2" /> value to the median of their sample (−3.4). The silicon abundance profile is based on Izotov & Thuan (1999). For all the other metals we used the values from the Cloudy ISM table, assuming they scale linearly with the oxygen abundance. The profiles used in our grid are given in Fig. A.1.

We apply the same scaling relative to the solar value to both the gas and stellar metallicity. In practice, since BPASS uses a slightly different solar value (Z☉ = 0.020) than the Cloudy one (Z☉ = 0.014), the metallicities are slightly offset. The helium abundance is fixed to the value provided in the ISM abundance set in Cloudy. The dust-to-gas mass ratio and abundance of polycyclic aromatic hydrocarbons (PAH) are computed following the metallicity-dependent prescription of Galliano et al. (2021) based on Bayesian dust SED fits of 798 galaxies (including our sample) using the code HerBIE (Galliano 2018) with the THEMIS grain properties (Jones et al. 2017). The dust-to-gas mass ratio hence follows a 4th degree polynomial relation at metallicity above 12+log(O/H) = 7.3 and scales linearly with metallicity below this threshold (see Eqs. (8)–(9) from Galliano et al. 2021).

Galliano et al. (2021) also provide an analytical prediction for the mass fraction of aromatic features emitting grains ($q_{AF}$), assuming those features are carried by small a-C:(H) grains. Here, we assume that such features are carried by PAHs instead and estimate the mass fraction of PAHs ($q_{PAH} = M_{PAH}/M_{dust}$) as $q_{PAH} \sim q_{AF}/2.2$. The abundance of PAHs is known to strongly vary under different physical conditions: they are especially sensitive to metallicity effects and to the strength of the interstellar radiation field (ISRF). Galliano et al. (2021) find that $q_{PAH}$ is primarily driven by metallicity effects while correlations with ISRF indicators are weaker. We adopt a single value of $q_{PAH}$ per metallicity bin, corresponding to the analytical fit from Galliano et al. (2021). We note that this prescription sets the total PAH abundance in our Cloudy models while the abundance profile follows the default one in Cloudy, that is, scaling as H⁰/H. This results in a profile where PAHs are mostly present in the PDR but completely destroyed in the H II regions and molecular zone. This assumption is consistent with observational studies that show that PAHs are absent in the ionized region (e.g., Relaño et al. 2018; Chastenet et al. 2019); however, PAHs can also be present in the molecular phase (e.g., Chastenet et al. 2019), which have not been considered in the Cloudy modeling. Considering that all the PAHs reside in the PDRs might in turn result in a slight overestimation of the PDR heating by photoelectric effect.

### 3.2.3. Geometry of H II regions

Each of our models consists in a spherical shell of gas around a central ionizing source, whose inner radius is set to match the input ionization parameter defined as:

$$U_{\text{in}} = \frac{Q(H^0)}{4\pi n c R_{\text{in}}^2}, \qquad (1)$$

where n is the hydrogen density at the inner radius, $U_{\text{in}}$ is the input ionization parameter, c is the speed of light, and $Q(H^0)$ is the total number ionizing photons produced by the central source defined as follows:

$$Q(H^0) = \int_{E_v=13.6 \text{ eV}}^{\infty} \frac{L_*(v) + L_X(v)}{\text{h}v} dv, \qquad (2)$$

where $L_*(v)$ is the stellar luminosity and $L_X(v)$ the luminosity of the X-ray source emitted with a given frequency v. We note that this definition of $U_{\text{in}}$ only constrains the ionization parameter at the inner irradiated edge and that the volume-averaged ionization parameter could be significantly different depending on the model parameters. In particular, we note that this definition is not appropriate to use for thick shells of gas for which the inner ionization parameter ($U_{\text{in}}$) at the illuminated front can be very different from the ionization parameter at the ionization front $U_{\text{out}}$. As an illustration, the relation between log $U_{\text{in}}$, log $U_{\text{out}}$ and the volume-averaged log $\langle U \rangle$ is provided in Fig. A.5. Since the inner radius is set automatically to match a given pair of input ionization parameter and gas density, the only possibility to change the geometry of the region is to change the incident luminosity. At fixed $U_{\text{in}}$, a large cluster luminosity result in a more shell-like geometry while a lower cluster luminosity result in a more filled-sphere geometry in which the gas lies closer to stars. Stasińska et al. (2015) show that such geometrical effects can affect the low-ionization to high-ionization line ratio (e.g., [O II]/[O III]) as the emission from outer regions is boosted in a compact configuration. Although this geometrical effect is only secondary for most lines arising from the H II regions, it has a

<footnote n="2">² The brackets indicate values normalized by the solar abundance ratio.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic view of the 17 cuts used to create sub-models.</figcaption></figure>

much stronger effect on lines emitted near the ionization front and in the PDR and molecular zone.

The bolometric luminosity of $10^9 L_\odot$ chosen by Cormier et al. (2019) corresponds to the case of a galaxy dominated by one single to a few giant H II regions. For the range of ionization parameters and stellar ages covered in our grid this corresponds to a thick shell geometry with log Hα between 39.8 and 40.7 (in erg s⁻¹). Such super-giant H II regions with Hα luminosities above $10^{39}$ erg s⁻¹ are numerous among BCD and starburst galaxies, but not always present (Youngblood & Hunter 1999). They populate the upper end of the observed Hα luminosity function derived from resolved galaxies (Bradley et al. 2006) and are likely to be optically thin (Pellegrini et al. 2012). We also consider a lower luminosity $L_{bol} = 10^7 L_\odot$ (37.8 ≤ log Hα (erg s⁻¹) ≤ 38.7), which corresponds to a typical H II region, similar to those that dominate the Hα luminosity function just before the cut-off value (log Hα = 38.6 erg s⁻¹) from Bradley et al. (2006). This case mimics a bursty star-formation where a few hundred to a few thousand compact clusters are responsible for the total emission.

We adopt the same density law as in Cormier et al. (2019) in which $n_H$ is nearly constant in the H II region and scales with column density above $10^{21}$ cm⁻². This law provides a simple first order prescription of a smoothly varying density, which can describe both the density profile expected in dynamically expanding H II regions (Hosokawa & Inutsuka 2005) and in the interior of turbulent molecular clouds (Wolfire et al. 2010).

One major change as compared to Cormier et al. (2019) is that luminosities are compiled in a cumulative fashion, meaning that we can access the intrinsic luminosity of each line at a given depth in the Cloudy model instead of considering only the total resulting luminosity, which obviously depends on the stopping criterion. This is a crucial step to study the escape fraction of ionizing photons as this parameter is sensitive to the stopping depth of the model, especially near the ionization front. As illustrated in Fig. 2, each initial Cloudy model is used to create 17 sub-models stopping at different $A_V$ controlled by the “cut” parameter. The original model is cut at the inner radius (cut = 0), ionization front (cut = 1), H₂ dissociation front (cut = 2), CO dissociation front (cut = 3) and outer radius (cut = 4). To sample the different phases (H II region, PDR, CO-dark H₂ region and CO-emitting H₂ region) defined by those cuts, three additional cuts are added between each integer i (cut = i + 0.25, i + 0.5, i + 0.75), equally spaced in $A_V$ between cut = i and cut = i + 1. We stress that stopping the model at a given cut and truncating it a posteriori are not strictly equivalent but our tests have shown that this is a secondary effect.

Our cut parameter is analogous to other parameters that were used to describe density-bounded models such as the Hβ fraction used in Stasińska et al. (2015) and Ramambason et al. (2020) or the zero-age optical depth to LyC photons ($\tau_\lambda$) used

in Lebouteiller et al. (2017) and Plat et al. (2019). Although the cut parameter does not have a physical meaning, it allows us to ensure a good sampling of the region near the ionization front and provides a simple characterization of density-bounded regions (cut < 1) vs. ionization-bounded models (cut > 1). Additionally, it is defined consistently regardless of the metallicity (as opposed to $A_V$).

Since each Cloudy model is cut a posteriori into sub-models, the stopping criterion is not crucial in this study. However, to incorporate the emission from different phases, models should be deep enough to include the neutral and molecular zone for each set of parameters. This transition is metallicity-dependent and requires going deep enough in $A_V$ for the lowest metallicity models. Most of our Cloudy models are computed until they reach a maximum $A_V$ = 10. Only the densest models cannot reach $A_V$ = 10 because their electronic temperature drops below 10 K before reaching this optical depth.

### 3.2.4. Escape fraction from H II regions

We calculate the escape fraction of ionizing photons from H II regions by using the ionizing continuum that is provided by Cloudy at each depth in the model. This continuum saves the number of photons per frequency bins, at a given depth, for all energies greater than 1 Rydberg. We compute the escape fraction in a given range of energy as follows:

$$f_{\text{esc,HII}}(E_1 \le h\nu \le E_2)(R) = \frac{L(E_1 \le h\nu \le E_2)(R)}{L(E_1 \le h\nu \le E_2)(R_{\text{in}})}, \quad (3)$$

where the numerator is the luminosity produced by the central sources reaching the radius R with energy between $E_1$ and $E_2$, and the denominator is the luminosity impinging the cloud at the inner radius $R_{\text{in}}$, within the same energy range. The number of ionizing photons at a given radius is calculated as:

$$L(E_1 \le h\nu \le E_2)(R) = 4\pi R^2 \int_{v_1}^{v_2} F_v(r)dv, \quad (4)$$

where $F_v$ in the flux of photons with frequency v at a given radius. In practice we consider $f_{\text{esc,HII}} = f_{\text{esc}}(1\text{Ryd} \le h\nu \le \infty)(R_{\text{cut}})$ that includes all ionizing photons even in the X-ray regime, with $R_{\text{cut}}$ the radius at which the sub-model is cut. This definition is equivalent to the ratio between the total observed emission below 912 Å and the intrinsic emission below 912 Å. We ensured that the sampling of the cut values results in a fine enough sampling for $f_{\text{esc,HII}}$ as well. This definition of $f_{\text{esc,HII}}$ consistently accounts for the absorption of photons in the gas and by dust. However, it assumes that the integrated emission of a galaxy is dominated by H II regions and that photons reaching the edge of a density-bounded region freely escape in the IGM. In practice, this tends to overestimate the global galactic escape fraction as some photons are likely to get reabsorbed by clumps or diffuse gas after having escaped from H II regions. This $f_{\text{esc,HII}}$ is to be considered as the escape fraction from H II regions and cautiously compared to other measurements, except for the two regions in NGC 4124 (see Sect. 2.1). While this definition corresponds to the true escape fraction from our model, accounting for the total ionizing luminosity below the Lyman edge, it differs from the definition commonly used in observational studies to measure $f_{\text{esc}}(\text{LyC})$, which most often relies on the flux ratio at 912 Å and 1500 Å. This difference should be taken into account when comparing with LyC measurements. This will be discussed in Sect. 7.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematic view of the ISM of a starburst galaxy and associated representative topology.</figcaption></figure>

## 4. MULTIGRIS runs

In this section we recall some of the most important features of the code used in the current study. A more detailed description of the code can be found in LR22.

### 4.1. Topological representation

The notion of topology introduced in Sect. 3.1 is especially important in the context of the escape fraction. Indeed, regardless of the nature of the physical mechanisms at the origin of escaping photons, it seems that the LyC escape fraction is tightly linked to the distribution of gas in the ISM. In particular, Gazagnes et al. (2020) have shown that the escape fractions of LyC and Ly$\alpha$ photons are quite sensitive to the distribution of neutral gas and dust, both at galactic scale and on specific lines of sight. This effect is even more pronounced for LyC photons that are easily absorbed by the neutral gas while Ly$\alpha$ photons scatter and are destroyed only by dust.

One of the main advantages of our grid of models compared to previous topological studies (see Sect. 1) is that it includes a free parameter that allows us to constrain the depth of each sector based on the observed line ratios (Sect. 3.2.3). This avoids having to fix an arbitrary stopping criteria that is not trivial to determine in studies that include PDRs and molecular regions (e.g., as in Cormier et al. 2019). The topological models aim to provide representations of the distribution of matter (number and covering factors of sectors) and phases (stopping depth) in a galaxy. The different phases we consider are represented in Fig. 3: an atomic ionized phase dominated by H II regions near star clusters, a neutral atomic phase dominated by a warm component around H II regions where hydrogen is photodissociated (PDR), and a neutral molecular phase concentrated in molecular clouds. Two phases are not accounted for in our Cloudy models: the DIG and the diffuse neutral gas (see discussion in Sect. 7).

In essence, topological models are not designed to determine the spatial distribution of the gas but only the relative contribution of each phase and sector. As shown in Fig. 3, we model a galaxy as a sum of representative star-forming regions that can be combined into one single representative model, assuming that the star clusters that dominate the emission share the same properties. In practice, we assume a single SED for each galaxy we model. This representation is especially well-suited for unre-

solved observations for which we precisely do not have access to spatial information. Additionally, most current photoionization codes cannot account for precise gas distribution that would require expensive 3D radiation transfer models. In our single cluster topological models, the emission lines are computed as a weighted linear combination of the emission from each sector as follows:

$$L = \sum_{i=1}^{N_{\text{sectors}}} w_i L_i, \qquad (5)$$

where $w_i$ is the mixing weight and $L_i$ the predicted luminosity of a given line in the ith sector.

This topological approach allows predictions for any other physical quantities available in Cloudy as long as it can be expressed as an analytical combination of the observables from individual sectors. We use the same formula for “extensive” quantities, meaning, that scale with luminosity (e.g., gas masses in the different phases, number of ionizing photons $Q$). In the current study, we consider configurations having a single cluster and different number of sectors (from 1 to 3). Under this hypothesis of a single cluster, the $f_{\text{esc,HII}}$ is also an extensive quantity and the averaged $f_{\text{esc,HII}}$ can be derived as follows:

$$f_{\text{esc,HII}} = \sum_{i=1}^{N_{\text{sectors}}} w_i f_{\text{esc,HII}}^i. \qquad (6)$$

For a configuration with $N$ sectors, the free parameters correspond to $N$ times the 8 free parameters in the Cloudy models from Table 2 determined for each sector plus the mixing weights $w_i$. In addition, the global luminosity scaling factor of the cluster is free. This amounts to a total of $9N+1$ parameters. The effective number of free parameters can be reduced by imposing priors that link some parameters together. To mimic the exposure of all sectors to a single cluster we impose that the stellar luminosity, X-ray luminosity, inner temperature of the of the X-ray-emitting accretion disk, metallicity, and age of the stellar burst are identical in all sectors. An additional constraint imposes that all $w_i$ sum up to 1 (no hole in the model). The effective number of free parameters thus reduces to $4N+5$ parameters.

The last constraint on the sum of $w_i$ translates into one major modeling assumption: we assume that the cluster is fully surrounded by the different sectors and that we can constrain its
total luminosity. If the cluster luminosity is unconstrained, the total luminosity scaling factor and the covering factors of each sector can be degenerate. In practice, this would result in quite different configurations that could all match a suite of emission lines but with different scaling factors. This caveat is particularly important for quantities that strongly vary with different configurations (in particular, the escape fraction) and will be discussed in Sect. 7. Hence, $L_{TIR}$, which is used as a constraint (see Table 1), is essential in our analysis to constrain the scaling factor and match the observed total luminosity.

## 4.2. Inference of the parameters

MULTIGRIS can be used with different Markov chain Monte-Carlo (MCMC) samplers, which are compared in LR22. We used the Sequential Monte Carlo Sampler (SMC) that is adapted to multidimensional grids with several likelihood peaks. Interpolating on all parameters is computationally expensive and raises some issues for models located at the edges of our grid. Instead, we use a nearest neighbor interpolation on all the parameters of the grid (described in Table 2) except for the metallicity, for which we perform a linear interpolation.

We define the likelihood of our data $p(O|\theta, M)$, where $O$ is the data and $\theta$ the set of parameters of a model $M$, by considering our suite of emission lines as independent identically distributed random variables (RV). Each RV is described as Gaussian distribution centered on the measured value and with a $\sigma$ corresponding to the uncertainty of the observation. Hence, the likelihood can be expressed as:

$$ \mathcal{L} = p(O|\theta, M) = \prod_{i=0}^{N_{\text{obs}}} N(\mu = O_i, \sigma^2 = U_i^2), \quad (7) $$

where $N$ is the number of emission lines with observed fluxes $O_i$ and uncertainties $U_i$. For undetected lines with instrumental upper limits, the Gaussian distribution is replaced by a half-Gaussian. We consider that upper limits correspond to a $2\sigma$ signal. To avoid possible biases due to lines detected with unrealistically small uncertainties, we force a minimal uncertainty of 10% for all lines.

In practice, single-sector models should be preferred if they are able to simultaneously reproduce all the emission lines of a given galaxy. However, galaxies having numerous lines are more likely to require a greater number of sectors. To which extent the addition of a supplementary sector improves the agreement with observations needs to be quantitatively studied. To compare models, one needs to estimate how well a model performs at predicting data. To do so, we estimate the marginal likelihood corresponding to each configuration by evaluating the model at the posterior distribution of the parameters. The marginal likelihood is defined as follows:

$$ \mathcal{L}_M = \int_{\theta} p(O|\theta, M)d\theta. \quad (8) $$

Because of the integration on all the free parameters $\theta$, this metric penalizes models that require a large number of sectors without significantly improving the agreement with the data. Although the code can be used with any number of components, those models become less and less likely to be selected by the marginal likelihood criterion. Additionally, the computing time required for the MCMC sampling step scales with the number of free parameters. We hence limit this study to combinations involving 1, 2, or 3 sectors. An example of a 4-sector model for

<figure><img src="image_5.png" /><figcaption>Fig. 4. Schematic view of the single cluster 3-sector configuration for He 2–10. The number of stars in the central cluster scales with the intensity of the radiation field, the proximity of the gas representing varying ionization parameters and the orange shading increases with density.</figcaption></figure>

1Zw 18, which includes additional constraints from the optical lines, is presented in LR22.

In Table A.1, we report the best (i.e., having the highest marginal likelihood) configuration for each galaxy, the corresponding accuracy at $3\sigma$ and the marginal likelihood values in logarithm, for a varying number of sectors. In practice, we combine the results obtained for configurations having 1, 2, or 3 sectors by considering a weighted mean of the 3 configurations with the weights are defined as follows:

$$ \alpha_{i \text{ sectors}} = \frac{\mathcal{L}_{M,i}}{\sum_{j=1}^{3 \text{ sectors}} \mathcal{L}_{M,j}}, \quad (9) $$

where $\mathcal{L}_{M,i}$ is the marginal likelihood associated with the configuration having $i$ sectors (with $1 \le i \le 3$). This weighted combination allows us to account for configurations having marginal likelihoods that are close to the best configuration. Since the weighted combination is performed directly on the MCMC draws, the uncertainties we derive for the parameters incorporate uncertainties on the combination of configurations, which can be considered as an uncertainty on the model itself. Other metrics could be used to define the weight, which are detailed in LR22. Using a different metric yields variations for a few individual objects in our sample but does not affect the global trends that we discuss in the following sections.

As an illustration, Fig. 4 shows the best configuration found by MULTIGRIS for He 2–10. This is one of the galaxies for which we have the most constraints available (22 detections and 6 upper limits, see Fig. A.4). For this galaxy, the 3-sector configuration is favored. The left-hand side plot represents the relative contribution of each of the 3-sectors and the right-hand side plot shows another version where the sectors have been randomly redistributed around the central source. We stress that the two topologies are completely equivalent in our study.

## 4.3. Probability density functions

One of the advantages of using a Bayesian method is that the ensemble of solutions can be represented by a continuous distribution, which is well-suited to identify possible degeneracies or multimodal solution. While frequentist methods like the $\chi^2$ only provide a group of the most likely values of a parameter, MULTIGRIS outputs a sample of draws from the posterior PDF $p(O|\theta, M)$ of any given parameter.

The PDF can be described using various estimators. The mean, mode, and median values of the parameters can be significantly different in the case of multimodal or asymmetric distribution. Ideally, one would like to represent the full posterior PDFs that are used in this study to derive trends and correlations.
tions. We also investigate the combined PDF concatenating all sources. To that purpose, we use kernel density estimate (KDE) plots, which provide a smoothed convolution (using a Gaussian kernel) of the 2D-PDF of our whole sample. Additionally, to represent 2D-PDFs of individual objects we adopt the skewed uncertainty ellipse (SUE) representation (see e.g., Appendix F from Galliano et al. 2021). A SUE represents the $1\sigma$ contour of a 2 dimensional split-normal distribution adjusted to have the same three first moments as the underlying PDF. The center of the SUE marks the location of the robust mean of the 2D-PDF. While the representation is convenient to locate the parameter space of highest probability (e.g., $1\sigma$) for a given object, SUEs are not always representative of the underlying PDF, especially when several modes are present. In the following plots we show either the KDE of the full sample or the individual SUEs corresponding to each galaxy. While the KDE of the full sample should always be used to study the statistical trends and correlations, the individual SUEs can also be useful to visualize the region of maximal likelihood associated with each object.

# 5. Consistency checks

The escape fraction of ionizing photons is a complex parameter with multiple dependences. Before diving into the interpretation of this complex observable (see Sect. 6), we perform consistency checks to see how the predictions from MULTIGRIS compare to classical diagnostics found in the literature. In particular, since the escape fraction of ionization photons is sensitive to both the gas and stellar content of a galaxy, we first examine some key parameters that may control the predicted values of $f_{esc,HI}$: the metallicity and the star-formation rate.

## 5.1. Metallicity estimates

To ensure that our metallicity estimates are not affected by the metallicity sampling of the model grid (see Sect. 3), we perform a linear interpolation on this parameter. For the other parameters, we take at each draw the nearest neighbor in the grid. Our approach differs from classical methods that usually rely on a single tracer or a combination of a few lines to derive the metallicity. Instead, our code uses the combined information from all emission lines to constrain model parameters, including the metallicity. Although this approach is more flexible (no need for a specific suite of emission lines), we need to rely on a particular abundance pattern because we do not have enough constraints on the abundance of all species. In this study, we use abundances tailored for compact, low-metallicity, dwarf galaxies, as described in Sect. 3.2.4.

The metallicity could be subject to multiple degeneracies. In particular, varying the cut parameter, that is, the $A_V$ at which each sector stops, impacts the emission of tracers used to estimate the metallicity. Previous works have reported that accounting for escaping ionizing photons tends to bias the metallicity estimates toward lower values (Xiao et al. 2018; Jiang et al. 2019). The metallicity is also degenerate with the stellar age parameter. As mentioned in Xiao et al. (2018), an old stellar population with no leakage can produce line ratios similar to a younger stellar population where part of the radiation is leaking through density-bounded regions.

To avoid such degeneracies between metallicity, stellar age and cut parameters, we allow the metallicity to vary under a weakly informative prior. This ensures the inference preferentially starts around the mean of the prior and prevents drawing too far from it. We use a Gaussian centered at the measured

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Predicted metallicity of MULTIGRIS vs. measured metallicity from the strong lines R23-P method from Pilyugin & Thuan (2005) based on optical oxygen lines ratio. The black dots represent the robust means of the posterior PDF. The $1\sigma$ contours of the SUEs are calculated by simulating mock data for the y-axis with a dispersion corresponding to the measured uncertainty. The vertical dashed lines correspond to the metallicity bins of the grid. The galaxies are labeled with the numbers reported in Table A.1. The numbers that are flagged with a star correspond to two pointings in NGC 4214, which are excluded from the KDE. The color code indicates whether the Humphreys $\alpha$ line (7–6) is detected or not.</figcaption></figure>

metallicity with a large width parameter ($\sigma = 0.1 + \sigma_{mes}$ in logarithmic scale, where $\sigma_{mes}$ is the uncertainty associated with the measured metallicity). We emphasize that although we use only constraints in the IR domain, the prior on metallicity incorporates information derived from optical lines since it is centered on the metallicity estimates derived from optical strong lines method. Galaxies that do not fall within the $\sigma = 0.1$ dex envelope are galaxies for which the information given by all emission lines have pushed the posterior PDF away from the prior distribution.

In Fig. 5, we compare this estimated metallicity to that measured from strong lines methods. The galaxies are labeled with the numbers reported in Table A.1. The measured metallicities from the literature come from Madden et al. (2013), and references therein, based mainly on the optical oxygen lines ratio (R23-P) method from Pilyugin & Thuan (2005). Their method is based on a two-parameter calibration involving the R23 ratio and an excitation parameter P, which corrects the estimates by accounting for the physical conditions in the H II regions. The method used in Madden et al. (2013) allows them to derive metallicity measurements for almost all the DGS galaxies except for three galaxies (HS 0017+1055, HS 0822+3542 and HS2351+2733) for which the R23-P method cannot be applied and the measured metallicity is derived using the $T_e$-direct method (Ugryumov et al. 2003; Izotov et al. 2006, 1994). We note that the latter method is more reliable than the R23-P method (e.g., Maiolino & Mannucci 2019) that was chosen only because it could be applied in a consistent way to almost all objects in the DGS sample.

The uncertainties derived ($1\sigma$ uncertainties of SUEs defined in Sect. 4.3) from our method are somewhat larger than those found in the literature, as they incorporate uncertainties on all
emission lines and on the modeling (e.g., in particular our pre-
scription for the abundance patterns vs. metallicity). When a
hydrogen recombination line (e.g., the Humphreys $\alpha$ line (7–6);
$H\alpha$) is available, the estimated metallicity tends to be more
tightly constrained as the code can directly infer elemental abun-
dances relative to hydrogen. This is illustrated in Fig. 5 where
galaxies that have $H\alpha$ detection tend to have smaller $1\sigma$
contours than the ones without. Interestingly, our code is also able to esti-
mate metallicities without hydrogen recombination line measure-
ments, in some cases with a precision comparable to galaxies with
$H\alpha$ measurements. As a consequence, the metallicity is inferred
by constraining metal-to-oxygen abundances, which depend on
metallicity in the grid (see Sect. 3). In this case the results we get
are strongly dependent on the assumed metal-to-oxygen profiles,
and in particular the N/O and C/O vs. O/H relations.

The metallicities are consistent with the corresponding
strong line measurements within 0.1 dex for 31 out 40 galaxies
and withing 0.2 dex for all galaxies except one (12: II Zw 40),
for which we predict a significantly smaller metallicity than that
derive with the R23-P method. We note, however, that our mea-
surement is consistent with the value measured for this galaxy
(8.09 ± 0.02) using the $T_e$-direct method from Izotov et al.
(2006). The scatter we obtain is somewhat smaller than that
derived for the DGS sample in Madden et al. (2013) when
comparing different calibration methods using the R23 ratio
(Pilyugin & Thuan 2005) and the $T_e$-direct method. (Izotov et al.
2006). We find that our predictions tend to systematically pre-
dict slightly higher metallicities than the R23-P method for
metallicities below 7.8. This might be due to the fact that we
find galaxies with numerous density-bounded regions (i.e., with
small cut parameters) among the lowest metallicities in our
sample (see Sect. 6.1). This effect might be linked to the fact
that the R23-P diagnostic relies on an empirical calibration based
on samples of galaxies, regardless of the presence or not of
density-bounded H II regions. If the calibration is dominated by
radiation-bounded H II regions, this diagnostic might slightly
underestimate the metallicity of galaxies in which numerous
density-bounded regions are present. Nevertheless, based on the
weak prior assumption and on the emission lines given as con-
straints, we find that our code infers metallicities in agreement
with previous measurements for all galaxies in our sample. This
is a necessary condition to derive parameters such as $f_{\text{esc,HII}}$
which is expected to be metallicity-sensitive.

5.2. Star-formation rate estimates

We now use the predictions provided by MULTIGRIS to infer
SFR estimates in our sample and compare them to classical diag-
nostics. Although the SFR is not one of the primary (used for
inference) or secondary (output from Cloudy) parameters, we
use classical SFR proxies ($H\alpha$, $L_{\text{TIR}}$ and $Q(H^0)$), which we con-
vert into SFR. We emphasize that $H\alpha$ observations are not used
as constraints in our models but our code can provide PDFs
of both observed and unobserved emission lines. In this work,
we use the constraints from IR lines that, to first order, are not
affected by extinction by dust, to estimate the number of ion-
izing photons produced by the central cluster, $Q(H^0)$. The pre-
dicted $H\alpha$ luminosity depends of the incident flux $Q(H^0)$ and
on the geometry of the gas in our representative model. The
main advantage of our method is that it takes into account infor-
mation coming from all the emission lines available. However,
like for other diagnostics, our estimates ultimately depend on the
assumption made regarding the stellar population.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Agreement between $H\alpha$ and SFR predictions with empirical measurements. **Top:** predicted intrinsic $H\alpha$ emission vs. measured $H\alpha$ corrected for Galactic extinction from (Rémy-Ruyer et al. 2015, and references therein). **Bottom:** predicted SFR vs. SFR($H\alpha+L_{\text{TIR}}$) from Rémy-Ruyer et al. (2015). The $1\sigma$ contours of the SUEs are calculated by simulating mock data for the y-axis with a dispersion corresponding to the measured uncertainty. For 4 galaxies (HS 0017+1055, HS 0052+2536, HS 1319+3224 and HS 2352+2733) no $H\alpha$ measurements are available. For 2 of those galaxies (5: HS 0052+2536 and 9: HS 1319+3224), their measured SFR is derived in using diagnostics based on FUV.</figcaption></figure>

We compare our predictions to a classical diagnostic based
on composite diagnostics that combine tracers of dust-obscured
(e.g., *L*<sub>TIR</sub>) and dust-unobscured (e.g., *H*<sub>α</sub>) star-formation (see
review from Calzetti et al. 2012). First, we compare the pre-
dicted *H*<sub>α</sub> luminosities to the values from the literature gathered
in Rémy-Ruyer et al. (2015) in Fig. 6. The latter correspond
to *H*<sub>α</sub> luminosities corrected for underlying stellar absorption,
[N II]6548,6584Å lines contamination, and foreground Galactic extinction. No correction has been applied to account for the intrinsic attenuation within galaxies. For consistency, we compare the corrected measurements to the predicted intrinsic $H\alpha$ emission i.e. what escapes from the galaxy without considering dust extinction.

The predicted $H\alpha$ luminosity is compatible with observations within 0.5 dex for all galaxies in our sample except for two galaxies (14: Mrk 153 and 2: Haro 3) for which our prediction is significantly larger. The fact that our models find $H\alpha$ values close to observed measurements is remarkable, especially since no optical lines were used as constraints. This means that the intrinsic values we predict are relatively robust even though we do not properly account for $H\alpha$ emitted by the DIG (see Sect. 7). It also indicates that taking into account the actual geometry of the gas by including density-bounded regions in our models only yields minor variations in the predicted intrinsic $H\alpha$ emission and remains compatible with observations using somewhat simple prescriptions to correct the observed $H\alpha$ emission for foreground Galactic extinction. We stress that the intrinsic fluxes predicted by our models do not account for internal attenuation within the galaxy. The fact that our predictions are in agreement with observations that were not corrected for internal attenuation suggests that this additional attenuation term has only a minor effect in low-metallicity galaxies. This could, however, be an important effect for more metal-rich sources and may also explain part of the scatter see in Fig. 6.

We then use our line predictions to estimate the SFRs. To allow the comparison with the estimates provided in Rémy-Ruyer et al. (2015), we use the same empirical calibration from Kennicutt et al. (2009) based on mixed tracers ($H\alpha$ and $L_{TIR}$). This diagnostic corresponds to an SFR conversion assuming a Kroupa IMF<fnref n="3" /> and with coefficients calibrated on the Spitzer-SINGS sample. For consistency, we use $L_{TIR}$ luminosities predicted by our code. We note, however, that the observed $L_{TIR}$ is used as a constraint and always well reproduced by our models within 0.1 dex. Similarly to the $H\alpha$ prediction, we find that our predictions for the SFR are consistent with measurements from Rémy-Ruyer et al. (2015) within 0.5 dex for all galaxies but one (14: Mrk 153), which is linked to our overestimation of $H\alpha$ emission<fnref n="4" />.

One advantage of our approach is that it allows us to account for the photon losses due to escape fractions in density-bounded sectors, which are ignored in empirical calibrations. In Fig. 7 we convert the incident flux of ionizing photon, $Q(H^0)$, into SFR using the analytical formula provided in Calzetti et al. (2012) assuming a Kroupa IMF. We then compare it to our prediction based on $H\alpha + L_{TIR}$. Our SFR estimates based on $Q(H^0)$ are in most cases shifted to higher values. If we consider SFR($Q(H^0)$) as the intrinsic SFR, the correction to apply to SFR($H\alpha$+TIR) follows the analytical curve $1/(1-f_{esc,HII})$ with small deviations, which can be attributed to several effects. First, although $H\alpha$ emission is directly proportional to the number of LyC photons absorbed by the gas, $(1-f_{esc,HII})Q(H^0)$, this conversion depends on the hydrogen recombination coefficient that varies with physical conditions (density and electronic temperature). Second, the

<figure><img src="image_8.png" /><figcaption>Fig. 7. Predicted SFR($Q(H^0)$)/SFR($H\alpha$+TIR) vs $f_{esc,HII}$. The symbols represent the robust means of the posterior PDF and the gray contours show the $1\sigma$ uncertainties of the SUEs. The dashed line corresponds to the analytical curve $1/(1-f_{esc})$.</figcaption></figure>

$L_{TIR}$ emission that is used to derive the SFR has more complex dependences; it is sensitive to the dust content of our models and can be partially powered by X-ray photons heating the PDR. Hence, the SFR estimates we derive from $H\alpha$+TIR should, in theory, be sensitive to variations of the average physical and chemical conditions in each galaxies of our sample. Nevertheless, we find a rather narrow dispersion around the $1/(1-f_{esc,HII})$ relation that may indicate that, to first order, such effects are small. Assuming that SFR($Q(H^0)$) is the intrinsic SFR accounting for all the photons produced by the central source, the correction to apply to SFR($H\alpha$+TIR) varies from a factor 1 to ~2.5 for well constrained galaxies having at least two sectors, and reaches a factor ~15 for one poorly-constrained galaxy (31: Tol 1214-277) for which the best solution is a single-sector, completely density-bounded model (see Fig. 7).

The difference between the predicted SFR($H\alpha$+TIR) and SFR($Q(H^0)$) can therefore be explained by the presence of density-bounded regions leading to escaping photons. Indeed, photons that escape from density-bounded regions are accounted for in $Q(H^0)$ but escape before ionizing the gas that produces $H\alpha$. Nevertheless, this effect is enhanced by our model assumptions (see Sect. 7) in which photons reaching the edge of density-bounded sectors can freely escape without further interaction. In principle, part of the leaking radiation should be reabsorbed and produce a DIG contributing to the global $H\alpha$ emission (Mathis 1986; Sembach et al. 2000; Wood et al. 2010; Belfiore et al. 2022). This effect would likely reduce the correction factors that we predict here in order to obtain the intrinsic SFR($Q(H^0)$). Assuming that all the ionizing photons escaping from H II regions are reabsorbed by the surrounding DIG and assuming a maximum value for the recombination coefficient of hydrogen (corresponding to densities typically found in H II region conditions), we can provide an upper limit on the additional $H\alpha$ emission, which might be powered by leaky H II regions. In Fig. 8, we show the predictions for the maximum additional $H\alpha$ emission powered by photon leakage,

<footnote n="3">³ We apply a corrective factor of 0.945 (assuming that all the luminosity comes from stars with masses above $1 M_\odot$ and that all the mass comes from stars with masses below $8 M_\odot$) to our predicted SFR to compare them to measurements that assumed a Kroupa IMF.</footnote>

<footnote n="4">⁴ There is also a good agreement with the measured SFR(UV+24 µm). However, this comparison should be considered with caution due to the bursty nature of star formation in local metal-poor dwarfs (e.g., De Looze et al. 2014).</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Ratio of the maximum $H\alpha$ emission associated with photon leakage from H II regions, assuming that all photons are reabsorbed and of the $H\alpha$ value predicted by our models vs. $f_{esc,HII}$. We assumed $Q(H^0)$-to-$H\alpha$ conversion coefficient of $7.31 \times 10^{11}$ for a case B recombination with $T_e = 10\,000$ K and $n = 100\,\text{cm}^{-3}$ from Kennicutt et al. (1995) to estimate $H\alpha(f_{esc,HII} \times Q(H^0))$. The $H\alpha_{pred}$ values are computed consistently in the Cloudy models. The dashed line corresponds to the analytical curve $f_{esc}/(1 - f_{esc})$.</figcaption></figure>

assuming a $Q(H^0)$-to-$H\alpha$ conversion coefficient of $7.31 \times 10^{11}$ for a case B recombination with $T_e = 10\,000$ K and $n = 100\,\text{cm}^{-3}$ (Kennicutt et al. 1995) to derive $H\alpha(f_{esc,HII} \times Q(H^0))$. This quantity follows the analytical curve $f_{esc,HII}/(1 - f_{esc,HII})$. The small variations around this relation come from the fact that we assume a fixed recombination coefficient to estimate $H\alpha(f_{esc,HII} \times Q(H^0))$ while $H\alpha_{pred}$ is calculated by integrating $H\alpha$ emission at each depth, assuming a varying recombination coefficient that corresponds to the local physical conditions. We note that this upper limit was calculated with an unrealistically high recombination coefficient for DIG conditions and only yields an upper limit on the additional $H\alpha$ emission. The range of plausible values corresponding to the reabsorption of leaky LyC-photons by DIG is represented by the gray shaded area. Unsurprisingly, the large $f_{esc,HII}$ we derive translate into large upper limits of the $H\alpha$ emission powered by DIG; this maximum emission exceeds 50% of the total $H\alpha$ value for the larger $f_{esc,HII}$ we infer (i.e., $f_{esc,HII} > 50\%$). Although a significant fraction of the $H\alpha$ emission is expected to be powered by the DIG (up to ~50–60%, see e.g., Oey et al. 2007; Lacerda et al. 2018), this is not sufficient to explain the $f_{esc,HII} > 50\%$ that are found for a few objects in our sample. Most importantly, this additional $H\alpha$ component can hardly be reconciled with the fact that our predicted $H\alpha$ values (which do not account for the DIG) are in good agreement with the observed values (see Fig. 7). Possible explanations of this discrepancy are further discussed in Sect. 7.2.2.

## 6. Escape fraction of ionizing photons

We have seen in Sect. 5 that our new code is able to predict several galactic observables that are well constrained by the suite of lines given as inputs and compatible with previous

observational studies. Those observables do not depend much on the inferred topology and can be estimated in a robust way, with little influence from, for example, the number of sectors considered (see LR22). The agreement of our predictions with empirical studies shows that simple prescriptions used in previous studies (e.g., single-component and no treatment of density-bounded regions) are sufficient, to first order, to estimate metallicities and SFR. Nevertheless, several other quantities can be predicted by our code, including more complex ones for which we have only a few, sometimes biased, observational proxies, and that strongly depend on the assumed topology. This possibility is especially interesting to study escape fractions for which direct observations are challenging and indirect proxies often yield different results (see Sect. 1). We now present the results that we obtain for $f_{esc,HII}$ and discuss the observed trends and their possible physical interpretations.

### 6.1. Escape fraction vs. metallicity

In Fig. 9 (upper panel) we show the relation between $f_{esc,HII}$ (defined in Sect. 3.2.4) and metallicity. The inferred mean $f_{esc,HII}$ in our sample ranges from 0 up to around ~60%. Since we used the ionizing continuum above 13.6 eV, it also includes high energy photons that may propagate unabsorbed throughout the neutral gas. Due to such photons, the inferred escape fraction is never exactly zero, although very low for some objects. In practice, even for completely radiation-bounded models, we predict escape fractions of the order of $10^{-10}-10^{-6}$. We find that 11 galaxies are best modeled with one single-sector. In most of those cases (9 out of 11), the best solution (meaning, having the highest marginal likelihood; see Sect. 4.2) is radiation-bounded. When such single-sector, radiation-bounded models are favored, the resulting escape fraction is biased toward low $f_{esc,HII}$ values.

Density-bounded, single-sector models have been proved to not reproduce well the optical emission lines from H II regions (e.g., Stasińska et al. 2015; Ramambason et al. 2020). They are even less likely to match the observational constraints in the IR domain where some tracers arise from the PDR. For two galaxies, however, our code favors single-sector solutions that are completely density-bounded (38: HS 2352+2733 and 31: Tol 1214-277). HS 2352+2733 has only one detection ([O III]) and an upper limit on [C II], while for Tol 1214-277 most of the detections are lines arising from the ionized region and only [C II] and [O I] constrains the PDR. However, the errorbars on its [C II] and [O I] detections are rather large, hence allowing the code to favor a completely density-bounded solution since there is not enough information to constrain the neutral gas properties in those galaxies. For such galaxies, the predicted LyC-leakage can reach high values since photons escape directly from one single density-bounded sector. Thus, for all the galaxies for which the best model has only a single-sector (whether completely radiation-bounded galaxies or completely density-bounded), the $f_{esc,HII}$ estimates that we derive are relatively uncertain, with a bias toward low values for radiation-bounded models and toward high values for density-bounded models. They can, however, be considered as lower and upper limits on the escape fraction. In the rest of the analysis, we focus on a subsample of 27 out of 40 galaxies for which the best models have at least two sectors and that correspond to galaxies covered in one single pointing to derive the main trends of $f_{esc,HII}$ with galactic parameters. In practice, we plot the KDEs and correlation coefficients for this subsample of galaxies while the individual SUEs are shown for
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Escape fraction dependence on metallicity. Top: escape fraction as a function of metallicity. The symbols and ellipses represent the robust means and associated uncertainties of the MCMC draws for individual galaxies, weighted by the marginal likelihood of each configuration (see Sect. 4.3). The symbols and colors indicate the configuration having the highest marginal likelihood for each galaxy. The underlying KDE represents the distribution of all the MCMC draws for galaxies whose best model has $N_{sectors} \ge 2$ only and $\rho$ the associated Spearman correlation coefficient. Orange triangles represent galaxies for which the best model is a single component either density-bounded (DB) or radiation-bounded (RB). Their values correspond respectively to upper and lower limits of $f_{esc,HII}$. The white symbols that are flagged with a star correspond to two pointings in NGC 4214, which are excluded from the KDE. Bottom: KDE per bins of $f_{esc,HII}$ for galaxies with at least two sectors.</figcaption></figure>

the whole sample, including single-sector models and the two pointings in NGC 4214. Although the KDE accounts for the full posterior probability distribution described by all the MCMC draws for each galaxy, we note that our sample remains stastically limited and that the trends derived should be taken with caution.

The predicted $f_{esc,HII}$ spans a wide range of values, which appear to be linked to the metallicity of the galaxy. By examining the KDE of the galaxies that are best constrained, we find a global trend of an increasing $f_{esc,HII}$ with decreasing metallicity. More specifically, in Fig. 9 (upper panel), the upper right

part of the KDE is completely unpopulated meaning that high-metallicity galaxies have very low probability of reaching high values of $f_{esc,HII}$. The dispersion increases with lower metallicity. While low-metallicity galaxies do not systematically exhibit leakage, if they do, they are more likely to reach high values of $f_{esc,HII}$. In Fig. 9 (bottom panel), we show 1D-KDE plots of the metallicity by bins of $f_{esc,HII}$, for galaxies having at least 2 sectors. While this representation is quite sensitive to the intrinsic metallicity distribution of our sample, it provides a clear visual representation of the $f_{esc,HII}$-metallicity relation. Although we find rather large distributions of metallicity for most of the $f_{esc,HII}$ bins, we see that the main peaks in metallicity are shifted toward lower metallicities when the escape fraction increases. We note that the current sample is poorly populated in the low-metallicity regime (only 3 galaxies below $12 + \log(O/H) = 7.5$) and does not allow us to derive a robust trend with metallicity. Nevertheless, we find enhanced values of $f_{esc,HII}$ in the low-metallicity regime, in which all the highest predicted $f_{esc,HII}$ are found. Future observations of galaxies in the low-metallicity regime may reveal objects with small $f_{esc,HII}$. The large scatter observed among galaxies that fall within a given metallicity bin may be linked to dependences in other parameters, which will be further examined in Sect. 6.2. The metallicity-$f_{esc,HII}$ trend we present here remains clearly visible when we force the number of sectors to be the same for all galaxies (i.e., no selection of the best configuration) as shown in Fig. A.3. By looking at the individual PDFs represented by the SUEs, we confirm that the highest mean $f_{esc,HII}$ are associated with the lowest metallicities in our sample. Specifically we infer $f_{esc,HII} > 40\%$ in 5 galaxies with metallicities below 7.8. and $f_{esc,HII} > 50\%$ in 2 galaxies with metallicities below 7.3. However, galaxies with $f_{esc,HII}$ below 10% are found at all metallicities above 7.5 and are not preferentially associated with high metallicities. We also note that some objects exhibit somewhat high escape fractions on the high-metallicity part of our plot (e.g., 39: NGC 4214-c and 40: NGC 4214-s). However, those two objects correspond to separate pointings within one galaxy instead of an integrated measurement. They will be further examined in Sect. 7.

The $f_{esc,HII}$-metallicity relation presented in this section is reminiscent of the PDR covering factor vs. metallicity relation presented in Cormier et al. (2019), that show that the PDR covering factor (calculated as a common scaling factor of PDR lines) correlates with metallicity. It is also globally consistent with observational studies in which LyC-leaking galaxies are predominantly found at low-metallicity, especially for the most extreme escape fractions. (e.g., Leitet et al. 2013; Izotov et al. 2018a; Flury et al. 2022b). In the recent Low-redshift Lyman Continuum Survey, which gathers 89 LyC measurements and associated metallicity estimates, Flury et al. (2022b) find that the strong LyC-emitters with measured $f_{esc} > 10\%$ are all found at metallicities below 8.1. However, they report a significant scatter in the $f_{esc}$-metallicity relation and do not observe any significant correlation for their whole sample. We further discuss possible explanations of this scatter in Sect. 7.

Several, possibly combined, effects might lead to the $f_{esc,HII}$-metallicity trend that we observe: first, the UV opacity of the diffuse gas in the ISM is strongly dependent on metallicity. This UV opacity is controlled, in particular, by the dust content of galaxies since dust grains have been shown to absorb a significant (<50%) part of the LyC radiation in H II regions (Inoue et al. 2001). Hence, in the low-metallicity regime where the dust-to-gas mass ratio drops steeply (Rémy-Ruyer et al. 2014; Galliano et al. 2021), the ISM becomes more transparent and photons can travel larger distances. Second, the escape of
LyC photons from H II regions also depends on the properties of the host molecular clouds in which stars form. Molecular clouds are optically thick in the UV at all metallicities. Therefore, the escape fraction strongly varies with their covering factor. Recent simulations (e.g., Kimm et al. 2019; Yoo et al. 2020) have shown that the neutral gas covering factor is sensitive to both metallicit y and stellar age effects with ionizing photons escaping more efficiently from their birth cloud when the metallicity is lower. Yoo et al. (2020) emphasize that the metallicity of the cloud in which stars form seems to be an important driver of the cloud disruption timescales. They find that in the low-metallicity runs, stars break out of their birth cloud within a few megayears only while at higher metallicities they tend to stay enshrouded for a longer period.

## 6.2. Dependences on other parameters

We further examine the $f_{esc,HII}$-metallicity relation by looking at possible secondary dependences that would explain the disper sion of this relation. In particular, we focus on examining the effect of various properties of the radiation sources: stellar mass, SFR, specific SFR (sSFR = SFR/$M_*$), age of the stellar burst, and X-ray-to-stellar luminosity ratio. We revisit the relation from Fig. 9 but split our sample by bins of values for secondary param- eters (see Fig. A.6). We note that the various parameters that we examine are all interdependent. Correlations between those fundamental parameters have been extensively studied in sam- ples of local and high-redshift galaxies (e.g., Nakajima & Ouchi 2014; Duarte Puertas et al. 2022) and, more specifically, among compact star-forming galaxies (Izotov et al. 2021; Flury et al. 2022b). In particular, the stellar mass, SFR, and sSFR are func- tions of the metallicity. Hence, we find that the highest $f_{esc,HII}$ that we infer are associated with low-metallicity, low stellar masses, and high sSFR. We find that those three parameters are the main drivers of large $f_{esc,HII}$ predicted in our sample. How- ever, our analysis is limited by the observational bias of our sam- ple in which only highly star-forming galaxies are populating the low-metallicity and low mass regime. With all those correlated parameters, it is not possible to infer which one has the strongest effect on $f_{esc,HII}$. We also find that stellar populations with ages above 3 Myr and an X-ray-to-stellar luminosity ratio above 0.5% (−2.3 in log) seem to be associated with enhanced escape frac- tions (see panels 4 and 5 from Fig. A.6). These results motivate the study of the $f_{esc,HII}$ dependences on galactic parameters other than metallicity, which we now examine in more details.

### 6.2.1. Stellar mass

It has been suggested that galaxies with shallow gravitational wells, associated with low stellar masses, would be good can- didates for potential LyC-leakage. In particular, this $f_{esc} - M_*$ trend has been observed in samples with masses down to $10^8 M_\odot$ (Leitet et al. 2013; Izotov et al. 2018b). However, a recent fol- low up study carried out in Izotov et al. (2021) targeting galax- ies with masses below $10^8 M_\odot$ finds no clear correlation between the global measured LyC escape fraction and the stellar masses.

In Fig. 10 (panel a), we find that large $f_{esc,HII}$ tend to be asso- ciated with low stellar masses, although with a large scatter. All the galaxies with masses above $10^9 M_\odot$ have mean escape frac- tions below 40%. On the other hand, galaxies with masses below $10^9 M_\odot$ do not always have high escape fractions. We note that this effect might be dominated by the mass-metallicity distribu- tion of our sample in which lower masses are associated with the lowest metallicity. The correlation between $f_{esc,HII}$ and metallic-

ity is, nevertheless, much clearer. Globally we observe a weak anti-correlation between $f_{esc,HII}$ and $M_*$, although with a large scatter, especially in the low-mass regime.

### 6.2.2. Star-formation rates

Recent simulations of photon leakage have shown that the process is a feedback driven mechanism (Trebitsch et al. 2017; Kimm et al. 2017, 2019; Kim et al. 2018, 2019, 2021; Kakiichi & Gronke 2021; Yoo et al. 2020). Although the main driver, whether stellar or supernovae (SN) feedback, remains elusive, the escape fraction is predicted to increase drastically when the star-formation efficiency is enhanced.

Surprisingly, we find no correlation between $f_{esc,HII}$ and SFR (see panel b from Fig. 10). The highest values of $f_{esc,HII}$ are asso- ciated with moderately star-forming galaxies in our sample. This absence of trend can be understood by looking at the KDE in Fig. A.6 (first panel); the highest SFRs in our sample correspond to the highest metallicities, which do not allow high amounts of photon leakage. On the other hand, we find that $f_{esc,HII}$ increases with increasing sSFR. We show on Fig. 10 the evolution of $f_{esc,HII}$ with respect to measured sSFR (panel c) and estimated sSFR from $Q(H^0)$ (panel f) as defined in Sect. 5.1. Although the trend is already visible when using the sSFR measurements from the literature, we find, not surprisingly, that applying a correction to account for the photon leakage through density-bounded sec- tors increases the dynamical range of the relation. Similarly to the $M_*$-$f_{esc,HII}$ relation and SFR-$f_{esc,HII}$ relation, the sSFR-$f_{esc,HII}$ relation is possibly driven by metallicity effects as the lower- metallicity sources are those with higher sSFR.

### 6.2.3. Stellar ages

Stellar evolution and the feedback mechanisms responsible for clouds dispersal may also impact the escape fraction. In our code, we only model very young bursts of ages below 10 Myr that represent the bulk of the H II regions that dominate the emis- sion. We expect both rotation and binary stars to affect our pre- diction for the H II region escape fractions by producing more LyC photons for a longer period of time (Choi et al. 2020, and reference therein). In particular, Choi et al. (2020) finds that at low-metallicity ($<1/3 Z_\odot$) models that include binary stars pro- duce 60% more ionizing flux compare to model without.

In Fig. 10 (panel d), we see that most of the stellar popula- tions in our sample are dominated by young starbursts with ages of 2 to 6 Myr. These timescales suggest that leakage might occur at very early stage of stellar evolution. We find no clear trend of $f_{esc,HII}$ with stellar age, although large $f_{esc,HII}$ above ~50% seem to be associated with ages greater than 3–4 Myr. We emphasize, however, that large ages are not necessarily associated with leak- ing galaxies. As previously mentioned in Sect. 5.1, we note that the age and metallicity are interdependent, with larger ages tend- ing to be associated with the lowest metallicities.

While our findings are consistent with scenarios in which photons might escape at very young ages (e.g., Kakiichi & Gronke 2021), we emphasize that our modeling assumptions based on a single burst stellar population prevent us from examining the effect of an older underlying population. Other sources of feedback (e.g., turbulence, SN feedback) might also play a role at later stellar ages. While at early ages, the feed- back is dominated by pressure and ionization mechanisms from young H II regions, older stars can also contribute through stel- lar winds, especially in the Wolf-Rayet phase. We note that this
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Escape fraction dependences on measured and inferred galactic parameters. **Top panel:** escape fraction vs. some measured galactic parameters: $M_*$, SFR and measured sSFR. **Bottom panel:** escape fraction vs. some inferred galactic parameters: stellar age, $L_x$, sSFR($Q(H^0)$). The underlying KDE represents the distribution of all the MCMC draws for galaxies with $n_{sectors} \ge 2$ only and $\rho$ the associated Spearman correlation coefficient. Orange triangles represent galaxies for which the best model is a single component either density-bounded (DB) or radiation-bounded (RB). The white symbols that are flagged with a star correspond to two pointings in NGC 4214, which are excluded from the KDE.</figcaption></figure>

stellar evolution also depends on metallicity effects with low-metallicity stars producing more ionizing photons over larger timescales (Eldridge et al. 2017; Xiao et al. 2018). Finally, older generation of stars might be responsible for creating low-density channels in the ISM, either due to mechanical feedback of the starburst (e.g., Ma et al. 2020) or through SN explosions. The presence of such low-density channels carved by previous generations of stars might favor leakage associated with young bursts that formed in regions where gas has already been partially cleared out (e.g., Ma et al. 2020). Nevertheless, it is worth noting that the stellar ages we derive are consistent with the disruption timescales found in simulations from Yoo et al. (2020). They are also very close to the feedback timescales of massive O/B stars measured in resolved giant molecular clouds (GMCs) from nearby spiral galaxies (e.g., Chevance et al. 2020a, 2022; Kim et al. 2022). Chevance et al. (2022) also point out that the efficiency at which the injected energy couples with the GMCs is low (a few tens of percent) and conclude that most of the energy and momentum produced by young stars escape outside of the host GMC, which is consistent with the high $f_{esc,HII}$ we derive.

### 6.2.4. X-ray sources

As discussed in Sect. 2.1, several bright X-ray sources have been observed in our sample. Due to our simple modeling assumption and the suite of lines used in this study, the X-ray luminosities we derive are not tightly constrained. One of the motivations to

include additional feedback mechanisms is the presence in our sample of extreme emission lines emitted by species with very high ionization potentials (>54 eV): [O IV] or [Ne V]. Such lines are, nevertheless, not numerous and not available for all our sample (see Fig. A.4). As a consequence, for most of the galaxies, the $1\sigma$ SUEs are large and the X-ray luminosity is not tightly constrained.

We note, however, that the X-ray-to-stellar luminosity ratio is also constrained by emission-lines in the neutral gas. This is due to the penetration of X-ray photons deep into the neutral and molecular gas because of their small photoionization cross-sections. Such photons are able to deposit their energy at large $A_V$, beyond the ionization front of H II regions. Because of this effect, our code is able to place an upper limit on the X-ray-to-stellar luminosity ratio based only on the emission lines arising from the neutral gas (e.g., [Fe II], [Si II], [O I]), even when no high ionization lines are detected. In particular, photoelectric heating on dust grains contributes to the heating of neutral gas in our models and is controlled by the assumed dust-to-gas mass ratio (fixed in our models) and total stellar luminosity (constrained by H II region emission lines). Hence, the additional X-ray contribution to neutral gas heating is limited by the emission lines in the neutral gas.

Although we derive uncertainties on the absolute log $L_X$ above 1 dex for most galaxies, we predict a relatively well-constrained X-ray luminosity for a few galaxies. Specifically, we obtain a total X-ray luminosity of log $L_X = 40.33^{+0.53}_{-0.49} \text{ erg s}^{-1}$
for IZw 18⁵. The value obtained for IZw 18 is slightly higher but compatible with the measurement of $1 \times 10^{40}$ erg s⁻¹ (intrinsic X-ray luminosity; 0.3–10 keV) reported in Kaaret & Feng (2013). In our case, the predictions are for the total X-ray luminosity integrated over the whole energy range covered by the multicolor blackbody assumed as an input. While the X-ray-to-stellar luminosity ratio may allow us to pinpoint some of the most extreme escape fractions in our sample, we note, however, that no global trend seems to be visible, with some galaxies having little to no contribution coming from X-ray being associated with large escape fractions.

### 6.3. Tracers of density-bounded regions and escaping photons

Our code allows us to disentangle the emission of tracers that emit throughout different phases of the ISM. This is in particular the case of the [C II]158 µm line, which is emitted both in the ionized and neutral gas phase. This impacts the [O III]88 µm/[C II]158 µm (O3C2) diagnostic as a tracer of $f_{esc}$, which is quite popular in high-redshift studies (e.g., Harikane et al. 2020; Katz et al. 2022b), as it can be detected with ALMA up to redshifts above 6. This O3C2 ratio has been extensively studied with results indicating that high O3C2 values could be explained by specific physical conditions such as a high burstiness parameter (Vallini et al. 2021) or a low covering factor of PDR (Harikane et al. 2020). Such a trend has also been reported in Chevance et al. (2016, 2020b) that find that O3C2 is sensitive to the PDR filling factor defined for each pixel in their resolved analysis of 30Dor in the LMC. Nevertheless, at galactic scale, Cormier et al. (2019) find little correlation of O3C2 with their estimated PDR covering factors in the DGS. We now discuss the complex origin of the [C II] line and how it might affect the predictive power of the O3C2 line ratio. We then discuss the potential use of other line ratios to trace the presence of density-bounded regions and of associated LyC-photon leakage from H II regions, both in the IR and optical domain.

#### 6.3.1. The complex origin of [C II] 158 µm emission

The [C II]158 µm line is an important diagnostic, which constrains in principle the neutral gas reservoir and the potential presence of density-bounded regions. However, the interpretation of this tracer is challenging as [C II] may be produced both throughout the neutral and ionized regions (e.g., Oberst et al. 2011; Ramos Padilla et al. 2021, 2022). In this section, we intend to disentangle the [C II] contribution from different phases by examining predictions for the fraction of [C II] luminosity arising from the ionized phase of the ISM ($f[\text{C II}]_{H+}$).

To do so, we weigh the cumulative [C II] line luminosity by the relative abundance profile of the hydrogen in each phase. For the ionized phase, we hence scale the cumulative line luminosity by the $H^+$ relative abundance profile and we integrate over the depth for each sector. The total [C II] luminosity arising from the ionized phase is obtained by summing over all the sectors of the model. We then derive $f[\text{C II}]_{H+}$ by computing the ratio of the prediction for [C II] in the ionized phase over the prediction for [C II] total luminosity.

In Fig. 11, we note that for one of the 1-sector, density-bounded galaxies (31: HS 2352+2733), $\log f[\text{C II}]_{H+} = 0$ since the gas in those models is completely ionized. For the other galaxies, $\log f[\text{C II}]_{H+}$ varies between -3 to -0.5 and its value correlates with the metallicity of the galaxy. This confirms previous

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Fraction of [C II]158 µm in the ionized phase vs. metallicity. The single-sector, density-bounded models represented by orange stars pointing down are PDR-deprived; all their [C II] emission arises from the ionized gas ($\log f[\text{C II}]_{H+} \sim 0$). The underlying KDE represents the distribution of all the MCMC draws for galaxies with $N_{\text{sectors}} \ge 2$ only and $\rho$ the associated Spearman correlation coefficient. The dashed line represents the relation derived by Cormier et al. (2019) for the DGS, assuming a fixed electron density $n_e = 30 \text{ cm}^{-3}$.</figcaption></figure>

results from Croxall et al. (2017) and Cormier et al. (2019) that show that the fraction of [C II] arising from the ionized phase decreases at low metallicities. We also find values in the same dynamical range as previously found in Cormier et al. (2019) but with a somewhat smaller dispersion. According to Croxall et al. (2017), this result might be linked to the presence of harder radiation fields in low-metallicity environments. As a result, the relative abundance of $C^{2+}$ in the ionized gas may increase with respect to the $C^+$ abundance, leading the [C II] emission to be dominated by the neutral phases. We discuss the implications on the use of [C II] as a potential tracer of the escape fraction in Sects. 6.3.2 and 7.

#### 6.3.2. IR line ratios

We now provide some observed or unobserved tracers of the escape fraction based on what is predicted by our multisector models. The MCMC algorithm used in our code (see Sect. 4.2) does not only sample the PDFs of primary parameters but can also infer the PDFs of secondary parameters and unobserved tracers. In particular, we can extract predictions for unobserved line fluxes and the associated uncertainties that we represent using the same formalism as for primary parameters (i.e., SUEs and KDE, see Sect. 4.3).

Line emission produced by ions with different ionization potentials allow us to trace different physical depths. Hence, ratios involving lines emitted closer to the stars, in highly excited gas and lines emitted further away in regions near or past the ionization front provide an efficient way to probe the structure of the gas around a cluster of stars. This method can help identify galaxies in which numerous density-bounded regions

<footnote n="5">⁵ The interval corresponds to the highest density interval at 94%.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** KDE showing $f_{esc,HII}$ as a function of IR ratios. The KDE represent the distribution of all the MCMC draws for galaxies with $N \ge 2$ only and $\rho$ the associated Spearman correlation coefficient. For this, we group the MCMC draws for our selected sample by bins of $f_{esc,HII}$ and fit each bin with a gaussian $N(x_0, \sigma^2)$. The blue dots and dashed lines respectively represent the gaussian mean $x_0$ and the gaussian width $\sigma$. The plain blue line shows a weighted linear regression applied to the gaussian centers, accounting for the widths of each bin. The corresponding fit values are provided in Table A.2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Individual SUEs showing $f_{esc,HII}$ as a function of IR ratios. We note that Tol 1214-277 (31) most often lies outside of the range shown in the x-axis. This is due to the lack of constraints on this galaxy for which the best model (one sector, completely density-bounded) predicts unrealistically low emission from ions with low ionisation potential.</figcaption></figure>

are present. We expect such regions, which are bounded by the lack of gas, to have a reduced emission of low ionization lines emitted further away from the stars. We show in Figs. 12–15 the line fluxes predicted by our models for all lines, regardless of whether or not they are detected. We note, however, than when a line is detected and used as a constraint, the prediction should be close to the observed value.

In Figs. 12–13 we show how classical ratios of IR lines relate with $f_{esc,HII}$. In particular, we examine the O3C2 ratio, which involves the [C II]158 µm whose complex origin was discussed in the previous subsection. In the framework of our multisector modeling, we do find a correlation between O3C2 and $f_{esc,HII}$. Although this correlation is not clearly seen by looking at individual SUEs, it is visible on the KDE representing the whole sample (i.e., the concatenated MCMC draws for all the galaxies with at least 2 sectors). We emphasize, however, that for most galaxies the dynamical range of the O3C2 in our sample is rather narrow (<2 dex), which do not allow us to easily identify leakers based on their O3C2 ratio.

Although the [N II] lines are fainter and more challenging to detect, we find that [O III]88 µm/[N II]122 µm and [O III]88 µm/[N II]205 µm (O3N2) are better tracers of $f_{esc,HII}$

both because of the tightness of the relation and of their wider dynamical range (over 3 orders of magnitudes), which makes them potentially interesting diagnostics for observation-based classification. These correlations might be explained by the fact that both [O III] and [N II] originate from H II regions but are emitted by ions with different ionization potentials (35.1 eV for $O^{2+}$ and 14.5 eV for $N^+$). This means that [O III] is predominantly emitted in a highly irradiated region close to the stars while [N II] arises from the outer part of the H II region, and hence their ratio depends on the depth of each sector. On the other hand, [C II] has a complex origin and arises both from H II regions and PDRs. As shown in Fig. 11, only a fraction of [C II] arises from the ionized gas. As one goes deeper in the neutral gas phase beyond the ionization front, the [C II] emission increases while $f_{esc,HII}$ remains unchanged. Hence, in low-metallicity environments in which [C II] is predominantly emitted in the neutral and molecular phases, the O3N2 ratio could be a better tracer of density-bounded H II regions as it is less contaminated by PDR emission than O3C2. We emphasize that, despite their interest for observational studies, the ratios presented here are strongly sensitive to the abundance profiles (C/O and N/O) assumed in this study. Other tracers, including ratios of the same
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** KDE showing $f_{\text{esc,HII}}$ as a function of line ratio involving ions with high ionization potentials.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Individual SUEs showing $f_{\text{esc,HII}}$ as a function of line ratio involving ions with high ionization potentials.</figcaption></figure>

element that are not sensitive to abundance profiles are provided in Table A.2 where we report the best potential lines ratio having a Spearman coefficient above 0.3.

In Figs. 14–15, we examine some ratios involving the lines emitted by ions with high ionization potentials: [Ne III]15 µm, [O IV]25 µm and [Ne V]14 µm. The luminosity of these lines strongly depends on the properties of the potential X-ray source and can span a large range of values. Preliminary tests performed with single-star population have shown that our code always needed an additional contribution of an X-ray component to match the observed values of such highly ionized species (as well as neutral gas lines). In the current study, which includes a contribution from binary stars with BPASS, an additional contribution from X-ray is not systematically needed and the X-ray-to-stellar luminosity ratios span several orders of magnitudes (see panel e from Fig. 10). We also note that typical ages found by MULTIGRIS are shifted to later ages when including effects from binary stars.

We find a correlation between $f_{\text{esc,HII}}$ and [Ne III]15 µm/ [Ne II]12 µm (Ne32), [O IV]25 µm/[C II]158 µm (O4C2) and [Ne V]14 µm/[Ne II]12 µm (Ne52). Nevertheless, we find that O4C2 and Ne52 cover much larger dynamical ranges. In Fig. 15, the large uncertainties of the ellipses reflect the uncertainty on the X-ray component, which is poorly constrained when no [O IV] and [Ne V] detections are provided. Although such ratios involving lines emitted by ions with high ionization potentials are subject to large uncertainties due to the unknown X-ray

source properties, they could become of high interest if more constraints on the X-ray spectrum are available. In particular, a more complete statistical sample would provide interesting insights on whether or not leakage is linked to the presence of an X-ray source. Compared to our current predictions, it would also help pinpoint some of the missing physics in our models and might provide insights on the nature of X-ray sources.

### 6.3.3. Optical line ratios

Although our study focuses on IR tracers, we also examined a few well-known optical and mixed tracers that have been discussed in the literature as possible tracers of LyC escape fractions (see Figs. 16–17). We find that the $L_{\text{TIR}}/L_{\text{bol}}$ ratio anti-correlates tightly with $f_{\text{esc,HII}}$. While $L_{\text{bol}}$ accounts for the total luminosity produced by the central stellar cluster (and the additional X-ray contribution), $L_{\text{TIR}}$ traces the photons that have been reprocessed by the surrounding gas and dust and re-emitted in the IR. Thus, this ratio traces well the escaping energetic photons that are produced by the stellar cluster but escape before being reprocessed by the gas. We note, however, that this method either requires modeling assumptions for the stellar population or fitting the continuum to derive $L_{\text{bol}}$.

Alternatively, we also put to test the [O III] $\lambda\lambda4959,5007\AA$/[O II] $\lambda\lambda3726, 3728\AA$ (O32) and [S III] $\lambda\lambda9068,9532\AA$/[S II] $\lambda\lambda6716,6731\AA$ (S32) optical line ratios. The O32 ratio has been claimed to efficiently pick up
Such tracers are mostly accessible in the IR domain. (e.g., [O I], [C II], [Si II]).

Previous IR studies have focused on quantifying a PDR covering factor (Cormier et al. 2012, 2019; Harikane et al. 2020), which has been used to describe the apparent porosity to ionizing radiation at galactic scales. This PDR covering factor was calculated using models where all the lines arising from the PDR were simultaneously scaled by a common factor between 0 and 1. This grouping is somewhat arbitrary as some lines emit throughout both H II regions and the PDRs. Additionally, these models did not account for possible density-bounded regions since the PDR zone (either total or partial) was always combined with a fully radiation-bounded H II region. Hence, the interpretation of the PDR covering factors is difficult as the simplistic geometry does not account well for complex gas distribution and cannot be linked to potentially escaping photons. In particular, this quantity is very sensitive to modeling assumptions such as the stopping criteria for models (e.g., depth or $A_V$). Among other factors, accounting for a nonunity filling factor of the ISM or for the physical depth of the photodissociated layers changes the values derived for the PDR covering factor: a very thin PDR with a large covering factor can produce the same emission lines as a thick PDR with a lower one.

In our topological models, we relax the assumption of a full PDR model scaled by a covering factor in favor of an agnostic combination of sectors (as described in Sect. 3) where both the weights and the depths of each sector are free parameters. The Bayesian framework described in Sect. 4.2 allows us to derive PDF and credibility intervals for secondary parameters and to calculate, for each sector, the amount of escaping photons. We then derive a global averaged $f_{esc,HII}$ (see Sect. 3.2.4) corresponding to the population of H II regions present in a galaxy. This integrated quantity directly relates to the ISM gas distribution as opposed to the previously used PDR covering factor.

Nevertheless, as mentioned in Sect. 4.1, our approach bears its own caveats. In particular, the interdependence between total cluster luminosity and escape fractions is not easy to solve. Nevertheless, although scaling the total cluster luminosity affects the absolute lines fluxes, it does not modify the line ratios that are tracing density-bounded regions (see Sect. 6.3). Constraining the cluster luminosity through $L_{TIR}$ is a crucial step that allows us to lift this degeneracy. Still, since we do not use any continuum emission, the UV luminosity we derive for our model depends on our modeling assumptions for the stellar population. The values we derive for the total bolometric luminosities and the escape fractions incorporate uncertainties on the stellar population. Another limit of our work comes from the assumed topology of galaxies; although the multisec-tor combination introduces a layer of complexity that is usually unaccounted for when using single-component photoionization models, this representation remains simplistic when compared to the actual ISM morphology.

## 7.2. From H II regions to galaxies

The H II region escape fractions that we derive in this study are hardly comparable with the direct UV detection of LyC. Such detections are, first of all, very challenging at low-redshift because the sensitivity of current observatories (e.g., HST/COS) drops at short wavelengths. The few detections available are remarkably small in the local universe and only one galaxy in our sample has a LyC detection (Haro 11, $f_{esc}$ = 3.2%). We note that this value is compatible with the prediction we derive

for Haro 11 of $f_{esc,HII}$ = 9.59$^{+9.80}_{-9.51}$%$^6$ Similarly, other observations of a few galaxies in the Local Universe detect little to no LyC leakage (e.g., Bergvall et al. 2006; Leitet et al. 2013; Borthakur et al. 2014; Leitherer et al. 2016). Those observations remain very sparse and suffer uncertainties due to line-of-sight variations. They are hard to reconcile with predictions from models (Fujita et al. 2003) and resolved observations of local star-bursts, which suggest that they should be leaking a substantial amount of photons. Quantifying the amount of ionizing photons that are actually escaping at galactic scales remains an elusive question to which we cannot answer within the scope of this study.

### 7.2.1. Resolution effects

Due to the simplified topology we consider (i.e., a galaxy being a weighted sum of star-forming regions), our predictions for $f_{esc,HII}$ represent an averaged value of the number of photons escaping from the H II regions that populate one galaxy. Although this observable is not easily linked to direct galactic-scale observations, it is representative of the ISM structure within star-forming regions and provides a metric to quantify how many density-bounded regions contribute to the integrated emission. This approach is an attempt to generalize previous analysis performed at local scales with resolved studies in which escape fractions could directly be mapped at the scale of H II regions.

Focusing on resolved H II regions in the local starburst galaxy IC10, Polles et al. (2019) introduced a method to take into account density-bounded regions by considering the cloud depth as a free parameter. While their approach is limited to single-sector models, they find that the observed spatial scale influences the predicted depth, with regions being almost all density-bounded at small scales (~200pc) but shifting to a radiation-bounded regime at large scales. At even smaller scales, hints of substantial escape fractions (above 80%) have been reported in Chevance et al. (2022) where they find very low feedback coupling efficiencies in GMCs from 9 nearby star-forming galaxies, meaning that most of the feedback energy injected by young stars is not dissipated and may escape.

Another technique based on ionization parameter mapping (IPM) using the [O III] $\lambda\lambda4959,5007\AA$/$[S II] \lambda\lambda6716,6731\AA$ line ratio was introduced by Pellegrini et al. (2012) to study the optical depth of H II regions in the LMC and SMC. This method has allowed them to classify the population of H II regions of both galaxies and derive rough estimates of the average luminosity weighted $f_{esc,HII}$ (≥40% in both galaxies). These estimates were then corrected using DIG measurements to derive an global average escape fraction of 4–9% in the SMC and 11–17% in the LMC. A similar technique using IPM based on S32 was used on local starbursts (including three galaxies in our sample) by Zastrow et al. (2011, 2013) to identify ionizing cones and gas features that could be associated with enhanced escape fractions. Despite the very good resolution of their data (~60 to 200 pc arcsec$^{-1}$), the analysis could not provide any quantitative estimate of the escape fractions since individual H II regions were not resolved.

Recent studies have been able to overcome this resolution limit and have provided quantitative estimates of the escape fraction for a few local galaxies. Using ~83 000 resolved stars, Choi et al. (2020) find that the escape fractions of star-forming regions within NGC 4214 exhibit a wide range of

<footnote>$^6$ The interval corresponds to the highest density interval at 94%.</footnote>
values going up to 40%. Their results yield a global escape fraction of about ~25% for the whole galaxy. Although our study focuses on two different regions of this galaxy, the values we derive (34.67<sup>+12.29</sup><sub>-13.27</sub>% for NGC 4214-c and 32.55<sup>+16.44</sup><sub>-15.40</sub>% for NGC 4214-s) are in good agreement with their predictions. Another approach based on H II region mapping was also made possible by the resolution of MUSE on the ELT. Weilbacher et al. (2018) mapped 386 H II regions from the central field of the Antennae nebula (having approximately a solar metallicity; Bastian et al. 2006; Lardo et al. 2015) and found 38 of them to exhibit large $f_{esc,HI}$ with a mean value of 72%. They estimate the overall escape fraction of the galaxy to be around 7%. Recently, Della Bruna et al. (2021) found an average escape fraction of 67% for the population of H II regions observed with MUSE in NGC 7793 (which has an averaged metallicity close to the solar value; Pilyugin et al. 2014). They find that the radiation leaking from H II regions is sufficient to explain by itself the measured fraction of DIG in this galaxy.

Our results corroborate the picture drawn by these recent studies that find that a substantial amount of ionizing photons might be leaking from H II regions. Connecting these escaping photons to the surrounding diffuse gas is a key element to understand whether or not photons might be leaking out at larger scales.

### 7.2.2. Missing diffuse component

In our models, a galaxy is represented as a combination of H II regions connected to PDR and molecular zone (see Fig. 3). Additionally, the potential presence of DIG is accounted for by allowing the code to add a sector with very low hydrogen density (as low as $1 \text{ cm}^{-3}$). This component is only added by our code if it is needed to reproduce the suite of constraints. In our solutions, however, such diffuse sectors are usually absent. Instead, the code favors configurations with hydrogen density typical of H II regions (around $100 \text{ cm}^{-3}$). This is understandable as the IR tracers that we use to constrain our models predominantly trace H II regions and PDRs. Additionally, our treatment of this DIG component is limited by the independent radiative transfers within each Cloudy model considered as a sector. Although the photon transfer is consistently computed throughout each phase, the different sectors are independent from one another and do not overlap, meaning that photons escaping from a given sector cannot be reabsorbed by a second component along the line of sight.

This unaccounted for DIG could be problematic to reproduce the optical emission lines (e.g., $H\alpha$ in Sect. 5 and optical line ratios in Sect. 6.3) we derive in our analysis. Indeed, contamination by DIG has been extensively studied in the optical domain and is known to significantly contribute to some line ratios such as $[S\II]/H\alpha$, $[N\II]/H\alpha$ and $[O\ III]/H\beta$ (Oey et al. 2007; Zhang et al. 2017; Kreckel et al. 2016; Vale Asari et al. 2019, 2020; Espinosa-Ponce et al. 2020; Della Bruna et al. 2021; Belfiore et al. 2022). Nevertheless, we have seen in Sect. 5.2 that our models predict $H\alpha$ fluxes in agreement with observations. This agreement is expected since our code reproduces both the metallicity (see Fig. 5) and the recombination line $H\alpha$ in the galaxies in which it is detected. Nevertheless, it appears at odds with the large $f_{esc,HI}$ we derive and could be explained by different scenarios. A first possibility is that the measured $H\alpha$ emission in the DGS galaxies is dominated by H II regions with a negligible contribution from DIG. This would be in line with observational studies that find that, as opposed to nonstarbursting galaxies, the $H\alpha$ emission in star-

bursts is largely dominated by H II regions in which the emission is boosted (Hanish et al. 2010). In that case, both our models and $H\alpha$ observations would be fairly independent of DIG. This would suggest that the set of observables used in this study is not sensitive to DIG and that additional tracers are needed to probe this diffuse component.

Alternatively, if part of the detected $H\alpha$ emission is powered by DIG, the fact that our models reproduce it without including a diffuse component means that our code somehow compensates for the missing DIG contribution, although assuming a simplified DIG-free geometry. As opposed to the previous scenario, this would only be possible if some of the IR tracers we use are sensitive to the DIG contribution. Although, to our knowledge, no systematic studies are available regarding DIG contamination in the IR domain, we do have some observational evidence that it may marginally contribute to the emission of some of the tracers considered in our study (Cormier et al. 2012, 2019). All in all, we find that our models based on combination of H II regions are sufficient to reproduce the observed $H\alpha$ emission, although they may be either missing or compensating for the DIG contribution. To disentangle those scenarios, further progress is needed to either correct the emission lines for DIG contamination before using them as constraints of MULTIGRIS or properly including in the models an additional DIG component. To that purpose, high spatial-resolution and sensitive telescopes in the far-IR would be useful to study the properties of the DIG.

Most importantly, we also do not account for the actual distribution of diffuse neutral gas. This is a major caveat considering that escaping photons has been shown to be very sensitive to the distribution of neutral hydrogen (e.g., Gazagnes et al. 2020). Although our models make the most of the few IR constraints tracing the PDRs, we are only accounting for the dense, highly irradiated neutral gas, which is seen in emission. This might partly explain why our estimates of the total H I masses (M(H I)), shown in Fig. A.2, tend to be significantly smaller than the measurements from Rémy-Ruyer et al. (2014) even though we find that those values agree within 0.5dex for 16 out of 40 galaxies. However, comparing the H I column densities of our sectors, we see that even the deepest sector of each model has a column density much lower than that measured in absorption (see Table A.3). This is one of the main limitations of our models: even when our models successfully recover the total H I gas mass, the distribution of sectors remains too simplistic. The actual distribution of the ISM likely resembles more the randomly distributed version presented in Fig. 3. In particular, introducing holes or very low column-density sight-lines in our topological models would naturally push the code to predict deeper sectors with higher column densities. Because of the aforementioned limitations regarding the diffuse ionized and neutral gas, we cannot provide estimates for the global galactic escape fraction in our sample. Instead, we focus on analysing the trends with various parameters and their physical interpretation.

## 7.3. Dependences of $f_{esc,HI}$

### 7.3.1. Correlation with galactic properties and observational limitations

We find a clear anti-correlation of $f_{esc,HI}$ with metallicity and a positive correlation with sSFR. Such trends are reminiscent of those observed in the local sample of few direct measurements from Leitet et al. (2013). They are, however, not clearly seen in larger surveys using UV detection of the LyC (e.g., Saxena et al. 2022; Pahl et al. 2022), although the recently
observed Low-redshift Lyman Continuum Survey (89 galaxies with LyC measurements) has revealed a similar trend of $f_{esc}$ with sSFR (Flury et al. 2022b).

The lack of observed correlation between galactic properties and measured LyC escape fractions might be due to a dilution effect of the localized H II regions from where photons escape into the galactic ISM (Saxena et al. 2022). This dilution may explain why we find clear trends of $f_{esc,HII}$ with several galactic parameters (see Sect. 6.2) while such correlations are unseen at galactic scales. Another crucial parameter that might strongly limit the observation of empirical correlation is the unknown inclination of unresolved galaxies. For example, Bassett et al. (2019) emphasize that different viewing angles can affect both the measured LyC values and the observed optical lines ratios (such as O32). Under the hypothesis of an anisotropic LyC leakage, Nakajima et al. (2020) suggest that the viewing angle of porous H II regions might be the main parameter controlling the detection or nondetection of LyC leakage in a sample of Lyα emitting galaxies from the LymAn Continuum Escape Survey (LACES). Although the dependence of emission line ratios on the viewing angle is somewhat mitigated by the use of optically thin tracers (e.g., IR lines), the exact impact of inclination on the observed emission from galaxies remains to be quantitatively examined. This effect was not examined in the current study since we cannot extract any information on the inclination in our sample of galaxies.

Further progress on this issue requires the use of complex 3D radiative transfer codes as this effect is not accounted for when using 1D radiative transfer codes such as Cloudy. For example, Mauerhofer et al. (2021) performed a detailed investigation of the viewing angle impact on simulated galaxies and find that the measured projected escape fractions are strongly impacted by orientation (with variations from 0 to 47%). This lack of correlation might also be due to propagation effects through the DIG, which are not taken into account in our H II region modeling. Those effects remain largely unconstrained since they depend on the H I gas 3D-distribution, which is unknown.

### 7.3.2. Potential impact from nonstellar feedback

Although the high $f_{esc,HII}$ in our sample are predominantly associated with low-metallicity, low mass, and high sSFR, a few objects in our sample exhibit surprisingly high $f_{esc,HII}$ with somewhat larger metallicities and moderate sSFR. We now discuss additional mechanisms that could explain the position of such objects on the $f_{esc,HII}$-metallicity relation. In particular, we examine if other feedback mechanisms such as the potential presence of an ULX or outflows could be responsible for the position of outliers in our sample. In Fig. 18 (first panel), we mark the galaxies having a high measured sSFR (log sSFR > −9.5). This criterion is common to all the galaxies with inferred $f_{esc,HII} > 40\%$ and shared by 10 out of 14 galaxies with $f_{esc,HII} > 20\%$. Nevertheless, 4 galaxies with $f_{esc,HII} > 20\%$ exhibit only a moderate sSFR. We now examine other potential mechanisms that may be responsible for their enhanced $f_{esc,HII}$.

On the middle panel from Fig. 18, we examine the impact of the X-ray-to-stellar luminosity ratio. Although it has not been shown that the presence of an X-ray source have a direct impact on escape fractions, high $L_X$ tends to be associated with of intense star-formation rates (Grimm et al. 2003) and the emission of high-energy photons, which travel deep into the ISM may enhance the global escape fraction of ionizing photons. We

find X-ray-to-stellar luminosity ratio above 0.5% in 8 out of 40 sources. Two out of the five reported potential ULX in our sample are associated with $L_X/L_* > 0.5\%$. We note that the two galaxies (13: IZw 18 and 25: SBS 0335-052) with the highest $f_{esc,HII}$ in our sample ($>50\%$) are associated with a known ULX. However, similar levels of X-ray contribution are also found in galaxies for which we infer little to no leakage. All in all, we find that large X-ray-to-stellar luminosities are not necessarily associated with enhanced leakage, although X-ray sources may contribute in some specific cases. Additionally, we note that a recent study from Marques-Chaves et al. (2022) finds no correlation between the LyC escape fractions measured in the Low-z Lyman Continuum Survey and the spectral hardness of the ionizing sources.

Finally, we expect dynamical effects disturbing the gas kinematics to also impact the H II region escape fractions. Outflows may remove gas from the galactic disk and create low-density channel from where photons can freely travel (Fujita et al. 2003). Recent studies have also examined how gas fragmentation induced by turbulence might play a role in producing the patchy structure of the ISM through which photons escape (e.g., Kakiichi & Gronke 2021; Hogarth et al. 2020). Although outflow signatures are only detected in a few objects (11/40), they offer a plausible explanation to interpret the position of some outliers having large predicted $f_{esc,HII}$ but modest sSFR and little to no contribution from an additional X-ray source. Detailed studies of the gas kinematics are available for a few objects in our sample (see Sect. 2.1 for references) based on either resolved observations or analysis of emission line profiles. Since photoionization and photodissociation models do not account for dynamical effects, we only provide a qualitative discussion on the potential role of gas outflows. In Fig. 18 (third panel), we show that the presence of outflows may explain the position of several galaxies with relatively high metallicities ($12+\log(O/H)>8$) but $f_{esc,HII} > 20\%$ (e.g., 2: Haro 3, 23: NGC 625, and 39–40: NGC 4214).

In summary, we find that although the presence of nonstellar feedback mechanisms such as a strong contribution of an X-ray source or outflowing gas may contribute to power photon leakage, such mechanisms seem to only have a secondary influence on the inferred $f_{esc,HII}$. Instead, we find that variations of $f_{esc,HII}$ appear to be mainly controlled by global galactic properties such as the metallicity, stellar mass and sSFR. However, once photons have escaped from their H II region, the additional feedback mechanisms that we have discussed in this section might play a role to convert $f_{esc,HII}$ into global galactic escape fractions by removing or ionizing the surrounding neutral gas.

### 7.4. Prospects for future observations

We provide in this study an observational framework to pinpoint galaxies with a substantial amount of $f_{esc,HII}$. The presence of leaking H II regions providing a substantial amount of ionizing photons is a necessary condition for a galaxy to be a potential LyC-leaker at larger scales. Hence, the tracers that we propose here could serve as a first basis to select galaxy samples with potential leakage.

We identify line ratios that could be used as tracers of the escape fraction. First, we find that among the classically used optical line ratio, O32 and S32 lines ratio correlate well with our predicted $f_{esc,HII}$. This finding is at odds with observational studies finding little correlations with such diagnostics, and large scatters. In particular, the clear O32-$f_{esc,HII}$ correlation is unexpected in the view of several recent studies that reported
<figure><img src="image_17.png" /><figcaption>**Fig. 18.** Impact of sSFR, X-ray sources and outflows on the escape fraction. **Left:** impact of the measured sSFR from Rémy-Ruyer et al. (2015) on the $f_{esc,HII}$-metallicity relation. No sSFR measurement are available for galaxies 37: HS 0017+1055 and 38: HS 2352+2733. **Middle:** impact of the X-ray-to-stellar luminosity ratio. The green annuli mark galaxies in which a potential ULX ($L_X > 10^{39}$ erg s$^{-1}$) have been reported. 1: Haro 2 (Otí-Floranes et al. 2012), 3: Haro 11 (Prestwich et al. 2015; Gross et al. 2021), 4: He 2–10 (Ott et al. 2005; Reines et al. 2011), 13: 1Zw 18(Thuan et al. 2004; Ott et al. 2005; Kaaret et al. 2011; Kaaret & Feng 2013; Brorby et al. 2014) and 25:SBS 0335-052 (Thuan et al. 2004; Prestwich et al. 2013). **Right:** impact of outflowing gas signatures. 1: Haro 2 Otí-Floranes et al. (2012), 3: Haro 11 (Menacho et al. 2019), 20: NGC 1569 (Sánchez-Cruces et al. 2015), 21: NGC 1705 (Zastrow et al. 2013), 22: NGC 5253 Zastrow et al. (2011), 23: NGC 625 (Cannon et al. 2004; McQuinn et al. 2018), 33: UM 461 (Carvalho & Plana 2018), 39-40: NGC 4214 (McQuinn et al. 2018).</figcaption></figure>

numerous caveats and limits in the use the O32 ratio for studying $f_{esc}$ (e.g., Stasińska et al. 2015; Plat et al. 2019; Bassett et al. 2019; Flury et al. 2022a,b). One possible explanation is that the $f_{esc,HII}$ we derive is less sensitive to viewing angle dependences than the measured escape fractions. Alternatively, this result might also be due to our sample selection, which focuses on highly star-forming galaxies. Additionally, both O32 and S32 ratios have been shown to correlate with the ionization parameter (Nakajima & Ouchi 2014; Kewley & Dopita 2002). Disentangling the effect of metallicity, ionization parameter, and escape fractions remains a difficult task, even within the Bayesian framework that we developed here.

Among the IR tracer, the well-known O3C2 ratio has been used in the analysis of high-z galaxies with ALMA (Harikane et al. 2020). Nevertheless, the interpretation of high O3C2 ratio is still uncertain and debated. Recently, Katz et al. (2022b) have shown that this ratio is sensitive to a wide variety of other galactic parameters including the C/O ratio, IMF, ionization parameter, gas density and CMB attenuation. In our study, assuming that the C/O in the local universe is observationally well constrained, we do find a trend of increasing H II region escape fraction associated with higher O3C2 ratios. However, we find that ratios involving the [N II] lines (122 and 205 µm) provide better diagnostics to trace $f_{esc,HII}$ since [N II] is less contaminated than [C II] by PDR emission. Although fainter, these lines also fall in the detection range of ALMA and have already been detected in a few high-redshift galaxies above $z \sim 4$. (e.g., De Breuck et al. 2019). Even for nondetections, upper limits on [N II] could also be useful to set lower limits on the resulting escape fractions.

For tracers such as [O IV] and [Ne V], associated with ions of high ionization potentials, our interpretation is limited by the little number of constraints that we have on the X-ray component. The JWST will open a new window for the observation of galaxies with signatures of high energy sources and their associated spectral signatures (e.g., [Mg IV]$\lambda$4.49 µm, [Ne VI]$\lambda$7.64 µm, and [Ne V]$\lambda$14.32 µm; Satyapal et al. 2021; Richardson et al. 2022). The combination of such lines with the emission of neutral gas lines (e.g., [Fe II]$\lambda$1.64 µm) boosted by X-ray pho-

tons may help disentangle the feedback from potential compact objects from the one induced by turbulence and shocks (Cresci et al. 2010, 2017). Those IR tracers could be used in combination with other tracers of high energy photons available in the optical range (e.g., [Ne V], [Fe V]) and in the UV range (e.g., C IV $\lambda$5808, He II $\lambda$1640, Wofford et al. 2021; Berg et al. 2019, 2021). A first example of the combination of optical and IR data is presented in LR22. Looking for counter-parts in the radio and X-ray domain with sensitive instruments (e.g., with the forthcoming *Athena* telescope), may also allow progress on the identification of the exact nature of the different sources producing high energy photons.

Finally, the numerous caveats mentioned in this section regarding the distribution of neutral gas might be lifted by incorporating into our models tracers of the neutral gas (e.g., the $Ly\alpha$ line, the hydrogen 21 cm line) or by using priors calibrated on absorption lines studies. Such observations, both in the local and high-redshift universe, will provide a more complex picture of the interplay between low-metallicity ISM and the intense sources of radiation expected in the early universe. In a subsequent publication we will explore and model the optical emission lines of galaxies from the Low-redshift Lyman Continuum Survey with MULTIGRIS (Marques-Chaves et al., in prep.).

## 8. Conclusion

This paper presents the first application of a new Bayesian code, MULTIGRIS (LR22), which opens new possibilities in terms of multiphase modeling of the ISM. We developed a framework to automatically combine photoionization and photodissociation models as different sectors representative of different physical and chemical ISM conditions, and including the presence of density-bounded regions from where ionizing photons can escape. We applied this new code to a sample of low-metallicity dwarf galaxies from the DGS (Madden et al. 2013) with extensive sets of IR emission line measurements and whose ISM properties may resemble those of galaxies in the early Universe. We examined the impact of using multisector models on
the predictions obtained for well-constrained physical quantities such as the metallicity and SFR. We then derived predictions for the H II region escape fraction of ionizing photons in the DGS and studied how this observable varies with various galactic parameters. Our main results are summarized below:

- Accounting for several, possibly density-bounded, sectors does not significantly affect the metallicity and SFR we derive. This confirms that empirical calibrations based on emission lines provide robust results that hold at galactic scales, even when accounting for a slightly more complex topology. We find, however, that not accounting for escaping photons may result in a moderate underestimation of the SFR. The correction factor we derive by estimating SFR based on $Q(H^0)$ tightly follows the analytical curve $1/(1-f_{esc,HII})$.

- We confirm previous results from Cormier et al. (2019) that show that the ISM structure and porosity strongly varies with metallicity. Specifically, we find that $f_{esc,HII}$ tends to increase when metallicity decreases. In the global KDE derived for the whole DGS sample, we find that the high-metallicity and high-$f_{esc,HII}$ region of the parameter space is completely unpopulated, meaning that high-metallicity galaxies have a low probability of reaching high values of $f_{esc,HII}$. Despite a somewhat large scatter around this relation, especially at low-metallicity, this trend is robust to changes in our modeling assumptions on the number of sectors and holds when forcing a fixed number of sectors. Our findings back up the picture of an increasingly porous ISM at low-metallicity, in which energetic photons can easily leak out of H II regions and ionize their surrounding gas.

- We also find that the fraction of [C II] emitted in the ionized phase drops by ~2 orders of magnitude when going from a solar metallicity down to a metallicity of ~1/35Z⊙. We explain this effect by the presence of enhanced radiation fields in low-metallicity environments, which favor the ionization of $C^{2+}$ ions with respect to $C^+$. Hence, we find that [C II] emission arises predominantly from the neutral atomic and molecular gas in low-metallicity environments. A more detailed study of the neutral and molecular gas mass distribution in the DGS, including the CO-dark component, will be carried out in a future work (Ramambason et al., in prep.).

- We examine secondary dependences of $f_{esc,HII}$ with galactic properties and find a weak anti-correlation with stellar mass and a clear correlation with sSFR. These findings fit in the picture of the escape fraction being a feedback-dominated process with complex dependences on the gas mass distribution and gas chemistry.

- Finally, we provide several line ratios that correlate with $f_{esc,HII}$, both in the IR domain (e.g., O3N2, O3C2), including some tracers with high ionization potential (e.g., Ne32, O4C2, and Ne52) and in the optical (e.g., O32, S32), as well as mixed tracers such as $L_{TIR}/L_{bol}$.

The modeling framework presented in this paper is highly sensitive to the wealth of tracers that are used as inputs, especially when they trace gas in different physical and chemical conditions. With that in mind, potentially more complex and realistic models could be constrained by making use of multiwavelength data combining emission from X-ray to IR.

*Acknowledgements.* The authors thank the anonymous referee for constructive feedback and useful comments. V. L., C. R. and L. R. thank the FACE Foundation Thomas Jefferson Fund (grant number TJF21_053). C. M. thanks the grant UNAM/PAPIIT - IN101220. M. C. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)

through an Emmy Noether Research Group (grant numbers KR4801/1-1 and CH2137/1-1), as well as from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement number 714907).

## References

- Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20

- Amorín, R. O., Pérez-Montero, E., & Vílchez, J. M. 2010, ApJ, 715, L128

- Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

- Bakx, T. J. L. C., Tamura, Y., Hashimoto, T., et al. 2020, MNRAS, 493, 4294

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

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

- Bastian, N., Nemslem, E., Kissler-Patig, M., & Maraston, C. 2006, A&A, 445, 471

- Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26

- Berg, D. A., Erb, D. K., Henry, R. B. C., Skillman, E. D., & McQuinn, K. B. W. 2019, ApJ, 874, 93

- Berg, D. A., Chisholm, J., Erb, D. K., et al. 2021, ApJ, 922, 170

- Bergvall, N., Zackrisson, B. G., & Andersson, B. G., et al. 2006, A&A, 448, 513

- Bian, F., Fan, X., McGreer, I., Cai, Z., & Jiang, L. 2017, ApJ, 837, L12

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

- Binder, B., Williams, B. F., Eracleous, M., et al. 2015, AJ, 150, 94

- Binette, L., Wilson, A. S., & Storchi-Bergmann, T. 1996, A&A, 312, 365

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

- Bradley, T. R., Knapen, J. H., Beckman, J. E., & Folkes, S. L. 2006, A&A, 459, L13

- Brorby, M., Kaaret, P., & Prestwich, A. 2014, MNRAS, 441, 2346

- Calzetti, D., Liu, G., & Koda, J. 2012, ApJ, 752, 98

- Cannon, J. M., McClure-Griffiths, N. M., Skillman, E. D., & Côté, S. 2004, ApJ, 607, 274

- Carniani, S., Maiolino, R., Pallottini, A., et al. 2017, A&A, 605, A42

- Carvalho, M. S., & Plana, H. 2018, MNRAS, 481, 122

- Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2019, ApJ, 876, 62

- Chevance, M., Madden, S. C., Lebouteiller, V., et al. 2016, A&A, 590, A36

- Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E., et al. 2020a, Space Sci. Rev., 216, 50

- Chevance, M., Madden, S. C., Fischer, C., et al. 2020b, MNRAS, 494, 5279

- Chevance, M., Kruijssen, J. M. D., Krumholz, M. R., et al. 2022, MNRAS, 509, 272

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

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

- Chisholm, J., Prochaska, J. X., Schaerer, D., Gazagnes, S., & Henry, A. 2020, MNRAS, 498, 2554

- Choi, Y., Dalcanton, J. J., Williams, B. F., et al. 2020, ApJ, 902, 54

- Cormier, D., Lebouteiller, V., Madden, S. C., et al. 2012, A&A, 548, A20

- Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53

- Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23

- Cresci, G., Vanzi, L., Sauvage, M., Santangelo, G., & van der Werf, P. 2010, A&A, 520, A82

- Cresci, G., Vanzi, L., Telles, E., et al. 2017, A&A, 604, A101

- Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731

- Croxall, K. V., Smith, J. D., Pellegrini, E., et al. 2017, ApJ, 845, 96

- de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51

- De Breuck, C., Weiß, A., Béthermin, M., et al. 2019, A&A, 631, A167

- De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62

- Della Bruna, L., Adamo, A., Bik, A., et al. 2020, A&A, 635, A134

- Della Bruna, L., Adamo, A., Lee, J. C., et al. 2021, A&A, 650, A103

- Duarte Puertas, S., Vilchez, J. M., Iglesias-Páramo, J., et al. 2022, A&A, 666, A186

- Eggen, N. R., Scarlata, C., Skillman, E., & Jaskot, A. 2021, ApJ, 912, 12

- Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058

- Espinosa-Ponce, C., Sánchez, S. F., Morisset, C., et al. 2020, MNRAS, 494, 1622

- Falkendal, T., Lehnert, M. D., Vernet, J., De Breuck, C., & Wang, W. 2021, A&A, 645, A120

- Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385

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

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

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

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

- Fujita, A., Martin, C. L., Mac Low, M.-M., & Abel, T. 2003, ApJ, 599, 50

- Galliano, F. 2018, MNRAS, 476, 1445

- Galliano, F., Nersessian, A., Bianchi, S., et al. 2021, A&A, 649, A18

- Gazagnes, S., Chisholm, J., Schaerer, D., et al. 2018, A&A, 616, A29
Gazagnes, S., Chisholm, J., Schaerer, D., Verhamme, A., & Izotov, Y. 2020, A&A, 639, A85

Grimm, H. J., Gilfanov, M., & Sunyaev, R. 2003, MNRAS, 339, 793

Gross, A. C., Prestwich, A., & Kaaret, P. 2021, MNRAS, 505, 610

Hanish, D. J., Oey, M. S., Rigby, J. R., de Mello, D. F., & Lee, J. C. 2010, ApJ, 725, 2029

Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93

Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71

Henry, A., Scarlata, C., Martin, C. L., & Erb, D. 2015, ApJ, 809, 19

Henry, A., Berg, D. A., Scarlata, C., Verhamme, A., & Erb, D. 2018, ApJ, 855, 96

Herenz, E. C., Hayes, M., Papaderos, P., et al. 2017, A&A, 606, L11

Hogarth, L., Amorín, R., Vílchez, J. M., et al. 2020, MNRAS, 494, 3541

Hosokawa, T., & Inutsuka, S.-I. 2005, ApJ, 623, 917

Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736

Inoue, A. K., Hirashita, H., & Kamaya, H. 2001, ApJ, 555, 613

Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559

Izotov, Y. I., & Thuan, T. X. 1999, ApJ, 511, 639

Izotov, Y. I., Thuan, T. X., & Lipovetsky, V. A. 1994, ApJ, 435, 647

Izotov, Y. I., Chaffee, F. H., & Green, R. F. 2001, ApJ, 562, 727

Izotov, Y. I., Noeske, K. G., Guseva, N. G., et al. 2004a, A&A, 415, L27

Izotov, Y. I., Papaderos, P., Guseva, N. G., Fricke, K. J., & Thuan, T. X. 2004b, A&A, 421, 539

Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955

Izotov, Y. I., Thuan, T. X., & Privon, G. 2012, MNRAS, 427, 1229

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

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

Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2017, MNRAS, 471, 548

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

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

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

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

James, B. L., Aloisi, A., Heckman, T., Sohn, S. T., & Wolfe, M. A. 2014, ApJ, 795, 109

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

Jiang, T., Malhotra, S., Rhoads, J. E., & Yang, H. 2019, ApJ, 872, 145

Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46

Kaaret, P., & Feng, H. 2013, AAS/High Energy Astrophys. Div., 13, 402.06

Kaaret, P., Schmitt, J., & Gorski, M. 2011, ApJ, 741, 10

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

Katz, H., Ďurovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164

Katz, H., Garel, T., Rosdahl, J., et al. 2022a, MNRAS, 515, 4265

Katz, H., Rosdahl, J., Kimm, T., et al. 2022b, MNRAS, 510, 5603

Kehrig, C., Vílchez, J. M., Guerrero, M. A., et al. 2018, MNRAS, 480, 1081

Kennicutt, R. C., Jr., Bresolin, F., Bomans, D. J., Bothun, G. D., & Thompson, I. B. 1995, AJ, 109, 594

Kennicutt, R. C., Jr, Hao, C.-N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35

Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2022, MNRAS, 516, 3006

Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68

Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2019, ApJ, 883, 102

Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021, ApJ, 911, 128

Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826

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

Kreckel, K., Blanc, G. A., Schinnerer, E., et al. 2016, ApJ, 827, 103

Kunth, D., Mas-Hesse, J. M., Terlevich, E., et al. 1998, A&A, 334, 11

Lacerda, E. A. D., Cid Fernandes, R., Couto, G. S., et al. 2018, MNRAS, 474, 3727

Lambert-Huyge, A., Madden, S. C., Lebouteiller, V., et al. 2022, A&A, 666, A112

Lardo, C., Davies, B., Kudritzki, R. P., et al. 2015, ApJ, 812, 160

Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506

Lebouteiller, V., & Ramambason, L. 2022, A&A, 667, A34

Lebouteiller, V., Péquignot, D., Cormier, D., et al. 2017, A&A, 602, A45

Lebouteiller, V., Cormier, D., Madden, S. C., et al. 2019, A&A, 632, A106

Lee, M. Y., Madden, S. C., Lebouteiller, V., et al. 2016, A&A, 596, A85

Lee, M. Y., Madden, S. C., Le Petit, F., et al. 2019, A&A, 628, A113

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Leitet, E., Bergvall, N., Hayes, M., Linné, S., & Zackrisson, E. 2013, A&A, 553, A106

Leitherer, C., Hernandez, S., Lee, J. C., & Oey, M. S. 2016, ApJ, 823, 64

Ma, X., Hopkins, P. F., Kasen, D., et al. 2016, MNRAS, 459, 3614

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

Madden, S. C., Rémy-Ruyer, A., Galametz, M., et al. 2013, PASP, 125, 600

Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141

Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3

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

Marques-Chaves, R., Schaerer, D., Amorín, R. O., et al. 2022, A&A, 663, L1

Martin, C. L. 1998, ApJ, 506, 222

Mathis, J. S. 1986, ApJ, 301, 423

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

McQuinn, K. B. W., Skillman, E. D., Heilman, T. N., Mitchell, N. P., & Kelley, T. 2018, MNRAS, 477, 3164

Menacho, V., Östlin, G., Bik, A., et al. 2019, MNRAS, 487, 3183

Menacho, V., Östlin, G., Bik, A., et al. 2021, MNRAS, 506, 1777

Meyer, R. A., Walter, F., Cicone, C., et al. 2022, ApJ, 927, 152

Mingozzi, M., Belfiore, F., Cresci, G., et al. 2020, A&A, 636, A42

Mirabel, I. F., Dijkstra, M., Laurent, P., Loeb, A., & Pritchard, J. R. 2011, A&A, 528, A149

Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741

Naidu, R. P., Oesch, P. A., Reddy, N., et al. 2017, ApJ, 847, 12

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

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

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

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

Nicholls, D. C., Sutherland, R. S., Dopita, M. A., Kewley, L. J., & Groves, B. A. 2017, MNRAS, 466, 4403

Oberst, T. E., Parshley, S. C., Nikola, T., et al. 2011, ApJ, 739, 100

Oey, M. S., Meurer, G. R., Yelda, S., et al. 2007, ApJ, 661, 801

Olivier, G. M., Berg, D. A., Chisholm, J., et al. 2021, ApJ, submitted [arXiv:2109.06725]

Otí-Floranes, H., Mas-Hesse, J. M., Jiménez-Bailón, E., et al. 2012, A&A, 546, A65

Ott, J., Walter, F., & Brinks, E. 2005, MNRAS, 358, 1453

Pahl, A. J., Shapley, A., Steidel, C. C., Chen, Y., & Reddy, N. A. 2021, MNRAS, 505, 2447

Pahl, A. J., Shapley, A., Steidel, C. C., Reddy, N. A., & Chen, Y. 2022, MNRAS, 516, 2062

Papaderos, P., Fricke, K. J., Thuan, T. X., & Loose, H. H. 1994, A&A, 291, L13

Pellegrini, E. W., Oey, M. S., Winkler, P. F., et al. 2012, ApJ, 755, 40

Péquignot, D. 2008, A&A, 478, 371

Pilyugin, L. S., & Thuan, T. X. 2005, ApJ, 631, 231

Pilyugin, L. S., Grebel, E. K., Zinchenko, I. A., & Kniazev, A. Y. 2014, AJ, 148, 134

Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978

Polles, F. L., Madden, S. C., Lebouteiller, V., et al. 2019, A&A, 622, A119

Prestwich, A. H., Tsantaki, M., Zezas, A., et al. 2013, ApJ, 769, 92

Prestwich, A. H., Jackson, F., Kaaret, P., et al. 2015, ApJ, 812, 166

Ramambason, L., Schaerer, D., Stasińska, G., et al. 2020, A&A, 644, A21

Ramos Padilla, A. F., Wang, L., Ploeckinger, S., van der Tak, F. F. S., & Trager, S. C. 2021, A&A, 645, A133

Ramos Padilla, A. F., Wang, L., van der Tak, F. F. S., & Trager, S. 2022, A&A, submitted [arXiv:2205.11955]

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

Reines, A. E., Sivakoff, G. R., Johnson, K. E., & Brogan, C. L. 2011, Nature, 470, 66

Relaño, M., De Looze, I., Kennicutt, R. C., et al. 2018, A&A, 613, A43

Rémy-Ruyer, A. 2013, Ph.D. Thesis

Richardson, C. T., Simpson, C., Polimera, M. S., et al. 2022, ApJ, 927, 165

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

Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71

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

Rosdahl, J., Katz, H., Blaizot, J., et al. 2018, MNRAS, 479, 994

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121

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

Sánchez-Cruces, M., Rosado, M., Rodríguez-González, A., & Reyes-Iturbide, J. 2015, ApJ, 799, 231

Satyapal, S., Kamal, L., Cann, J. M., Secrest, N. J., & Abel, N. P. 2021, ApJ, 906, 35

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

Schaerer, D., Contini, T., & Pindao, M. 1999, A&AS, 136, 35

Schaerer, D., Fragos, T., & Izotov, Y. I. 2019, A&A, 622, L10

Sembach, K. R., Howk, J. C., Ryans, R. S. I., & Keenan, F. P. 2000, ApJ, 528, 310

Senchyna, P., Stark, D. P., Mirocha, J., et al. 2020, MNRAS, 494, 941

Senchyna, P., Stark, D. P., Charlot, S., et al. 2021, MNRAS, 503, 6112

Senchyna, P., Stark, D. P., Charlot, S., et al. 2022, ApJ, 930, 105
Shapley, A. E., Steidel, C. C., Strom, A. L., et al. 2016, ApJ, 826, L24

Simmonds, C., Schaerer, D., & Verhamme, A. 2021, A&A, 656, A127

Stanway, E. R., & Eldridge, J. J. 2019, A&A, 621, A105

Stark, D. P., Walth, G., Charlot, S., et al. 2015, MNRAS, 454, 1393

Stasińska, G., Izotov, Y., Morisset, C., & Guseva, N. 2015, A&A, 576, A83

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

Sutherland, R., Dopita, M., Binette, L., & Groves, B. 2018, Astrophysics Source Code Library [record ascl:1807.005]

Thuan, T. X., Bauer, F. E., Papaderos, P., & Izotov, Y. I. 2004, ApJ, 606, 213

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

Ugryumov, A. V., Engels, D., Pustilnik, S. A., et al. 2003, A&A, 397, 463

Umeda, H., Ouchi, M., Nakajima, K., et al. 2022, ApJ, 930, 37

Vale Asari, N., Couto, G. S., Cid Fernandes, R., et al. 2019, MNRAS, 489, 4721

Vale Asari, N., Wild, V., de Amorim, A. L., et al. 2020, MNRAS, 498, 4205

Vallini, L., Ferrara, A., Pallottini, A., Carniani, S., & Gallerani, S. 2021, MNRAS, 505, 5543

Vanzella, E., de Barros, S., Castellano, M., et al. 2015, A&A, 576, A116

Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41

Vanzella, E., Nonino, M., Cupani, G., et al. 2018, MNRAS, 476, L15

Vanzella, E., Caminha, G. B., Calura, F., et al. 2020, MNRAS, 491, 1093

Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7

Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13

Walter, F., Riechers, D., Novak, M., et al. 2018, ApJ, 869, L22

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

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

Weilbacher, P. M., Monreal-Ibero, A., Verhamme, A., et al. 2018, A&A, 611, A95

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Westmoquette, M. S., Smith, L. J., & Gallagher, J. S. 2008, MNRAS, 383, 864

Wofford, A., Vidal-García, A., Feltre, A., et al. 2021, MNRAS, 500, 2908

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191

Wood, K., Hill, A. S., Joung, M. R., et al. 2010, ApJ, 721, 1397

Xiao, L., Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 477, 904

Xu, X., Henry, A., Heckman, T., et al. 2022, ApJ, 933, 202

Yoo, T., Kimm, T., & Rosdahl, J. 2020, MNRAS, 499, 5175

Youngblood, A. J., & Hunter, D. A. 1999, ApJ, 519, 55

Zastrow, J., Oey, M. S., Veilleux, S., McDonald, M., & Martin, C. L. 2011, ApJ, 741, L17

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

Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217

Zurita, A., Beckman, J. E., Rozas, M., & Ryder, S. 2002, A&A, 386, 801
# Appendix A: Additional notes and figures

In this section we gather additional figures that complement the main plots presented in this article. Those figures are mentioned throughout the article to illustrate some of the discussions but do not display essential results. We recall in the captions the corresponding sections to which they are related.

<figure><img src="image_18.png" /><figcaption>**Fig. A.1.** Abundance profile used to compute the grid of Cloudy models. This plot complements the description of the grid in Sect. 3 and illustrate the prescription described in Sect. 3.2.4.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. A.2.** Predicted vs. measured H I masses. As discussed in Sect. 7.2.2, the predicted H I masses for the galaxies in our sample are consistent with the observed H I masses from Rémy-Ruyer et al. (2014). This agreement is at odds with the underprediction of H I column densities (see Table A.3) and is discussed in Sect. 7.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. A.3.** Metallicity-$f_{esc,HI}$ relation for a fixed number of sectors for all galaxies. These plots complement the Fig. 9 presented in Sect. 6.1 where the optimal number of sectors for a galaxy is chosen based on the score criteria presented in Sect. 4. We show here that the $f_{esc,HI}$-metallicity trend presented in Sect. 6.1 is robust and holds when we force an arbitrary number of sectors. We note that for single sector models, most solutions are completely radiation-bounded leading to escape fractions of 0%.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. A.4.** Detected lines (red dots) and upper limits (blue triangle) used in this paper. This figure illustrates the description of the available constraints for each galaxy in our sample that is presented in Sect. 2. The total number of detections and upper limits used for each object are reported on the right-hand side y-axis.</figcaption></figure>
<table><caption>**Table A.1.** Best configuration selected for each galaxy and corresponding values of some indicators used to check model accuracy.</caption><thead><tr><th>#</th><th>Galaxy</th><th>N<sub>sectors,best</sub></th><th>p-3σ [%]<sup>1</sup></th><th>L<sub>M,best</sub><sup>2</sup></th><th>ΔL<sub>M,1</sub><sup>3</sup></th><th>ΔL<sub>M,2</sub><sup>3</sup></th><th>ΔL<sub>M,3</sub><sup>3</sup></th></tr></thead><tbody><tr><td>1</td><td>Haro2</td><td>3.0</td><td>87.01</td><td>-45.49</td><td>-12.98</td><td>-4.41</td><td>0.0</td></tr><tr><td>2</td><td>Haro3</td><td>2.0</td><td>74.67</td><td>-60.51</td><td>-14.58</td><td>0.0</td><td>-1.28</td></tr><tr><td>3</td><td>Haro11</td><td>3.0</td><td>86.35</td><td>-51.95</td><td>-16.09</td><td>-1.32</td><td>0.0</td></tr><tr><td>4</td><td>He2-10</td><td>3.0</td><td>75.03</td><td>-93.15</td><td>-17.21</td><td>-2.59</td><td>0.0</td></tr><tr><td>5</td><td>HS0052+2536</td><td>2.0</td><td>86.21</td><td>-31.15</td><td>-0.17</td><td>0.0</td><td>-0.48</td></tr><tr><td>6</td><td>HS0822+3542</td><td>1.0</td><td>84.76</td><td>-13.19</td><td>0.0</td><td>-1.4</td><td>-3.62</td></tr><tr><td>7</td><td>HS1222+3741</td><td>1.0</td><td>88.61</td><td>-10.64</td><td>0.0</td><td>-0.96</td><td>-1.14</td></tr><tr><td>8</td><td>HS1304+3529</td><td>1.0</td><td>78.53</td><td>-18.49</td><td>0.0</td><td>-0.36</td><td>-0.75</td></tr><tr><td>9</td><td>HS1319+3224</td><td>3.0</td><td>94.67</td><td>-3.02</td><td>-1.25</td><td>-0.88</td><td>0.0</td></tr><tr><td>10</td><td>HS1330+3651</td><td>1.0</td><td>89.63</td><td>-24.35</td><td>0.0</td><td>-2.07</td><td>-0.64</td></tr><tr><td>11</td><td>HS1442+4250</td><td>3.0</td><td>74.18</td><td>-21.04</td><td>-0.23</td><td>-0.38</td><td>0.0</td></tr><tr><td>12</td><td>IIZw40</td><td>2.0</td><td>96.07</td><td>-39.4</td><td>-13.23</td><td>0.0</td><td>-1.6</td></tr><tr><td>13</td><td>IZw18</td><td>3.0</td><td>87.91</td><td>-31.75</td><td>-10.68</td><td>-2.24</td><td>0.0</td></tr><tr><td>14</td><td>Mrk153</td><td>2.0</td><td>97.54</td><td>-17.36</td><td>-1.42</td><td>0.0</td><td>-0.72</td></tr><tr><td>15</td><td>Mrk209</td><td>2.0</td><td>88.04</td><td>-24.83</td><td>-2.81</td><td>0.0</td><td>-0.64</td></tr><tr><td>16</td><td>Mrk930</td><td>1.0</td><td>98.09</td><td>-26.55</td><td>0.0</td><td>-0.67</td><td>-2.18</td></tr><tr><td>17</td><td>Mrk1089</td><td>3.0</td><td>80.25</td><td>-43.25</td><td>-10.14</td><td>-1.52</td><td>0.0</td></tr><tr><td>18</td><td>Mrk1450</td><td>3.0</td><td>87.57</td><td>-45.83</td><td>-9.31</td><td>-0.07</td><td>0.0</td></tr><tr><td>19</td><td>NGC1140</td><td>3.0</td><td>76.21</td><td>-73.47</td><td>-0.72</td><td>-3.14</td><td>0.0</td></tr><tr><td>20</td><td>NGC1569</td><td>2.0</td><td>74.77</td><td>-62.13</td><td>-5.98</td><td>0.0</td><td>-0.5</td></tr><tr><td>21</td><td>NGC1705</td><td>2.0</td><td>74.48</td><td>-51.53</td><td>-1.65</td><td>0.0</td><td>-1.45</td></tr><tr><td>22</td><td>NGC5253</td><td>2.0</td><td>76.91</td><td>-51.4</td><td>-4.36</td><td>0.0</td><td>-0.84</td></tr><tr><td>23</td><td>NGC625</td><td>3.0</td><td>74.99</td><td>-57.43</td><td>-4.6</td><td>-1.31</td><td>0.0</td></tr><tr><td>24</td><td>Pox186</td><td>3.0</td><td>92.09</td><td>-12.49</td><td>-3.46</td><td>-1.45</td><td>0.0</td></tr><tr><td>25</td><td>SBS0335-052</td><td>3.0</td><td>80.58</td><td>-34.49</td><td>-2.05</td><td>-0.34</td><td>0.0</td></tr><tr><td>26</td><td>SBS1159+545</td><td>2.0</td><td>78.39</td><td>-22.8</td><td>-0.14</td><td>0.0</td><td>-1.57</td></tr><tr><td>27</td><td>SBS1211+540</td><td>1.0</td><td>89.23</td><td>-16.74</td><td>0.0</td><td>-1.5</td><td>-1.94</td></tr><tr><td>28</td><td>SBS1249+493</td><td>2.0</td><td>89.72</td><td>-8.14</td><td>-2.76</td><td>0.0</td><td>-0.66</td></tr><tr><td>29</td><td>SBS1415+437</td><td>1.0</td><td>77.33</td><td>-36.14</td><td>0.0</td><td>-1.13</td><td>-0.01</td></tr><tr><td>30</td><td>SBS1533+574</td><td>2.0</td><td>77.95</td><td>-32.51</td><td>-0.83</td><td>0.0</td><td>-0.05</td></tr><tr><td>31</td><td>Toll1214-277</td><td>1.0</td><td>61.86</td><td>-21.46</td><td>0.0</td><td>-4.99</td><td>-4.56</td></tr><tr><td>32</td><td>UM448</td><td>3.0</td><td>75.5</td><td>-67.58</td><td>-2.34</td><td>-4.63</td><td>0.0</td></tr><tr><td>33</td><td>UM461</td><td>2.0</td><td>84.32</td><td>-40.5</td><td>-10.49</td><td>0.0</td><td>-0.47</td></tr><tr><td>34</td><td>VIIZw403</td><td>2.0</td><td>83.39</td><td>-34.71</td><td>-1.59</td><td>0.0</td><td>-0.05</td></tr><tr><td>35</td><td>UGC4483</td><td>2.0</td><td>75.98</td><td>-11.14</td><td>-0.81</td><td>0.0</td><td>-0.14</td></tr><tr><td>36</td><td>UM133</td><td>1.0</td><td>99.15</td><td>-2.47</td><td>0.0</td><td>-1.09</td><td>-0.99</td></tr><tr><td>37</td><td>HS0017+1055</td><td>1.0</td><td>96.12</td><td>-2.77</td><td>0.0</td><td>-0.05</td><td>-0.09</td></tr><tr><td>38</td><td>HS2352+2733</td><td>1.0</td><td>96.6</td><td>-1.08</td><td>0.0</td><td>-0.17</td><td>-0.27</td></tr><tr><td>39</td><td>NGC4214-c</td><td>3.0</td><td>90.66</td><td>-17.95</td><td>-9.85</td><td>-0.39</td><td>0.0</td></tr><tr><td>40</td><td>NGC4214-s</td><td>2.0</td><td>95.49</td><td>-11.49</td><td>-19.3</td><td>0.0</td><td>-2.45</td></tr></tbody></table>

(1): Probability that the best model predictions are within $3\sigma$ of the observations, for all the lines used as constraints. (2): Absolute value of the logarithm of the marginal likelihood for the best configuration corresponding to $N_{sectors,best}$. (3): $\Delta L_{M,i} = \ln(L_{M,i}/L_{M,best})$
<figure><img src="image_22.png" /><figcaption>**Fig. A.5.** Evolution of the volume-averaged ionization parameter log 〈U〉 as a function log U_in and log U_out in the model grid. We show only a subsample corresponding to models with no X-ray source and stopping at the ionization front (cut=1). This figure complements the discussion on the possible caveats induced by the use of U_in as a primary parameter in our grid (see Sect. 3.2.3). We note that U_in only sets an upper limit on the actual volume averaged log 〈U〉. Variations in terms of metallicity and initial density n_0 yield different geometries, which cover a large range of volume averaged log 〈U〉 for the same U_in.</figcaption></figure>
<figure><img src="image_23.png" /></figure>

**Fig. A.6.** Impact of measured and inferred galactic parameters on the escape fraction vs metallicity relation. **Panels 1 to 3:** Impact of measured galactic parameters on the escape fraction vs metallicity relation: M*, SFR, sSFR. **Panels 3 to 6:** Impact of inferred galactic parameters on the escape fraction vs metallicity relation: stellar age, Lx, sSFR(Q0). This figure complements the results presented in Sect. 6.1. It provides a reinterpretation of Fig. 9 but accounting for the additional dependences on galactic parameters by splitting the sample by bins of values for secondary parameters. Those secondary parameters are further examined in Sect. 6.2.