© D. Cormier et al. 2019

1. Introduction

The interstellar medium (ISM) of galaxies is a complex environment composed of phases with inhomogeneous structures and different physical conditions, evolving through dynamical effects such as feedback from massive stars or turbulence (e.g., Hennebelle & Chabrier 2008; Hopkins et al. 2012; Naab & Ostriker 2017). Characterizing the ISM in galaxies is necessary to understand how dense structures are formed and how the environmental properties (metal content, activity, stellar mass) may affect the star-formation process and thereby galaxy evolution. Space and airborne missions in the infrared, with IRAS, ISO, KAO, Spitzer, Herschel, and SOFIA, have been fundamental for unveiling a large part the energy budget of galaxies and enabling for a better census of the cooling of the ISM. Advances in sensitivity have allowed detection of emission lines probing regions related to star formation in low-metallicity galaxies for the first time (e.g., Hunter et al. 2001; Madden et al. 2006; Wu et al. 2006; Hunt et al. 2010; Cormier et al. 2015; Cigan et al. 2016). These galaxies stand out compared to more metal-rich galaxies. In particular, the prominent [O III] line at 88 μm in low-metallicity star-forming regions highlights the presence of an extended ionized gas phase, while the lack of CO emission suggests an extremely clumpy dense gas distribution (e.g., Indebetouw et al. 2013; Cormier et al. 2014; Lebouteiller et al. 2017; Turner et al. 2017; Vallini et al. 2017; Jameson et al. 2018). Those observations provide evidence of a possibly different, porous ISM structure in those galaxies as opposed to more metal-rich galaxies. Establishing whether there is such a metallicity-dependent change in the ISM structure is very important for understanding the characteristics of low-mass galaxies in the high-redshift Universe. It is particularly relevant for galaxies at the epoch of reionization because enhanced porosity could be a precondition that facilitates the escape of ionizing photons from those galaxies (Stark 2016). However, a quantitative and systematic characterization of the topology and composition of the ISM in galaxies has yet to be achieved. Because of the complexity of the ISM and the fact that current observations do not resolve individual ISM components in external galaxies, such a characterization can currently only be reliably achieved via multiphase modeling.

Emission lines in the mid-infrared (MIR) and far-infrared (FIR) are powerful diagnostics of the ISM properties (e.g., Díaz-Santos et al. 2013; Cormier et al. 2015; Fernández-Ontiveros et al. 2016; Lapham et al. 2017; Spinoglio et al. 2017; Ucci et al. 2017; Herrera-Camus et al. 2018a; Polles et al. 2019; van der Tak et al. 2018) because, contrary to optical lines, they probe a wider range of physical conditions (temperatures, densities, pressure, etc.), including more embedded regions. In addition, their observed intensities are less affected by extinction. The various lines that are accessible in the MIR and FIR may arise from different ionized and neutral gas phases that can be heated by several mechanisms such as ultraviolet photoionization from young stars, photoelectric effect on dust grains, X-rays, and shocks (Tielens 2005). Among those lines is the [C II] at 157 μm, which is one of the brightest cooling lines in galaxies, used as a star-formation rate (SFR) and mass tracer (e.g., De Looze et al. 2011, 2014; Herrera-Camus et al. 2015; Vallini et al. 2015; Jameson et al. 2018; Zanella et al. 2018), but with uncertainties linked in part to the possibility of [C II] being emitted in diffuse ionized gas as well as in photodissociation regions (PDR) with moderate densities.

In a first paper, Cormier et al. (2015), we present observations of the FIR fine-structure lines observed with the PACS instrument onboard Herschel (Poglitsch et al. 2010; Pilbratt et al. 2010) in the galaxies of the “Dwarf Galaxy Survey” (DGS; Madden et al. 2013). We found strikingly higher FIR line ratios of [O III]₈₈/[N II]₁₂₂, [N III]₅₇/[N II]₁₂₂, [C II]₁₅₇/[N II]₁₂₂, and [O III]₈₈/[O I]₆₃ in dwarf galaxies compared to more metal-rich galaxies. Along with observations from the Spitzer IRS spectrograph (Werner et al. 2004; Houck et al. 2004), we interpret those ratios in terms of radiation field hardness and high filling factor of ionized gas relative to the neutral gas. However, our analysis in that paper is focused on the average properties of the dwarf galaxy sample and the model assumptions are quite simplified (for example, we assume a single set of abundances). The DGS sample is also biased toward dwarf galaxies with high specific star-formation rate (sSFR) for detectability reasons. This bias does not always make it easy to infer whether differences with metal-rich galaxies are driven by their higher sSFR or the lower metallicity. Our goal in this paper is to perform a more detailed modeling of the ISM properties (physical conditions, porosity, etc.) of the DGS survey by modeling each galaxy of the sample and by tuning some of the model parameters to more appropriate values (from previous studies or ancillary data) whenever possible. We aim to investigate correlations of parameters characteristic of the ISM structure and conditions with a galaxy metallicity (Z) and star-formation activity, and, ultimately, to infer what drives changes in the ISM of galaxies, which is important to understand how we can expect galaxies in the early Universe to be affected by their low metallicities.

The outline of the paper is as follows. In Sect. 2 we present the galaxy sample and dataset used. Section 3 lays out the modeling strategy, grid assumptions, and the method to compare models and observations. We present the model results in Sect. 4 and their sensitivity to some of the choices in the modeling strategy in Sect. 5. Finally, we discuss the origin of [C II] emission as well as the ISM porosity in Sect. 6. Our results are summarized in Sect. 7.

2. Data

2.1. Sample selection

Our sample is drawn from the Herschel Dwarf Galaxy Survey (Madden et al. 2013) which consists of 48 nearby galaxies at distances between 0.5 and 200 Mpc. In the closest galaxies, only specific regions were observed with the PACS spectrometer while the more distant galaxies were observed in their entirety. We have focussed our study on the galaxies of the compact subsample (43 galaxies) that were observed in full (thus excluding NGC 2366, NGC 4861, and UM 311), resulting in 40 galaxies. Among those, 37 have ancillary Spitzer IRS spectroscopic data (all but HS 0017+1055, UGC 4483, and UM 133). We further require that our galaxies are detected in at least three spectral lines (IRS or PACS) to proceed with the modeling (thus excluding HS 1236+3937 and HS 2352+2733). This leaves us with 38 galaxies. Details on the data reduction can be found in Cormier et al. (2015). Most galaxies were observed in both high-resolution (HR) and low-resolution (LR) modes with the IRS instrument. Since many lines are better detected in HR, we used mainly fluxes from this mode. In cases where the LR observation is deeper and provides a detection in LR while the line is undetected in HR, we use the LR flux. We also verify that fluxes in LR and HR are consistent. Several lines are close to each other in wavelength and blended spectrally in the LR mode (H₂12.28 μm and Huα 12.36 μm, [O IV] 25.89 μm (not used) and [Fe II] 25.99 μm). For those, we consider only the HR fluxes. For convenience, we provide a table with all the PACS and IRS line fluxes used for the modeling at the CDS.

We note that NGC 4214 (D = 3 Mpc) is quite extended and has not been observed in full in spectroscopy. The physical properties of its two main star-forming regions have been studied in Dimaratos et al. (2015) using both IRS and PACS data. Here we have included those two regions in our sample because they have (rare) observations of the [N II] 122 μm line as well as of the [N II] 205 μm line from the SPIRE FTS instrument. Those are important constraints for the origins of [C II] emission discussed in Sect. 6.1. We refer the reader to Dimaratos et al. (2015) for fluxes and details on the modeling of the ionized and PDR gas in this object.

Continuum measurements for our sample have been obtained with Spitzer and Herschel. We used the total infrared luminosities (L_TIR), derived by modeling the dust spectral energy distributions (SEDs), from Rémy-Ruyer et al. (2015). We also used bolometric luminosities, L_BOL, as measured from the continuum SEDs presented in Rémy-Ruyer et al. (2015), extrapolating the curve in the UV/optical range with a slope of −0.6 (in νL_ν versus λ) and integrating the SED from 0.1 to 1000 microns.

2.2. Properties of the main observed lines

The standard ionic lines observed by the IRS and PACS are relatively bright in low-Z galaxies. Their properties are summarized in Table 1. The main ionic lines used in this study to constrain models are: [Ar II] 6.99 μm, [Ar III] 8.99 μm, [S IV] 10.51 μm, Huα 12.36 μm, [Ne II] 12.81 μm, [Ne III] 15.56 μm, [S III] 18.71 μm, [Ar III] 21.83 μm, [S III] 33.48 μm, [N III] 57.32 μm, [O III] 88.36 μm, and [N II] 121.9 μm. These lines are found in the ionized gas only, and usually trace rather dense H II regions. The ratios of [S IV]/[S III], [Ne III]/[Ne II], [Ar III]/[Ar II], [N III]/[N II] are indicators of the hardness and strength of the radiation field. They are independent of elemental abundances and less dependent on density than optical lines. The ratios of the two [S III] lines or [Ar III] lines measure the electron density in the H II regions for densities in the range 100 − 10 000 cm⁻³ and 3000 − 300 000 cm⁻³ respectively as they are not very sensitive to the electron temperature nor to extinction (Houck et al. 1984; Keenan & Conlon 1993; Osterbrock & Ferland 2006). The H I recombination line Huα is very useful because it is a good tracer of the total ionized gas mass. It can be used to calibrate elemental abundances from MIR lines or to measure the extinction in the optical by comparison to Hα (e.g., Hummer & Storey 1987; Lebouteiller et al. 2008).

Other ionic lines such as [Fe III] (22.93 μm) and [Fe II] (17.93 μm and 25.99 μm) are present in the IRS spectral range but their detection rate is lower compared to the main lines described above. The abundance of iron in the gas phase is enhanced in the presence of shocks. As a result, we have not used those iron lines to constrain the models but we compared a posteriori predictions from the models to observations in Sect. 4.4.1. Many ionic lines are also detected in the optical. However, optical observations often do not match our infrared spatial coverage. Optical lines are also sensitive to extinction and generally do not probe embedded regions as well as the infrared lines. Again, we have not used those optical lines to constrain our models. Instead, we compared optical line ratio predictions from the models to observations a posteriori in Sect. 4.4.3.

The main lines observed by the IRS and PACS tracing the neutral ISM are: [Si II] 34.82 μm, [O I] 63.18 μm, [O I] 145.5 μm, and [C II] 157.7 μm. The ionization potentials of Si⁰ and C⁰ being below that of hydrogen, both [Si II] and [C II] can be found outside of H II regions, in the neutral phase. [O I] is only found in neutral gas and usually arises from warm dense regions. [O I] tracks the high-A_V layer of PDRs better than [C II]. This is because the transition of C⁺ into C⁰ (∼1−3 mag) occurs before the formation of CO occurs (∼5 mag), and because at high-A_V, almost all the carbon is converted into molecules while some oxygen remains in atomic form. The ratio of the two [O I] lines is an indicator of the gas temperature for T ≃ 300 K. It is a density tracer for high temperatures and high densities (Hollenbach & Tielens 1999). In our modeling, the main constraints for the PDR are the [C II], [O I] lines and the infrared continuum radiation. We have not used [Si II] as a constraint since its abundance and emission can be linked to shocks, but we compared [Si II] observations to model predictions a posteriori. Mid-IR H₂ rotational lines are also detected in nine of our galaxies. However, models are very sensitive to the treatment of H₂ and shocks can also contribute to this emission. Hence, we also compared H₂ observations to model predictions a posteriori (see Sect. 4.4.1).

3. Methodology

Herschel and Spitzer spectral data provide powerful diagnostics on the ISM physical conditions. The observables that we have in hand trace various conditions and phases of the ISM, including diffuse ionized or neutral gas, compact H II regions, and dense PDRs. In order to interpret them correctly, a self-consistent multiphase modeling is required. In particular, we are interested in fitting the ionic and neutral gas lines simultaneously to account for the contribution from all phases modeled to the line emission and to better constrain conditions at the ionized-neutral phase boundary when few spectral lines are observed. Our goal is to identify a parameter space (density, ionization parameter, neutral gas covering factor) that fits all our tracers, thus requiring one or two model components (see Sect. 4).

3.1. Grid initial conditions

We used the spectral synthesis code Cloudy v.17.00 (Ferland et al. 2017) to model the ionized and neutral gas phases of our galaxies. The geometry is 1D plane-parallel. The main model input parameters are explained below and listed in Table 2.

Since the modeling results are sensitive to the gas and dust abundances, we built five grids of models with the following metallicity values: Z_bin = [0.05, 0.1, 0.25, 0.5, 1]. Z_bin = 1 corresponds to the default set of ISM abundances in Cloudy, for which 12 + log(O/H) = 8.5. For the other metallicities, we scaled the gas and dust abundances by Z_bin. The abundance of several elements (C, N, Ne, S, Si, Ar, Fe, Cl) were further scaled to measured patterns in low-metallicity galaxies according to Izotov et al. (2006). Table 2 reports the adopted abundances of those elements. For the dust, we used grain properties similar to those of the SMC and included polycyclic aromatic hydrocarbons (PAHs) which are important for neutral gas heating by photoelectric effect (Abel et al. 2008). The dust-to-gas ratio (DGR) is 5 × 10⁻³ for Z_bin = 1 and decreases linearly with metallicity. The abundance of PAHs was further scaled with Z_bin to the power 1.3, as determined by Rémy-Ruyer et al. (2015). The implications of these assumptions on the DGR and PAH abundance are discussed in Sect. 5.5.

The source of radiation is chosen as a young starburst that we simulated with Starburst99 (Leitherer et al. 2010). We opted for a continuous star formation SED computed out to an age of 10 Myr. The stars were selected from the standard Geneva tracks (without rotation) and distributed on a Kroupa IMF (slope −1.3 between 0.1 and 0.5 M_⊙ and −2.3 between 0.5 and 100 M_⊙; Kroupa 2001). The Starburst99 spectra are shown in Fig. 1 for the different Z_bin values. We also added a small soft X-ray component (blackbody peaking at ≃1.5 keV) to the input radiation field and adopt a cosmic ray rate of 6 × 10⁻¹⁶ s⁻¹ (see Sect. 5.2.2 for a justification of those parameters). We adopted a density profile that is roughly constant in the H II region and increases linearly with the hydrogen column density in the neutral gas (see Fig. 2, and Table 2). This density profile is an intermediate choice between the more extreme, common assumptions of a constant density or of constant pressure (that we test in Sect. 5). In the following, unless stated otherwise, n_H refers to the initial density. It is the density in the H II region. All models were stopped at an A_V of 5 mag to ensure going deep enough in the cloud for complete PDR predictions. We note that this A_V corresponds to different N(H) values for different Z values. Stopping the models at a fixed N(H) value instead of a fixed A_Vxs value would not change our results as long as N(H) is chosen large enough so that the PDR lines ([C II], [O I]) are fully produced. The transitions from the ionized to the atomic and from the atomic to the molecular regions are defined where the fractions of e⁻ and H⁰, and H⁰ and H₂ respectively, are equal to 0.5. We note that in environments with high ionization parameters, radiation pressure could become dynamically important if it dominates the other pressure terms. This could be considered an unstable model. For our models, we have verified that this is not the case.

Fig. 1. Spectra from Starburst99 used as input for Cloudy. Spectra are normalized and shown for each metallicity bin value of our grid. There are more hard photons at low metallicity.

Fig. 2. Adopted hydrogen density profile in the default grid of models as a function of visual extinction, shown for Zbin = 0.5 (solid line) and Zbin = 0.1 (dot-dash line). We also show the density profile for the constant pressure test performed in Sect. 5.4 (gray solid line). Model parameters are set to nH = 102 cm−3 and log U = −2. The phase transitions from H+-dominated to H I-dominated and from H I-dominated to H2-dominated are indicated with vertical lines. We note that AV of 1 mag corresponds to a hydrogen column density that is 5× larger at Zbin = 0.1 than at Zbin = 0.5.

We generated a grid of models for each of the five metallicity bins by varying the following parameters: the initial density (n_H) from 10^0.5 to 10^3.5 cm⁻³, and the inner radius such that we cover an ionization parameter (U) from 10⁻⁴ to 10⁻¹. Values of n_H and U are those at the start of the calculation. Since the PDR conditions are set by the H II region model (by continuity), there is little degree of freedom on the PDR parameters. The radiation field is set by that in the H II region and the density is set by the adopted density profile. As a free parameter we chose the covering of the PDR with respect to the H II region which always has a covering factor of unity. The PDR covering factor (cov_PDR) is not varied within the grid but a posteriori, from zero to unity, because it just corresponds to a linear scaling of the PDR intensities. We note that with this setup, the H II region is always radiation bounded, in other words no ionizing radiation escapes. In the case of a low PDR covering factor, segments of the H II region that are not covered by neutral gas could be stopped at lower A_V, before the ionization front is reached (making the model matter bounded with ionizing radiation escaping). However, this stopping parameter can only be constrained by specific lines tracing the ionization front (e.g., [Ar II], [Ne II]), and it introduces degeneracies with other parameters (see Polles et al. 2019, for a detailed analysis with this method). Hence we prefer not to add an additional free parameter and to keep all models radiation bounded. For the following discussion, we use the term porosity when referring to a reduced covering factor of the dense neutral gas (PDR) relative to the ionized gas as determined by the modeling. For convenience, we provide line luminosities from the grid of models predicted at the ionization front and at an A_V of 5 mag at the CDS.

Most of the choices for the model parameters described above are motivated by observations or set to default values. However, uncertainties in those choices can remain. In Sect. 5, we change some of the input parameters (namely the radiation field, the presence of cosmic rays and X-rays, the turbulent velocity, the density law, and the dust abundance) and we discuss their influence on the results.

3.2. Comparison to observations

3.2.1. Observed and adopted model abundances

Elemental abundances have been measured from optical spectroscopy in several of our galaxies. Their values are reported in Table A.1. Oxygen abundances are estimated with the Pilyugin & Thuan (2005) calibration. We refer the reader to Rémy-Ruyer et al. (2014) for further details and references on abundances and to Madden et al. (2013) and De Vis et al. (2017) for a comparison of metallicities obtained with other methods in the DGS galaxies. Since abundances in the model grid are fixed, a correction needs to be applied when comparing models to observations of a given galaxy. For small abundance variations, abundances, and intensities scale linearly. For each galaxy, we select the grid of models with metallicity (oxygen abundance) closest to the observed value (tabulated in Tables 3 and 4 as Z_bin). Then, for elements with measured abundances, we multiplied the predicted line intensities by the offset between the observed abundances (Table A.1) and the adopted model abundances (Table 2). For example, for a galaxy with log(O/H) = − 4.40, we select the grid in the metallicity bin Z_bin = 0.1 and multiply the oxygen line intensities by 10^{−4.40 + 4.50} = 1.26. Oxygen abundances are measured for all of our galaxies. For the other elements, if not measured, their abundances are still assumed to follow the patterns in Table 2, hence we interpolate values in that table considering the observed log(O/H) value. The abundance values quoted here can be uncertain and vary within galaxies or according to the datasets and methods used to measure them. For our galaxies, literature studies indicate that those variations are on the order of ∼0.2 dex. We see similar offsets between the measured abundances and the abundance patterns adopted in Table 2. Hence we took into account an uncertainty from abundances that we set to 50% of the observed line intensities.

3.2.2. Observed and modeled luminosities

To compare the absolute fluxes observed to model predictions, we needed to scale the luminosity of the model to that of the source. We calculated a nominal scaling factor as the ratio of the average observed luminosity to the average predicted luminosity of the main ionic lines to fit. In addition to reproducing the line luminosities, the successful scaled model should also reproduce the continuum luminosity (L_TIR) because this corresponds to the amount of stellar light reprocessed by dust in the PDR. Contrary to L_TIR which is used as a formal constraint on the models, the bolometric luminosity of the galaxies (L_BOL) is not used as a direct constraint. In principle, the input luminosity of the models should be equal to the bolometric luminosity. However, the correct value of the input luminosity in case of partial PDR coverage is difficult to assess. By energy conservation, the input luminosity of the models must be equal to L_TIR (if cov_PDR = 1) or greater than L_TIR (in case of escaping photons, cov_PDR <  1). It should be, at most, the bolometric luminosity of the galaxy (i.e., all energy is associated with the model stellar cluster). L_BOL is compared a posteriori to the input luminosity of the models to identify large mismatches that could, for example, indicate sources of heating other than stars. We identify a limited number of such cases (see Sect. 4.4.2).

3.2.3. Goodness of fit

We searched for the best-fitting models by minimizing the χ² defined as

(1)

where M_i and O_i are the modeled and observed intensities, and σ_i is the fractional error on the observed flux (uncertainty/flux) with calibration and abundance uncertainties added in quadrature to the measured uncertainties. A line is considered not detected when its flux is lower than three times its measured uncertainty. For non-detections with a 1σ upper limit σ_ul, we set O_i to M_i if M_i ≤ 3σ_ul (i.e., no contribution of the upper limit to the χ²) and we set O_i to 3σ_ul and σ_i to 1 if M_i >  3σ_ul. The number of free parameters in the models is three (density, ionization parameter, and PDR covering factor). The reduced χ² is given by where ν is the difference between the number of constraints and the number of free parameters.

3.2.4. Reliability of the derived physical parameters

The χ² minimization gives an idea on the models that are closest to the observations. However, it is possible that those models do not converge to a single value of the free parameters that we aim to constrain. To assess how well the free physical parameters are constrained, we compute their probability distribution functions (PDFs). The marginalized PDF of a given parameter is calculated at each fixed value of this parameter. It is taken as the sum of the probabilities (e^−χ2/2) of all the models having this parameter set to its fixed value and all possible values for the other parameters. The PDFs are then normalized such that the sum of probabilities of all models in the grid is equal to one. We estimated uncertainties on the best-fit parameters by considering the average and standard deviation of the values of the free parameters taken by models with χ² values within a certain range of . For instance, for a number of degrees of freedom of four, and a 1σ confidence level (probability of 68.3%), this range is (Press et al. 1992).

4. Results

We present general grid results in Sect. 4.1 and results from fitting the grid models to the observations in the following subsections.

4.1. Grid results

Figure 3 shows the parameter space covered by the models and the observations for Z_bin = 0.1 and Z_bin = 0.5. In the left panels, we see the values of G₀ produced when varying the ionization parameter and the density, which are input parameters of the model grid. We remind the reader that G₀ is the intensity of the ultraviolet radiation field at the PDR front, given in units of the background Habing (1968) value, 1.6 × 10⁻³ erg cm⁻² s⁻¹, as defined by Tielens & Hollenbach (1985). For a given density, U and G₀ scale linearly (except at the high-U, low-n_H end due to geometric dilution of the field; see e.g., Cormier et al. 2015). The middle panels of Fig. 3 show diagnostics of the H II region. The [S III] line ratio traces intermediate densities (n_H ≃ 10² − 10^3.5 cm⁻³), independently of Z. Model predictions for the [S IV]/[S III] line ratio are slightly shifted up at Z_bin = 0.1 compared to Z_bin = 0.5, especially for low-U values, because the starlight radiation field is harder. Observations span values between 0.1 and 10 for the [S IV]₁₀/[S III]₁₈ line ratio and between 0.4 and 1.2 for the [S III]₁₈/[S III]₃₃ line ratio. Both grids at Z_bin = 0.1 and Z_bin = 0.5 cover the observations quite well. Observations are rather compatible with models of low/ intermediate-densities and high-U values. The right panels of Fig. 3 show a diagnostic linking the H II region and the PDR. We see that predictions of the [O III]/[C II] and [C II]/[O I] ratios vary significantly with Z. The main effects at work are the C/O abundance ratio that decreases by a factor of two from Z_bin = 0.5 to Z_bin = 0.1, the harder stellar radiation field at Z_bin = 0.1 making more [O III], and the fact that the PDR structure is very sensitive to Z. In particular, with decreasing Z, the PDR lines are formed up to a larger column density and their emission is boosted relative to the ionized gas metal lines (i.e., predictions going down in Fig. 3, right panels). Observations span values between 0.3 and 10 for the [O III]₈₈/[C II]₁₅₇ line ratio and between 0.3 and 5 for the [C II]₁₅₇/[O I]₆₃ line ratio. At Z_bin = 0.5, models cover quite well the observed ranges, while at Z_bin = 0.1, we notice that model predictions are at the edge of the observed ratio ranges. This suggests that a combination of models and/or low covering factor of the PDR may be needed to match better the observations.

Fig. 3. Grid results: predictions of selected line intensities and model parameters, for Zbin = 0.5 (top panels) and Zbin = 0.1 (bottom panels). Color scheme: increasing initial density from 100.5 cm−3 (light green) to 103.5 cm−3 (dark purple). Symbol size: large symbols correspond to small inner radius (large U) and small symbols to large inner radius (small U). Dotted lines connect models with the same inner radius, solid lines connect models with the same density. All observations are overplotted in orange and galaxies in the metallicity bin plotted are highlighted in red. We note that the error bars are large because they include abundance uncertainties (assumed to be 50% of the flux). The PDR covering factor is set to unity. Top right panel: quantitative indication on how the C/O abundance ratio affects the grid predictions from Zbin = 0.5 to Zbin = 0.1. Other effects due to metallicity (e.g., stellar spectra, column density, etc.) are discussed in the text.

We remind the reader that normal, star-forming, metal-rich galaxies have similar [S III]₁₈/[S III]₃₃ and [C II]₁₅₇/[O I]₆₃ ratio values than the DGS galaxies, but their [S IV]₁₀/[S III]₁₈ and [O III]₈₈/[C II]₁₅₇ values are respectively about seven times and four times lower (Brauher et al. 2008; Dale et al. 2009; Cormier et al. 2015). Those values are covered by our grid of models and compatible with models of low or intermediate-densities and intermediate-U values (for Z_bin = 1).

4.2. Best-fitting single H II+PDR models

Parameters of the best-fitting single models as well as their corresponding χ² and reduced χ² () are indicated in Tables 3 and 4. The number of detected lines in each galaxy varies from three (model under-constrained) to 16. Figure 4 shows how the observed and predicted line intensities compare for the best-fitting model of He 2-10, as well as the PDFs of n_H, log(U), and cov_PDR. Figures for the other galaxies of the sample are shown in Appendix D. The ionization parameters vary between −2.9 and −0.3 with median value −2.4, and the PDR covering factor varies from 0.2 to one with a median value of 0.4. Densities of the best-fitting models vary from 10^0.5 to 10^3.0 cm⁻³, with median at 10^2.0 cm⁻³. The medians are very close to the values log(U) = − 2.5 and n_H = 10^2.0 cm⁻³ found in Cormier et al. (2015) from average line ratios of the DGS sample. The classical MIR density tracers ([S III] and [Ar III] line ratios) are not particularly sensitive to those low densities, and although we can recover such densities by modeling a suite of lines, future determinations of electron densities in dwarf galaxies would benefit from FIR observations of line doublets such as the [O III]₅₂ and [O III]₈₈ or [N II]₁₂₂ and [N II]₂₀₅ lines.

Fig. 4. Model results, shown only for the example galaxy He 2-10 here and for all galaxies in Appendix D. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

About half of the galaxies (16/39) have minimum below two and their observed lines are satisfyingly reproduced by a single model. Inspection of the results indicates that galaxies with larger usually have one or more lines for which predictions are off by a factor of more than two. The lines that are most difficult to reproduce simultaneously with the other lines are [S IV] and [Ne II].

4.3. Mixing two model components

4.3.1. Motivation

It is possible that the minimum χ² for single models as described above does not provide a satisfying fit to the data because our model description is too simple. In addition to single models, we explored the possibility of mixing two models of the grids: one dense high-U model (model #1), which should represent the material that gave birth to the starburst, and one more diffuse lower-U model (model #2), which should represent a hypothetical more extended ionized phase. This two-component modeling approach is supported by earlier work on those DGS galaxies suggestive of a large filling factor, diffuse, ionized gas phase (Cormier et al. 2012, 2015; Dimaratos et al. 2015). For this, we have generated a new set of models by combining two models for which parameters are varied within the parameter space of our original grid, and their relative contribution f_c (in terms of luminosity) varies from 10% to 90% in steps of 10%. That is, for a given line, the intensity prediction of the combined model is the sum of the prediction of the first model multiplied by f_c and of the prediction of the second model times 1 − f_c. Given the small number of constraints for the PDR, we consider that only the dense H II model produces a PDR. In that case, we have as free parameters for the mixed models: n_H (but n_H, 1 ≥ n_H, 2), U (but U₁ ≥ U₂), and cov_PDR, 1 (with cov_PDR, 2 = 0). The number of free parameters is five. Figure 5 shows an illustration of the best-fitting model for one example galaxy and figures for each galaxy of the sample can be found in Appendix B. We stress that for mixed models, f_c and cov_PDR are different quantities. f_c corresponds to the contribution of a given model to the final (2-component) model, it is therefore translated into a covering factor of the H II region in Fig. 5. cov_PDR is the covering factor of the PDR relative to that of the H II region, hence the covering factor of the final model is cov_PDR × f_c.

Fig. 5. 2D representation of the 1D best-fitting (mixed) model for the example galaxy He 2-10. Schematics for the other galaxies of the sample are shown in Appendix B. The model parameters are the same as those reported in Table 4. Color coding corresponds to density. The central stellar cluster is represented in green. The distance to the inner black shell is proportional to the inner radius, and the thickness of the inner black shell is proportional to the ionization parameter U, while the thickness of the H II and PDR regions are kept fixed. U is a function of the input stellar spectrum, density, and inner radius. The angle shown by dotted lines indicates the contribution from each component (i.e., parameter fc in Table 4).

For galaxies with fewer than six lines detected (for which a mixed model would be under-constrained) or with from the single models lower than two (or lower than with the mixed models), we consider the single best-fit model to be the best representation. There are 23 galaxies that satisfy this condition (results reported in Table 3). For the other 16 galaxies, we consider the mixed best-fit model to be better (results reported in Table 4).

4.3.2. Results for mixed models

In order to reproduce all observed line intensities in galaxies where a single model fails (), we search for a combination of two model components. We stress that the second component is not added to the single model but a new set of two model components is found (with some constraints on those parameters as explained in the previous subsection). Parameters of the best-fitting mixed models as well as their corresponding χ² are indicated in Tables 3 and 4. Similar figures as for the single model case are shown in Fig. 4 and in Appendix D for all galaxies.

In general, the lines that were not satisfyingly reproduced with single models are well fit by mixed models. The low-U component contributes mostly to the [N II], [Ne II], [Ar II], and [S III] lines. We note that the two U-values of the mixed models always bracket the U-value of the single models, while the n_H-values of the mixed models are sometimes lower than the n_H-value of the single models (for example, NGC 625, SBS 1533+574). For those galaxies, the single models tend to under-predict the [C II]₁₅₇/[O I]₆₃ ratio by a factor of two to three while the lower n_H-values of the mixed models reproduce this ratio better. Overall, our observations seem to be more sensitive to variations in the ionization parameter, although there might be also some degeneracy at play between density and covering factor for the PDR lines. The PDFs of the parameters of the mixed models are sometimes broader than for the single models, especially for the lower-density component, indicating that it is not very well constrained by the observations. In most cases, the denser component contributes more than the lower-density component to the final mixed model (f_c >  0.5).

For most galaxies for which single models do not provide a good fit, the mixed models perform better (lower ). All galaxies for which mixed models are preferred have . For several galaxies though, mixed models do not improve the results over the single models but provide similar (II Zw 40, Mrk 209, Pox 186, SBS 0335−052, SBS 1159+545, SBS 1415+437) or worse (Tol 1214−277). For two of them (SBS 1159+545, Tol 1214−277), remains above three. For those two galaxies, we find that the combination of a low-U, high-n_H component going in the PDR and a high-U, low-n_H component would actually fit the data better (with close to unity). For the galaxies where a single model was already giving a good fit, the mixed models provide lower χ² but are often under-constrained. The parameters of the two model components are sometimes very similar to the parameters of the single model or they bracket them (especially for U and less so for n_H, as noted above).

Comparison to previous, similar studies. At this stage, it is worth comparing our findings to two previous studies that we carried out to model MIR and FIR lines using a similar approach. Results for Haro 11 are in global agreement with previous modeling done in Cormier et al. (2012). For the H II region, we find similar densities (10^2.8 cm⁻³ and 10^1.0 cm⁻³ versus 10^3.0 cm⁻³ and 10^1.5 cm⁻³ here). The U-values are somewhat different (−2.5 and −1.7 versus −1.9 and −3.9). This is mainly driven by the different stellar spectra used in input (instantaneous burst and single star versus continuous star-formation scenario here). For the PDR, we find a larger PDR covering factor because the PDR density is lower than in Cormier et al. (2012) where the density was set by pressure equilibrium between the H II region and the PDR. For the two star-forming regions of NGC 4214, we also find very similar results for single models to Dimaratos et al. (2015), with lower density, ionization parameter, and PDR covering factor in the central region compared to the southern region. Mixed models (not included in Dimaratos et al. 2015) are however better suited to reproduce the [N II] emission in NGC 4214.

4.4. A posteriori checks

4.4.1. Mid-infrared lines

We have not used the silicon, iron, and H₂ lines from the Spitzer IRS instrument to constrain the models due to the ambiguous origin (possibly shocks) of those lines and uncertain gas-phase abundances. However, we have compared predictions from the best models a posteriori. We report on the performance of the single and mixed models in predicting those secondary lines in Table 5.

The [Si II] line at 34 μm can originate in the ionized gas and in the neutral ISM. It is quite bright and detected in many of the galaxies of our sample. We find that the [Si II] line is systematically under-predicted by a factor of two to six by the best-fitting models. Mixed models perform slightly better than single models although [Si II] remains under-predicted in most cases. Since our models include the H II region and a diffuse, low-excitation ionized phase in the case of mixed models, any contribution of ionized gas to [Si II] emission is largely taken into account by our models. Thus, reasons for the under-prediction of [Si II] could be additional excitation of [Si II] by shocks or uncertainty in the adopted gaseous abundances. Given the difficulty in modeling iron (see below), we cannot really compare our results for [Si II] with those for [Fe II] to rule out shock excitation.

The MIR iron lines are faint and detected in only a few galaxies of our sample. The [Fe II] line at 25.99 μm is generally well reproduced while the [Fe III] line at 22.93 μm is over-predicted by a factor of two to three for half of the sources with detection, even when adding those lines in the list of lines to fit for. We note that collisional excitation of [Fe II] by H and H₂ is not included in the version of Cloudy used in this work, hence an additional contribution from PDRs to the [Fe II] emission (not included in our models) would then over-predict the observations. However, the main effect affecting the line prediction is again the adopted abundances. Although measured for many of our galaxies (Table A.1), abundances can be uncertain because of unknown ionization corrections, depletion of iron in the cold dense medium (Savage & Sembach 1996), or enhancement of iron from supernovae-driven shocks (O’Halloran et al. 2008). In our case, adopted abundances may be systematically too high. At solar metallicity, Kaufman et al. (2006) consider an iron abundance ten times lower than our adopted abundance.

The H₂ rotational lines are also faint and detected in few galaxies in our sample. The S(1) transition at 17.03 μm is sometimes over-predicted and sometimes under-predicted by a factor of two to fifteen, while the S(2) transition at 12.28 μm is usually under-predicted by a factor of three to ten. The ratio of those two transitions is mostly sensitive to density for G₀ >  10² (Kaufman et al. 2006). When fitting for the H₂ lines, the densities required for the PDR are indeed larger (typically ≥10² cm⁻³ at the PDR front) but increases a lot (models performing worse for some H II lines, [C II], [O I]). More H₂ detections are required to really investigate the presence and importance of shocks.

4.4.2. Bolometric luminosity

As explained in Sect. 3.2.2, the total luminosity of the models is not an easy parameter to constrain. Here we just verify that the total luminosity of the best-fitting models does not exceed the observed bolometric luminosities. Table 5 shows that this condition is satisfied for all galaxies in the case of single models. In the case of mixed models, for 12/37 galaxies, the luminosity required by the models is larger than the observed L_BOL by a factor of three to four in most cases and up to 30 for one galaxy. For all those 12 galaxies, the diffuse ionized component contributes to the majority of the total luminosity, and for ten of those galaxies, the single models are actually preferred. For the two galaxies where mixed models are preferred (HS 0052+2536 and Mrk 1089), L_BOL is over-predicted by a factor of three. When interpreting this over-prediction, one should not forget that there is a range of acceptable values for each model parameter (PDFs in Fig. 4 and in Appendix D). Indeed, a density of 10 cm⁻³ (value often within the acceptable range) instead of 3 cm⁻³ for the diffuse component of those two galaxies would reconcile better the total observed and modeled luminosities. We also note that we calculated the observed bolometric luminosity relatively coarsely.

4.4.3. Optical lines

Since the optical lines are very sensitive to extinction and to the distribution of dust in the models, as opposed to the infrared lines, and their observations often cover a smaller field-of-view, we do not include them in finding the best-fitting models. A posteriori comparisons indicate that the intensity ratios of the optical lines relative to Hβ are generally well reproduced by the models (within a factor of three). To assess how the best-fitting models would perform for the total fluxes (i.e., not normalized to Hβ), we also compare the ratio of total observed Hα luminosity (references in Rémy-Ruyer et al. 2015, mainly from Gil de Paz et al. 2003) to the [O III]₈₈ luminosity. We find that the observed Hα/[O III]₈₈ ratio is matched by the models (within a factor of two) for more than half of the galaxies of the sample, and it is over-predicted for the other galaxies by a factor of three to eight. The discrepancies with optical observations are likely caused by extinction/geometry considerations and reproducing them, simultaneously with IR lines, is not trivial (and neither is it the goal of this study; see also Polles et al. 2019).

4.5. Contributions to the [C II] and [Si II] emission

The use of the [C II] line at 157 μm as a reliable neutral gas/PDR tracer has been questioned since it can also originate from diffuse (n_{crit, e⁻, [CII]}  ∼  50 cm⁻³), low-excitation ionized gas. To quantify such contribution to the observed [C II] intensity, it is often compared to the [N II] lines at 122 μm or 205 μm, which have low ionization potentials and low critical densities. In metal-rich environments, the fraction of [C II] emission arising in the H II region ranges from 5% to 60% (e.g., Malhotra et al. 2001; Pineda et al. 2013; Parkin et al. 2014; Hughes et al. 2015; Abdullah et al. 2017; Croxall et al. 2017; Díaz-Santos et al. 2017; Lapham et al. 2017). The low [N II]₁₂₂/[C II]₁₅₇ values measured in the dwarfs (Cormier et al. 2015) suggest little contribution of the H II region to the observed C⁺ emission. Similarly, [Si II] can arise from low-excitation ionized gas but it has higher critical density with electrons than [C II] (n_{crit,e⁻,[Si II]} ∼ 1000 cm⁻³). The H II region produces significant [Si II] emission, more than the PDR, at moderate densities and low-U (see also Abel et al. 2005). With the physical conditions constrained in most of the DGS galaxies, we can quantify more precisely how much [C II] and [Si II] emission is predicted to arise from the ionized gas. The definition that we use to measure this is to take the intensity predictions at the phase transition from ionized to atomic gas (i.e., where the fractions of e⁻ and H⁰ are equal to 0.5).

In the case of single models, we find that the fractions of [C II] and [Si II] predicted in the H II region are small, typically less than 10% for [C II] and less than 40% for [Si II]. In the case of mixed models, those fractions generally increase slightly but remain below 40% for [C II] and 50% for [Si II]. Regardless of density, it is the low-U component which contributes most to those fractions. The fractions of [C II] and [Si II] emission from the ionized gas are generally low in the single models because most ionic lines that we observed or detected require high-U values. A possible additional low-U component will increase these fractions but the parameters of this component can only be constrained when [N II] and, to a lesser extent, [Ne II] and [Ar II] are available. We further discuss the fractions of [C II] from the ionized gas with other galaxy parameters and line diagnostics in Sect. 6.1.

5. Influence of some modeling choices

In the methodology described in Sect. 3, we have made choices regarding the set of spectral lines to fit as well as the set of model parameters to fix or vary. Not all of the input parameters can be varied and fit by our observations. We have taken care to use the best input fixed parameters and the motivation of their choices for several of them is given below. In this section, we further explore: (1) how the inclusion or exclusion of specific lines can change the best-fitting results (Sect. 5.1); and (2) the influence of varying one of the fixed parameters at a time (Sects. 5.2–5.5). For (2), we focus on the metallicity bins Z_bin = 0.5 and Z_bin = 0.1, and we make comparison to galaxies in those bins (including Haro 11 which we study in detail in Cormier et al. 2012).

5.1. Set of spectral lines to fit

As discussed in the previous section, spectral lines that are known to be problematic for the adopted modeling strategy were excluded from the fit (for instance if shock heating or geometry effects can be important). Now, within the default set of fitted lines, we aim to determine the effect of excluding some specific lines on the best-fitting results and to identify possible biases. For this, we have considered excluding, one at a time, the [C II]₁₅₇ line which can arise from more diffuse gas than [O I], the [O I]₆₃ line which can be optically thick, the [Ne III]₁₅ line which can be problematic to model from population synthesis models (see Sect. 5.2), and the [S IV]₁₀ line which is usually well detected but one of the lines that is least well reproduced by our models. We also consider including the H₂ lines in the fit and fitting a small, fixed set of lines ([C II]₁₅₇, [O III]₈₈, [Ne III]₁₅, and L_TIR) since several of the low-Z galaxies are only detected in those lines.

The results are shown in Fig. 6. We plot histograms of the best-fitting parameters and values for both single and mixed models. We note that all galaxies are included and uncertainties on the best-fitting parameters as well as quality of the fits are not conveyed on the plots. For single models, seems to improve quite significantly when excluding [S IV]. For mixed models, the distributions of remain similar and the mean values are only slightly lower when excluding a specific line. Compared to the default case, the distribution of model parameters for the tests performed have relatively similar shapes, and individual galaxies typically move to the neighboring bin at most. Excluding [O I]₆₃ or [Ne III]₁₅ does not have a significant impact on the best-fitting parameters, while the distribution of those parameters change noticeably when excluding [S IV] or [C II]. Indeed, by excluding [S IV], there are fewer galaxies with best-fitting models that have very low density (n_H = 10^0.5 cm⁻³), hence the mean density of the mixed, low-density component is larger. The mean ionization parameters of the single and mixed, high-density components are also slightly lower, probably because [S IV] has one of the highest ionization potentials. Similarly, when considering the small, fixed set of fitted lines (which only allows for single models, not mixed models), the distribution of the best-fitting solutions does not change significantly but the uncertainties on the model parameters are larger. There are fewer galaxies at the highest densities, probably because the [S III] lines are not available to constrain further the density and the mean ionization parameter is also slightly lower because [S IV] is excluded from the fit. We want to stress that the trends of best-fitting parameters with global galaxy properties (SFR, Z) discussed in Sect. 6 persist with this test. By excluding [C II], the distributions of densities shift to larger (factor of three on average) values and the distribution of PDR covering factors also shifts to slightly lower values. This can be understood as the PDR solution is driven by L_TIR and the [O I] lines which generally trace the denser, lower filling factor regions of PDRs than [C II]. Including the warm H₂ lines in the fit does not affect the overall distributions and mean values of the model parameters because those lines are detected in only eight galaxies of our sample. Those lines are difficult to reproduce though, yielding larger values on average.

Fig. 6. Sensitivity of the model results to the fitted lines. We show histograms of the reduced χ2 () and best-fitting parameters (nH, log U, covPDR) for the entire sample of galaxies. The results from the default grid of models and the default set of lines are in black, and the results are in color if the set of fitted lines is changed, that is, when one spectral line is excluded or when the H2 lines included or when we consider a small, fixed set of lines to fit ([Ne III]15, [O III]88, [C II]157, and LTIR). The vertical dotted lines indicate the mean values. Each histogram is offset vertically by a value of 20. See Sect. 5.1 for details.

5.2. Choice of input radiation field

The sources of radiation giving rise to the observed line emission are multiple and mainly clustered in a few star-forming regions. Given the limited geometry of Cloudy, we consider only one central source per model and need to choose its radiation field. Several input SEDs were tested (Fig. 7). In this section, we discuss the effects of different starburst spectra and heating by X-rays and cosmic rays on the line predictions and best-fitting models.

Fig. 7. Cloudy input radiation fields tested. We show the input (dashed) and output (solid) spectral energy distributions for: a continuous star-formation scenario of 10 Myr (with an X-ray component), a Starburst99 instantaneous burst of 3.5 Myr, a Popstar burst of 3.5 Myr, a Magphys fit to the optical-IR photometry. All have model parameters: nH = 102 cm−3, rin = 1020.4 cm, and Zbin = 0.5.

5.2.1. Starburst spectrum

We tried as input the Starburst99 stellar spectrum both from an instantaneous burst with an age between 3 and 6 Myr, and from a continuous (up to an age of 10 Myr) star-formation scenario. In the instantaneous case, we find that the emission of lines tracing the radiation field hardness is too sensitive to the step size for sampling the burst age, making it easy to miss the “best” age (see also Polles et al. 2019). It is also probably not the best representation of the cluster age distribution in galaxies. Moreover, the PDR lines are better fit with a continuous star-formation scenario, in part because G₀ is too high with an instantaneous burst (fewer lower-energy photons). We also tried an instantaneous burst using the Popstar evolutionary synthesis models (Mollá et al. 2009). A different treatment of the stellar atmosphere can produce significant differences in the stellar spectrum for photons with energy > 40 eV, hence affecting predictions of lines such as [Ne III] (Morisset & Georgiev 2009). However, similar problems persist, as noted above for the Starburst99 instantaneous case.

For a more complex star-formation scenario, we also tried to constrain the shape of the incident spectrum by taking the optical photometry into account. Focusing on Haro 11, we used the SED code Magphys (da Cunha et al. 2008) to obtain an unattenuated spectrum and used that spectrum as input in Cloudy. However, for the lines that we fit (H II region and PDR) the results are not especially improved ( was not lower) than with the continuous scenario of Starburst99. Moreover, since the available optical photometry is limited for some of our galaxies, we cannot use Magphys for all galaxies of our sample. It is important to have a consistent method in order to study trends. Therefore, we opt for a continuous star-formation scenario for our grid of models.

5.2.2. Heating of the PDR and [O I] line ratio

The standard model contains a soft X-ray component of luminosity ∼10⁻⁴ times the bolometric luminosity and cosmic rays. This is primarily needed for understanding the [O I] line emission. The [O I] 145 μm line is clearly detected in 10 galaxies and the [O I] 63 μm line in 30 galaxies of the sample, with [O I]₁₄₅/[O I]₆₃ luminosity ratios ranging from 0.06 to 0.09 (Cormier et al. 2015). The [O I] 63 μm line is known to suffer from optical depth effects or self-absorption in regions of the Milky Way and dusty galaxies (e.g., Liseau et al. 2006; Abel et al. 2007; Rosenberg et al. 2015), giving rise to [O I] line ratios above 0.1. The emergent line ratio is somewhat sensitive to the depth (A_V) of the models (Abel et al. 2007). Other heating sources than the stellar radiation field can also affect the PDR chemistry and line emission, such as cosmic rays or X-ray or shock heating. We explored the effect of high column density (high A_V), cosmic rays, and soft X-ray heating in Fig. 8.

Fig. 8. Model predictions for the [O I]145/[O I]63 luminosity ratio as a function of visual extinction in the cloud. The default model (solid curve) has parameters: nH = 103 cm−3, G0 = 103, Zbin = 0.5. Increasing the cosmic ray rate by a factor of 10 (dashed curve) or adding a soft X-ray component to the input radiation field (dotted curve) both increase the PDR heating and predicted [O I] line ratios. The gray shading indicates the observed range of ratios.

With no additional heating source than the stellar spectrum, the observed range of [O I] line ratios requires systematically high A_V, at which the CO optical depths are also very large, and this seems unlikely. Indeed, the mean A_V is expected to decrease with decreasing mean metallicity, and CO emission is known to be clumpy at low Z, suggesting that the [O I] emission integrated over entire galaxies is probably little affected by high-A_V (> 10 mag) sightlines.

Soft X-ray emission (0.2−2 keV) is detected in several galaxies of our sample, with an observed luminosity about 10⁴ times lower than L_TIR (except in I Zw 18 which is strong in X-rays; see Lebouteiller et al. 2017). The intrinsic (unabsorbed) emission should be higher, however we do not know how much of that emission goes into the heating of the dense ISM (as opposed to diffuse isotropic emission). Hence we decide to add a soft X-ray component to the input radiation field of our models that we take as a blackbody spectrum of temperature 5 × 10⁶ K (to peak around 1 keV) and luminosity in the soft X-ray range 10⁻⁴ L_BOL. Adding the X-rays creates a hybrid XDR/PDR model in Cloudy. X-rays generally do not affect the predictions of the inner ionized gas cloud but increases the heating at the ionization front and in the PDR (e.g., Meijerink & Spaans 2005; Lebouteiller et al. 2017). In the presence of strong X-ray flux, the deeper penetrating X-rays produce more Ne⁺ and Ar⁺, while the increased ionization has an effect on CO, suppressing its formation due to reactions with He⁺ and (e.g., Meijerink & Spaans 2005; Abel et al. 2009). This produces more free oxygen and a higher gas temperature. In our case, the X-ray component has low luminosity. The main effect is that the PDR temperature increases roughly from 150 K to 200 K. The emission of the PDR lines increases slightly, as well as the [Ne II] and [Ar II] emission for models with low-density or high-U. In the case tested in Fig. 8, the predicted [O I] line ratios go up by a factor of 1.2.

Regarding the cosmic rays, we show in Fig. 8 the effect of increasing the standard cosmic ray rate (2 × 10⁻¹⁶ s⁻¹; Indriolo et al. 2007) by an order of magnitude. Similarly to the case with soft X-rays, the PDR heating and line emission are enhanced. The cosmic ray rate varies significantly from sightline to sightline in the Galactic ISM (e.g., Dalgarno 2006; Neufeld & Wolfire 2017), and the dwarfs of our survey have undergone recent star formation, hence it seems reasonable to adopt a different value than the standard rate. To match the observed [O I] line ratios, we decided to increase the cosmic ray rate by a factor of three.

We note that cosmic rays and X-rays have similar effects (ionization and heating) and distinguishing between the two with no other constraints is not possible. One way to constrain the cosmic ray flux would be to study its effect on the ionization and coupling of the magnetic field and the gas in the PDR (Padovani & Galli 2011).

5.3. Turbulent velocity

The micro-turbulent velocity affects the shielding and pumping of lines. It is included as a pressure term in the equation of state. Hence, we expect that it impacts mainly predictions of PDR and high optical depth lines. We set the turbulent velocity v_turb to 1.5 km s⁻¹ in the default grid of models, as in Kaufman et al. (1999), and here, we explored the cases v_turb = 0 km s⁻¹ and v_turb = 5 km s⁻¹ for the grid of models in the metallicity bin Z_bin = 0.5. We find that, by increasing the turbulent velocity from 0 km s⁻¹ to 5 km s⁻¹, the PDR temperature reduces marginally by about 10%, the [O I] line intensities are down by 15%, the [C II] intensity is up by 20%, and the H₂ lines are down by ∼30%. Overall, the turbulent velocity does not have a significant impact on our resulting best-fitting models and therefore on our trend analysis.

5.4. Density law

One of the model choices with the largest consequence is the density profile to adopt. As default, we adopted a density profile that increases smoothly with depth (see also Hosokawa & Inutsuka 2005; Wolfire et al. 2010). Here we have tested, for the grid of models in the metallicity bin Z_bin = 0.5, two other often-used options: a constant density throughout the H II region and the PDR, and hydrostatic equilibrium. Hydrostatic equilibrium assumes that the total pressure (including gas, turbulent, radiation pressure terms) is constant and implies a jump in density at the transition between the H II region and the PDR (of almost two orders of magnitude when gas pressure dominates, see Fig. 2). The adopted density law affects mostly the PDR predictions with the constant density models predicting less [O I] emission and the constant pressure models predicting more [O I] emission than the default grid. For Haro 11, the dense component of the mixed model prefers higher initial densities in constant density and lower initial densities in constant pressure. However, we have tested the constant density and constant pressure hypotheses on eight other galaxies at Z_bin = 0.5 (Haro 2, Haro 3, II Zw 40, Mrk 1089, NGC 625, NGC 1705, NGC 5253, UM 448), and those trends are not systematic. In general, we find that the constant density/pressure single models do not provide better fits to the observations. With mixed models, the constant pressure case provides (very marginally) lower than the default grid for eight of those nine galaxies. Overall, the smoothly increasing density law (adopted as default) is a reasonable choice for single models. It also provides satisfying fits for mixed models, although constant pressure models could also be considered but there is no clear evidence for the need to switch to constant pressure.

5.5. Abundance of PAHs and grains

Grains and PAHs are important for the energy balance, in particular for the cooling of the neutral ISM and heating via the photoelectric effect. The adopted abundances of grains and PAHs are motivated by recent analysis of the Spitzer and Herschel data (Rémy-Ruyer et al. 2014, 2015). However, there is both large scatter at all metallicities investigated (0.03 ≤ Z ≤ 3) and a limited number of observations at Z ≤ 0.2. For example, the relation linking the DGR with metallicity seems to become steeper at Z <  0.2 than the relation calibrated on Z ≥ 0.2 galaxies (e.g., Rémy-Ruyer et al. 2014). Concerning the PAHs, their abundance in the PDRs of dwarf galaxies is very uncertain as they are scarcely detected for Z <  0.2 with the Spitzer IRS instrument and we expect a large variation of the PAH-to-dust ratio between H II regions and PDRs. We note that, if dust is present in the ionized gas but PAHs are weak there, a lower PDR covering factor at low Z (see Sect. 6.3) would significantly influence the globally observed PAH-to-dust emission ratios. Constraints on PAHs from the James Webb Space Telescope (JWST) will be very valuable. Departure from the default values and thus the effect on the results will be larger at low Z. Hence we carried out two tests for Z_bin = 0.1 (and not Z_bin = 0.5 as in the previous subsections): (a) we adopt a PAH abundance 10× lower; (b) we adopt a DGR 10× lower than in Table 2.

Results for the galaxies Mrk 209, Pox 186, SBS 1159+545, SBS 1415+437, Tol 1214−277, UM 461, VII Zw 403 (all at Z_bin = 0.1) are inspected. Test (a) does not change any of the results, probably because the PAH abundance was already low by default hence the PAHs contribute very little to the energy balance of those galaxies. Dust grains are responsible for the heating. Test (b) affects the energy balance noticeably. By reducing the DGR by a factor of ten and therefore the dust grain abundance, the photoelectric effect remains the main heating source in the PDR but it is less important, and the PDR temperature goes down from ≃150 K to 80 K. The cloud is also more transparent and the stopping criterion of A_V of 5 mag corresponds to a larger depth. The emission of the [C II] and [O I] lines increases by a factor of ∼2.5 and 1.4, respectively. This grid with lower DGR at Z_bin = 0.1 still does not cover the highest observed [C II]/[O I]₆₃ ratios (lower, right panel of Fig. 3. We find that the values are lower, especially for SBS 1159+545 and Tol 1214−277, and marginally for three other galaxies. The best-fitting models have slightly higher densities for four of the seven galaxies in this metallicity bin, but the best-fitting solutions do not change much and there are no clear systematic trends for all galaxies (in particular for the PDR covering factor which is discussed later).

In conclusion, we find that additional heating sources are not required to reproduce the PDR emission, as opposed to what was found by Lebouteiller et al. (2017) for I Zw 18, where X-rays are mainly responsible for the PDR heating and not grains. For this particular galaxy, our adopted DGR of ≃2 × 10⁻⁴ might be too high since Lebouteiller et al. (2017) adopt a value of ≃7 × 10⁻⁶. We also conclude that the abundance of PAHs is not an important parameter for the models (as long as the dust grain abundance is larger). The dust grain abundance (and DGR) does, however, change the temperature structure of the cloud. For our study, this does not have a noticeable impact on the results, but additional PDR observations (such as MIR H₂ lines for which line ratios are expected to change) and better constraints on the DGR would be needed for a more thorough study of the PDR properties.

6. Discussion

6.1. Origin of [C II] emission

All of the ionic lines used in this study have critical densities for collisions with electrons larger than that of C⁺, making it difficult to constrain a low-density, low-excitation ionized phase from which [C II] emission could arise. Optical lines such as [S II] are often used as density tracers, but they probe a similar range of densities as [S III], i.e., 100–10 000 cm⁻³, as well as shocks (Osterbrock & Ferland 2006). [N II]₁₂₂ is the best line that we have in hand but, in principle, [N II]₂₀₅ is better suited (see Table 1 for critical densities). Croxall et al. (2017) and Díaz-Santos et al. (2017) carried out studies on the origin of [C II] emission in nearby late-type spiral and luminous galaxies based on [N II]₁₂₂ and [N II]₂₀₅. In those galaxies, [N II] line ratios that are indicative of electron densities between 1 and 300 cm⁻³ with median values at 30 and 41 cm⁻³, (Herrera-Camus et al. 2016 and Díaz-Santos et al. 2017, respectively; see also Herrera-Camus et al. 2018a).

Because of the issue of constraining the low-density ionized gas, it is interesting to compare our results on densities and fractions of [C II] emission from the ionized gas with those obtained from theoretical calculations based on the [N II] lines. We have also explored the sensitivity of those parameters to a fixed-density ionized gas component in the models.

Inferred electron densities from the [N II] line ratio only. [N II]₁₂₂ is observed in seven of our galaxies and [N II]₂₀₅ in only three of our galaxies: Haro 11, He 2-10, and NGC 4214-c. The [N II]₂₀₅/[N II]₁₂₂ ratios are observed to be 0.67, 0.35, and 0.4, respectively and are matched by models with densities between 30 and 80 cm⁻³.

Comparison of observed [N II] emission with models. The single models often tend to under-predict [N II]₁₂₂ while the mixed models match [N II]₁₂₂ observations better. In Haro 11, He 2-10, and NGC 4214-c, we also fit for the [N II]₂₀₅ line and the best mixed models reproduce successfully both [N II] lines, which emission arises from the low-U (log(U)≃ − 3.9), low-density component (n_H of 30, 30 and 100 cm⁻³, respectively, values close to those inferred from the [N II] line ratio only).

Forcing a low-density component in the models. For galaxies without [N II] observations, we still would like to estimate how much [C II] emission could arise from such low-density, low-excitation ionized phase. Here, we assume n_e ≃ 30 cm⁻³. We derived the fraction of ionized [C II] by forcing the single models or one component in our mixed models to have a density of 30 cm⁻³ and an ionization parameter lower than that of the other component.

The results on the [C II] emission arising from the ionized gas are presented in Fig. 9. For the three galaxies with [N II]₂₀₅ detected, we first estimated the fraction of [C II] emission from the ionized gas from theoretical calculations, based on statistical equilibrium of collisionally excited gas (see, e.g., Oberst et al. 2006 and Croxall et al. 2017). Those fractions range from 13% to 17% (triangles in Fig. 9) and are close to those derived from our models (circles and crosses in Fig. 9), though on the high side for NGC 4214-c. For all galaxies, we then compared estimates from our models with free n_H and with fixed, low n_H. In all cases, we find that the fraction of [C II] emission coming from the ionized phase is at most 50%. In galaxies where [N II] is observed, the models predict lower fractions of [C II] in the ionized phase than predictions from Croxall et al. (2017) because our observed [N II]/[C II] ratios are on the low side of most of their observations. Predictions based on assuming a fixed low-density component (n_H = 30 cm⁻³) can be unreliable because they sometimes yield much higher χ² than when n_H is let free. Those predictions are not plotted in Fig. 9 if they yield . Based on our models, we also see a trend of decreasing fraction of [C II] emission arising in the ionized gas with decreasing Z. The trend is similar to that reported by Croxall et al. (2017) and extends to lower Z values. For Z ≤ 0.2, our observations agree better with predictions from Kaufman et al. (2006, their Eq. A9 for a low-density medium) than with predictions from Croxall et al. (2017) (extrapolated to lower Z). Following Kaufman et al. (2006), we fit the DGS data only as well as the DGS and Croxall et al. (2017) data with a function of the form

(2)

Fig. 9. Fraction of [C II] emission coming from the ionized gas phase. We show predictions for each DGS galaxy based on: (1) the best single or mixed models reported in Tables 3 and 4 (crosses), and (2) the best single or mixed models with one model having a fixed low density of 30 cm−3 (circles); predictions yielding (all models are bad) are not plotted. Red symbols indicate that [N II]122 is detected and orange symbols indicate that both [N II]122 and [N II]205 are detected. For the three galaxies where [N II]205 is detected, we also show the predicted values from simple theoretical calculations, as done in Croxall et al. (2017). We overlaid the data (diamonds) and fit from Croxall et al. (2017) as well as the prediction from PDR models from Kaufman et al. (2006) in gray. Fits to the DGS data and to the DGS (free density) and Croxall et al. (2017) data are shown in blue.

where Z is the metallicity normalized to the solar value (i.e., 12 + log([O/H]_⊙) = 8.7). The metallicities of both datasets are based on the Pilyugin & Thuan (2005) calibration. For the DGS galaxies and with n_e ≃ 30 cm⁻³, we find coefficients (a, b)≃(0.3, −2.3) (blue dotted curve in Fig. 9). For the DGS and DGS+Croxall et al. (2017) galaxies with free n_e, we find coefficients (a, b)≃(1.2, −1.4) (solid blue and dash-dotted blue curves in Fig. 9).

In Fig. 10, we present correlations of the fractions of [C II] emission from the ionized gas, f[CII]_{ionized gas}, obtained from our models with n_H free, with metallicity, sSFR and several combinations of observed line ratios to see if it is possible to predict the [C II] fractions more easily (from observables rather than modeling). We quantified correlations using the IDL procedure r_correlate.pro which outputs the Spearman’s rank correlation coefficient and the two-sided significance of its deviation from zero. Here, we vary quantities within their uncertainties a 1,000 times to obtain a distribution of correlation coefficients. The plots report the median value ± the standard deviation of the distribution, as well as the median value of the significance distribution in parenthesis. As seen in Fig. 10 (and already noted in Fig. 9), f[CII]_{ionized gas} is correlated with metallicity. f[CII]_{ionized gas} is not correlated with sSFR. Regarding observed line ratios, we find the strongest trends with [C II]/[Ne II], [Ne III]/[Ne II], and [S IV]/[S III]₁₈. We also find strong trends of f[CII]_{ionized gas} with [C II]/[N II]₁₂₂ and [O III]/[N II]₁₂₂ when taking into account predictions of [N II]₁₂₂ in addition to observations. The trends with [Ne III]/[Ne II], [S IV]/[S III]₁₈, and [O III]/[N II]₁₂₂ can be understood as those ratios trace the average radiation field strength and, indirectly, metallicity. Regarding the [O III]/[N II]₁₂₂ ratio, we overlay in Fig. 10 the fit of Díaz-Santos et al. (2017) for local luminous galaxies. We find a similar but offset trend, the dwarfs probing higher ratios of [O III]/[N II]₁₂₂. We do not find any trend of f[CII]_{ionized gas} with the [S III] line ratio, which is a density tracer, though sensitive to higher densities than the critical density of [C II] with electrons. From our model results, we find a weak anticorrelation between electron densities and f[CII]_{ionized gas}, but with a lot of scatter. Hence, as in Díaz-Santos et al. (2017), f[CII]_{ionized gas} correlates more strongly with the radiation field intensity than with the electron density (conveyed in the color coding of the last two panels of Fig. 10, and see also Fig. C.1).

Fig. 10. Correlation of the fraction of [C II] emission from the ionized gas with: metallicity, specific star-formation rate, and line ratios for which the strongest trends are found. Values for the y-axis are taken from the models in boldface in Tables 3 and 4 (triangles for single models, filled circles for mixed models). Color coding corresponds to metallicity, except for the last two panels which are color coded by electron density and ionization parameter (of component #2 in case of mixed models). Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner (ρ). In the last three panels, star symbols are the ratios of observed [C II]157 or [O III]88 intensity divided by the predicted [N II]122 intensity. The correlation coefficients ρall include those star symbols.

Accurso et al. (2017) and Olsen et al. (2017) investigate the [C II] emission based on simulations of galaxies. In general, they find larger fractions of [C II] in the ionized gas than we find here. Equation (16) of Accurso et al. (2017) provides a prescription of the fraction of [C II] associated with PDRs that depends (only) on a galaxy’s sSFR. This prescription predicts that 55–75% of the [C II] emission in our galaxies should arise from the PDR, which is somewhat coherent with our results. However, their Equation 12, which also depends on the dust mass fraction, Z, and electron density that we set to 30 cm⁻³, brings those fractions below 50%. Olsen et al. (2017) also find that the [C II] emission from PDRs is below 50% in simulations of high-redshift star-forming galaxies. Those fractions of [C II] emission arising from the PDR are low compared to our results and to Croxall et al. (2017). This discrepancy might be related to the fact that a two-phase model is a simplistic view of a galaxy, or alternatively that prescriptions or simulations are not representative of our population of galaxies, each galaxy being dominated by a few localized star-forming events.

Finally, our results may have implications for the use of [C II] as a SFR indicator. Most of the [C II] emission arises in PDRs at low metallicities, which indicates that the [C II]₁₅₇ line could still be a good tracer of the SFR in the low-Z galaxies at high redshift (see also Stacey et al. 2010; Vallini et al. 2015; Lagache et al. 2018). De Looze et al. (2014) found an offset in the [C II]–SFR relation between nearby dwarf galaxies, and normal star-forming spirals and IR-luminous galaxies. This offset might be linked to contamination of the ionized gas to the observed [C II] emission, resulting in higher [C II] level of emission for a given level of SFR in metal-rich galaxies. Given our results and those of Croxall et al. (2017) on higher-Z galaxies, the contamination of the ionized gas to the observed [C II] emission remains however moderate. Another important factor to consider in the [C II]–SFR relation is the [C II] line deficit (with respect to L_TIR) which varies by orders of magnitude in galaxies. This line deficit is closely related to the surface density of star formation and therefore probably to local variations of the dense gas conditions (Smith et al. 2017).

6.2. Correlation of modeled physical conditions, metallicity, and SFR

We aim to see if we can link the physical conditions, namely the density and radiation field in the H II regions and PDRs, to more global parameters such as Z and sSFR. Figure 11 shows the model parameters of each galaxy with color codes corresponding to Z and sSFR.

Fig. 11. Distribution of model parameters. Top row: ionization parameter U and density nH for the best-fitting single models (circles) or mixed models (component #1: squares and component #2: triangles). For mixed models, the black lines link the two model components and the symbol size is proportional to the scaling factor fc of each component. Bottom row: same for the radiation field G0 impinging the PDR and the average PDR density nH, PDR. Color coding corresponds to Z (left panels) and sSFR (right panels). Galaxies with no sSFR value are displayed with open symbols.

Regarding the model parameters of the grid (top panels in Fig. 11), we find that there are no obvious systematics in the distribution of mixed models. For several galaxies, the mixed models highlight a need to have components with different U rather than different n_H. There is a tendency for the low-U models to also have smaller scaling factors. Galaxies with higher Z tend to have lower ionization parameters than galaxies with lower Z (see also Fig. C.1). Low-to-moderate ionization parameters (log(U) <  − 2) are indeed found in more metal-rich nearby galaxies (e.g., Herrera-Camus et al. 2018b). It is in the dwarf galaxies with lower Z that we find the highest U values (log(U) >  − 2). Finally, galaxies with higher sSFR tend to have higher n_H than galaxies with lower sSFR.

Regarding the PDR parameters (bottom panels in Fig. 11), we find a clear correlation between G₀ and n_H, PDR. We also find a tendency for lower-Z galaxies to have higher G₀ and n_H, PDR values, but that tendency is clearer for high sSFR galaxies. Malhotra et al. (2001) conducted a study on 60 nearby normal star-forming galaxies. They also found a relation between G₀ and n, with G₀ ∝ n^1.4, and explain that such scaling is expected from H II regions (Strömgren spheres) surrounded by PDRs. The ratio of the two, G₀/n_H, PDR, drives the heating of the PDR (e.g., Tielens & Hollenbach 1985) and, interestingly, it is roughly constant in the galaxies of our sample. We find a median value of 0.5 and a standard deviation of 0.4 dex. This is in the same range of values found by Díaz-Santos et al. (2017) for IR-luminous galaxies. High sSFR galaxies also have slightly higher G₀/n_H, PDR values than the median, indicating more compact/intense star-forming regions, consistent with the increase of G₀/n_H, PDR at high IR surface densities (though for more moderate values of G₀/n_H, PDR and Σ_IR in our sample) reported by Díaz-Santos et al. (2017). Analytically, given our adopted density law, G₀/n_H, PDR is proportional to U × Z. Since we find that lower-Z galaxies tend to have higher-U values, this naturally leads to a roughly constant G₀/n_H, PDR ratio. However, a constant ratio is not purely a consequence of model assumptions because the PDR parameters have to fit the (PDR) observations. A different density law would lead to a similar relation.

6.3. Phase filling factors

One main result from analyzing Spitzer and Herschel observations of both resolved regions and entire galaxies is that the ionized phase is extended and fills a large volume of the ISM of star-forming dwarf galaxies (Cormier et al. 2012; Chevance et al. 2016; Polles et al. 2019). With the modeling that we have performed for each individual galaxy of the Herschel DGS, we can quantify further the relative filling factor of the ionized and PDR phases. First, we investigated the relative volume filling factor of the two ionized gas components in the case of mixed models. We could expect a diffuse, low-excitation ionized medium to fill a larger volume than denser H II regions. This is indeed the case for galaxies where the two components have very different densities (for example, Haro 11 and NGC 5253 for which the diffuse phase fills a volume that is about ten and one hundred times larger than the dense phase, respectively). However, as discussed in the previous section (Fig. 11, top panels), we do not find a systematic contrast in density for the mixed models, even in galaxies where [N II]₁₂₂ is detected. In particular, the [O III] emission, often found extended in studies that have sufficient spatial resolution (Chevance et al. 2016; Polles et al. 2019), is well reproduced by the high-U (not necessarily low density) models.

Second, we investigate the covering factor of the PDR phase relative to that of the ionized phase. Figure 12 shows the PDR covering factors of the models as a function of metallicity, sSFR, and selected line ratios. The covering factor is simply the fraction of solid angle covered by PDR gas as seen from the center of the ionized region. For the mixed models, it is the product of the covering factor of the dense model multiplied by the scaling factor (the covering factor of the low-density model being set to zero by default). Correlating the PDR covering factor with several line ratios, such as [O III]₈₈/[C II]₁₅₇, [O III]₈₈/L_TIR, or the optical line ratio of [O III]5007 Å -to- [O II]3727 Å, does not yield strong trends. Although the dwarf galaxies are offset in the observed parameter space of those ratios compared to metal-rich galaxies (Cormier et al. 2015), we do not find a robust line ratio predictor of the PDR covering factor (i.e., ISM porosity). It seems that the emission of the lines under study is too dependent on the other physical conditions and only the combination of many constraints provided by various lines can predict such parameter. However, we find that there is a correlation between the PDR covering factor and both the broadband luminosity ratio of L_TIR/L_FUV and Z, with large scatter. The trend indicates that the porosity of the ISM (relative covering factors of ionized and PDR gas) increases at low Z.

Fig. 12. Correlations of the PDR covering factor with metallicity (top left panel), specific star-formation rate (top right panel), LTIR/LFUV ratio (bottom left panel), and [O III]88/[C II] ratio (bottom right panel). Galaxies are color coded by Z. Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner. Lower-Z galaxies tend to have smaller covering factors but we do not find strong trends with observed line ratios. The apparent floor of 0.2 in the PDR covering factor is due to the chosen steps for this parameter in the grid.

A decrease of the PDR covering factor may facilitate the escape of ionizing photons from the star-forming regions. Although our models by design cannot be used to quantify this escape, it is still a parameter worth discussing. In the Magellanic Clouds (with Z = 1/2 and 1/5 Z_⊙, where 12 + log(O/H)_⊙ = 8.69; Asplund et al. 2009), the fraction of ionizing photons escaping H II regions is of the order of 45% (Pellegrini et al. 2012). Moreover, Gazagnes et al. (2018) show for galaxies with known Lyman continuum leakage that the escape of ionizing photons is related to the covering factor of the H I gas (rather than low H I column densities). Recent modeling of the Local Group dwarf galaxy IC 10 also shows that on scales of individual star-forming regions, clouds tend to be matter-bounded, allowing for the possibility of photons escaping the regions, while on larger spatial scales, clouds tend to be more radiation-bounded (Polles et al. 2019). At large enough scales the escape fraction is expected to be near zero. However, given that our models are radiation-bounded (i.e., we stopped them further than the ionization front) and that several nearby dwarf galaxies have extended H I halos that we have not modeled, the PDR covering factors derived here cannot be translated into galaxy escape fractions. Such escape fractions have been measured to be up to a few percent in nearby dwarf galaxies (e.g., Leitet et al. 2011, 2013; see also Zastrow et al. 2013), but could be larger in dwarfs at higher redshift (Razoumov & Sommer-Larsen 2010; Yajima et al. 2011). Leitet et al. (2013) find that the escape fraction is somewhat larger in galaxies with lower metallicity, lower stellar mass and higher sSFR for the energetics of the star-forming regions to affect the whole galaxy (Clarke & Oey 2002). Here, we find correlations between the PDR covering factor and metallicity (as well as with stellar mass and SFR) but not with sSFR. A lower PDR covering factor in nearby low-Z galaxies might indicate that the escape fraction is likely to be higher for low-Z galaxies in the high-redshift Universe.

7. Conclusions

We have modeled the MIR and FIR emission of the ISM of galaxies from the Herschel Dwarf Galaxy Survey with the spectral synthesis code Cloudy. The goal of this study is to characterize the gas physical conditions in the H II region and PDR of those galaxies (namely density, ionization parameter, covering factor of the PDR), and to understand how or if they correlate with a galaxy’s metallicity or star-formation activity. We focus our discussion on the phase origin of [C II] emission and on the ISM porosity. Our results are summarized as follows.

• Several galaxies are well reproduced by single models, but for other galaxies we need mixed models, especially models with a different ionization parameter (rather than a different gas density) to reproduce satisfyingly the observed MIR-FIR emission. We find values of ionization parameters in the range log(U)≃ − 3.0 to −0.3, and electron densities in the range n_e ≃ 10^0.5 − 10^3.0 cm⁻³, corresponding to PDR densities roughly ten times larger (given the adopted density profile). The highest U values are higher than in metal-rich galaxies and are found in the lower-Z galaxies of our sample. Radiation fields and densities in the PDR are also found higher for the lowest-Z galaxies. Densities are higher in galaxies with higher specific star-formation rates.

• We evaluated the effect of the choice of spectral lines to fit and of individual input parameters of the models on the results. The ionization parameter is mainly sensitive to the choice of input radiation field library, while choices on the density profile and dust-to-gas ratio can affect the density in the PDR. We find the need for a low-luminosity soft X-ray component and/or increased cosmic ray rate to reproduce the observed [O I]₁₄₅/[O I]₆₃ line ratios.

• The contribution of the ionized phase to the [C II] emission is low, typically < 30%, making [C II] a good tracer of PDRs and of the SFR at low metallicities (and potentially at high redshift). As found in Croxall et al. (2017), there is a trend of increasing [C II] fractions from the ionized gas with metallicity, though with scatter. This fraction seems more driven by the ionization parameter (and indirectly metallicity) than by density.

• We find that the covering factor of the neutral gas relative to that of the ionized gas decreases with decreasing metallicity and TIR-to-FUV luminosity ratio. This provides evidence for a change in the ISM porosity, which we conjecture may facilitate the propagation and escape of ionizing photon in some systems. If such ISM porosity is present in low-metallicity galaxies at high redshift, this might have implications for the fraction of escaping photons and cosmic reionization.

The DGS galaxies span more than an order of magnitude in metallicity and sSFR with a bias toward high sSFR values. Metallicity (rather than sSFR) seems to be the main parameter driving variations in the ISM properties of the DGS galaxies. However, extending the modeling approach of this paper to galaxies that cover a larger parameter space in terms of metallicity, SFR, stellar mass, amongst others, will be important to map the evolution of ISM properties through cosmic times.

Our suite of models makes predictions for many more lines than observed and discussed in this paper. These predictions will be used to guide follow-up observations that will allow us to exploit unique PDR coolants, such as [C I] (available with ALMA), MIR H₂ and PAH features (available with the JWST), and to test the importance of specific heating mechanisms, such as shocks. Spatially resolving the main ISM phases would also be extremely useful for a more direct confirmation of the phase origin of the emission lines, respective filling factors, ISM porosity and potentially, the escape of ionizing photons.

Acknowledgments

The authors thank P. De Vis for providing the Magphys SED of Haro 11 and M. Krumholz for useful discussion. We also thank the referee for the valuable comments which helped to improve the manuscript. DC is supported by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 702622. NPA acknowledges funding from NASA and SOFIA through grant program SOF 05-0084. SH is supported by the DFG program HO 5475/2-1. IDL gratefully acknowledges the support of the Research Foundation – Flanders (FWO). We also acknowledge support from the DAAD/PROCOPE projects 57210883/35265PE and from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

References

Appendix A: Observed abundances for the DGS galaxies

In Table A.1, we report the observed abundances from the literature for the galaxies of the Herschel Dwarf Galaxy Survey (Madden et al. 2013) studied here.

Appendix B: Atlas of the H II region+PDR best-fitting model schematics

For each galaxy, we show a 2D schematic of the best-fitting single or mixed model, based on model results presented in Tables 3 and 4. Galaxies are sorted by metallicity (from highest to lowest). Color coding corresponds to density. The central stellar cluster is represented in green. The distance to the inner black shell varies with the inner radius, and the thickness of the inner black shell is proportional to the ionization parameter U, while the thickness of the H II and PDR regions are kept fixed. U is a function of the input stellar spectrum, density, and inner radius. In the case of mixed models, the angle shown by dotted lines indicates the contribution from each component (i.e., parameter f_c in Tables 3 and 4).

Fig. B.1. Best-fitting single or mixed model, based on model results presented in Tables 3 and 4.

Fig. B.1. continued.

Fig. B.1. continued.

Appendix C: Additional correlations for the galaxy sample

In Fig. C.1, we show correlations of the fraction of [C II] emission from the ionized gas deduced from the models as a function of the model density and ionization parameter, as well as the correlation between the galaxy metallicity and the model ionization parameter.

Fig. C.1. Top: correlation of the fraction of [C II] emission from the ionized gas with the model parameters: the electron density (left panel) and the ionization parameter (right panel). Values of component #2 are plotted in case of mixed models. Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner. Bottom: distribution of ionization parameter values from the models with the observed galaxy metallicity. A horizontal line links the two U-values in the case of mixed models. In all panels, the data points are color coded by metallicity.

Appendix D: Results of the H II region+PDR models for individual galaxies

We show the results of the H II region+PDR modeling of individual galaxies with three main types of plots: (a) Comparison of the observed and predicted line intensities; (b) Contribution of each gas component to the total line prediction in case of mixed models; (c) Probability density functions (PDFs) of the fitted parameters. Only the results of the galaxy He 2-10 are shown in the main text of the paper (Fig. 4).

Fig. D.1. Results for Haro 2. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.2. Results for Haro 3. See below for caption description.

Fig. D.3. Results for Haro 11. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.4. Results for He 2-10. See below for caption description.

Fig. D.5. Results for HS 0052+2536. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.6. Results for HS 0822+3542. See below for caption description.

Fig. D.7. Results for HS 1222+3741. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.8. Results for HS 1304+3529. See below for caption description.

Fig. D.9. Results for HS 1319+3224. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.10. Results for HS 1330+3651. See below for caption description.

Fig. D.11. Results for HS 1442+4250. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.12. Results for II Zw 40. See below for caption description.

Fig. D.13. Results for I Zw 18. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.14. Results for Mrk 153. See below for caption description.

Fig. D.15. Results for Mrk 209. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.16. Results for Mrk 930. See below for caption description.

Fig. D.17. Results for Mrk 1089. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.18. Results for Mrk 1450. See below for caption description.

Fig. D.19. Results for NGC 1140. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.20. Results for NGC 1569. See below for caption description.

Fig. D.21. Results for NGC 1705. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.22. Results for NGC 5253. See below for caption description.

Fig. D.23. Results for NGC 625. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.24. Results for Pox 186. See below for caption description.

Fig. D.25. Results for SBS 0335-052. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.26. Results for SBS 1159+545. See below for caption description.

Fig. D.27. Results for SBS 1211+540. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.28. Results for SBS 1249+493. See below for caption description.

Fig. D.29. Results for SBS 1415+437. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.30. Results for SBS 1533+574. See below for caption description.

Fig. D.31. Results for Tol 1214-277. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.32. Results for UM 448. See below for caption description.

Fig. D.33. Results for UM 461. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.34. Results for VII Zw 403. See below for caption description.

Fig. D.35. Results for UGC 4483. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.36. Results for UM 133. See below for caption description.

Fig. D.37. Results for HS 0017+1055. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

Fig. D.38. Results for NGC 4214-c. See below for caption description.

Fig. D.39. Results for NGC 4214-s. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

All Tables

All Figures

	Fig. 1. Spectra from Starburst99 used as input for Cloudy. Spectra are normalized and shown for each metallicity bin value of our grid. There are more hard photons at low metallicity.
In the text

Fig. 2. Adopted hydrogen density profile in the default grid of models as a function of visual extinction, shown for Zbin = 0.5 (solid line) and Zbin = 0.1 (dot-dash line). We also show the density profile for the constant pressure test performed in Sect. 5.4 (gray solid line). Model parameters are set to nH = 102 cm−3 and log U = −2. The phase transitions from H+-dominated to H I-dominated and from H I-dominated to H2-dominated are indicated with vertical lines. We note that AV of 1 mag corresponds to a hydrogen column density that is 5× larger at Zbin = 0.1 than at Zbin = 0.5.

In the text

Fig. 3. Grid results: predictions of selected line intensities and model parameters, for Zbin = 0.5 (top panels) and Zbin = 0.1 (bottom panels). Color scheme: increasing initial density from 100.5 cm−3 (light green) to 103.5 cm−3 (dark purple). Symbol size: large symbols correspond to small inner radius (large U) and small symbols to large inner radius (small U). Dotted lines connect models with the same inner radius, solid lines connect models with the same density. All observations are overplotted in orange and galaxies in the metallicity bin plotted are highlighted in red. We note that the error bars are large because they include abundance uncertainties (assumed to be 50% of the flux). The PDR covering factor is set to unity. Top right panel: quantitative indication on how the C/O abundance ratio affects the grid predictions from Zbin = 0.5 to Zbin = 0.1. Other effects due to metallicity (e.g., stellar spectra, column density, etc.) are discussed in the text.

In the text

Fig. 4. Model results, shown only for the example galaxy He 2-10 here and for all galaxies in Appendix D. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

Fig. 5. 2D representation of the 1D best-fitting (mixed) model for the example galaxy He 2-10. Schematics for the other galaxies of the sample are shown in Appendix B. The model parameters are the same as those reported in Table 4. Color coding corresponds to density. The central stellar cluster is represented in green. The distance to the inner black shell is proportional to the inner radius, and the thickness of the inner black shell is proportional to the ionization parameter U, while the thickness of the H II and PDR regions are kept fixed. U is a function of the input stellar spectrum, density, and inner radius. The angle shown by dotted lines indicates the contribution from each component (i.e., parameter fc in Table 4).

In the text

Fig. 6. Sensitivity of the model results to the fitted lines. We show histograms of the reduced χ2 () and best-fitting parameters (nH, log U, covPDR) for the entire sample of galaxies. The results from the default grid of models and the default set of lines are in black, and the results are in color if the set of fitted lines is changed, that is, when one spectral line is excluded or when the H2 lines included or when we consider a small, fixed set of lines to fit ([Ne III]15, [O III]88, [C II]157, and LTIR). The vertical dotted lines indicate the mean values. Each histogram is offset vertically by a value of 20. See Sect. 5.1 for details.

In the text

	Fig. 7. Cloudy input radiation fields tested. We show the input (dashed) and output (solid) spectral energy distributions for: a continuous star-formation scenario of 10 Myr (with an X-ray component), a Starburst99 instantaneous burst of 3.5 Myr, a Popstar burst of 3.5 Myr, a Magphys fit to the optical-IR photometry. All have model parameters: nH = 102 cm−3, rin = 1020.4 cm, and Zbin = 0.5.
In the text

Fig. 8. Model predictions for the [O I]145/[O I]63 luminosity ratio as a function of visual extinction in the cloud. The default model (solid curve) has parameters: nH = 103 cm−3, G0 = 103, Zbin = 0.5. Increasing the cosmic ray rate by a factor of 10 (dashed curve) or adding a soft X-ray component to the input radiation field (dotted curve) both increase the PDR heating and predicted [O I] line ratios. The gray shading indicates the observed range of ratios.

In the text

Fig. 9. Fraction of [C II] emission coming from the ionized gas phase. We show predictions for each DGS galaxy based on: (1) the best single or mixed models reported in Tables 3 and 4 (crosses), and (2) the best single or mixed models with one model having a fixed low density of 30 cm−3 (circles); predictions yielding (all models are bad) are not plotted. Red symbols indicate that [N II]122 is detected and orange symbols indicate that both [N II]122 and [N II]205 are detected. For the three galaxies where [N II]205 is detected, we also show the predicted values from simple theoretical calculations, as done in Croxall et al. (2017). We overlaid the data (diamonds) and fit from Croxall et al. (2017) as well as the prediction from PDR models from Kaufman et al. (2006) in gray. Fits to the DGS data and to the DGS (free density) and Croxall et al. (2017) data are shown in blue.

In the text

Fig. 10. Correlation of the fraction of [C II] emission from the ionized gas with: metallicity, specific star-formation rate, and line ratios for which the strongest trends are found. Values for the y-axis are taken from the models in boldface in Tables 3 and 4 (triangles for single models, filled circles for mixed models). Color coding corresponds to metallicity, except for the last two panels which are color coded by electron density and ionization parameter (of component #2 in case of mixed models). Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner (ρ). In the last three panels, star symbols are the ratios of observed [C II]157 or [O III]88 intensity divided by the predicted [N II]122 intensity. The correlation coefficients ρall include those star symbols.

In the text

Fig. 11. Distribution of model parameters. Top row: ionization parameter U and density nH for the best-fitting single models (circles) or mixed models (component #1: squares and component #2: triangles). For mixed models, the black lines link the two model components and the symbol size is proportional to the scaling factor fc of each component. Bottom row: same for the radiation field G0 impinging the PDR and the average PDR density nH, PDR. Color coding corresponds to Z (left panels) and sSFR (right panels). Galaxies with no sSFR value are displayed with open symbols.

In the text

Fig. 12. Correlations of the PDR covering factor with metallicity (top left panel), specific star-formation rate (top right panel), LTIR/LFUV ratio (bottom left panel), and [O III]88/[C II] ratio (bottom right panel). Galaxies are color coded by Z. Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner. Lower-Z galaxies tend to have smaller covering factors but we do not find strong trends with observed line ratios. The apparent floor of 0.2 in the PDR covering factor is due to the chosen steps for this parameter in the grid.

In the text

	Fig. B.1. Best-fitting single or mixed model, based on model results presented in Tables 3 and 4.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

Fig. C.1. Top: correlation of the fraction of [C II] emission from the ionized gas with the model parameters: the electron density (left panel) and the ionization parameter (right panel). Values of component #2 are plotted in case of mixed models. Spearman’s rank correlation coefficient and significance of its deviation from zero in parenthesis are indicated in the top-right corner. Bottom: distribution of ionization parameter values from the models with the observed galaxy metallicity. A horizontal line links the two U-values in the case of mixed models. In all panels, the data points are color coded by metallicity.

In the text

Fig. D.1. Results for Haro 2. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.2. Results for Haro 3. See below for caption description.
In the text

Fig. D.3. Results for Haro 11. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.4. Results for He 2-10. See below for caption description.
In the text

Fig. D.5. Results for HS 0052+2536. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.6. Results for HS 0822+3542. See below for caption description.
In the text

Fig. D.7. Results for HS 1222+3741. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.8. Results for HS 1304+3529. See below for caption description.
In the text

Fig. D.9. Results for HS 1319+3224. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.10. Results for HS 1330+3651. See below for caption description.
In the text

Fig. D.11. Results for HS 1442+4250. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.12. Results for II Zw 40. See below for caption description.
In the text

Fig. D.13. Results for I Zw 18. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.14. Results for Mrk 153. See below for caption description.
In the text

Fig. D.15. Results for Mrk 209. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.16. Results for Mrk 930. See below for caption description.
In the text

Fig. D.17. Results for Mrk 1089. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.18. Results for Mrk 1450. See below for caption description.
In the text

Fig. D.19. Results for NGC 1140. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.20. Results for NGC 1569. See below for caption description.
In the text

Fig. D.21. Results for NGC 1705. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.22. Results for NGC 5253. See below for caption description.
In the text

Fig. D.23. Results for NGC 625. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.24. Results for Pox 186. See below for caption description.
In the text

Fig. D.25. Results for SBS 0335-052. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.26. Results for SBS 1159+545. See below for caption description.
In the text

Fig. D.27. Results for SBS 1211+540. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.28. Results for SBS 1249+493. See below for caption description.
In the text

Fig. D.29. Results for SBS 1415+437. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.30. Results for SBS 1533+574. See below for caption description.
In the text

Fig. D.31. Results for Tol 1214-277. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.32. Results for UM 448. See below for caption description.
In the text

Fig. D.33. Results for UM 461. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.34. Results for VII Zw 403. See below for caption description.
In the text

Fig. D.35. Results for UGC 4483. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.36. Results for UM 133. See below for caption description.
In the text

Fig. D.37. Results for HS 0017+1055. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text

	Fig. D.38. Results for NGC 4214-c. See below for caption description.
In the text

Fig. D.39. Results for NGC 4214-s. Panel a: comparison of the observed (losange) and predicted intensities for our best-fitting single (blue bar) and mixed (red bar) H II region+PDR models. Lines are sorted by decreasing energy. The fitted lines have labels in boldface. Panel b: for mixed models: respective contributions from component #1 and component #2 to the line prediction. For PDR lines ([C II], [O I], and [Si II]), only the contribution from the ionized gas is shown; the PDR contribution is one minus the sum of the contributions from the plot. Contributions are expected to be dominated by one or the other component if their model parameters are noticeably different. Panel c: probability density functions of the model parameters (nH, U, covPDR/scaling factor) for single (left panel) and mixed (right panel) models. The scaling factor corresponds to the proportion in which we combine two models in the mixed model case. It is equal to unity otherwise. It is shown with dashed lines in the same panels as the PDR covering factor. Vertical dotted lines show values of the best-fitting model. Black asterisks show the range and step of values of the grid parameters.

In the text