The Herschel Dwarf Galaxy Survey: II. Physical conditions, origin of [CII] emission, and porosity of the multiphase low-metallicity ISM

The sensitive infrared telescopes, Spitzer and Herschel, have been used to target low-metallicity star-forming galaxies, allowing us to investigate the properties of their interstellar medium (ISM) in unprecedented detail. Interpretation of the observations in physical terms relies on careful modeling of those properties. We have employed a multiphase approach to model the ISM phases (HII region and photodissociation region) with the spectral synthesis code Cloudy. Our goal is to characterize the physical conditions (gas densities, radiation fields, etc.) in the ISM of the galaxies from the Herschel Dwarf Galaxy Survey. We are particularly interested in correlations between those physical conditions and metallicity or star-formation rate. Other key issues we have addressed are the contribution of different ISM phases to the total line emission, especially of the [CII]157um line, and the characterization of the porosity of the ISM. We find that the lower-metallicity galaxies of our sample tend to have higher ionization parameters and galaxies with higher specific star-formation rates have higher gas densities. The [CII] emission arises mainly from PDRs and the contribution from the ionized gas phases is small, typically less than 30% of the observed emission. We also find correlation - though with scatter - between metallicity and both the PDR covering factor and the fraction of [CII] from the ionized gas. Overall, the low metal abundances appear to be driving most of the changes in the ISM structure and conditions of these galaxies, and not the high specific star-formation rates. These results demonstrate in a quantitative way the increase of ISM porosity at low metallicity. Such porosity may be typical of galaxies in the young Universe.


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 Tables of observed line fluxes and predicted line luminosities from the models are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc. u-strasbg.fr/viz-bin/qcat?J/A+A/626/A23 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 metalrich 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

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

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

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

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 . The dust-to-gas ratio (DGR) is 5 × 10 −3 for Z bin = 1 and decreases linearly with metallicity. The abundance of PAHs was further scaled with Z bin to the power A23, page 3 of 44 A&A 626, A23 (2019) Table 2. Input parameters and abundances for the model grids.

Radiation field
Starburst99 continuous SF of 10 Myr, . . . Geneva stellar track, scaled to Z, . . . Kroupa IMF, . . . log(L/L ) = 9 Blackbody T = 5 × 10 6 K (X-ray), . . . log(L/L ) = 5.2 Cosmic rays background 0.5 log CMB Notes. Other metals (not shown here) have abundances that are based on the ISM abundance pattern, and are then scaled to the metallicity (Z bin ) of the grid. We note that the luminosity of the different radiation field components are the reference values used for the model computation, but models are scaled afterwards to the observed luminosities of each source.
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 −16 s −1 (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 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 Z bin = 0.5 (solid line) and Z bin = 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 n H = 10 2 cm −3 and log U = −2. The phase transitions from H + -dominated to H i-dominated and from H i-dominated to H 2 -dominated are indicated with vertical lines. We note that A V of 1 mag corresponds to a hydrogen column density that is 5× larger at Z bin = 0.1 than at Z bin = 0.5. at a fixed N(H) value instead of a fixed A V 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 0 , and H 0 and H 2 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.
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 −3 , and the inner radius such that we cover an ionization parameter (U) from 10 −4 to 10 −1 . 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.

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

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

Goodness of fit
We searched for the best-fitting models by minimizing the χ 2 defined as 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 χ 2 ) 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 χ 2 is given by χ 2 ν = χ 2 /ν where ν is the difference between the number of constraints and the number of free parameters.

Reliability of the derived physical parameters
The χ 2 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 A&A 626, A23 (2019) Notes. This table reports results of galaxies for which single models are preferred, and Table 4 reports results for galaxies where mixed models are preferred. The following notes apply to both tables. (a) N c is the number of observational constraints (number of detected lines plus 1 for the luminosity constraint). (b) n H,HII corresponds to the initial density. (c) n H,PDR corresponds to the density of the PDR at the H i-H 2 transition (i.e., is the contribution of the first (second) model to the mixed model. We note that the values of the covering factor are those of the initial grid (with full coverage of the H ii region). For 2-component models, cov PDR needs to be further multiplied by f c to obtain the effective covering factor. Galaxies are sorted by decreasing metallicity. For one-component models, uncertainties on the free parameters are given in upper and lower scripts and are calculated as the average ± the standard deviation of models within the 1σ confidence level (see Sect. 3.2 for details). The number of free parameters n free is three for the one-component model, and five for the two-component model. The degrees of freedom are given by ν = N c − n free . The reduced χ 2 is given by χ 2 ν = χ 2 /ν.
average and standard deviation of the values of the free parameters taken by models with χ 2 values within a certain range of χ 2 min . For instance, for a number of degrees of freedom of four, and a 1σ confidence level (probability of 68.3%), this range is ∆χ 2 = χ 2 − χ 2 min = 4.72 (Press et al. 1992).

Results
We present general grid results in Sect. 4.1 and results from fitting the grid models to the observations in the following subsections. 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 0 produced when varying the ionization parameter and the density, which are input parameters of the model grid. We remind the reader that G 0 is the intensity of the ultraviolet radiation field at the PDR front, given in units of the background Habing (1968) value, 1.6 × 10 −3 erg cm −2 s −1 , as defined by Tielens & Hollenbach (1985). For a given density, U and G 0 scale linearly (except at the high-U, low-n H end due to geometric dilution of the field; see e.g., Cormier Notes. See caption of Table 3  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] 88 /[C ii] 157 line ratio and between 0.3 and 5 for the [C ii] 157 /[O i] 63 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.

Grid results
We remind the reader that normal, star-forming, metal-rich galaxies have similar [ (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).

Best-fitting single H ii+PDR models
Parameters of the best-fitting single models as well as their corresponding χ 2 and reduced χ 2 (χ 2 ν ) 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 bestfitting 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 −3 , with median at 10 2.0 cm −3 . The medians are very close to the values log(U) = −2.5 and n H = 10 2.0 cm −3 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 A&A 626, A23 (2019) Grid results: predictions of selected line intensities and model parameters, for Z bin = 0.5 (top panels) and Z bin = 0.1 (bottom panels). Color scheme: increasing initial density from 10 0.5 cm −3 (light green) to 10 3.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 Z bin = 0.5 to Z bin = 0.1. Other effects due to metallicity (e.g., stellar spectra, column density, etc.) are discussed in the text. About half of the galaxies (16/39) have minimum χ 2 ν below two and their observed lines are satisfyingly reproduced by a single model. Inspection of the results indicates that galaxies with larger χ 2 ν 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].

Motivation
It is possible that the minimum χ 2 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. 2012Dimaratos 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 1 ≥ U 2 ), 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 .
For galaxies with fewer than six lines detected (for which a mixed model would be under-constrained) or with χ 2 ν 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).   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 f c in Table 4).

Results for mixed models
In order to reproduce all observed line intensities in galaxies where a single model fails (χ 2 ν > 2), 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 bestfitting mixed models as well as their corresponding χ 2 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] 157 /[O i] 63 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).
A23, page 9 of 44 A&A 626, A23 (2019) For most galaxies for which single models do not provide a good fit, the mixed models perform better (lower χ 2 ν,min ). All galaxies for which mixed models are preferred have χ 2 ν,min ≤ 2. For several galaxies though, mixed models do not improve the results over the single models but provide similar χ 2 ν,min (II Zw 40, Mrk 209, Pox 186, SBS 0335−052, SBS 1159+545, SBS 1415+437) or worse χ 2 ν,min (Tol 1214−277). For two of them (SBS 1159+545, Tol 1214−277), χ 2 ν,min 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 χ 2 ν,min close to unity). For the galaxies where a single model was already giving a good fit, the mixed models provide lower χ 2 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 −3 and 10 1.0 cm −3 versus 10 3.0 cm −3 and 10 1.5 cm −3 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.

Mid-infrared lines
We have not used the silicon, iron, and H 2 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. Notes. Columns 1-2: line label and number of galaxies in which the line is detected. Columns 3-4: number of galaxies in which the ratio of predicted over observed intensity is low, ok, and high, where the threshold is set to a factor of two. We report numbers for single and mixed models.
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 2 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 2 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 0 > 10 2 (Kaufman et al. 2006). When fitting for the H 2 lines, the densities required for the PDR are indeed larger (typically ≥10 2 cm −3 at the PDR front) but χ 2 ν increases a lot (models performing worse for some H ii lines, More H 2 detections are required to really investigate the presence and importance of shocks.

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 A23, page 10 of 44 D. Cormier et al.: Modeling the multiphase ISM of star-forming dwarf galaxies 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 −3 (value often within the acceptable range) instead of 3 cm −3 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.

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] 88 luminosity. We find that the observed Hα/[O iii] 88 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).

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 −3 ), 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   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 −3 ). 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 0 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.

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

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] 157 line which can arise from more diffuse gas than [O i], the [O i] 63 line which can be optically thick, the [Ne iii] 15 line which can be problematic to model from population synthesis models (see Sect. 5.2), and the [S iv] 10 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 2 lines in the fit and fitting a small, fixed set of lines ([C ii] 157 , [O iii] 88 , [Ne iii] 15 , 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 χ 2 ν 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, χ 2 ν seems to improve quite significantly when excluding [S iv]. For mixed models, the distributions of χ 2 ν remain similar and the mean χ 2 ν 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] 63 or [Ne iii] 15 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 −3 ), 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 A&A 626, A23 (2019) Fig. 6. Sensitivity of the model results to the fitted lines. We show histograms of the reduced χ 2 (χ 2 ν ) and best-fitting parameters (n H , log U, cov PDR ) 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 H 2 lines included or when we consider a small, fixed set of lines to fit ( 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 2 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 χ 2 ν values on average.

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.

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 A23, page 12 of 44 D. Cormier et al.: Modeling the multiphase ISM of star-forming dwarf galaxies 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 0 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 (χ 2 ν,min 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.  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.

Heating of the PDR and [O i] line ratio
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 4 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 6 K (to peak around 1 keV) and luminosity in the soft X-ray range 10 −4 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 H + 3 (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 −16 s −1 ; 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).

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 −1 in the default grid of models, as in Kaufman et al. (1999), and here, we explored the cases v turb = 0 km s −1 and v turb = 5 km s −1 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 −1 to 5 km s −1 , 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 2 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.

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 χ 2 ν,min 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.

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. 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] 63 ratios (lower, right panel of Fig. 3. We find that the χ 2 ν,min 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 −4 might be too high since Lebouteiller et al. (2017) adopt a value of 7 × 10 −6 . 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 A23, page 14 of 44 D. Cormier et al.: Modeling the multiphase ISM of star-forming dwarf galaxies 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) 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.
does not have a noticeable impact on the results, but additional PDR observations (such as MIR H 2 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.

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 −3 , as well as shocks (Osterbrock & Ferland 2006). [N ii] 122 is the best line that we have in hand but, in principle, [N ii] 205 is better suited (see Table 1  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. 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 −3 . 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 −3 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] 205 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 andCroxall 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 lowdensity component (n H = 30 cm −3 ) can be unreliable because they sometimes yield much higher χ 2 than when n H is let free. Those predictions are not plotted in Fig. 9 if they yield χ 2 ν > 3. 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 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 −3 , 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 A&A 626, A23 (2019) 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.
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 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] 122 . 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). 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 −3 , 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 starforming 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 twophase 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] 157 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 starforming spirals and IR-luminous galaxies. This offset might be Fig. 11. Distribution of model parameters. Top row: ionization parameter U and density n H 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 f c of each component. Bottom row: same for the radiation field G 0 impinging the PDR and the average PDR density n H,PDR . Color coding corresponds to Z (left panels) and sSFR (right panels). Galaxies with no sSFR value are displayed with open symbols. 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 ).

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.
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 0 and n H,PDR . We also find a tendency for lower-Z galaxies to have higher G 0 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 0 and n, with G 0 ∝ 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 0 /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 0 /n H,PDR values than the median, indicating more compact/intense starforming regions, consistent with the increase of G 0 /n H,PDR at high IR surface densities (though for more moderate values of G 0 /n H,PDR and Σ IR in our sample) reported by Díaz-Santos et al. (2017). Analytically, given our adopted density law, G 0 /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 0 /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.

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] 122 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 Although the dwarf galaxies are offset in the observed parameter space of those ratios compared to metalrich galaxies , 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.
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 radiationbounded (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. 2011Leitet et al. , 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 highredshift Universe.

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 −3 , corresponding to PDR densities roughly ten times larger (given the adopted density profile). The highest U values are higher than in metalrich 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  (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 2 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.