Calibrating mid-infrared emission as a tracer of obscured star formation on HII-region scales in the era of JWST

Measurements of the star formation activity on cloud scales are fundamental to uncovering the physics of the molecular cloud, star formation, and stellar feedback cycle in galaxies. Infrared (IR) emission from small dust grains and polycyclic aromatic hydrocarbons (PAHs) are widely used to trace the obscured component of star formation. However, the relation between these emission features and dust attenuation is complicated by the combined effects of dust heating from old stellar populations and an uncertain dust geometry with respect to heating sources. We use images obtained with NIRCam and MIRI as part of the PHANGS--JWST survey to calibrate dust emission at 21$\rm \mu m$, and the emission in the PAH-tracing bands at 3.3, 7.7, 10, and 11.3$\rm \mu m$ as tracers of obscured star formation. We analyse $\sim$ 20000 optically selected HII regions across 19 nearby star-forming galaxies, and benchmark their IR emission against dust attenuation measured from the Balmer decrement. We model the extinction-corrected H$\alpha$ flux as the sum of the observed H$\alpha$ emission and a term proportional to the IR emission, with $a_{IR}$ as the proportionality coefficient. A constant $a_{IR}$ leads to extinction-corrected H$\alpha$ estimates which agree with those obtained with the Balmer decrement with a scatter of $\sim$ 0.1 dex for all bands considered. Among these bands, 21$\rm \mu m$ emission is demonstrated to be the best tracer of dust attenuation. The PAH-tracing bands underestimate the correction for bright HII regions, since in these environments the ratio of PAH-tracing bands to 21$\rm \mu m$ decreases, signalling destruction of the PAH molecules. For fainter HII regions all bands suffer from an increasing contamination from the diffuse infrared background.


Introduction
Interstellar dust biases our view of star formation, obscuring the ultraviolet (UV) and optical emission of newly formed massive stars.The energy absorbed by dust at short wavelengths is reradiated in the IR, making observations of dust emission a powerful probe of embedded star formation.Multi-wavelength studies, at both low and high redshift, have highlighted the importance of the IR to provide a full accounting of star formation rates (SFRs) in galaxies (Calzetti et al. 2007;Wuyts et al. 2011b;Kennicutt & Evans 2012;Gruppioni et al. 2013;Madau & Dickinson 2014;Calzetti 2020).
The thermal dust emission spectrum peaks in the far-IR.In the limit where the dust heating is dominated by young stellar populations, the far-IR emission can be directly related to the embedded SFR using energy balance (Kennicutt 1998).However, such a bolometric measure of dust luminosity is sensitive to heating by older stellar populations (Cortese et al. 2008).This heating term contributes progressively more to the cooler, longer wavelength far-IR emission, leading to diffuse 'IR cirrus'.On the other hand, the mid-IR (MIR) region around 24 µm, corresponding to the IRAS 25µm, the Spitzer MIPS 24 µm band, or the WISE W4 (22 µm) filter, traces emission of small, hot grains.While not directly relatable to the bolometric dust luminosity without further assumptions on the shape of the dust SED, emis-⋆ e-mail: francesco.belfiore@inaf.itsion at 24 µm has the advantage of tracing a hotter dust component, and therefore being less sensitive to IR cirrus.
Several authors have advocated the use of 24 µm fluxes, in combination with either UV or Hα recombination line emission, to account for both the attenuated and un-attenuated component of the SFR (Hirashita et al. 2003;Kennicutt et al. 2007;Calzetti et al. 2007;Kennicutt et al. 2009;Wuyts et al. 2011a;Hao et al. 2011;Leroy et al. 2012;Kennicutt & Evans 2012;Jarrett et al. 2013;Catalán-Torrecilla et al. 2015;Leslie et al. 2018).These hybrid recipes have been demonstrated to be successful for both integrated galaxies and kpc-scale sub-galactic regions in local galaxies (Calzetti et al. 2007;Leroy et al. 2012;Belfiore et al. 2023).However, dust heating from old stellar populations can still constitute a significant fraction of a galaxy's integrated emission at 24 µm (Bendo et al. 2012;Groves et al. 2012;Crocker et al. 2013;Lu et al. 2015;Boquien et al. 2016;Viaene et al. 2017) and complicates the association of the MIR with star formation in observations that have low spatial resolution (Leroy et al. 2012).
In addition to 24µm, several authors have studied the potential of using the flux in PAH features or the brightness of PAHtracing bands (most commonly the IRAC 8 µm band) to probe obscured star formation (Peeters et al. 2004;Shipley et al. 2016;Maragkoudakis et al. 2018;Whitcomb et al. 2020).These calibrations offer significant advantages with respect to the longer mid-IR wavelengths in terms of angular resolution and sensitivity (Calzetti 2020;Li 2020).Moreover, PAH features remain ac-cessible and are readily identifiable with Spitzer, and now JWST, out to high redshifts (Lutz et al. 2005;Riechers et al. 2014).Studies of nearby galaxies at high spatial resolution have demonstrated, however, that PAH emission is suppressed relative to 24 µm (or total IR) emission at the position of H ii regions (Helou et al. 2004;Bendo et al. 2006;Bolatto et al. 2007;Lebouteiller et al. 2011;Chastenet et al. 2019Chastenet et al. , 2023b;;Egorov et al. 2023).In particular, PAH emission is found to generate bright rings around H ii regions, probably due to PAH photo-destruction, therefore complicating the relation between PAHs and star formation (Galliano et al. 2018).
In this work, we test the use of both the MIR dust continuum at 21 µm and bands dominated by PAH emission to trace embedded star formation on the scale of individual H ii regions (∼40−110 pc) in a sample of 19 nearby disc galaxies observed with the MIRI and NIRCam instruments onboard JWST as part of the PHANGS (Physics at High Angular resolution in Nearby Galaxies) JWST program (Lee et al. 2023).We compare the attenuation estimates obtained using literature calibrations based on the MIR/Hα ratios with matched-resolution Balmer decrement (Hα/Hβ) measurements obtained from the PHANGS-MUSE integral field spectroscopy (IFS) mapping with the VLT/MUSE instrument (Emsellem et al. 2022).
Using this JWST data we calibrate photometric measurements of the MIR continuum and PAH emission to accurately trace the obscured component of star formation on scales of individual H ii regions.Our analysis builds on the work of Leroy et al. (2023b) and Hassani et al. (2023), who used early data for a subset of four galaxies in our sample to demonstrate that attenuation-corrected Hα emission correlates well with MIR emission at 21 µm, albeit non-linearly, for regions that are bright in the MIR.
In Sec. 2 we present our data, paying special attention to the new JWST images, and summarise the methods for inferring the dust attention.In Sec. 3 we present our results concerning the use of the IR bands in correcting Hα for dust attenuation in H ii regions, and discuss the recommended recipes and potential future work in Sec. 4. We summarise the results in Sec. 5.

PHANGS-JWST data
We analysed 19 galaxies observed as part of the cycle 1 PHANGS-JWST Treasury program (ID 02107, PI: J. Lee).The survey consists of MIRI and NIRCam images for galaxies observed by the PHANGS-MUSE program (Emsellem et al. 2022).Key properties of our target galaxies are summarised in Table 1.
For each target we obtained one or more NIRCam and MIRI tiles, covering roughly the same area on sky and maximising the overlap with existing MUSE observations.MIRI images were obtained in the F770W, F1000W, F1130W and F2100W filters, centred at 7.7, 10, 11.3, and 21 µm, respectively.With NIRCam we obtained images in the F200W, F300M, F335M, and F360M filters, centred at 2.0, 3.0, 3.37, and 3.6 µm, respectively.NIR-Cam and MIRI observations were obtained in sequence, in order to obtain a MIRI background pointing in parallel with the NIRCam observations.However, for three galaxies (IC 5332,NGC 4303,NGC 4321) all or part of the NIRCam observations were not obtained due to guide star failures.In the case of IC 5332 NIRCam data was obtained during a subsequent visit, while for NGC 4303 and NGC 4321 only one NIRCam tile was obtained out of the intended two-tile mosaic, therefore covering roughly half of the area imaged with MIRI.Because of the failure of the first NIRCam visit for IC 5332, to process its MIRI data the background frame for NGC 7496 (taken on the same day) was used.
The observing strategy and data reduction workflow are presented in Lee et al. (2023) and Williams et al. (in preparation), respectively.Here we briefly discuss specific reduction steps of particular importance for this work, and which represent an update with respect to the early reduction process described in Lee et al. (2023).We used the latest version of the CRDS context (jwst_1070) and science calibration pipeline (v1.10.0).Secondly, we have implemented a custom background matching technique to minimise the per-pixel difference in each pair of overlapping tiles.The JWST pipeline assumes a constant background across the image, which is not the case for our observations (in many cases, one half of the image has a significantly lower 'background' level than the other).Our routine avoids this by calculating the median per-pixel difference in each overlapping pair, and finding values that minimise this across the tile set.Thirdly, we have improved the absolute astrometry for the long wavelength MIRI tiles, by basing the astrometric solution of the longer bands on that of the shortest MIRI band (F770W).As each tile in each filter is observed close in time, and using the same guide stars, this produces significantly better astrometric solutions than attempting to match each band separately.Finally, as discussed in Leroy et al. (2023b) the background level of the overall image mosaics are re-scaled to match those of existing Spitzer or WISE data.This step is necessary because we generally do not have empty sky regions that can used for background subtraction within our mosaics.
We convolved the JWST images to match the MUSE point spread function (PSF), which is different for each target as discussed in the next subsection.The required kernels were generated using the methodology of Aniano et al. (2011) and the JWST PSFs were obtained via the WebbPSF tool (Perrin et al. 2014).Error maps produced by the JWST reduction pipeline were also convolved with the same kernels and the errors are propagated analytically through the convolution process.For one galaxy, NGC 4535, the MUSE PSF full width at half maximum (FWHM) (0.56 ′′ ) is smaller than the F2100W PSF.We therefore convolve the JWST data with the smallest suitable kernel (an 'aggressive' kernel in the definition of Aniano et al. 2011) to generate an image with a Gassian PSF.The resulting PSF FWHM is of 0.71 ′′ FWHM, and the MUSE maps for this galaxy are smoothed to the same PSF.
In the NIRCam wavelength range, the PAH emission feature at 3.3 µm falls in the F335M medium-band filter.As discussed in Sandstrom et al. (2023a), we performed a bespoke continuum subtraction procedure using the adjoining F300M and F360M medium-band filters in order to remove stellar continuum emission and minimise the impact of the PAH emission contaminating the F360M side-band.Our procedure is based on and improves upon the scheme first suggested by Lai et al. (2020).In this work we refer to the resulting continuum-subtracted emission as where F335M cont is the continuum contamination estimated using equations 10 and 11 from Sandstrom et al. (2023a).For the rest of this work we only work with the continuum-subtracted F335M PAH band.
The data for a few galaxies necessitated further by-hand masking of prominent diffraction spikes (NGC 1365, NGC 1566, NGC 7496) caused by bright AGN, evident in the F2100W images presented in Fig. 1, further described below.

MUSE data and the optical H ii region catalog
We compare the JWST images with maps of Hα and Hβ obtained with IFS from the PHANGS-MUSE program (Emsellem et al. 2022), and the associated H ii region masks (Groves et al. 2023).In short, individual MUSE pointings (field of view of 1 ′ × 1 ′ ) were convolved with a kernel chosen to generate a common Gaussian PSF across each galaxy mosaic which is also constant as a function of wavelength.We refer to the resulting PSF as convolved-optimised or 'copt'.We adopt this copt FWHM as the target resolution for the JWST images, except in the case of NGC 4535 already noted above.Table 1 summarises the adopted FWHM resolution for each galaxy.
H ii region masks were derived with the hiiphot algorithm (Thilker et al. 2000), which separates nebulae from the surrounding emission by first defining seed regions around bright peaks, and extending them until a termination gradient in the surface brightness profile is reached.Line fluxes and key nebular properties (e.g.metallicity, ionisation parameter) were derived within each H ii region integrating the spectra inside the mask following the procedure detailed in Groves et al. (2023).
Fluxes in each JWST band were extracted within our H ii region masks by reprojecting the convolved JWST images onto the MUSE astrometric grid.Regions only partially covered by the footprint of the MIRI data or including masked or saturated pixels were discarded.Since the MIRI and NIRCam coverage in each target is not identical, and some targets lack NIRCam data coverage, some H ii regions have MIRI fluxes but lack NIRCam ones, or vice versa.In general, for each section of the analysis we use all H ii regions with valid JWST data, leading to a slightly different sample of regions when different bands are considered.We have checked that restricting the analysis to H ii regions that have valid fluxes in all JWST bands does not alter the results.
No local background subtraction was performed in our fiducial analysis because of the substantial crowding of sources and the resulting uncertainty intrinsic to a local background definition.However, we test the impact of this choice by performing a simple local background subtraction for each H ii region.We define the background as the mean surface brightness in an annulus around each region two MUSE pixels (each pixel is 0.2 ′′ ) distant from the region boundary and three pixels in width.If such an annulus overlaps with other regions, we remove the overlap area from the background mask.This procedure is simplistic and leads to fairly large background levels.For example, at 21µm 25% of the H ii regions in our sample have fluxes consistent with the local background with their error, and the average contrast between the flux in the H ii regions and in the background is 1.5.Because of the fairly arbitrary choices involved in the definition of a background level, we focus instead on the analysis of H ii region fluxes without background subtraction throughout this work, but we comment in the text and tables on the effect of the background subtraction, where relevant.
When considering the properties of H ii regions we further applied cuts in the BPT (Baldwin-Phillips-Terlevich, Baldwin et al. 1981;Phillips et al. 1986) Kauffmann et al. (2003), and Kewley et al. (2001), respectively, to minimise contamination from supernova remnants and planetary nebulae (Groves et al. 2023).This excluded 25% of the starting sample of nebulae.Finally, we discarded regions overlapping with foreground stars (0.3 %), or that have a signal-to-noise ratio lower than three in Hα or Hβ (0.4%).We obtain a final sample of 19 901 regions with signal-to-noise ratio larger than three in F2100W.Their distribution across the sample is summarised in Table 1.

Balmer decrement and extinction correction
The attenuation correction using the Balmer decrement is computed assuming Case B recombination, temperature T = 10 4 K and density n e = 10 2 cm −3 , leading to I Hα,corr /I Hβ,corr = 2.86 (Osterbrock & Ferland 2006).Under these assumptions where k Hα and k Hβ are the value of the assumed reddening curve at the wavelengths of the two emission lines.We assume the attenuation law of O'Donnell (1994).
In order to minimise the effect of distance (and correlated distance uncertainties) across our sample, we present results for H ii regions in terms of luminosity per unit area, and in the following I Hα , I Hα,corr and I band refer to the luminosity per unit area of observed Hα emission, Balmer decrement corrected Hα, and one of the JWST bands, all in units of erg s −1 kpc −2 .
The range of distances covered by the sample (from 5.2 to 19.6 Mpc) may affect the amount of background contamination, although within our small sample we cannot identify a distancerelated change in the relations we study in the rest of this work.Values of luminosity per unit area are corrected by a factor of cos i, where i is the galaxy inclination, to account for projection effects (Lang et al. 2020, Table 1).Since our sample excludes highly inclined galaxies, not performing this correction has minimal impact on the results.

Other ancillary data
Throughout this work we use a number of additional physical quantities to characterise H ii regions and their local environments.We summarise the definitions, data sets leveraged, and the publications where each of these physical quantities were first derived in Table 2.
In particular, we define a specific star-formation rate (sSFR) at the location of our H ii regions from the attenuation-corrected Hα assuming full sampling of the IMF, to allow for a direct comparison with the work of Belfiore et al. (2023).We define log (sSFR/yr −1 ) = C Hα I Hα,corr /Σ ⋆ , where C Hα is the conversion factor derived for a fully-sampled initial mass function (IMF) by Calzetti et al. (2007), and Σ ⋆ is the stellar mass surface density estimated from the MUSE data at the location of the H ii region, corrected downwards by 0.09 dex to bring it into agreement with the stellar masses obtained by Salim et al. (2018) and Leroy et al. (2019) based on SED fitting (see discussion in Belfiore et al. 2023).This stellar mass surface density estimate does not refer to the mass of the cluster ionising the H ii region, which cannot be resolved in the MUSE data, but is derived in larger spatial bins, of size comparable to the copt PSF FWHM (Emsellem et al. 2022).
region masks the small mismatch in resolution does not affect our results.

Using 21 µm emission as a tracer of obscured star formation
We aim to provide attenuation-corrected measures of the ionising photon luminosity from massive stars to trace recent star formation in individual H ii regions.To do so, we follow the well-established approach of using the MIR to correct for the obscured part of the Hα luminosity (Hirashita et al. 2003;Calzetti et al. 2007;Kennicutt et al. 2009;Belfiore et al. 2023).Taking the MIRI F2100W band centred at 21 µm as our fiducial band, we can write the attenuation-corrected Hα luminosity per unit area as where all luminosities are expressed in units of erg s −1 kpc −2 (i.e., I 21 = ν 21 I ν ) and a 21 is a dimensionless coefficient.This coefficient is generally calibrated empirically because of its complex dependence on star-formation history, star-dust geometry, shape of the dust attenuation law, and potential non-linearity in the dust emissivity with radiation field strength.The dust attenuation at the wavelength of Hα can therefore be related to the I 21 /I Hα ratio via A Hα = 2.5 log (1 + a 21 I 21 /I Hα ). ( 5) An equivalent formalism can be used with the PAH-tracing bands (Sec.3.2).As a starting point for our analysis we consider the numerical value of the coefficient obtained by Calzetti et al. (2007)   for a sample of star-forming complexes in nearby galaxies using Spitzer 24 µm MIPS data: a C07 24 = 0.031.To apply this to JWST data we take into account the different filter bandpasses between MIPS 24 µm and MIRI F2100W using the average flux ratio observed across our sample and tabulated by Leroy et al. (2023b), ⟨I 21 /I 24 ⟩ = 0.91, leading to a C07 21 = 0.034.

21 µm emission in Hii regions and the diffuse component
In this section we present the relative spatial distribution of 21 µm and Hα emission.In Fig. 1 we show three-colour images combining Hα from MUSE (red) with 21 µm emission from MIRI at copt resolution (green) and stellar mass surface density from MUSE full-spectral fitting (blue).Galaxies are ordered by increasing total stellar mass.The figure highlights the variety of stellar and ISM morphologies present in our sample.Throughout the sample we observe a strong spatial correlation between bright 21 µm sources and bright Hα emission, as expected if both components are powered by young stars (and dust attenuation is moderate, see discussion in Leroy et al. 2023b).The fainter 21 µm emission outside H ii regions appears filamentary, suggesting that such emission traces primarily changes in the column density of cold gas illuminated by the diffuse interstellar radiation field (Leroy et al. 2023b;Sandstrom et al. 2023b).The presence of diffuse Hα, on the other hand, correlates spatially with the positions of H ii region, and can largely be attributed to leaking ionising radiation (Belfiore et al. 2022).
To discuss the relation between Hα and 21 µm emission more quantitatively we focus on one of the galaxies in our sample, NGC 628, a grand-design spiral galaxy with log(M ⋆ /M ⊙ ) = 10.3.Fig. 2 (top row) presents a comparison between the 21 µm and Hα maps for NGC 628 at matched (copt) resolution.To quantify the spatial differences between Hα and 21 µm emission we compute the 21 µm flux within and outside the optically defined H ii region masks (blue contours in Fig. 2, bottom left).We find that in NGC 628, 54% of the 21 µm emission, but only 36% of the Hα emission, lies outside the H ii region masks.Across the full sample the fraction of 21 µm flux outside H ii regions is on average ∼ 60%.This fraction is higher than the fraction of diffuse Hα, which averages to 38%.
The large amount of 21 µm emission outside H ii regions is not due to a population of fully embedded star-forming regions, which may have eluded our optical catalogs.To demonstrate this, we match the MUSE H ii region catalog with the catalog of 21 µm compact point sources from Hassani et al. (2023) for the four galaxies in common (IC 5332, NGC 628, NGC 1365, NGC 7496).For this sub-sample we identify 21 µm sources which do not overlap MUSE H ii regions.We allow a maximum overlap of 10% in area and require sources to have MIR colours typical of the ISM, in order to exclude background galaxies and dusty stars (the selection criteria are discussed in Hassani et al. 2023).While the vast majority of 21 µm point sources overlap with H ii regions, in NGC 628 we find 53 non-overlapping sources, accounting for 0.7% of the total 21 µm emission (shown as red circles in Fig. 2, bottom left).These sources do not overlap with supernova remnants or planetary nebulae in the Groves et al. ( 2023) catalog, so we consider them likely fully embedded star-forming regions.Over the sub-sample of four galaxies in common with Hassani et al. (2023) we find that the fraction of 21 µm emission in such embedded point sources is of the order of a few percent.Their contribution is therefore negligible in explaining the diffuse 21 µm emission.
The difference in diffuse fraction between Hα and 21 µm emission implies that areas that are faint in Hα show a higher I 21 /I Hα ratio.In the bottom-right panel of Fig. 2 we show a ratio map of I 21 /I Hα with red contours corresponding to regions of bright Hα emission1 .The higher I 21 /I Hα in the DIG with respect to H ii regions does not imply, however, that the diffuse component suffers from a higher dust attenuation, as may be expected from a naive application of Eq. 5.In fact, the 21 µm and Hα emission are powered by different components of the radiation field, and the observed change in the ratio likely reflects differences in the radiation field rather than dust attenuation.The diffuse Hα is mostly the result of the leakage of ionising photons from H ii regions, with a small (few percent) contribution from old stars (Belfiore et al. 2022).The 21 µm emission, on the other hand, is powered by UV and optical light escaping star-forming regions but also originating from the old stellar population, with this ratio changing as a function of local environment (De Looze et al. 2014;Williams et al. 2019).In this work we therefore focus enitirely on the emission properties of H ii regions, where both dust and Hα should be powered by the same sources.The contamination from the diffuse component is, however, likely to affect the measured properties of the fainter H ii regions.A more detailed study of the diffuse infrared emission and its dependence on local environment will instead be presented in a future publication.For the remainder of this work we focus on calibrating the dust correction measures for H ii regions.

Hii region scaling relations at 21 µm
We test the framework for dust attenuation correction presented in Sect.3.1 on the scale of H ii regions using the JWST F2100W band.We start by presenting the relation between the observed I Hα and I 21 in H ii regions (Fig. 3).The two quantities show an excellent correlation.A linear fit to log quantities is performed by using the Bayesian framework described in Kelly (2007) which takes into account the errors on both variables.This leads to a sub-linear slope (0.92) and residual scatter around the bestfit relation of 0.27 dex.
If one corrects the Hα emission for dust attenuation via the Balmer decrement, we obtain a relation with lower scatter (0.22 dex), but still a sub-linear slope of 0.86 (Fig. 3, middle panel).In this space, however, the best-fit relation (blue line) deviates from the observed data for regions of high surface brightness.This result may be compared with Leroy et al. (2023b), who analyse most of the area in each galaxy, therefore including diffuse emission, and obtain a slope of 0.77 (inverse of their reported slope of 1.29 fitting the I Hα as a function of I 21 ).For H ii regions, at log(I Hα /erg s −1 kpc −2 ) > 39.0 the relation approaches a linear slope (black line).Since bright H ii regions also show higher values of dust attenuation, we interpret this as indicating that, when obscured SFR is dominant, both I Hα,corr and I 21 are good proxies for the total ionising photon production rate.
Finally, we show the relation between I Hα,corr and the hybrid 21 µm and Hα tracer (Eq.4) using the constant a C07 21 value.In this case the quantities on the two axes show correlated errors and we take the non-diagonal elements of the covariance matrix into account to perform the fit, as described in Kelly (2007).The effect of covariance between the x-and y-axes is small in our case, with an average correlation coefficient of 0.02.Performing a fit to the log quantities, we obtain a power law slope of 0.85 and scatter of the residuals of 0.10 dex.The hybrid tracer therefore provides an excellent match to the Balmer decrement-corrected Hα, but the resulting relation is still not linear.In particular, comparing the best-fit line (dashed blue) with the one-to-one line (black) the calibration leads to an overestimate of the dust correction at low Hα surface brightness and a slight underestimate for the high end.
A local background subtraction leads to relations that are closer to linear, but suffer from larger scatter (e.g. a slope of 1.03 for the relation between I Hα,corr and the hybrid 21 µm and Hα, and a scatter of 0.37 dex).The deviation from linearity at low Hα surface brightness is less evident, but still present.Because of the difficulty to provide a robust measure of the background, especially for faint regions, we do not further pursue this calculation.

Secondary dependencies of the calibration coefficient for 21 µm emission
In this section we test whether an additional physical parameter drives the deviation from linearity observed using a constant a 21 coefficient.We assume that the Balmer decrement traces dust attenuation across our sample (i.e. that we do not see any fully embedded emission) and calculate the a 21 coefficient using equation 5.For ease of benchmarking we consider the ratio between the inferred a 21 and the C07-equivalent value of a C07 21 .Fig. 4 shows a 21 /a C07 21 as a function of several H ii region properties derived either in this work or in Groves et al. (2023).In several of the panels in Fig. 4 the quantities on the x-and y-axis are not statistically independent.We have checked, however, that analytical propagation of the error correlations implies a degree of correlation much smaller (on average 2 dex smaller) than that observed in the data: i.e. the observed trends cannot be explained by covariance among the errors alone.
We quantify the correlations between each set of quantities via the Spearman rank correlation coefficient.In order to estimate the error on this statistic we perform a Monte Carlo simulation by adding Gaussian noise to the data.We take into account the correlations between the errors on the x-and y-axis by sampling from a 2D Gaussian with the appropriate covariance matrix for each data point.The median and standard deviations of 1000 Monte Carlo runs are taken as our best estimates of the Spearman rank correlation coefficient and its error, and are shown in the top-left for each panel in Fig. 4. Given the relatively high signal-to-noise ratio in H ii regions and the large dynamic range spanned in most quantities with respect to their errors, the errors on the correlation coefficients are found to be small (∼ 0.01 in all panels).
The inferred a 21 coefficient demonstrates significant positive correlations with EW(Hα) (ρ = 0.65, panel a), sSFR estimated from the attenuation-corrected Hα (ρ = 0.63, panel b), attenuation-corrected Hα surface brightness (I Hα,corr , ρ = 0.60, panel c), E(B−V) computed from the Balmer decrement (ρ = 0.54, panel d), and the ratio of 21 µm to 2.0 µm surface brightness, where the surface brightness at 2.0 µm is taken from the NIRCam imaging data in the F200W filter (ρ = 0.39, panel e).EW(Hα), sSFR, and I 21 /I 2.0 trace the luminosity of the H ii region with respect to that of the old stars.E(B−V) is known to correlate with both attenuation-corrected Hα (e.g.Emsellem et al. 2022), and EW(Hα) (Groves et al. 2023).Overall, these trends suggest that a 21 is slightly higher than the Calzetti et al. (2007) value for the brightest (and dustiest H ii regions), while it should be significantly lower, by up to 1.0 dex for faint, less dusty regions.All these correlations flatten for the brightest H ii regions.
Substantially weaker correlations (with |ρ| < 0.3) are found as a function of galactrocentric radius (panel f), stellar mass surface density (panel g), molecular gas surface density (panel h), gas-phase metallicity (panel i), or ionisation parameter (the ratio between the ionising photon flux and the hydrogen density in H ii regions, panel j).
Finally, the relation between a 21 and [S ii]λλ6717,31/Hα (panel k) or [O i]λ6300/Hα (panel l) is flat at low values of these line ratios, but shows a strong negative correlation (ρ = −0.46 and ρ = −0.45,respectively) for higher ones.These line ratios have been studied in our sample by Belfiore et al. (2022) and are found to increase with decreasing Hα surface brightness, corresponding to the transition between H ii regions and diffuse ionised gas (Haffner et al. 2009;Zhang et al. 2017).Even though we are focusing on H ii regions here, we are likely seeing the increased contribution of the diffuse ionised gas background to the  2007) for H ii regions as a function of several physical properties.Grey points are individual H ii regions, while red hexagons represent the median values as a function of the quantity on the x-axis.In each panel we report the values of the Spearman rank correlation coefficient (and its error).The 16 th -84 th percentiles are shown as a red shaded area.The quantities considered on the x-axis are a) equivalent width of Hα in emission, b) sSFR surface density, defined as the ratio between the SFR surface density obtained from dust attenuation corrected Hα (assuming a fully sampled IMF) and the stellar mass surface density measured from MUSE data via spectral fitting (Σ ⋆ ), c) Hα luminosity per unit area corrected for dust attenuation with the Balmer decrement (I Hα,corr ), d) E(B − V) of the gas obtained from the Balmer decrement, e) ratio of the luminosity per unit area in F2100W to that in F200M (I 21 /I 2.0 ), f) luminosity per unit area in F2100W, g) MUSE stellar mass surface density from fits to the stellar population to optical spectra, h) molecular gas surface density, derived from ALMA CO(2-1), i) gas-phase metallicity, j) ionisation parameter, k) [S ii]λλ6717,31/Hα ratio, and l) [O i]λ6300/Hα ratio.
flux within the masks corresponding to faint H ii regions going hand in hand with an increase in the relative significance of the 21 µm background.
We have tested the effect of defining a local background subtraction, which leads to an increase in the scatter and a decrease in the strength of all correlation coefficients, but does not change the shape and sign of the strongest ones.

Using PAH-tracing bands to measure obscured star formation
In this section we test the ability of PAH-dominated bands to trace the obscured component of star formation.We consider the continuum-subtracted F335M NIRCam band (F335M PAH ), and the MIRI bands centred at 7.7 µm (F770W), 10 µm (F1000W), and 11.3 µm (F1130W).Emission in the 7.7 µm and 11.3 µm bands is dominated by features generally associated with PAHs (Smith et al. 2007;Draine & Li 2007), while the star-light contamination in the F335M band is removed as described in Sec.
2.1 and Sandstrom et al. (2023b).The nature of the continuum emission underlying the PAH features and dominating the 10 µm band emission remains unclear (Smith et al. 2007;Li 2020).Leroy et al. (2023b) find that the 10 µm continuum closely follows the PAH-dominated bands, rather than the 21 µm emission.In this work we therefore refer collectively to F335M PAH , F770W, F1000W, and F1130W as 'PAH-tracing' bands.Table 3. Median (and 16 th and 84 th percentiles) of the ratio between the F335M PAH , F770W, F1000W, F1130W to the F2100W luminosity (in units of erg s −1 kpc −2 ) for the H ii regions in our sample.The dependence of the I band /I 21 on EW(Hα) is shown in Fig. 7.

PAH-tracing band ratios in Hii regions
We adopt an empirical strategy to set fiducial values of the hybridisation coefficients for each band, a band .We calculate these coefficients, which we refer to as 'C07-equivalent', by scaling the value of the Calzetti et al. (2007) 21 µm coefficient using the average band ratio between PAH-tracing bands and F2100W.In detail, we define log (a C07 band ) ≡ log (a C07 21 ) − ⟨log (I band /I 21 )⟩, (6) where ⟨log (I band /I 21 )⟩ is the median band ratio between the bands of interest (3.3, 7.7, 10, and 11.3 µm) and 21 µm in our H ii region sample.We report the median, 16 th , and 84 th percentiles of the distribution of I band /I 21 in Table 3.The scatter in the I band /I 21 lies in the range of 0.2 dex (F770W) to 0.4 dex (F335M PAH ).Part of this larger scatter for F335M PAH may be attributed to residuals in the continuum subtraction.

Hii region scaling relations for PAH-tracing bands
Fig. 5 shows the relations between Hα corrected using the Balmer decrement and the Hα hybridised with each of the PAHdominated bands using the a C07 band factor defined in Eq. 6.For all bands considered we obtain a best-fit power law with sub-linear slopes (ranging from 0.82 to 0.89) and ∼0.10 dex scatter (except for F335M PAH for which we obtain a marginally larger scatter of 0.13 dex).
These results are remarkably consistent with those obtained with F2100W in Sec.3.1.2.In fact, Fig. 5 demonstrates that the hybrid recipe overestimates the dust correction at low surface brightness levels (log(I Hα,corr /erg s −1 kpc −2 ) < 39), while underestimating it at high surface brightness levels, as already observed for F2100W.

Secondary relations of a band for PAH-tracing bands
In Fig. 6 we show the median trends of the a band coefficients for different bands as a function of EW(Hα).The values of a band are scaled logarithmically and normalised, so that the zero on the y-axis corresponds to the C07-equivalent value.The trend for F2100W (purple) was already presented in Fig. 4, but here we compare it directly with the behaviour of the other bands.For log EW(Hα)/Å < 1.5 (or log(sSFR/yr −1 ) < −10) all bands follow a closely matching trend of increasing a band with EW(Hα), reaching the C07-equivalent value at log EW(Hα)/Å ∼ 1.5.
At high EW(Hα), on the other hand, the behaviour of different bands diverges.The dependence of a 21 on EW(Hα) flattens for log EW(Hα)/Å > 2, plateauing to a constant value of a 21 that is 0.16 dex higher than the C07-equivalent one.The F335M PAH band also shows a flattening in the slope of the relation but with no plateau, while F770W, F1000W, and F1130W follow increasing trends.
The difference between the PAH-tracing bands reflects the relative brightness of different features in the MIR spectrum.For example, Chastenet et al. (2023b) find that the ratio of (I 7.7 + I 11.3 )/I 21 , which traces PAH abundance (Draine et al. 2021), is relatively constant in the diffuse ISM, but decreases in bright H ii regions.Egorov et al. (2023) find that this ratio in H ii regions anti-correlates with the ionisation parameter, which is tracing the density of the extreme UV photons in the region.Moreover, Chastenet et al. (2023a) find that the F335M PAH band is enhanced and the F1130W band is suppressed in regions of high Hα surface brightness.These trends indicate that the abundance of PAHs is lowered in H ii regions, due to their destruction by extreme UV photons, but also imply changes in the average properties of the PAH population, namely a decrease in average size and an increased abundance of charged grains and/or hotter grains.
We show the change in the I band /I 21 ratios as a function of EW(Hα) in Fig. 7.The ratio for each band is normalised by subtracting the median value for the sample, allowing for easier comparison of the trend across bands.Fig. 7 shows the expected decrease in the ratios of I 7.7 /I 21 , I 10 /I 21 , and I 11.3 /I 21 as a func- Finally, our targets are relatively metal-rich, so we do not probe the low-metallicity regime typical of dwarf galaxies where the PAH fraction is observed to be lower (Madden et al. 2006; EW(Hα) remains the most significant secondary correlation for the PAH bands, and we therefore do not consider additional ones in this work.

Discussion
This work presents an unprecedented set of resolved IR data for star-forming regions at ∼ 100 pc scales.Our results show the promise, and caveats, involved in the use of hybrid SFR indicators on the scales of H ii regions.In particular, we use IR observations with JWST at 3.3,7.7,10,11.3,and 21 µm to correct Hα emission from H ii regions for dust attenuation using a hybrid recipe.Using a scaled version of the C07 coefficient (Eq.6) leads to hybrid SFR estimations which agree with the SFR derived via the Balmer decrement to within ∼0.1 dex, but show systematic deviations for both low and high SFR.These deviations correlate most notably with quantities tracing sSFR (e.g., EW(Hα)).

Recommendation for using constant a band
We report in Table 4 the median values of the a band coefficients for our full sample of 19 901 H ii regions, together with the C07equivalent values.Our median a band are in excellent agreement with the C07-equivalent values, but the scatter around these values (quantified as the 16 th and 84 th percentiles) are large.
In Table 4 we also present the median values of the a band coefficients after performing a local background subtraction.These values are ∼ 0.2 to 0.4 dex larger than our nominal ones.The background subtraction more strongly affects the F770W, The agreement between our values and the C07-equivalent ones is not immediately expected.Calzetti et al. (2007) estimate the background in large sub-galactic regions, but not using local annuli, therefore likely resulting in an intermediate estimate of the value of the background level.However, the median a band will also depend on the distribution of EW(Hα) in the sample of regions considered and on the median spatial resolution.
In fact, Fig. 6 demonstrates that the plateau observed using kpc-resolution IR data (WISE 22µm in Belfiore et al. 2023) is a resolution effect, due to the blending of different H ii regions with each other and with the diffuse medium.The trend observed by Belfiore et al. (2023) using kpc, 15 ′′ -resolution WISE W4 22 µm data and the full PHANGS-MUSE galaxy sample, is shown as the black dashed line in Fig. 6.JWST F2100W observations, providing a factor of ∼ 15 higher spatial resolution than WISE W4, reveal that the relation between a 21 and EW(Hα) plateaus at higher EW(Hα) (and a 21 ) than seen in the low-resolution data.Since our observations resolve the inter-cloud distance (∼ 100 pc, Kim et al. 2022), and the contamination from the diffuse background is minimal for the bright regions, we expect the trend uncovered by JWST to persist at even higher spatial resolution.This hypothesis can be tested with upcoming JWST observations of more nearby (or Local Group) galaxies.
Finally, we note that our background subtraction strategy does not remove the trend of increasing a band with EW(Hα).Future work may consider more advanced approaches to model the background and test whether such contamination can explain the behaviour observed in Fig. 6 at low EW(Hα).

Modelling the impact of old stellar population on a band
The trend of increasing a band with EW(Hα) likely reflects the increasing contribution from dust heated by old stellar popula-tions in fainter regions (Cortese et al. 2008;Leroy et al. 2012;Boquien et al. 2016;Nersesian et al. 2019), and, for the PAHtracing bands, the roughly constant band ratios in diffuse regions (Chastenet et al. 2019).
Modelling the behaviour of a band as a function of EW(Hα), or one of the other parameters studied in Sec.3.1.3,would provide a more accurate estimate of the dust-attenuated star formation, and a prescription more readily applicable to lowerresolution data.Several authors have presented calibrations including such additional terms, for example colours, sSFR, or M ⋆ (Boquien et al. 2016;Zhang & Ho 2023;Belfiore et al. 2023;Kouroumpatzakis et al. 2023).
We aim to provide such a second-order correction applicable to nearby galaxies data.The functional form that this correction should take is not evident a priori.For the sake of simplicity, here we add an additional linear term to the conventional hybrid star formation recipe presented in Eq. 4, leading to where I IR is one of the JWST IR dust-tracing bands, and I Q is a second physical quantity appropriately chosen to best reproduce the extinction correction from the Balmer decrement.As before, we express all surface brightness measurements in units of erg s −1 kpc −2 .We first consider the task of deriving an estimate for the obscured SFR using two JWST bands: a dust-tracing band (F2100W or any of the PAH-tracing bands), and a second band tracing stellar continuum emission (F200W, F300M or F360M) to account for the dust heating from old stars.To determine the best combination of bands for such a calibration we run a random forest regression, using the algorithm implementation in scikit-learn (Pedregosa et al. 2011).In particular, we use the random forest model to determine the relative importance of input features (Bluck et al. 2022;Baker et al. 2022).This approach is conceptually similar to more traditional principal component analysis, but allows the generation of more complex, non-linear models.
The random forest algorithm is a form of supervised machine learning that builds a non-parametric model using multiple decision trees, trained to decrease the mean square error at each split.We use E(B − V) as the target variable since the goal is to understand how best to perform the dust correction, and I band /I Hα , where 'band' is one of the JWST bands available to us (F200W, F300M, F360M, F335M PAH , F770W, F1000W, F1130W, F2100W), as input features.We also consider an array of random numbers as an additional input feature to test the performance of the algorithm.The dataset is subdivided into a training (80% of the full sample) and testing (remaining 20%) subset of H ii regions.Only H ii regions covered by both NIR-Cam and MIRI and with signal-to-noise ratio grater than 3 in all JWST bands are used for this computation.After training, the mean square error (MSE) of the test set is compared with that of the training sample in order to check the quality of the resulting model and avoid overfitting.
We find that the combination of I 21 /I Hα and I 2.0 /I Hα accounts for nearly all the variance in E(B−V) (Fig. 8).Once these two variables are taken into account, the remaining ones have residual importance consistent with random within the errors.The MSE of the training and test set are comparable (reported on the top right of Fig. 8), demonstrating that the algorithm produced a successful model.This analysis confirms two of the key insights from this work: 1) F2100W is the best band to use for hybridisation among the bands tracing dust emission available in our PHANGS-JWST filter set, and 2) there exists an evident secondary dependence on the stellar mass-tracing NIRCam bands, of which the best is F200W.
Having established I 21 and I 2.0 as the two best variables, we fit the extinction-corrected Hα surface brightness with a twoparameter regression model, using the Bayesian formalism discussed in Hogg et al. (2010).We fit for the α and β parameters in Eq. 7 and for intrinsic scatter in the relation, using the Bayesian Monte Carlo sampler emcee, and obtain best-fit as relation We also fit using the F300M and F360M NIRCam bands, obtaining respectively α = 0.025, β = 8 × 10 −4 , and α = 0.025, β = 7 × 10 −4 .
Our approach here is similar to the one of Kouroumpatzakis et al. ( 2023), who use spectral energy distribution models generated with the code cigale to study how best to determine SFR from JWST bands.In their work, however, Kouroumpatzakis et al. (2023) do not consider I Hα , but model the SFR as a linear combination of I 21 and I 2.0 .Their results are therefore not directly comparable.We leave a direct comparison of our results with spectral energy distribution models to future work.
Finally, we repeat our two-parameter fit using stellar mass surface density Σ ⋆ , which may be used in the absence of JWST NIRCam data.We obtain where surface brightnesses are in units of erg −1 s −1 kpc −2 and Σ ⋆ is in units of M ⊙ kpc −2 .

Does the IR emission trace embedded SFR or cold gas?
The MIR, and in particular the bands containing strong PAH features, show excellent correlations with the cold molecular gas content, mapped by CO emission (Regan et al. 2006;Chown et al. 2021;Whitcomb et al. 2023;Leroy et al. 2023a).Leroy et al. (2021b) and Leroy et al. (2023a), for example, found a nearly linear relation between the intensity in the WISE W3 band centred at 12 µm and CO(2-1) emission.This correlation is found to have lower scatter than that between CO and IR emission at ∼ 24 µm.Whitcomb et al. (2023) similarly argued using correlation analysis that emission around 12 µm traces CO emission better than SFR, while IR emission at 24 µm traces SFR better than CO.Leroy et al. (2023b) use the early JWST data for four galaxies from our sample and compare the MIR emission with both Hα and CO emission, finding that MIR emission correlates well with both extinction-corrected Hα and CO.Such a result is expected, since dust emission depends on both the dust surface density and the intensity of the radiation field heating it.The convolution of these effects leads to MIR emission being a useful tracer of both I Hα,corr and CO.Within H ii regions, however, one expects a tighter relation between MIR emission and star formation than between MIR emission and cold gas content.To test this hypothesis we evaluate the relative ability of each of the JWST bands considered in this work to trace extinction-corrected Hα, molecular gas (as traced by CO(2-1)), and stellar mass surface density, focusing our analysis on H ii regions specifically.The analysis is performed by building a random forest regression model for the surface brightness in each JWST band and using I Hα,corr , I CO , and Σ ⋆ , in addition to a vector of random numbers, as input features.The relative importance of each of the features in reproducing the brightness in each JWST band considered is presented in Fig. 9.
As expected, both F300M and F200W trace mostly stellar mass surface density (importance ∼ 90%).For F360M Hα takes up part of the importance, probably reflecting the increasing importance of hot dust and residual contamination from PAHs in this band (Sandstrom et al. 2023a).F335M PAH mostly traces Hα, highlighting the success of the continuum subtraction strategy for this band.The MIR bands trace mostly Hα, but with 20% importance attributed also to CO for F770W, F1000W, F1130W.The contribution of stellar mass is small, and it is largest at F1000W, where is amounts to 5% of the overall importance.This band contains the largest contamination from point sources, potentially evolved stars, whose density correlates with the overall stellar mass surface density.For the canonical F2100W band the relative importance of Hα and CO are 83% and 16% respectively.In summary, both 3.3 µm and the MIR bands trace mostly extinction-corrected Hα within H ii regions, with a smaller (<20%) secondary dependence on molecular gas, as traced by CO.The small dependence on CO is consistent with the finding that many H ii regions exists without associated molecular gas, in the sense of gas within the H ii region boundary (Zakardjian et al. 2023).These conclusions are likely to differ within the more diffuse medium, where the relative importance of CO in determining the MIR emission is expected to increase (Leroy et al. 2023b), but we leave an exploration of these trends to future work.

Conclusions
We use NIRCam and MIRI images in combination with IFS from MUSE for 19 nearby galaxies to calibrate hybrid recipes to correct Hα for dust extinction on the scales of individual H ii regions (∼ 100 pc).We examine hybrid recipes that consider a linear combination of Hα with the F2100W MIRI band, dominated by emission from small grains, and bands tracing PAH emission, namely the continuum-subtracted NIRCam F335M band (F335M PAH ), and the F770W, F1000W, and F1130W MIRI bands.We summarise our main results below.
-The ratio between 21 µm and Hα emission is systematically higher in the diffuse ISM, outside optically defined H ii regions.This fact reflects a change in the dominant source of dust heating -young stars in H ii regions and the old stellar population in the DIG -and not an increase in dust attenuation.Fully embedded regions make a negligible contribution to the diffuse 21 µm emission.-Focusing on H ii regions, we calibrate the coefficient in front of the IR emission term (I Hα,corr = I Hα +a band I band ) by benchmarking the hybrid dust correction recipe with the correction obtained using the Balmer decrement measured from the MUSE data.We adopt fiducial hybridisation coefficients for the different dust-tracing JWST bands (a C07 band ) by rescaling the coefficient derived by Calzetti et al. (2007) for Spitzer 24 µm emission by the average band ratio ⟨I 24 /I band ⟩.These C07-equivalent coefficients are in good general agreement with the median of the measured a band coefficient across our sample of H ii regions for all bands considered (F335M PAH , F770W, F1000W, F1130W, F2100W).A local background subtraction causes an increase of the cofficients associated with F770W, F1000W, and F1130W bands of ∼ 0.4 dex, while it has a smaller effect on F2100W and F335M PAH (a increases by ∼ 0.2 dex).Our recommended coefficients are summarised in Table 4.
-For both 21 µm and the PAH-tracing bands, using a hybrid recipe with a constant coefficient a C07 band leads to extinctioncorrected Hα fluxes that correlate well (scatter ∼ 0.1 dex), albeit sub-linearly (slope 0.82-0.89),with those inferred from the Balmer decrement.The sub-linear slope implies that a constant a band overestimates the dust correction for faint H ii regions, and it underestimates it for bright regions.
-We tested for secondary dependencies of a band on a number of additional physical parameters.The strongest correlations are with physical parameters that trace the ratio between young and old stellar components, e.g.EW(Hα), sSFR, or indirectly correlate with them (e.g.Hα surface brightness, or dust attenuation measured from the Balmer decrement).
Strong correlations are also found with line ratios that trace the change in the ionisation condition between H ii regions and DIG ([O i]/Hα, [S ii]/Hα).These trends indicate that dust heating from the old stellar component contributes substantially to the IR emission for the fainter H ii regions in our sample.a band correlates with EW(Hα) for all bands considered.For EW(Hα) < 30 Å, the trends are similar for all bands and agree with the recent kpc-resolution study of Belfiore et al. (2023).At higher EW(Hα), a 21 plateaus to a constant value, while the coefficients for the F770W, F1000W, and F1130W keep increasing.The behaviour of the MIR PAH-tracing bands reflects the relative decrease in PAH emission in bright H ii regions, in agreement with the literature predicting the destruction of PAH molecules in such harsh environments.The F335M PAH shows an intermediate behaviour between F2100W and the other MIR bands, which can be explained by the relative increase in 3.3 µm emission with respect to the other PAH-tracing bands in H ii regions, possibly tracing a decrease in the average size of PAH grains or an increase in the hardness of the UV radiation field.-A random forest analysis confirms that F2100W is the preferred band to consider in combination with Hα to infer dust attenuation, confirming results from previous literature using WISE 22 µm or Spitzer 24 µm bands.-The main secondary dependence of a band can be modelled by explicitly adding an additional linear term proportional to stellar mass surface density or 2 µm emission to the hybrid recipe.We provide coefficients for these corrections in Eq. 8 and 9.
-Finally, we analyse the relative importance of molecular gas (traced by CO(2-1)), stellar mass surface density, and extinction-corrected Hα in determining the surface brightness for H ii regions within our JWST bands.The MIR bands trace mostly Hα but with a smaller (15-20%) importance attributed to CO.

Fig. 1 .
Fig. 1.Three-colour images for the galaxies in our sample showing a combination of Hα emission from MUSE IFS (red), 21 µm emission from JWST MIRI F2100W images (green), and stellar mass surface density from MUSE full spectral fitting (blue).The JWST data and MUSE data are shown at matched resolution.Galaxies are shown in order of increasing stellar mass.Foreground stars and the AGN in NGC 1566 and NGC 1365 are masked.The diffraction spikes due to bright AGN in NGC 4535, NGC 7496 and NGC 1365 are shown here but masked in subsequent analysis.The jagged edges of the images are caused by the intersection between the MUSE and MIRI image coverage.

Fig. 2 .
Fig. 2. Comparison of the Hα and 21µm emission for NGC 0628.Top-left, 21µm surface brightness image obtained using the F2100W filter on MIRI, convolved to the MUSE copt resolution (0.92 ′′ ).Top-right, map of the Hα surface brightness from MUSE (Emsellem et al. 2022).White contours enclose regions of bright 21 µm emission (I 21 > 10 40.5 erg s −1 kpc −2 ).Bottom-left: Same as top-left, but with the boundaries of the optically selected H ii regions from Groves et al. (2023) in blue and the 21 µm sources which do not overlap with optical H ii regions from Hassani et al. (2023) in red.Bottom-right: Ratio map of the Hα to F2100W surface brightness.Red contours are bright Hα regions, I Hα > 10 38.8 erg s −1 kpc −2 .White regions within the mapped area represent areas with S/N<3 in F2100W.The boundaries of the MUSE (solid black) and MIRI (dashed black) mosaics are shown to aid in visualising the data overlap.

Fig. 3 .
Fig. 3. Left: I 21 as a function of I Hα for the H ii regions (grey points) in our sample.The white dots represent the median relation, while the blue dashed line is the best-fit power law obtained by fitting the data considering the errors associated with the quantities on both axes.The solid black line shows a linear relation normalised to match the data at high surface brightness.The slope and intercept of the best-fit power law together with the scatter with respect to the best model are reported in the top left corner.The meaning of the symbols and lines is the same across the three panels, except that the solid black line the right panel represents a one-to-one relation.Middle: I 21 as a function of I Hα,corr .Right: I Hα + a C07 21 I 21 as a function of I Hα,corr .

Fig. 4 .
Fig.4.The measured coefficient a 21 normalised to the (rescaled) value fromCalzetti et al. (2007) for H ii regions as a function of several physical properties.Grey points are individual H ii regions, while red hexagons represent the median values as a function of the quantity on the x-axis.In each panel we report the values of the Spearman rank correlation coefficient (and its error).The 16 th -84 th percentiles are shown as a red shaded area.The quantities considered on the x-axis are a) equivalent width of Hα in emission, b) sSFR surface density, defined as the ratio between the SFR surface density obtained from dust attenuation corrected Hα (assuming a fully sampled IMF) and the stellar mass surface density measured from MUSE data via spectral fitting (Σ ⋆ ), c) Hα luminosity per unit area corrected for dust attenuation with the Balmer decrement (I Hα,corr ), d) E(B − V) of the gas obtained from the Balmer decrement, e) ratio of the luminosity per unit area in F2100W to that in F200M (I 21 /I 2.0 ), f) luminosity per unit area in F2100W, g) MUSE stellar mass surface density from fits to the stellar population to optical spectra, h) molecular gas surface density, derived from ALMA CO(2-1), i) gas-phase metallicity, j) ionisation parameter, k) [S ii]λλ6717,31/Hα ratio, and l) [O i]λ6300/Hα ratio.

Fig. 5 .
Fig. 5. Comparison of Hα hybridised with F335M PAH , F770W, F1000W and F1130W and the attenuation-corrected Hα obtained using the Balmer decrement.The dashed blue line is the best linear fit while the white circles represent a binned average.The slope and intercept of the best-fit power law together with the scatter with respect to the best model are reported on the top left corner.The solid black line represents the one-to-one relation.

Fig. 6 .Table 4 .
Fig. 6. log a band normalised to the C07-equivalent values as a function of EW(Hα) and sSFR (alternative x-axis on top) for H ii regions in our sample, where band = [F335M PAH (blue), F770W (orange), F1000W (green), F1130W (red), F2100W (purple)].The coloured solid lines represent the median relations, while the coloured dashed lines show the 16 th and 84 th percentiles of the distribution.The black dashed line is the best-fit to the low-resolution (kpc-scale) Hα + WISE W4 data from Belfiore et al. (2023).The grey histogram (top) shows the distribution of EW(Hα) in the H ii regions used in this work (with 16 th , 50 th , 84 th percentiles marked and labelled).The purple histogram (right) shows the distribution of log a 21 /a C07 21 for the region in our sample, with the 16 th , 50 th , 84 th percentiles marked and labelled.

Fig. 7 .
Fig. 7.The ratio of I band /I 21 normalised to the sample mean (⟨I band /I 21 ⟩) as a function of EW(Hα), for band= [F335M PAH (blue), F770W (orange), F1000W (green), F1130W (red)].The solid lines represent the median, while the dashed lines correspond to the 16 th and 84 th percentiles.The alternative x-axis on top shows the sSFR associated with each EW(Hα) value.The median trends show that I 3.3 /I 21 does not show the same decrease with EW(Hα) as the ratios of the other PAH bands.

Fig. 8 .
Fig. 8. Relative importance of different band ratios in predicting the E(B−V) of H ii regions using a random forest model.The error-bars are obtained by bootstrap random sampling.The mean squared errors (MSE) are shown for both the testing and the training set.I 21 /I Hα and I 2.0 /I Hα account for most of the variance in the data.

Fig. 9 .
Fig. 9. Relative importance of extinction-corrected Hα (I Hα,corr ), CO(2-1) surface density, and stellar mass surface density (Σ ⋆ ) in predicting the surface brightness of H ii regions in a set of JWST bands.

Table 1 .
Distances Leroy et al. (2021b)ompilation ofAnand et al. (2021), stellar mass and SFR are fromLeroy et al. (2021b), inclinations fromLang et al. (2020)where available, otherwise fromLeroy et al. (2021b), and PSF FWHM of the copt MUSE data are from Emsellem et al. (2022) ( ⋆ for this galaxy the copt FWHM is smaller than the F2100W PSF, so we used an 'aggressive' kernel to convolve the JWST data to a Gaussian PSF).θ bm gives the size of the common resolution element with no inclination correction.N HII gives the number of HII regions fromGroves et al. (2023)overlapping with the MIRI footprint, detected in F2100W, and meeting our other selection criteria (see text).The notes include information on which galaxies have partial coverage with NIRCam.CO I CO(2−1) /R 21 , a ALMALeroy et al. (2021b)E(B − V)

Table 2 .
Summary of additional data and physical quantities used in this work.
Notes.a: assuming a Milky Way conversion factor α CO