© L. Ramambason et al. 2022

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₂ 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; 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 ionizing 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_esc(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⁸ M_⊙ and extreme escape fractions (the “oligarchs”).

Regardless of which galaxy population drives the reionization process, the inclusion of binary stars, which provide ionizing photons at later stellar evolution stages than single stars, might also play a crucial role in providing energetic photons over large timescales. Simulations of Ma et al. (2016) and Rosdahl et al. (2018) find that this leads to significantly higher time-averaged escape fractions of ionizing photons. In low-metallicity environments hosting very massive stars, stellar-mass black holes might also contribute significantly 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_esc(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 observation of potentially LyC-leaking galaxies directly within the EoR. This observational barrier makes it difficult to perform quantitative 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 sensitive 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 absorption 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 selection 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 indicators of LyC leakage. In particular, line ratios involving ions with different ionization potentials, produced at different depths, have been proposed as indicators to discriminate between radiation-bounded and density-bounded H II regions. Radiation-bounded regions correspond to ionized spheres delimited by their Strömgren radii set by the equilibrium between production of photons by stars and ionization of the surrounding gas. Density-bounded regions are instead delimited by the lack of matter, which sets their outer radius before the Strömgren radius. Hence, they allow part of the LyC-photons produced 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 successfully used to select LyC-leaking candidates but no strong correlation was found with the measured values of escape fraction (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 ionization 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 resolution. A large fraction of those escaping photons is, however, likely reabsorbed within galactic scales and contributes to create 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. Aditionally, 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 I Zw 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, I Zw 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 χ² 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 χ² 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, HII).

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 spectrocopic 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-σ of 0.112 × 10¹⁶W m⁻². We also include the total IR luminosities (L_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_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 χ² 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α 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 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³⁹ erg s⁻¹: Haro 2 (Otí-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¹ metallicities ranging from 12+log(O/H) = 7.14 (∼1/35 Z_⊙) up to 8.43 (∼1/2 Z_⊙), 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⁷ to ∼10¹¹ M_⊙ (Rémy-Ruyer et al. 2014) and specific gass masses (M_HI/M_*) 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²¹ cm⁻² (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; 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 χ² 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α emission drops. We stress, however, that for the lowest metallicity models, the contribution from binary stars delays this drop in Hα 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 and constant star-formation rate models corresponding to similar ages (Chisholm et al. 2019).

Fig. 1. Stellar (left) and X-ray (right) incident SEDs of Cloudy models for a solar metallicity and bolometric luminosity of 109 L⊙. QE shows the number of photons produced above a given energy Eν. The vertical dashed lines represent the ionization potentials of H+ (13.6 eV), O2+ (35.1 eV), O3+ (54.9 eV) and Ne4+ (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%.

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_⊙, which is the default IMF in Eldridge et al. (2017). Although very massive stars with masses above 100 M_⊙ 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_⊙.

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³ K and an inner temperature varying from 10⁵ to 10⁷ 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]² 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 ∼ 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 primarly 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:

$$\begin{array}{c}\hfill U_{\mathrm{in}}=\frac{Q({\mathrm{H}}^0)}{4\pi nc{R}_{\mathrm{in}}^2},\end{array}$$(1)

where n is the hydrogen density at the inner radius, U_in is the input ionization parameter, c is the speed of light, and Q(H⁰) is the total number ionizing photons produced by the central source defined as follows:

$$\begin{array}{c}\hfill Q({\mathrm{H}}^0)={\displaystyle {\int }_{E_{\nu }=13.6eV}^{\infty }}\frac{L_{\ast }(\nu )+L_X(\nu )}{h\nu }\mathrm{d}\nu ,\end{array}$$(2)

where L_*(ν) is the stellar luminosity and L_X(ν) the luminosity of the X-ray source emitted with a given frequency ν. We note that this definition of U_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_in) at the illuminated front can be very different from the ionization parameter at the ionization front U_out. As an illustration, the relation between log U_in, log U_out and the volume-averaged log ⟨U⟩ 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_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 I]/[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 much stronger effect on lines emitted near the ionization front and in the PDR and molecular zone.

The bolometric luminosity of 10⁹ L_⊙ 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³⁹ erg s⁻¹ are numerous amoung 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⁷ L_⊙ (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²¹ 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 the 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.

Fig. 2. Schematic view of the 17 cuts used to create sub-models.

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 (τ_λ) 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 follow:

$$\begin{array}{c}\hfill f_{\mathrm{esc},\mathrm{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_{\mathrm{in}})},\end{array}$$(3)

where the numerator is the luminosity produced by the central sources reaching the radius R with energy between E₁ and E₂, and the denominator is the luminosity impinging the cloud at the inner radius R_in, within the same energy range. The number of ionizing photons at a given radius is calculated as:

$$\begin{array}{c}\hfill L(E_1\le h\nu \le E_2)(R)=4\pi R^2{\displaystyle {\int }_{{\nu }_1}^{{\nu }_2}}{F}_{\nu }(r)\mathrm{d}\nu ,\end{array}$$(4)

where F_ν in the flux of photons with frequency ν at a given radius. In practice we consider f_esc, HII = f_esc(1Ryd ≤ hν  ≤  ∞)(R_cut) that includes all ionizing photons even in the X-ray regime, with R_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_esc, HII as well. This definition of f_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_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_esc(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.

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α 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α 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).

Fig. 3. Schematic view of the ISM of a starburst galaxy and associated representative topology.

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 unresolved 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:

$$\begin{array}{c}\hfill L={\displaystyle \sum _{\mathrm{i}=1}^{N_{\mathrm{sectors}}}}{w}_{\mathrm{i}}{L}_{\mathrm{i}},\end{array}$$(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_esc, HII is also an extensive quantity and the averaged f_esc, HII can be derive as follows:

$$\begin{array}{c}\hfill f_{\mathrm{esc},\mathrm{HII}}={\displaystyle \sum _{\mathrm{i}=1}^{N_{\mathrm{sectors}}}}{w}_{\mathrm{i}}{f}_{\mathrm{esc},\mathrm{HII}}^{\mathrm{i}}.\end{array}$$(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|θ, ℳ), where O is the data and θ the set of parameters of a model ℳ, 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 σ corresponding to the uncertainty of the observation. Hence, the likelihood can be expressed as:

$$\begin{array}{c}\hfill \mathcal{L}=p(\mathbf{O}|\mathit{\theta },\mathcal{M})={\displaystyle \prod _{i=0}^{N_{\mathrm{obs}}}}\mathcal{N}(\mu =O_{\mathrm{i}},{\sigma }^2=U_{\mathrm{i}}^2),\end{array}$$(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σ 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:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathcal{M}}={\displaystyle {\int }_{\theta }}p(\mathbf{O}|\mathit{\theta },\mathcal{M})\mathrm{d}\mathit{\theta }.\end{array}$$(8)

Because of the integration on all the free parameters θ, 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 I Zw 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σ 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:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{i}\mathrm{sectors}}=\frac{{\mathcal{L}}_{\mathcal{M},\mathrm{i}}}{{\sum }_{\mathrm{j}=1}^{3\mathrm{sectors}}{\mathcal{L}}_{\mathcal{M},\mathrm{j}}},\end{array}$$(9)

where ℒ_ℳ,i is the marginal likelihood associated with the configuration having i sectors (with 1 ≤ i ≤ 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.

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.

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 χ² only provide a group of the most likely values of a parameter, MULTIGRIS outputs a sample of draws from the posterior PDF p(O|θ, ℳ) 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 asymetric distribution. Ideally, one would like to represent the full posterior PDFs that are used in this study to derive trends and correlations. 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σ 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σ) 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, HII: 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 metallicity with a large width parameter (σ = 0.1 + σ_mes in logarithmic scale, where σ_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 σ = 0.1 dex envelop 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.

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σ 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 α line (7–6) is detected or not.

The uncertainties derived (1σ 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 prescription for the abundance patterns vs. metallicity). When a hydrogen recombination line (e.g., the Humphreys α line (7–6); Huα) is available, the estimated metallicity tends to be more tightly constrained as the code can directly infer elemental abundances relative to hydrogen. This is illustrated in Fig. 5 where galaxies that have Huα detection tend to have smaller 1σ contours than the ones without. Interestingly, our code is also able to estimate metallicities without hydrogen recombination line measurements, in some cases with a precision comparable to galaxies with Huα 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 measurement 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 predict 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 constraints, 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_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 diagnostics. Although the SFR is not one of the primary (used for inference) or secondary (output from Cloudy) parameters, we use classical SFR proxies (Hα, L_TIR and Q(H⁰)), which we convert into SFR. We emphasize that Hα 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 ionizing photons produced by the central cluster, Q(H⁰). The predicted Hα luminosity depends of the incident flux Q(H⁰) and on the geometry of the gas in our representative model. The main advantage of our method is that it takes into account information coming from all the emission lines available. However, like for other diagnostics, our estimates ultimately depend on the assumption made regarding the stellar population.

We compare our predictions to a classical diagnostic based on composite diagnostics that combine tracers of dust-obscured (e.g., L_TIR) and dust-unobscured (e.g., Hα) star-formation (see review from Calzetti et al. 2012). First, we compare the predicted Hα luminosities to the values from the literature gathered in Rémy-Ruyer et al. (2015) in Fig. 6. The latter correspond to Hα 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α emission i.e. what escapes from the galaxy without considering dust extinction.

Fig. 6. Agreement between Hα and SFR predictions with empirical measurements. Top: predicted intrinsic Hα emission vs. measured Hα corrected for Galactic extinction from (Rémy-Ruyer et al. 2015, and references therein). Bottom: predicted SFR vs. SFR(Hα+LTIR) from Rémy-Ruyer et al. (2015). The 1σ 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α 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.

The predicted Hα 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α 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α 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α emission and remains compatible with observations using somewhat simple prescriptions to correct the observed Hα 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α and L_TIR). This diagnostic corresponds to an SFR conversion assuming a Kroupa IMF³ 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α 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: Mkr 153), which is linked to our overestimation of Hα emission⁴.

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⁰), 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α + L_TIR. Our SFR estimates based on Q(H⁰) are in most cases shifted to higher values. If we consider SFR(Q(H⁰)) as the intrinsic SFR, the correction to apply to SFR(Hα+TIR) follows the analytical curve 1/(1-f_esc, HII) with small deviations, which can be attributed to several effects. First, although Hα emission is directly proportional to the number of LyC photons absorbed by the gas, (1 − f_esc, HII)Q(H⁰), this conversion depends on the hydrogen recombination coefficient that varies with physical conditions (density and electronic temperature). Second, the 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α+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⁰)) is the intrinsic SFR accounting for all the photons produced by the central source, the correction to apply to SFR(Hα+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).

Fig. 7. Predicted SFR(Q(H0))/SFR(Hα+TIR) vs fesc, HII. The symbols represent the robust means of the posterior PDF and the gray contours show the 1σ uncertainties of the SUEs. The dashed line corresponds to the analytical curve 1/(1–fesc).

The difference between the predicted SFR(Hα+TIR) and SFR(Q(H⁰)) 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⁰) but escape before ionizing the gas that produces Hα. 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α 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⁰)). 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α emission, which might be powered by leaky H II regions. In Fig. 8, we show the predictions for the maximum additional Hα emission powered by photon leakage, assuming a Q(H⁰)-to-Hα conversion coefficient of 7.31 × 10¹¹ for a case B recombination with T_e = 10 000 K and n = 100 cm⁻³ (Kennicutt et al. 1995) to derive Hα(f_esc, HII× Q(H⁰)). 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α(f_esc, HII×Q(H⁰)) while Hα_pred is calculated by integrating Hα 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α 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α emission powered by DIG; this maximum emission exceeds 50% of the total Hα value for the larger f_esc, HII we infer (i.e., f_esc, HII > 50%). Although a significant fraction of the Hα 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α component can hardly be reconciled with the fact that our predicted Hα 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.

Fig. 8. Ratio of the maximum Hα emission associated with photon leakage from H II regions, assuming that all photons are reabsorbed and of the Hα value predicted by our models vs. fesc, HII. We assumed Q(H0)-to-Hα conversion coefficient of 7.31 × 1011 for a case B recombination with Te = 10 000 K and n = 100 cm−3 from Kennicutt et al. (1995) to estimate Hα(fesc, HII× Q(H0)). The Hαpred values are computed consistently in the Cloudy models. The dashed line correspond to the analytical curve fesc/(1 − fesc).

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⁻⁶. 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.

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 Nsectors ≥ 2 only and ρ 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 fesc, 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 fesc, HII for galaxies with at least two sectors.

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 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 metallicity 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 dispersion 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 parameters (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 samples 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 functions 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. However, our analysis is limited by the observational bias of our sample 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 fractions (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 candidates for potential LyC-leakage. In particular, this f_esc − M_* trend has been observed in samples with masses down to 10⁸ M_⊙ (Leitet et al. 2013; Izotov et al. 2018b). However, a recent follow up study carried out in Izotov et al. (2021) targeting galaxies with masses below 10⁸ M_⊙ 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 associated with low stellar masses, although with a large scatter. All the galaxies with masses above 10⁹ M_⊙ have mean escape fractions below 40%. On the other hand, galaxies with masses below 10⁹ M_⊙ do not always have high escape fractions. We note that this effect might be dominated by the mass-metallicity distribution of our sample in which lower masses are associated with the lowest metallicity. The correlation between f_esc, HII and metallicity 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.

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, Lx, sSFR(Q(H0)). The underlying KDE represents the distribution of all the MCMC draws for galaxies with nsectors ≥ 2 only and ρ 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.

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 associated 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⁰) (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 sectors 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 emission. We expect both rotation and binary stars to affect our prediction 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_⊙) models that include binary stars produce 60% more ionizing flux compare to model without.

In Fig. 10 (panel d), we see that most of the stellar populations 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 leaking galaxies. As previously mentioned in Sect. 5.1, we note that the age and metallicity are interdependent, with larger ages tending 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 feedback is dominated by pressure and ionization mechanisms from young H II regions, older stars can also contribute through stellar winds, especially in the Wolf-Rayet phase. We note that this 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σ 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.49}^{+0.53}$ erg s⁻¹ for I Zw 18⁵. The value obtained for I Zw 18 is slightly higher but compatible with the measurement of 1 × 10⁴⁰ 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[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[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[C II]_H+ = 0 since the gas in those models is completely ionized. For the other galaxies, log f[C II]_H+ varies between −3 to −0.5 and its value correlates with the metallicity of the galaxy. This confirms previous 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²⁺ 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.

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[C II]H+ ∼ 0). The underlying KDE represents the distribution of all the MCMC draws for galaxies with Nsectors ≥ 2 only and ρ 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 ne = 30 cm−3.

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 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.

Fig. 12. KDE showing fesc, HII as a function of IR ratios. The KDE represent the distribution of all the MCMC draws for galaxies with N ≥ 2 only and ρ the associated Spearman correlation coefficient. For this, we group the MCMC draws for our selected sample by bins of fesc, HII and fit each bin with a gaussian 𝒩(x0, σ2). The blue dots and dashed lines respectively represent the gaussian mean x0 and the gaussian width σ. 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.

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.

Fig. 13. Individual SUEs showing fesc, HII as a function of IR ratios. We note that Tol 1214-277 (31) most often lies outside of the range shown in 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.

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²⁺ 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 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.

Fig. 14. KDE showing fesc, HII as a function of line ratio involving ions with high ionization potentials.

Fig. 15. Individual SUEs showing fesc, HII as a function of line ratio involving ions with high ionization potentials.

We find a correlation between f_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_TIR/L_bol ratio anti-correlates tightly with f_esc, HII. While L_bol accounts for the total luminosity produced by the central stellar cluster (and the additional X-ray contribution), L_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_bol.

Fig. 16. KDE showing fesc, HII as a function of IR and optical line ratios. The black dashed line on the second panel shows the relation O32-fesc(LyC) from Izotov et al. (2018a), assuming that [O III]λ 5007 Å/[O III]λ 4959 Å ≈ 3.

Fig. 17. Individual SUEs showing fesc, HII as a function of IR and optical line ratios.

Alternatively, we also put to test the [O III] λλ4959,5007 Å/[O II] λλ3726, 3728 Å (O32) and [S III] λλ9068,9532 Å/[S II] λλ6716,6731 Å (S32) optical line ratios. The O32 ratio has been claimed to efficiently pick up some of the most extreme leakers (Nakajima & Ouchi 2014; Izotov et al. 2016a,b, 2018a,b; Flury et al. 2022b). In particular, Flury et al. (2022b) find a significant correlation of the escape fraction detected in the UV with O32, although with large scatter. Nevertheless, the O32-f_esc relation is not straightforward; Stasińska et al. (2015) have examined galaxies with enhanced O32 ratio and find that a large fraction of them exhibit no sign of leakage. Additionally, recent simulations from Barrow et al. (2020) find that the O32 ratio fluctuates over timescales different than the dynamical time associated with escape fractions. They conclude that although high O32 tend to be associated with high f_esc, no systematic correlation is expected. Although less classically used than the O32 ratio, the S32 ratio has been shown to be a good tracer of the ionization parameter (Kewley & Dopita 2002; Mingozzi et al. 2020). Because it involves the [S II] λλ6716,6731 Å lines emission, it should also be sensitive to the presence of density-bounded regions with reduced [S II] λλ6716,6731 Å emission (Wang et al. 2019, 2021; Ramambason et al. 2020). This ratio has been used, for example, in Zastrow et al. (2011, 2013) as a diagnostic of H II region optical depths. Still, this trend has not yet been confirmed by any observational study since [S III] λλ9068,9532 Å emission has only been detected in a few LyC-leaking galaxies. Based on our models constrained by IR emission lines, we find that f_esc, HII correlates well with the predicted O32 and S32. In particular, the fit we derive for the O32-f_esc, HII relation is remarkably close to the relation from Izotov et al. (2018b) derived from fitting 11 LyC-leaking galaxies, and compatible with the new sample of 89 galaxies (detection and nondetection) observed in the Low-z Lyman Continuum Survey (Flury et al. 2022b). We further discuss the reasons that may explain those trends, and why they might be associated with a large scatter in observational studies, in the following section.

7. Discussion

7.1. From gas covering factors to H II region escape fraction

While it has been shown that the escape fraction of ionizing photons is tightly related to the structure of the neutral surrounding gas (e.g., Gazagnes et al. 2018, 2020; Chisholm et al. 2018; Saldana-Lopez et al. 2022), quantitatively linking the gas distribution to the amount of escaping photons remains a complex task. Absorption lines are powerful proxies to probe the properties of the gas on a given line of sight. In particular, H I absorption lines that saturate above a given column density have been used to derive H I covering fractions, which describe the distribution of gas on global galactic scales. Due to the gas kinematics, the covering fractions that are derived are always a lower limit to the real geometric covering fraction (e.g., Gazagnes et al. 2020). Part of this neutral gas component can also be probed in emission where specific tracers emit in dense and irradiated PDRs. 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 multisector 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_{\mathrm{esc},\mathrm{HII}}=9.{59}_{-9.51}^{+9.80}$%⁶ 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 obervations of local starbursts, 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] λλ4959,5007 Å/[S II] λλ6716,6731 Å 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⁻¹), 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 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${}_{-13.27}^{+12.29}$% for NGC 4214-c and 32.55${}_{-15.40}^{+16.44}$% 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, HII 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 cm⁻³). 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 cm⁻³). 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α 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α, [N II]/Hα and [O III]/Hβ (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α fluxes in agreement with observations. This agreement is expected since our code reproduces both the metallicity (see Fig. 5) and the recombination line Huα in the galaxies in which it is detected. Nevertheless, it appears at odds with the large f_esc, HII we derive and could be explained by different scenarios. A first possibility is that the measured Hα 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α emission in starbursts 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α 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α 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α 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 that 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, HII

7.3.1. Correlation with galactic properties and observational limitations

We find a clear anti-correlation of f_esc, HII 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.

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 fesc, 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 (LX > 1039 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: I Zw 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).

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: I Zw 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 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-know 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 ∼ 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]λ4.49 μm, [Ne VI]λ7.64 μm, and [Ne V]λ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]λ1.64 μm) boosted by X-ray photons 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λ5808, He IIλ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α 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⁰) 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/35 Z_⊙. We explain this effect by the presence of enhanced radiation fields in low-metallicity environments, which favor the ionization of C²⁺ 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.

Acknowledgments

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

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.

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.

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.

Fig. A.3. Metallicity-fesc, HII 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 fesc, HII-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%.

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.

Fig. A.5. Evolution of the volume-averaged ionization parameter log ⟨U⟩ as a function log Uin and log Uout 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 Uin as a primary parameter in our grid (see Sect. 3.2.3). We note that Uin only sets an upper limit on the actual volume averaged log ⟨U⟩. Variations in terms of metallicity and initial density n0 yield different geometries, which cover a large range of volume averaged log ⟨U⟩ for the same Uin.

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.

All Tables

All Figures

Fig. 1. Stellar (left) and X-ray (right) incident SEDs of Cloudy models for a solar metallicity and bolometric luminosity of 109 L⊙. QE shows the number of photons produced above a given energy Eν. The vertical dashed lines represent the ionization potentials of H+ (13.6 eV), O2+ (35.1 eV), O3+ (54.9 eV) and Ne4+ (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%.

In the text

	Fig. 2. Schematic view of the 17 cuts used to create sub-models.
In the text

	Fig. 3. Schematic view of the ISM of a starburst galaxy and associated representative topology.
In the text

	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.
In the text

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σ 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 α line (7–6) is detected or not.

In the text

Fig. 6. Agreement between Hα and SFR predictions with empirical measurements. Top: predicted intrinsic Hα emission vs. measured Hα corrected for Galactic extinction from (Rémy-Ruyer et al. 2015, and references therein). Bottom: predicted SFR vs. SFR(Hα+LTIR) from Rémy-Ruyer et al. (2015). The 1σ 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α 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.

In the text

	Fig. 7. Predicted SFR(Q(H0))/SFR(Hα+TIR) vs fesc, HII. The symbols represent the robust means of the posterior PDF and the gray contours show the 1σ uncertainties of the SUEs. The dashed line corresponds to the analytical curve 1/(1–fesc).
In the text

Fig. 8. Ratio of the maximum Hα emission associated with photon leakage from H II regions, assuming that all photons are reabsorbed and of the Hα value predicted by our models vs. fesc, HII. We assumed Q(H0)-to-Hα conversion coefficient of 7.31 × 1011 for a case B recombination with Te = 10 000 K and n = 100 cm−3 from Kennicutt et al. (1995) to estimate Hα(fesc, HII× Q(H0)). The Hαpred values are computed consistently in the Cloudy models. The dashed line correspond to the analytical curve fesc/(1 − fesc).

In the text

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 Nsectors ≥ 2 only and ρ 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 fesc, 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 fesc, HII for galaxies with at least two sectors.

In the text

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, Lx, sSFR(Q(H0)). The underlying KDE represents the distribution of all the MCMC draws for galaxies with nsectors ≥ 2 only and ρ 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.

In the text

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[C II]H+ ∼ 0). The underlying KDE represents the distribution of all the MCMC draws for galaxies with Nsectors ≥ 2 only and ρ 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 ne = 30 cm−3.

In the text

Fig. 12. KDE showing fesc, HII as a function of IR ratios. The KDE represent the distribution of all the MCMC draws for galaxies with N ≥ 2 only and ρ the associated Spearman correlation coefficient. For this, we group the MCMC draws for our selected sample by bins of fesc, HII and fit each bin with a gaussian 𝒩(x0, σ2). The blue dots and dashed lines respectively represent the gaussian mean x0 and the gaussian width σ. 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.

In the text

	Fig. 13. Individual SUEs showing fesc, HII as a function of IR ratios. We note that Tol 1214-277 (31) most often lies outside of the range shown in 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.
In the text

	Fig. 14. KDE showing fesc, HII as a function of line ratio involving ions with high ionization potentials.
In the text

	Fig. 15. Individual SUEs showing fesc, HII as a function of line ratio involving ions with high ionization potentials.
In the text

	Fig. 16. KDE showing fesc, HII as a function of IR and optical line ratios. The black dashed line on the second panel shows the relation O32-fesc(LyC) from Izotov et al. (2018a), assuming that [O III]λ 5007 Å/[O III]λ 4959 Å ≈ 3.
In the text

	Fig. 17. Individual SUEs showing fesc, HII as a function of IR and optical line ratios.
In the text

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 fesc, 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 (LX > 1039 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: I Zw 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).

In the text

	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.
In the text

	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.
In the text

Fig. A.3. Metallicity-fesc, HII 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 fesc, HII-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%.

In the text

	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.
In the text

Fig. A.5. Evolution of the volume-averaged ionization parameter log ⟨U⟩ as a function log Uin and log Uout 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 Uin as a primary parameter in our grid (see Sect. 3.2.3). We note that Uin only sets an upper limit on the actual volume averaged log ⟨U⟩. Variations in terms of metallicity and initial density n0 yield different geometries, which cover a large range of volume averaged log ⟨U⟩ for the same Uin.

In the text

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.

In the text