11institutetext: INAF — Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5,
I-50125, Florence, Italy 11email<EMAIL_ADDRESS>22institutetext:
Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus, OH 43210, USA 33institutetext: Sub-department of Astrophysics,
Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
44institutetext: European Southern Observatory, Karl-Schwarzschild Straße 2,
D-85748 Garching bei München, Germany 55institutetext: Argelander-Institut für
Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany
66institutetext: Centro de Astronomía (CITEVA), Universidad de Antofagasta,
Avenida Angamos 601, Antofagasta, Chile 77institutetext: Max-Planck-Institut
für extraterrestrische Physik (MPE), Giessenbachstrasse 1, D-85748 Garching,
Germany 88institutetext: Sterrenkundig Observatorium, Ghent University,
Krijgslaan 281-S9, 9000 Gent, Belgium 99institutetext: European Southern
Observatory (ESO), Alonso de Córdova 3107, Casilla 19, Santiago 19001, Chile
1010institutetext: Department of Physics and Astronomy, University of Wyoming,
Laramie, WY 82071, USA 1111institutetext: Astronomisches Rechen-Institut,
Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14,
D-69120 Heidelberg, Germany 1212institutetext: Univ Lyon, Univ Lyon1, ENS de
Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-
Genis-Laval France 1313institutetext: Universität Heidelberg, Zentrum für
Astronomie, Institut für theoretische Astrophysik, Albert-Ueberle-Straße 2,
D-69120, Heidelberg, Germany 1414institutetext: International Centre for Radio
Astronomy Research, University of Western Australia, 7 Fairway, Crawley, 6009,
WA, Australia 1515institutetext: Universität Heidelberg, Interdisziplinäres
Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, D-69120
Heidelberg, Germany 1616institutetext: Max-Planck-Institute for Astronomy,
Königstuhl 17, D-69117 Heidelberg, Germany 1717institutetext: Observatorio
Astronómico Nacional (IGN), C/Alfonso XII, 3, E-28014 Madrid, Spain
1818institutetext: Department of Physics, University of Alberta, Edmonton,
Alberta, T6G 2E1, Canada 1919institutetext: Departamento de Física de la
Tierra y Astrofísica, Universidad Complutense de Madrid, E-28040 Madrid, Spain
2020institutetext: Instituto de Física de Partículas y del Cosmos IPARCOS,
Facultad de CC Físicas, UCM, E-28040, Madrid, Spain 2121institutetext:
Department of Astronomy & Astrophysics, University of California, San Diego,
9500 Gilman Drive, San Diego, CA 92093, USA 2222institutetext: Department of
Physics and Astronomy, McMaster University, 1280 Main St. West, Hamilton ON
L8S 4M1 Canada 2323institutetext: Canadian Institute for Theoretical
Astrophysics (CITA), University of Toronto, 60 St George Street, Toronto, ON
M5S 3H8, Canada

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

Francesco Belfiore, Calibrating mid-infrared emission as a tracer of obscured
star formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST<EMAIL_ADDRESS>Adam K. Leroy Calibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWSTCalibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWST Thomas G. Williams
Calibrating mid-infrared emission as a tracer of obscured star formation on
Hii-region scales in the era of JWSTCalibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of JWST
Ashley T. Barnes Calibrating mid-infrared emission as a tracer of obscured
star formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST Frank Bigiel Calibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWST Médéric Boquien Calibrating mid-infrared emission as
a tracer of obscured star formation on Hii-region scales in the era of
JWSTCalibrating mid-infrared emission as a tracer of obscured star formation
on Hii-region scales in the era of JWST Yixian Cao Calibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWSTCalibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWST Jérémy Chastenet Calibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWST Enrico Congiu
Calibrating mid-infrared emission as a tracer of obscured star formation on
Hii-region scales in the era of JWSTCalibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of JWST
Daniel A. Dale Calibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST Oleg V. Egorov Calibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWST Cosima Eibensteiner Calibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWSTCalibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWST Eric Emsellem Calibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWST Simon C. O.
Glover Calibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST Brent Groves Calibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWST Hamid Hassani Calibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of
JWSTCalibrating mid-infrared emission as a tracer of obscured star formation
on Hii-region scales in the era of JWST Ralf S. Klessen Calibrating mid-
infrared emission as a tracer of obscured star formation on Hii-region scales
in the era of JWSTCalibrating mid-infrared emission as a tracer of obscured
star formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWSTCalibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWST Kathryn Kreckel Calibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWST Lukas Neumann
Calibrating mid-infrared emission as a tracer of obscured star formation on
Hii-region scales in the era of JWSTCalibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of JWST
Justus Neumann Calibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST Miguel Querejeta Calibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWST Erik Rosolowsky Calibrating mid-infrared emission as
a tracer of obscured star formation on Hii-region scales in the era of
JWSTCalibrating mid-infrared emission as a tracer of obscured star formation
on Hii-region scales in the era of JWST Patricia Sanchez-Blazquez Calibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWSTCalibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWST Karin
Sandstrom Calibrating mid-infrared emission as a tracer of obscured star
formation on Hii-region scales in the era of JWSTCalibrating mid-infrared
emission as a tracer of obscured star formation on Hii-region scales in the
era of JWST Eva Schinnerer Calibrating mid-infrared emission as a tracer of
obscured star formation on Hii-region scales in the era of JWSTCalibrating
mid-infrared emission as a tracer of obscured star formation on Hii-region
scales in the era of JWST Jiayi Sun Calibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of
JWSTCalibrating mid-infrared emission as a tracer of obscured star formation
on Hii-region scales in the era of JWSTCalibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of
JWSTCalibrating mid-infrared emission as a tracer of obscured star formation
on Hii-region scales in the era of JWST Jessica Sutter Calibrating mid-
infrared emission as a tracer of obscured star formation on Hii-region scales
in the era of JWSTCalibrating mid-infrared emission as a tracer of obscured
star formation on Hii-region scales in the era of JWST Elizabeth J. Watkins
Calibrating mid-infrared emission as a tracer of obscured star formation on
Hii-region scales in the era of JWSTCalibrating mid-infrared emission as a
tracer of obscured star formation on Hii-region scales in the era of JWST

(Received XXX; accepted XXX)

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 the IR 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$ 16000 optically
selected H ii regions across 16 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 an
extinction-corrected H$\alpha$ estimate which agrees 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 H ii 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 H ii regions all bands suffer from an increasing
contamination from the diffuse infrared background. We present calibrations
that take this effect into account by adding an explicit dependence on 3
$\rm\mu m$ emission or stellar mass surface density.

###### Key Words.:

Galaxies: ISM – Galaxies: star formation – ISM: general

##  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 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 $\mu$m, corresponding to the IRAS 25$\rm\mu m$, the Spitzer
MIPS 24 $\mu$m band, or the WISE W4 (22 $\mu$m) filter, traces emission of
small, hot grains transiently heated by the local radiation field. While not
directly relatable to the bolometric dust luminosity without further
assumptions on the shape of the dust SED, emission at 24 $\rm\mu 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 $\rm\mu m$ fluxes, in combination
with either UV or H$\alpha$ 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 $\rm\mu 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$\rm\mu m$, several authors have studied the potential of
using the flux in PAH features or the brightness of PAH-tracing bands (most
commonly the IRAC 8 $\rm\mu m$ band) to probe obscured star formation (Peeters
et al., 2004; Shipley et al., 2016; Maragkoudakis et al., 2018). 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 accessible 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 $\rm\mu 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, 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 $\rm\mu m$
and bands dominated by PAH emission to trace embedded star formation on the
scale of individual H ii regions (${\sim}40{-}110$ pc) in a sample of 16
nearby disc galaxies observed with the MIRI and (partially) with 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$\alpha$ ratios with matched-resolution Balmer decrement
(H$\alpha$/H$\beta$) 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$\alpha$ emission correlates well with MIR emission at 21 $\rm\mu
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$\alpha$ 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.

##  2 Data and Methods

###  2.1 PHANGS–JWST data

We analysed 16 galaxies observed as part of the cycle 1 PHANGS–JWST Treasury
program (ID 02107, PI: J. Lee). The full survey will consist of MIRI and
NIRCam images for the 19 nearby galaxies observed by the PHANGS–MUSE survey
(Emsellem et al., 2022). The data used in this work corresponds to the subset
of the PHANGS–JWST program available at the time of writing. Key properties of
our 16 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 $\rm\mu 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 $\rm\mu 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 no NIRCam data was
obtained, 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. In the case of IC 5332 no MIRI background pointing is
available, so 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), while Leroy et al. (2023b) discuss the determination and anchoring of
the background level for the MIRI imaging. We briefly discuss updates with
respect to the reduction procedure implemented in those works. Firstly, 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.
Further details of the reduction strategy will be described in an upcoming
publication from the PHANGS–JWST team.

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 15 of the 16 targets the MUSE PSF is larger than that
of the lowest-resolution JWST band in our filter set (F2100W). 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 $\rm\mu 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

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

where $\rm F335M_{\rm 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 $\rm F335M_{\rm 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.

Galaxy | D | $\log(\rm M_{\star})$ | $\log(\rm SFR)$ | $i$ | Copt PSF | $\theta_{\rm bm}$ | $\rm N_{HII}$ | Notes
---|---|---|---|---|---|---|---|---
| [Mpc] | $[\rm M_{\odot}]$ | $[\rm M_{\odot}\leavevmode\nobreak\ yr^{-1}]$ | [deg] | [arcsec] | [pc] |  |
IC 5332 | 9.0 | 9.67 | $-0.39$ | 26.9 | 0.87 | 38 | 546 | no NIRCam
NGC 0628 | 9.8 | 10.34 | 0.24 | 8.9 | 0.92 | 44 | 1711 | grand-design spiral
NGC 1087 | 15.9 | 9.93 | 0.12 | 42.9 | 0.92 | 71 | 910 |
NGC 1300 | 19.0 | 10.62 | 0.07 | 31.8 | 0.89 | 82 | 946 | nuclear ring
NGC 1365 | 19.6 | 10.99 | 1.23 | 55.4 | 1.16 | 109 | 555 | AGN, nuclear ring
NGC 1385 | 17.2 | 9.98 | 0.32 | 44.0 | 0.77 | 64 | 912 |
NGC 1433 | 18.6 | 10.87 | 0.05 | 28.6 | 0.91 | 82 | 849 | nuclear ring
NGC 1512 | 18.8 | 10.71 | 0.11 | 42.5 | 1.25 | 114 | 477 | nuclear ring
NGC 1566 | 17.7 | 10.78 | 0.66 | 29.5 | 0.80 | 69 | 1561 | AGN
NGC 1672 | 19.4 | 10.73 | 0.88 | 42.6 | 0.96 | 90 | 934 | AGN, nuclear ring
NGC 3627 | 11.3 | 10.83 | 0.58 | 57.3 | 1.05 | 58 | 973 | AGN
NGC 4303 | 17.0 | 10.52 | 0.73 | 23.5 | 0.78 | 64 | 2042 | AGN, nuclear ring, NIRCam covers half
NGC 4321 | 15.2 | 10.75 | 0.55 | 38.5 | 1.16 | 86 | 1239 | nuclear ring, NIRCam covers half
NGC 4535 | 15.8 | 10.53 | 0.33 | 44.7 | 0.71⋆ | 54 | 1101 | AGN
NGC 5068 | 5.2 | 9.40 |  $-$0.56 | 35.7 | 1.04 | 26 | 1182 |
NGC 7496 | 18.7 | 10.00 | 0.35 | 35.9 | 0.89 | 81 | 536 | AGN
Table 1: 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). $\theta_{\rm bm}$ gives the size of
the common resolution element with no inclination correction. $\rm N_{HII}$
gives the number of HII 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 missing or
partial coverage with NIRCam.

###  2.2 MUSE data and the optical H ii region catalog

We compare the JWST images with maps of H$\alpha$ and H$\beta$ 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${}^{\prime}\times 1^{\prime}$) 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 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$\rm\mu 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)
diagrams [N ii]/H$\alpha$ versus [O iii]/H$\beta$ and [S ii]/H$\alpha$ versus
[O iii]/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 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$\alpha$ or H$\beta$ (0.4%). We obtain a final
sample of 16 474 regions with signal-to-noise ratio larger than three in
F2100W. Their distribution across the sample is summarised in Table 1.

Quantity | Definition | Derived from dataset | Reference
---|---|---|---
$\Sigma_{\rm H_{2}}$ | $\alpha_{\rm CO}I_{\rm CO(2-1)}/R_{21},^{\rm a}$ | ALMA | Leroy et al. (2021b)
E($B-V$) | Eq. 2 | MUSE | Groves et al. (2023)
$\Sigma_{\star}$ | spectral fitting | MUSE | Emsellem et al. (2022)
sSFR |  $C_{\rm H\alpha}\leavevmode\nobreak\ I_{\rm H\alpha,corr}/\Sigma_{\star}$,b | MUSE | Belfiore et al. (2023)
EW(H$\alpha$) | moments of spectrac | MUSE | Groves et al. (2023)
12+log(O/H) | Scal calibrationd | MUSE | Groves et al. (2023)
log(U) | from $\dfrac{[\mathrm{S}\textsc{iii}]\lambda\lambda 9069,9532}{[\mathrm{S}\textsc{ii}]\lambda\lambda 6717,31}$,e | MUSE | Groves et al. (2023)
Table 2: Summary of additional data and physical quantities used in this work.
Notes. a: assuming a Milky Way conversion factor $\alpha_{\rm
CO}=4.35\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\
pc^{-2}(K\leavevmode\nobreak\ km\leavevmode\nobreak\ 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_{\rm H\alpha}/(\mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}/(erg\leavevmode\nobreak\ 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).

###  2.3 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_{\rm
e}=10^{2}\,\mathrm{cm}^{-3}$, leading to $I_{\rm H\alpha,corr}/I_{\rm
H\beta,corr}=2.86$ (Osterbrock & Ferland, 2006). Under these assumptions

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

where $k_{\rm H\alpha}$ and $k_{\rm H\beta}$ are the value of the assumed
reddening curve at the wavelengths of the two emission lines. We assume a
foreground attenuation law following 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_{\rm H\alpha}$, $I_{\rm
H\alpha,corr}$ and $I_{\rm band}$ refer to the luminosity per unit area of
observed H$\alpha$ emission, Balmer decrement corrected H$\alpha$, and one of
the JWST bands, all in units of $\rm erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ 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$\alpha$ assuming full
sampling of the IMF, to allow for a direct comparison with the work of
Belfiore et al. (2023). We define

$\log{\mathrm{(sSFR/yr^{-1})}}=C_{\rm H\alpha}\leavevmode\nobreak\ I_{\rm
H\alpha,corr}/\Sigma_{\star},$ (3)

where $C_{\rm H\alpha}$ 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 Sec. 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 program (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 ($\sim$ 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.

Figure 1: Three-colour images for the 16 galaxies in our sample showing a
combination of H$\alpha$ emission from MUSE IFS (red), 21 $\rm\mu 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.

##  3 Results

Figure 2: Comparison of the H$\alpha$ and 21$\rm\mu m$ emission for NGC 0628.
Top-left, 21$\rm\mu 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$\alpha$ surface brightness from MUSE (Emsellem et al., 2022). White
contours enclose regions of bright 21 $\rm\mu m$ emission
($I_{21}>10^{40.5}\leavevmode\nobreak\ \mathrm{erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ 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 $\rm\mu 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$\alpha$ to F2100W surface brightness. Red contours are bright H$\alpha$
regions, $I_{\rm H\alpha}>10^{38.8}\leavevmode\nobreak\
\mathrm{erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ 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. Figure 3: Left: $I_{21}$ as a function
of $I_{\rm 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_{\rm H\alpha,corr}$. Right: $I_{\rm
H\alpha}+a_{21}^{\rm C07}I_{21}$ as a function of $I_{\rm H\alpha,corr}$.

###  3.1 Using 21 $\rm\mu 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$\alpha$ 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 $\rm\mu m$ as our fiducial
band, we can write the attenuation-corrected H$\alpha$ luminosity per unit
area as

$I_{\mathrm{H\alpha,corr}}^{\rm
hybrid}=I_{\mathrm{H\alpha}}+a_{\mathrm{21}}\leavevmode\nobreak\
I_{\mathrm{21}},$ (4)

where all luminosities are expressed in units of $\rm erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ kpc^{-2}$ (i.e., $I_{\mathrm{21}}=\nu_{21}I_{\nu}$)
and $a_{\mathrm{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_{\rm H\alpha}$ ratio via

$A_{\rm H\alpha}=2.5\log{(1+a_{\mathrm{21}}\leavevmode\nobreak\ I_{21}/I_{\rm
H\alpha})}.$ (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 $\rm\mu m$ MIPS data: $a_{\rm
24}^{\rm C07}=0.031$. To apply this to JWST data we take into account the
different filter bandpasses between MIPS 24 $\rm\mu m$ and MIRI F2100W using
the average flux ratio observed across our sample and tabulated by Leroy et
al. (2023b), $\langle I_{21}/I_{24}\rangle=0.91$, leading to $a_{\rm 21}^{\rm
C07}=0.034$.

####  3.1.1 21 $\rm\mu$m emission in Hii regions and the diffuse component

In this section we present the relative spatial distribution of 21 $\rm\mu m$
and H$\alpha$ emission. In Fig. 1 we show three-colour images combining
H$\alpha$ from MUSE (red) with 21 $\rm\mu 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 $\rm\mu 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 $\rm\mu 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$\alpha$, 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$\alpha$ and 21 $\rm\mu m$ emission more
quantitatively we focus on one of the galaxies in our sample, NGC 628, a
grand-design spiral galaxy with $\rm log(M_{\star}/M_{\odot})=10.3$. Fig. 2
(top row) presents a comparison between the 21 $\rm\mu m$ and H$\alpha$ maps
for NGC 628 at matched (copt) resolution. To quantify the spatial differences
between H$\alpha$ and 21 $\rm\mu m$ emission we compute the 21 $\rm\mu 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 $\rm\mu m$
emission, but only 36% of the H$\alpha$ emission, lies outside the H ii region
masks. Across the full sample the fraction of 21 $\rm\mu m$ flux outside H ii
regions is on average $\sim$ 60%. This fraction is higher than the fraction of
diffuse H$\alpha$, which averages to 38%.

The large amount of 21 $\rm\mu 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 $\rm\mu 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 $\rm\mu 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 $\rm\mu 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 $\rm\mu 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
$\rm\mu 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 $\rm\mu m$ emission.

The difference in diffuse fraction between H$\alpha$ and 21 $\rm\mu m$
emission implies that areas that are faint in H$\alpha$ show a higher
$I_{21}/I_{\rm H\alpha}$ ratio. In the bottom-right panel of Fig. 2 we show a
ratio map of $I_{21}/I_{\rm H\alpha}$ with red contours corresponding to
regions of bright H$\alpha$ emission111$\log{(I_{\rm
H\alpha}/\mathrm{erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\
kpc^{-2}}})>38.8$, chosen to correspond to the 16th percentile of the surface
brightness distribution of H ii regions in the PHANGS-MUSE sample from
Belfiore et al. (2022).. The higher $I_{21}/I_{\rm H\alpha}$ 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 $\rm\mu 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 H ii regions, with a small (few percent)
contribution from old stars (Belfiore et al., 2022). The 21 $\rm\mu 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$\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 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.

Figure 4: The 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 ($\rm\Sigma_{\star}$), c) H$\alpha$ luminosity per
unit area corrected for dust attenuation with the Balmer decrement ($I_{\rm
H\alpha,corr}$), d) $\rm 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.

####  3.1.2 Hii region scaling relations at 21 $\rm\mu 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_{\rm H\alpha}$ and $I_{\rm
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.91) and residual scatter
around the best-fit relation of 0.26 dex.

If one corrects the H$\alpha$ emission for dust attenuation via the Balmer
decrement, we obtain a relation with lower scatter (0.21 dex), but still a
sub-linear slope of 0.85 (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_{\rm H\alpha}$ as a function of $I_{\rm 21}$). For H ii
regions, at $\log(I_{\rm H\alpha}/\mathrm{erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ 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_{\rm H\alpha,corr}$ and $I_{\rm 21}$ are good proxies for
the total ionising photon production rate.

Finally, we show the relation between $I_{\rm H\alpha,corr}$ and the hybrid 21
$\rm\mu m$ and H$\alpha$ tracer (Eq. 4) using the constant $a_{\rm 21}^{C07}$
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$\alpha$, 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$\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_{\rm H\alpha,corr}$ and the hybrid 21 $\rm\mu 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
$\rm\mu 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_{21}^{\rm C07}$.

Fig. 4 shows $a_{21}/a_{21}^{\rm 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: 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 ($\sim 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.64$, panel b), attenuation-
corrected H$\alpha$ surface brightness ($I_{\rm H\alpha,corr}$, $\rho=0.60$,
panel c), E(B$-$V) computed from the Balmer decrement ($\rho=0.55$, panel d),
and the ratio of 21 $\rm\mu m$ to 2.0 $\rm\mu m$ surface brightness, where the
surface brightness at 2.0 $\rm\mu m$ is taken from the NIRCam imaging data in
the F200W filter ($\rho=0.37$, 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 $I_{21}$ (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\lambda$6717,31/H$\alpha$ (panel k) or [O i]$\lambda$6300/H$\alpha$
(panel l) is flat at low values of these line ratios, but shows a strong
negative correlation ($\rho=-0.48$ and $\rho=-0.50$, 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 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 $\rm\mu 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 (F335MPAH), and the MIRI bands centred at 7.7 $\rm\mu m$
(F770W), 10 $\rm\mu m$ (F1000W), and 11.3 $\rm\mu m$ (F1130W). Emission in the
7.7 $\rm\mu m$ and 11.3 $\rm\mu 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 $\rm\mu m$ band emission remains unclear
(Smith et al., 2007; Li, 2020). Leroy et al. (2023b) find that the 10 $\rm\mu
m$ continuum closely follows the PAH-dominated bands, rather than the 21
$\rm\mu m$ emission. In this work we therefore refer collectively to F335MPAH,
F770W, F1000W, and F1130W as ‘PAH-tracing’ bands.

$\log{(I_{\rm band}/I_{21})}$
---
Band | Median | 16th perc. | 84th perc.
F335MPAH |  $-$0.38 |  $-$0.62 |  $-$0.19
F770W | 0.65 | 0.53 | 0.73
F1000W | 0.09 |  $-$0.04 | 0.17
F1130W | 0.59 | 0.44 | 0.68
Table 3: Median (and 16th and 84th percentiles) of the ratio between the
F335MPAH, F770W, F1000W, F1130W to the F2100W luminosity (in units of $\rm
erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ kpc^{-2}$) for the H ii
regions in our sample. The dependence of the $I_{\rm band}/I_{21}$ on
EW(H$\alpha$) is shown in Fig. 7. Figure 5: Comparison of H$\alpha$ hybridised
with F335MPAH, F770W, F1000W and F1130W and the attenuation-corrected
H$\alpha$ 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.

####  3.2.1 PAH-tracing band ratios in Hii regions

We adopt an empirical strategy to set fiducial values of the hybridisation
coefficients for each band, $a_{\rm band}$. We calculate these coefficients,
which we refer to as ‘C07-equivalent’, by scaling the value of the Calzetti et
al. (2007) 21 $\rm\mu m$ coefficient using the average band ratio between PAH-
tracing bands and F2100W. In detail, we define

$\log{(a_{\rm band}^{\rm C07})}\equiv\log{(a_{21}^{\rm
C07})}-\langle\log{(I_{\rm band}/I_{21})}\rangle,$ (6)

where $\langle\log{(I_{\rm band}/I_{21})}\rangle$ is the median band ratio
between the bands of interest (3.3, 7.7, 10, and 11.3 $\rm\mu m$) and 21
$\rm\mu m$ in our H ii region sample. We report the median, 16th, and 84th
percentiles of the distribution of $I_{\rm band}/I_{21}$ in Table 3. The
scatter in the $I_{\rm band}/I_{21}$ lies in the range of 0.2 dex (F770W) to
0.4 dex (F335MPAH). Part of this larger scatter for F335MPAH may be attributed
to residuals in the continuum subtraction.

####  3.2.2 Hii region scaling relations for PAH-tracing bands

Fig. 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_{\rm band}^{\rm 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.87, but all consistent with each other within the errors) and
$\sim$0.10 dex scatter (except for F335MPAH for which we obtain a marginally
larger scatter of 0.14 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 ($\rm
log(I_{H\alpha,corr}/erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\
kpc^{-2})<39$), while underestimating it at high surface brightness levels, as
already observed for F2100W.

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

Figure 6: $\log{a_{\rm band}}$ normalised to the C07-equivalent values as a function of EW(H$\alpha$) and sSFR (alternative x-axis on top) for H ii regions in our sample, where band = [F335MPAH (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$\alpha$ \+ WISE W4 data from Belfiore et al. (2023). The grey histogram (top) shows the distribution of EW(H$\alpha$) 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_{\rm 21}/a_{\rm 21}^{\rm C07}}$ for the region in our sample, with the 16th, 50th, 84th percentiles marked and labelled.  |  | F335MPAH | F770W | F1000W | F1130W | F2100W
---|---|---|---|---|---|---
constant coeff. C07 | $a_{\rm band}^{\rm C07}$ | 0.078 | 0.0077 | 0.027 | 0.0087 | 0.034
median measured | $a_{\rm band}$ | 0.079${}^{+0.08}_{-0.05}$ | 0.0072${}^{+0.007}_{-0.004}$ | 0.026${}^{+0.028}_{-0.015}$ | 0.0082${}^{+0.009}_{-0.005}$ | 0.031${}^{+0.023}_{-0.017}$
median measured (backg. subtracted) | $a_{\rm band}^{\rm bkg\leavevmode\nobreak\ sub}$ | 0.13${}^{+0.20}_{-0.08}$ | 0.017${}^{+0.03}_{-0.01}$ | 0.06${}^{+0.1}_{-0.04}$ | 0.020${}^{+0.04}_{-0.014}$ | 0.051${}^{+0.07}_{-0.03}$
Table 4: Hybridisation coefficients for various JWST bands. The top row shows
the value of the Calzetti et al. (2007) coefficient, originally derived from
Spitzer 24$\rm\mu m$, re-scaled using the median band ratios (see text). The
second and third rows show the median $a_{\rm band}$ obtained in our sample of
H ii regions before and after applying a local background subtraction.

In Fig. 6 we show the median trends of the $a_{\rm band}$ coefficients for
different bands as a function of EW(H$\alpha$). The values of $a_{\rm 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{\rm EW(\textup{H}\alpha)/\AA}<1.5$ (or
$\log(\rm sSFR/yr^{-1})<-10$) all bands follow a closely matching trend of
increasing $a_{\rm band}$ with EW(H$\alpha$), reaching the C07-equivalent
value at $\log{\rm EW(\textup{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{\rm
EW(\textup{H}\alpha)/\AA}>2$, plateauing to a constant value of $a_{21}$ that
is 0.16 dex higher than the C07-equivalent one. The F335MPAH 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 F335MPAH 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.

Figure 7: The ratio of $\rm I_{band}/I_{21}$ normalised to the sample mean
($\rm\langle I_{band}/I_{21}\rangle$) as a function of EW(H$\alpha$), for
band= [F335MPAH (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.

We show the change in the $I_{\rm 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. 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 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, with only a
small decrease in the median ratio for $\rm\log{(EW(H\alpha/\AA))}>1.8$.

##  4 Discussion

This work presents and unprecedented set of resolved IR data for star-forming
regions at $\sim$ 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
$\rm\mu 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 $\sim$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_{\rm band}$

We report in Table 4 the median values of the $a_{\rm band}$ coefficients for
our full sample of 16 474 H ii regions, together with the C07-equivalent
values. Our median $a_{\rm 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_{\rm band}$
coefficients after performing a local background subtraction. These values are
0.22 to 0.38 dex larger than our nominal ones. The background subtraction more
strongly affects the F770W, F1000W, and F1130W bands ($a$ increases by $\sim$
0.4 dex), while it has a smaller effect on F2100W and F335MPAH ($a$ increases
by $\sim$ 0.2 dex).

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_{\rm band}$ will also depend on the distribution of EW(H$\alpha$) 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$\rm\mu 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 $\rm\mu 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 $\sim$ 15 higher spatial resolution than WISE W4, reveal
that the relation between $a_{21}$ and EW(H$\alpha$) plateaus at higher
EW(H$\alpha$) (and $a_{21}$) than seen in the low-resolution data. Since our
observations resolve the inter-cloud distance ($\sim$ 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_{\rm band}$ with EW(H$\alpha$). 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$\alpha$).

Figure 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_{\rm H\alpha}$ and
$I_{2.0}/I_{\rm H\alpha}$ account for most of the variance in the data.

###  4.2 Modelling the impact of old stellar population on $\rm a_{band}$

The trend of increasing $a_{\rm band}$ with EW(H$\alpha$) 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_{\rm band}$ as a function of EW(H$\alpha$), 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 lower-resolution data. Several authors have
presented calibrations including such additional terms, for example colours,
sSFR, or $\rm 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_{\rm H\alpha,corr}=I_{\rm H\alpha}+\alpha\leavevmode\nobreak\ I_{\rm
IR}+\beta\leavevmode\nobreak\ I_{Q},$ (7)

where $I_{\rm 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 $\rm erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ 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 Gini impurity 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_{\rm band}/I_{\rm H\alpha}$, where ‘band’ is one of the
JWST bands avaiable to us (F200W, F300M, F360M, F335MPAH, 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 NIRCam 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_{\rm H\alpha}$ and $I_{3.0}/I_{\rm
H\alpha}$ 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 mass-tracing NIRCam bands, of which the best is F300M.

Having established $I_{21}$ and $I_{3.0}$ as the two best variables, we fit
the extinction-corrected H$\alpha$ surface brightness with a two-parameter
regression model, using the Bayesian formalism discussed in Hogg et al.
(2010). We fit for the $\alpha$ and $\beta$ parameters in Eq. 7 and for
intrinsic scatter in the relation, using the Bayesian Monte Carlo sampler
emcee, and obtain best-fit as relation

$I_{\rm H\alpha,corr}=I_{\rm H\alpha}+0.026\leavevmode\nobreak\
I_{21}-0.0011\leavevmode\nobreak\ I_{3},$ (8)

We also fit using the F200W and F360M NIRCam bands, obtaining respectively
$\alpha=0.026$, $\beta=4.2\leavevmode\nobreak\ 10^{-4}$, and $\alpha=0.025$,
$\beta=-0.001$.

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_{\rm H\alpha}$, 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
$\Sigma_{\star}$, which may be used in the absence of JWST NIRCam data. We
obtain

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

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

###  4.3 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., 2022; 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 $\rm\mu m$ and CO(2-1) emission. This correlation is found to
have lower scatter than that between CO and IR emission at $\sim$ 24 $\rm\mu
m$. Whitcomb et al. (2022) similarly argued using correlation analysis that
emission around 12 $\rm\mu m$ traces CO emission better than SFR, while IR
emission at 24 $\rm\mu 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$\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_{\rm
H\alpha,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$\alpha$,
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_{\rm H\alpha,corr}$, $I_{\rm 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 $\sim$ 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). F335MPAH
mostly traces H$\alpha$, highlighting the success of the continuum subtraction
strategy for this band. The MIR bands trace mostly H$\alpha$, but with 15-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
7$\%$ 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 81% and 16% respectively. In
summary, both 3.3 $\rm\mu 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.

Figure 9: Relative importance of extinction-corrected H$\alpha$ ($\rm
I_{H\alpha,corr}$), CO(2-1) surface density, and stellar mass surface density
($\rm\Sigma_{\star}$) in predicting the surface brightness of H ii regions in
a set of JWST bands.

##  5 Conclusions

We use NIRCam and MIRI images in combination with IFS from MUSE for 16 nearby
galaxies to calibrate hybrid recipes to correct H$\alpha$ for dust extinction
on the scales of individual H ii regions ($\sim$ 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 (F335MPAH), and the F770W,
F1000W, and F1130W MIRI bands. We summarise our main results below.

  * •

The ratio between 21 $\rm\mu m$ and H$\alpha$ 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 $\rm\mu m$ emission.

  * •

Focusing on H ii regions, we calibrate the coefficient in front of the IR
emission term ($I_{\rm H\alpha,corr}=I_{\rm H\alpha}+a_{\rm
band}\leavevmode\nobreak\ I_{\rm 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_{\rm band}^{\rm C07}$) by rescaling
the coefficient derived by Calzetti et al. (2007) for Spitzer 24 $\rm\mu m$
emission by the average band ratio $\langle I_{24}/I_{\rm band}\rangle$. These
C07-equivalent coefficients are in good general agreement with the median of
the measured $a_{\rm band}$ coefficient across our sample of H ii regions for
all bands considered (F335MPAH, F770W, F1000W, F1130W, F2100W).

  * •

For both 21 $\rm\mu m$ and the PAH-tracing bands, using a hybrid recipe with a
constant coefficient $a_{\rm band}^{\rm C07}$ leads to extinction-corrected
H$\alpha$ fluxes that correlate well (scatter $\sim$ 0.1 dex), albeit sub-
linearly (slope 0.82-0.87), with those inferred from the Balmer decrement. The
sub-linear slope implies that a constant $a_{\rm band}$ overestimates the dust
correction for faint H ii regions, and it underestimates it for bright
regions.

  * •

We tested for secondary dependencies of $a_{\rm 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$\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 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_{\rm 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 F335MPAH shows an intermediate
behaviour between F2100W and the other MIR bands, which can be explained by
the relative increase in 3.3 $\rm\mu 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.

  * •

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 $\rm\mu m$ or Spitzer 24
$\rm\mu m$ bands.

  * •

The main secondary dependence of $a_{\rm band}$ can be modelled by explicitly
adding an additional linear term proportional to stellar mass surface density
or 3 $\rm\mu 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$\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.

###### 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 program 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). FB acknowledges support from the INAF Fundamental Astrophysics
program 2022. MB acknowledges support from FONDECYT regular grant 1211000 and
by the ANID BASAL project FB210003. HAP acknowledges support by the National
Science and Technology Council of Taiwan under grant 110-2112-M-032-020-MY3.
JMDK 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. MC gratefully acknowledges funding from the DFG through an
Emmy Noether Research Group (grant number CH2137/1-1). RSK and SCOG
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 50OO2206). RSK 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. JC acknowledges support from ERC starting grant #851622
DustOrigin. MQ acknowledges support from the Spanish grant
PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. KK, OVE and
EJW 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). JN acknowledges funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No. 694343). ER and HH
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). AKL 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. KS acknowledges funding support from grant support by
JWST-GO-02107.006-A

## References

  * Anand et al. (2021) Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621
  * Aniano et al. (2011) Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218
  * Baker et al. (2022) Baker, W. M., Maiolino, R., Bluck, A. F. L., et al. 2022, 7, 1
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Belfiore et al. (2023) Belfiore, F., Leroy, A., Sun, J., et al. 2023, A&A, 670, A67
  * Belfiore et al. (2022) Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26
  * Bendo et al. (2012) Bendo, G. J., Boselli, A., Dariush, A., et al. 2012, MNRAS, 419, 1833
  * Bendo et al. (2006) Bendo, G. J., Dale, D. A., Draine, B. T., et al. 2006, ApJ, 652, 283
  * Bluck et al. (2022) Bluck, A. F., Maiolino, R., Brownson, S., et al. 2022, A&A, 659, A160
  * Bolatto et al. (2007) Bolatto, A. D., Simon, J. D., Stanimirović, S., et al. 2007, ApJ, 655, 212
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Boquien et al. (2016) Boquien, M., Kennicutt, R., Calzetti, D., et al. 2016, A&A, 591, A6
  * Calzetti (2020) Calzetti, D. 2020, Nat. Astron., 4, 437
  * Calzetti et al. (2007) Calzetti, D., Kennicutt, R. C., Engelbracht, C. W., et al. 2007, ApJ, 666, 870
  * Catalán-Torrecilla et al. (2015) Catalán-Torrecilla, C., Gil de Paz, A., Castillo-Morales, A., et al. 2015, A&A, 584, A87
  * Chastenet et al. (2019) Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2019, ApJ, 876, 62
  * Chastenet et al. (2023a) Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023a, ApJL, 944, L12
  * Chastenet et al. (2023b) Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023b, ApJL, 944, L11
  * Chown et al. (2021) Chown, R., Li, C., Parker, L., et al. 2021, MNRAS, 500, 1261
  * Cortese et al. (2008) Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157
  * Crocker et al. (2013) Crocker, A. F., Calzetti, D., Thilker, D. A., et al. 2013, ApJ, 762, 79
  * De Looze et al. (2014) De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, 1
  * Den Brok et al. (2021) Den Brok, J. S., Chatzigiannakis, D., Bigiel, F., et al. 2021, MNRAS, 504, 3221
  * Diaz et al. (1991) Diaz, A. I., Terlevich, E., Vflchez, J. M., et al. 1991, MNRAS, 253, 245
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Draine et al. (2021) Draine, B. T., Li, A., Hensley, B. S., et al. 2021, ApJ, 917, 3
  * Egorov et al. (2023) Egorov, O. V., Kreckel, K., Sandstrom, K. M., et al. 2023, ApJ, 944, L16
  * Emsellem et al. (2022) Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191
  * Galliano et al. (2018) Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673
  * Groves et al. (2012) Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892
  * Groves et al. (2023) Groves, B., Kreckel, K., Santoro, F., et al. 2023, MNRAS, 520, 4902
  * Gruppioni et al. (2013) Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23
  * Haffner et al. (2009) Haffner, L. M., Dettmar, R.-J. J., Beckman, J. E., et al. 2009, Rev. Mod. Phys., 81, 969
  * Hao et al. (2011) Hao, C. N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124
  * Hassani et al. (2023) Hassani, H., Rosolowsky, E., Leroy, A. K., et al. 2023, ApJL, 944, L21
  * Helou et al. (2004) Helou, G., Roussel, H., Appleton, P., et al. 2004, ApJS, 154, 253
  * Hirashita et al. (2003) Hirashita, H., Buat, V., & Inoue, A. K. 2003, A&A, 410, 83
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010 [arXiv:1008.4686]
  * Jarrett et al. (2013) Jarrett, T. H., Masci, F., Tsai, C. W., et al. 2013, AJ, 145, 6
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
  * Kennicutt (1998) Kennicutt, R. C. 1998, ARA&A, 36, 189
  * Kennicutt et al. (2007) Kennicutt, R. C., Calzetti, D., Walter, F., et al. 2007, ApJ, 671, 333
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531
  * Kennicutt et al. (2009) Kennicutt, R. C., Hao, C. N., Calzetti, D., et al. 2009, ApJ, 703, 1672
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kim et al. (2022) Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2022, MNRAS, 516, 3006
  * Kouroumpatzakis et al. (2023) Kouroumpatzakis, K., Zezas, A., Kyritsis, E., Salim, S., & Svoboda, J. 2023, arXiv, 1
  * Lai et al. (2020) Lai, T. S.-Y., Smith, J. D. T., Baba, S., Spoon, H. W. W., & Imanishi, M. 2020, ApJ, 905, 55
  * Lang et al. (2020) Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122
  * Lebouteiller et al. (2011) Lebouteiller, V., Bernard-Salas, J., Whelan, D. G., et al. 2011, ApJ, 728
  * Lee et al. (2023) Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJL, 944, L17
  * Leroy et al. (2012) Leroy, A. K., Bigiel, F., de Blok, W. J. G., et al. 2012, AJ, 144, 3
  * Leroy et al. (2023a) Leroy, A. K., Bolatto, A. D., Sandstrom, K., et al. 2023a, ApJ, 944, L10
  * Leroy et al. (2021a) Leroy, A. K., Hughes, A., Liu, D., et al. 2021a, ApJS, 255, 19
  * Leroy et al. (2023b) Leroy, A. K., Sandstrom, K., Rosolowsky, E., et al. 2023b, ApJL, 944, L9
  * Leroy et al. (2019) Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24
  * Leroy et al. (2021b) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021b, ApJS, 257, 43
  * Leslie et al. (2018) Leslie, S. K., Schinnerer, E., Groves, B., et al. 2018, A&A, 616, 1
  * Li (2020) Li, A. 2020, Nat. Astron., 4, 339
  * Lu et al. (2015) Lu, Y., Blanc, G. A., & Benson, A. 2015, ApJ, 808, 129
  * Lutz et al. (2005) Lutz, D., Valiante, E., Sturm, E., et al. 2005, ApJ, 625, L83
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415
  * Maragkoudakis et al. (2018) Maragkoudakis, A., Ivkovich, N., Peeters, E., et al. 2018, MNRAS, 481, 5370
  * Nersesian et al. (2019) Nersesian, A., Xilouris, E. M., Bianchi, S., et al. 2019, A&A, 624, A80
  * O’Donnell (1994) O’Donnell, J. E. 1994, ApJ, 422, 158
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825
  * Peeters et al. (2004) Peeters, E., Spoon, H. W. W., & Tielens, A. G. G. M. 2004, ApJ, 613, 986
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Sp. Telesc. Instrum. 2014 Opt. Infrared, Millim. Wave, 9143, 91433X
  * Phillips et al. (1986) Phillips, M. M., Jenkins, E. B., Dopita, M. A., Sadler, E. M., & Binette, L. 1986, AJ, 91, 1062
  * Pilyugin & Grebel (2016) Pilyugin, L. S. & Grebel, E. K. 2016, MNRAS, 457, 3678
  * Regan et al. (2006) Regan, M. W., Thornley, M. D., Vogel, S. N., et al. 2006, ApJ, 652, 1112
  * Riechers et al. (2014) Riechers, D. A., Pope, A., Daddi, E., et al. 2014, ApJ, 786, 31
  * Salim et al. (2018) Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11
  * Sandstrom et al. (2023a) Sandstrom, K. M., Chastenet, J., Sutter, J., et al. 2023a, ApJ, 944, L7
  * Sandstrom et al. (2023b) Sandstrom, K. M., Koch, E. W., Leroy, A. K., et al. 2023b, ApJL, 944, L8
  * Shipley et al. (2016) Shipley, H. V., Papovich, C., Rieke, G. H., Brown, M. J. I., & Moustakas, J. 2016, ApJ, 818, 60
  * Smith et al. (2007) Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007, ApJ, 656, 770
  * Thilker et al. (2000) Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2000, AJ, 120, 3070
  * Viaene et al. (2017) Viaene, S., Sarzi, M., Baes, M., Fritz, J., & Puerari, I. 2017, MNRAS, 472, 1286
  * Westfall et al. (2019) Westfall, K. B., Cappellari, M., Bershady, M. A., et al. 2019, AJ, 158, 57
  * Whitcomb et al. (2022) Whitcomb, C. M., Sandstrom, K., Leroy, A., & Smith, J. D. T. 2022, ArXiv:2212.00180v1 [arXiv:2212.00180]
  * Williams et al. (2019) Williams, T. G., Baes, M., Looze, I. D., et al. 2019, MNRAS, 487, 2753
  * Wuyts et al. (2011a) Wuyts, S., Forster Schreiber, N. M., Lutz, D., et al. 2011a, ApJ, 738, 106
  * Wuyts et al. (2011b) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011b, ApJ, 742, 96
  * Zakardjian et al. (2023) Zakardjian, A., Pety, J., Herrera, C. N., et al. 2023, arXiv, arXiv:2305.03650
  * Zhang et al. (2017) Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217
  * Zhang & Ho (2023) Zhang, L. & Ho, L. C. 2023, ApJ, 943, 60