Disentangling emission from star-forming regions in the Magellanic Clouds: Linking [OIII]88 micron and 24 micron

This study explores the link between the [OIII]88mu emission, a well-known tracer of HII regions, and 24mu continuum, often used to trace warm dust in the ionized phases of galaxies. We investigate the local conditions driving the relation between those tracers in the Magellanic Clouds, comparing observations with Cloudy models consisting of an HII region plus a photodissociation region (PDR) component, varying the stellar age, the initial density (at the illuminated edge of the cloud), and the ionization parameter. We introduce a new parameter, cPDR, to quantify the proportion of emission arising from PDRs and that with an origin in HII regions along each line of sight. We use the ratio ([CII]+[OI])/[OIII] as a proxy for the ratio of PDR versus HII region emission, and compare it to [OIII]/24mu. The use of [OIII]/24mu and [OIII]/70mu together allows us to constrain the models most efficiently. We find a correlation over at least 3 orders of magnitude in [OIII]88mu and 24mu continuum in spatially resolved maps of the Magellanic Cloud regions as well as unresolved galaxy-wide low metallicity galaxies of the Dwarf Galaxy Survey. Most of the regions have low proportions of PDRs along the lines of sight (<12%), while a limited area of some of the mapped regions can reach 30 to 50%. For most lines of sight within the star-forming regions we have studied in the Magellanic Clouds, HII regions are the dominant phase. We propose the use of the correlation between the [OIII]88mu and 24mu continuum as a new predictive tool to estimate, for example, the [OIII]88mu emission when the 24mu continuum is available or inversely. This can be useful to prepare for ALMA observations of [OIII]88mu in high-z galaxies. This simple and novel method may also provide a way to disentangle different phases along the line of sight, when other 3D information is not available.


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. 2017Matthee et al. , 2019Carniani et al. 2018Carniani et al. , 2020Hashimoto 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 con-ditions 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).
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: The most luminous farinfrared (FIR) line observed in lowmetallicity 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. 2019Sutter et al. , 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(Cormier et al. , 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. 2007Bendo et al. , 2008Bendo et al. , 2012Croxall 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 (50 kpc 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 1997Pagel , 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 Section 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 Section 3. In Section 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 Section 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 Section 6, and we conclude our work in Section 7.

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.

Spitzer data
We used the Spitzer IRAC and MIPS maps available on the IRSA Infrared Science Archive 1 and selected subregions covering our 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, that 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.

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) andPACSman (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 −3 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, repectively 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).

Relation between [O iii] 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.

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 χ 2 method, leading to the following relation: where the surface brightness, Σ, is expressed in W m −2 sr −1 . 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 (Figure 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 (Figure 2), for example, the slope begins to steepen when the 24µm surface brightness is less than ∼10 −6 W m −2 sr −1 . These particular low-surface brightness pixels lie on the outskirts of the star-forming peak in N11C (see Figure 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 ( Figure 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 starforming regions from which the 24µm originates. More studies on other regions in the SMC would be necessary to conclude on this point. 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 Figures 1 and 2, corresponding to values of 24µm lower than 10 6 W/m 2 /sr) showing a different behavior from that of the global relation.
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 Section 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 pointsource removal methods: a) effect of different fitting methods: comparison of a simple, nonweighted χ 2 , and a Monte Carlo simulation to determine uncertainties and the weighted χ 2 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 Figure 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 Section 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.

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 Figure 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 χ 2 method on those star-forming regions, we obtained the relation of Equation 2 for the integrated Magellanic Cloud regions: log(L [OIII] 88 ) = 0.948(±0.086) × log(L 24 ) − 0.622(±0.360), (2) with L [OIII] 88 being the luminosity of [O iii] and L 24 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 (Equation 1) within the uncertainties. We also investigate the relation with the surface luminosity.  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 Section 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. 2003Galliano et al. , 2005Madden et al. 2006;Cormier et al. 2015Cormier et al. , 2019Madden et al. 2020;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

Predictive power of the [O iii] -24µm calibration
We used our two calibrations from Equations 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 val- We first compared the [O iii] predictions we derive from the observed 24µm using surface brightness, Equation 1, to the observed [O iii] maps of the Magellanic Cloud regions. As examples, we show in Figures 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 −8 W/m 2 /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 Figure 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 Figure 7   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 Section 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 Section 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 ( Figure 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 Section 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 nonnegligible spread in some cases, around the correlations derived from Figures 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.

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 0 . 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 stelar 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 0 , 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 0 ), impacting the inner surface of the cloud of radius, r, the speed of light, c and n 0 : n 0 is varied in log(cm −3 ) 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 2 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.
The ionization parameter, U, is defined in Equation 3 Fig. 9. Example from our Cloudy models. The cumulative emission of different bands and emission lines are presented (upper panel), 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, A V . The emission of all of the tracers were 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(n 0 )= 2.5 and log(U) = −3. 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 −2 in our example). The density is low in the H ii region, and increases in the PDR (up to a value of ∼ log(n 0 )= 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 Section 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 can not 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 2 , 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 equation 4: where I H II is the intensity of the tracer in the H ii region only in W m −2 sr −1 ; 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 HII+PDR ) is similar to a luminosity scaling, but we prefered 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 Equation 4. We note that this mixture is made only from models having the same initial n 0 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 Equation 4 are illustrated in Figure 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 Equations 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 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 Equation 4, the 0% c PDR case corresponds to no PDR in the model, hence fully ionized gas; 50% c PDR corresponds to a single Cloudy model, with a connected H ii region and PDR (Equation A.2); and 100% c PDR 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.
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 (Section 5.1), n 0 and U (Section 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 ( We first investigated the physical conditions of the different regions, and especially the parameters we used to create our grid of models: age, n 0 and U. Then, we used Equation 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 = ( , that we used as a proxy for the ratio of PDR over H ii region emission.

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 Myrs, and compared the observations to the predicted line ratio [O iii]/24µm for models that cover a range of n 0 , U and c PDR (Figure 11).
From Figure 11, we see that the predicted [O iii]/24µm ratio is roughly constant from 1 to ∼4 Myrs, 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 Myrs, 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 Myrs), 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 Figure 11, we can also begin to investigate the likelihood of some values of the density and ionization parameter. The models with log(n 0 ) = 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 0 ) = 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 0 ) ∼ 2) is more consistent for the bulk of the data. From Figure 11, it is not possible to conclude on the ionization parameter or the c PDR value.

Determining the ISM conditions
We then considered the model parameters, n 0 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 .

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 Figure 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 0 agreeing with all of the observations. The models that best agree with the entire data set have log(n 0 )=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 3 .
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.

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, that 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 0 (Table 5). It is probable that the density we report here corresponds to a rather diffuse ionized gas, that 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 5 -10 6 cm −3 at high A V , corresponding to the high-density clumps observed. Figure 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 0 value: log(n 0 ) = 1.5 while the other sources have log(n 0 ) 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 −3 ) 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 Figure 13, that other models with with lower U (log(U) = −3) and higher densities (log(n 0 ) = 2.5to3) 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 −3 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 0 ) = 2.5 to 3 agree well with the higher density structures of swept up material in the region.

N44 (
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 Figure 13 Figure 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  −3, 2 -As the observations trace a wide range in parameter space, other models of our grid can be consistent with each region (although less than our best model), or with a subpart of the region. Those models are indicated in the third column of the Table. of c PDR , to have a better understanding of how it can be linked to physical interpretation.

Insight on the distribution of ionized and neutral gas
For each region, assuming the best fitting model determined in Section 5.2, we calculated c PDR needed to match the observed  bach (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]   ter 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;Requena-Torres et al. 2016;Fahrion et al. 2017;Okada et al. 2019;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) (Figures 14, 15 and C.1 to C.10). We see (Fig. 16) that R C II/PDR is correlated with c PDR . When R C II/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 C II/PDR would be consistent with the c PDR value we derive: when R C II/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 C II/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.

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 Table 5, we see that only models with intermediate density (log(n 0 ) = 2 − 2.5) are relevant for those two regions, where models with lower density (log(n 0 ) ≤ 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 Figure 18, where we notice some spread of the data relating to the PDR proportion. The pixels with similar c PDR 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 4 different bins of PDR proportion with a weighted χ 2 . The results from the fitting are presented in Table  6. The data and the relation fitted for the different bins are shown in Figure 19. Table 6. Fit parameters of the relation between [O iii] and 24µm for different bins of log(c PDR ), with the values fit on the full sample noted in the last row, and 1σ uncertainty on the slope and intercept. Those fits are represented for each c PDR bin in Figure 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 Figure 1 (Section 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 Equation 1 in Section 3.1). On the other hand, the relation for the bin of highest c PDR values (log(c PDR ) =1.5 to 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 c PDR 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.

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 0 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 parame-ters 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. 2020), or investigated the geometry and properties of the ISM with velocity-resolved studies (see e.g., Okada et al. 2015Okada et al. , 2019Requena-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 starformation 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, that 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. 2020and associated studies, Carniani et al. 2020, or even at z∼7 or higher(Hashimoto 2020Bouwens et al. 2021). Likewise, the [O iii] line has also been detected at z∼7 (e.g., Hashimoto 2020; Bouwens et al. 2021). 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 6 cm −3 ) 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.

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: (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 0 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 ( Dale, D. A., Sandstrom, K., et al. 2021, Monthly Notices of the Royal Astronomical Society, 503, 911 Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, Annual Review of Astronomy and Astrophysics, 58, 157 Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, The Astrophysical Journal, 874, 27 Tarantino, E., Bolatto, A. D., Herrera-Camus, R., et al. 2021 , 338, 687 Vutisalchavakul, N. & Evans, Neal J., I. 2013, ApJ, 765, 129 Walker, A. R. 2012, Astrophysics and Space Science, 341, 43 Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, The Astrophysical Journal Supplement Series, 154, 1 Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, The Astrophysical Journal, 716, 1191Yan, L., Sajina, A., Fadda, D., et al. 2007, The Astrophysical Journal, 658, 778 Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, The Astrphysical Journal Letters, 749, L17 Zanella, A., Daddi, E., Magdis, G., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 1976   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. 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.  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 χ 2 , 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 χ 2 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 χ 2 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 2 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 2 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 χ 2 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 2 parameters are very similar for the two χ 2 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 2 values is smaller than 0.1.
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 χ 2 and the Monte Carlo method shows that the uncertainties on the fitted parameters driven by the fitting method are much larger than the un-certainties driven by data uncertainties. The weighted χ 2 , 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 χ 2 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 2 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.  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 Original" is for noncorrected data, "tiled" and "zoomed diffuse emission" are the two calculation of the diffuse emission from the software used to extract the point sources, and "point source subtracted" is for the original MIPS 24µm maps with point sources subtracted.
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 N11I 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 Article number, page 23 of 35 A&A proofs: manuscript no. Dust_and_gas_in_HII_regions_journal_v4  Table B.3 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.
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 2 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 resolu-