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

Francesco Belfiore, Adam K Leroy, Thomas G Williams, Ashley T Barnes, Frank Bigiel, Médéric Boquien, Yixian Cao, Jérémy Chastenet, Enrico Congiu, Daniel A Dale, et al.

► To cite this version:

Francesco Belfiore, Adam K Leroy, Thomas G Williams, Ashley T Barnes, Frank Bigiel, et al. Calibrating mid-infrared emission as a tracer of obscured star formation on HII-region scales in the era of JWST. Astronomy & Astrophysics - A&A, 2023, 678, pp.A129. 10.1051/0004-6361/202347175. hal-05018645

HAL Id: hal-05018645

https://hal.science/hal-05018645v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Calibrating mid-infrared emission as a tracer of obscured star formation on HII-region scales in the era of JWST

Francesco Belfiore¹,², Adam K. Leroy², Thomas G. Williams³, Ashley T. Barnes⁴, Frank Bigiel⁵, Médéric Boquien⁶,
Yixian Cao⁷, Jérémy Chastenet⁸, Enrico Congiu⁹, Daniel A. Dale¹⁰, Oleg V. Egorov¹¹, Cosima Eibensteiner⁵,
Eric Emsellem⁴,¹², Simon C. O. Glover¹³, Brent Groves¹⁴, Hamid Hassani¹⁸, Ralf S. Klessen¹³,¹⁵, Kathryn Kreckel¹¹,
Lukas Neumann⁵, Justus Neumann¹⁶, Miguel Querejeta¹⁷, Erik Rosolowsky¹⁸, Patricia Sanchez-Blazquez¹⁹,²⁰,
Karin Sandstrom²¹, Eva Schinnerer¹⁶, Jiayi Sun²²,²³, Jessica Sutter²¹, and Elizabeth J. Watkins¹¹

*(Affiliations can be found after the references)*

Received 13 June 2023 / Accepted 1 September 2023

## ABSTRACT

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) is 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 used images obtained with NIRCam and MIRI as part of the PHANGS–JWST survey to calibrate the IR emission at 21 µm, and the emission in the PAH-tracing bands at 3.3, 7.7, 10, and 11.3 µm as tracers of obscured star formation. We analysed ~20 000 optically selected H II regions across 19 nearby star-forming galaxies, and benchmarked their IR emission against dust attenuation measured from the Balmer decrement. We modelled the extinction-corrected H α flux as the sum of the observed H α emission and a term proportional to the IR emission, with $a_{IR}$ as the proportionality coefficient. A constant $a_{IR}$ leads to an extinction-corrected H α estimate that agrees with those obtained with the Balmer decrement with a scatter of ~0.1 dex for all bands considered. Among these bands, 21 µm emission is demonstrated to be the best tracer of dust attenuation. The PAH-tracing bands underestimate the correction for bright H II regions, since in these environments the ratio of PAH-tracing bands to 21 µm decreases, signalling destruction of the PAH molecules. For fainter H II regions, all bands suffer from an increasing contamination from the diffuse IR background. We present calibrations that take this effect into account by adding an explicit dependence on 2 µm emission or stellar mass surface density.

**Key words.** dust, extinction – galaxies: ISM – galaxies: star formation – infrared: ISM

## 1. 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 re-radiated in the infrared (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 account of star formation rates (SFRs) in galaxies (Calzetti et al. 2007; Wuyts et al. 2011a; 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 spectral energy distribution (SED), emission at 24 µm has the advantage

of tracing a hotter dust component, and therefore being less sensitive to IR cirrus.

Several authors have advocated for 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, 2009; Calzetti et al. 2007; Wuyts et al. 2011b; 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 kiloparsec-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 polycyclic aromatic hydrocarbons (PAH) features or the brightness of PAH-tracing 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
<table><caption>**Table 1.** Key properties of the PHANGS–MUSE sample of galaxies used in this work.</caption><thead><tr><th>Galaxy</th><th><i>D</i><br>[Mpc]</th><th>log(<i>M</i><sub>⋆</sub>)<br>[<i>M</i><sub>☉</sub>]</th><th>log(SFR)<br>[<i>M</i><sub>☉</sub> yr<sup>-1</sup>]</th><th><i>i</i><br>[deg]</th><th>Copt. PSF<br>[arcsec]</th><th>θ<sub>bm</sub><br>[pc]</th><th><i>N</i><sub>HII</sub></th><th>Notes</th></tr></thead><tbody><tr><td>IC 5332</td><td>9.0</td><td>9.67</td><td>-0.39</td><td>26.9</td><td>0.87</td><td>38</td><td>545</td><td></td></tr><tr><td>NGC 0628</td><td>9.8</td><td>10.34</td><td>0.24</td><td>8.9</td><td>0.92</td><td>44</td><td>1712</td><td>Grand-design spiral</td></tr><tr><td>NGC 1087</td><td>15.9</td><td>9.93</td><td>0.12</td><td>42.9</td><td>0.92</td><td>71</td><td>912</td><td></td></tr><tr><td>NGC 1300</td><td>19.0</td><td>10.62</td><td>0.07</td><td>31.8</td><td>0.89</td><td>82</td><td>944</td><td>Nuclear ring</td></tr><tr><td>NGC 1365</td><td>19.6</td><td>10.99</td><td>1.23</td><td>55.4</td><td>1.16</td><td>109</td><td>550</td><td>AGN, nuclear ring</td></tr><tr><td>NGC 1385</td><td>17.2</td><td>9.98</td><td>0.32</td><td>44.0</td><td>0.77</td><td>64</td><td>913</td><td></td></tr><tr><td>NGC 1433</td><td>18.6</td><td>10.87</td><td>0.05</td><td>28.6</td><td>0.91</td><td>82</td><td>851</td><td>Nuclear ring</td></tr><tr><td>NGC 1512</td><td>18.8</td><td>10.71</td><td>0.11</td><td>42.5</td><td>1.25</td><td>114</td><td>477</td><td>Nuclear ring</td></tr><tr><td>NGC 1566</td><td>17.7</td><td>10.78</td><td>0.66</td><td>29.5</td><td>0.80</td><td>69</td><td>1561</td><td>AGN</td></tr><tr><td>NGC 1672</td><td>19.4</td><td>10.73</td><td>0.88</td><td>42.6</td><td>0.96</td><td>90</td><td>934</td><td>AGN, nuclear ring</td></tr><tr><td>NGC 2835</td><td>12.2</td><td>10.00</td><td>0.09</td><td>41.3</td><td>1.15</td><td>68</td><td>732</td><td></td></tr><tr><td>NGC 3351</td><td>10.0</td><td>10.36</td><td>0.12</td><td>45.1</td><td>1.05</td><td>51</td><td>690</td><td>Nuclear ring</td></tr><tr><td>NGC 3627</td><td>11.3</td><td>10.83</td><td>0.58</td><td>57.3</td><td>1.05</td><td>58</td><td>970</td><td>AGN</td></tr><tr><td>NGC 4254</td><td>13.1</td><td>10.42</td><td>0.49</td><td>34.4</td><td>0.89</td><td>57</td><td>2022</td><td></td></tr><tr><td>NGC 4303</td><td>17.0</td><td>10.52</td><td>0.73</td><td>23.5</td><td>0.78</td><td>64</td><td>2043</td><td>AGN, nuclear ring, NIRCam covers half</td></tr><tr><td>NGC 4321</td><td>15.2</td><td>10.75</td><td>0.55</td><td>38.5</td><td>1.16</td><td>86</td><td>1239</td><td>Nuclear ring, NIRCam covers half</td></tr><tr><td>NGC 4535</td><td>15.8</td><td>10.53</td><td>0.33</td><td>44.7</td><td>0.71*</td><td>54</td><td>1097</td><td>AGN</td></tr><tr><td>NGC 5068</td><td>5.2</td><td>9.40</td><td>-0.56</td><td>35.7</td><td>1.04</td><td>26</td><td>1174</td><td></td></tr><tr><td>NGC 7496</td><td>18.7</td><td>10.00</td><td>0.35</td><td>35.9</td><td>0.89</td><td>81</td><td>535</td><td>AGN</td></tr></tbody></table>

Notes. Distances are taken from the compilation of Anand et al. (2021), stellar mass and SFR are from Leroy et al. (2021b), inclinations from Lang et al. (2020) where available, otherwise from Leroy 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). θ<sub>bm</sub> gives the size of the common resolution element with no inclination correction. N<sub>HII</sub> gives the number of H II regions from Groves 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.

in terms of angular resolution and sensitivity (Calzetti 2020; Li 2020). Moreover, PAH features remain accessible and are readily identifiable with Spitzer, and now the James Webb Space Telescope (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. 2019, 2023a; 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 url (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 Sect. 2 we present our data, paying special attention to the new JWST images, and summarise the methods for inferring the dust attention. In Sect. 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 Sect. 4. We summarise the results in Sect. 5.

## 2. Data and methods

### 2.1. PHANGS–JWST data

We analysed 19 galaxies observed as part of the cycle 1 PHANGS–JWST Treasury programme (ID 02107, PI: J. Lee). The survey consists of MIRI and NIRCam images for galaxies observed by the PHANGS–MUSE programme (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. NIRCam 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 prep.), 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 be 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 Gaussian 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

$$F335M_{\text{PAH}} = F335M - F335M_{\text{cont}}, \quad (1)$$

where $F335M_{\text{cont}}$ is the continuum contamination estimated using Eqs. (10) and (11) from Sandstrom et al. (2023a). For the rest of this work we only work with the continuum-subtracted $F335M_{\text{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. This is described further below.

## 2.2. MUSE data and the optical HII region catalogue

We compared the JWST images with maps of $H\alpha$ and $H\beta$ obtained with IFS from the PHANGS-MUSE programme (Emsellem et al. 2022), and the associated HII 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 adopted 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.

HII 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 HII 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 HII 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 HII regions have MIRI fluxes but lack NIRCam ones, or vice versa. In general, for each section of the analysis we used all HII 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 HII 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 tested the impact of this choice by performing a simple local background subtraction for each HII region. We defined 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 overlapped with other regions, we removed 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 HII regions in our sample have fluxes consistent with the local background with their error, and the average contrast between the flux in the HII 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 HII 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 HII regions we further applied cuts in the BPT (Baldwin-Phillips-Terlevich, Baldwin et al. 1981; Phillips et al. 1986) diagrams [NII]/$H\alpha$ versus [OIII]/$H\beta$ and [SII]/$H\alpha$ versus [OIII]/$H\beta$, corresponding to the dividing lines of Kauffmann et al. (2003) and Kewley et al. (2001), respectively, to minimise contamination from supernova remnants and planetary nebulae (Groves et al. 2023). This
<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>
<table><caption>**Table 2.** Summary of additional data and physical quantities used in this work.</caption><thead><tr><th>Quantity</th><th>Definition</th><th>Derived from dataset</th><th>Reference</th></tr></thead><tbody><tr><td>Σ<sub>H2</sub></td><td>α<sub>CO</sub>I<sub>CO(2-1)</sub>/R<sub>21</sub><sup>(a)</sup></td><td>ALMA</td><td>Leroy et al. (2021b)</td></tr><tr><td>E(B-V)</td><td>Equation (2)</td><td>MUSE</td><td>Groves et al. (2023)</td></tr><tr><td>Σ<sub>★</sub></td><td>Spectral fitting</td><td>MUSE</td><td>Emsellem et al. (2022)</td></tr><tr><td>sSFR</td><td>C<sub>Hα</sub> I<sub>Hα,corr</sub>/Σ<sub>★</sub><sup>(b)</sup></td><td>MUSE</td><td>Belfiore et al. (2023)</td></tr><tr><td>EW(Hα)</td><td>Moments of spectra<sup>(c)</sup></td><td>MUSE</td><td>Groves et al. (2023)</td></tr><tr><td>12 + log(O/H)</td><td>Scal. calibration<sup>(d)</sup></td><td>MUSE</td><td>Groves et al. (2023)</td></tr><tr><td>log(U)</td><td>From [SIII]λλ9069, 9532 / [SII]λλ6717, 31<sup>(e)</sup></td><td>MUSE</td><td>Groves et al. (2023)</td></tr></tbody></table>

**Notes.** (a) Assuming a Milky Way conversion factor $α_{CO} = 4.35 M_{\odot} pc^{-2}$ (K km s$^{-1}$)$^{-1}$ (Bolatto et al. 2013) and a CO(2–1)-to-CO(1–0) ratio $R_{21} = 0.65$ (Den Brok et al. 2021). (b) $log C_{Hα}/(M_{\odot} yr^{-1}/(erg s^{-1})) = -41.26$, from Calzetti et al. (2007). (c) Definition in Westfall et al. (2019). (d) Scal. calibration from Pilyugin & Grebel (2016). (e) Ionisation parameter calculated using the calibration from Diaz et al. (1991).

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

## 2.3. Balmer decrement and extinction correction

The attenuation correction using the Balmer decrement was 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

$$E(B-V) = \frac{2.5}{k_{H\beta} - k_{H\alpha}} \log_{10} \left[ \frac{L_{H\alpha}/L_{H\beta}}{2.86} \right], \quad (2)$$

where $k_{Hα}$ and $k_{Hβ}$ are the value of the assumed reddening curve at the wavelengths of the two emission lines. We assumed 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 distance-related 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.

## 2.4. 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 (\text{sSFR}/\text{yr}^{-1}) = C_{H\alpha} I_{H\alpha,\text{corr}}/\Sigma_{\star}, \quad (3)$$

where $C_{Hα}$ is the conversion factor derived for a fully sampled initial mass function (IMF) by Calzetti et al. (2007), and $\Sigma_{\star}$ 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).

In addition, in Sect. 3.1.3 we use the molecular gas mass surface density, as traced by CO(2–1) maps obtained as part of the PHANGS–ALMA programme (Leroy et al. 2021a,b), at the positions of our optically defined H II regions. The ALMA data has slightly lower angular resolution than copt (~1.3′′) for each of our galaxy targets, but considering the integration over the H II region masks the small mismatch in resolution does not affect our results.

# 3. Results

## 3.1. 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

$$I_{H\alpha,\text{corr}}^{\text{hybrid}} = I_{H\alpha} + a_{21} I_{21}, \quad (4)$$

where all luminosities are expressed in units of erg s$^{-1}$ kpc$^{-2}$ (i.e. $I_{21} = v_{21}I_v$) 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\alpha$ can therefore be related to the $I_{21}/I_{H\alpha}$ ratio via

$$A_{H\alpha} = 2.5 \log(1 + a_{21} I_{21}/I_{H\alpha}). \qquad (5)$$

An equivalent formalism can be used with the PAH-tracing bands (Sect. 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_{24}^{C07} = 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 tabulated by Leroy et al. (2023b), $\langle I_{21}/I_{24} \rangle = 0.91$, leading to $a_{21}^{C07} = 0.034$.

### 3.1.1. 21 µm emission in HII regions and the diffuse component

In this section we present the relative spatial distribution of 21 µm and $H\alpha$ emission. In Fig. 1 we show three-colour images combining $H\alpha$ 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\alpha$ 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 HII 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\alpha$, on the other hand, correlates spatially with the positions of HII region, and can largely be attributed to leaking ionising radiation (Belfiore et al. 2022).

To discuss the relation between $H\alpha$ 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_{\star}/M_{\odot}) = 10.3$. Figure 2 (top row) presents a comparison between the 21 µm and $H\alpha$ maps for NGC 628 at matched (copt) resolution. To quantify the spatial differences between $H\alpha$ and 21 µm emission we compute the 21 µm flux within and outside the optically defined HII 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\alpha$ emission, lies outside the HII region masks. Across the full sample the fraction of 21 µm flux outside HII regions is on average ~60%. This fraction is higher than the fraction of diffuse $H\alpha$, which averages to 38%.

The large amount of 21 µm emission outside HII regions is not due to a population of fully embedded star-forming regions, which may have eluded our optical catalogues. To demonstrate this, we match the MUSE HII region catalogue with the catalogue 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 HII 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 HII 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) catalogue, 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\alpha$ and 21 µm emission implies that areas that are faint in $H\alpha$ show a higher $I_{21}/I_{H\alpha}$ ratio. In the bottom-right panel of Fig. 2 we show a ratio map of $I_{21}/I_{H\alpha}$ with red contours corresponding to regions of bright $H\alpha$ emission<fnref n="1" />. The higher $I_{21}/I_{H\alpha}$ in the DIG with respect to HII 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\alpha$ 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\alpha$ is mostly the result of the leakage of ionising photons from HII 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 entirely on the emission properties of HII regions, where both dust and $H\alpha$ should be powered by the same sources. The contamination from the diffuse component is, however, likely to affect the measured properties of the fainter HII regions. A more detailed study of the diffuse IR 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 HII regions.

### 3.1.2. HII region scaling relations at 21 µm

We test the framework for dust attenuation correction presented in Sect. 3.1 on the scale of HII regions using the JWST F2100W band. We start by presenting the relation between the observed $I_{H\alpha}$ and $I_{21}$ in HII 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 best-fit relation of 0.27 dex.

If one corrects the $H\alpha$ 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\alpha}$ as a function of $I_{21}$). For HII regions, at $log(I_{H\alpha}/erg s^{-1} kpc^{-2}) > 39.0$ the relation approaches a linear slope (black line). Since bright HII regions

<footnote n="1">¹ $log(I_{H\alpha}/erg s^{-1} kpc^{-2}) > 38.8$, chosen to correspond to the 16th percentile of the surface brightness distribution of HII regions in the PHANGS-MUSE sample from Belfiore et al. (2022).</footnote>
<figure><img src="image_3.png" /></figure>

**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<sub>21</sub> > 10<sup>40.5</sup> erg s<sup>-1</sup> kpc<sup>-2</sup>). 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<sub>Hα</sub> > 10<sup>38.8</sup> erg s<sup>-1</sup> kpc<sup>-2</sup>. 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.

also show higher values of dust attenuation, we interpret this as
indicating that, when obscured SFR is dominant, both I<sub>Hα,corr</sub> and
I<sub>21</sub> are good proxies for the total ionising photon production rate.

Finally, we show the relation between *I*<sub>Hα,corr</sub> and the hybrid 21 µm and Hα tracer (Eq. (4)) using the constant *a*<sub>21</sub><sup>C07</sup> 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,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Scaling relations between $I_{21}$ and $I_{H\alpha}$ for H II regions. Left: $I_{21}$ as a function of $I_{H\alpha}$ 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\alpha,corr}$. Right: $I_{H\alpha} + a_{21}^{C07} I_{21}$ as a function of $I_{H\alpha,corr}$.</figcaption></figure>

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\alpha$ 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\alpha,corr}$ and the hybrid 21 µm and $H\alpha$, and a scatter of 0.37 dex). The deviation from linearity at low $H\alpha$ 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.

### 3.1.3. 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 Eq. (5). For ease of benchmarking we consider the ratio between the inferred $a_{21}$ and the C07-equivalent value of $a_{21}^{C07}$.

Figure 4 shows $a_{21}/a_{21}^{C07}$ 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: that is to say 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\alpha$) ($\rho$ = 0.65, panel a), sSFR estimated from the attenuation-corrected $H\alpha$ ($\rho$ = 0.63, panel b), attenuation-corrected $H\alpha$ surface brightness ($I_{H\alpha,corr}$, $\rho$ = 0.60, panel c), $E(B-V)$ computed from the Balmer decrement ($\rho$ = 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 ($\rho$ = 0.39, panel e). EW($H\alpha$), 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\alpha$ (e.g. Emsellem et al. 2022), and EW($H\alpha$) (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 $|\rho| < 0.3$) are found as a function of galactocentric 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] $\lambda\lambda6717,31/H\alpha$ (panel k) or [O I] $\lambda6300/H\alpha$ (panel l) is flat at low values of these line ratios, but shows a strong negative correlation ($\rho$ = -0.46 and $\rho$ = -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\alpha$ 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Measured coefficient $a_{21}$ normalised to the (rescaled) value from Calzetti 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 16th–84th percentiles are shown as a red shaded area. The quantities considered on the x-axis are (a) equivalent width of $H\alpha$ in emission, (b) sSFR surface density, defined as the ratio between the SFR surface density obtained from dust attenuation corrected $H\alpha$ (assuming a fully sampled IMF) and the stellar mass surface density measured from MUSE data via spectral fitting ($\Sigma_*$), (c) $H\alpha$ luminosity per unit area corrected for dust attenuation with the Balmer decrement ($I_{H\alpha,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]$\lambda\lambda$6717,31/H$\alpha$ ratio, and (l) [O I]$\lambda$6300/H$\alpha$ ratio.</figcaption></figure>

to the 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.

## 3.2. 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<sub>PAH</sub>), 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 Sect. 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<sub>PAH</sub>, F770W, F1000W, and F1130W as ‘PAH-tracing’ bands.

### 3.2.1. PAH-tracing band ratios in H II regions

We adopt an empirical strategy to set fiducial values of the hybridisation coefficients for each band, $a_{band}$. We calculate
<table><caption>**Table 3.** Median (and 16th and 84th percentiles) of the ratio between the $F335M_{\text{PAH}}$, $F770W$, $F1000W$, $F1130W$ to the $F2100W$ luminosity (in units of erg s$^{-1}$ kpc$^{-2}$) for the H II regions in our sample.</caption><thead><tr><th rowspan="2">Band</th><th colspan="3">log (I<sub>band</sub>/I<sub>21</sub>)</th></tr><tr><th>Median</th><th>16th perc.</th><th>84th perc.</th></tr></thead><tbody><tr><td>F335M<sub>PAH</sub></td><td>-0.41</td><td>-0.62</td><td>-0.24</td></tr><tr><td>F770W</td><td>0.65</td><td>0.54</td><td>0.73</td></tr><tr><td>F1000W</td><td>0.10</td><td>-0.04</td><td>0.18</td></tr><tr><td>F1130W</td><td>0.60</td><td>0.45</td><td>0.68</td></tr></tbody></table>

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_{\text{band}}^{\text{C07}}) \equiv \log(a_{21}^{\text{C07}}) - \langle\log(I_{\text{band}}/I_{21})\rangle, \quad (6) $$

where $\langle\log(I_{\text{band}}/I_{21})\rangle$ 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, 16th, and 84th percentiles of the distribution of $I_{\text{band}}/I_{21}$ in Table 3. The scatter in the $I_{\text{band}}/I_{21}$ lies in the range of 0.2 dex (F770W) to 0.4 dex (F335M<sub>PAH</sub>). Part of this larger scatter for F335M<sub>PAH</sub> may be attributed to residuals in the continuum subtraction.

### 3.2.2. H II region scaling relations for PAH-tracing bands

Figure 5 shows the relations between $H\alpha$ corrected using the Balmer decrement and the $H\alpha$ hybridised with each of the PAH-dominated bands using the $a_{\text{band}}^{\text{C07}}$ 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<sub>PAH</sub> for which we obtain a marginally larger scatter of 0.13 dex).

These results are remarkably consistent with those obtained with F2100W in Sect. 3.1.2. In fact, Fig. 5 demonstrates that the hybrid recipe overestimates the dust correction at low surface brightness levels ($\log(I_{H\alpha,corr}/\text{erg s}^{-1}\text{kpc}^{-2}) < 39$), while underestimating it at high surface brightness levels, as already observed for F2100W.

### 3.2.3. Secondary relations of $a_{\text{band}}$ for PAH-tracing bands

In Fig. 6 we show the median trends of the $a_{\text{band}}$ coefficients for different bands as a function of EW($H\alpha$). The values of $a_{\text{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 \text{EW}(H\alpha)/\AA < 1.5$ (or $\log(\text{sSFR}/\text{yr}^{-1}) < -10$) all bands follow a closely matching trend of increasing $a_{\text{band}}$ with EW($H\alpha$), reaching the C07-equivalent value at $\log \text{EW}(H\alpha)/\AA \sim 1.5$.

At high EW($H\alpha$), on the other hand, the behaviour of different bands diverges. The dependence of $a_{21}$ on EW($H\alpha$) flattens for $\log \text{EW}(H\alpha)/\AA > 2$, plateauing to a constant value of $a_{21}$ that is 0.16 dex higher than the C07-equivalent one. The F335M<sub>PAH</sub> 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. (2023a) 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. (2023b) find that the F335M<sub>PAH</sub> band is enhanced and the F1130W band is suppressed in regions of high $H\alpha$ 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_{\text{band}}/I_{21}$ ratios as a function of EW($H\alpha$) 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. Figure 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 function of EW($H\alpha$). $I_{7.7}/I_{21}$ shows a slightly shallower slope than $I_{10}/I_{21}$ and $I_{11.3}/I_{21}$. The $I_{3.3}/I_{21}$ ratio, on the other hand, shows a much flatter (and even increasing) trend.

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; Draine & Li 2007; Khramtsova et al. 2013). Previous work, however, focused on studying entire galaxies or kpc-size regions. On the scale of H II regions in our metallicity range the PAH fraction is found not to depend on metallicity (Egorov et al. 2023), but rather on the intensity of hydrogen-ionising radiation that probably regulates the PAH destruction in H II regions. Overall EW($H\alpha$) remains the most significant secondary correlation for the PAH bands, and we therefore do not consider additional ones in this work.

## 4. 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\alpha$ 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\alpha$)).

### 4.1. Recommendation for using constant $a_{\text{band}}$

We report in Table 4 the median values of the $a_{\text{band}}$ coefficients for our full sample of 19 901 H II regions, together with the C07-equivalent values. Our median $a_{\text{band}}$ are in excellent agreement with the C07-equivalent values, but the scatter around these values (quantified as the 16th and 84th percentiles) are large.

In Table 4 we also present the median values of the $a_{\text{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, F1000W, and F1130W bands (a increases by ~0.4 dex), while it has a smaller effect on F2100W and F335M<sub>PAH</sub> (a increases by ~0.2 dex).
<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of Hα hybridised with F335M<sub>PAH</sub>, 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.</figcaption></figure>

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

## 4.2. 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 populations in fainter regions (Cortese et al. 2008; Leroy et al. 2012; Boquien et al. 2016; Nersesian et al. 2019), and, for the PAH-tracing 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 Sect. 3.1.3, would provide a more accurate estimate of the dust-attenuated star formation, and a prescription more readily applicable to lower-resolution data. Several authors have presented calibrations including such additional terms, for example colours, sSFR, or $M_\star$ (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

$$I_{H\alpha,corr} = I_{H\alpha} + \alpha I_{IR} + \beta I_{Q}, \qquad (7)$$

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⁻¹ kpc⁻².

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 ran a random forest regression, using the algorithm implementation in SCIKIT-LEARN (Pedregosa et al. 2011). In particular, we used 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 used 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\alpha}$, where 'band' is one of the JWST bands available to us (F200W, F300M, F360M, F335M<sub>PAH</sub>, F770W, F1000W, F1130W, F2100W), as input features. We also considered an array of random numbers as an additional input feature to test the performance of the algorithm. The dataset was 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 NIRCam and MIRI and with signal-to-noise ratio grater than 3 in all JWST bands were used for this computation. After training, the mean square error (MSE) of the test set was compared with that of the training sample in order to check the quality of the resulting model and avoid overfitting.

We found that the combination of $I_{21}/I_{H\alpha}$ and $I_{2,0}/I_{H\alpha}$ accounts for nearly all the variance in E(B - V) (Fig. 8). Once
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** log *a*<sub>band</sub> 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<sub>PAH</sub> (blue), F770W (orange), F1000W (green), F1130W (red), F2100W (purple)]. The coloured solid lines represent the median relations, while the coloured dashed lines show the 16th and 84th 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 16th, 50th, 84th percentiles marked and labelled). The purple histogram (right) shows the distribution of log *a*<sub>21</sub>/*a*<sub>21</sub><sup>C07</sup> for the region in our sample, with the 16th, 50th, 84th percentiles marked and labelled.</figcaption></figure>

<table><caption>**Table 4.** Hybridisation coefficients for various JWST bands.</caption><thead><tr><th></th><th></th><th>F335M<sub>PAH</sub></th><th>F770W</th><th>F1000W</th><th>F1130W</th><th>F2100W</th></tr></thead><tbody><tr><td>Constant coeff. C07</td><td>a<sup>C07</sup><sub>band</sub></td><td>0.078</td><td>0.0077</td><td>0.027</td><td>0.0087</td><td>0.034</td></tr><tr><td>Median measured</td><td>a<sub>band</sub></td><td>0.087<sup>+0.07</sup><sub>-0.05</sub></td><td>0.0072<sup>+0.007</sup><sub>-0.004</sub></td><td>0.026<sup>+0.027</sup><sub>-0.015</sub></td><td>0.0084<sup>+0.009</sup><sub>-0.005</sub></td><td>0.031<sup>+0.023</sup><sub>-0.017</sub></td></tr><tr><td>Median measured (backg. subtracted)</td><td>a<sup>bkg sub</sup><sub>band</sub></td><td>0.13<sup>+0.16</sup><sub>-0.08</sub></td><td>0.016<sup>+0.03</sup><sub>-0.01</sub></td><td>0.06<sup>+0.1</sup><sub>-0.04</sub></td><td>0.020<sup>+0.03</sup><sub>-0.014</sub></td><td>0.051<sup>+0.07</sup><sub>-0.03</sub></td></tr></tbody></table>

**Notes.** The top row shows the value of the Calzetti et al. (2007) coefficient, originally derived from Spitzer 24 µm, re-scaled using the median band ratios (see text). The second and third rows show the median *a*<sub>band</sub> obtained in our sample of H II regions before and after applying a local background subtraction.

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 fitted the extinction-corrected Hα surface brightness with a two-parameter regression model, using the Bayesian formalism discussed in Hogg et al. (2010). We fitted for the α and β parameters in Eq. (7) and for intrinsic scatter in the relation, using the Bayesian Monte Carlo sampler EMCEE, and obtained best-fit as
<figure><img src="image_8.png" /><figcaption>Fig. 7. Ratio of $I_{\text{band}}/I_{21}$ normalised to the sample mean ($I_{\text{band}}/I_{21}$) as a function of EW(H$\alpha$), for band = [F335M$_{\text{PAH}}$ (blue), F770W (orange), F1000W (green), F1130W (red)]. The solid lines represent the median, while the dashed lines correspond to the 16th and 84th percentiles. The alternative x-axis on top shows the sSFR associated with each EW(H$\alpha$) value. The median trends show that $I_{3.3}/I_{21}$ does not show the same decrease with EW(H$\alpha$) as the ratios of the other PAH bands.</figcaption></figure>

relation

$$I_{H\alpha,\text{corr}} = I_{H\alpha} + 0.025 I_{21} - 3.0 \times 10^{-4} I_{2}. \quad (8)$$

We also fitted using the F300M and F360M NIRCam bands, obtaining $\alpha = 0.025$, $\beta = 8 \times 10^{-4}$ and $\alpha = 0.025$, $\beta = 7 \times 10^{-4}$, respectively.

Our approach here is similar to the one of Kouroumpatzakis et al. (2023), who used 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) did not consider $I_{H\alpha}$, but modelled 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 repeated our two-parameter fit using stellar mass surface density $\Sigma_{\star}$, which may be used in the absence of JWST NIRCam data. We obtained

$$I_{H\alpha,\text{corr}} = I_{H\alpha} + 0.028 I_{21} - 0.085 \times 10^{31} \Sigma_{\star}, \quad (9)$$

where surface brightnesses are in units of $\text{erg}^{-1}\text{s}^{-1}\text{kpc}^{-2}$ and $\Sigma_{\star}$ is in units of $M_{\odot}\text{kpc}^{-2}$.

## 4.3. IR emission as a tracer 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, 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 was 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) used the early JWST data for four galaxies from our sample and compare the MIR emission with both H$\alpha$ and CO emission, finding that MIR emission correlates well with both extinction-corrected H$\alpha$ 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\alpha,\text{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 evaluated the relative ability of each of the JWST bands considered in this work to trace extinction-corrected H$\alpha$, molecular gas (as traced by CO(2–1)), and stellar mass surface density, focusing our analysis on H II regions specifically. The analysis was performed by building a random forest regression model for the surface brightness in each JWST band and using $I_{H\alpha,\text{corr}}$, $I_{\text{CO}}$, and $\Sigma_{\star}$, 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$\alpha$ 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$_{\text{PAH}}$ mostly traces H$\alpha$, highlighting the success of the continuum subtraction strategy for this band. The MIR bands trace mostly H$\alpha$, 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$\alpha$ and CO are 83% and 16%, respectively. In summary, both 3.3 µm and the MIR bands trace mostly extinction-corrected H$\alpha$ 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.

# 5. Conclusions

We use NIRCam and MIRI images in combination with IFS from MUSE for 19 nearby galaxies to calibrate hybrid recipes to correct H$\alpha$ for dust extinction on the scales of individual H II regions (~100 pc). We examine hybrid recipes that consider a linear combination of H$\alpha$ with the F2100W MIRI band, dominated by emission from small grains, and bands tracing PAH emission, namely the continuum-subtracted NIRCam F335M band (F335M$_{\text{PAH}}$), and the F770W, F1000W, and F1130W MIRI bands. We summarise our main results below.

- The ratio between 21 µm and H$\alpha$ emission is systematically higher in the diffuse ISM, outside optically defined H II
<figure><img src="image_9.png" /><figcaption>**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<sub>21</sub>/I<sub>Hα</sub> and I<sub>2.0</sub>/I<sub>Hα</sub> account for most of the variance in the data.</figcaption></figure>

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\alpha,corr} = I_{H\alpha} + 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_{band}^{C07}$) by rescaling the coefficient derived by Calzetti et al. (2007) for Spitzer 24 µm emission by the average band ratio $\langle I_{24}/I_{band} \rangle$. 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 coefficients 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_{band}^{C07}$ leads to extinction-corrected $H\alpha$ 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, for example EW($H\alpha$), sSFR, or indirectly correlate with them (e.g. $H\alpha$ 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

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Relative importance of extinction-corrected $H\alpha$ ($I_{H\alpha,corr}$), CO(2–1) surface density, and stellar mass surface density ($\Sigma_{\star}$) in predicting the surface brightness of H II regions in a set of JWST bands.</figcaption></figure>

between H II regions and DIG ([O I]/H$\alpha$, [S II]/H$\alpha$). 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\alpha$) for all bands considered. For EW($H\alpha$) < 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\alpha$), $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\alpha$ 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 Eqs. (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\alpha$ in determining the surface brightness for H II regions within our JWST bands. The MIR bands trace mostly $H\alpha$, but with a smaller (15–20%) importance attributed to CO.

The analysis presented in this work extends the validity of previous widely used hybrid SFR calibrations to the ~100 pc scale and demonstrates the transformational ability of JWST to unveil the dust-obscured component of star formation in galaxies.

**Acknowledgements.** This work is carried out as part of the PHANGS Collaboration. Based on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The observations are associated with JWST programme 2107 (PI: J. Lee) based on observations collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), 1100.B-0651 (PHANGS-MUSE; PI: Schinnerer), as well as 094.B-0321
(MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). F.B. acknowledges support from the INAF Fundamental Astrophysics programme 2022. M.B. acknowledges support from FONDECYT regular grant 1211000 and by the ANID BASAL project FB210003. H.A.P. acknowledges support by the National Science and Technology Council of Taiwan under grant 110-2112-M-032-020-MY3. J.M.D.K. gratefully acknowledges funding from the ERC under the European Union's Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant number 714907). COOL Research DAO is a Decentralized Autonomous Organization supporting research in astrophysics aimed at uncovering our cosmic origins. M.C. gratefully acknowledges funding from the DFG through an Emmy Noether Research Group (grant number CH2137/1-1). R.S.K. and S.C.O.G. acknowledge financial support from the ERC via the ERC Synergy Grant 'ECOGAL' (project ID 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) 'STRUCTURES', and from the German Ministry for Economic Affairs and Climate Action in project 'MAINN' (funding ID 50002206). R.S.K. also thanks for computing resources provided bwHPC and DFG through grant INST 35/1134-1 FUGG and for data storage at SDS@hd through grant INST 35/1314-1 FUGG. J.C. acknowledges support from ERC starting grant #851622 DustOrigin. M.Q. acknowledges support from the Spanish grant PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. K.K., O.V.E. and E.J.W. gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI: Kreckel). J.N. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694343). E.R. and H.H. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2022-03499, and the support of the Canadian Space Agency (22ASTALBER). A.K.L. gratefully acknowledges support by grants 1653300 and 2205628 from the National Science Foundation, by award JWST-GO-02107.009-A, and by a Humboldt Research Award from the Alexander von Humboldt Foundation. K.S. acknowledges funding support from grant support by JWST-GO-02107.006-A.

## References

Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621

Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218

Baker, W. M., Maiolino, R., Bluck, A. F. L., et al. 2022, MNRAS, 510, 3622

Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5

Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26

Belfiore, F., Leroy, A., Sun, J., et al. 2023, A&A, 670, A67

Bendo, G. J., Dale, D. A., Draine, B. T., et al. 2006, ApJ, 652, 283

Bendo, G. J., Boselli, A., Dariush, A., et al. 2012, MNRAS, 419, 1833

Bluck, A. F., Maiolino, R., Brownson, S., et al. 2022, A&A, 659, A160

Bolatto, A. D., Simon, J. D., Stanimirović, S., et al. 2007, ApJ, 655, 212

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, A&A, 51, 207

Boquien, M., Kennicutt, R., Calzetti, D., et al. 2016, A&A, 591, A6

Calzetti, D. 2020, Nat. Astron., 4, 437

Calzetti, D., Kennicutt, R. C., Engelbracht, C. W., et al. 2007, ApJ, 666, 870

Catalán-Torrecilla, C., Gil de Paz, A., Castillo-Morales, A., et al. 2015, A&A, 584, A87

Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2019, ApJ, 876, 62

Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023a, ApJ, 944, L11

Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023b, ApJ, 944, L12

Chown, R., Li, C., Parker, L., et al. 2021, MNRAS, 500, 1261

Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157

Crocker, A. F., Calzetti, D., Thilker, D. A., et al. 2013, ApJ, 762, 79

De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, A69

Den Brok, J. S., Chatzigiannakis, D., Bigiel, F., et al. 2021, MNRAS, 504, 3221

Diaz, A. I., Terlevich, E., Vflchez, J. M., et al. 1991, MNRAS, 253, 245

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Draine, B. T., Li, A., Hensley, B. S., et al. 2021, ApJ, 917, 3

Egorov, O. V., Kreckel, K., Sandstrom, K. M., et al. 2023, ApJ, 944, L16

Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892

Groves, B., Kreckel, K., Santoro, F., et al. 2023, MNRAS, 520, 4902

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Haffner, L. M., Dettmar, R.-J. J., Beckman, J. E., et al. 2009, Rev. Mod. Phys., 81, 969

Hao, C. N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124

Hassani, H., Rosolowsky, E., Leroy, A. K., et al. 2023, ApJ, 944, L21

Helou, G., Roussel, H., Appleton, P., et al. 2004, ApJS, 154, 253

Hirashita, H., Buat, V., & Inoue, A. K. 2003, A&A, 410, 83

Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv e-prints [arXiv:1008.4686]

Jarrett, T. H., Masci, F., Tsai, C. W., et al. 2013, AJ, 145, 6

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Kelly, B. C. 2007, ApJ, 665, 1489

Kennicutt, R. C. 1998, ARA&A, 36, 189

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kennicutt, R. C., Calzetti, D., Walter, F., et al. 2007, ApJ, 671, 333

Kennicutt, R. C., Hao, C. N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

Khramtsova, M. S., Wiebe, D. S.,Boley, P. A., & Pavlyuchenkov, Y. N. 2013, MNRAS, 431, 2006

Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2022, MNRAS, 516, 3006

Kouroumpatzakis, K., Zezas, A., Kyritsis, E., Salim, S., & Svoboda, J. 2023, A&A, 673, A16

Lai, T. S.-Y., Smith, J. D. T., Baba, S., Spoon, H. W. W., & Imanishi, M. 2020, ApJ, 905, 55

Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122

Lebouteiller, V., Bernard-Salas, J., Whelan, D. G., et al. 2011, ApJ, 728, 45

Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17

Leroy, A. K., Bigiel, F., de Blok, W. J. G., et al. 2012, AJ, 144, 3

Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24

Leroy, A. K., Hughes, A., Liu, D., et al. 2021a, ApJS, 255, 19

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021b, ApJS, 257, 43

Leroy, A. K., Bolatto, A. D., Sandstrom, K., et al. 2023a, ApJ, 944, L10

Leroy, A. K., Sandstrom, K., Rosolowsky, E., et al. 2023b, ApJ, 944, L9

Leslie, S. K., Schinnerer, E., Groves, B., et al. 2018, A&A, 616, A157

Li, A. 2020, Nat. Astron., 4, 339

Lu, Y., Blanc, G. A., & Benson, A. 2015, ApJ, 808, 129

Lutz, D., Valiante, E., Sturm, E., et al. 2005, ApJ, 625, L83

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Madden, S. C., Galliano, F., Jones, A. P., & Sauvage, M. 2006, A&A, 446, 877

Maragkoudakis, A., Ivkovich, N., Peeters, E., et al. 2018, MNRAS, 481, 5370

Nersesian, A., Xilouris, E. M., Bianchi, S., et al. 2019, A&A, 624, A80

O'Donnell, J. E. 1994, ApJ, 422, 158

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Sausalito: University Science Books)

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825

Peeters, E., Spoon, H. W. W., & Tielens, A. G. G. M. 2004, ApJ, 613, 986

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X

Phillips, M. M., Jenkins, E. B., Dopita, M. A., Sadler, E. M., & Binette, L. 1986, AJ, 91, 1062

Pilyugin, L. S., & Grebel, E. K. 2016, MNRAS, 457, 3678

Regan, M. W., Thornley, M. D., Vogel, S. N., et al. 2006, ApJ, 652, 1112

Riechers, D. A., Pope, A., Daddi, E., et al. 2014, ApJ, 786, 31

Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11

Sandstrom, K. M., Chastenet, J., Sutter, J., et al. 2023a, ApJ, 944, L7

Sandstrom, K. M., Koch, E. W., Leroy, A. K., et al. 2023b, ApJ, 944, L8

Shipley, H. V., Papovich, C., Rieke, G. H., Brown, M. J. I., & Moustakas, J. 2016, ApJ, 818, 60

Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007, ApJ, 656, 770

Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2000, AJ, 120, 3070

Viaene, S., Sarzi, M., Baes, M., Fritz, J., & Puerari, I. 2017, MNRAS, 472, 1286

Westfall, K. B., Cappellari, M., Bershady, M. A., et al. 2019, AJ, 158, 57

Whitcomb, C. M., Sandstrom, K., Murphy, E. J., & Linden, S. 2020, ApJ, 901, 47

Whitcomb, C. M., Sandstrom, K., Leroy, A., & Smith, J. D. T. 2023, ApJ, 948, 88

Williams, T. G., Baes, M., Looze, I. D., et al. 2019, MNRAS, 487, 2753

Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011a, ApJ, 742, 96

Wuyts, S., Förster Schreiber, N. M., Lutz, D., et al. 2011b, ApJ, 738, 106

Zakardjian, A., Pety, J., Herrera, C. N., et al. 2023, ArXiv e-prints [arXiv:2305.03650]

Zhang, L., & Ho, L. C. 2023, ApJ, 943, 60

Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217

<footnote n="1">¹ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Florence, Italy e-mail: francesco.belfiore@inaf.it</footnote>

<footnote n="2">² Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus OH 43210, USA</footnote>

<footnote n="3">³ Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="4">⁴ European Southern Observatory, Karl-Schwarzschild Straße 2,</footnote>
85748 Garching bei München, Germany

5 Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

6 Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile

7 Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessen-bachstrasse 1, 85748 Garching, Germany

8 Sterrenkundig Observatorium, Ghent University, Krijgslaan 281-S9, 9000 Gent, Belgium

9 European Southern Observatory (ESO), Alonso de Córdova 3107, Casilla 19, Santiago 19001, Chile

10 Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA

11 Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, 69120 Heidelberg, Germany

12 Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France

13 Universität Heidelberg, Zentrum für Astronomie, Institut für theoretische Astrophysik, Albert-Ueberle-Straße 2, 69120 Heidelberg, Germany

14 International Centre for Radio Astronomy Research, University of Western Australia, 7 Fairway, Crawley 6009, WA, Australia

15 Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany

16 Max-Planck-Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

17 Observatorio Astronómico Nacional (IGN), C/Alfonso XII, 3, 28014 Madrid, Spain

18 Department of Physics, University of Alberta, Edmonton AB T6G 2E1, Canada

19 Departamento de Física de la Tierra y Astrofísica, Universidad Complutense de Madrid, 28040 Madrid, Spain

20 Instituto de Física de Partículas y del Cosmos IPARCOS, Facultad de CC Físicas, UCM, 28040 Madrid, Spain

21 Department of Astronomy & Astrophysics, University of California, San Diego, 9500 Gilman Drive, San Diego, CA 92093, USA

22 Department of Physics and Astronomy, McMaster University, 1280 Main St. West, Hamilton, ON L8S 4M1, Canada

23 Canadian Institute for Theoretical Astrophysics (CITA), University of Toronto, 60 St George Street, Toronto, ON M5S 3H8, Canada