# Topological models to infer multiphase interstellar medium properties

Vianney Lebouteiller, Lise Ramambason

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

## To cite this version:

Vianney Lebouteiller, Lise Ramambason. Topological models to infer multiphase interstellar medium properties. Astronomy & Astrophysics - A&A, 2022, 10.1051/0004-6361/202243865. hal-03747406

HAL Id: hal-03747406

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

Submitted on 16 Mar 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Topological models to infer multiphase interstellar medium properties

V. Lebouteiller and L. Ramambason

Université Paris-Saclay, Université Paris-Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France
e-mail: vianney.lebouteiller@cea.fr

Received 26 April 2022 / Accepted 1 July 2022

## ABSTRACT

**Context.** Spectroscopic observations of high-redshift galaxies slowly reveal the same complexity of the interstellar medium (ISM) as expected from resolved observations in nearby galaxies. While providing, in principle, a wealth of diagnostics concerning galaxy evolution, star formation, or the nature and influence of compact objects, such high-z spectra are often spatially and spectrally unresolved, and inferring reliable diagnostics represents a major obstacle. Bright, nearby, unresolved galaxies observed in the optical and infrared domains provide many constraints to design methods to infer ISM properties, but they have so far been limited to deterministic methods and/or with simple topological assumptions (e.g., single 1D model).

**Aims.** It is urgent to build upon previous ISM multiphase and multicomponent methods by using a probabilistic approach that makes it possible to derive probability density functions for relevant parameters while also enabling a large number of free parameters with potential priors. The goal is to provide a flexible statistical framework that is agnostic to the model grid and that considers either a few discrete components defined by their parameter values and/or statistical distributions of parameters. In this paper, we present a first application with the objective to infer probability distributions of several physical parameters (e.g., the mass of H⁰, H₂, escape fraction of ionizing photons, and metallicity) for the star-forming regions of the metal-poor dwarf galaxy IZw 18 in order to confirm the low molecular gas content and high escape fraction of ionizing photons from H II regions.

**Methods.** We present a Bayesian approach to model a suite of spectral lines using a sequential Monte Carlo method provided by the Python package PyMC which combines several concepts such as tempered likelihoods, importance sampling, and independent Metropolis-Hastings chains. The algorithm, provided by the associated code MULTIGRIS, accepts a few components which can be represented as sectors around one or several stellar clusters, or continuous (e.g., power-law, normal) distributions for any given parameter. We applied this approach to a grid of models calculated with the photoionization and photodissociation code Cloudy in order to produce topological models of IZw 18.

**Results.** The statistical framework we present makes it possible to consider a large number of spectroscopic tracers, with the extinction and systematic uncertainties as potential additional random variables. We applied this technique to the galaxy IZw 18 in order to reproduce and go beyond previous topological models specifically tailored to this object. While our grid is designed for global properties of low-metallicity star-forming galaxies, we were able to calculate accurate values for the metallicity, number of ionizing photons, masses of ionized and neutral hydrogen, as well as the dust mass and the dust-to-gas mass ratio in IZw 18. We find a relatively modest amount of H₂ (~10⁵ M⊙) which is predominantly CO-dark and traced by C⁺ rather than C⁰. Nevertheless, more than 90% of the [C II] emission is associated with the neutral atomic gas. Our models confirm the necessity to include an X-ray source with an inferred luminosity in good agreement with direct X-ray observations. Finally, we investigate the escape fraction of ionizing photons for different energy ranges. While the escape fraction for the main H II region lies around 50–65%, we show that most of the soft X-ray photons are able to escape and may play a role in the ionization and heating of the circumgalactic or intergalactic medium.

**Conclusions.** Multicomponent ISM models associate a complex enough distribution of matter and phases with a simple enough topological description to be constrained with probabilistic frameworks. Despite ignoring effects such as reflected light, the diffuse radiation field, or ionization by several non-cospatial sources, they remain well adapted to individual H II regions and to star-forming galaxies dominated by one or a few H II regions, and the improvement due to the combination of several components largely compensates for other secondary effects.

**Key words.** HII regions – ISM: general – ISM: structure – galaxies: ISM – galaxies: individual: IZw18 – methods: numerical

## 1. Introduction

Spectroscopic observations of very high redshift ($z \ge 7$) galaxies have become routinely available (e.g., Harikane et al. 2018; Wagg et al. 2020; Lee et al. 2021), with an increasing number of tracers arising from different phases of the interstellar medium (ISM), enabling numerous diagnostics (e.g., star formation rate, active galactic nuclei, and molecular gas content) that had been so far limited to galaxy surveys around or after the cosmic noon ($z \le 2–3$). The integrated emission-line spectrum of galaxies, in particular, makes it possible to infer chemical

abundances, physical conditions, and relative contributions of energetic sources at work in galaxies, which in turn provide useful constraints on cosmological evolution of galaxies (e.g., Spinoglio et al. 2017) or to examine specific processes such as the escape of ionizing photons (e.g., Harikane et al. 2020) and the nature and influence of compact objects in low-metallicity environments (e.g., Lebouteiller et al. 2017).

The integrated spectrum reflects a combination of (1) many interstellar components with different properties (e.g., ionization parameter, abundances, column density, and incident radiation field) that need to be disentangled, and (2) tracers with
different emission conditions (e.g., critical density and excita-
tion temperature) that need to be consistently accounted for. It is
often assumed either that integrated measurements can be linked
to some kind of loosely defined “average” parameter or that
some regions with specific properties (e.g., H II regions) dom-
inate the galaxy emission. The modeling of integrated spectra
is thus often limited to a single cloud hypothesis (i.e., assuming
that the physical conditions of a single cloud correspond to the
average conditions of the galaxy). Specific diagnostics involving
only a few tracers are sometimes considered to isolate a specific
parameter (e.g., density measurements with optical [S II] lines,
ionization parameter, or metallicity; e.g., Kewley et al. 2019) or
to isolate a specific region or phase (e.g., low-density ionized gas
with far-IR [N II] lines; e.g., Lee et al. 2019). Nevertheless, the
derived conditions may not always represent a meaningful aver-
age galaxy value due to multiple nonlinear effects that produce
integrated emission lines.

Possible ways forward include global approaches treating all
observables coherently, such as what is done for the full spec-
tral energy distribution (SED) with, for example, MAGPHYS
(da Cunha et al. 2008) or BEAST (Gordon et al. 2016). While
many tools are available to infer physical conditions from the
SED (see review in Walcher et al. 2011), the treatment of emis-
sion lines is, however, not systematic and emission lines are
sometimes included only to decontaminate observed photome-
try bands for full SED modeling (Burgarella et al. 2005). Using
a large number of emission lines as constraints is a difficult task
due to the complex, multiphase nature of the ISM. For instance,
efforts are underway to consider photodissociation regions or
molecular gas in CIGALE (Boquien et al. 2019). Concerning
emission lines models specifically, efforts to treat many lines
at once have been mostly focused on the ionized gas, includ-
ing with modern statistical frameworks (e.g., BOND, Vale Asari
et al. 2016; NebulaBayes, Thomas et al. 2019; and WARPFIELD,
Kang et al. 2022). Still, some excitation mechanisms are often
not considered (e.g., shocks and X-rays) and a single-sector
approach is often preferred due to the quickly overwhelming
number of free parameters and the difficulty in defining the
distribution of parameters.

It is thus urgent to continue and build upon such efforts
especially since multiwavelength and multiphase ISM studies in
external galaxies have been facilitated by the advent of infrared
space instruments enabling access to neutral gas tracers such
as the Infrared Spectrograph (IRS; Houck et al. 2004) onboard
Spitzer (Werner et al. 2004), the Photodetector Array Cam-
era and Spectrometer (PACS; Poglitsch et al. 2010) onboard
Herschel (Pilbratt et al. 2010) and airborne ones such as FIFI-LS
onboard SOFIA (Fischer et al. 2018), as well as millimeter obser-
vatories that provide multiphase tracers at high-z (e.g., ALMA
and IRAM/NOEMA).

Several projects have been undergone over recent years to
infer physical parameters from IR and optical spectra of nearby
 galaxies. Studies focusing on low-metallicity dwarf galaxies
 from the Dwarf Galaxy Survey (Madden et al. 2013) pioneered
 the multiphase models of the ISM in external galaxies. We
 modeled simultaneously the ionized and neutral (atomic and
 molecular) gas with the photoionization and photodissociation
code Cloudy (Ferland et al. 2013) while considering increasingly
 complex arrangements of regions with different characteris-
tics. Since a realistic description of the ISM structure for a
given galaxy is impossible to infer, simplifications are neces-
sary. For instance, one can assume that the ionizing sources are
cospatial. Assuming galaxies are dominated by one or a few
massive H II regions, our approach has then been to consider a

stellar cluster surrounded by a small number of sectors defined
by their covering factor (fraction of the total solid angle seen
by the cluster), with the incident radiation field propagating into
each – independent – sector. The various iterations are described
in detail in Sect. 2. This approach had been devised by Péquignot
(2008), who showed that difficulties in explaining $T_e$([O III]) in
the dwarf galaxy I Zw 18 can be lifted using radiation-bounded
condensations embedded in a low-density matter-bounded dif-
fuse medium. Following the approach of Péquignot (2008), the
principle is to constrain an ISM topology, which represents the
relative contribution to the total integrated emission from dis-
crete sectors having different physical conditions. In other words,
same results will be obtained for a given topology represent-
ing different geometries (i.e., the way in which these sectors are
distributed, intertwined, or split; Sect. 4.3).

In most cases, the number of available tracers and their dif-
ferent origins (e.g., ionized/neutral and diffuse/dense phases)
require at least two sectors with different properties even when
the integrated spectrum appears to be dominated by a sin-
gle H II region. In our studies, the reflected and diffuse light
were not accounted for but the associated uncertainties are
largely compensated for by improvements enabled by a mul-
tiphase/sector approach. While the ISM topology description
above is adapted to H II region-dominated galaxies, several
efforts are underway to extend the sectors into a diffuse com-
ponent and to combine several H II regions. Our approach has
been so far to combine a rather small number of discrete sectors,
but another approach considers that the integrated emission-line
spectrum effectively corresponds to an vast ensemble of clouds
whose properties (i.e., distribution of density and distance to
source) are described and linked through a specific parameter-
ized law (e.g., power law; Richardson et al. 2014, 2016). The
observed emission is then the result of strong selection effects
due to the fact that some lines emit preferentially under some
physical conditions.

As such models become increasingly complex, it is interest-
ing to realize how they complement 3D simulations. While sim-
ulations do consider self-consistent and realistic ISM structures,
the treatment of the radiative transfer is often too calculation-
intensive. On the other hand, state-of-the-art radiative transfer
codes (e.g., Ercolano et al. 2005; Bisbas et al. 2012, 2015; Jin
et al. 2022) are difficult to scale to 3D (Morisset 2011, 2013) or
even to adapt to multiple components since the exact distribution
of matter and ionizing/energetic sources has to be predefined.
Simulations and models could eventually meet halfway pro-
vided the latter are able to consider complex geometries. As
the complexity increases (e.g., by increasing the number of
sectors, stellar populations, or excitation mechanisms), how-
ever, it becomes important to implement a statistically robust
method to identify high-probability regions in the parameter
space. Previous efforts to model galaxies with multiple sec-
tors and phases have mostly relied on $\chi^2$ methods, with several
limitations described in Sect. 2. We have thus designed a new
statistical framework to evaluate solutions in a multidimensional
grid and to quantify confidence intervals in a multiphase/sector
approach. The approach is close to NebulaBayes (Thomas et al.
2019) but here adapted to complex configurations with several
sectors with varying escape fraction of ionizing photons, that is
with high dimensionality and several continuous variables that
can be implemented in a customized environment. The tool we
present is agnostic to the grid of models used as input (i.e., it is
possible to use grids from different photoionization codes or any
grid for which sets of input parameters and output values can
be tabulated). We present here ISM applications, with a focus
on the low-metallicity dwarf galaxy IZw 18. The latter has been extensively studied and open questions include in particular the amount and distribution of the elusive molecular gas.

We summarize the various steps and results using topological models of galaxies in Sect. 2. We then describe the new probabilistic method in Sect. 3 and present specific applications to ISM studies in Sect. 4. Benchmark results are presented in Sect. 5. Finally, we illustrate the use of the new method with the dwarf galaxy IZw 18 in Sect. 6.

## 2. Modeling approach and previous results

In order to motivate the general approach presented in this study, we illustrate below recent progresses in the modeling approach with multiple phases and sectors (summarized in Fig. 1). The first step was to study spatially-unresolved emission lines in the dwarf galaxy Haro 11 (Cormier et al. 2012). We were able to construct a model with Cloudy that matches the flux of 17 IR emission lines. The full model implements the propagation of the radiation field into a single sector from a dense ionized gas into the photodissociation region (PDR). Since the model includes only one sector, the modeled “PDR covering factor” is one by definition. However, the model is constrained using only unambiguous H II region lines and the predictions for the PDR lines [C II] and [O I] are compared to the observations a posteriori to calculate an effective PDR covering factor. The model suggests that the PDR covering factor is about 10% in this galaxy. A diffuse ionized/neutral phase was included a posteriori to fully explain the emission of [N II], [Ne II], and [C II]. Structure-wise, the volume filling factors of the diffuse component, H II region, and the PDR are ≥90%, ~0.2%, and ≤0.01% respectively, hinting at a “porous” ISM (small filling factor of the dense phase) which could be linked to the relatively low-metallicity environment.

In Cormier et al. (2019), the technique was improved and applied to most of the DGS targets which were fully covered by the IRS and PACS spectrograph apertures (38 spatially-unresolved galaxies; Madden et al. 2013). The IR emission lines are modeled with two radiation-bounded sectors distinguished by their ionization parameter and densities. The PDR is assumed to be associated with the high-ionization parameter sector (i.e., the PDR is clipped from the low-ionization parameter sector). Possible PDR covering factor values, corresponding to a scaling factor of all PDR lines, are then tabulated along with other physical parameters (e.g., density and ionization parameter). Contrary to the iterative process to characterize the components needed in Haro 11, all sector parameter values are determined simultaneously. The main result shows that the PDR covering factor increases with metallicity, strengthening the result that the ISM porosity increases at low metallicity. However, it was not possible to convert this parameter into a fraction of escaping ionizing photons.

In Polles et al. (2019), we studied IR spectral cubes of the nearby (≈700 kpc) galaxy IC 10. The H II region lines are modeled with a single sector, as in Haro 11, but the sector is allowed to be matter-bounded. As such, the physical depth of the sector is another free parameter in addition to cluster age, density etc. This new parameter allowed us to determine that regions are more and more radiation-bounded as the spatial aperture over which the emission lines are integrated increases. In other words, H II regions are systematically matter-bounded and the full optically body scale is mostly radiation-bounded. While this should translate into a similar conclusion for the escape of ionizing photons, the calculation of the latter was unfortunately not possible.

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

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

c) Model: single sector with varying depth

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

d) Model #1: HII region only Model #2: PDR only

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

e) Model #1: multiple sectors with varying depth

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

**Fig. 1.** Illustration of the topological configurations considered in previous studies. For each approach, we show the relative contributions of sectors each corresponding to a specific Cloudy model. These contributions are equivalent to the fraction of the total solid angle in 3D. The light blue shell corresponds to the ionized gas, the dark gray shell to the neutral gas, and the red arc to the photoionization front. In Cormier (2012) (a), the model used an H II region connected with a PDR, with the PDR covering factor scaled a posteriori (cropped shell in white), and an additional model is considered a posteriori to account for the diffuse ionized gas. In Cormier et al. (2019) (b), two sectors are considered with the PDR associated only with the high-ionization parameter sector. In Polles et al. (2019) (c), models use a single sector with a varying depth. In Lambert-Huyghe et al. (2022) (d), observations are compared to a mix of pure H II region and PDR components drawn from the same Cloudy radiation-bounded model. In Lebouteiller et al. (2017) (e), up to four sectors were considered with different depths.

We have also examined IR emission lines in star-forming regions in the Magellanic Clouds (Lambert-Huyghe et al. 2022). The observations are spatially resolved and individual pixels are modeled with Cloudy with a single sector that can be matter-bounded. The particularity of the approach is that each Cloudy model is split between its H II region and PDR counterparts and each pixel is modeled with a mix of these. As an illustration, a
pixel can be a naked H II region (matter-bounded) or a pure PDR with no associated ionized gas emission (reflecting the fact that the ionizing sources may be located in another pixel). Any combination in between is possible, but, for the sake of this specific project, the H II region and PDR counterparts for any given pixel originate from the same Cloudy model.

Finally, the study of I Zw 18 (Lebouteiller et al. 2017) was inspired by the work of Péquignot (2008), with a conversion of the Nebu radiative transfer code to Cloudy, and with the inclusion of new IR emission lines, X-ray constraints, and dust treatment. Several sectors are used, some of them matter-bounded, and both optical and IR emission lines were used. While the model is relatively more complex compared to other efforts, the topology has not been inferred statistically but rather from a partly-manual and specific convergence process. Structure-wise, the model shows that the covering factor of the sector reaching into the molecular gas must be lower than ~0.1%, corresponding to subparsec size clumps, thereby confirming the low porosity of the ISM in an extremely metal-poor environment.

The studies described above allowed us to make considerable progress for ISM studies. As far as the technique is concerned, calibration uncertainties and other correlated uncertainties have been successfully accounted for with a covariance matrix. However, several limitations have been identified, concerning in particular the $\chi^2$ minimization technique to evaluate the solution and associated error bars. Firstly, upper limits were often ignored or else treated in a binary way (i.e., either the model is above or below the upper limit). Secondly, the number of degrees of freedom is difficult to estimate, even for linear models (Andrae et al. 2010), and some observations may trace the same physical parameter while a given observation may constrain several parameters (this difficulty increases when priors need to be considered). Degeneracy and correlation between parameters are generally difficult to account for. Due to the large number of free parameters and to the difficulty in assessing similar or degenerate solutions, it was difficult to consider combinations of more than two components (i.e., either dealing with more sectors and/or more stellar populations). Thirdly, the type of constraints (ratios and/or absolute values) and choice of tracers (ignoring or adding some tracers) often produce erratic probability density functions with solutions that are not compatible with each other within errors. Fourthly, the $\chi^2$ distribution is difficult to translate robustly into a confidence/credible interval. Finally, priors could not be used to inform the model of known quantities with uncertainties.

In the present study, we present a new probabilistic tool called MULTIGRIS<fnref n="1" /> to solve most of the issues above. The underlying objective is to derive automatically a configuration as complex as the one illustrated in Fig. 1e without the need to tabulate some parameters (in particular the mixing weight of components (i.e., their relative contributions) which can be treated as continuous variables instead. For a more general discussion of deterministic and Bayesian approaches in astrophysics we refer to Galliano (2022).

# 3. Methodology

## 3.1. General principle

The goal is to evaluate the posterior distribution of parameters from a multidimensional grid of individual precalculated models, with a series of observations (measured values with

uncertainties) as constraints. Primary parameters that define a unique set of models are used for the inference, with each parameter value drawn from a Monte-Carlo Markov Chain (MCMC). Mixing weights may be used as additional inference parameters to combine model prediction sets. Predicted observations are then computed as a linear combination of prediction sets and compared to measured values accounting for their error bars (and potentially upper/lower limits). MULTIGRIS is completely agnostic to the grid of precomputed models used as long as predicted observations can be expressed as a linear combination of prediction sets. In practice, the code can therefore perform on any kind of grid for which input (primary) parameters and output values (predictions) are tabulated.

Although one would prefer to independently compute predicted values on-the-fly for each set of drawn parameter values, the computing time may be too long for a time-efficient MCMC sampling. In our approach, the sampler explores the grid of precomputed models either with nearest neighbor or linear interpolation and with several additional continuous variables (e.g., mixing weights for the relative contribution of components or any nuisance variables). The sampling process, which is required to draw continuous random variables such as the mixing weight or any nuisance variable, is much slower than algorithms such as NebulaBayes (Thomas et al. 2019), for which likelihoods are evaluated at a predefined number of grid points, with or without interpolation. Such “brute-force” algorithms remain ideal when the dimensionality is relatively low and/or when all combination of parameters can be precalculated and tabulated. The main advantage of the sampling is to consider continuous or discrete variables without the need to tabulate them and the ability to implement new on-the-fly distributions using the same input grid.

## 3.2. Inference

The Bayesian inference is performed with PyMC<fnref n="2" /> which implements and enables gradient-based MCMC algorithms in an intuitive and relatively abstract way (Salvatier et al. 2016). The predicted value for each observation *i* is defined as:

$$M_i \sim s \sum_{c=1}^{N_{\text{comps}}} w_c f(\theta_c), \quad (1)$$

where *f* is a function that retrieves the predicted observation for any set of primary parameters from a grid using a given interpolation method, $w_c$ are the mixing weights of each component, $\theta_c$ are the primary parameters uniquely identifying a precomputed model in the grid, and the scaling factor *s* is described as a Normal distribution with mean $\mu$ and standard deviation $\sigma$:

$$s \sim N(\mu, \sigma^2), \quad (2)$$

with

$$\mu = \frac{1}{N_{\text{obs}}} \sum_{i=1}^{N_{\text{obs}}} O_i - \bar{M}_i, \quad (3)$$

where $N_{\text{obs}}$ is the total number of observations, $O_i$ is the measured value of observation *i* and $\bar{M}_i$ is the median of the predicted values for the corresponding observation, and

$$\sigma = \sqrt{\frac{1}{N_{\text{obs}} - 1} \sum_{i=1}^{N_{\text{obs}}} [O_i - \bar{M}_i - \mu]^2}. \quad (4)$$

<footnote n="2">² https://www.pymc.io/</footnote>

<footnote n="1">¹ https://gitlab.com/multigris</footnote>
<table><caption>**Table 1.** Summary of inference model parameters.</caption><thead><tr><th colspan="2">Inference model parameters</th></tr></thead><tbody><tr><td>Primary parameters for each component <i>c</i></td><td><i>θ</i><sub>c</sub></td></tr><tr><td>Mixing weight for each component <i>c</i></td><td><i>w</i><sub>c</sub></td></tr><tr><td>Scaling factor</td><td><i>s</i></td></tr></tbody></table>

By default, the model is scaled globally to match the set of observed values, making the scaling factor one of the inferred parameters (see Table 1 for the list of parameters used for inference). The use of the median of predicted values ensures that quantities do not reach large numbers while optimizing the computing time for sampling random variables.

The mixing weights are described as a Dirichlet distribution:

$$w_c \sim \mathcal{D}\left(\sum_{c=1}^{N_{\text{comps}}} w_c = 1\right). \qquad (5)$$

Primary parameters are described with Normal distributions truncated to the minimum/maximum values in the grid:

$$\theta_j \sim \min(\theta) < \mathcal{N}(\mu_j, \sigma_j^2) < \max(\theta), \qquad (6)$$

with $\mu_j$ the average parameter value in the grid and $\sigma_j^2$ spans the entire parameter value range by default for a weakly informative prior when no user-provided prior is set. The model likelihood is calculated as the product of asymmetric Student-T distributions $S$ centered on the observation for each draw of model values $M_i$:

$$\mathcal{L} = \prod_{i=1}^{N_{\text{obs}}} S(v, \mu = O_i, \sigma^2 = U_i^2), \qquad (7)$$

with $v$ the normality parameter (with a default value corresponding to a Normal distribution) and $U_i$ the asymmetric uncertainty on the observed value. The normality parameter $v$ may be changed to account for outliers, but if there are no outliers, decreasing the normality parameter effectively increases the number of high-probability regions. The model is then allowed to scan wider ranges and find that other solutions are acceptable.

The posterior probability distribution for a given model $M$ is then defined as

$$p(\theta|O, M) = \frac{p(O|\theta, M)p(\theta|M)}{p(O|M)}, \qquad (8)$$

with $p(O|\theta, M)$ the likelihood, $p(\theta, M)$ the prior probability, and $p(O, M)$ the marginal likelihood which integrates all parameter combinations:

$$p(O|M) = \int_{\theta} p(O|\theta, M)p(\theta|M)d\theta. \qquad (9)$$

For lower and upper limits we use a half-Student-T likelihood with the same normality parameter $v$ as for Eq. (7), with a default value corresponding to a Normal distribution. Systematic uncertainties on observations are treated as nuisance variables encompassing a set of one or several observations, distributed by default as a Normal distribution.

Priors on, or between, parameters are meant to be implemented in an intuitive way. Built-in priors currently allow constraining parameter values between components (e.g., value in component 1 lower than in component 2), constraining absolute parameter values as Normal distributions, upper or lower limits, and using literal equations binding parameter values to each other.

### 3.3. Input data grid

A model table is required with each row defining a unique pre-computed model prediction. The model table is then transformed to produce a grid spanning all possible combination of parameters. If some parameter combinations are not available, the invalid data is accounted for during the inference process with zero likelihood. A mask is used to track the regions where the modeled values are not finite, either due to grid completion (but also due to model upper limits).

A preprocessing script may be used to manipulate the grid before the inference step. This is useful, for instance to create new observations based on existing ones in order to be later used as constraints.

It is possible to upsample the grid, for instance with a linear interpolation. The inference (Sect. 3.2) enables on-the-fly interpolation for each MCMC draw but nonlinear interpolation that are calculation-intensive would need to be performed beforehand in order to generate a refined grid. The default grid interpolation method uses nearest neighbors, as it allows a relatively quick solution to be found. A multidimensional linear interpolation is also implemented that can be used for either all parameters or a subset. For performance issues the code samples fractional grid indices by default with nearest neighbor or linear interpolation but it is possible to force sampling real values instead.

The conversion of the initial model table into a multidimensional grid may produce many nonvalid values. For this reason, the linear interpolation often deals with an n-dimensional cell with nonfinite vertices. The current performance-driven workaround is to perform a discrete interpolation (upsample) in order to force the sampler to probe vertices instead of exploring around it. This is done by introducing a categorical random variable.

It is also expected that the grid is sampled well enough in each dimension so that the variation of the likelihood function across each grid node is accurately described with a linear interpolation (see an illustration and the discussion in Galliano et al. 2021). Potential solutions to better sample the likelihood function include computing new models possibly thanks to an adaptative grid, using a quadratic interpolation to force smooth variations of the likelihood function, or to perform a linear regression using Machine Learning techniques (Morisset et al. in prep.).

As for all grid search algorithms, there can be some edge effects. The inference model needs to probe around the high probability regions, and if the latter happen to lie on the edge of the grid or on the edge of an invalid data region (which is more likely to occur as the number of dimensions increases), the posterior mean will be slightly offset from the edge even if the maximum likelihood lies in principle at the edge. Edge effects can be mitigated in part by using the median or the mode of the chain.

Finally, the inference is performed on the primary parameters themselves by default, but it is possible to consider power-law, broken power-law, normal, or double-normal distributions instead. In such cases, the slopes and pivot point (power-law
distributions) or means and standard deviations (normal distributions) as well as lower- and upper-boundaries are described as random variables.

### 3.4. Sampling methods

Several MCMC step methods are implemented in PyMC, some of which we have benchmarked for our purpose. The Sequential Monte Carlo (SMC) sampler is a method that progresses by a series of successive annealed sequences from the prior to the posterior. The PyMC implementation is based on Ching & Chen (2007) and Minson et al. (2013). The No U-turn sampler (NUTS) works on continuous variables with Hamiltonian mechanics (Hoffman & Gelman 2011). Other methods include Metropolis-Hastings (e.g., Robert 2015) and Adaptive Differential Evolution Metropolis, the latter which uses the past chain values to inform and generate jumps (ter Braak & Vrugt 2008). Other methods may be investigated in the future.

Random walkers such as NUTS need long enough chains to authorize excursions either to or from the high-probability region. In comparison, SMC identifies the high-likelihood regions by running in parallel a large number of chains starting from random values sampled from the prior distribution. Among the many advantages of SMC, it can sample from distributions with multiple peaks and does not have a burn-in phase. Random walkers are more likely to be stuck around one peak. While several methods have been included in MULTIGRIS, the envisioned use for ISM studies implies that the probability distribution is expected to be multipeaked and we thus describe the SMC method further below.

SMC uses a combination of various techniques, with in particular tempered likelihoods and importance sampling. The tempered posterior is proportional to the product of the tempered likelihood and the prior. For a given model $M$ this translates as

$$p(\theta|O, M)_{\beta} \propto p(O|\theta, M)^{\beta} p(\theta). \quad (10)$$

SMC first samples the prior ($\beta = 0$), calculates a new $\beta$ value to match a predefined effective sample size (by default half the number of draws), computes importance weights using the tempered likelihood $p(O|\theta, M)^{\beta}$, computes a new set of samples by resampling according to importance weights, computes the mean/covariance of the proposed distribution, runs relatively short (on the order of a few to a few dozens) independent Metropolis-Hastings chains from the proposed distribution to explore the tempered posterior. The importance weights increasingly account for the likelihood as the stages increase. The stages continue until the true posterior is sampled ($\beta \ge 1$). The overall process is illustrated in Fig. 2.

Since SMC samples a large number of chains exploring the entire parameter space instead of a single chain, it is important to have enough samples to probe the prior parameter space. Furthermore, since each individual chain walks with independent Metropolis-Hastings, it is in principle possible to have only few samples per parameter but better results are obtained when using a large number of chains per parameter combination to either obtain an average probability density function (PDF) if the chains have significantly different distributions and/or to evaluate whether all chains have converged to the same distribution. While the number of samples may seem appropriate for the considered prior distribution, they also need to be numerous enough to probe the posterior distribution, which may show multiple peaks. If a small number of samples is considered, this may lead to significant stochasticity between independent model

<figure><img src="image_8.png" /><figcaption>Fig. 2. Illustration of Sequential Importance Resampling Particle Filter. From top to bottom, the distribution is sampled from the prior, reweighted and resampled according to the tempered posterior density. Surviving samples seed new Markov chains and walk through a number of Metropolis steps.</figcaption></figure>

runs. SMC runs several jobs in parallel, and in the following we always show the combined chain result. The inference result contains the posterior distribution of the parameters, which are then processed for further analysis.

### 3.5. Component identification

When the configuration requires more than one component, an important postsampling step is to identify them. By default, drawn values of primary parameters (used for the inference) are not assigned to any specific component (i.e., components are not explicitly identified or tied to any given parameter or parameter set). This is meant to assert whether inferred parameters have indeed significantly different values between each component. If parameter values are not significantly different between
<figure><img src="image_9.png" /><figcaption>**Fig. 3.** Posterior distribution of parameters for a model with two components. One parameter (age) has been purposefully forced to be the same of the two components. The raw posterior distribution is shown on the *left*, for which the two components may switch during inference, and the final result is shown on the *right*, in which the individual component parameter values have been redistributed a posteriori. For visual purposes only a subset of the parameters are shown.</figcaption></figure>

components, the components are thus allowed to switch back and forth while the solution should remain stable (specifically, for the NUTS/Metropolis samplers this means that draws may switch with time, while for the SMC sampler this means that each draw at any given stage may probe any component).

If the main objective is to derive PDFs of secondary parameters (those not used for inference and calculated in the post-processing step; Sect. 3.8), there is no need to identify explicitly each component. Nevertheless, if the parameter value for each component is difficult to disentangle, some samples may be wasted to switch from one component to the other and the probability distributions may not be as smooth as possible.

In some cases, it may be interesting to examine the posterior distribution of primary parameters for each component individually. MULTIGRIS performs a processing step to identify and characterize individual components from the chain, by redistributing components through a minimization of the standard deviation for all parameters for individual chains and components (Fig. 3). This is an iterative process, redistributing components for each draw along the entire chain several times. This redistribution bears some uncertainties and the individual component properties may still have some level of degeneracy. The global solution (combination of the components) is, however, not affected by the redistribution of components since all parameters are switched together.

Depending on the solution, it may be difficult to redistribute components a posteriori, which may indicate a need for more constraints. For this reason, it may be useful to force the identification of components a priori. In MULTIGRIS, components can be sorted using a prior constraint on any given parameter. Sorting with a prior constraint will obviously affect the inference, and the result will be different depending on which parameter is used for sorting.

## 3.6. Diagnostics

PyMC automatically returns useful statistics for each parameter (e.g., effective sample size, autocorrelation, Gelman–Rubin convergence test using multiple chains, and marginal likelihood for SMC). Several diagnostics are available in particular to check the model convergence. To calculate the final posterior distribution, we draw a random set of values from the chains (and, depending on the step method, removing the tuning steps). We also compute posterior predictive p-values for each

<figure><img src="image_10.png" /><figcaption>**Fig. 4.** Example of results for a 2-component model. The kernel density estimate (nonparametric method to estimate a PDF) is shown on top and a corner plot for one of the two components on the bottom. We note that some parameters have been forced to have identical values for the two components.</figcaption></figure>

parameter and for all parameters combined following the method in Galliano et al. (2021). Figures 4–5 provide some illustrations of some available diagnostic plots (see also Sect. 6).

## 3.7. Model selection and comparison

Model selection methods in PyMC include WAIC (widely-applicable Information Criterion; Watanabe 2010) and LOO (Leave-one-out cross-validation; Vehtari et al. 2017). LOO cross-validation, in particular, repeatedly extracts a training set to fit
<figure><img src="image_11.png" /><figcaption>Fig. 5. Example of the posterior distributions for the primary parameters (two leftmost columns) and for the observations (two rightmost columns). For each quantity, the left plot shows the PDF (with the vertical line indicating the reference measured value for observations), and the right plot shows the sequence of samples. For these plots individual jobs are combined. We note that some parameters have been forced to have identical values for the two components.</figcaption></figure>

the model and evaluates the fit using the remainder data. Approximations based on importance sampling make it possible to perform the iterations without having to refit the data.

The SMC sampler also provides the marginal likelihood $p(O|M)$ (Eq. (8)), that is the probability of the observed data given the model (where parameters are integrated out). While the $\chi^2$ asymptotically reaches a minimum value when the number of components increases, reflecting the fact that the fit does not improve significantly anymore (the reduced $\chi^2$ is difficult to estimate; Sect. 2), the marginal likelihood reaches a maximum and decreases as the number of potential solutions increases through the larger number of components (see details in Sect. 5.2). Even though increasing the number of components increases the number of potential solutions, such models become less and less likely as the overall number of parameter sets increases (spreading the prior), implying that the probability of the data given the model decreases.

To compare models, we are eventually interested in computing for each model:

$$p(M|O) \propto p(O|M)p(M). \qquad (11)$$

One can then compare models directly with different priors or models with a different number of components by using the marginal likelihood $p(O|M)$ as long as models have the same a priori probability $p(M)$ under all circumstances.

Model comparison may be used to select one preferred model or to combine models through averaging:

$$p(\theta|O) = \sum_{i=0}^{K} p(\theta|O, M_i)p(M_i|O), \qquad (12)$$

where $K$ is the number of models. An alternative to weighing through the marginal likelihood is to use information criteria methods such as WAIC or LOO with the Akaike weight for a model $i$ as

$$w_i = \frac{e^{-\frac{1}{2}dIC_i}}{\sum_{j=0}^{K} e^{-\frac{1}{2}dIC_j}}, \qquad (13)$$

with $dIC = -2 \times \text{elpd}_{\text{LOO,WAIC}}$ where elpd is the expected log pointwise predictive density for either LOO or WAIC.

### 3.8. Postprocessing

Primary parameters, that uniquely define precomputed model predictions in the grid, are the ones used for inference. From the posterior distribution of primary parameters, we can obtain the corresponding distribution for any other parameter in the grid. The postprocessing step uses the PDF of primary parameters in order to deduce secondary parameters that completely depend
on the primary parameters, or new observations that were not
used as constraints. The semantic distinction between new obser-
vations to predict and secondary parameters is only due to the
assumption that secondary parameters can never be observed.

4. Application to ISM studies

While MULTIGRIS is agnostic to the grid of precomputed mod-
els, it has been initially designed for the purpose of combining
ISM components corresponding to sectors around stellar clus-
ters or to several clusters surrounded with one or several sectors.
In that vein, the mixing weight described in Sect. 3.2 is inter-
preted as the integrated covering factor of lines of sight (or
equivalently as the fraction of the total solid angle) correspond-
ing to the same 1D model. More complex configurations are
described in Sect. 4.3. If a single cluster is considered with sev-
eral ISM components, we refer to the components as sectors in
the following.

Other configurations can also be considered, for instance
considering a large number of clouds whose parameters
are linked with a distribution (e.g., power-law and normal;
Richardson et al. 2014). In the latter cases, the distribution
parameters (e.g., slope and boundaries) need to be inferred
(instead of primary parameters for a discrete number of
clouds/components). It is of course possible to simply tabulate
the statistical cloud distribution parameters as a preprocessing
step but the ability to infer these parameters on-the-fly through
the sampling process of continuous random variables makes it
easier to consider many dimensions (e.g., varying the lower- and
upper-boundaries).

Although MULTIGRIS can consider power-law or normal
distributions, we illustrate here the applications for a multi-
sector approach building on previous works using the same
approach. Following our previous efforts, we consider that star-
forming dwarf galaxies are dominated by one or a few giant
H II regions. The ISM topology is then defined as a group of
sectors around a given cluster including thermal and possibly
nonthermal sources.

MULTIGRIS includes by default part of the publicly
available BOND database (Vale Asari et al. 2016) and a
specific model database calculated with the photoioniza-
tion/photodissociation code Cloudy (SFGX; Ramambason et al.
2022). The SFGX database is especially adapted to low-
metallicity ISM studies and includes predictions for UV to
IR emission lines and X-ray sources. This model database con-
sists of a table with each row corresponding to a specific model
outcome computed with Cloudy (Ferland et al. 2017). Primary
parameters defined in MULTIGRIS are the luminosity and age
of the stellar cluster, the luminosity and temperature of a poten-
tial X-ray source, the ionization parameter, the gas density, the
metallicity, the dust-to-gas mass ratio, and the physical depth of
the cloud. Secondary parameters include in particular masses of
dust and of various gas phases as well as the escape fraction of
ionizing photons. We briefly describe the database below with
the details available in Ramambason et al. (2022).

4.1. Cluster properties

Stellar clusters are defined by their age and luminosity, as well as
the temperature/luminosity of a potential X-ray component. The
latter is described by a multicolor blackbody parameterized by
the temperature at the inner radius of the accretion disk and by
the total luminosity.

We consider Cloudy runs with a 10<sup>7</sup> or 10<sup>9</sup> L<sub>☉</sub> cluster. A single scaling factor is used to match the observed fluxes (the scaling factor may use all observations like in Cormier et al. 2019 or a single observation, such as the total IR emission). The scaling factor together with the model luminosity provide the effective luminosity. This means that a 10<sup>10</sup> L<sub>☉</sub> galaxy is effectively considered to be a sum of a thousand 10<sup>7</sup> L<sub>☉</sub> or ten 10<sup>9</sup> L<sub>☉</sub> H II regions rather than a 10<sup>10</sup> L<sub>☉</sub> H II region. These two reference luminosities are considered to avoid variations due to shape and degeneracies between luminosity and scaling factor (filled sphere vs. thin shell; see Ramambason et al. 2022).

There is no degeneracy between the stellar luminosity and
the scaling factor used to reproduce the lines, since the relative
line fluxes do change with luminosity while the scaling factor
simply scales up or down all lines from the inferred luminos-
ity. Line ratios are not affected by the scaling factor, contrary to
topological factors like the depth of each sector (Sect. 4.2).

4.2. ISM properties

The ISM properties for each sector are defined by the density,
ionization parameter, and cloud depth. The density parameter is
set by the density at the illuminated front. The density profile is
then constant until the ionization front and then scales with the
total H column density.

An ISM parameter of interest to study the neutral gas, low-
ionization H II region species, or the escape of ionizing photons
is the physical depth of the shells. We need a general definition
for the depth cuts that applies to all models, in particular because
the scaling of parameters such as visual extinction A_V as a func-
tion of depth varies depending on metallicity. We introduce in
the SFGX model database a “cut” parameter whose variation
depends on observables, described as

- cut = 0: inner radius (illuminated edge).

- cut = 1: H⁺-H⁰ transition (i.e., ionization front; electron fraction 0.5).

- cut = 2: H⁰-H₂ transition (photodissociation front; molecular gas fraction 0.5).

- cut = 3: C⁰-CO transition (photodissociation front; molecular gas fraction 0.5).

- cut = 4: full depth corresponding to stopping criterion A_V = 10.

Fractional cut values are considered by sampling the extinc-
tion A_V profile. With this definition, radiation-bounded models
have cut ≥1. The last cut value is somewhat less trivial as
the corresponding stopping criterion is the maximum model
A_V = 10 which vary dramatically between models with different
metallicities.

4.3. ISM topology

Possible configurations enabled by MULTIGRIS include any
combination of components, either as a discrete number of sec-
tors around several clusters with the inferred mixing weight
defining the relative contribution of each sector or as a statistical
distribution of clouds parameterized by, for instance, a power-
law or a normal distribution. In the following, we focus on a
discrete number of sectors.

The contexts enabled by MULTIGRIS make it possible to
use predefined sets of priors, for instance, one cluster (defined by
the age of the stellar population and by the X-ray source proper-
ties) surrounded with *n* sectors (possibly matter-bounded), a set
of *n* clusters each associated with one sector (in practice, we mix
<figure><img src="image_12.png" /><figcaption>**Fig. 6.** Configuration of inferred sectors for an example model. Top: the left plot shows the integrated covering factors while the right plot shows the same topology but with a different, random, geometry. The color scales from yellow to red as the gas volume density increases. Different inner radii indicate different ionization parameter (U) values, with dotted arcs showing U = 0, 1, .... The gray arcs show the photoionization and photodissociation fronts. The depth is not drawn to scale. Bottom: sectors are distributed randomly in 3D according to their relative covering factors and then projected in 2D. Screening/extinction effects are ignored. Stars are added for visual purposes only.</figcaption></figure>

models with different stellar population ages and/or different X-ray source properties), or a combination of several clusters and several sectors around each cluster. When several sectors are considered, their covering factor is not tabulated a priori, contrary to Cormier et al. (2019). It is instead inferred on-the-fly through the mixing weight continuous random variable.

Figure 6 shows an illustration of inferred mixing weights. The default representation is simply a pie chart indicating the contributions of each sector that is the fraction of solid angle around the cluster that corresponds to a given sector. Another representation randomly distributes the sectors around the cluster according to their covering factor and reflects the fact that the ISM topology alone is constrained and that all the random distributions of sectors (different geometries) are thus equivalent. Finally, it is possible to illustrate a 3D distribution by creating a spherical random map to distribute all sectors based on their covering factor and by populating a line emissivity cube with each line of sight from the cube center toward a given angle corresponding to a given sector. The projected result is shown in Fig. 6.

A caveat exists when dealing with two sectors including one with negligible flux contribution for all lines (e.g., very small depth). There is then a degeneracy between the covering factor

of the main sector and the model scaling factor. This is important only for parameters that heavily depend on the topology (see also Sect. 5.3). This problem is exacerbated if using a single observation that arises mostly from one sector. For these reasons, it is advised to check whether a second sector is in fact necessary through model selection methods (Sect. 3.7).

As discussed in Ramambason et al. (2022), the overall topology and the depth of each sector may have a significant impact on the metallicity and SFR determinations, and an even stronger impact on the escape fraction of ionizing photons (see also Sect. 5.3). The escape fraction is constrained by line ratios from tracers originating at various depths and also depends on the total cluster luminosity determination (Sect. 4.1).

### 4.4. Observations

Any spectral line or photometric band can be used as observational constraint as long as Cloudy is able to compute it. The SFGX database includes several hundred tracers mostly in the optical and IR domains. While extinction-corrected fluxes can be used as constraints, MULTIGRIS may include extinction or attenuation law parameters as random variables, therefore consistently matching observations and models. This is especially useful when comparing optical and IR datasets because optical datasets only consider extinction within the layer from which optical lines arise while IR datasets may probe more embedded regions.

Nuisance variables can be added to reflect systematic uncertainties on observations. For instance, a systematic uncertainty can be considered for elemental abundances. On first order, line fluxes scale with the elemental abundance for small variations that do not modify significantly the physical conditions (e.g., cooling rates; see Cormier et al. 2019) and a systematic factor can thus be used to scale all species of a given element. Other nuisance variables may consider calibration errors between several groups of tracers.

Several observation sets can be defined with some tracer potentially appearing in several sets (e.g., observed with different instruments and requiring a different scaling factor). As new observations become available, it is therefore possible to update the posterior distributions.

## 5. Benchmark

MULTIGRIS is designed to run on standard multicore personal computers, with the ability to use large grids. The SFGX grid currently contains more than a billion values and can be read in less than a second.

With the SMC sampler, relatively short individual chains are sampled with independent Metropolis-Hastings (Sect. 3.4). Depending on the parameter variation across the grid, the minimum number of chains per parameter may range from two to a few. From our benchmark studies, good results are obtained with the SFGX context with about two chains per parameter. For a configuration with a single cluster with two or three matter-bounded sectors, this translates into a few $10^{3-4}$ samples. The inference time scales with the number of parameters to explore, and, to a somewhat lesser extent, with the number of observations due to the likelihood calculation.

### 5.1. Multimodal solutions

A common issue with ISM grids is the existence of multiple potential solutions across the grid. Random MCMC walkers are
<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Benchmark model run with a degeneracy artificially introduced for the age (idx_age_0) and metallicity (idx_Z_0). The SMC step method is shown on the *left* (each individual line corresponding to a given job) and NUTS on the *right* (each line corresponding to a given chain). The x-axes in the two columns have different ranges.</figcaption></figure>

not particularly adapted because the sampler may not be able to probe multiple peaks if the latter are too distant from each other, unless the chain is sufficiently long. On the other hand, the SMC step method initially samples from the prior and is able to probe multiple peaks through the convergence process to the posterior distribution (Sect. 3.4).

We show in Fig. 7 a model run with the SMC and NUTS samplers with similar execution time. The grid was modified to introduce a degeneracy for the age and metallicity with two equivalent solutions. SMC is able to explore the two peaks in the posterior distribution while individual chains in NUTS may be stuck to one or the other peak.

## 5.2. Number of sectors

Another common issue with multisector models is to estimate the most likely, or at least minimum, number of sectors required. Model comparison metrics help in making a choice to estimate how many components are necessary (Sect. 3.7). One can either choose a frequentist tests by computing likelihood ratios (testing how well the models explains the data using the most likely parameter set) or hypothesis tests by computing the marginal likelihood (testing how much belief one should have in each model given the data and the priors). In Cormier et al. (2019), the decision was based on the $\chi^2$ value but this only indicated the

minimum number of sectors to obtain a satisfactory agreement between observations and models, that is from a pure frequentist point of view.

In principle, the marginal likelihood reaches a maximum for the most likely number of sectors given the data and the priors. In other words, a higher number of sectors is less likely to reproduce the data because the latter is overfitted. In that vein, the Bayes factor, which computes the marginal likelihood ratio is a useful metrics, but the actual probability we are interested in to compare models is $p(M|O)$, which supposes the knowledge of the overall a priori probability of the model, $p(M)$ (Sect. 3.7).

Another model comparison method is to run a model with $n$ sectors and infer whether up to $n-1$ sectors can be ignored by setting their covering factor to zero during inference. The $n-1$ sectors that can be ignored are described with a categorical random variable. This allows for an on-the-fly estimate of the most likely number of components, with the advantage of inferring a probability for each number of components. On important caveat though is that this method is adapted only in cases when potential additional sectors do not change significantly the properties of the other sectors in order to prevent erratic jumps during inference.

In Fig. 8, we show a benchmark test for estimating the number of sectors. We start from a model with a given number of sectors (from one to three) and run the MCMC model fixing the number of sectors from one to three. We also show the results when letting the number of sectors free (but limited to two or three). One can see that all metrics (LOO, WAIC, SMC marginal likelihood) reach a maximum for the expected number of sectors. When the number of sectors is itself allowed to vary, the expected number is inferred (obviously as long as the maximum number of sectors is enough). We note that such methods are reliable only if (1) all models have passed the burn-in phase (for random walkers), (2) all chains have converged toward an identical solution, and (3) there are significantly more observations that parameters. Furthermore, we made the implicit assumption that there is no intrinsic reason to prefer any specific configuration for all possible combinations of parameters and data (i.e., $p(M)$ is the same for all models).

While the benchmark test performs as expected, the application for ISM components is much less straightforward. Models are supposed to represent any given number of sectors around a single cluster or possibly even several clusters (Sect. 4.3). However, depending on the complexity of the dataset (e.g., observations coming from different phases) and the complexity of the studied object (e.g., single H II region or a full galaxy) the model may not be fully adapted a priori, resulting in posterior probabilities that need to be regarded with caution. Moreover, models may not have the same a priori probability $p(M)$. For a multisector approach, $p(M)$ can be somewhat subjective and cannot be computed rigorously since we know the model does not correspond to the reality of a galaxy or a stellar cluster. At the very least, $p(M)$ is expected to be low for a single sector with a given density/pressure profile.

The systematic uncertainties related to the model description as a few sectors may motivate the use of more statistically appropriate models with a large number of clouds whose physical conditions can be parameterized. In that vein, it is important to consider and compare various types of distributions (Sect. 6).

## 5.3. Influence of topology on parameters

While it is possible to compare models with different number of sectors based on the available observations (Sect. 5.2), it must be
<figure><img src="image_14.png" /><figcaption>**Fig. 8.** Model comparison to estimate the number of sectors required. We show models with an input of one, two, and three sectors from left to right. In the top panel we show the log values of LOO, WAIC, and the SMC marginal likelihood ("lml") for MCMC runs with a number of sectors from one to three. In the middle panel we show the likelihood weights (not to be confused with the mixing weight) of each sector when the number of sectors is allowed to vary. Each horizontal line corresponds to a model with the number of sectors allowed to vary up to two or three. For each horizontal line, the filled circles show the inferred weight and the circled one highlights the number of sectors with the largest weight. In the bottom panel we show the fraction of draws that agree with the observed value within 2- and 3-σ.</figcaption></figure>

noted that some parameters are particularly sensitive to the num-
ber of sectors used. First and foremost, it is generally assumed
here that all sectors contribute significantly to the emission of
at least one tracer. This is to avoid having to consider an empty
or unconstrained sector (model with holes), which would effec-
tively change the covering factors of the other sectors. If we were
to assume empty sectors, it is possible to calculate the effec-
tive total covering factor of nonempty sectors provided the input
stellar cluster luminosity (i.e., the bolometric luminosity) is con-
strained independently (Sect. 4.3). In the following we always
assume only nonempty sectors.

Figure 9 shows the variation of inferred parameters as a func-
tion of the number of sectors used for the same series of typical
benchmark models as in Sect. 5.2. When the true solution cor-
responds to a single sector, using a model with two or three
sectors does not change significantly the parameter values except
for the escape fraction of ionizing photons. However, such mod-
els become more and more unlikely as the number of sectors
increases. Similarly, using three sectors for a true solution with
two sectors leads to a good agreement with expected results, this
time even for the escape fraction. More generally, using a larger
than necessary number of sectors usually does not impact signifi-
cantly the results except for the escape fraction in some cases. On
the other hand, using a lower than necessary number of sectors
may result in significant differences for all parameters, but keep-
ing in mind that the model is likely unsatisfactory (Fig. 8). In
summary and unsurprisingly, the parameters which depend the
most on the topology are the ones which depend on the depth of
each sector, that is the escape fraction of ionizing photons but
also the mass of H₂.

5.4. *Selection of observations*

In Figs. 10–11 we show a benchmark test analyzing the evolu-
tion of the solution with an increasing number of lines. For this

illustration, the series of lines we consider is sorted based on
their detection or potential detection in high-z sources. The first
set of three lines includes [C II], [O III] 88 µm, and [O I] 63 µm.
We then add [N II] 205 µm and [O I] 145 µm, then [Ne II]
12.8 µm, [Ne III] 15.6 µm, [N III] 57 µm and [Ar II] 7 µm, then
[Ar III] 9 µm, H₂ 17 µm, [Si II] 34 µm, and [Ne V] 14 µm. As the
figures show, the solution naturally improves with the number
of lines considered. In the radiation-bounded case (eight param-
eters; Fig. 10), a satisfactory solution is found even with three
lines. For the matter-bounded case (ten parameters), for which
there is an added degeneracy due to the inclusion of the depth
parameter, a satisfactory solution is found with nine lines.

While it may seem that considering as many lines as pos-
sible would be the best option, there are several important
caveats to take into account when dealing with actual obser-
vations. Since one is looking for an approximate/simplified
topological model, it is may be somewhat illusory to search a
model that will reproduce many observations at once. There is
a risk of diluting the actual solution if the configuration is not
adapted and/or if the systematic uncertainties (e.g., instrumen-
tal, atomic data, and abundances) are not properly accounted for.
For instance, a model may be controlled by the likelihoods of
a subset of many observations corresponding to lines arising in
similar physical conditions – though in slight disagreement – at
the expense of other few observations constraining potentially
interesting/complementary parameters. Furthermore, outliers or
observations with abnormally small error bars may produce a
pathological posterior with local probability peaks where the
sampler gets stuck (either preventing excursions for random
walkers or affecting the computation of important samples for
SMC). Unfortunately, this prevents the sampler from explor-
ing other regions in the parameter space and, as a result,
the sampler may be stuck in a global (all observations com-
bined) low-probability region. Inversely, this also means that a
global high-probability region may be dominated by an outlier
<figure><img src="image_15.png" /><figcaption>**Fig. 9.** Output parameters for the benchmark model comparison. We show models with an input of one, two, and three sectors and for radiation-bounded ("rb") or matter-bounded ("mb") template models. Each plot shows a standard boxplot configuration with the median value in red, with the box spanning from the first to third quartile, and the horizontal gray bars extending to 1.5 times the inter-quartile range. Faint open circles show the data that extend beyond the horizontal gray bars. In each panel we show the results forcing the number of components from one to three. *The two top panels* show the marginal likelihood ("lml") and the likelihood while the *other panels below* show several parameters (*from top to bottom* the number of ionizing photons, the mass of H⁺, the mass of H₂ traced by CO, the mass of H₂ traced by C⁺, the X-ray luminosity, the metallicity, and the escape fraction of ionizing photons). The horizontal lines show the expected value, except for the marginal likelihood where it indicates the highest value.</figcaption></figure>

or and observation with an abnormally small uncertainty. The $\hat{r}$ diagnostic (comparing means and variances of chains; Vehtari et al. 2021) may be useful to examine whether one observation dominates the likelihood calculation.

Since these issues intrinsically depend on which observations are considered and with what uncertainties, the Bayesian approach by itself does not perform necessarily better than the $\chi^2$ approach and great care must be taken to adapt the model

configuration to the set of observations and to account for all possible systematic uncertainties. Nevertheless, we can use the PDF to analyze the dependencies between parameters and predicted observations. Once the sampler converges, the values drawn according to the primary parameter PDFs make it possible to investigate correlations between parameters, between parameters and predicted observations, and between predicted observations themselves. With this in mind, it may be useful to
<figure><img src="image_16.png" /></figure>

**Fig. 10.** Evolution of the solution with the number of lines considered for a benchmark run. Each row corresponds to a parameter, and each column represents a different number of lines from three to 13 ([C II]+[O III] 88 µm+[O I] 63 µm, +[N II] 205 µm+[O I] 145 µm, +[Ne II] 12.8 µm+[Ne III] 15.6 µm+[N III] 57 µm+[Ar II] 7 µm, +[Ar III] 9 µm+H₂ 17 µm+[Si II] 34 µm+[Ne V] 14 µm. The horizontal black lines indicate the expected solution. The model used here is radiation-bounded.

include all available observations at first and, (1) identify observations that are not well reproduced in order to consider potential missing systematic uncertainties and/or incorrect model assumptions (e.g., number of parameters and specific configuration), (2) check if some observations effectively constrain at least one parameter, and inversely, if every parameter is constrained by at least one observation, and (3) check if some observations constrain the same parameter(s), in which case redundant observations may be removed to prevent that the likelihood is controlled only by few parameters without changing significantly the solution.

Figure 12 shows the map of the correlation coefficients between predicted observations for a given run. For each set of two observations, we compare the correlation coefficients with each parameter for all MCMC draws (i.e., around the solution). The final correlation coefficient between the two predicted observations is zero if the two predicted observations constrain entirely different parameters and one if the two predicted observations constrain precisely the same parameters. The correlations do not account for the observed values/uncertainties and are calculated using the input grid values. Nevertheless, they may help in choosing an optimal set of observations.

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

**Fig. 11.** Same as Fig. 10, but for a matter-bounded model.

**6. Modeling of the low-metallicity star-forming dwarf galaxy IZw 18**

The low-metallicity star-forming dwarf galaxy IZw 18 has been extensively studied with topological models in Péquignot (2008, hereafter P08) and Lebouteiller et al. (2017, hereafter L17). Among the most pressing questions concerning this object are the amount and distribution of the molecular gas and the influence of the known bright X-ray source on the ISM properties.

We use the SFGX grid (Sect. 4) with a linear interpolation for the metallicity and luminosity parameters. The interpolation is particularly useful for the luminosity which is not well sampled with only two values (Sect. 4.1). Output parameter values are provided as the median and highest density intervals with 94% confidence level. The full list of observations used is shown in Table 2. We note that the observed values for IR lines have been completely revised in L17 compared to P08, but these lines had not been used for the convergence process in the latter. Furthermore, the fluxes used here are tailored to the study of the NW H II region in IZw 18 and differ slightly from the uniform dataset used for the DGS sample in Ramambason et al. (2022).
<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Map of the correlation coefficients between predicted observations for a benchmark run. Here one can see that, according to the input grid, most H lines trace the same parameters around the solution (the same applies to He or [S II]).</figcaption></figure>

**6.1. Reproducing the modeling strategy of Péquignot (2008) with optical lines: Model $M_{\text{I,opt}}$**

We first build model $M_{\text{I,opt}}$ to illustrate how MULTIGRIS may reproduce the results in P08 with the use of a general photoionization grid. P08 used an iterative and partly-manual convergence process specifically tailored to IZw 18. The observed spectrum used in P08 is close to Izotov et al. (1999) but includes a correction from $E(B-V) = 0.08$ to 0.04, which in turn provides a different value and a constraint on the $H\alpha/H\beta$ ratio as compared to Izotov et al. (1999).

While P08 used 12.97 Mpc for the distance, we use here 18.2 Mpc (Aloisi et al. 2007) to be consistent with L17. Predicted physical quantities in P08 are modified accordingly here. Contrary to P08 we use error bars as they are naturally required to produce random draws in the MCMC process. In order to prevent the sampler from being stuck in local high-probability regions (Sect. 5.4), we introduce a default minimum 10% error on all lines.

Instead of using the [S II] line ratio as a specific observation to constrain the density around the ionization front as in P08, we use the absolute fluxes of all lines including both [S II] lines. The lines used for inference depend on the densities at various depths, which all translate into a constraint for the density at the illuminated front, which is the only parameter defining the density (Sect. 4.2). While we could use the known X-ray luminosity (Kaaret & Feng 2013) and dust-to-gas mass ratio (Galliano et al. 2021) in IZw 18 as observed values together with the spectral line fluxes, we chose to predict them instead and compare them with observations a posteriori.

The helium abundance in the SFGX database is $He/H=0.098$ and does not scale with metallicity. We tentatively scale all helium line intensities using a nuisance variable with the assumption that (1) on first order all line fluxes scale with the abundance as long as the scaling factor is small enough

<figure><img src="image_19.png" /><figcaption>**Fig. 13.** Comparison between observed and modeled fluxes for the ionized gas model $M_{\text{I,opt}}$ of IZw 18. Plotted values in red are the difference between observations and models normalized by the model value. MULTIGRIS produces this plot in cartesian and polar coordinates, we show here the latter for illustration and because it is better adapted to a large number of observations. The black circle, black dots, and orthogonal black lines indicate the model and the red points and red error bars indicate the observations. The color of the sector is green if observations and models agree within errors and yellow if the model agrees within the observed error bar.</figcaption></figure>

and (2) the chemical and physical conditions of the model are not significantly impacted by the abundance variation. In all results below, we find He/H between 0.07 and 0.08 for all runs, compared to 0.08 used in P08. Helium is the only element for which we force the abundance. The other elemental abundances scale with O/H depending on several prescriptions described in Ramambason et al. (2022).

Since our goal is to compare to P08, we first use a similar three-sector configuration with an identical inner radius for all sectors, with no prior. We allow the cloud depth to reach slightly beyond the ionization front in order to fully consider the ionized layer. The primary parameters that we use for inference as well as predicted (secondary) parameters are listed in Table 3.

We show the comparison between the observations and the predictions of model $M_{\text{I,opt}}$ in Table 2 and Fig. 13. Results for parameters are shown in Table 4. P08 derived one sector as a matter-bounded “veil” with covering factor (CF) 44% to account for the fact that absorbing gas is expected for any line of sight from the cluster, best corresponding to our sector #1, one sector as a matter-bounded sector with CF = 30%, best corresponding to our sector #2, and one sector as a radiation-bounded sector with CF = 26% best corresponding to our sector #3. The inner radius in model $M_{\text{I,opt}}$ is $\approx 6_{-5}^{+25} \times 10^{20}$ cm as compared to a fixed value of $5.6 \times 10^{20}$ cm in P08 (18.2 Mpc). Sector #1 contributes little to the total line emission but contributes significant fraction of He II 4686Å, which is the brightest line in this sector (Fig. 14), as in P08.

We find good agreements for the number of ionizing photons $Q_0$ and the ionized gas mass $M(H^+)$ compared to P08. The escape fraction is calculated using the ionizing continuum output from Cloudy and is also in good agreement, although it is by far the parameter with the largest relative uncertainty. In P08,
<table><caption>**Table 2.** Observations and predictions for the ionized and neutral gas models of IZw 18.</caption><thead><tr><th>Observations</th><th>Observed value</th><th>M<sub>Iopt</sub></th><th>M<sub>Iopt+IR</sub></th><th>M<sub>NA</sub></th><th>M<sub>NM</sub></th></tr></thead><tbody><tr><td>S26716.44A</td><td>22.50 ± 7.00</td><td>16.06<sup>+4.62</sup><sub>-3.71</sub></td><td>15.29<sup>+4.49</sup><sub>-3.08</sub></td><td>17.48<sup>+5.50</sup><sub>-3.38</sub></td><td>18.27<sup>+3.86</sup><sub>-4.89</sub></td></tr><tr><td>S26730.82A</td><td>16.90 ± 7.00</td><td>12.61<sup>+2.85</sup><sub>-2.70</sub></td><td>11.88<sup>+2.39</sup><sub>-3.07</sub></td><td>13.66<sup>+4.18</sup><sub>-3.04</sub></td><td>14.29<sup>+3.05</sup><sub>-3.16</sub></td></tr><tr><td>H14101.73A</td><td>266.00 ± 26.60</td><td>250.68<sup>+41.85</sup><sub>-25.57</sub></td><td>249.06<sup>+22.45</sup><sub>-27.66</sub></td><td>250.26<sup>+22.62</sup><sub>-27.61</sub></td><td>241.69<sup>+22.96</sup><sub>-26.11</sub></td></tr><tr><td>H14340.46A</td><td>461.00 ± 46.10</td><td>487.67<sup>+43.61</sup><sub>-60.95</sub></td><td>486.39<sup>+57.38</sup><sub>-40.49</sub></td><td>487.33<sup>+51.74</sup><sub>-48.84</sub></td><td>477.37<sup>+51.21</sup><sub>-45.53</sub></td></tr><tr><td>H14861.33A</td><td>1000.00 ± 100.00</td><td>1019.23<sup>+87.81</sup><sub>-125.72</sub></td><td>1015.64<sup>+119.14</sup><sub>-108.84</sub></td><td>1017.26<sup>+113.61</sup><sub>-102.38</sub></td><td>1002.37<sup>+104.68</sup><sub>-101.70</sub></td></tr><tr><td>H16562.81A</td><td>2860.00 ± 286.00</td><td>2901.65<sup>+278.02</sup><sub>-366.10</sub></td><td>2958.23<sup>+286.61</sup><sub>-342.18</sub></td><td>3089.73<sup>+344.61</sup><sub>-357.65</sub></td><td>2993.75<sup>+351.15</sup><sub>-292.05</sub></td></tr><tr><td>He14471.49A</td><td>21.40 ± 2.14</td><td>24.60<sup>+3.06</sup><sub>-3.64</sub></td><td>23.94<sup>+3.26</sup><sub>-3.30</sub></td><td>23.25<sup>+3.78</sup><sub>-2.43</sub></td><td>24.08<sup>+3.27</sup><sub>-3.10</sub></td></tr><tr><td>He15875.64A</td><td>67.70 ± 6.77</td><td>68.73<sup>+7.09</sup><sub>-10.57</sub></td><td>66.95<sup>+7.86</sup><sub>-8.92</sub></td><td>65.13<sup>+8.88</sup><sub>-7.84</sub></td><td>66.42<sup>+7.29</sup><sub>-9.12</sub></td></tr><tr><td>He16678.15A</td><td>25.30 ± 2.53</td><td>18.60<sup>+2.44</sup><sub>-2.73</sub></td><td>18.01<sup>+2.69</sup><sub>-2.19</sub></td><td>17.43<sup>+2.29</sup><sub>-2.53</sub></td><td>17.74<sup>+2.16</sup><sub>-2.28</sub></td></tr><tr><td>He17065.22A</td><td>24.40 ± 2.44</td><td>22.76<sup>+4.29</sup><sub>-4.10</sub></td><td>22.10<sup>+4.63</sup><sub>-3.99</sub></td><td>22.82<sup>+5.47</sup><sub>-3.79</sub></td><td>22.70<sup>+4.82</sup><sub>-3.73</sub></td></tr><tr><td>He24685.64A</td><td>36.80 ± 3.68</td><td>35.59<sup>+8.56</sup><sub>-10.48</sub></td><td>38.49<sup>+7.11</sup><sub>-10.18</sub></td><td>34.85<sup>+9.22</sup><sub>-10.53</sub></td><td>35.54<sup>+8.53</sup><sub>-9.30</sub></td></tr><tr><td>N26548+84A</td><td>9.20 ± 2.00</td><td>9.81<sup>+2.95</sup><sub>-2.83</sub></td><td>9.18<sup>+2.88</sup><sub>-2.62</sub></td><td>9.06<sup>+4.16</sup><sub>-2.53</sub></td><td>9.93<sup>+2.84</sup><sub>-2.66</sub></td></tr><tr><td>O23726+9A</td><td>238.00 ± 23.80</td><td>217.80<sup>+36.13</sup><sub>-4.13</sub></td><td>213.20<sup>+44.43</sup><sub>-39.32</sub></td><td>202.23<sup>+51.98</sup><sub>-48.84</sub></td><td>207.29<sup>+43.47</sup><sub>-43.02</sub></td></tr><tr><td>O35007+4959A</td><td>2683.00 ± 268.30</td><td>2541.64<sup>+419.60</sup><sub>-372.86</sub></td><td>2382.20<sup>+314.02</sup><sub>-274.36</sub></td><td>2407.67<sup>+498.47</sup><sub>-277.28</sub></td><td>2446.83<sup>+293.20</sup><sub>-348.72</sub></td></tr><tr><td>Ne33869+968A</td><td>191.00 ± 19.10</td><td>204.78<sup>+37.24</sup><sub>-29.07</sub></td><td>192.99<sup>+21.14</sup><sub>-25.22</sub></td><td>195.72<sup>+20.34</sup><sub>-26.47</sub></td><td>194.82<sup>+26.52</sup><sub>-22.43</sub></td></tr><tr><td>Ar37135+751A</td><td>23.50 ± 2.35</td><td>22.90<sup>+3.82</sup><sub>-4.10</sub></td><td>23.54<sup>+2.87</sup><sub>-2.86</sub></td><td>22.97<sup>+2.32</sup><sub>-4.18</sub></td><td>22.21<sup>+2.89</sup><sub>-3.07</sub></td></tr><tr><td>Fe34658.01A</td><td>4.50 ± 2.00</td><td>7.08<sup>+1.78</sup><sub>-1.86</sub></td><td>7.12<sup>+1.54</sup><sub>-1.90</sub></td><td>6.70<sup>+2.00</sup><sub>-1.44</sub></td><td>6.30<sup>+1.46</sup><sub>-1.79</sub></td></tr><tr><td>O388.3323m</td><td>200.00 ± 30.00</td><td>(187.20)<sup>+41.45</sup><sub>60.46</sub></td><td>174.74<sup>+41.43</sup><sub>-30.17</sub></td><td>175.54<sup>+31.36</sup><sub>-31.22</sub></td><td>166.47<sup>+31.80</sup><sub>-41.51</sub></td></tr><tr><td>H112.3684m</td><td>8.50 ± 3.00</td><td>(6.38)<sup>+0.82</sup><sub>-0.99</sub></td><td>6.30<sup>+0.60</sup><sub>-0.76</sub></td><td>6.25<sup>+0.65</sup><sub>-0.73</sub></td><td>6.00<sup>+0.68</sup><sub>-0.67</sub></td></tr><tr><td>O425.8832m</td><td>38.00 ± 8.00</td><td>(35.36)<sup>+18.04</sup><sub>-21.47</sub></td><td>27.63<sup>+6.75</sup><sub>-9.73</sub></td><td>27.58<sup>+8.68</sup><sub>-9.89</sub></td><td>27.51<sup>+8.95</sup><sub>-8.54</sub></td></tr><tr><td>Ne212.8101m</td><td>6.00 ± 3.00</td><td>(2.36)<sup>+0.83</sup><sub>-0.85</sub></td><td>2.31<sup>+0.69</sup><sub>-0.56</sub></td><td>4.62<sup>+3.30</sup><sub>-1.93</sub></td><td>5.03<sup>+3.36</sup><sub>-2.35</sub></td></tr><tr><td>Ne315.5509m</td><td>47.00 ± 10.00</td><td>(51.59)<sup>+8.90</sup><sub>-7.10</sub></td><td>49.24<sup>+4.96</sup><sub>-6.39</sub></td><td>49.82<sup>+5.60</sup><sub>-5.72</sub></td><td>49.27<sup>+6.10</sup><sub>-5.29</sub></td></tr><tr><td>Ne336.0036m</td><td>&lt;7.00</td><td>(4.67)<sup>+0.87</sup><sub>-0.70</sub></td><td>4.43<sup>+0.45</sup><sub>-0.58</sub></td><td>4.48<sup>+0.49</sup><sub>-0.53</sub></td><td>4.43<sup>+0.55</sup><sub>-0.48</sub></td></tr><tr><td>Ne514.3228m</td><td>&lt;2.30</td><td>(3.02)<sup>+8.75</sup><sub>-3.02</sub></td><td>0.93<sup>+1.14</sup><sub>-0.92</sub></td><td>1.50<sup>+1.42</sup><sub>-1.17</sub></td><td>1.58<sup>+1.17</sup><sub>-1.13</sub></td></tr><tr><td>Ne524.2065m</td><td>&lt;4.00</td><td>(3.61)<sup>+10.13</sup><sub>-3.61</sub></td><td>1.10<sup>+1.33</sup><sub>-1.10</sub></td><td>1.83<sup>+1.59</sup><sub>-1.41</sub></td><td>1.87<sup>+1.26</sup><sub>-1.39</sub></td></tr><tr><td>S318.7078m</td><td>28.00 ± 7.00</td><td>(37.01)<sup>+6.77</sup><sub>-6.25</sub></td><td>37.98<sup>+4.39</sup><sub>-4.75</sub></td><td>37.85<sup>+3.69</sup><sub>-5.15</sub></td><td>37.12<sup>+4.41</sup><sub>-4.44</sub></td></tr><tr><td>S333.4704m</td><td>40.00 ± 7.00</td><td>(50.64)<sup>+13.24</sup><sub>-9.58</sub></td><td>50.15<sup>+8.59</sup><sub>-6.49</sub></td><td>51.27<sup>+8.50</sup><sub>-7.77</sub></td><td>49.90<sup>+8.26</sup><sub>-6.55</sub></td></tr><tr><td>S410.5076m</td><td>49.00 ± 10.00</td><td>(65.65)<sup>+29.61</sup><sub>-21.61</sub></td><td>53.05<sup>+10.92</sup><sub>-9.87</sub></td><td>56.33<sup>+9.74</sup><sub>-11.53</sub></td><td>56.86<sup>+11.20</sup><sub>-10.97</sub></td></tr><tr><td>Ar26.98337m</td><td>&lt;5.00</td><td>(0.22)<sup>+0.13</sup><sub>-0.11</sub></td><td>0.18<sup>+0.06</sup><sub>-0.06</sub></td><td>0.47<sup>+0.45</sup><sub>-0.24</sub></td><td>0.57<sup>+0.42</sup><sub>-0.33</sub></td></tr><tr><td>Ar321.8253m</td><td>&lt;8.00</td><td>(0.60)<sup>+0.11</sup><sub>-0.11</sub></td><td>0.60<sup>+0.07</sup><sub>-0.08</sub></td><td>0.58<sup>+0.09</sup><sub>-0.09</sub></td><td>0.56<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>Ar38.98898m</td><td>&lt;10.00</td><td>(8.25)<sup>+1.62</sup><sub>-1.27</sub></td><td>8.45<sup>+0.76</sup><sub>-1.37</sub></td><td>8.19<sup>+1.15</sup><sub>-1.24</sub></td><td>7.92<sup>+1.07</sup><sub>-1.09</sub></td></tr><tr><td>Fe322.9190m</td><td>&lt;6.00</td><td>(3.75)<sup>+0.89</sup><sub>-1.23</sub></td><td>3.53<sup>+1.05</sup><sub>-1.00</sub></td><td>3.35<sup>+1.30</sup><sub>-0.78</sub></td><td>3.12<sup>+0.96</sup><sub>-0.87</sub></td></tr><tr><td>LTIR500.500m</td><td>2808.40 ± 2808.40</td><td>(5425.82)<sup>+36065.1</sup><sub>-2509.9</sub></td><td>(16977.40)<sup>+184385.2</sup><sub>-13675.2</sub></td><td>28190.57<sup>+7133.6</sup><sub>-7615.3</sub></td><td>26628.15<sup>+9958.7</sup><sub>-7175.0</sub></td></tr><tr><td>C2157.636m</td><td>76.00 ± 20.00</td><td>(1.28)<sup>+2.86</sup><sub>-1.13</sub></td><td>(1.95)<sup>+1.88</sup><sub>-1.51</sub></td><td>27.32<sup>+19.42</sup><sub>-23.75</sub></td><td>29.30<sup>+30.18</sup><sub>-24.53</sub></td></tr><tr><td>O163.1679m</td><td>60.00 ± 20.00</td><td>(0.25)<sup>+0.26</sup><sub>-0.25</sub></td><td>(1.36)<sup>+26.86</sup><sub>-1.05</sub></td><td>77.21<sup>+44.25</sup><sub>-41.66</sub></td><td>85.29<sup>+46.83</sup><sub>-41.94</sub></td></tr><tr><td>Si234.8046m</td><td>65.00 ± 12.00</td><td>(7.69)<sup>+6.00</sup><sub>-4.40</sub></td><td>(4.79)<sup>+4.95</sup><sub>-1.64</sub></td><td>52.70<sup>+27.57</sup><sub>-17.67</sub></td><td>54.34<sup>+27.82</sup><sub>-20.43</sub></td></tr><tr><td>Fe217.9314m</td><td>&lt;3.00</td><td>(0.38)<sup>+0.39</sup><sub>-0.36</sub></td><td>(0.05)<sup>+0.51</sup><sub>-0.04</sub></td><td>0.79<sup>+1.10</sup><sub>-0.55</sub></td><td>0.94<sup>+0.74</sup><sub>-0.78</sub></td></tr><tr><td>Fe225.9811m</td><td>13.00 ± 8.00</td><td>(0.80)<sup>+0.92</sup><sub>-0.77</sub></td><td>(0.11)<sup>+0.94</sup><sub>-0.08</sub></td><td>6.35<sup>+5.62</sup><sub>-2.57</sub></td><td>6.61<sup>+3.63</sup><sub>-3.71</sub></td></tr><tr><td>CO2600.05m</td><td>1.50 ± 0.50 × 10<sup>-04</sup></td><td>(7.79)<sup>+3.2×10<sup>-11</sup></sup><sub>7.78</sub></td><td>(1.19)<sup>+4×10<sup>-03</sup></sup><sub>1.12</sub></td><td>5.06<sup>+1.4×10<sup>-03</sup></sup><sub>-5.05</sub></td><td>7.05<sup>+×10<sup>-04</sup></sup><sub>-7.05</sub></td></tr><tr><td>H22.12099m</td><td>&lt;2.00</td><td>(2.02)<sup>+1.2×10<sup>-14</sup></sup><sub>-1.89</sub></td><td>(1.16)<sup>+18.77</sup><sub>-0.65</sub></td><td>4.52<sup>+10.76</sup><sub>-4.52</sub></td><td>0.10<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>H228.2130m</td><td>&lt;7.00</td><td>(2.02)<sup>+1.2×10<sup>-14</sup></sup><sub>-1.89</sub></td><td>(1.16)<sup>+18.77</sup><sub>-0.65</sub></td><td>1.29<sup>+128.18</sup><sub>-1.29</sub></td><td>0.28<sup>+5.37</sup><sub>-0.28</sub></td></tr><tr><td>H212.2752m</td><td>&lt;7.50</td><td>(2.02)<sup>+1.2×10<sup>-14</sup></sup><sub>-1.89</sub></td><td>(1.16)<sup>+21.41</sup><sub>-0.60</sub></td><td>0.14<sup>+0.87</sup><sub>-0.14</sub></td><td>0.59<sup>+0.91</sup><sub>-0.59</sub></td></tr><tr><td>H217.0300m</td><td>&lt;7.00</td><td>(2.02)<sup>+1.2×10<sup>-14</sup></sup><sub>-1.89</sub></td><td>(1.16)<sup>+18.90</sup><sub>-0.65</sub></td><td>0.17<sup>+1.56</sup><sub>-0.17</sub></td><td>0.71<sup>+2.61</sup><sub>-0.71</sub></td></tr><tr><td>O16300+63A</td><td>8.50 ± 2.00</td><td>(1.62)<sup>+1.99</sup><sub>-1.57</sub></td><td>(0.19)<sup>+1.91</sup><sub>-0.16</sub></td><td>(3.54)<sup>+3.39</sup><sub>-2.14</sub></td><td>(4.37)<sup>+2.26</sup><sub>-3.33</sub></td></tr><tr><td>O27320+30A</td><td>6.30 ± 2.00</td><td>(4.91)<sup>+1.63</sup><sub>-1.39</sub></td><td>(4.98)<sup>+2.29</sup><sub>-0.90</sub></td><td>(4.76)<sup>+1.51</sup><sub>-1.62</sub></td><td>(4.86)<sup>+1.27</sup><sub>-1.02</sub></td></tr><tr><td>O34363.21A</td><td>65.90 ± 6.59</td><td>(58.86)<sup>+15.18</sup><sub>-15.15</sub></td><td>(50.96)<sup>+8.26</sup><sub>-8.92</sub></td><td>(53.71)<sup>+9.86</sup><sub>-8.24</sub></td><td>(55.66)<sup>+7.90</sup><sub>-11.46</sub></td></tr><tr><td>S24068.60A</td><td>3.70 ± 1.00</td><td>(1.97)<sup>+0.61</sup><sub>-0.49</sub></td><td>(1.91)<sup>+0.71</sup><sub>-0.34</sub></td><td>(2.02)<sup>+0.80</sup><sub>-0.37</sub></td><td>(2.13)<sup>+0.45</sup><sub>-0.43</sub></td></tr><tr><td>S36312.06A</td><td>6.70 ± 2.00</td><td>(9.96)<sup>+1.78</sup><sub>-1.46</sub></td><td>(9.96)<sup>+1.16</sup><sub>-1.06</sub></td><td>(9.83)<sup>+0.89</sup><sub>-1.42</sub></td><td>(9.56)<sup>+1.16</sup><sub>-1.33</sub></td></tr><tr><td>Ar44740.12A</td><td>4.50 ± 3.00</td><td>(4.66)<sup>+2.70</sup><sub>-2.37</sub></td><td>(3.00)<sup>+1.16</sup><sub>-0.91</sub></td><td>(3.49)<sup>+1.45</sup><sub>-1.01</sub></td><td>(3.74)<sup>+1.36</sup><sub>-1.12</sub></td></tr><tr><td>Fe34987.33A</td><td>7.40 ± 3.00</td><td>(0.26)<sup>+0.09</sup><sub>-0.07</sub></td><td>(0.25)<sup>+0.09</sup><sub>-0.06</sub></td><td>(0.25)<sup>+0.06</sup><sub>-0.06</sub></td><td>(0.23)<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>Fe54227.19A</td><td>1.80 ± 0.70</td><td>(5.41)<sup>+4.26</sup><sub>-3.12</sub></td><td>(4.04)<sup>+1.69</sup><sub>-1.70</sub></td><td>(4.47)<sup>+1.85</sup><sub>-1.46</sub></td><td>(4.19)<sup>+1.02</sup><sub>-1.44</sub></td></tr></tbody></table>

Notes. Fluxes normalized to $H\beta = 1000$ need to be scaled by $4.56 \times 10^{36}$ to convert to luminosity in erg s$^{-1}$. Values in parentheses correspond to predicted observations. Model $M_{Iopt}$ considers only optical ionized gas observations, $M_{Iopt+IR}$ considers infrared lines as additional observations, $M_{NA}$ considers all lines tracing the ionized and neutral atomic gas (CO is only used as an upper limit), and $M_{NM}$ considers the neutral molecular gas.
<figure><img src="image_20.png" /><figcaption>**Fig. 14.** Origin of each line emission in the MULTIGRIS model for the ionized gas model *M*<sub>L opt</sub> (**left**), the neutral atomic gas model *M*<sub>NA</sub> (**middle**), and the neutral molecular gas model *M*<sub>NM</sub> (**right**).</figcaption></figure>
<table><caption>**Table 3. ISM model parameters.**</caption><thead><tr><th>Parameter</th><th>Description</th><th>Values/range</th></tr></thead><tbody><tr><td colspan="3">Primary parameters with distinct values for each component (sector)</td></tr><tr><td>w</td><td>Component (sector) mixing weight</td><td>[0–1]</td></tr><tr><td>log U<sub>0</sub></td><td>Ionization parameter <sup>(a)</sup></td><td>[-4, -3, -2, -1, 0]</td></tr><tr><td>log n<sub>0</sub> (cm<sup>-3</sup>)</td><td>Gas density <sup>(a)</sup></td><td>[0, 1, 2, 3, 4]</td></tr><tr><td>cut</td><td>Depth into the cloud <sup>(b)</sup></td><td>[0 – -4] step 0.25</td></tr><tr><td colspan="3">Primary parameters with a single value for all components (sectors)</td></tr><tr><td>log Z</td><td>Metallicity</td><td>[-2.19, -1.889, -1.19, -0.889, -0.667, -0.491, -0.19, 0.111]</td></tr><tr><td>log Z<sub>dust</sub></td><td>Dust-to-gas mass ratio <sup>(c)</sup></td><td>[-0.68, 0, +0.70]</td></tr><tr><td>log L<sub>X</sub> (L<sub>☉</sub>)</td><td>X-ray source luminosity</td><td>[-∞, -3, -2, -1]</td></tr><tr><td>log T<sub>X</sub> (K)</td><td>X-ray source temperature</td><td>[0, 5, 6, 7]</td></tr><tr><td>log L* (L<sub>☉</sub>)</td><td>Stellar cluster luminosity</td><td>[7, 9]</td></tr><tr><td>Age (Myr)</td><td>Burst age</td><td>[1, 2, 3, 4, 5, 6, 8, 10]</td></tr><tr><td colspan="3">Secondary (predicted) parameters in this study</td></tr><tr><td>log Q(H) (s<sup>-1</sup>)</td><td>Number of ionizing photons</td><td></td></tr><tr><td>f<sub>esc</sub>(1 – ∞) (%)</td><td>Escape fraction of photons with energy range [1 – ∞] Ryd</td><td></td></tr><tr><td>f<sub>esc</sub>(1–1.8) (%)</td><td>Escape fraction of photons with energy range [1–1.8] Ryd</td><td></td></tr><tr><td>f<sub>esc</sub>(1.8–4) (%)</td><td>Escape fraction of photons with energy range [1.8–4] Ryd</td><td></td></tr><tr><td>f<sub>esc</sub>(4–20) (%)</td><td>Escape fraction of photons with energy range [4–20] Ryd</td><td></td></tr><tr><td>f<sub>esc</sub>(20 – ∞) (%)</td><td>Escape fraction of photons with energy range [20 – ∞] Ryd</td><td></td></tr><tr><td>M<sub>H+</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sup>+</sup></td><td></td></tr><tr><td>M<sub>H0</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sup>0</sup></td><td></td></tr><tr><td>M<sub>H2</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sub>2</sub></td><td></td></tr><tr><td>M<sub>H2,C+</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sub>2</sub> traced by C<sup>+</sup></td><td></td></tr><tr><td>M<sub>H2,C0</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sub>2</sub> traced by C<sup>0</sup></td><td></td></tr><tr><td>M<sub>H2,CO</sub> (log M<sub>☉</sub>)</td><td>Mass of H<sub>2</sub> traced by CO</td><td></td></tr><tr><td>M<sub>dust</sub> (log M<sub>☉</sub>)</td><td>Dust mass</td><td></td></tr><tr><td>f<sub>CII,H+</sub></td><td>Fraction of [C II] associated with H<sup>+</sup></td><td></td></tr><tr><td>f<sub>CII,H0</sub></td><td>Fraction of [C II] associated with H<sup>0</sup></td><td></td></tr><tr><td>f<sub>CII,H2</sub></td><td>Fraction of [C II] associated with H<sub>2</sub></td><td></td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Calculated at the illuminated face of the cloud. <sup>(b)</sup>See definition in Sect. 4.2. <sup>(c)</sup>The range [-1, 1] corresponds to the 95% envelop (in dex) around the median relationship in Galliano et al. (2021).

the fraction of photons absorbed within the H II regions, $Q_{abs}/Q$, was estimated a priori to be 0.37 ± 0.03 and at least 0.30 based on the observed vs. expected $H\alpha$ emission and assuming that all ionizing photons are eventually absorbed in the galaxy. This leads to an expected escape fraction $f_{esc} \approx 63\%$ and at most 70% for the H II regions and the best models in P08 agree with these values. The escaping photons are responsible for the diffuse $H\alpha$ halo emission in the galaxy.

While the dust-to-gas mass ratio has little impact on the predicted emission lines in the H II region, it is interesting that the model is still able to recover the expected value. The range of [-1, 1] for $Z_{dust}$ corresponds to the 95% envelope around the median relationship of $Z_{dust}$ vs. $Z_{gas}$ in Galliano et al. (2021). I Zw 18 lies only slightly above the median relationship in that study while the SFGX models are allowed to reach up to five times below or above (Table 3).

The X-ray luminosity we infer is difficult to compare to P08. In P08, a stellar radiation field including an ad-hoc EUV component was used and an X-ray blackbody was added to explain the neutral gas emission lines. The stellar model BPASS (Eldridge et al. 2017) does not have enough energy in the EUV range to explain the He II emission so there is a need to include the X-ray multicolor blackbody model from the SFGX database.

Relatively, He II is the main line emission in sector #1 (Fig. 14). The X-ray luminosity we derive, $39.40^{+1.14}_{-0.33}$ (log erg s$^{-1}$), agrees with the value measured by Thuan et al. (2004), ≈39.48, with Chandra in the 0.5–10 keV range. It is somewhat lower than the XMM-Newton value in Kaaret & Feng (2013) of 40.15 (log erg s$^{-1}$) but still within uncertainties.

The fact that we obtain globally similar results to P08 is remarkable considering that we use different pressure profiles and a general prescription for the abundance patterns vs. metallicity. However, the automatic approach obviously cannot pretend to reach the same level of fine-tuning as specific studies.

Not surprisingly, the predicted [C II] in model $M_{l, opt}$, $1.6^{+2.2}_{-1.2}$, is much lower than the observed value ($76^{+20}_{-20}$) because [C II] is expected to emit deeper while we restricted ourselves to models reaching the ionization front. This is an indication that the fraction of observed [C II] in the ionized gas is small, which is confirmed by models accounting for the neutral gas (Sect. 6.2).

## 6.2. New topological models

We now attempt to update the results in L17 who converted the P08 modeling approach to Cloudy. The use of the latter
<table><caption>**Table 4.** I Zw 18 model parameters.</caption><thead><tr><th>Parameter</th><th><i>M</i><sub>L opt</sub></th><th><i>M</i><sub>L opt+IR</sub></th><th><i>M</i><sub>NA</sub></th><th><i>M</i><sub>NM</sub></th><th>P08/L17</th></tr></thead><tbody><tr><td><i>w</i></td><td>0.32<sup>+0.49</sup><sub>-0.30</sub><br>0.24<sup>+0.48</sup><sub>-0.24</sub><br>0.28<sup>+0.50</sup><sub>-0.27</sub></td><td>0.52<sup>+0.36</sup><sub>-0.49</sub><br>0.20<sup>+0.53</sup><sub>-0.19</sub><br>0.16<sup>+0.68</sup><sub>-0.14</sub></td><td>0.51<sup>+0.40</sup><sub>-0.41</sub><br>0.39<sup>+0.39</sup><sub>-0.37</sub><br>0.07<sup>+0.11</sup><sub>-0.06</sub></td><td>0.64<sup>+0.24</sup><sub>-0.36</sub><br>0.16<sup>+0.28</sup><sub>-0.14</sub><br>0.10<sup>+0.26</sup><sub>-0.10</sub><br>0.07<sup>+0.11</sup><sub>-0.05</sub></td><td></td></tr><tr><td><i>U</i><sub>0</sub></td><td>-2.08<sup>+1.63</sup><sub>-1.51</sub><br>-3.03<sup>+1.16</sup><sub>-0.97</sub><br>-1.72<sup>+1.72</sup><sub>-1.44</sub></td><td>-3.06<sup>+0.54</sup><sub>-0.44</sub><br>-2.46<sup>+1.70</sup><sub>-1.22</sub><br>-2.89<sup>+1.49</sup><sub>-1.02</sub></td><td>-2.82<sup>+1.21</sup><sub>-0.68</sub><br>-2.98<sup>+1.35</sup><sub>-0.54</sub><br>-1.78<sup>+1.23</sup><sub>-1.55</sub></td><td>-3.00<sup>+0.44</sup><sub>-0.43</sub><br>-2.25<sup>+1.39</sup><sub>-1.13</sub><br>-2.82<sup>+1.10</sup><sub>-0.65</sub><br>-2.27<sup>+1.41</sup><sub>-1.21</sub></td><td></td></tr><tr><td><i>n</i><sub>0</sub> (log cm<sup>-3</sup>)</td><td>1.14<sup>+1.13</sup><sub>-1.09</sub><br>2.20<sup>+1.29</sup><sub>-1.15</sub><br>0.72<sup>+1.22</sup><sub>-0.71</sub></td><td>1.01<sup>+0.91</sup><sub>-0.91</sub><br>1.53<sup>+1.96</sup><sub>-1.35</sub><br>2.34<sup>+1.18</sup><sub>-1.70</sub></td><td>1.21<sup>+1.12</sup><sub>-0.69</sub><br>1.49<sup>+1.41</sup><sub>-0.92</sub><br>2.17<sup>+1.01</sup><sub>-0.64</sub></td><td>1.44<sup>+0.95</sup><sub>-0.82</sub><br>1.85<sup>+1.26</sup><sub>-1.24</sub><br>1.95<sup>+1.38</sup><sub>-1.13</sub><br>2.02<sup>+1.04</sup><sub>-0.75</sub></td><td></td></tr><tr><td>cut <sup>(a)</sup></td><td>0.42<sup>+0.45</sup><sub>-0.35</sub><br>0.59<sup>+0.52</sup><sub>-0.45</sub><br>0.93<sup>+0.20</sup><sub>-0.51</sub></td><td>0.65<sup>+0.33</sup><sub>-0.39</sub><br>0.25<sup>+0.44</sup><sub>-0.23</sub><br>0.76<sup>+0.35</sup><sub>-0.47</sub></td><td>0.36<sup>+0.38</sup><sub>-0.35</sub><br>0.72<sup>+0.38</sup><sub>-0.31</sub><br>1.32<sup>+0.51</sup><sub>-0.19</sub></td><td>0.59<sup>+0.28</sup><sub>-0.24</sub><br>0.52<sup>+0.58</sup><sub>-0.40</sub><br>1.01<sup>+0.66</sup><sub>-0.65</sub><br>1.56<sup>+0.83</sup><sub>-0.52</sub></td><td></td></tr><tr><td>Z (log Z<sub>☉</sub>)</td><td>-1.45<sup>+0.10</sup><sub>-0.08</sub></td><td>-1.46<sup>+0.10</sup><sub>-0.09</sub></td><td>-1.44<sup>+0.07</sup><sub>-0.07</sub></td><td>-1.41<sup>+0.07</sup><sub>-0.07</sub></td><td>-1.46 (P08; O/H)</td></tr><tr><td>Z<sub>dust</sub> <sup>(b)</sup></td><td>-0.17<sup>+0.66</sup><sub>-0.46</sub></td><td>-0.04<sup>+0.57</sup><sub>-0.50</sub></td><td>0.10<sup>+0.45</sup><sub>-0.41</sub></td><td>0.34<sup>+0.27</sup><sub>-0.52</sub></td><td></td></tr><tr><td><i>L</i><sub>X</sub> (log erg s<sup>-1</sup>)</td><td>39.40<sup>+1.14</sup><sub>-0.33</sub></td><td>39.99<sup>+0.59</sup><sub>-0.74</sub></td><td>39.87<sup>+0.35</sup><sub>-0.67</sub></td><td>39.96<sup>+0.36</sup><sub>-0.60</sub></td><td></td></tr><tr><td><i>T</i><sub>X</sub> (log K)</td><td>6.08<sup>+0.66</sup><sub>-0.58</sub></td><td>6.16<sup>+0.47</sup><sub>-0.64</sub></td><td>5.92<sup>+0.47</sup><sub>-0.37</sub></td><td>6.04<sup>+0.38</sup><sub>-0.50</sub></td><td></td></tr><tr><td><i>Q</i>(H) (log s<sup>-1</sup>)</td><td>52.24<sup>+0.42</sup><sub>-0.24</sub></td><td>52.28<sup>+0.39</sup><sub>-0.25</sub></td><td>52.21<sup>+0.21</sup><sub>-0.21</sub></td><td>52.13<sup>+0.26</sup><sub>-0.16</sub></td><td>52.1 (P08; 18.2 Mpc)</td></tr><tr><td><i>f</i><sub>esc</sub>(1 &rarr; &infin;) (%)</td><td>56.07<sup>+37.20</sup><sub>-34.01</sub></td><td>59.20<sup>+35.85</sup><sub>-19.19</sub></td><td>54.17<sup>+16.90</sup><sub>-21.74</sub></td><td>47.56<sup>+29.14</sup><sub>-14.80</sub></td><td>&approx;63, [60, 70] (P08)</td></tr><tr><td>- <i>f</i><sub>esc</sub>(1&ndash;1.8) (%)</td><td>48.12<sup>+35.72</sup><sub>-31.57</sub></td><td>51.89<sup>+30.48</sup><sub>-26.93</sub></td><td>47.02<sup>+20.89</sup><sub>-27.18</sub></td><td>41.73<sup>+30.10</sup><sub>-22.32</sub></td><td></td></tr><tr><td>- <i>f</i><sub>esc</sub>(1.8&ndash;4) (%)</td><td>68.53<sup>+31.07</sup><sub>-42.82</sub></td><td>78.05<sup>+20.58</sup><sub>-18.82</sub></td><td>71.55<sup>+23.42</sup><sub>-18.14</sub></td><td>71.26<sup>+21.88</sup><sub>-18.39</sub></td><td></td></tr><tr><td>- <i>f</i><sub>esc</sub>(4&ndash;20) (%)</td><td>20.38<sup>+35.63</sup><sub>-16.27</sub></td><td>26.31<sup>+33.29</sup><sub>-15.58</sub></td><td>24.14<sup>+33.70</sup><sub>-8.43</sub></td><td>22.07<sup>+14.30</sup><sub>-7.92</sub></td><td></td></tr><tr><td>- <i>f</i><sub>esc</sub>(20 &rarr; &infin;) (%)</td><td>79.99<sup>+19.09</sup><sub>-17.88</sub></td><td>87.37<sup>+11.78</sup><sub>-8.55</sub></td><td>81.34<sup>+7.60</sup><sub>-11.75</sub></td><td>78.07<sup>+9.72</sup><sub>-9.59</sub></td><td></td></tr><tr><td><i>M</i><sub>H+</sub> (log <i>M</i><sub>☉</sub>)</td><td>6.46<sup>+0.91</sup><sub>-0.92</sub></td><td>6.47<sup>+1.05</sup><sub>-0.94</sub></td><td>6.42<sup>+0.12</sup><sub>-0.91</sub></td><td>6.29<sup>+0.24</sup><sub>-0.79</sub></td><td>6.5 (P08; 18.2 Mpc)</td></tr><tr><td><i>M</i><sub>H0</sub> (log <i>M</i><sub>☉</sub>)</td><td>4.25<sup>+1.32</sup><sub>-1.41</sub></td><td>3.82<sup>+1.10</sup><sub>-1.19</sub></td><td>7.28<sup>+0.56</sup><sub>-0.78</sub></td><td>7.23<sup>+0.60</sup><sub>-0.95</sub></td><td>&approx;[6.9&ndash;7.3] (L17)</td></tr><tr><td><i>M</i><sub>H2</sub> (log <i>M</i><sub>☉</sub>)</td><td>-5.67<sup>+2.78</sup><sub>-4.38</sub></td><td>-7.22<sup>+3.25</sup><sub>-3.03</sub></td><td>3.82<sup>+2.51</sup><sub>-2.50</sub></td><td>5.39<sup>+1.67</sup><sub>-4.17</sub></td><td></td></tr><tr><td>- <i>M</i><sub>H2,C+</sub> (log <i>M</i><sub>☉</sub>)</td><td>-5.81<sup>+2.79</sup><sub>-5.25</sub></td><td>-7.40<sup>+3.10</sup><sub>-4.04</sub></td><td>3.82<sup>+2.41</sup><sub>-2.50</sub></td><td>5.38<sup>+1.50</sup><sub>-4.16</sub></td><td></td></tr><tr><td>- <i>M</i><sub>H2,C</sub> (log <i>M</i><sub>☉</sub>)</td><td>-9.23<sup>+3.15</sup><sub>-5.07</sub></td><td>-10.41<sup>+4.04</sup><sub>-3.68</sub></td><td>0.38<sup>+3.24</sup><sub>-2.98</sub></td><td>2.28<sup>+3.05</sup><sub>-4.90</sub></td><td></td></tr><tr><td>- <i>M</i><sub>H2,CO</sub> (log <i>M</i><sub>☉</sub>)</td><td>-19.04<sup>+3.31</sup><sub>-11.29</sub></td><td>-22.34<sup>+6.66</sup><sub>-7.35</sub></td><td>-3.32<sup>+5.17</sup><sub>-4.95</sub></td><td>-0.19<sup>+4.37</sup><sub>-7.74</sub></td><td></td></tr><tr><td><i>M</i><sub>dust</sub> (log <i>M</i><sub>☉</sub>)</td><td>1.16<sup>+1.09</sup><sub>-1.08</sub></td><td>1.16<sup>+1.20</sup><sub>-1.09</sub></td><td>2.23<sup>+0.52</sup><sub>-0.63</sub></td><td>2.45<sup>+0.50</sup><sub>-0.91</sub></td><td>2.2 (L17), 2.86 &pm; 0.24 <sup>(c)</sup></td></tr><tr><td><i>f</i><sub>CII,H+</sub></td><td>0.95<sup>+0.05</sup><sub>-0.03</sub></td><td>1.00<sup>+0.00</sup><sub>-0.05</sub></td><td>0.06<sup>+0.16</sup><sub>-0.05</sub></td><td>0.04<sup>+0.15</sup><sub>-0.04</sub></td><td></td></tr><tr><td><i>f</i><sub>CII,H</sub></td><td>0.05<sup>+0.03</sup><sub>-0.05</sub></td><td>0.00<sup>+0.05</sup><sub>-0.00</sub></td><td>0.94<sup>+0.06</sup><sub>-0.16</sub></td><td>0.93<sup>+0.05</sup><sub>-0.20</sub></td><td></td></tr><tr><td><i>f</i><sub>CII,H2</sub></td><td>0.00<sup>+0.00</sup><sub>-0.00</sub></td><td>0.00<sup>+0.00</sup><sub>-0.00</sub></td><td>0.00<sup>+0.03</sup><sub>-0.00</sub></td><td>0.01<sup>+0.13</sup><sub>-0.01</sub></td><td></td></tr></tbody></table>

**Notes.** Models $M_{L, opt}$ and $M_{L, opt+IR}$ consider ionized gas observations assuming P08's modeling strategy without and with infrared lines respectively. $M_{NA}$ adds neutral atomic gas observations ([C II], [O I], [Si II], and [Fe II]). $M_{NM}$ adds the CO marginal detection of Zhou et al. (2021) and upper limits on near-IR and mid-IR H$_2$. The escape fractions are given for several energy ranges in Rydberg. The reported values correspond to the median of the posterior distribution together with the confidence interval. Combined quantities (e.g., $M_{H2}$) are calculated for each draw and the final median value can therefore deviate from the combination of median values. (a) 1: ionization front (H$^+$–H$^0$ transition), 2: PDR front (H$^0$–H$_2$ transition), 3: C$^0$–CO transition. (b) Dust-to-gas mass ratio. See Table 3 for more details. (c) Value from Rémy-Ruyer et al. (2015) assuming graphite composition.

was motivated by the need to account for dust in order to compare the relative importance of the photoelectric effect heating and X-ray photoionization in the neutral gas. To this effect L17 introduced several neutral atomic and molecular lines in the model but parameters were adjusted manually and there was no convergence test. Apart from an exploratory model with four

sectors including a deep sector with constant pressure to investigate the molecular gas, the sector configuration in L17 was precisely the same as in P08.

For all the models described below, in addition to new IR constraints probing the ionized or neutral gas, we also relax some of the constraints on the topology, notably by letting the
inner radius free for each sector. Furthermore, while we fixed the number of components in model $M_{\text{Iopt}}$ to reproduce the results of P08 in Sect. 6.1, we now test different number of components. All model predictions are provided for convenience in the appendix (Figs. A.1–A.3; Table A.1).

### 6.2.1. Ionized gas: Model $M_{\text{Iopt+IR}}$

For model $M_{\text{Iopt+IR}}$ we consider the ionized gas again but this time considering IR lines in addition to optical ones. Results are shown in Table 2 and Fig. 13. Some newly added observations arise from the same physical conditions as the optical lines and thus provide further constraints to the model, but some other observations provide new constraints, such as [O III] 88 µm for the rather diffuse highly-ionized gas Fig. 15 confirms that the optical and IR [O III] lines do not trace the same physical conditions.

The model comparison metrics indicate that the most likely number of sectors is one but this does not account for the fact that the prior odds for a model with a single sector is expected to be low (Sect. 3.7). It turns out, however, that the predictions are in good agreement with models with more than one sector (Fig. A.1). We also tested various combinations with the ionization parameter $U$ and gas density $n$ as power-law distributions and the cut parameter as either a power-law or normal distribution (model $M_{\text{I, power-law}}$). Results are shown in Fig. A.1 and Table A.1. Globally, even though those latter models have a somewhat lower probability, the predictions for all secondary parameters are in agreement with the values inferred with a small, discrete, number of components.

Using the model $M_{\text{Iopt+IR}}$ as an illustration (i.e., with three sectors), the X-ray luminosity is now slightly larger compared to $M_{\text{I, opt}}$. The predicted value, $L_X = 39.99^{+0.59}_{-0.74}$ (log erg s⁻¹), is now in better agreement with the XMM-Newton value in Kaaret & Feng (2013) of 40.15 (log erg s⁻¹).

Other notable differences with model $M_{\text{I, opt}}$ include the mass of neutral atomic hydrogen $M(\text{H}^0)$ and the dust mass $M_{\text{dust}}$ which become somewhat lower. However, the values have large uncertainties reflecting the fact that model $M_{\text{I, opt+IR}}$ holds little information on the neutral gas.

### 6.2.2. Neutral atomic gas: Model $M_{\text{NA}}$

For the neutral atomic gas model $M_{\text{NA}}$, we include [C II], [O I], [Si II], and [Fe II] as additional constraints. We also include the total IR emission (which, according to L17's models originates from dust emission in the H II region, a moderate contribution from dust in the H I region, and a small but significant contribution from free-free emission). We allow the cloud depth to reach slightly beyond the H/H₂ transition in order to fully consider the neutral atomic layer.

The model comparison metrics (Sect. 3.7) indicate that model with a single sector is highly unlikely (Fig. A.2). Results with more than one sector are in good agreement. As for the ionized gas model ($M_{\text{I, opt+IR}}$; Sect. 6.2.1), we also tested various combinations with $U$ and $n$ as power-law distributions and the cut parameter as a broken power-law (model $M_{\text{NA, power-law}}$). The use of a broken power-law instead of a single power-law is motivated by the fact that the extrapolation of the distribution of depth values for matter-bounded regions does not have to hold for the radiation-bounded regions. Not surprisingly, we find a pivot cut around one, that is at the ionization front, with significantly different slopes below and above (Table A.1). The corresponding model $M_{\text{NA, power-law}}$ has a relatively lower

probability as compared to the models with a small, discrete, number of components, indicating that more complex distributions may be needed. Nevertheless, the predictions are globally in good agreement with model $M_{\text{NA}}$ (Fig. A.1).

Parameter values for model $M_{\text{NA}}$ are shown in Tables 2 and 4. Most global parameters (e.g., metallicity and mass) are unchanged when compared to model $M_{\text{I, opt+IR}}$ with only ionized gas observations. However, the mixing weights and depth of the sectors are significantly different compared to the ionized gas model, which explains the change in the escape fraction values, which strongly depends on the topology (Sect. 5.3). The escape fraction somewhat decreased compared to $M_{\text{I, opt+IR}}$, and, while the difference is relatively small, this could imply that some of the neutral gas emission arises on somewhat larger spatial scales than the H II region.

The mass of atomic hydrogen is now much larger compared to $M_{\text{I, opt+IR}}$ and in agreement with expectations (Table 4). We find a mass of molecular gas of about $\approx 10^{3.8} M_{\odot}$. L17 had evaluated the molecular gas mass assuming the same topology as $M_{\text{I, opt}}$, that is with the same covering factors as in P08. This effectively corresponded to a rather diffuse molecular medium as opposed to clumps that were considered in a specific, separate, four-sector model. The H₂ mass distributed in the large-scale radiation-bounded sector in L17 was extremely low, on the order of $1 M_{\odot}$. We find a much larger value in our model $M_{\text{NA}}$ because the topology is different (sector mixing weight and different inner radii) and the H₂ mass is furthermore controlled by the somewhat arbitrary maximum depth criterion (which we fix to slightly beyond the H/H₂ transition for $M_{\text{NA}}$). This rather small H₂ mass of $\approx 10^{3.8} M_{\odot}$ exists in the C⁺ zone and is therefore CO-dark with a CO-dark H₂ fraction of $\approx 100\%$ (but the overall [C II] emission is still dominated by the neutral atomic gas). The dust mass is also close to the observed value and the dust-to-gas mass ratio is now slightly larger compared to $M_{\text{I, opt+IR}}$ and in better agreement with Galliano et al. (2021), which is not surprising since the dust mass mostly comes from the neutral phase. However, by comparing models $M_{\text{I, opt}}$ (or $M_{\text{I, opt+IR}}$) and $M_{\text{NA}}$ we find that $\approx 2/3$ of the TIR emission comes from the ionized phase, which is in agreement with L17 who found that TIR can be decomposed as $\approx 20\%$ free-free emission, $\approx 15\%$ dust in the H I region, and $\approx 65\%$ dust in the H II region.

The [C II] emission is much larger compared to $M_{\text{I, opt+IR}}$ (Table 2) confirming that most of the predicted [C II] emission arises in the neutral atomic gas. We find a fraction of [C II] arising in the neutral atomic gas, $f_{\text{CII,H}}$ of about 94%. The deep narrow sector #3 is responsible for most of the [C II], [O I], [Si II], [Fe II], and H₂ line emission and also dominates the [Ar II] and [Ne II] emission (Fig. 14; Table 4). We note that the relative fluxes predicted for [C II], [O I], and [Si II] differ somewhat compared to observations, which is likely due to the density profile assumed in the models. The model $M_{\text{NA}}$ predicts an essentially null CO emission.

### 6.2.3. Neutral molecular gas: model $M_{\text{NM}}$

We now add molecular gas observations as constraints for model $M_{\text{NM}}$. For CO, we use the recently claimed CO(2–1) detection of Zhou et al. (2021), corresponding to a CO(1–0) luminosity of $0.18 L_{\odot}$ ($\approx 3.5\sigma$ detection). Even though the signal is quite uncertain, the new observation is much deeper than previous studies and provides at the very least a useful upper limit. We also include several upper limits for near-IR and mid-IR H₂ lines (see L17).
<figure><img src="image_21.png" /><figcaption>**Fig. 15.** Map of the correlation coefficients between predicted values for the ionized gas model $M_{\text{Iopt}}$ of I Zw 18 (see Fig. 12 for the plot description). We can see for instance that [O III] 88 µm does not correlate with the same parameters compared to other lines.</figcaption></figure>

The model comparison metrics (Sect. 3.7) indicate that the model with a single sector is highly unlikely (Fig. A.3). As for model $M_{\text{NA}}$, we also tested various combinations with *U* and *n* as power-law distributions and the cut parameter as a broken power-law. The corresponding model $M_{\text{NM,power-law}}$ reproduces better the observed CO but has a much lower probability as compared to the models with a small, discrete, number of components, indicating that more complex distributions of the cut parameter throughout the ionized, neutral atomic, and neutral molecular gas may be needed (Fig. A.1). In the following we discuss the model $M_{\text{NM}}$ composed of 4 sectors.

In Table 4, we can see that the narrowest sector now reaches somewhat deeper, but with a most likely value not quite reaching the H⁰-H₂ transition. Otherwise the other primary parameters are mostly unchanged compared to $M_{\text{NA}}$.

In Fig. 16, we see that the most likely values for I(CO) are much lower than the observed value even within uncertainties. CO being the only observation constraining such depths and the detection level being so low, the model is mostly constrained by other lines, in particular those arising mostly in the neutral atomic gas as well as TIR. In other words, the CO detection is in fact unlikely within the model topological assumptions in model $M_{\text{NM}}$. There is a solution agreeing with the observed CO intensity, corresponding to an H₂ mass ≈10⁷ $M_{\odot}$, but with a relatively low probability of ≤20%. This large value coincides with the upper limit on the H₂ mass of L17 based on the observed CO(1–0) upper limit of <8 × 10⁻³ (normalized to Hβ = 1000) from Leroy et al. (2007).

When the observed CO flux is well reproduced, the H⁰ mass becomes somewhat too large, being significantly above 10⁷ $M_{\odot}$

(Fig. 16) and the total IR emission starts deviating from the observed value (Fig. 17). Adding CO tends to make a deeper narrower sector as the CO emission is better and better reproduced (though these models are unlikely). The mixing weight of the deep sector ought to be lower otherwise it would produce too much [O I] and H₂ in particular. Even with four sectors (model $M_{\text{NM}}$, the CO emission always come from the same sector as [C II] and [O I]. In other words, if we need to account for CO then we need a deep narrow sector and this sector produces a significant/dominant emission of [C II] and [O I]. This implies that the [C II] and [O I] emission should be probing a mostly neutral atomic medium associated to the CO clumps if they do exist, since the maximum fraction of [C II] associated with H₂ is ≈15% (Fig. 16).

A more stringent observational constraint would be useful either to confirm that the CO emission is in fact much lower or to constrain a better topological model. Nevertheless, the predicted H₂ mass is about 10⁵.⁴ $M_{\odot}$ and anticorrelates with the covering factor of the dense sector (top panel of Fig. 16). As for model $M_{\text{NA}}$, the small H₂ mass is better traced by C⁺ than C⁰ or CO (Table 4), and the CO-dark H₂ gas fraction is therefore ≈100%. However, most of the [C II] emission is not associated with H₂ but with H⁰, in agreement with L17.

The CO luminosity we predict is ≈3500 times lower than Zhou et al. (2021), implying $\alpha_{\text{CO}} = M(\text{H}_2)/L(\text{CO}) = 9 \times 10^4 M_{\odot} (\text{K km s}^{-1})^{-1}$. For the metallicity of I Zw 18, together with the reference solar value for $\alpha_{\text{CO}} = 4.3 M_{\odot} \text{ pc}^{-2} (\text{K km s}^{-1})^{-1}$ (e.g., Bolatto et al. 2013), this translates into a slope of $y = -3.1$ with $\alpha_{\text{CO}} \propto (Z/Z_{\odot})^y$. This coefficient is on the high end of metallicity-dependent
<figure><img src="image_22.png" /><figcaption>**Fig. 16.** Properties of the molecular sector in model $M_{\text{NM}}$. From top to bottom, we show as a function of the total predicted H$_2$ mass, the mixing weight of the molecular sector, the mass of H$^0$, of dust, the CO(1–0) intensity (scaled to $4.56 \times 10^{36}$ erg s$^{-1}$; with the horizontal lines showing the observed range in Zhou et al. (2021), and the fraction of [C II] originating from the molecular sector.</figcaption></figure>

$\alpha_{\text{CO}}$ studies but quite close to the value derived by Madden et al. (2020) of –3.39 and the extrapolation of Schruba et al. (2012; Fig. 18). If we consider instead the less likely solution at $M(H_2) \sim 10^7 M_\odot$ corresponding to the CO observation in Zhou et al. (2021), we obtain a $\alpha_{\text{CO}} \approx 2200 M_\odot \text{ pc}^{-2} (\text{K km s}^{-1})^{-1}$, implying $y \approx -2$. The calculation and discussion of $\alpha_{\text{CO}}$ for the full DGS sample will be discussed in Ramambason et al. (in prep.).

### 6.2.4. Escape fraction of ionizing photons

In all models, more than half of the global escape fraction is due to sector #1 which is always matter-bounded, and with the largest mixing weight (Fig. 19). The observed fluxes correspond to the NW H II region of 1Zw 18 and the global escape fraction therefore also corresponds to that region, while the galaxy escape fraction is probably significantly lower (P08). We refer to

<figure><img src="image_23.png" /><figcaption>**Fig. 17.** Variation of the observed fluxes as a function of I(CO). From top to bottom we show the IR flux, the [O I] flux, and the H$_2$ S(1) line flux. Fluxes are scaled to $4.56 \times 10^{36}$ erg s$^{-1}$. Horizontal and vertical lines show the observed uncertainties or the upper limit (for H$_2$).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 18.** Variation of $\alpha_{\text{CO}}$ as a function of metallicity compared to Schruba et al. (2012); Sandstrom et al. (2013); Amorín et al. (2016); Madden et al. (2020). The red data points and upper limits are from Madden et al. (2020). We also show on top the most likely solution (empty black circle with the gray circle indicating the standard deviation) with $M(H_2) = 10^{5.4} M_\odot$ and CO 3500 times fainter than in Zhou et al. (2021) and on bottom the secondary, less likely, solution with $10^7 M_\odot$ and CO in agreement with Zhou et al. (2021).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 19.** Sector configurations for models $M_{\text{I,opt+IR}}$ (a), $M_{\text{NA}}$ (b), and $M_{\text{NM}}$ (c) for 1Zw 18. See Fig. 6 for the plot description. The arrows indicate the escape fraction $(1 - \infty \text{Ryd})$ for each sector.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 20.** Absorption of the ionizing radiation throughout the gas slab. Top: shape of the radiation field spectrum for each zone calculated in Cloudy for the radiation-bounded sector of model $M_{\text{NM}}$, from the illuminated edge of the cloud (blue) to the last zone (green). The vertical lines show the 1, 1.8, 4, 8, 20, and 2000 Ryd energy values. The spectrum is heavily and rapidly absorbed right above 4 Ryd due to $\text{He}^{2+}$ absorption. Bottom: evolution of the integrated flux across the cloud depth (normalized to the value at the illuminated edge) for different energy ranges in Rydberg. The vertical lines show the ionization front, PDR, and C/CO transitions from left to right. The radiation field >20 Ryd is able to penetrate in the neutral gas (top panel) but this contributes little to the total escape fraction (bottom panel). The kink for the 1.8–4 Ryd range is due to secondary ionization following He recombination.</figcaption></figure>

which traces the escape of soft to hard X-ray photons is ≈87% in the model. Even though soft X-ray photons deposit their energy in the neutral gas (L17), $f_{\text{esc}}(20 - \infty)$ remains large even in model $M_{\text{NA}}$ (≈80%), which is due to the fact that the mixing weight of the radiation-bounded sector is relatively low. Most X-ray photons therefore escape the H II region, but the global X-ray

Ramambason et al. (2022) for a more general discussion of the escape fractions in the Dwarf Galaxy Survey.

Our models include the presence of an X-ray source associated with the H II region (L17). For the first time we are now able to calculate the escape fraction for different photon energy ranges and evaluate the escape fraction of He-ionizing photons, which have received much less attention than H-ionizing photons (Benson et al. 2013). The escape fractions of photons in our models is saved in the energy ranges 1 – ∞, 1–1.8, 1.8–4 (He⁰-ionizing photons), 4–20 (He⁺-ionizing photons), and 20 – ∞ Ryd (Table 2). In model $M_{\text{I,opt+IR}}$ (Sect. 6.2.1), the escape fraction $f_{\text{esc}}(1 - \infty) \approx 59\%$ lies between $f_{\text{esc}}(1-1.8) \approx 52\%$ and $f_{\text{esc}}(1.8-4) \approx 78\%$, that is in the range where most of the input radiation field energy lies. The escape fraction $f_{\text{esc}}(4-20)$, on the other hand, is remarkably lower (≈26%), which is due to the absorption by $\text{He}^{2+}$ close to the illuminated edge of the cloud (Fig. 20). Nevertheless, the $f_{\text{esc}}(4-20)$ values remains significant due to the fact that the medium gradually becomes more transparent to photons significantly above 4 Ryd. As a matter of fact, $f_{\text{esc}}(20 - \infty)$,
escape fraction reaching the IGM ultimately depends on the dis-
tribution of diffuse H I gas. Since low-metallicity star-forming
dwarf galaxies are known to harbor bright X-ray sources (e.g.,
Basu-Zych et al. 2016; Ponnada et al. 2019; Lehmer et al. 2020),
this makes such sources interesting suspects for the IGM heating
at z > 10 (Mirabel et al. 2011).

7. Conclusions

In this study, we present a new statistical framework to infer
probability density functions (PDFs) of parameters constrained
by observations with uncertainties. While the method and asso-
ciated code, MULTIGRIS, are agnostic to the input grid used,
the method was developed in particular for ISM studies with the
objective of inferring a topology representing either a discrete
number of ISM components around one or several stellar clusters
or a continuous distribution of clouds whose physical conditions
are parameterized by power- or normal-laws.

The inference is performed using a sequential Monte Carlo
method that is well adapted to multipeaked posterior distribu-
tions. Primary parameters are sampled using either a nearest
neighbor or linear interpolation on the grid values. Built-in pri-
ors are available to use predefined configurations flexible enough
to various combinations of components.

We illustrate the method with a grid of models calculated
with Cloudy. Primary parameters describe the stellar clusters and
potential X-ray sources as well as the surrounding ISM distribu-
tion. The grid, described in Ramambason et al. (2022) includes
a free parameter defining the cloud depth. While the inference is
performed on these primary parameters, the algorithm is able to
calculate PDFs of secondary parameters, that is any parameter
not considered as an observation and which Cloudy can com-
pute. A large number of random variables and observations can
be accommodated. Observations may be gathered in observation
sets with associated systematic uncertainties, and upper/lower
limits are fully considered. For ISM studies, the extinction may
be included as one additional random variable.

We first apply this method to the well-known low-metallicity
star-forming dwarf galaxy I Zw 18, whose ISM has been modeled
extensively using optical and IR lines in particular. Our objective
is to reproduce previous results with a general grid that is not
tailored to any specific object and to provide PDFs for the most
important parameters. We first consider only ionized gas obser-
vations, then add neutral atomic and finally neutral molecular
gas observations in order to test the robustness of the diagnos-
tics adding/removing some observations while examine physical
parameters related to specific ISM phases.

We are able to reproduce the metallicity, number of ionizing
photons, as well as the mass of H+ and H0, the dust mass and
the dust-to-gas mass ratio. We confirm that an X-ray source with
luminosity ≈10⁴⁰ erg s⁻¹ is necessary to explain the ionized and
neutral gas emission, with properties similar to those observed in
specific X-ray studies. This implies interesting prospects for con-
straining the presence, properties, and nature of compact objects
from their signatures in the ISM. Such a method is only possible
if the topology is well constrained by many observations, with at
least some of them themselves influenced by the X-ray source.

We examine the H₂ mass, using in particular the mid-IR H₂
upper limits as well as the recently claimed detection of
CO(1–0) by Zhou et al. (2021). With the topology used, we
find that the CO(1–0) emission is likely much lower than the
observed value. It is possible, however, that CO exists in small
dense clumps that we are not able to model. Nevertheless,
we derive a most likely mass of H₂ of ≈10⁵.⁴ M☉, leading to

$\alpha_{CO} = 9 \times 10^4 M_{\odot} \text{ pc}^{-2} (\text{K km s}^{-1})^{-1}$ and to a metallicity coefficient of $y = -3.1$ with $\alpha_{CO} \propto (Z/Z_{\odot})^y$, in excellent agreement with Madden et al. (2020) and Schruba et al. (2012). The calculation and discussion of $\alpha_{CO}$ for the full DGS sample will be discussed in Ramambason et al. (in prep.).

The H₂ mass is better traced by C⁺ than C⁰ or CO and the
CO-dark H₂ mass is therefore ≈100%. However, the fraction of
[C II] arising from the molecular gas is small and most of the
[C II] emission is associated with H⁰. Assuming CO clumps
exist, the predicted [C II] and [O I] emission arise from the same
sector as CO, implying that [C II] and [O I] lines probe a mostly
neutral atomic medium directly associated with the hypothetical
CO clumps.

Finally we investigate the escape fraction of ionizing photons
and find values around 50–65% in agreement with expectations.
We decompose the photons by energy range and find that most
(≥80%) soft X-ray photons are able to escape. This has potential
implications for the ionization and heating of the IGM as such
photons can escape the diffuse H I gas halo.

*Acknowledgements.* This work has greatly benefited from the environment pro-
vided by the FACE Foundation Thomas Jefferson Fund (award TJF2I_053). We
would like to thank Frédéric Galliano for insightful and inspiring discussions
about Bayesian statistics, and the anonymous referee for constructive feedback
and encouraging comments.

References

Aloisi, A., Clementini, G., Tosi, M., et al. 2007, ApJ, 667, L151

Amorín, R., Muñoz-Tuñón, C., Aguerri, J. A. L., & Planesas, P. 2016, A&A, 588, A23

Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, ArXiv e-prints
[arXiv:1012.3754]

Basu-Zych, A. R., Lehmer, B., Fragos, T., et al. 2016, ApJ, 818, 140

Benson, A., Venkatesan, A., & Shull, J. M. 2013, ApJ, 770, 76

Bisbas, T. G., Bell, T. A., Viti, S., Yates, J., & Barlow, M. J. 2012, MNRAS, 427, 2100

Bisbas, T. G., Haworth, T. J., Barlow, M. J., et al. 2015, MNRAS, 454, 2828

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Ching, J., & Chen, Y.-C. 2007, J. Eng. Mech., 133, 816

Cormier, D. 2012, Ph.D. Thesis, Université Paris-Diderot, France

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

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

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

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

Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038

Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137

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

Fischer, C., Beckmann, S., Bryant, A., et al. 2018, J. Astron. Instrum., 7, 1840003

Galliano, F. 2022, Habilitation Thesis, Université Paris-Saclay, France

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

Gordon, K. D., Fouesneau, M., Arab, H., et al. 2016, ApJ, 826, 104

Harikane, Y., Ouchi, M., Shibuya, T., et al. 2018, ApJ, 859, 84

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

Hoffman, M. D., & Gelman, A. 2011, ArXiv e-prints [arXiv:1111.4246]

Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

Izotov, Y. I., Chaffee, F. H., Foltz, C. B., et al. 1999, ApJ, 527, 757

Jin, Y., Kewley, L. J., & Sutherland, R. 2022, ApJ, 927, 37

Kaaret, P., & Feng, H. 2013, ApJ, 770, 20

Kang, D. E., Pellegrini, E. W., Ardizzone, L., et al. 2022, MNRAS, 512, 617

Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511

Lambert-Huyghe, A., Madden, S. C., Lebouteiller, V., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202243353

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

Lee, M. M., Nagao, T., De Breuck, C., et al. 2019, ApJ, 883, L29

Lee, M. M., Nagao, T., De Breuck, C., et al. 2021, ApJ, 913, 41

Lehmer, B. D., Ferrell, A. P., Doore, K., et al. 2020, ApJS, 248, 31
Leroy, A., Bolatto, A., Stanimirovic, S., et al. 2007, ApJ, 658, 1027

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

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

Minson, S. E., Simons, M., & Beck, J. L. 2013, Geophys. J. Int., 194, 1701

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

Morisset, C. 2011, Cloudy_3D: Quick Pseudo-3D Photoionization Code

Morisset, C. 2013, pyCloudy: Tools to manage astronomical Cloudy photoionization code

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

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, A1

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, A2

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

Ponnada, S., Brorby, M., & Kaaret, P. 2019, MNRAS, 2537

Ramambason, L., Lebouteiller, V., Bik, A., et al. 2022, A&A, in press https://doi.org/10.1051/0004-6361/202243866

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

Richardson, C. T., Allen, J. T., Baldwin, J. A., Hewett, P. C., & Ferland, G. J. 2014, MNRAS, 437, 2376

Richardson, C. T., Allen, J. T., Baldwin, J. A., et al. 2016, MNRAS, 458, 988

Robert, C. P. 2015, ArXiv e-prints [arXiv:1504.01896]

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55

Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5

Schruba, A., Leroy, A. K., Walter, F., et al. 2012, AJ, 143, 138

Spinoglio, L., Alonso-Herrero, A., Armus, L., et al. 2017, PASA, 34, e057

ter Braak, C. J. F., & Vrugt, J. A. 2008, Stat. Comput., 18, 435

Thomas, A. D., Kewley, L. J., Dopita, M. A., et al. 2019, ApJ, 874, 100

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

Vale Asari, N., Stasińska, G., Morisset, C., & Cid Fernandes, R. 2016, MNRAS, 460, 1739

Vehtari, A., Gelman, A., & Gabry, J. 2017, Stat. Comput., 27, 1413

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. 2021, Bayesian Anal., 16, 667

Wagg, J., Aravena, M., Brisbin, D., et al. 2020, MNRAS, 499, 1788

Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1

Watanabe, S. 2010, J. Mach. Learn. Res., 11, 3571

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Zhou, L., Shi, Y., Zhang, Z.-Y., & Wang, J. 2021, A&A, 653, A10
## Appendix A: Other models

In this appendix we show the results for the ionized, neutral-atomic, and neutral-molecular gas models with a varying number of sectors and with a power- or normal-law distribution for some parameters.

<table><caption>**Table A.1.** I Zw 18 model parameters for other considered models.</caption><thead><tr><th>Parameter</th><th>M<sub>I,normal</sub></th><th>M<sub>I,power-law</sub></th><th>M<sub>NA,power-law</sub></th><th>M<sub>NM,power-law</sub></th></tr></thead><tbody><tr><td>U<sub>0</sub></td><td>[-3.34<sup>+0.58</sup><sub>-0.61</sub>, -0.62<sup>+0.51</sup><sub>-2.57</sub>]<br>α<sub>U</sub> = -1.45<sup>+2.15</sup><sub>-2.16</sub></td><td>[-3.31<sup>+0.67</sup><sub>-0.43</sub>, -1.11<sup>+0.90</sup><sub>-1.52</sub>]<br>α<sub>U</sub> = -1.57<sup>+1.51</sup><sub>-1.34</sub></td><td>[-2.98<sup>+0.46</sup><sub>-0.54</sub>, -0.42<sup>+0.34</sup><sub>-1.70</sub>]<br>α<sub>U</sub> = -2.23<sup>+3.59</sup><sub>-1.10</sub></td><td>[-3.04<sup>+0.33</sup><sub>-0.68</sub>, -0.54<sup>+0.35</sup><sub>-0.55</sub>]<br>α<sub>U</sub> = -2.66<sup>+0.77</sup><sub>-0.75</sub></td></tr><tr><td>n<sub>0</sub> [log cm<sup>-3</sup>]</td><td>[1.11<sup>+0.75</sup><sub>-1.03</sub>, 3.30<sup>+0.61</sup><sub>-0.81</sub>]<br>α<sub>n</sub> = -3.30<sup>+2.28</sup><sub>-1.67</sub></td><td>[1.01<sup>+1.32</sup><sub>-0.91</sub>, 3.13<sup>+0.79</sup><sub>-0.78</sub>]<br>α<sub>n</sub> = -3.48<sup>+1.74</sup><sub>-1.29</sub></td><td>[0.85<sup>+1.12</sup><sub>-0.48</sub>, 3.30<sup>+0.49</sup><sub>-1.35</sub>]<br>α<sub>n</sub> = -2.72<sup>+1.98</sup><sub>-1.29</sub></td><td>[1.89<sup>+0.42</sup><sub>-0.38</sub>, 3.36<sup>+0.46</sup><sub>-0.76</sub>]<br>α<sub>n</sub> = -2.63<sup>+1.14</sup><sub>-0.75</sub></td></tr><tr><td>cut<sup>a</sup></td><td>[0.38<sup>+0.46</sup><sub>-0.28</sub>, 0.74<sup>+0.27</sup><sub>-0.44</sub>]<br>μ<sub>cut</sub> = 0.38<sup>+0.90</sup><sub>-0.47</sub><br>σ<sub>cut</sub> = 0.77<sup>+1.75</sup><sub>-1.62</sub></td><td>[0.20<sup>+0.16</sup><sub>-0.19</sub>, 1.03<sup>+0.20</sup><sub>-0.20</sub>]<br>α<sub>cut</sub> = -0.81<sup>+1.33</sup><sub>-1.29</sub></td><td>[0.20<sup>+0.33</sup><sub>-0.15</sub>, 1.19<sup>+0.61</sup><sub>-0.55</sub>]<br>α<sub>cut</sub> = -0.37<sup>+2.15</sup><sub>-2.20</sub>, -2.35<sup>+3.62</sup><sub>-4.17</sub><br>x<sup>b</sup><sub>cut</sub> = 0.92<sup>+0.56</sup><sub>-0.58</sub></td><td>[0.14<sup>+0.21</sup><sub>-0.10</sub>, 1.46<sup>+0.76</sup><sub>-0.33</sub>]<br>α<sub>cut</sub> = -0.57<sup>+2.71</sup><sub>-2.97</sub>, -2.62<sup>+1.34</sup><sub>-1.05</sub><br>x<sup>b</sup><sub>cut</sub> = 0.74<sup>+0.34</sup><sub>-0.29</sub></td></tr><tr><td>Z [log Z<sub>☉</sub>]</td><td>-1.46<sup>+0.06</sup><sub>-0.08</sub></td><td>-1.45<sup>+0.08</sup><sub>-0.12</sub></td><td>-1.46<sup>+0.11</sup><sub>-0.06</sub></td><td>-1.40<sup>+0.06</sup><sub>-0.09</sub></td></tr><tr><td>Z<sub>dust</sub><sup>b</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>L<sub>X</sub> [log erg s<sup>-1</sup>]</td><td>40.06<sup>+0.21</sup><sub>-0.93</sub></td><td>40.07<sup>+0.33</sup><sub>-0.14</sub></td><td>40.12<sup>+1.37</sup><sub>-0.24</sub></td><td>40.13<sup>+0.53</sup><sub>-0.19</sub></td></tr><tr><td>T<sub>X</sub> [log K]</td><td>5.92<sup>+0.94</sup><sub>-1.25</sub></td><td>6.04<sup>+0.45</sup><sub>-2.54</sub></td><td>5.90<sup>+0.52</sup><sub>-0.33</sub></td><td>5.80<sup>+0.70</sup><sub>-1.69</sub></td></tr><tr><td>Q(H) [log s<sup>-1</sup>]</td><td>52.24<sup>+0.40</sup><sub>-0.23</sub></td><td>52.26<sup>+0.44</sup><sub>-0.28</sub></td><td>52.25<sup>+0.28</sup><sub>-0.15</sub></td><td>52.18<sup>+0.68</sup><sub>-0.18</sub></td></tr><tr><td>f<sub>esc</sub>(1 &minus; &infin;) [%]</td><td>54.04<sup>+31.41</sup><sub>-33.20</sub></td><td>55.33<sup>+31.06</sup><sub>-28.47</sub></td><td>58.10<sup>+19.49</sup><sub>-21.28</sub></td><td>54.56<sup>+39.26</sup><sub>-17.41</sub></td></tr><tr><td>- f<sub>esc</sub>(1 &minus; 1.8) [%]</td><td></td><td></td><td></td><td></td></tr><tr><td>- f<sub>esc</sub>(1.8 &minus; 4) [%]</td><td></td><td></td><td></td><td></td></tr><tr><td>- f<sub>esc</sub>(4 &minus; 20) [%]</td><td></td><td></td><td></td><td></td></tr><tr><td>- f<sub>esc</sub>(20 &minus; &infin;) [%]</td><td></td><td></td><td></td><td></td></tr><tr><td>M<sub>H+</sub> [log M<sub>☉</sub>]</td><td>6.54<sup>+0.82</sup><sub>-0.93</sub></td><td>6.52<sup>+1.03</sup><sub>-1.04</sub></td><td>6.54<sup>+0.92</sup><sub>-1.01</sub></td><td>5.54<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>M<sub>H0</sub> [log M<sub>☉</sub>]</td><td>3.95<sup>+1.04</sup><sub>-1.11</sub></td><td>4.43<sup>+3.67</sup><sub>-1.54</sub></td><td>7.42<sup>+0.90</sup><sub>-3.69</sub></td><td>7.48<sup>+0.49</sup><sub>-0.29</sub></td></tr><tr><td>M<sub>H2</sub> [log M<sub>☉</sub>]</td><td>-7.32<sup>+4.49</sup><sub>-1.20</sub></td><td>-4.83<sup>+12.08</sup><sub>-4.66</sub></td><td>5.09<sup>+2.38</sup><sub>-13.36</sub></td><td>5.70<sup>+1.54</sup><sub>-0.98</sub></td></tr><tr><td>- M<sub>H2,C+</sub> [log M<sub>☉</sub>]</td><td>-7.87<sup>+4.91</sup><sub>-1.66</sub></td><td>-4.99<sup>+11.68</sup><sub>-5.45</sub></td><td>5.09<sup>+1.83</sup><sub>-14.35</sub></td><td>5.65<sup>+1.33</sup><sub>-0.92</sub></td></tr><tr><td>- M<sub>H2,C</sub> [log M<sub>☉</sub>]</td><td>-10.48<sup>+4.71</sup><sub>-2.24</sub></td><td>-7.98<sup>+14.80</sup><sub>-5.33</sub></td><td>2.01<sup>+4.93</sup><sub>-14.33</sub></td><td>3.04<sup>+2.49</sup><sub>-3.49</sub></td></tr><tr><td>- M<sub>H2,CO</sub> [log M<sub>☉</sub>]</td><td>-24.56<sup>+72.41</sup><sub>-4.41</sub></td><td>-15.85<sup>+23.40</sup><sub>-11.95</sub></td><td>-1.33<sup>+25.71</sup><sub>-14.02</sub></td><td>1.22<sup>+3.47</sup><sub>-3.00</sub></td></tr><tr><td>M<sub>dust</sub> [log M<sub>☉</sub>]</td><td>1.30<sup>+0.77</sup><sub>-0.93</sub></td><td>1.27<sup>+1.59</sup><sub>-0.88</sub></td><td>2.42<sup>+0.62</sup><sub>-1.23</sub></td><td>2.42<sup>+0.42</sup><sub>-0.27</sub></td></tr><tr><td>f<sub>CII,H+</sub></td><td>0.99<sup>+0.00</sup><sub>-0.05</sub></td><td>0.96<sup>+0.04</sup><sub>-0.85</sub></td><td>0.08<sup>+0.52</sup><sub>-0.07</sub></td><td>0.01<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>f<sub>CII,H</sub></td><td>0.01<sup>+0.05</sup><sub>-0.00</sub></td><td>0.04<sup>+0.78</sup><sub>-0.04</sub></td><td>0.86<sup>+0.12</sup><sub>-0.86</sub></td><td>0.98<sup>+0.01</sup><sub>-0.06</sub></td></tr><tr><td>f<sub>CII,H2</sub></td><td>0.00<sup>+0.00</sup><sub>-0.00</sub></td><td>0.00<sup>+0.07</sup><sub>-0.00</sub></td><td>0.00<sup>+0.07</sup><sub>-0.00</sub></td><td>0.01<sup>+0.06</sup><sub>-0.01</sub></td></tr></tbody></table>

**Notes.** For each parameter *U*, *n*, and *cut*, the first line shows the lower and upper boundaries. $Z_{dust}$ was fixed to zero (i.e., corresponding to the median relationship in Galliano et al. 2021).

<footnote>(a) 1: ionization front (H⁺–H⁰ transition), 2: PDR front (H⁰–H₂ transition), 3: C⁰–CO transition. (b) Pivot point for the broken power law.</footnote>
<figure><img src="image_27.png" /><figcaption>**Fig. A.1.** Results for model *M*<sub>Iopt+IR</sub>. See Figure 9 for the plot description.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. A.2. Results for model $M_{NA}$. See Figure 9 for the plot description.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. A.3.** Results for model *M*<sub>NM</sub>. See Figure 9 for the plot description.</figcaption></figure>