© A. Lambert-Huyghe et al. 2022

1 Introduction

The interstellar medium (ISM) harbors many processes that play key roles in the evolution of galaxies. It is composed of material that has been accreted from the intergalactic reservoir as well as the heavy elements (metals heavier than hydrogen or helium) produced in stars and returned to the ISM at the end of their lives to feed the next generation of stars. In return, this star formation process and the interaction between stars and the ISM evolve with time, as the maturing stars are cycling their newly made elements to the ISM, enriching the latter ISM reservoir with metals from each generation. Despite recent progress in high-z surveys (e.g., Maiolino et al. 2015; Inoue et al. 2016; Knudsen et al. 2016; Pentericci et al. 2016; Matthee et al. 2017, 2019; Carniani et al. 2018, 2020; Hashimoto et al. 2019b,a,c; Hashimoto 2020; Le Fèvre et al. 2020; Izumi et al. 2021), following the process of metal enrichment as galaxies evolve through time, from the earliest galaxies to the present day, has not been an easy observational endeavor with the current telescope sensitivities. Thus, many questions are still open regarding low metallicity ISM conditions and the formation of stars at the earliest epochs. A robust understanding of the origins of the tracers of star-forming regions and insight into their diagnostic capabilities are some of the necessary steps to guide the interpretation of high-z ISM and star-forming conditions.

While not completely mimicking the earliest galaxies, local Universe dwarf galaxies are often used as laboratories to understand the physical properties of the gas and dust and their interplay with star formation in early Universe environments. The proximity of nearby dwarf galaxies makes many different tracers accessible, some of which are too faint to be detected in more distant galaxies. Moreover, the diverse collection of local Universe dwarf galaxies exhibit a wide range of metallicities and star formation rates (SFR).

Local dwarf galaxies were the focus of large Herschel and Spitzer surveys (e.g., The Dwarf Galaxy Survey, DGS; Madden et al. 2006). Studies on both resolved and integrated-galaxy scales have highlighted some distinctively unique observational signatures of star-forming low-metallicity dwarf galaxies. A non-linear relation of the dust-to-gas mass (D/G) with metallicity is observed, with extremely low dust masses observed for the lowest metallicity galaxies (Z ≤ 0.1 Z_⊙; Herrera-Camus et al. 2012; Fisher et al. 2014; Rémy-Ruyer et al. 2015; Galliano et al. 2018, 2021; Cigan et al. 2021). Furthermore, the hard radiation fields in star-forming dwarf galaxies, along with their lower dust abundance, result in extended ionized gas regions prominent on global galaxy scales (Hunter et al. 2011; Cormier et al. 2015, 2019). The consequence is often a largely photodissociated molecular phase, existing in clumps which are difficult to observe with the usual molecular gas tracer, CO (1-0) (Cormier et al. 2014; Hunt et al. 2015; Accurso et al. 2017b), beckoning the presence of a CO-dark molecular gas phase (Grenier et al. 2005; Röllig et al. 2006; Wolfire et al. 2010; Glover & Clark 2012; Bolatto et al. 2013; Accurso et al. 2017a; Madden et al. 2020). Other emission lines, however, such as the far-infrared [C ii]λl58 µm line, are strikingly enhanced on global scales in dwarf galaxies (e.g., Cormier et al. 2015, 2019; Cigan et al. 2016; Lebouteiller et al. 2017; Jameson et al. 2018), making the [C ii]λl58 µm line a potential tool for tracing star formation activity (Malhotra et al. 2001; Papadopoulos et al. 2007; Pineda et al. 2014; De Looze et al. 2014; Olsen et al. 2015; Herrera-Camus et al. 2015, Herrera-Camus et al. 2018; Carniani et al. 2018; Matthee et al. 2019; Izumi et al. 2021; Fujimoto et al. 2021) and a tracer of the total H₂ in galaxies, near and far (Poglitsch et al. 1995; Wolfire et al. 2010; Pineda et al. 2013; Nordon & Sternberg 2016; Fahrion et al. 2017; Accurso et al. 2017b; Zanella et al. 2018; Madden et al. 2020; Schaerer et al. 2020; Tacconi et al. 2020).

In this paper, we explore the use of different infrared (IR) tracers, in particular the [O iii]λ88 µm line together with the midinfrared (MIR) continuum at 24 µm, as a new method to understand the evolution of the warm dust and gas phases in galaxies. This new tool can also be used to prepare future observations, both at high redshift and in local galaxies. For the rest of the paper we make use of the following shorthand notation: [O iii] = [O iii]λ88 µm; [C ii] = [C ii]λl58 µm; and [O i]λ63 µm = [O i].

The most luminous farinfrared (FIR) line observed in low-metallicity star-forming galaxies is the [O iii] line, even brighter than the classically used [C ii] line, which is the brightest FIR line in metal-rich galaxies (Stacey et al. 1991; Malhotra et al. 2001; Luhman et al. 2003; Brauher et al. 2008; Smith et al. 2017; Díaz-Santos et al. 2017; Croxall et al. 2017; Sutter et al. 2019, 2021). The energy needed to create O⁺⁺ (35 eV) makes it a direct tracer of the ionized gas around young stars and thus an important coolant of the ionized medium and a SFR tracer (e.g., De Looze et al. 2014; Herrera-Camus et al. 2015). The predominance of the [O iii] line over full dwarf galaxy scales (Cormier et al. 2015, 2019) is an indication of an extensive filling factor of ionized gas in low metallicity environments, and can explain the potential for large scale photodissociation of the CO molecule in dwarf galaxies. The [O iii] line has been observed in surveys of local galaxies of a wide range of metallicities with Herschel (e.g., DGS; KINGFISH: Kennicutt et al. 2011; Madden et al. 2013; Lamarche et al. 2017; Pereira-Santaella et al. 2017; Fernández-Ontiveros et al. 2016) and is becoming a line frequently observed in high-z galaxies (e.g., Ferkinhoff et al. 2010; Matthee et al. 2017; Carniani et al. 2017; Tamura et al. 2019; Hashimoto et al. 2019b,c; Hashimoto 2020).

Often associated with warm dust in the ionized phases in galaxies, the MIR continuum at 24 µm, has been mapped with Spitzer/MIPS in many galaxies (e.g., Draine & Li 2007; Bendo et al. 2007, 2008, 2012; Croxall et al. 2012; Vutisalchavakul & Evans 2013; Kendall et al. 2015; Boquien et al. 2016) including full maps of the nearby Magellanic Clouds (e.g., Meixner et al. 2006; Gordon et al. 2011). Assumed to be a tracer of dust in H ii regions, it is often used as a SFR tracer (see e.g., Kennicutt & Evans 2012, and references therein), but it can also originate in dust heated by older stellar populations (e.g., Leroy et al. 2012).

Therefore, most of the [O iii] and 24 µm emission is associated with recent star formation in galaxies, suggesting that there should be an observational relationship between them. This paper investigates the link between [O iii] and 24 µm for the first time and explores the origins of the contributions as well as the deviations of the relationship between these tracers. We investigate the potential to predict the [O iii] line from the 24µm band, with spatially resolved and integrated data (also to be used to predict 24 µm from [O iii] observations). We also investigate the accuracy of the prediction and what drives deviations from the global behavior.

In order to understand what those relations on global scales, it is important to obtain a precise understanding of the properties and relation to the local environment of gas and dust phases at a resolved scale, and the consequences on the observed tracers. Only in this way can one be confident in employing certain diagnostic tracers for unresolved high-z studies as well. The closest low metallicity galaxies of the Milky way are the Large and Small Magellanic Clouds (respectively LMC and SMC). Their proximity (50kpc for the LMC, Walker 2012, 61 kpc for the SMC, Hilditch et al. 2005) gives us access to details of the properties of their star-forming regions, and allows us to disentangle the different phases of the ISM under low metallicity environments (Z_LMC = 1/2 Z_⊙, Pagel 1997, 2003, Z_SMC = 1/5 Z_⊙, Russell & Dopita 1992; Cioni et al. 2000). Thus, the Magellanic Clouds are the ideal laboratories to investigate the relation between [O iii] and 24 µm in detail.

This paper is structured as follows. In Sect. 2, we present the Herschel and Spitzer data used in the study. A direct relation between [O iii] and the MIPS 24 µm band is studied in Sect. 3. In Sect. 4, we present the models prepared for our analysis, and the specific technique we use to disentangle the emission from H ii regions and PDRs. The dependency with different physical parameters is analyzed in Sect. 5, where we also present a new method developed for this study allowing investigation of the effects of various physical parameters, including the proportion of PDRs needed in the models to reproduce the observations. This method can be seen as a step toward understanding the 3D geometry of the ISM. The method and the analyses are discussed in Sect. 6, and we conclude our work in Sect. 7.

2 Data

The regions studied here are driven by the available Herschel/PACS spectroscopy (Pilbratt et al. 2010; Poglitsch et al. 2010) from the DGS survey, where maps of some FIR fine structure lines were observed for numerous star-forming regions of the LMC and SMC (Table 1; Cormier et al. 2015). To this we added Spitzer 24µm photometry (Werner et al. 2004; Houck et al. 2004; Rieke et al. 2004) from the SAGE project (Meixner et al. 2006) which mapped the LMC and SMC completely.

2.1 Spitzer data

We used the Spitzer IRAC and MIPS maps available on the IRSA Infrared Science Archive¹ and selected subregions covering our Herschel spectroscopy maps. Two regions present some saturation in the MIPS 24 µm band: 30 Doradus and N160. For 30 Doradus, the MIPS 24 µm data were replaced with Spitzer/IRS observations, integrated over the transmission profile of the MIPS band (Chevance et al. 2016). As N160 was not observed with the IRS, the part of the map affected by saturation was simply masked.

Visual inspection showed that the 24 µm maps are strongly affected by point source emission, often corresponding to stellar sources. As we are interested in the ISM properties, we removed the point sources from the MIPS 24 µm map. Positions of the point sources were first selected from dedicated catalogs (Seale et al. 2012; Gruendl & Chu 2009; Jones et al. 2017); among these, sources to be removed were selected as those also detected as point sources in the IRAC 8 µm band. Verification was made by checking the presence of point source emission at the selected positions in the IRAC maps. A 2D point spread function (PSF) was fitted on the point sources, and removed from the 24 µm maps. Then, the background was fitted, assumed to be flat on tiles with a size depending on the PSF and the intensity of the point-sources. The aim was to avoid degeneracy in the extraction of the sources. For two regions we did not make this correction. For N160 the saturated part of the map was already masked, and all the point sources lie within this masked area. The N66 region contains a large number of point-sources, which lead to numerous artifacts when removing them, making the map impossible to use. Thus, we decided to keep the MIPS 24 µm map of N66 without correction for the point-sources. Keeping the 24 µm emission corresponding to point sources in this region does not significantly affect the global results in our study, but increases the uncertainty on the examined relationships.

The limiting spatial resolution for the correlations investigated in this study is due to the MIPS 70 µm observations, which have a beam full width at half maximum (FWHM) of 18″. Thus, the maps of all the tracers were convolved to 18 with a Gaussian kernel and resampled to a grid with pixels of 12″ width in order to avoid oversampling of the beam.

2.2 Herschel data

The Herschel/PACS spectroscopic maps of LMC and SMC regions were initially presented in Cormier et al. (2015), along with other spectroscopic observations of the DGS sample. The processing of the data is fully described in Cormier et al. (2015). In this study, data were reprocessed with the latest version of HIPE (v.15, Ott 2010) and PACSman (v3.63, Lebouteiller et al. 2012). The maps were also convolved with a Gaussian kernel to match the MIPS 70 µm resolution and similarly resampled to match the 24 µm maps.

For further comparison of our work with other dwarf galaxies, we used DGS (Madden et al. 2013), which targets 48 dwarf galaxies of the local Universe, covering a large range of conditions, including metallicity and SFR. Indeed, the metallicity of the DGS sample ranges between 1/50 Z_⊙ and 1/2 Z_⊙, with the LMC and SMC being at values of 1/2 Z_⊙ and 1/5 Z_⊙, respectively. The SFR of the DGS varies widely, from 10⁻³ M_⊙ yr⁻¹ to 100 M_⊙ yr⁻¹, with values of ~0.2 M_⊙ yr⁻¹ and ~0.04 M_⊙ yr⁻¹ for the LMC and SMC respectively. The distances range between 50kpc and 60kpc for the LMC and SMC, respectively to nearly 200 Mpc for the furthest target of the DGS sample (see Meixner et al. 2013; Madden et al. 2013, and references therein for further details).

3 Relation between [OIII] and MIPS 24µm

We first investigated if there is a relation between surface flux density of [O iii] and 24 µm. In general, this comparison could be done for any line unambiguously directly tracing H ii regions. We used [O iii] here since this is a line frequently observed at both high and low redshift, especially with Herschel and ALMA surveys. We also looked into the accuracy of the relationship between these tracers and, in cases where there may be deviations from the global behavior, into what drives the differences.

3.1 Spatially resolved relation

We took advantage of the excellent spatial resolution of Herschel/PACS and Spitzer/MIPS together with the proximity of the Magellanic Clouds to investigate the relation between the [O iii] and 24 µm at ~5 pc scale. Figure 1 shows the correlation between the observed 24 µm continuum and the [O iii] line at a spatial resolution of 12″ pixels in all of our sources for which the [O iii] line was detected at more than 5σ. We see that the different sources together exhibit a wide range of values, spanning around 3 orders of magnitude in 24 µm and 4 orders of magnitude in [O iii]. We found a clear linear relation between those two tracers in log-log space. This correlation is fitted with a weighted χ² method, leading to the following relation: $$\log \left({\Sigma }_{\left[\text{OIII}\right]88}\right)=0.879\left(\pm 0.012\right)\times \log \left({\Sigma }_{24}\right)-1.446\left(\pm 0.065\right),$$(1)

where the surface brightness, Σ, is expressed in Wm⁻² sr⁻¹. We see that the spread varies significantly with the flux. It is tight (around 0.5 dex) for the brightest regions, especially 30 Doradus, N160, and the brightest regions in N159, but it can reach almost 2 dex deviations in other regions. This spread is partly due to differences between the local physical properties of the regions (e.g., between N44 and N11C), and intrinsic variations within each region (e.g., N159).

If we consider each region separately (Fig. 2), we see that the slopes of most of the individual regions are close to the global slope which seems to be driven by 30Dor. All except a few positions in N44, N11I and N159 fall within the 95% confidence level found from the global relation, with the deviant pixels often being low brightness pixels. In N11C (Fig. 2), for example, the slope begins to steepen when the 24 µm surface brightness is less than ~ 10⁻⁶ W m⁻² sr⁻¹. These particular low-surface brightness pixels lie on the outskirts of the star-forming peak in N11C (see Fig. 3). Thus, we see that the [O iii]-24um relation depends on the surface brightness. Moreover, we see that the regions are distributed along the relation, some of them showing strong emission of both [O iii] and 24 µm (e.g., 30 Doradus, N159), and some showing weak emission of both tracers (e.g., N11I, N180). This variation could be linked to the age of the ionizing source. A younger cluster will produce more ionizing photons, which increases the emission of the [O iii] line, and it will be surrounded by a more compact shell of dust, resulting in an increase in the 24 µm emission. The two tracers are probably not affected similarly, thus the relation may be affected by the age of the star-formation burst. On the other hand, 24 µm emission can be powered by an older population of stars, especially in regions that have sustained multiple episodes of star formation (Leroy et al. 2012). Thus, part of the spread of the data around the relation could be due to 24 µm emission arising from older stars, especially in more quiescent regions such as N11I or N180.

N66 is the only region in the SMC in this study. Both the [O iii] and 24 µm may be affected by the lower metallicity, which can globally decrease the amount of dust and [O iii] in the gas phase simply owing to the lower metal abundance. On the other hand, as the dust abundance is diminished, the ISM may be subjected to an overall harder radiation field, which may increase the 24 µm emission by shifting the peak of the SED to shorter wavelengths (see e.g., Galliano et al. 2018, and references therein). Despite those effects, the [O iii]-24 µm relation in N66 does not differ much from that of the LMC regions. The slope of the region appears flatter (Fig. 2) than that of the global relation pointing to the possibility that [O iii] and 24 µm may not be as closely linked in the SMC as in the LMC regions. For example, the [O iii] in N66 may originate from a more extended, diffuse component that is not completely associated with compact star-forming regions from which the 24 µm originates. More studies on other regions in the SMC would be necessary to conclude on this point.

We note that the fits to the [O iii] and 24 µm correlations were made using convolved and resampled maps, and the MIPS 24 µm maps were treated to remove contamination from point sources (cf. Sect. 2.1). Any of these steps could possibly influence the fitting results. In Appendix B, we detail the effects of different fitting methods and different data treatment and point-source removal methods: (a) effect of different fitting methods: comparison of a simple, nonweighted χ², and a Monte Carlo simulation to determine uncertainties and the weighted χ² fit; (b) effect of resolution and sampling of the maps on the fits; (c) tests of the method for point source subtraction and corrections. We conclude from the different tests that the fit remains robust in these different cases.

The spatially resolved data in Fig. 1 show a clear, linear relation between [O iii] and 24 µm in log-log space, which can be used, for example, to predict [O iii] line emission from maps of the 24 µm band, or inversely. The predictive power of this relation and the effect of the spread around the relation are studied in detail in Sect. 3.3. However, the LMC and SMC are the closest neighbors of the Milky Way, and most of the local Universe galaxies will not provide such spatially detailed observations. We thus explored if the relation also holds for unresolved regions, to see if the correlation can be applied in a more general case and to more distant galaxies.

Fig. 1 Relation between [O iii] and 24 µm in all of our SMC and LMC regions, only including pixels above 5σ in the [O iii]. The black solid line represents the fitted relation, linear in log-log space, and the dashed lines and gray filling represent the 95% confidence interval of the fit. The colored symbols refer to the individual regions labeled in the legend to the right of the plot.

Fig. 2 Relation between surface brightness of [O iii] and 24 µm for all of the pixels of our individual star-forming regions detected at more than 5σ. The black solid line represent the fitted linear relation in loglog space for the total sample (Fig. 1), and the dashed lines and gray filling represents the 95% confidence interval associated with the total sample.

Fig. 3 Spatial distribution of the low emission pixels in N11C. Left panel: [O iii] map for N11C; right panel: MIPS 24 µm map for N11C. Both maps are convolved to 12″ and resampled to 12″ pixels. The black crosses correspond to pixels that have low surface brightness (on the left of the bulk of the regions in Figs. 1 and 2, corresponding to values of 24 µm lower than 106 W m−2 sr−1) showing a different behavior from that of the global relation.

3.2 Integrated relation and comparison with the DGS

In order to mimic unresolved observations, we used the spatial maps of the regions at their original resolution and sampling, without correction from point sources. The surface brightness of the complete regions were then calculated over the zones where the observations of the two tracers, [O iii] and 24 µm, overlapped. To take into account the different distances of the Magellanic Clouds and the DGS galaxies, we converted the flux values into luminosities, normalized to solar luminosity, L_⊙. The uncertainties of the integrated surface brightness values were calculated via a quadratic sum, and converted in the same way as the integrated luminosity values. We see in Fig. 4 that the integrated luminosities of the different regions of the Magellanic Clouds also show a linear trend in log-log space. Using a weighted χ² method on those star-forming regions, we obtained the relation of Eq. (2) for the integrated Magellanic Cloud regions: $$\log \left(L_{\left[\text{OIII}\right]88}\right)=0.948\left(\pm 0.086\right)\times \log \left(L_{24}\right)-0.622\left(\pm 0.360\right),$$(2)

with L_[OIII]88 being the luminosity of [O iii] and L₂₄ being the luminosity of the 24 µm continuum, both expressed in L_⊙. We see that the slope obtained for the integrated luminosities of the Magellanic Clouds is consistent with that obtained for surface brightness values (Eq. (1)) within the uncertainties. We also investigate the relation with the surface luminosity.

In Fig. 4, we also compare the [O iii] and 24 µm luminosities of the integrated Magellanic Cloud regions with those of the integrated DGS targets. They are normalized by the solar luminosity (upper panel) or by the surface (lower panel). We examined three different fits, one on our Magellanic Cloud integrated sample, one on the DGS sample and one on all the data (Table 2). We see that the different fits give close values, agreeing with each other considering the uncertainties. However, in the upper panel, some of the DGS targets lie significantly below the fit on the Magellanic Clouds regions, and the fit on the DGS sample has a somewhat shallower slope. As discussed in Sect. 6, this may indicate a small enhancement of 24 µm emission compared to [O iii] on the scale of total galaxies. This effect is more evident in the lower panel, where the luminosity per surface area is shown. Indeed, [O iii] arises from energetic ionized media heated by young O and B stars, while 24 µm emission can be powered by older stellar populations as well, more prominently when integrating over full galaxy scales. Leroy et al. (2012) studied the proportion of 24 µm not associated with recent star formation in a sample of 30 local disk galaxies. They found a median value of 19% of the total 24 µm emission that is linked to older stellar populations, and possibly up to 40% in some parts of the galaxies. The DGS sample we use here is a collection of star-forming dwarf galaxies, known to have a bursty star formation history, in contrast to Leroy et al. (2012) sample of disk galaxies for which a more continuous star formation is invoked. The effect of prominent recent SF and lower dust abundance in this dwarf galaxy sample (e.g., Rémy-Ruyer et al. 2014; Galliano et al. 2018) leads to a more porous ISM for hard UV photons globally (see e.g., Galliano et al. 2003, 2005; Madden et al. 2006, 2020; Cormier et al. 2015, 2019; Ramambason et al. 2022). Thus, the proportion of 24 µm emission not associated with recent star formation is probably lower than that found in disk galaxies. Some of 24 µm not arising from the neighborhood of young star-forming regions may thus be linked to the presence of non-star-forming regions within the beam. However, this difference is small, suggesting that most of the 24 µm emission is still associated with star-forming regions, even when considering the full emission of a galaxy (e.g., Contursi et al. 2002; Calzetti et al. 2005; Alonso-Herrero et al. 2006; Pérez-González et al. 2006; Prescott et al. 2007). Our small maps of Magellanic Cloud regions are centered on the brightest zones of star formation, largely filled by ionized gas (see e.g., Kurt & Dufour 1998; Nazé et al. 2001; Barbá et al. 2003; Lebouteiller et al. 2012; Chevance et al. 2016; Gordon et al. 2017; McLeod et al. 2019). If we could observe more extended regions within the Magellanic Clouds, similar to full galaxy scales, we might collect more PDR gas within our beam and possibly some 24µm emission from PDRs as well as older stellar populations, which would then further enhance the 24 µm emission compared to [O iii] emission. Despite this caveat, this relation holds over more than three orders of magnitude for both [O iii] and 24 µm, and is still valid for numerous galaxies at different distances from the Milky Way, giving our calibration a potentially powerful predictive potential, that we demonstrate in the following section.

Fig. 4 Upper panel: relation between the integrated [O iii] and 24 µm luminosity, in units of solar luminosity (L⊙). Lower panel: integrated luminosity per surface area of [O iii] and 24 µm (Wkpc−2). The DGS sample corresponds to the gray circles, and the integrated regions from the Magellanic Cloud sample are coded with colors, as indicated in the inserted legend in the upper panel. The fit on the Magellanic Clouds alone is indicated as MCs, the fit on DGS sample is indicated as DGS, and the fit on the whole data is indicated as DGS+MCs.

3.3 Predictive power of the [O iii] – 24 µm calibration

We used our two calibrations from Eqs. (1) and (2) to explore the predictive potential of the observed correlation between 24 µm and [O iii], and we compared them to the resolved and unresolved observations. We calculated the uncertainties of the predicted values with a Monte Carlo simulation, using the uncertainties on the fitted parameters. Those uncertainties are compared to the difference between the observed and predicted values of [O iii].

We first compared the [O iii] predictions we derive from the observed 24 µm using surface brightness, Eq. (1), to the observed [O iii] maps of the Magellanic Cloud regions. As examples, we show in Figs. 5 and 6 the comparison for the 30Dor and N11I regions, respectively. The predicted and observed maps for the other regions are shown in the Appendix C.

The predicted and observed [O iii] emission agree within a factor of 3 for 90% of the full data set, with a global discrepancy not higher than 1 dex, and for all of the regions but two, the predicted and observed [O iii] agree within this factor of 3 for more than 80% of the data set (see Table 3). We also point out that the intensity of [O iii] does not seem to strongly affect the goodness of the prediction down to 3 × 10⁻⁸ Wm⁻² sr⁻¹. Only N11I, which is the faintest region in both 24 µm and [O iii], exhibits [O iii] emission under this value. This is also one of the two regions for which only ~50% of the data shows agreement between observations and predictions within a factor of 3, and shows a large discrepancy between the observed spatial distributions of [O iii] and 24 µm (see Fig. 6). This region seems to be more quiescent than the rest of the sample, and [O iii] and 24 µm emission do not seem to be tightly linked, leading to different spatial distributions.

Other regions with significant differences between predictions and observations are N11B and N11C (see Figs. 7, C.2 and C.3). These regions lie slightly above the global relation in Fig. 2, although they have similar slopes. Table 3 also shows that those regions exhibit the highest values of [O iii]_obs/[O iii]_pred. Likewise, we see from Fig. 2 that N44 and N159 have numerous pixels with observed 24 µm and [O iii] values lying below the global relation. In particular, only 53% of the N44 observed emission agrees within a factor of 3 with the predictions. The exploration of the physical conditions in the different regions in Sect. 5 shows that the physical parameters are significantly different in this region compared to the rest of the sample. The offsets of some regions from the observed global correlation (e.g., N11B) leads to under or overestimation of the [O iii] emission by a factor of a few over most maapped areas. It can be due to differences in physical conditions within the individual regions, which is investigated in Sect. 5.2.2. Nevertheless, our correlation allows us to make predictions of [O iii] from the observed 24 µm which are accurate up to a factor 3.

We also investigated the predictions for the integrated values of our sample and the DGS values, based on the luminosities. Here again, predicted and observed [O iii] luminosities agree within a factor of 3 for the Magellanic Cloud sample and for all the DGS targets (Fig. 8) except four. The discrepancies between predictions and observations are similar to those found for the resolved data, with a factor of a few for most of the regions, and a couple of extreme DGS galaxies for which the disagreement can be up to almost 1 dex. To summarize, the calibrations we derived here allowed to predict [O iii] within a factor 3 for most of the data, in both surface brightness and luminosity, for resolved and unresolved data. They also hold for galaxies with different physical conditions and distances from the Milky Way see Sect. 2.2 for more details. Thus, the calibrations we derive here may be useful to prepare observations at both low and high redshift, from the scale of spatially resolved, individual star-forming regions to integrated galaxies. We point out that this calibration can be used in both ways, to also predict 24 µm emission, adding potentially useful continuum information for high-redshift galaxies where the [O iii] is detected.

However, we also showed that there is a non-negligible spread in some cases, around the correlations derived from Figs. 1 and 4, that can be due to differences in the overall physical conditions between the individual regions of our Magellanic Clouds sample, as well as locally within an individual region. In order to understand the spread around this relation, we built models to study the effect of some physical parameters.

Fig. 5 Predicted 30 Dor map of [O iii] (top panel) given the observed 24 µm and Eq. (1). For comparison, the actual observed map of [O iii] (middle panel), and the ratio of observed over predicted emission (lower panel). The two maps of [O iii] are represented on the same scale and are displayed in Wm−2 sr−1.

Fig. 6 Similar to Fig. 5, but for N11I.

Fig. 7 Spatially resolved ratios of observed ([O iii]obs) over predicted ([O iii]pred) emission of [O iii]. The value of 1, corresponding to good agreement of the predictions and the observations, is indicated with a horizontal dashed line. The horizontal gray lines indicate differences of a factor 3 from agreement.

Fig. 8 Predicted [O iii] luminosity, from the observed 24 µm for integrated luminosities of the Magellanic Clouds sample and DGS galaxies compared to the observed luminosity values. The horizontal dashed line indicates when predicted and observed values are equal. The horizontal gray lines indicate differences of a factor 3 from agreement.

4 Modeling strategy

We used the spectral synthesis code, Cloudy (version 17.01, Ferland et al. 2017) to interpret the observations. Cloudy is a 1D spherical photodissociation and photoionization code that selfconsistently solves the thermal, chemical, and radiation equilibrium of a multiphase medium exposed to a radiation source and produces a complete spectral energy distribution of the gas and dust. Input parameters include the spectrum of the central source of the radiation field and the total initial hydrogen density at the illuminated face of the cloud, n₀. We chose the stellar population model PopStar (Mollá et al. 2009), which simulates the SED of a full stellar population based on a Chabrier initial mass function (Chabrier 2003) at different ages, for a single burst of star formation. We did not test the assumption of a continuous star formation history because our observations are centered on small regions presenting recent star formation, thus it is likely that the ionizing clusters in our maps have similar ages. This is also consistent with previous studies on the stellar content of our regions (e.g., Bica et al. 1996; Heydari-Malayeri et al. 2000; Caulet et al. 2008; Carlson et al. 2012; Ochsendorf et al. 2017). In our models, we supposed that the ISM is in pressure equilibrium, thus we set constant total pressure through the cloud. We also ran a set of models with density kept constant through the cloud, but it was not possible to reproduce all the observed tracers. We also noticed that characteristic tracers from PDRs were better reproduced by models with a higher density at the inner radius than those reproducing the characteristic tracers of H ii regions, which is consistent with constant pressure models, and is consistent with previous works comparing the different assumptions (e.g., Tielens & Hollenbach 1985; Hollenbach & Tielens 1999; Bron et al. 2018; Cormier et al. 2019).

We then varied two parameters: n₀, the density at the inner radius, and the ionization parameter (U), which describes the effect of the ionizing source on the gas and dust. U is the dimensionless parameter relating the number of hydrogen-ionizing photons emitted by the central source, Q(H⁰), impacting the inner surface of the cloud of radius, r, the speed of light, c and $$U=\frac{Q\left({\text{H}}^0\right)}{4\pi r^2{n}_{0}c}.$$(3)

n₀ is varied in log(cm⁻³) from −1 to 3 with steps of 0.5, and U is varied in log space from −4 to −1 with steps of 1. The age of the simulated cluster is varied between 1 and 10 Myr, with steps of 0.2 in log space. The default dust abundance and polycyclic aromatic hydrocarbons (PAH) properties in Cloudy are those from the Milky Way. We applied a multiplicative factor to the abundances to match the SMC and LMC metallicities. In order to completely model the PDR up to the molecular region, the models were stopped at an A_v of 2.5, which ensures that the formation of the H₂ molecules in all our models is reached, which is considered as the beginning of the molecular region. We initially computed models for two metallicities, 1/2 Z_⊙ and 1/5 Z_⊙, corresponding to the conditions in the LMC and SMC respectively However, the difference between the two sets of models in the outputs used in our work were not significant. Thus, we kept one set of models with the LMC metallicity, as most of our regions are from the LMC. Table 4 summarizes the range of the parameters and the steps of each parameter used to create the grids.

Figure 9 is an example of the models we created with Cloudy, showing the evolution of different physical parameters and the cumulative emission of the tracers we use in our following study with the depth of the cloud (presented with A_V values). The limit between the H ii region and the neutral gas is fixed at the point where the electronic density drops below half its maximum value (A_V ~ 3 × 10⁻² in our example). The density is low in the H ii region, and increases in the PDR (up to a value of ~ log(n₀) = 5.6 in our example, to be compared to the initial value at the illuminated edge of the cloud of log(n) = 2.5).

Concerning the lines and bands studied here, two immediate features are noteworthy. Firstly, we noticed that all of the [C ii] emission arises from the PDR in our models, not the H ii region. This is consistent with the global studies of low-metallicity dwarf galaxies where most of the [C ii] arises from PDRs (e.g., Cormier et al. 2019). Secondly, we found in our models that a significant fraction of the 24 µm continuum emission arises from PDRs, while it is often implicitly assumed to trace warm dust in H ii regions (see e.g., Draine & Li 2007; Croxall et al. 2012, for this assumption). These two points will be important for the comparison with data in Sect. 5. The ISM observed within our beam is composed of a mixture of different phases with different physical conditions (diffuse ionized regions, PDRs, denser H ii regions, etc.), and cannot always be modeled by a single Cloudy model (i.e., 1D model simultaneously solving one H ii region plus one PDR).

We saved the cumulative emission as a function of the depth inside the cloud to examine the variation of some tracers with the cloud depth. Then, we separated the emission into two parts: that arising from the H ii region, until the ionization front, and that arising from the PDR², beyond the ionization front and up to A_v = 2.5. We considered the total emission from the H ii region as a pure H ii region model and the total emission from the PDR as a pure PDR model. We then mixed the two components to mimic a physical mixture of ionized region and PDR along the line of sight for each pixel in the map, using Eq. (4): $$\{\begin{array}{ccc}{I}_{\text{mixed\hspace{0.17em}model}}\hfill & =\hfill & c_{\text{H\hspace{0.17em}II}}{I}_{\text{H\hspace{0.17em}II}}+c_{\text{PDR}}{I}_{\text{PDR}}\hfill \\ \hfill & =\hfill & c_{\text{H\hspace{0.17em}II}}{I}_{\text{H\hspace{0.17em}II}}+\left(c_{\text{PDR}}-c_{\text{H\hspace{0.17em}II}}\right)I_{\text{PDR}}\hfill \\ c_{\text{H\hspace{0.17em}II}}+c_{\text{PDR}}\hfill & =\hfill & 1,\hfill \end{array}$$(4)

where I_{H ii} is the intensity of the tracer in the H ii region only in W m⁻² sr⁻¹; I_PDR is the intensity of the tracer in the PDR only and I_{mixed model} is the total intensity of a tracer for the mixed model along the line of sight. c_{H ii} and c_PDR are the scaling factors of the H ii region and PDR respectively. As we keep their sum equal to 1 in our strategy, we solve for c_PDR in the following work. The multiplicative factor c_{H ii} for the full model (I_{H ii+PDR}) is similar to a luminosity scaling, but we preferred to use ratios between lines and continuum bands.

With this method, we illustrated different physical configurations for the distribution of H ii regions and PDRs by changing c_PDR, the proportion of PDR (and by construct, H ii regions) in Eq. (4). We note that this mixture is made only from models having the same initial n₀ and U. We used this simple approach to avoid degeneracy between the parameters, as a first step. A more complex use of our model will be left for future research.

Some of the cases that can be described by Eq. (4) are illustrated in Fig. 10. For example, when there is 0% of PDR, the model corrspond to a naked H ii region. Similarly, a proportion of 100% PDR corresponds to a model of an isolated PDR. We underline that in reality, the model of an isolated PDR component corresponds to an observation with only PDR emission along the line of sight, but the PDR can be a part of a larger structure that also contains H ii regions not seen along this line of sight. A proportion of 50% c_PDR corresponds to a full Cloudy model (1D connected H ii region and PDR).

The extreme cases of 0% and 100% c_PDR are easy to represent, and all of the other cases are based on the interpretation of Eqs. (A.1) to (A.3) detailed in Appendix A. The c_PDR we used here does not correspond to the amount of PDRs in the medium, but to the amount of PDR emission from the model we used. A c_PDR greater than 50% corresponds to a covering factor of the PDR on the H ii region greater than 1, with more PDRs along the line of sight than a simple model, or in other words there is more than one cloud along the line of sight, with potentially different intercepts between the cloud and the line of sight. Thus, it can be a way of interpreting the distribution of matter in the third dimension.

In reality, the distribution of the matter in the ISM can be much more complex. Around a young star, neither the H ii region nor the PDR may fully cover the stellar source. The H ii region can be reduced to a clump of ionized gas near the star, and the covering factor of the PDR can be lower than unity, leading to a configuration with part of the cloud surrounding the star being matter-bounded H ii regions, where ionizing photons escape (no or little PDR beyond the ionized gas), and part of it being radiation bounded H ii regions where no photons with energies greater than 13.6 eV escape (with a full PDR beyond the ionized gas). This means that if along the line of sight there is a thick PDR covering only a part of the H ii region, we model it as a thinner PDR which is fully covering the H ii region. Higher spatial resolution observations, of course, would help to constrain a more detailed modeling scheme.

This method allowed us to constrain the input physical parameters of our Cloudy models: the age of the burst (Sect. 5.1), n₀ and U (Sect. 5.2). Based on those parameters, it was possible to determine c_PDR for each pixel in the map. To constrain our mixed model, we used a ratio of observed PDR and H ii region tracers, such as ([C ii]+[O i])/[O iii], to determine c_PDR needed in the model to reproduce the observed ratio, using Eq. (5) as follows. $$\begin{array}{ccc}{R}_{\text{obs}}& =& \frac{I_{\text{CII,obs}}\text{+}{I}_{\text{OI,obs}}}{I_{\text{OIII,obs}}}\\ & =& \frac{\left(c_{\text{HII}}{I}_{\text{CII,HII}}+c_{\text{PDR}}{I}_{\text{CII,PDR}}\right)+\left(c_{\text{HII}}{I}_{\text{OI,HII}}+c_{\text{PDR}}{I}_{\text{OI,PDR}}\right)}{c_{\text{HII}}{I}_{\text{OII,HII}}+c_{\text{PDR}}{I}_{\text{OIII,PDR}}}\\ c_{PDR}& =& {\left(1-\frac{I_{\text{CII,PDR}}+I_{\text{OI,PDR}}-R_{\text{obs}}{I}_{\text{OIII,PDR}}}{I_{\text{CII,HII}}+I_{\text{OI,HII}}-R_{obs}\times I_{\text{OIII,HII}}}\right)}^{-1},\end{array}$$(5)

I_CII,_obs the observed line intensity of [C ii]; I_{OI obs} is the observed line intensity of [O i] and I_OIII,obs is the observed line intensity of [O iii], all along the same line of sight. We used Eq. (4) in the ratio R_obs for each tracer in order to compare the models to the observations (Sect. 5), where I_{CII, HII} is the intensity of the [C ii] line in the H ii region model, with similar definitions for the other lines, and I_CII,PDR is the intensity of the [C ii] line in the PDR model, with similar definitions for the other lines. Equation (5) provides c_PDR which gives the final formula for the PDR proportion along each sight line.

Fig. 9 Example from our Cloudy models. The cumulative emission of different bands and emission lines are presented (upperpanel), together with some physical parameters (density, electronic density and temperature, lower panel) as a function of the depth in the cloud, measured with the parameter, Ay. The emission of all of the tracers was normalized to a maximum of 1 to allow a direct comparison. The age of the central ionizing source is 4 Myr, and the initial values of physical parameters are log(n0)= 2.5 and log(U) = −3.

Fig. 10 Sketch of different PDR proportions (cPDR) in our mixed models. The yellow regions represent H ii regions, and PDRs are represented in gray. The blue arrow indicates the direction of increasing cPDR. The black arrows show the lines of sight through the represented clouds (each line of sight representing a given pixel in our maps). As per Eq. (4), the 0% cPDR case corresponds to no PDR in the model, hence fully ionized gas; 50% cPDR corresponds to a single Cloudy model, with a connected H ii region and PDR (Eq. (A.2)); and 100% cPDR represents a full PDR model, without an H ii region contribution. We note that 50% of PDR does not correspond to a medium with half of the volume being filled by PDR, but corresponds to a mix that can be modeled by a Cloudy model with a full PDR layer.

5 Disentangling dust emission from the lonized phase and PDRs

We saw in Sect. 3 that the relation between [O iii] and 24µm emission shows a non-negligible dispersion. Based on models presented in Sect. 4, we investigated the possibility that part of the 24 µm may not arise from the same phase as [O iii]. On the other hand, some regions are globally offset from the global relation, thus different physical conditions between the regions could also explain part of the spread around the relation. In order to examine those two potential effects and constrain our models, we used the [O iii] line as a tracer of the H ii regions, and the [C ii]λ158 µm line and [O i]λ63 µm line as tracers of the PDR. In higher metallicity galaxies, the [C ii] line can arise from both H ii region and PDRs, but in low-metallicity dwarf galaxies, it mostly originates from PDRs (20–30% at most coming from H ii regions, see e.g., Cormier et al. 2019; Tarantino et al. 2021 and it is not greater than 5% in the SMC according to Jameson et al. 2018). We thus considered in this paper that the observed [C ii] emission comes from PDRs. The [O i] line has a clearer origin, and is used to trace PDRs with high UV illumination and/or high density (see e.g., Tielens & Hollenbach 1985; Tielens 2005; Osterbrock & Ferland 2006).

We first investigated the physical conditions of the different regions, and especially the parameters we used to create our grid of models: age, n₀ and U. Then, we used Eq. (5) to estimate in each pixel the contribution from PDRs and the ionized phase to the 24 µm emission. We used ratios of the tracers to be able to compare the different regions without requiring a calibration of the total luminosity, especially the ratio [O iii]/24 µm, to explore the contribution of each phase to the 24 µm emission, and the ratio R_PDR/HII = ([C ii]+[O i])/[O iii], that we used as a proxy for the ratio of PDR over H ii region emission.

5.1 Age of the burst

There are many studies on the star-forming regions on the Magellanic Clouds, and the age of recent star formation is frequently studied, mainly through the presence of O and B stars, or Young Stellar Objects (YSOs; see e.g., Heydari-Malayeri & Testor 1986; Heydari-Malayeri et al. 2002; Heydari-Malayeri & Selier 2010; Bolatto et al. 2000; Nazé et al. 2001; Tsamis et al. 2003; Gruendl & Chu 2009; Lebouteiller et al. 2012; Seale et al. 2012; McLeod et al. 2019; Okada et al. 2019). As we previously noted, the [O iii] line is a characteristic tracer of the massive young stellar population, while the 24 µm emission can be produced by other stellar types. Thus, the ratio [O iii]/24 µm can be used to trace the overall age of the stellar population. We sampled ages between 1 and 10 Myr, and compared the observations to the predicted line ratio [O iii]/24 µm for models that cover a range of n0, U and c_PDR (Fig. 11).

From Fig. 11, we see that the predicted [O iii]/24 µm ratio is roughly constant from 1 to ~4 Myr, and drops after this age, as the younger, more massive stars die off. Our observations show values that are consistent with the early age models. This is also consistent with several previous studies on these individual regions, proposing ages younger than 10 Myr, often found to be ages of a few Myrs (see e.g., Lucke & Hodge 1970; Oey & Massey 1995; Bica et al. 1996; Parker et al. 1996; Heydari-Malayeri et al. 2000; Caulet et al. 2008; Carlson et al. 2012; Ochsendorf et al. 2017; Bestenlehner et al. 2020). For the rest of the study, we used the model with an age of 10^6.6 (~4 Myr), noting that it is not possible with our tracers to constrain ages younger than this. We also point out that the final results of our study are not affected if younger ages are chosen.

Based on Fig. 11, we can also begin to investigate the likelihood of some values of the density and ionization parameter. The models with log(n₀) = 3 seem less consistent with the bulk of the data, and only those with a lower c_PDR agree with the data. On the other hand, models with low density (log(n₀) = 1) and low c_PDR (<50%) are also not very consistent with most of the data. A deeper investigation is done in the following sections, but the global trend from this figure shows that an intermediate density (log(n₀) ~ 2) is more consistent for the bulk of the data. From Fig. 11, it is not possible to conclude on the ionization parameter or the c_PDR value.

Fig. 11 [O iii]/24 µm for different models as a function of burst age. Each model is represented with a color for n0, and a line style for U, and the lines are following the proportion of PDR emission, cPDR, with the markers representing the cPDR of 0%, 50 % and 99%. The blue solid horizontal line corresponds to the median of the spatially resolved observations and the dashed blue lines to the median absolute deviation. The distribution of our observations is shown on the right as a histogram.

5.2 Determining the ISM conditions

We then considered the model parameters, n₀ and U. We inspected the ratio [O iii]/24 µm as a function of R_PDR/HII. A high value of R_PDR/HII indicates that the emission originating from PDRs is important. As the MIPS 70 µm band emission generally traces cooler dust than that of the 24 µm band, and is emitted mainly by PDRs in our models, we also compared the ratio [O iii]/70 µm to R_PDR/HII

5.2.1 Global results for the sample

The observed [O iii]/24 µm and [O iii]/70 µm ratios are represented as functions of R_PDR/HII in Fig. 12 with constant pressure models overlaid, for an age of 4 Myrs. We see that the modeled ratios with the PDR proportions agree with the data, with a relatively flatter trend at low PDR proportions (50% or less), and a decreasing trend for higher proportions (over 50%). The use of [O iii]/24 µm and [O iii]/70 µm together also allows us to better constrain the models, and find a value of U and n₀ agreeing with all of the observations. The models that best agree with the entire data set have log(n₀)= 2.5, log(U) = −3 and various values for c_PDR. Notice that the c_PDR fitting the data is below 50% for the most part. This means that the emission can be modeled with a single connected H ii region and PDR model and a covering factor of the PDR between 0 and 1³.

We see that the values of the physical parameters seems quite homogeneous for our sample. However, there is a non-negligible dispersion of the data, thus we investigated the physical conditions for each region individually.

Fig. 12 Representation of the data of the full sample and models with constant pressure for [O iii] 88 µm/MIPS 24 µm as a function of ([C іі]+[O i])/[O iii] (top panel), and [O iii] 88 µm/MIPS 70 µm as a function of ([C ii]+[O i])/[O iii] (lower panel). Each model is represented with a color for the initial density, n0, and a line style for U, and the lines follow the cPDR, with the markers representing the proportions of 0%, 50% and 99%. The models displayed are calculated for an age of star-forming burst of 4 Myr. The gray points are the observed values for each pixel of the regions.

5.2.2 Individual regions

The different regions have been investigated individually (see figures in Appendix C). We see for each region that the models that best fit the data are close to the best model found for the whole sample, but differences can be seen from one region to another. We summarized the best fitting models for the total sample and the individual regions in Table 5.

In addition to this best model, there are other models in our grid that agree with parts of some regions. Those models are indicated in the second column of Table 5. The possibility for multiple models to agree with part of the data can be understood as our regions are not homogeneous: there are clumps and structures that are spatially resolved in our maps, thus having different physical conditions, which correspond to different models. The best model for each region is then the most representative model of the global conditions in the region, despite possible local variations.

We see that when multiple models agree, an increase of U corresponds to a decrease of n₀ (Table 5). It is probable that the density we report here corresponds to a rather diffuse ionized gas, which is ubiquitous in our maps. Although there are dense clumps observed in almost all of the regions in our sample, the subregions with the highest densities are hidden by two effects: the dilution of dense clumps by the surrounding lower density medium, due to the resolution and sampling of our data, and the fact that we report the initial densities at the illuminated edge of the cloud, where the actual densities in our models can rise to a few 10⁵–10⁶ cm⁻³ at high A_V, corresponding to the high-density clumps observed.

N44 (Fig. 13) is the region that deviates most from the global distribution of the parameters, being the only region with the highest ionization parameter, log(U) = −2, while all of the other regions and the global value of N44 indicate log(U) = −3. N44 also has a lower n₀ value: log(n₀) = 1.5 while the other sources have log(n₀) values between 2 and 2.5. McLeod et al. (2019) found that a large bubble (45.7 pc radius) was carved out in the N44 region by the central star cluster, leading to a large zone with very low density (n ≤ 100 cm⁻³) with multiple regions of star formation triggered around this main bubble, each of them carving out the ISM resulting in different measured densities. Our PACS region is similar to the subregion of triggered star formation N44C for both the position and the size of the map, although our maps are slightly shifted toward the main bubble. Since the very low density medium inside the bubbles cover almost half of our PACS map, this morphological observation may explain the global low density we found, which leads to a higher ionization parameter for a given radius. It is a condition unique to N44, as the other regions in our sample do not show such a large excavation of the ISM. However, we also see from Table 5 and Fig. 13, that other models with lower U (log(U) = −3) and higher densities (log(n₀) = 2.5–3) can be consistent with a subregion of N44, which is in agreement with the density found by McLeod et al. (2019) for the N44C region. Here they measure an electron density of 152 ± 42 cm⁻³ with a much smaller aperture (~8″) than that of our PACS observations. Thus, best model values are consistent with globally very low density, where the other models with log(U) = −3 and log(n₀) = 2.5–3 agree well with the higher density structures of swept up material in the region.

Our best model is selected by picking the model that agrees best with the trend delineated by the observations. However, in some regions, the observations are consistent with different models in our grid. In Fig. 13, the model grids flatten at low values of R_PDR/HII with an almost constant ([C ii]+[O i])/[O iii], inconsistent with the observations which show continually rising values of R_PDR/HII. In the case of N44, we see on Fig. 13, that models with a lower U and a higher n₀ are more consistent with the observations at low values of R_PDR/HII. In this region, it can be interpreted as the edge of the bubble described in McLeod et al. (2019).

While the PDR tracers, [O i] and [C ii] are prominently emitted in all of the star-forming regions studied, most regions have quite low proportions of PDR components along the lines of sight throughout their regions (see Fig. 10 for visual representation of different PDR proportions). All but 5 regions have c_PDR,< 12% throughout their maps. The N11C region shows pixels with c_PDR at least 30% along their lines of sight and only three regions exhibit c_PDR > 50% for some of the pixels of the regions: N44, N11I and N159. Both N159 and N44 show complex morphologies, with arcs, bubbles, and clumps. They also show ongoing star formation, with associated molecular clumps. On the other hand, N11I harbors a more quiescent and extended molecular clump. The higher c_PDR values exhibited toward parts of the maps could be linked to the recent stellar formation activity and the birth cloud that may not be totally dispersed. For N11I, it could also be linked to ongoing star formation inside the quiescent molecular cloud. As most of the regions have low c_PDR, this implies that most of the pixels are largely dominated by H ii regions, and there is little neutral gas that can be observed along the lines of sight compared to the more prominent ionized gas. On the other hand, pixels where the PDR proportion is higher than 50% correspond to an excess of PDRs considering a simple model with a connected H ii region plus PDR. This may mean that we have multiple clouds along the line of sight, or that we are looking through the external part of a cloud that is spatially resolved. In all of the cases, we see that the PDR proportion varies significantly within an individual region. In the following we thus investigated the spatially resolved distribution of c_PDR, to have a better understanding of how it can be linked to physical interpretation.

Fig. 13 Similar to Fig. 12, but for the region N44.

5.3 Insight on the distribution of ionized and neutral gas

For each region, assuming the best fitting model determined in Sect. 5.2, we calculated c_PDR needed to match the observed ([C ii]+[O i])/[O iii] ratio (Eq. (5)). Those maps are displayed for 30 Doradus (Fig. 14) and N11I (Fig. 15), along with other emission maps and ratio maps for comparison. More region maps are displayed in Appendix C.

From the comparison of the [C ii], [O i] and [O iii] emission (Fig. 14; Fig. 15), we see that the distribution of the c_PDR is, perhaps surprisingly, not closely correlated spatially with that of the PDR tracers, [C ii] and [O i]. The spatial distribution of the two PDR tracers is quite similar while they peak somewhat offset from the [O iii] emission. The lines of sight for the lowest c_PDR values are often found in close proximity to the peaks of the H ii region tracers. On the other hand, significant [O iii] emission can be found toward the [C ii] and [O i] peaks, thus a substantial amount of ionized gas is also present in those pixels and mixed with the PDR tracers. We concluded that strong emission of PDR tracers in one line of sight does not mean that PDRs are dominant along the line of sight.

We investigated how the physical conditions in the PDRs, in terms of the characteristic properties of the gas traced by [C ii] and [O i] influence the c_PDR. The cooling of the gas associated with dense and strongly-illuminated PDRs is normally dominated by the [O i] emission line (e.g., Tielens & Hollenbach 1985). On the other hand, [C ii] is the dominant coolant in PDRs that may be less dense, or subjected to lower irradiation. Projection effects could complicate this simple view; it is difficult to access the 3D distribution of matter. This problem arises in particular when trying to determine the density which is a key parameter for a study using [O i]. In that vein, Chevance et al. (2016) applied PDR models to FIR lines in 30 Doradus, including [C ii] and [O i] and continuum emission, to constrain the density of the different clouds and investigate the 3D matter distribution. The addition of [O i] leads to an unexpected distribution of matter around the R136 cluster at the center of the region, showing the strong benefits of including that line. Velocity-resolved observations of the lines can also bring valuable information on the ISM structure (e.g., Okada et al. 2015, 2019; Requena-Torres et al. 2016; Fahrion et al. 2017; Tarantino et al. 2021), helping to disentangle the presence of multiple clumps inside the ISM, to compare the different phases, and to investigate the possibility of self-absorption.

Okada et al. (2019) found that [O i] 63 µm shows probable self-absorption features in 30 Doradus and N159 in the LMC, perhaps accounting for ~ 10 – 20% of the total [O i] emission toward these regions, in some cases, higher than the PACS calibration uncertainties (which is around 10%, see Cormier et al. 2015). If the [O i] self-absorption is as high as 20% in our maps, it would slightly increase the ([C ii]+[O i])/[O iii] ratio. In our analysis, it may increase the derived c_PDR value, but this effect would be small or insignificant and the global results of our study would not be strongly changed.

We compared c_PDR maps to the ratio [C ii]/([C ii]+[O i]), which we call, R_CII/PDR (Figs. 14, 15 and C.1 to C.10). We see (Fig. 16) that R_CII/PDR is correlated with c_PDR. When R_CII/PDR is a low value, c_PDR is also low, and H ii regions dominate (see e.g., Fig. 14). When the ratio is high, c_PDR is higher and tends toward PDRs dominating (see e.g., Fig. 15 or C.5). The [C ii] line is the major tracer in PDRs that are less dense or less irradiated than those in which the [O i] line is the dominant tracer, thus it would be tempting to interpret [O i]-dominated PDRs as regions that are closer to the illumination source. This interpretation of R_CII/PDR would be consistent with the c_PDR value we derive: when R_CII/PDR is low, it would indicate that the PDR is close to the ionizing source, which is consistent with the low c_PDR value and that the H ii region emission dominates along the line of sight. Further from the sources, we see less ionized gas along the line of sight. c_PDR increases, and the cooling of PDR gas becomes dominated by [C ii] emission, leading to an increase of R_CII/PDR.

These results show that our approach is useful in investigating the spatial distribution of gas in a 3D sense. It may provide a rough idea of the proximity of the cloud to the ionizing sources, by combining the [O iii], [C ii] and [O i] lines with the MIPS bands at 24 µm and 70 µm.

Fig. 14 Distribution of [C ii] and [O i] emission (upper left and right panels respectively) in 30 Doradus; the ratio of ([C ii]+[O i])/[O iii] (middle left panel), the distribution of cPDR (middle right panel) with the colorbar indicating the percentage of PDR along the line of sight and the ratio [C ii]/([C ii]+[O i]) (lower panel). The contours of [O iii] are overlaid on all the panels for comparison, with the color intensity ranging from dark gray for lower emission to light gray for higher [O iii] emission.

Fig. 15 Similar to Fig. 14, but for the region N11I.

5.4 Physical parameters driving the dispersion of the relation between [O iii] and 24 µm

We now wish to know how the physical conditions we are studying affect the relation between [O iii] and 24 µm, and especially their predictive power. As we saw, the solutions for ionization parameters and densities are similar for all of the regions, except for N44. Thus, these parameters do not seem to have a significant effect on the relation and its predictive power.

Figure 17 shows that pixels with c_PDR of 50% or higher have over-predicted values of [O iii], often higher than the observation by a factor of 3 or more. On the other hand, the pixels that show under-predicted [O iii] values also exhibit some of the lowest PDR proportions found (a few percent or less). Both effects could be explained by the fact that the relation was fitted on the whole sample, including regions where the PDR emission dominates the line of sight. This behavior is not surprising. Indeed, pixels with a higher c_PDR are those for which [O iii] and 24 µm emission are less tightly associated, as there is a larger part of 24 µm arising from PDRs. This leads to enhanced values of 24 µm compared to the [O iii] emission, and thus to over-prediction of the [O iii] emission when using the calibration presented here. For pixels with low c_PDR values that show under-predicted [O iii], this can be partly linked to the use of data with high c_PDR in the fit. The pixels with high c_PDR values can influence the global fit, which is thus less accurate for the low c_PDR pixels. However, not all of the pixels with low c_PDR value show an underestimation of [O iii]. Comparing Fig. 7 and Fig. 17, the pixels with low c_PDR values and under-predicted [O iii] emission mostly come from the regions N11B and N11C. From Table 5, we see that only models with intermediate density (log(n₀) = 2−2.5) are relevant for those two regions, where models with lower density (log(n₀) ≤ 1.5) are appropriate for small parts in the other regions. Thus, differences between the physical conditions of the regions can have a significant effect on the prediction. All of those effects make the predictions less accurate for pixels with extreme PDR proportions.

We also investigated the effect of PDR proportion on the direct relation between [O iii] and 24 µm, shown in Fig. 18, where we notice some spread of the data relating to the PDR proportion. The pixels with similar cPDR values seem to follow a shallower trend than the global relation, with an offset depending on the value of c_PDR. To test this point, we fit the relation between [O iii] and 24 µm for four different bins of PDR proportion with a weighted;χ². The results from the fitting are presented in Table 6. The data and the relation fitted for the different bins are shown in Fig. 19.

Table 6 shows that fits to the data of the bins of lower c_PDR have a shallower slope than the global fit of Fig. 1 (Sect. 3.1). While the slopes for the three lowest bins are consistent with each other within the uncertainties, they are offset from each other, and not consistent with the global relation fitted earlier (see Eq. (1) in Sect. 3.1). On the other hand, the relation for the bin of highest c_PDR values (log(c_PDR) = 1.5–2.0) has a slope closer to, but somewhat steeper than that fitted on the full sample. This suggests that the fit on the complete sample is strongly affected by the pixels with high c_PDR, despite their low numbers.

While the calibration is greatly improved by this separation based on the c_PDR, we recommend using the initial calibration carried out on the global sample for general predictive use. Indeed, the calibration for low cPDR may provide a better prediction for the sight lines dominated by the H ii region emission, but it would overestimate more often for regions dominated by PDR emission. Thus, using the calibrations based on c_PDR bins would require some prior knowledge of the distribution between the phases of the ISM.

Fig. 16 Relation between the PDR proportion and the ratio [C ii]/([C ii]+[O i]), for all of the pixels in all of the regions studied here, with the same region color code as in Fig. 1. The 1% and 10% cPDR values are indicated by the two dashed horizontal lines. The vertical gray line separates regions where cooling is dominated by [C ii] (right) and regions where the cooling is dominated by [O i] (left).

Fig. 17 Ratio of [O iii]obs/[O iii]pred represented as a function of [O iii]obs, similar to Fig. 7, except that the color code corresponds to the log(cPDR).

Fig. 18 Relation between [O iii] and MIPS 24 µm in all of the pixels in all of our star-forming regions for those pixels above 5σ in the [O iii] line. This figure is similar to Fig. 1, but the points are color-coded with the value of log(cPDR) here.

6 Discussion and caveats

This new method estimating the distribution of the ISM phases along sight lines in the studied regions, is based on a few tracers only: [O iii] to trace the ionized gas, [C ii] and [O i] to trace the PDR gas, and MIPS 24 µm and 70 µm tracing the warm and cooler dust emission. We first show that 24 µm can have a significant contribution from PDR components, not only from the H ii regions as is normally assumed in the literature. The possible 24 µm emission from PDRs cannot be neglected. Then, we determined some properties of the star-forming regions, such as U and the n₀ at the illuminated edge of the cloud, which is important to constrain the age of the burst of star formation. Our method also gives access to a new parameter, c_PDR, allowing us to investigate the distribution of gas between H ii regions and PDRs along the line of sight in a new and original way.

It is always difficult to investigate the distribution of matter, whether it is spatial distribution or contribution from the different phases of the ISM. A number of studies have been carried out during the past decades to try to disentangle the geometry and distribution effects. These include, but are not limited to, multisector models with different densities, ionization parameters or radiation fields (Cormier et al. 2015); the combination of radiation-bounded and matter-bounded models (Cormier et al. 2015; Polles et al. 2019); the addition of other sources of excitation beyond optical/UV stellar light, such as X-rays (Lebouteiller et al. 2017); as well as the use of a PDR covering factor parameter applied to a standard Cloudy simulation with H ii regions and PDRs (Cormier et al. 2019). Other works attempted to make 3D simulations of turbulent star-forming clouds with a chemical network and self-consistent physical calculations, sometimes by coupling different codes (see e.g., Röllig et al. 2007; Glover & Mac Low 2011; Glover et al. 2015; Glover & Clark 2016; Offner et al. 2014; Accurso et al. 2017a; Ramos Padilla et al. 2021), or investigated the geometry and properties of the ISM with velocity-resolved studies (see e.g., Okada et al. 2015, 2019; Requena-Torres et al. 2016; Bisbas et al. 2019; Fahrion et al. 2017; Tarantino et al. 2021). Despite all of these efforts, there are still degeneracies related to the conditions and geometry of star-forming regions. Detailed observations with increased spatial and spectroscopic resolution indeed reveal very high levels of complexity in galactic ISM regions. It helps to understand the global picture when the different phases are studied as an ensemble. However, the tracers used can raise different issues, like absorption (optical/UV, but also self-absorption in some conditions, such as the [O i] line discussed above, when viewed through the Milky Way disk, or edge-on galaxies), and extinction effects. Those issues are not significant for our data, but should be investigated before using the methods developed. The possibility that a single tracer can arise from multiple phases, under some galactic environments, such as [C ii] potentially originating in both PDRs and H ii regions, also needs to be accounted for.

As it is usually associated with star-forming regions, the 24 µm continuum is often used directly as a tracer of star-formation regions and of the SFR (e.g., Hao et al. 2011; Kennicutt & Evans 2012). However, for studies that require disentangling the different phases of the ISM in star-forming regions, the origin of 24 µm emission can also lead to an ambiguous interpretation, as we have seen. It is often implicitly considered to originate from the ionized phase only (e.g., Malhotra et al. 2001; Croxall et al. 2012), but our study demonstrated that it is not the case, and that the origin of the 24µm emission from PDRs should also be considered for a more accurate interpretation. By comparing the ratios [O iii]/24 µm and ([C ii]+[O i])/[O iii], it should be possible to derive a correction factor (which can also depend on other parameters, such as density) to apply to 24µm to separate the relative contribution from H ii regions and PDRs.

One of the interests of this work is that our method relies on IR tracers, which are not normally extinguished nor absorbed in most cases. In particular [C ii] and [O iii] are among the brightest FIR emission lines in galaxies, and are popular tracers being observed in high-redshift galaxies with ALMA. [C ii] observations have been carried out up to a redshift of z ~ 6 in the large ALPINE sample of Le Fèvre et al. (2020) and associated studies, (Carniani et al. 2020), or even at z ~ 7 or higher (Hashimoto 2020; Bouwens et al. 2022). Likewise, the [O iii] line has also been detected at z ~ 7 (e.g., Hashimoto 2020; Bouwens et al. 2022). The rest-frame 24 µm band has also been observed for galaxies at different redshifts with Spitzer and Herschel telescopes (up to z ~ 3.5, see e.g., Yan et al. 2007; Oliver et al. 2012; Liu et al. 2021; Nagaraj et al. 2021; Schouws et al. 2022, among others), and is accessible with current facilities, such as FIFI-LS (Fischer et al. 2018), GREAT (Heyminck et al. 2012; Risacher et al. 2016) and HAWC+ (Harper et al. 2018) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA, Young et al. 2012). In the future, the Fred Young Submillimeter Telescope (FYST) may be able to access 24 µm in galaxies near the epoch of reionization.

The origin of the [C ii] line can be uncertain, more so in high metallicity galaxies. It can be emitted by both neutral and ionized gas as it requires 11.3 eV to create C⁺. However, particularly in low metallicity galaxies the contribution of ionized gas to [C ii] is often small (typically less than 30%, Cormier et al. 2019).

Considering the high redshift observations, our new method would thus certainly be convenient to derive a picture of the physical conditions in distant low-metallicity galaxies, where these few tracers already exist. There are, however, additional caveats to keep in mind when using this method to probe PDR and H ii region distributions.

First, we assume that in a given pixel the different phases (H ii regions, PDRs and neutral gas) have the same physical conditions, in particular the same density at the illuminated edge of the cloud and the same ionization parameter. Considering that the ISM is clumpy, it is not likely that all of the clumps along the line of sight have the same density, or the same illumination conditions. The assumption of constant pressure through all the models can also be questioned. The density profile may be more complex than these simple assumptions, and other effects can disturb the assumed equilibrium (winds, steps in pressure, turbulence, flows, etc.).

Another caveat is that our models are computed up to A_V of 2.5. This can lead to very high densities (up to a few 10⁶ cm⁻³) in the PDR, which can induce strong [C ii] and [O i] emission. It is probable that not all of the PDRs attain those high depths, and an investigation of different final cuts of the model have to be taken into account eventually.

Bringing more diagnostic tracers into the modeling scheme may also help to refine the results and make them more robust. Spitzer has observed numerous MIR H ii region tracers such as [Ne ii] and [Ne iii], [Si ii], iron lines and molecular hydrogen lines in some regions, potentially helping to further constrain the conditions in PDRs and H ii regions mixed within the telescope beam. JWST will sure bring us a rich array of these tracers. Additionally, velocity resolved observations to separate different clumps and phases, would also bring further information on the phase distributions. Despite all of these challenges, the method presented here brings a new and simple angle to study and disentangle the emission from the ionized and neutral gas in star-forming regions, with the possibility of application to other unresolved low metallicity galaxies, even at high redshift.

Fig. 19 Relation between [O iii] and MIPS 24 µm. The points are color-coded with the value of log(cPDR). The different panels correspond to the four bins applied on log(cPDR), namely [-1,0] (upper left panel), [0,1] (upper right panel), [1,1.5] (lower left panel), [1.5,2] (lower right panel). The black solid line for each panel correspond to the fitted relation on the bin, the black dashed lines represent the 95% confidence interval, and the gray solid line is the relation fitted in Fig. 1 (Sect. 3.1).

7 Summary

We examined the relation between the Herschel [O iii] 88 µm and the Spitzer 24 µm continuum band for a sample of star-forming regions that have been mapped in the Magellanic Clouds, using a novel and simple method to spatially disentangle different phases along the line of sight. We also compared our spatially-resolved results with those of the unresolved observations of the Herschel Dwarf Galaxy Survey. Out findings include the following:

• The [O iii] and Spitzer 24 µm continuum show a clear correlation covering three orders of magnitude for 24 µm and four orders of magnitude for [O iii]. We find, in surface brightness: $$\log \left({\Sigma }_{\left[\text{OIII}\right]88}\right)=0.879\left(\pm 0.012\right)\times \log \left({\Sigma }_{24}\right)-1.446\left(\pm 0.065\right),$$(6)

  with (Σ_[OIII]88) and Σ₂₄ in Wm⁻² sr⁻¹) and with L_[OIII] and L_[O24] in luminosity (L_⊙): $$\log \left(L_{\left[\text{OIII}\right]88}\right)=0.948\left(\pm 0.086\right)\times \log \left(L_{24}\right)-0.622\left(\pm 0.360\right).$$(7)

• To investigate the local conditions driving this relation, we compared our observations with a grid of 1D Cloudy models containing an H ii region component and a PDR component, and vary the stellar age, the initial density (at the illuminated edge of the cloud), n₀ and the ionization parameter, U. We also created a new parameter, the proportion of PDR, c_PDR, to determine the proportion of emission arising from PDRs and that with an origin in H ii regions along each line of sight.

• We used the ratio ([C ii]+[O i])/[O iii] as a proxy for the ratio of PDR vs. H ii region emission, and compared it to the [O iii]/24 µm. The use of [O iii]/24 µm and [O iii]/70 µm together allows us to constrain the models most efficiently.

• While the 24 µm continuum is often considered to trace dust specifically in H ii regions, we showed that it can also arise in PDRs, at least in low metallicity environments investigated here. Thus its use as a probe of H ii regions only should be treated with caution.

• The PDR proportion, c_PDR, is spatially anticorrelated with the [O iii] line emission, emitted only by H ii regions, and not necessarily correlated with the two PDR tracers, [C ii] and [O i]. This suggests that c_PDR is not tracing the simple presence of PDRs. Instead, the c_PDR parameter indicates whether PDRs or H ii regions are the dominant phase along the line of sight.

• We found that c_PDR correlates with the ratio [C ii]/([C ii]+[O i]). This suggests that regions of high c_PDR values, where PDRs dominate the line of sight emission, originate in environments of relatively lower illumination, thus further away from the young stars that ionize their surroundings.

We demonstrated the predictive power of [O iii] and 24 µm for both integrated and spatially resolved Magellanic Cloud regions and for the unresolved DGS data. Areas where the predicted [O iii] is significantly higher than the actual observations are dominated by PDRs, with higher c_PDR values. On the other hand, areas with the lower c_PDR values give predictions closely resembling the observations (observed and predicted values agree within a factor 3 for 90% of the data), although the predicted values have to be taken with caution for very low [O iii] fluxes (≤3 × 10⁻⁸ Wm⁻² sr⁻¹), and for places where the ionized gas may not be dominant.

Acknowledgements

We would like to thank the anonymous referee for useful comments that greatly helped to improve the quality of the paper. This work is based in part on Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

Appendix A Analytic Examples of the Modeling Strategy

Following on the description of the modeling strategy described in Section 4, we elaborate further on other examples of mixing different proportions of H ii regions plus PDR Cloudy models, related to the illustrations of Figure 10 which correspond to other physical configurations along the line of sight.

We rewrite Equation 4 in terms of I_single, which denotes a single, full Cloudy model. We then scale the full single Cloudy model along with the fraction of PDR intensity (I_PDR) necessary to match the balance of H ii region plus PDR.

Appendix A.0.1 Example: 90% H ii Region and 10% PDR

For a mixed line of sight dominated by an H ii region component, Equation 4 gives: $$\{\begin{array}{ccc}{I}_{\text{single}}& =& I_{\text{HII+PDR}}\hfill \\ 0.9I_{\text{HII}}+0.1I_{\text{PDR}}& =& 0.9I_{\text{HII}}+0.9I_{\text{PDR}}-0.8I_{\text{PDR}}\hfill \\ & =& 0.9I_{\text{single}}-0.8I_{\text{PDR}}\hfill \\ & =& 0.9\left(I_{\text{single}}-\frac{0.8}{0.9}\times I_{\text{PDR}}\right),\hfill \end{array},$$(A.1)

In Equation A.1, we see that a proportion of 10% of PDR can be simulated as a full Cloudy model with part of the PDR suppressed, assuming all of the modeled ISM have identical conditions (same ionizing sources, same initial density and same initial ionization parameter). The physical interpretation of a single Cloudy model, with part of the PDR removed, may represent a matter-bounded region, where the PDR is decreased due to lack of gas and dust.

Appendix A.0.2 Example: 50% H ii Region and 50% PDR

For an equal mix of H ii region and PDR, Equation 4 gives: $$\begin{array}{ccc}0.5I_{\text{H\hspace{0.17em}II}}+0.5I_{\text{PDR}}& =& 0.5I_{\text{single}}\end{array}.$$(A.2)

We see in Equation A.2 that a mixture with 50% of PDR (Figure 10) corresponds to a single Cloudy model, with all of the PDR emission taken into account, still under the assumption that the H ii region and the PDR component have the same initial density and ionization parameter. We highlight the point that 50% of PDR does not correspond to a medium with half of the volume being filled by PDR, but corresponds to a mix that can be modeled by a Cloudy model with a full PDR layer. This configuration can be physically interpreted as lines of sight through the different zones of a star-forming cloud with a young star cluster, and may represent a radiation-bounded cloud.

Appendix A.0.3 Example: 10% H ii Region and 90% PDR

For a mix dominated by PDR region, Equation 4 gives: $$\begin{array}{lll}0.1I_{\text{H\hspace{0.17em}II}}+0.9I_{\text{PDR}}\hfill & =\hfill & 0.1I_{\text{H\hspace{0.17em}II}}+0.1I_{\text{PDR}}+0.8I_{\text{PDR}}\hfill \\ \hfill & =\hfill & 0.1I_{\text{single}}+0.8I_{\text{PDR}}\hfill \\ \hfill & =\hfill & 0.1I_{\text{single}}+0.8I_{\text{PDR}},\hfill \end{array}$$(A.3)

In Equation A.3, we see that 90% of PDR in the mixture can be represented by a full Cloudy model with a further addition of eight times more PDR in this case. Those PDRs are connected to H ii regions lying outside the line of sight. Such a model also represents a radiation-bounded case along the corresponding line of sight between the stellar source and the cloud. Physically, this configuration can not be modeled by a single Cloudy model, as it would overestimate the emission of tracers emitted by the H ii region.

Appendix B Tests on the Strength of the Fit

We investigate the robustness of the relation between [O iii] 88µm and MIPS 24µm, studied in Section 3.1 when changing various aspects of the analysis: the pixel size, the method used for the fitting, the correction for point sources and the resolution at which we analyze the maps.

Appendix B.1 Test of the Fitting Method

Here we compare the effect of the fitting method of the parameters of the linear relation of [O iii] and 24µm discussed in Section 3. A simple, non weighted χ², method is made to investigate the uncertainties due to the fitting method. Another method uses a Monte Carlo simulation to calculate the uncertainties, which can then be considered to be directly linked to the observation uncertainties. Finally, the weighted χ² used in Section 3.1, takes into account uncertainties from the data and the method of fit. The fitted slopes and intercepts for each method, and their uncertainties, are shown in Table B.1.

The Monte Carlo simulation used the observed sample as a base to run 100 simulated samples, with the same number of pixels, each varying inside the limits of 1σ, following a normal distribution. The fit is then calculated with a nonweighted χ² on each of the simulated samples, and the final value for the parameters are the median from all of the resulting fittings, with the median absolute deviation being the uncertainty. We added the value of the R² coefficient of the fitted linear regression, to estimate the difference in fitting accuracy. The three fitted relations are illustrated with the data in Figure B.1.

The parameters fitted with different methods are very similar and roughly agreeing within the uncertainties (Table B.1). We also see that the R² values are high for all the methods, between 80 and 90%, which correspond to a high probabilities of a linear positive correlation. The different methods show very similar relations when compared to the fitted data (Figure B.1).

The uncertainties derived by the two χ² methods are almost identical, and much larger than the uncertainties derived by the Monte Carlo method. It indicates that the uncertainties derived by the fitting procedure and due to the spread of the data are dominant compared to the effect of observation uncertainties.

We can see that the R² parameters are very similar for the two χ² methods, but it is a bit higher for the Monte Carlo simulation, indicating that the Monte-Carlo fit is a bit more representative of the data, but this effect is very small, as the difference between the R² values is smaller than 0.1.

Fig. B.1 Comparison of the linear fits between [O iii] and 24µm for the three fitting methods. The 95% confidence intervals are represented with dashed lines and gray shading for the simple and weighted χ2 methods. The confidence interval is not determined for the Monte Carlo method. The upper right panel, corresponding to the weighted χ2, is identical to Figure 1 in Section 3.1, showing the fit of the correlation.

We can conclude from this test that the different methods give similar results, both on the fitted value and the accuracy of the fit. The comparison of the nonweighted χ² and the Monte Carlo method shows that the uncertainties on the fitted parameters driven by the fitting method are much larger than the uncertainties driven by data uncertainties. The weighted χ², which is the method finally adopted for our study, has the advantage of taking into account both types of uncertainties, although the uncertainties driven by the fitting method are dominant.

Appendix B.2 Test on the Resolution and Sampling of the Maps

We compare different spatial resolutions and pixel sampling sizes of the maps, to investigate the potential effects it can have on the fit. We use a weighted χ² method to be consistent with that used in Section 3.1, and we compare two different resolutions with the same sampling, and then two different samplings with the same resolution. The resolutions we use are 12", to match the PACS emission lines, and 18" resolution, to match the MIPS 70µm. For both resolutions, we resample the maps to 12" pixel width, to avoid oversampling of the beam. For the spatial sampling, we compare our final pixel size of 12" to the native sampling of PACS spectroscopy, 3".1. It also allows us to investigate the effect of oversampling the beam, as we use 12" resolution for the tests on sampling size (Figure B.2).

The different resolutions and sampling sizes of the maps return very similar results for the fits of the parameters of the relation (Table B.2), which roughly agree within the uncertainties. The R² parameters are also very similar for each of the data settings.

The derived uncertainties seem to be affected by the pixel width. For the largest pixel size (12"), the uncertainties on the fitted parameters are 3 times larger than those for the larger pixel width (3".1). Thus, for comparable dispersion, the fit can be more precise where it is carried out on smaller pixels.

While the fits from Figure B.2 (Table B.2) do not differ much, the confidence intervals are affected by both the resolution and the pixel size: at the low resolution and low pixel width (12" resolution and 3".1 pixel width), the confidence intervals are highly irregular, especially at the low brightness end of the relation, but they become more regular when increasing the pixel size (12" resolution and 12" pixel width) and when increasing the beam size (18" resolution and 12" pixel width). It corresponds to a smoothing of the data.

We can conclude from this test that the effect of resolution and pixel size is negligible on the fitted parameters of the relation between [O iii] 88µm and MIPS 24µm, even when the beam is oversampled. Although using a larger number of pixels, leads to a somewhat better precision of the fitted values.

Appendix B.3 Test on Method for Point Source Correction

The MIPS 24µm continuum band contains emission from stellar point sources as well as dust emission. We investigate the effect of the point source removal, and the effect of the method used, on the fitting of the relation between [O iii] and MIPS 24µm.

Fig. B.2 The fitted linear relation and the corresponding data for the different resolution and sampling of the maps. The black solid line is the fit, the dashed black lines and gray filling illustrate the 95% confidence interval. The upper right panel, corresponding to 12" resolution and a pixel size of 12", is identical to the Figure 1 in Section 3.1, showing the fit of the correlation.

We developed a dedicated software to extract the emission of the point sources from the maps. The position of the point sources are taken from catalogs. A point source is considered to be a 2D Gaussian. The emission from each point of the map is then decomposed in two components: the background, and the contribution from the point source. The map is divided into tiles and the background is considered to be flat on each tile. The size of the tiles depends on the size of the total map and the spread of the point source emission. The tiles roughly correspond to the 12" pixel except for the 30 Doradus region. The background value is estimated on each tile, in order to fit the 2D Gaussian function for each point source referenced. The software produces three maps as outputs: the tiled background map, a "zoom" background map, which is based on the tiled background and reprocessed to improve the resolution, and a map of the point source emission, at the resolution of the initial map. This point source map is subtracted from the initial map, to retrieve the extended, diffuse emission only. This method comes with some caveats, especially in the case of overestimation of the point source emission. It can lead to some artifacts, or some pixels with negative emission, which are artificially set to zero, but a detailed examination of our sample showed that this type of artifact is not present or is negligible in the studied regions.

In order to study the effect of the method used to correct maps from point source emission, we fitted the relation between [O iii] and MIPS 24µm by using four 24µm emission maps: the original MIPS 24µm maps, convolved and resampled, without removing the point sources (later called original data), the two types of diffuse emission maps created by the software, namely the tiled background and zoomed background, and the original data with the point source maps subtracted, which is referred to as point source subtracted data or simply subtracted data. In region N160, the 24µm data show a large saturated zone. The saturated pixels were masked, and when trying to apply a point source correction, no point source was found in nonsaturated pixels. Thus the N160 24µm masked map is used in all of the four procedures, without attempts to correct from point source emission. Table B.3 summarizes the parameters fitted for the four 24µm maps, and the results are illustrated in Figure B.3.

We notice that the global dispersion is not affected by the method used to correct the 24µm maps from point source emission (Figure B.3), indicating that it is linked to the discrepancies between the different regions, as we have already noticed in Section 3.1. However, the internal dispersion for each region is affected by the removal of point sources. The two cases of tiled and zoomed background emission reduce the spread for some regions, like N iii or N159, for the MIPS 24µm data (roughly from 1 dex or a bit more to 0.5 dex); it naturally leaves untouched the spread in [O iii] emission. Some regions, on the contrary, exhibit a spread a bit larger with the tiles and zoomed background than with the original MIPS 24µm data, such as 30 Doradus. This increase of the spread is probably due to the tiling of the diffuse emission used to calculate the point source emission, and thus is identical for the two diffuse emission calculations. The effect is most important for the 30 Doradus region, probably due to the fact that this is the only region for which the tiling is done on a scale much larger than the resampling of the maps at 12" pixel width.

Fig. B.3 Results of the fit with different treatment of point-sources in the 24µm maps. The original MIPS 24µm maps (upper left panel), tiled (lower left panel) and smoothed background (lower right panel) calculated by the software extracting the point sources, and the original maps with emission from the extracted point sources subtracted (upper right panel). The two first maps of the diffuse emission are outputs from the software, where the last one is the software estimation of the point sources emission subtracted from the initial map. The black solid line represents the fitted relation between the [O iii] emission line and the 24µm band and the black dashed lines and the gray shading indicate the 95% confidence intervals. Fitted parameters for the different cases are in Table B.3

The diffuse emission obtained by subtracting the point source emission from the original maps also does not influence the global spread, but only the spread of each individual region. In this case, we observe the same reduction of the spread for the regions affected by the other point source correction methods, but there is no increase of the spread in the other regions.

Removing the point source emission has little effect on the fitted parameters, independent of the method used (Table B.3). The slopes and intercepts are very similar, and completely agreeing within the uncertainties. The R² parameter is somewhat affected: the method of original data, with no point sources subtracted, gives the least accurate fit, the method using two different diffuse emission calculations gives a result a bit more accurate, and the best fitted relation is achieved for the maps with the point source emission subtracted from the original maps.

From this test, we conclude that taking into account the contamination of the 24µm emission by point source emission is important, especially to investigate the behavior and conditions in the star-forming regions of the individual regions, but it does not influence much the fitting of the parameters of the overall relation between [O iii] emission and MIPS 24µm, nor the uncertainties on the fitted parameters. To avoid an increase of the spread in some individual regions, and the loss of spatial resolution due to the tiling used to calculate the point source emission, we use the MIPS 24µm data with point sources subtracted.

Appendix C Maps and specific plots for the individual regions

Fig. C.1 For the region 30 Dor, the predicted map of [O iii] given the observed 24µm (panel a) is compared to the actual map of [O iii] observed with Herschel/PACS (panel b), and the ratio of observed over predicted emission (panel c). The distribution of [C ii] and [O i] emission (panels d and e respectively) are overlaid with the [O iii] line contours. The ratio of RCII/PDR (panel f), RPDR/HII (panel g) and cPDR with associated uncertainties (panel h and i respectively) are also presented with the [O iii] line contours. The observed ratios [O iii]/24µm and [O iii]/70µm presented as a function of RPDR/HII are compared to our grid of models (panelj and k respectively). The large cross is the value integrated on the global region.

Fig. C.2 Similar to Figure C.1, but for the region N11B.

Fig. C.3 Similar to Figure C.1, but for the region N11C.

Fig. C.4 Similar to Figure C.1, but for the region N11I.

Fig. C.5 Similar to Figure C.1, but for the region N44.

Fig. C.6 Similar to Figure C.1, but for the region N158.

Fig. C.7 Similar to Figure C.1, but for the region N159.

Fig. C.8 Similar to Figure C.1, but for the region N160.

Fig. C.9 Similar to Figure C.1, but for the region N180.

Fig. C.10 Similar to Figure C.1, but for the region N66.

References

All Tables

All Figures

	Fig. 1 Relation between [O iii] and 24 µm in all of our SMC and LMC regions, only including pixels above 5σ in the [O iii]. The black solid line represents the fitted relation, linear in log-log space, and the dashed lines and gray filling represent the 95% confidence interval of the fit. The colored symbols refer to the individual regions labeled in the legend to the right of the plot.
In the text

	Fig. 2 Relation between surface brightness of [O iii] and 24 µm for all of the pixels of our individual star-forming regions detected at more than 5σ. The black solid line represent the fitted linear relation in loglog space for the total sample (Fig. 1), and the dashed lines and gray filling represents the 95% confidence interval associated with the total sample.
In the text

Fig. 3 Spatial distribution of the low emission pixels in N11C. Left panel: [O iii] map for N11C; right panel: MIPS 24 µm map for N11C. Both maps are convolved to 12″ and resampled to 12″ pixels. The black crosses correspond to pixels that have low surface brightness (on the left of the bulk of the regions in Figs. 1 and 2, corresponding to values of 24 µm lower than 106 W m−2 sr−1) showing a different behavior from that of the global relation.

In the text

Fig. 4 Upper panel: relation between the integrated [O iii] and 24 µm luminosity, in units of solar luminosity (L⊙). Lower panel: integrated luminosity per surface area of [O iii] and 24 µm (Wkpc−2). The DGS sample corresponds to the gray circles, and the integrated regions from the Magellanic Cloud sample are coded with colors, as indicated in the inserted legend in the upper panel. The fit on the Magellanic Clouds alone is indicated as MCs, the fit on DGS sample is indicated as DGS, and the fit on the whole data is indicated as DGS+MCs.

In the text

	Fig. 5 Predicted 30 Dor map of [O iii] (top panel) given the observed 24 µm and Eq. (1). For comparison, the actual observed map of [O iii] (middle panel), and the ratio of observed over predicted emission (lower panel). The two maps of [O iii] are represented on the same scale and are displayed in Wm−2 sr−1.
In the text

	Fig. 6 Similar to Fig. 5, but for N11I.
In the text

	Fig. 7 Spatially resolved ratios of observed ([O iii]obs) over predicted ([O iii]pred) emission of [O iii]. The value of 1, corresponding to good agreement of the predictions and the observations, is indicated with a horizontal dashed line. The horizontal gray lines indicate differences of a factor 3 from agreement.
In the text

	Fig. 8 Predicted [O iii] luminosity, from the observed 24 µm for integrated luminosities of the Magellanic Clouds sample and DGS galaxies compared to the observed luminosity values. The horizontal dashed line indicates when predicted and observed values are equal. The horizontal gray lines indicate differences of a factor 3 from agreement.
In the text

Fig. 9 Example from our Cloudy models. The cumulative emission of different bands and emission lines are presented (upperpanel), together with some physical parameters (density, electronic density and temperature, lower panel) as a function of the depth in the cloud, measured with the parameter, Ay. The emission of all of the tracers was normalized to a maximum of 1 to allow a direct comparison. The age of the central ionizing source is 4 Myr, and the initial values of physical parameters are log(n0)= 2.5 and log(U) = −3.

In the text

Fig. 10 Sketch of different PDR proportions (cPDR) in our mixed models. The yellow regions represent H ii regions, and PDRs are represented in gray. The blue arrow indicates the direction of increasing cPDR. The black arrows show the lines of sight through the represented clouds (each line of sight representing a given pixel in our maps). As per Eq. (4), the 0% cPDR case corresponds to no PDR in the model, hence fully ionized gas; 50% cPDR corresponds to a single Cloudy model, with a connected H ii region and PDR (Eq. (A.2)); and 100% cPDR represents a full PDR model, without an H ii region contribution. We note that 50% of PDR does not correspond to a medium with half of the volume being filled by PDR, but corresponds to a mix that can be modeled by a Cloudy model with a full PDR layer.

In the text

Fig. 11 [O iii]/24 µm for different models as a function of burst age. Each model is represented with a color for n0, and a line style for U, and the lines are following the proportion of PDR emission, cPDR, with the markers representing the cPDR of 0%, 50 % and 99%. The blue solid horizontal line corresponds to the median of the spatially resolved observations and the dashed blue lines to the median absolute deviation. The distribution of our observations is shown on the right as a histogram.

In the text

Fig. 12 Representation of the data of the full sample and models with constant pressure for [O iii] 88 µm/MIPS 24 µm as a function of ([C іі]+[O i])/[O iii] (top panel), and [O iii] 88 µm/MIPS 70 µm as a function of ([C ii]+[O i])/[O iii] (lower panel). Each model is represented with a color for the initial density, n0, and a line style for U, and the lines follow the cPDR, with the markers representing the proportions of 0%, 50% and 99%. The models displayed are calculated for an age of star-forming burst of 4 Myr. The gray points are the observed values for each pixel of the regions.

In the text

	Fig. 13 Similar to Fig. 12, but for the region N44.
In the text

Fig. 14 Distribution of [C ii] and [O i] emission (upper left and right panels respectively) in 30 Doradus; the ratio of ([C ii]+[O i])/[O iii] (middle left panel), the distribution of cPDR (middle right panel) with the colorbar indicating the percentage of PDR along the line of sight and the ratio [C ii]/([C ii]+[O i]) (lower panel). The contours of [O iii] are overlaid on all the panels for comparison, with the color intensity ranging from dark gray for lower emission to light gray for higher [O iii] emission.

In the text

	Fig. 15 Similar to Fig. 14, but for the region N11I.
In the text

Fig. 16 Relation between the PDR proportion and the ratio [C ii]/([C ii]+[O i]), for all of the pixels in all of the regions studied here, with the same region color code as in Fig. 1. The 1% and 10% cPDR values are indicated by the two dashed horizontal lines. The vertical gray line separates regions where cooling is dominated by [C ii] (right) and regions where the cooling is dominated by [O i] (left).

In the text

	Fig. 17 Ratio of [O iii]obs/[O iii]pred represented as a function of [O iii]obs, similar to Fig. 7, except that the color code corresponds to the log(cPDR).
In the text

	Fig. 18 Relation between [O iii] and MIPS 24 µm in all of the pixels in all of our star-forming regions for those pixels above 5σ in the [O iii] line. This figure is similar to Fig. 1, but the points are color-coded with the value of log(cPDR) here.
In the text

Fig. 19 Relation between [O iii] and MIPS 24 µm. The points are color-coded with the value of log(cPDR). The different panels correspond to the four bins applied on log(cPDR), namely [-1,0] (upper left panel), [0,1] (upper right panel), [1,1.5] (lower left panel), [1.5,2] (lower right panel). The black solid line for each panel correspond to the fitted relation on the bin, the black dashed lines represent the 95% confidence interval, and the gray solid line is the relation fitted in Fig. 1 (Sect. 3.1).

In the text

Fig. B.1 Comparison of the linear fits between [O iii] and 24µm for the three fitting methods. The 95% confidence intervals are represented with dashed lines and gray shading for the simple and weighted χ2 methods. The confidence interval is not determined for the Monte Carlo method. The upper right panel, corresponding to the weighted χ2, is identical to Figure 1 in Section 3.1, showing the fit of the correlation.

In the text

	Fig. B.2 The fitted linear relation and the corresponding data for the different resolution and sampling of the maps. The black solid line is the fit, the dashed black lines and gray filling illustrate the 95% confidence interval. The upper right panel, corresponding to 12" resolution and a pixel size of 12", is identical to the Figure 1 in Section 3.1, showing the fit of the correlation.
In the text

Fig. B.3 Results of the fit with different treatment of point-sources in the 24µm maps. The original MIPS 24µm maps (upper left panel), tiled (lower left panel) and smoothed background (lower right panel) calculated by the software extracting the point sources, and the original maps with emission from the extracted point sources subtracted (upper right panel). The two first maps of the diffuse emission are outputs from the software, where the last one is the software estimation of the point sources emission subtracted from the initial map. The black solid line represents the fitted relation between the [O iii] emission line and the 24µm band and the black dashed lines and the gray shading indicate the 95% confidence intervals. Fitted parameters for the different cases are in Table B.3

In the text

Fig. C.1 For the region 30 Dor, the predicted map of [O iii] given the observed 24µm (panel a) is compared to the actual map of [O iii] observed with Herschel/PACS (panel b), and the ratio of observed over predicted emission (panel c). The distribution of [C ii] and [O i] emission (panels d and e respectively) are overlaid with the [O iii] line contours. The ratio of RCII/PDR (panel f), RPDR/HII (panel g) and cPDR with associated uncertainties (panel h and i respectively) are also presented with the [O iii] line contours. The observed ratios [O iii]/24µm and [O iii]/70µm presented as a function of RPDR/HII are compared to our grid of models (panelj and k respectively). The large cross is the value integrated on the global region.

In the text

	Fig. C.2 Similar to Figure C.1, but for the region N11B.
In the text

	Fig. C.3 Similar to Figure C.1, but for the region N11C.
In the text

	Fig. C.4 Similar to Figure C.1, but for the region N11I.
In the text

	Fig. C.5 Similar to Figure C.1, but for the region N44.
In the text

	Fig. C.6 Similar to Figure C.1, but for the region N158.
In the text

	Fig. C.7 Similar to Figure C.1, but for the region N159.
In the text

	Fig. C.8 Similar to Figure C.1, but for the region N160.
In the text

	Fig. C.9 Similar to Figure C.1, but for the region N180.
In the text

	Fig. C.10 Similar to Figure C.1, but for the region N66.
In the text