11institutetext:  Max-Planck-Institute for Astronomy, Königstuhl 17, D-69117
Heidelberg, Germany 22institutetext: Astronomisches Rechen-Institut, Zentrum
für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, D-69120
Heidelberg, Germany 33institutetext: INAF — Osservatorio Astrofisico di
Arcetri, Largo E. Fermi 5, I-50125, Florence, Italy 44institutetext: Research
School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia 55institutetext: Departamento de Astronomía,
Universidad de Chile, Santiago,Chile 66institutetext: Department of Physics
and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
77institutetext: Observatories of the Carnegie Institution for Science,
Pasadena, CA, USA 88institutetext: Sterrenkundig Observatorium, Universiteit
Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium 99institutetext: Universität
Heidelberg, Zentrum für Astronomie, Institut für theoretische Astrophysik,
Albert-Ueberle-Straße 2, D-69120, Heidelberg, Germany 1010institutetext:
Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany
1111institutetext: Max-Planck-Institute for extraterrestrial Physics,
Giessenbachstraße 1, D-85748 Garching, Germany 1212institutetext: Observatorio
Astronómico Nacional (IGN), C/Alfonso XII, 3, E-28014 Madrid, Spain
1313institutetext: European Southern Observatory, Karl-Schwarzschild-Straße 2,
85748 Garching, Germany 1414institutetext: Univ Lyon, Univ Lyon1, ENS de Lyon,
CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-
Laval France 1515institutetext: Sydney Institute for Astronomy, School of
Physics, Physics Road, The University of Sydney, Darlington 2006, NSW,
Australia 1616institutetext: Argelander-Institut für Astronomie, Universität
Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany 1717institutetext: Centro de
Astronomía (CITEVA), Universidad de Antofagasta, Avenida Angamos 601,
Antofagasta, Chile 1818institutetext: Department of Physics and Astronomy,
University of Wyoming, Laramie, WY 82071, USA 1919institutetext: Gemini
Observatory/NSF’s NOIRLab, 950 N. Cherry Avenue, Tucson, AZ, 85719, USA
2020institutetext: Department of Astronomy, The Ohio State University, 140
West 18th Avenue, Columbus, OH 43210, USA 2121institutetext: Department of
Physics, Tamkang University, No.151, Yingzhuan Rd., Tamsui Dist., New Taipei
City 251301, Taiwan 2222institutetext: Department of Physics, University of
Alberta, Edmonton, AB T6G 2E1, Canada 2323institutetext: Department of
Physics, General Studies, College of Engineering, Nihon University, 1
Nakagawara, Tokusada, Tamuramachi, Koriyama, Fukushima, 963-8642, Japan
2424institutetext: National Astronomical Observatory of Japan, 2-21-1 Osawa,
Mitaka, Tokyo, 181-8588, Japan 2525institutetext: Departamento de Física de la
Tierra y Astrofísica, Universidad Complutense de Madrid, E-28040 Madrid, Spain

# PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies

Francesco Santoro<EMAIL_ADDRESS>the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxies Kathryn Kreckel PHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxies Francesco Belfiore
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Brent Groves PHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxies Enrico Congiu PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies David A. Thilker
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Guillermo A. Blanc PHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxies Eva Schinnerer PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies I-Ting Ho PHANGS–MUSE: the
H ii region luminosity function of local star-forming galaxiesPHANGS–MUSE: the
H ii region luminosity function of local star-forming galaxies J. M. Diederik
Kruijssen PHANGS–MUSE: the H ii region luminosity function of local star-
forming galaxiesPHANGS–MUSE: the H ii region luminosity function of local
star-forming galaxies Sharon Meidt PHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxies Ralf S. Klessen PHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxies Andreas Schruba
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Miguel Querejeta PHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxies Ismael Pessa PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies Mélanie Chevance
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Jaeyeon Kim PHANGS–MUSE: the H ii region luminosity function of local
star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of local
star-forming galaxies Eric Emsellem PHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxies Rebecca McElroy PHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxies Ashley T. Barnes
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Frank Bigiel PHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxies Médéric Boquien PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies Daniel A. Dale PHANGS–MUSE:
the H ii region luminosity function of local star-forming galaxiesPHANGS–MUSE:
the H ii region luminosity function of local star-forming galaxies Simon C. O.
Glover PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Kathryn Grasha PHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxies Janice Lee PHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity
function of local star-forming galaxies Adam K. Leroy PHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii
region luminosity function of local star-forming galaxies Hsi-An Pan
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Erik Rosolowsky PHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function of
local star-forming galaxies Toshiki Saito PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies Patricia Sanchez-Blazquez
PHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxiesPHANGS–MUSE: the H ii region luminosity function of local star-forming
galaxies Elizabeth J. Watkins PHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxiesPHANGS–MUSE: the H ii region luminosity function
of local star-forming galaxies Thomas G. Williams PHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxiesPHANGS–MUSE: the H ii region
luminosity function of local star-forming galaxies

(Received August 30, 2021; accepted November 11, 2021)

We use an unprecedented sample of about $23\,000$ H ii regions detected at an
average physical resolution of $67$ pc in the PHANGS–MUSE sample to study the
extragalactic H ii region H$\alpha$ luminosity function (LF). Our observations
probe the star-forming disk of 19 nearby spiral galaxies with low inclination
and located close to the star formation main sequence at $z=0$. The mean LF
slope $\alpha$ in our sample is $=1.73$ with a $\sigma$ of $0.15$. We find
that $\alpha$ decreases with the galaxy’s star formation rate surface density
$\Sigma_{\mathrm{SFR}}$and argue that this is driven by an enhanced clustering
of young stars at high gas surface densities. Looking at the H ii regions
within single galaxies we find that no significant variations occur between
the LF of the inner and outer part of the star-forming disk, whereas the LF in
the spiral arm areas is shallower than in the inter-arm areas for six out of
the 13 galaxies with clearly visible spiral arms. We attribute these
variations to the spiral arms increasing the molecular clouds arm–inter-arm
mass contrast and find suggestive evidence that they are more evident for
galaxies with stronger spiral arms. Furthermore, we find systematic variations
in $\alpha$ between samples of H ii regions with high and low ionization
parameter $q$ and argue that they are driven by the aging of H ii regions.

###### Key Words.:

ISM: H ii regions – ISM: structure – Galaxies: ISM – Galaxies: star formation
– Galaxies: evolution – Galaxies: spirals

##  1 Introduction

H$\alpha$ emission is one of the most effective tracers of young stars that
formed within the last $10$ Myr (e.g. Kennicutt & Evans, 2012a; Haydon et al.,
2020). In star-forming regions, it originates from massive OB-type stars which
energetic (hard UV) radiation is able to heat and ionize the surrounding gas
to form an H ii region. This is why H ii regions have long been considered as
the optimal probes of massive star formation in galaxies (e.g. Kennicutt et
al., 1989; Thilker et al., 2000; Lawton et al., 2010), even though we should
not forget that star formation currently ongoing within dust enshrouded dense
cores may remain hidden for about $3$ Myr (Kim et al., 2021). Since the
H$\alpha$ luminosity of an H ii region is directly related to the amount of
ionizing radiation emitted by the OB stars at its center, the H ii region
luminosity function (LF) allows us to constrain the mass function (MF) of
young stellar regions.

Similarities between the spatial distribution of stars, stellar associations
and clusters on different physical scales indicate that star formation is a
scale-free process (e.g. Efremov & Elmegreen, 1998; Bastian et al., 2009;
Kruijssen, 2012; Hopkins, 2013; Krumholz, 2014, and references therein). For
this reason, the mass/luminosity function of star-forming regions is expected
to follow a power-law with slope of about $-2$ (Elmegreen & Falgarone, 1996).
Studies probing the H ii region LFs using different tracers largely agree with
this expectation, measuring LF slopes close to $-2$ with minor variations of
$\pm 0.2{-}0.5$ among galaxies (e.g. UV: Cook et al., 2016, H$\alpha$:
Kennicutt et al., 1989, Pa$\alpha$: Liu et al., 2013, infrared and radio:
Mascoop et al., 2021). It is worth noticing that the main physical driver
behind this is hierarchical growth under the influence of gravity.

The accessibility of the H$\alpha$ emission line to ground-based observations
favored the study of the so-called nebular LF (i.e. built using the H$\alpha$
luminosity of H ii regions) in numerous extra-galactic H$\alpha$ surveys (e.g.
Knapen, 1998; Thilker et al., 2002; Kennicutt et al., 2008). Thanks to studies
in the Milky Way (MW) and Local Group, we know that H ii regions can be
extremely diverse in their nature, ranging from the least bright and sub-
parsec-sized regions (e.g. ultra-compact H ii regions observed within the MW;
see e.g. Hoare, 2005) to the most bright and ${\sim}400$ pc sized superbubbles
(e.g. resembling 30 Doradus in the Large Magellanic Cloud; see Oey, 1996;
Pellegrini et al., 2010). Studies of nebular LFs in nearby external galaxies
commonly probe H ii regions ionized by star clusters and associations while
only in a few cases does the sensitivity of the observations allow probing H
ii regions ionized by single OB stars (e.g. Azimlu et al., 2011).

Variations in the slope of the LF can unveil whether global properties of
galaxies, as well as local parameters such as chemical abundance, dust
content, gas dynamics (e.g. spiral arm perturbations) influence the star
formation process. Only a handful of studies have attempted to examine the
relation between the LF slope and global galaxy parameters mostly finding weak
or statistically insignificant correlations (Kennicutt et al., 1989; Elmegreen
& Salzer, 1999; Youngblood & Hunter, 1999; van Zee, 2000; Thilker et al.,
2002; Liu et al., 2013; Cook et al., 2016). On the other hand, studies of
individual star-forming galaxies with larger H ii region samples mainly
focused on variations of the LF between the inner and outer disk and between
the spiral arm and inter-arm areas. These studies showed that in the inter-arm
areas and in the outer disk the H ii region LF becomes steeper, however such
variations have not been found to be universal (see, e.g. Cepa & Beckman,
1990; Rand, 1992; Banfi et al., 1993; Knapen et al., 1993; Rozas et al., 1996;
Knapen, 1998; Lelièvre & Roy, 2000; Thilker et al., 2000; Helmboldt et al.,
2005; Gutiérrez et al., 2011; Scoville et al., 2001; Azimlu et al., 2011).
Quantifying variations of the LF slope with the galaxy’s SFR, morphology, or
local environmental properties can elucidate whether e.g. environment affects
the demographics of young stellar populations and may provide new clues for
our understanding of the star formation process.

New observations with integral field units (IFUs) are now allowing us to
characterize the physical properties of H ii regions in external galaxies in
greater detail, with higher sensitivity and spatial resolution, and
unprecedented statistics (see e.g. Rousseau-Nepton et al., 2018; Kreckel et
al., 2019; Espinosa-Ponce et al., 2020). The
PHANGS–MUSE111http://www.phangs.org data set (Emsellem et al., 2021) is now
starting to unveil its full potential in resolving and studying H ii region
properties and connecting them to galactic structure and galaxy evolution (see
e.g. Ho et al., 2019; Kreckel et al., 2019, 2020). With its 19 nearby star-
forming galaxies and the detection of tens of thousands of H ii regions at a
spatial resolution of about $70$ pc, the PHANGS–MUSE data are ideal to study
the H ii region LF, marking a turning point in terms of statistics and ability
to deblend/resolve H ii regions. Furthermore, the availability of PHANGS–ALMA
(Leroy et al., 2021) and PHANGS–HST (Lee et al., 2021) observations opens up
the possibility to compare the H ii region LFs to the MFs of giant molecular
clouds (GMCs) and young stellar regions (Wei et al., 2020; Rosolowsky et al.,
2021; Thilker et al., 2021; Whitmore et al., 2021) and to investigate how the
demographics of substructure changes over the course of the star formation and
feedback process.

This work is focused on the study of nebular LFs and their variations as a
function of galaxy global properties (e.g. galaxy morphology, star formation
rate, stellar mass, and gas-phase metallicity) and galactic environment (e.g.
spiral arm vs. inter-arm) across the entire PHANGS–MUSE sample. The paper is
structured as follows: Section 2 briefly describes the PHANGS–MUSE data set,
the data reduction strategy, and the ancillary data that have been used. In
Section 3, we describe the source identification technique and present our H
ii region catalogs. In Section 4, we present the nebular LF of our galaxies
and the fitting technique (Section 4.1), and we investigate variations of the
LF slope as a function of global properties across the sample (Section 4.2) or
for different populations of H ii regions within individual galaxies (Section
4.3). In Section 5, we assess the effects of completeness and crowding on our
results. In Section 6, we discuss our results and what drives variations in
the LF slope. We finally summarize our main findings in Section 7 and
highlight future prospects in Section 8.

##  2 Observations and data

For a complete and detailed description of the PHANGS–MUSE sample and data we
redirect the reader to Emsellem et al. (2021), while in this section we
summarize the overall properties of the sample and provide essential details
of the data reduction and analysis pipelines. In the last part of this
section, we also describe the ancillary data products that we use in this
work.

###  2.1 PHANGS–MUSE sample

The PHANGS–MUSE sample includes 19 nearby ($D<20$ Mpc) star-forming galaxies
(Sa–Sc Hubble morphological type) with relatively low inclination ($i\lesssim
55$ degrees). We summarize their main properties in Table 3. The galaxies are
located close to the so-called star formation main sequence and span a stellar
mass range of $\log(M_{\star}/M_{\odot})=9.4{-}11$. The galactocentric radii
used in this paper have been deprojected according to the inclinations and
position angles reported in Table 3, which are taken from Lang et al. (2020)
who performed modeling of the CO(2–1) kinematics using the PHANGS–ALMA data.
The area covered by the MUSE observations mainly samples the star-forming disk
of the galaxies. The mean maximum galactocentric radius covered across the
sample is $0.86R_{25}$, the exact coverage for each galaxy is reported in
Table 3.

###  2.2 PHANGS–MUSE data reduction and analysis

The PHANGS–MUSE data set (Emsellem et al., 2021) is the result of an extensive
(${\sim}170$ h) observational campaign (PI E. Schinnerer) with the MUSE IFU
(Bacon et al., 2010) mounted on the Unit 4 telescope (UT4) at ESO’s VLT, with
the addition of archival observations for NGC0628 and the centers of 5
galaxies in the sample. Each galaxy has been covered by several pointings
(from 3 up to 15) to obtain a contiguous mosaic of the star-forming disk. All
the observations have been performed in wide-field mode (WFM with
$1\arcmin\times 1\arcmin$ field of view), either in natural seeing (WFM-noAO,
for 11 galaxies) or taking advantage of the GALACSI adaptive optics system
(WFM-AO, for eight galaxies).

The data have been reduced using a version of the pymusepipe python
package,222Available from https://github.com/emsellem/pymusepipe a wrapper of
the esorex MUSE reduction recipes, tailored to the PHANGS–MUSE program
(Emsellem et al., 2021). Broad-band images taken with the ESO/MPG 2.2m Wide
Field Imager and with the Du Pont DirectCCD camera (Razza et al. in
preparation) were used as a reference to astrometrically align and flux
calibrate the MUSE data. The complete data reduction provides a sky-
subtracted, flux-calibrated, mosaicked data cube for each galaxy with a spaxel
size of $0.2\arcsec$. The point spread function (PSF) of the observations was
determined for each pointing, and for each galaxy the largest value (reported
in Table 3) was used to create PSF-homogenized data cubes. We refer to these
as the copt data cubes to distinguish them from the native data cubes
maintaining the native resolution of the observations.

The data analysis of the mosaicked data cubes has been performed with a data
analysis pipeline (DAP) based on the gist (Galaxy IFU Spectroscopy Tool;
Bittner et al. 2019) software package and tailored to the PHANGS–MUSE program
which provides maps of stellar kinematics, stellar population properties
(including stellar ages and mass surface density), and emission line fluxes
and kinematics (Emsellem et al., 2021). The DAP employs the penalized pixel
fitting method via the pPXF package (Cappellari, 2017) to extract information
on the stellar population and the ionized gas from the MUSE spectra over the
wavelength range $4850{-}7000$ Å taking into account the MUSE spectral
resolution as parametrized in Bacon et al. (2017). First, the stellar
kinematics and the stellar population information are extracted in two
subsequent steps using binned spectra with a continuum S/N target of 35, and
E-MILES simple stellar population models (Vazdekis et al., 2016) generated
with a Chabrier (2003) initial mass function, BaSTI isochrones (Pietrinferni
et al., 2004), eight ages ($0.15{-}14$ Gyr), and four metallicities
($\textrm{[Z/H]}=[-1.5,-0.35,0.06,0.4]$). Finally, the ionized gas analysis is
performed by leveraging the previous steps and simultaneously fitting the
stellar continuum and the emission lines using single spaxel spectra. Emission
lines are modeled with a single Gaussian function and emission line flux maps
are corrected for foreground Galactic extinction, using the O’Donnell (1994)
extinction law and the ${E(B-V)}$ values from Schlafly & Finkbeiner (2011). We
refer the reader to the PHANGS-MUSE survey paper (Emsellem et al., 2021) for
further details on the choice of models and the spectral fitting procedure. In
this paper, we will make use of the copt H$\alpha$ maps to localize H ii
regions and the native MUSE mosaicked data cubes to extract their integrated
spectra (see Section 3).

###  2.3 Environmental and foreground stars masks

Morphological masks delimiting stellar structures by using Spitzer NIR
$3.6\leavevmode\nobreak\ \mu$m imaging are available for the PHANGS–MUSE
galaxies from Querejeta et al. (2021). We use these masks to define galaxy
centers, bars, spiral arms, and inter-arms regions. In the work of Querejeta
et al. (2021), spiral arms have been defined only when they are dominant
features across the galaxy disk (omitting flocculent galaxies) by fitting a
log-spiral function to bright regions in the NIR images along each spiral arm.
These modeled log-spiral curves are then assigned an empirically determined
width based on the overlap with CO(2–1) emission from PHANGS–ALMA data and as
a last step the starting and ending azimuth of each spiral segment is adjusted
to define the final spiral arm footprint. The outer edges of bars have been
defined as ellipses where the bar length, ellipticity, and position angle come
from a compilation of measurements from the literature based on NIR imaging
(mostly from Herrera-Endoqui et al., 2015).

We make use of the Querejeta et al. (2021) environmental masks to define five
environments – namely galaxy center, bar, spiral arm, inter-arm, and disk –
according to the criteria described below. Our galaxy “centers” include small
bulges, inner star-forming rings or nuclei. We define “bars” in the same way
they are defined in the Querejeta et al. (2021) masks. For the galaxies with
spiral arms, we define as “spiral arms” the regions flagged as spiral arms
with the addition of the bar ends regions (when spiral arms overlap with bar
ends). All remaining regions are defined as “inter-arm” for the galaxies with
spiral arms and “disc” for the galaxies without spiral arms. These five
environments can be visualized for NGC4321 in Fig. 1 and for the remaining
galaxies in Appendix A, in Fig. 1 we also show a sketch of the different
environments to guide the reader.

There are a number of foreground stars across the field of view (FoV) of our
observations which need to be excluded from our analysis. Masks of foreground
stars for the PHANGS–MUSE galaxies are described in Emsellem et al. (2021).
They have been generated by matching the positions of Gaia point sources. To
avoid masking compact H ii regions and galactic nuclei, which may be included
as point sources in the Gaia catalog, a further check was performed to
identify the rest-frame Ca ii triplet absorption features at $8498$, $8542$,
and $8662$ Å in the MUSE spectra.

Table 1: Main properties of the galaxies in the PHANGS–MUSE sample.

Galaxy RA Dec Morph1 Distance2 $i^{3}$ PA3 Log $M_{*}^{4}$ SFR4
$\Delta_{\mathrm{SFMS}}^{4}$ $R_{25}^{1}$ $\mathrm{FWHM}_{\mathrm{PSF}}^{5}$
pc/′′ rmax ∘ ∘ Mpc ∘ ∘ $M_{\odot}$ $M_{\odot}\leavevmode\nobreak\
\mathrm{yr}^{-1}$ ′′ ′′ IC5332 353.61453 -36.10108 SABc 9.01 26.9 74.4 9.67
0.41 0.01 182 0.87 43.6 0.68 NGC0628 24.173855 15.783643 Sc 9.84 8.9 20.7
10.34 1.75 0.18 297 0.92 47.6 0.51 NGC1087 41.60492 -0.498717 Sc 15.85 42.9
359.1 9.93 1.31 0.33 89 0.92 76.5 1.37 NGC1300 49.920815 -19.411114 Sbc 18.99
31.8 278.0 10.62 1.17 -0.18 178 0.89 91.6 0.96 NGC1365 53.40152 -36.140404 Sb
19.57 55.4 201.1 10.99 16.90 0.72 361 1.15 94.4 0.73 NGC1385 54.369015
-24.501162 Sc 17.22 44.0 181.3 9.98 2.09 0.50 102 0.67 83.1 1.21 NGC1433
55.506195 -47.221943 SBa 18.63 28.6 199.7 10.87 1.13 -0.36 186 0.91 89.9 1.05
NGC1512 60.975574 -43.348724 Sa 18.83 42.5 261.9 10.71 1.28 -0.21 253 1.25
90.8 0.62 NGC1566 65.00159 -54.93801 SABb 17.69 29.5 214.7 10.78 4.54 0.29 217
0.8 85.3 0.73 NGC1672 71.42704 -59.247257 Sb 19.4 42.6 134.3 10.73 7.60 0.56
185 0.96 93.5 1.00 NGC2835 139.47044 -22.35468 Sc 12.22 41.3 1.0 10.00 1.24
0.26 192 1.15 59.0 0.69 NGC3351 160.99065 11.70367 Sb 9.96 45.1 193.2 10.36
1.32 0.05 217 1.05 48.2 0.76 NGC3627 170.06252 12.9915 Sb 11.32 57.3 173.1
10.83 3.84 0.19 308 1.05 54.7 0.66 NGC4254 184.7068 14.416412 Sc 13.1 34.4
68.1 10.42 3.07 0.37 151 0.89 63.3 1.28 NGC4303 185.47888 4.473744 Sbc 16.99
23.5 312.4 10.52 5.33 0.54 207 0.78 82.0 0.66 NGC4321 185.72887 15.822304 SABb
15.21 38.5 156.2 10.75 3.56 0.21 183 1.16 73.4 0.97 NGC4535 188.5846 8.197973
Sc 15.77 44.7 179.7 10.53 2.16 0.14 244 0.56 76.1 0.56 NGC5068 199.72807
-21.038744 Sc 5.2 35.7 342.4 9.40 0.28 0.02 224 1.04 25.2 0.89 NGC7496
347.44702 -43.42785 Sb 18.72 35.9 193.7 10.00 2.26 0.53 100 0.89 90.3 1.06

1 HyperLEDA, Makarov et al. (2014), 2 Anand et al. (2021), 3 Lang et al.
(2020), 4 Leroy et al. (2021), 5 Emsellem et al. (2021)

333The table reports the galaxy name (col 1), galaxy center sky coordinates
(col 2–3), Hubble morphological type (col 4), distance (col 5), inclination
(col 6), position angle (col 7), total stellar mass (col 8), total star
formation rate (col 9), offset from the star formation main sequence (col 10),
$R_{25}$ $B$-band isophotal radius (col 11), angular resolution of the copt
data (col 12), the parsec per arcsec scaling factor (col 13), and the maximum
deprojected galactocentric radius in units of $R_{25}$ (col 14).

##  3 Data analysis

###  3.1 Ionized nebulae identification

To isolate H ii regions across our sample, we make use of the PHANGS–MUSE copt
H$\alpha$ emission line maps provided by the DAP procedure summarized in Sec.
2.2. We chose the copt rather than the native resolution maps to mitigate the
effect of PSF spatial variations across the FoV of single targets. Ionized
nebulae are identified using an adaptation of the HIIphot software (Thilker et
al., 2000) originally implemented to work with H$\alpha$ narrow-band images.
The software first detects what are called “seed regions” and then grows them
up to a given contrast (i.e. the termination gradient) that is selected by the
user. We set the background level of the H$\alpha$ maps to the mean H$\alpha$
flux of the pixels with surface brightness below
$\Sigma_{\mathrm{H\alpha}}<1\\!\times\\!10^{-17}$ erg s-1 arcsec-2 and the
detection threshold to $3\sigma$ above the background, with $\sigma$ being the
standard deviation of the background pixels. Once the seed regions have been
identified, we perform further cleaning based on visual inspection to avoid
artifacts due to noise in the H$\alpha$ map. We consider the total integrated
flux of the seed region and impose a S/N cut of $50$ (using the H$\alpha$
emission line error maps to assess the noise) and a
$\Sigma_{\mathrm{H\alpha}}$ cut above $3\sigma$ of the background. To avoid
the detection of regions with unphysical sizes, we find that the best solution
is to limit the spatial smoothing operated by HIIphot on the input maps for
the iterative search of seed regions to three iterations, each time increasing
the starting PSF by $10$%.

Once the seed regions have been identified, HIIphot starts to grow them.
Different methods have been used in the literature to define H ii region
boundaries: emission line ratios, H$\alpha$ equivalent width, H$\alpha$
spatial gradients, or a combination thereof (see e.g. Blanc et al., 2009;
Thilker et al., 2000; Sánchez et al., 2012; Zhang et al., 2017). In this work,
we rely on HIIphot, which uses a user-selected termination gradient of the
H$\alpha$ surface brightness in units of emission measure (EM), to stop the
growth of seed regions – with a lower termination gradient the seed regions
are allowed to grow more and vice versa (see Thilker et al., 2000, for an
example). In principle, in the classical scenario of an H ii region
represented as a Strömgren sphere, closer galaxies and/or better seeing
conditions lead to the detection of steeper gradients and vice versa (see e.g.
Oey et al., 2007). Considering the distances of the galaxies in our sample
($D<20$ Mpc) and the seeing of our observations
($\mathrm{FWHM}_{\mathrm{PSF}}\leq 1.2\arcsec$), after visual inspection we
select a termination gradient of $2.43\\!\times\\!10^{-16}$ erg s-1 arcsec-2
pc-1 (i.e. corresponding to $5$ EM pc-1) to also ensure coherence with results
presented in the literature (see Oey et al., 2007). In general a higher/lower
termination gradient results in a steeper/shallower LF. We checked that
varying the termination gradient in the interval $1.5{-}10$ EM pc-1 gives LF
with slopes that are, within the errors, fully compatible with what we discuss
in Sec. 4.1 of this paper. Excluding the regions overlapping with foreground
stars (see Sec. 2.3), we detect a total number of $31\,399$ ionized nebulae
across the 19 galaxies of the PHANGS–MUSE sample. We use the nebulae spatial
masks for each galaxy to extract the nebula integrated spectra from the MUSE
native mosaicked data cubes. These spectra are then fitted using the DAP of
the PHANGS–MUSE data as described in Sec. 2.2. The only difference at this
level is that we fit the MUSE spectra across their entire wavelength range so
to include also the modeling of the [S iii]$\lambda 9069$Å emission line. It
is worth noting that in this work we do not correct the emission line fluxes
for the diffuse background (i.e. the emission of the diffuse ionized gas DIG
component permeating star-forming disks, see Haffner et al., 2009). However,
we checked that such a correction starts to be more relevant only for the
H$\alpha$ emission of the fainter H ii regions and, therefore, does not affect
the LF slopes presented in Sect. 4. Depending on the galaxy, the H$\alpha$
emission outside the ionized nebulae footprints can represents $20{-}50\%$ of
the total H$\alpha$ emission (see Belfiore et al., 2021).

###  3.2 The final H ii region catalog

Figure 1: H ii regions and environments for NGC4321. The figure shows the
H$\alpha$ emission in the background, color coded according to the color
scheme on the right, with overlaid the borders of the H ii regions in our
catalog. The centers of the nebulae which have been discarded by our selection
criteria are marked with crosses. In the lower left corner, the black circle
indicates the PSF of the MUSE observations while the black line marks a
physical scale corresponding to $1$ kpc. Both the H ii regions and the
discarded nebulae are color coded according to our definition of environments
as outlined by the color scheme at the bottom and shown in the bottom-right
sketch.

In order to obtain our fiducial H ii regions catalog we exclude:

  * •

Nebulae located in galaxy centers (as defined in Sec. 2.3); this is meant to
avoid areas with a high SFR surface density (e.g. starburst rings in NGC1300,
NGC1512, NGC1672, NGC3351, NGC4303, and NGC4321), where the deblending of H ii
regions is not possible at our resolution (i.e. typically few parsec in
radius; see e.g. Barnes et al., 2020), and contamination by active galactic
nuclei (AGN) (i.e. NGC1365, NGC1512, NGC1566, and NGC1672 show evidence of
hosting low-luminosity AGN when looking at emission line maps). This cut
comprises about $1.5$% of the nebulae.

  * •

Nebulae located outside the “H ii region” area in at least one of the
classical Baldwin–Phillips–Terlevich (BPT) diagrams (Baldwin et al., 1981);
more specifically, these are the regions located above the Kauffmann et al.
(2003) line in the [O iii]/H$\beta$ vs. [N ii]/H$\alpha$ diagnostic, and above
the Kewley et al. (2001) line in the [O iii]/H$\beta$ vs. [S ii]/H$\alpha$, as
well as the [O iii]/H$\beta$ vs. [O i]/H$\alpha$ diagnostic. To build the BPT
diagrams we require a $\mathrm{S/N}>3$ for the H$\beta$, [O iii]$\lambda
5007$Å, [O i]$\lambda 6300$Å, H$\alpha$, [N ii]$\lambda 6584$Å, [S ii]$\lambda
6716$Å, and [S ii]$\lambda 6731$Å emission lines; this is meant to ensure a
robust estimate of both the dust correction based on the H$\alpha$/H$\beta$
Balmer decrement and the quantities relying on line ratios (e.g. ionization
parameter and gas-phase metallicity) used in this work. This cut comprises
about $23$% of the nebulae.

  * •

Nebulae with velocity dispersion $\sigma>400$ km s-1 for all three emission
line groups (i.e. hydrogen, high ionization and low ionization lines) fitted
separately by the DAP (see Emsellem et al., 2021) in order to avoid cases
where the emission line fitting performs poorly and to exclude objects that
are potentially supernovae remnant (SNR) candidates. This cut comprises about
$1$% of the nebulae.

  * •

Regions whose geometric center matches within $0.5\arcsec$ the positions of
planetary nebulae (PNe) identified using [O iii] emission line maps
(Scheuermann et al. in preparation). This cut comprises about $0.6$% of the
nebulae.

  * •

Regions where the distance of the geometric center from the edges of the MUSE
mosaic FoV is less than one $\mathrm{FWHM}_{\mathrm{PSF}}$. This cut comprises
about $2$% of the nebulae.

We note that the S/N of the H$\alpha$ and H$\beta$ lines is usually high, more
specifically, $99$% of the region spectra have a $\mathrm{S/N}\gtrsim 32$ for
the H$\alpha$ line and ${\gtrsim}12$ for the H$\beta$ line. All together, the
cuts remove about $26$% of the original nebulae and leave us with a final
joined catalog of $23\,301$ H ii regions. A discussion on how selection
effects may influence our results is presented in Appendix C.

In Fig. 1, we show the spatial masks of the H ii regions found in NGC4321,
highlighting also the different galactic environments defined in Sec. 2.3.
Analogous figures for the remaining galaxies in our sample are presented in
Appendix A.

We assume a screen geometry and use pyneb444https://pypi.org/project/PyNeb/
(Luridiana et al., 2015) to correct line fluxes for dust extinction via the
H$\alpha$/H$\beta$ ratio, adopting the O’Donnell (1994) reddening law with
$R_{V}=3.1$ and a theoretical H$\alpha$/H$\beta=2.86$. The extinction
corrected emission line luminosities of the H ii regions are then computed
using the distances reported in Table 3. For every H ii region in our catalog,
we also estimate the gas-phase metallicity and the gas ionization parameter by
using extinction corrected emission line fluxes. The gas-phase metallicity
$\mathrm{O/H}$ is calculated using the Pilyugin & Grebel (2016) S-calibration
(Scal hereafter). This calibration relies on three diagnostic line ratios
(i.e. [N ii]/H$\beta$, [S ii]/H$\beta$, and [O iii]/H$\beta$) and provides an
empirical calibration against H ii regions that have direct constraints on
their nebular temperatures and hence their abundances. This calibration is
relatively insensitive to changes in gas pressure or ionization parameter and
we adopted it as fiducial approach in this paper (see Kreckel et al., 2019,
for a discussion). In addition, for each galaxy we fit the radial metallicity
gradient by using an unweighted least-square linear fitting of the trend
between $12+\log(\mathrm{O/H})$ and the deprojected galactocentric radius (see
Fig. 30 in Appendix A). The gas ionization parameter represents the ratio
between the ionizing photon flux density and the gas hydrogen density. In this
paper, we express the ionization parameter as $q=U\times c=Q(H^{0})/4\pi
R^{2}n$, where $c$ is the speed of light, $U$ is the dimensionless ionization
parameter, $Q(H^{0})$ is the number of hydrogen ionizing photons ($E>13.6$ eV)
emitted per second, $R$ is the empty (wind-blown) radius of the H ii region,
and $n$ its hydrogen density. The ionization parameter is ultimately defined
by the structure of an H ii region (e.g. size, gas density, filling factor)
and the properties of its ionizing source. Photoionization models show that it
can be extracted via different diagnostic lines (Kewley & Dopita, 2002; Dors
et al., 2011). In this paper, we use the calibration proposed by Diaz et al.
(1991) based on the [S iii](9069+9532)/[S ii](6717+6713) line ratio. As the [S
iii]$\lambda 9532$Å line falls outside the wavelength range covered by MUSE,
we assume that $\text{{[S\,{iii}]}}\lambda
9532$Å$=2.47\text{{[S\,{iii}]}}\lambda 9069$Å according to default atomic data
in pyneb (Luridiana et al., 2015). It should be noted that for about $3000$ H
ii regions we are not able to estimate the ionization parameter due to lack of
detection of the [S iii]$\lambda 9532$Å emission line.

The H ii regions catalog containing the host galaxy name, the sky coordinates
of the geometrical center, the deprojected galactocentric distance from the
host galaxy center, the H$\alpha$ and H$\beta$ observed fluxes, the H$\alpha$
extinction corrected flux, the ionization parameter, and an environmental flag
for all the H ii regions discussed in the current paper is only available in
electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr
(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.

##  4 Results

###  4.1 The H ii region luminosity function

Figure 2: Best-fitting models of the H ii region LF for the galaxies of our
sample. Galaxies are ordered by increasing stellar mass from top left to
bottom right and their names are indicated within each panel. The dotted black
line indicates the empirical PDF, plotted as $p(L)$ versus the logarithm of
the H ii regions H$\alpha$ luminosity expressed in erg s-1. The vertical
dashed line indicates the estimated completeness limit. The vertical solid
black line and the shaded gray area indicate the best-fitting
$L_{\mathrm{min}}$ and its uncertainty. The LF best-fitting model is marked by
the solid red line while the associated slope $\alpha$ and its uncertainty are
indicated in the upper right corner of each panel.

We make use of the extinction corrected H$\alpha$ luminosities of our H ii
regions to obtain the nebular luminosity function (LF) of each galaxy in our
sample. H ii region LFs have long been known to be well-described by a power-
law (Kennicutt et al., 1989). The key parameter of the LF is its slope
$\alpha$. Several different methods for deriving this from the empirical data
have been employed. Past studies have often used least-squares linear
regression of binned data (i.e. histograms) in log–log space. This method,
despite being widespread, can give incorrect results under relatively common
conditions, as discussed in Clauset et al. (2009). In the context of LFs, this
is mainly related to the fact that linear regression assumes Gaussian noise in
the dependent variable while the noise of the logarithm of a histogram is not
Gaussian. Furthermore, the choice of the binning scheme (e.g. equal-luminosity
bins versus equal-number bins) and the definition of the bin center can affect
the LF modeling and add uncertainties to $\alpha$ (see e.g. Maíz Apellániz &
Úbeda, 2005; Cook et al., 2016).

Table 2: Properties of the H ii region LF in the PHANGS–MUSE galaxies.

Galaxy | Slope $\alpha$ | $\log L_{\mathrm{min}}$ | $\log L_{\mathrm{compl}}$ | O7V eq | $\log L_{\mathrm{max}}$ | $\log L_{\mathrm{tot}}$ | $N_{\mathrm{tot}}$ | $N_{\mathrm{LFfit}}$ | $N_{\mathrm{DB}}$ | $\mathrm{FWHM}_{\mathrm{seeing}}$
---|---|---|---|---|---|---|---|---|---|---
|  | $\mathrm{erg\leavevmode\nobreak\ s^{-1}}$ | $\mathrm{erg\leavevmode\nobreak\ s^{-1}}$ |  | $\mathrm{erg\leavevmode\nobreak\ s^{-1}}$ | $\mathrm{erg\leavevmode\nobreak\ s^{-1}}$ |  |  |  | $\mathrm{pc}$
NGC5068 | 1.58$\pm$0.09 | 35.92$\pm$0.06 | 35.86 | 0.06 | 39.26 | 40.38 | 1465 | 1206 | 9 | 26.18
IC5332 | 1.82$\pm$0.13 | 36.25$\pm$0.19 | 35.71 | 0.12 | 38.62 | 39.71 | 612 | 357 | 1 | 37.91
NGC1087 | 1.71$\pm$0.10 | 37.49$\pm$0.10 | 36.88 | 2.04 | 40.01 | 41.24 | 891 | 439 | 101 | 70.38
NGC1385 | 1.52$\pm$0.08 | 37.29$\pm$0.17 | 36.68 | 1.29 | 40.11 | 41.47 | 914 | 556 | 157 | 55.67
NGC7496 | 1.69$\pm$0.10 | 37.16$\pm$0.09 | 36.64 | 0.95 | 39.70 | 40.75 | 547 | 300 | 39 | 80.35
NGC2835 | 1.76$\pm$0.11 | 37.08$\pm$0.10 | 36.67 | 0.78 | 39.56 | 40.72 | 819 | 432 | 35 | 67.90
NGC0628 | 1.71$\pm$0.10 | 36.63$\pm$0.08 | 36.24 | 0.28 | 39.34 | 40.87 | 2230 | 1273 | 42 | 43.77
NGC3351 | 1.98$\pm$0.16 | 36.84$\pm$0.13 | 36.06 | 0.45 | 38.77 | 40.14 | 784 | 369 | 3 | 50.56
NGC4254 | 1.61$\pm$0.08 | 37.36$\pm$0.12 | 36.77 | 1.49 | 39.93 | 41.62 | 2536 | 1430 | 333 | 56.32
NGC4303 | 1.56$\pm$0.09 | 37.26$\pm$0.15 | 36.92 | 1.18 | 40.03 | 41.75 | 2208 | 1567 | 353 | 63.94
NGC4535 | 1.63$\pm$0.09 | 36.79$\pm$0.09 | 36.13 | 0.40 | 39.87 | 41.00 | 1444 | 934 | 65 | 42.63
NGC1300 | 1.77$\pm$0.12 | 37.08$\pm$0.25 | 36.26 | 0.79 | 39.42 | 40.84 | 1169 | 632 | 49 | 81.50
NGC1512 | 2.04$\pm$0.17 | 37.41$\pm$0.10 | 36.48 | 1.71 | 39.07 | 40.43 | 472 | 209 | 15 | 113.51
NGC1672 | 1.67$\pm$0.10 | 37.55$\pm$0.12 | 37.03 | 2.32 | 40.03 | 41.40 | 1051 | 570 | 152 | 89.80
NGC4321 | 1.95$\pm$0.15 | 37.95$\pm$0.27 | 36.93 | 5.82 | 39.78 | 41.31 | 1385 | 382 | 111 | 85.17
NGC1566 | 1.57$\pm$0.09 | 37.14$\pm$0.14 | 36.77 | 0.91 | 40.53 | 41.58 | 1655 | 1026 | 186 | 68.27
NGC3627 | 1.57$\pm$0.09 | 37.43$\pm$0.23 | 36.84 | 1.78 | 40.10 | 41.56 | 1007 | 657 | 188 | 57.44
NGC1433 | 1.96$\pm$0.16 | 37.32$\pm$0.14 | 36.20 | 1.37 | 39.41 | 40.69 | 1258 | 407 | 28 | 81.76
NGC1365 | 1.74$\pm$0.11 | 37.82$\pm$0.20 | 37.04 | 4.31 | 40.76 | 41.44 | 854 | 347 | 119 | 108.51

555The table reports the galaxy name (col 1), the LF slope (col 2) and
$L_{\mathrm{min}}$ (col 3), the LF completeness limit (col 4), the number of
equivalent O7 V stars corresponding to $L_{\mathrm{min}}$ (col 5), the
luminosity of the brightest H ii region (col 6), the integrated H ii regions
luminosity (col 7), the total number of H ii regions (col 8), the number of H
ii regions used to fit the LF slope (col 9), the number of “density bounded” H
ii regions with $L_{\mathrm{H\alpha}}>10^{38.6}$ erg s-1 (col 10), and the
physical resolution of the H$\alpha$ maps used to detect H ii regions in
parsecs (col 11).

The slope of the LF is commonly constrained by fitting the bright end of the
LF or, more specifically, H ii regions with a luminosity above a given
$L_{\mathrm{min}}$. This is the luminosity below which the LF starts to
flatten – it is often referred to as the “turnover point” and, depending on
the sensitivity of the observations, can arise due to incompleteness of the
data. $L_{\mathrm{min}}$ is an important parameter for determining the LF
slope and in many studies it has been, quite arbitrarily, fixed to the
luminosity where the histogram of the LF peaks or, alternatively, where the
bin count number steadily decreases (e.g. Cook et al., 2016; Azimlu et al.,
2011). As described in Clauset et al. (2009), by adopting a maximum likelihood
estimation (MLE) method, in combination with the Kolmogorov–Smirnov (KS)
statistic, it is possible to obtain both $L_{\mathrm{min}}$ and $\alpha$ of a
given LF in a statistically robust way (see Sec. 5 for further discussion).

In this paper, we adopt this method to fit an empirical LF: it does not depend
on any binning scheme and, as discussed in Clauset et al. (2009), provides a
statistically robust way to fit data following heavily tailed distributions.
Whitmore et al. (2014) compared the performance of this method to classical
ones based on histograms by fitting the LF of star clusters observed with HST
and found an overall good agreement, with MLE fits giving steeper slopes for
steeper LFs. To perform the LF fits, we make use of the Python package
powerlaw (Alstott et al., 2014) which follows the prescription given in
Clauset et al. (2009) and models the probability distribution function (PDF)
$p(L)$ connected to the empirical LF using a power-law of the form:

$p(L)=(\alpha-1)\,L_{\mathrm{min}}^{\alpha-1}\,L^{-\alpha}\quad\text{with}\quad
L\geq L_{\mathrm{min}}\leavevmode\nobreak\ .$ (1)

The algorithm performs an MLE fit by recursively fixing $L_{\mathrm{min}}$
equal to each empirical data point. For a given $L_{\mathrm{min}}$, the
algorithm tries a range of different slopes and selects as the best-fitting
slope the one maximizing the likelihood estimator. At this point we have a set
of models that includes the best-fitting model for each fixed
$L_{\mathrm{min}}$. The KS statistic is then used to compute the maximum
distance between this set of models and the empirical cumulative distribution
function. The final best-fitting model is chosen as the one that minimizes the
KS statistic.

We limit the search of $L_{\mathrm{min}}$ to an interval of $\pm 1\sigma$
around the median of the distribution (i.e. where, visually, the LFs start to
flatten) and the LF slope to the interval of $\alpha=1{-}3$. The error on
$L_{\mathrm{min}}$ is taken to be the standard deviation of the
$L_{\mathrm{min}}$ values found by rerunning the fitting procedure on $1000$
mock LFs extracted via bootstrapping from a given LF, after making sure that
the error converges to a stable value as the number of bootstraps approaches
1000. For the error on $\alpha$, we consider the trend of the likelihood
estimator as a function of $\alpha$ for the models with the best
$L_{\mathrm{min}}$. This trend can be well represented by a Gaussian that
peaks at the best $\alpha$ and we compute the error on $\alpha$ as the
standard deviation of the best-fitting Gaussian function modeling the data
(similar to what has been done by Whitmore et al., 2014). We do not run any
test to check for the existence of an upper cutoff in the LFs.

The best-fitting models of our LF are shown in Fig. 2 and reported in Table 5.
The slopes we find are in the interval of $\alpha=1.5{-}2$ with a mean value
of $1.73$ and a $\sigma$ of $0.15$, in agreement with what has been discussed
in the literature for Sa–Sc spiral galaxies (see e.g. Kennicutt et al., 1989;
Elmegreen & Salzer, 1999; Whitmore et al., 2014). We note that for our sample,
the minimum number of H ii regions used to fit the LF is $300$ (i.e. for
NGC7496) and that our fits can be considered robust against biases due to low
number statistics (see Sec. 5 for further discussion). In Fig. 2, we also mark
the completeness limit of our LFs that will be discussed in Sec. 5. In
Appendix C, we discuss how the results would change if H$\alpha$ luminosities
were not corrected for dust extinction (e.g. in narrow-band imaging). We find
that neglecting the dust extinction correction causes a slight steepening of
the LFs that, within their errors, remain compatible with the slopes we report
in this section.

Before describing further the details of our data analysis, it is worth
discussing the overall properties of our LFs. As can be seen in Fig. 2, our H
ii regions span a luminosity range of $\log(L_{\mathrm{H\alpha}}\
[\mathrm{erg\leavevmode\nobreak\ s^{-1}}])\sim 35{-}40$. H ii regions with
luminosities less than a few times $10^{37}$ erg s-1 are typically ionized by
single O- or B-type stars, whereas at higher luminosities and up to about
$10^{39}$ erg s-1 we enter the regime where H ii regions are ionized by
stellar associations or clusters containing multiple OB stars. In this regime,
the H$\alpha$ luminosity of an H ii region is expected to be roughly
proportional to its number of ionizing stars. H ii regions with luminosities
higher than a few $10^{39}$ erg s-1 are referred to as giant or supergiant H
ii regions (e.g. 30 Doradus in the LMC) and are typical of late-type galaxies
(Sc and Irr Hubble types; see e.g. Kennicutt et al., 1989; Elmegreen et al.,
1996).

The ionizing luminosities corresponding to $L_{\mathrm{min}}$ in each galaxy
of our sample are reported in Table 5 in terms of the number of equivalent O7
V type stars, assuming that the number of ionizing photons per second for an
O7 V star along the zero-age main sequence (ZAMS) is
$\log(Q_{0}^{\mathrm{O7\,V}}\ [\mathrm{photons\leavevmode\nobreak\
s^{-1}}])=49.05$ following Vacca (1994). The number of equivalent O7 V stars
is used here as a first order indicator of the number of ionizing stars in an
H ii region and at $L_{\mathrm{min}}$ it is in the range of $0.05{-}6$ in our
sample. This testifies that our LF fits, extending ${\sim}2{-}3$ orders of
magnitude above $L_{\mathrm{min}}$, are mainly probing H ii regions ionized by
star clusters and associations, and only for a few targets (i.e. NGC5068,
IC3352 with smaller distances) are we probing the regime of H ii regions
ionized by single massive stars.

Different studies in the literature argue that the LF steepens at luminosities
higher than $\log(L_{\mathrm{H\alpha}}\ [\mathrm{erg\leavevmode\nobreak\
s^{-1}}])=38.6\pm 0.1$ (i.e. the so-called LF cutoff or break) and is better
reproduced by a double power-law (i.e. type II LF, reported in individual
galaxy studies; e.g. Kennicutt et al., 1989; Rand, 1992; Rozas et al., 1996;
Beckman et al., 2000; Thilker et al., 2000; Gutiérrez et al., 2011; Lee et
al., 2011). Beckman et al. (2000) also observed a local sharp peak (adopting
their nomenclature, we refer to this feature as “glitch”) in the LF of their
galaxies at $\log(L_{\mathrm{H\alpha}}\ [\mathrm{erg\leavevmode\nobreak\
s^{-1}}])=38.6$ and suggested that it marks the transition of H ii regions
from ionization bounding (the central star/star cluster only ionizes the gas
within the H ii region) at low luminosities to density bounding (a large
amount of Lyman continuum photons from the central ionizing source escapes the
H ii region and ionizes the diffuse surrounding medium) at higher
luminosities.

The number of putative density bounded H ii regions in our galaxies (i.e.
$L_{\mathrm{H\alpha}}\geq 10^{38.6}$ erg s-1) varies from a few up to about
$350$ and is reported in Table 5. Neither the number of these regions nor
their fraction with respect to the total number of H ii regions correlates
with galaxy morphological type (two targets with Sa morphology have a few tens
of these regions while Sb and Sc galaxies show comparable numbers and scatter)
or the physical resolution of the observations. This suggests that in our
sample the detection of this kind of H ii region is not purely an effect of
blending due to limited spatial resolution (as suggested by the study of Lee
et al., 2011).

Figure 3: Best-fitting model of the upper end of the H ii region LF for the
sub-sample of nine galaxies with more than $100$ density bounded H ii regions.
Galaxies are ordered by increasing stellar mass from top left to bottom right
and their names are indicated within each panel. The dotted black line
indicates the empirical PDF, plotted as $p(L)$, versus the logarithm of the H
ii regions H$\alpha$ luminosity expressed in erg s-1. The vertical dashed line
indicates the estimated completeness limit. The vertical solid black line and
the solid red line indicate the best-fitting $L_{\mathrm{min}}$ and $\alpha$
of the global LF as shown in Fig. 2. The orange solid line marks the best-
fitting model for the regions with $L_{\mathrm{H\alpha}}\geq 10^{38.5}$ erg
s-1. The LF slopes of the two best-fitting models are indicated in the upper
right corner of each panel following the same color coding as the models.

For the nine galaxies in which we detect more than $100$ density bounded H ii
regions, we perform an additional fit of the LF only at
$L_{\mathrm{H\alpha}}\geq 10^{38.5}$ erg s-1 and show the results in Fig. 3.
We find that for these galaxies the slope of the upper end is indeed steeper
with respect to the slope measured globally; only for two of the nine targets
(i.e. NGC1365 and NGC4321) the variations are within the errors. However,
looking at the LFs of our sample in Fig. 2, the steepening at
$\log(L_{\mathrm{H\alpha}}\ [\mathrm{erg\leavevmode\nobreak\ s^{-1}}])>38.6\pm
0.1$ and the glitch at $\log(L_{\mathrm{H\alpha}}\
[\mathrm{erg\leavevmode\nobreak\ s^{-1}}])=38.6\pm 0.1$ are subtle features,
contrary to what has been shown by Beckman et al. (2000). Assessing which LFs
are better described by a double power-law is behind the scope of the current
paper. Overall, a single power-law fit gives a good representation of the
empirical data in our sample and is it thus adopted as our fiducial LF model.
In the following, we will take the advantage of this simpler parametrization
of the LF to directly compare the LF slopes of different galaxies and of H ii
region sub-samples within individual galaxies.

###  4.2 LF variations among galaxies

Figure 4: Trends of the LF slope $\alpha$ with global galaxy properties. From
top left to bottom right we show $\alpha$ as a function of the galaxy
morphological T class, the Scal metallicity at the median galactocentric
radius covered by our observations $12+\log(\mathrm{O/H})_{\langle
r_{\mathrm{gal}}\rangle}$, the total stellar mass $M_{*}$, the total star
formation rate SFR, the specific star formation rate sSFR, the SFR surface
density $\Sigma_{\mathrm{SFR}}$, the offset from the resolved star formation
main sequence $\Delta$ rSFMS, the offset from the resolved Kennicutt–Schmidt
relation $\Delta$ rKS, and the offset from the resolved molecular gas main
sequence $\Delta$ rMGMS. For each panel we report the Spearman correlation
coefficient $\rho$ and $p$-value of the plotted quantities. In the middle
right panel the solid red line shows the best-fitting relation found by Cook
et al. (2016) while the dashed gray and solid back horizontal lines
respectively mark the $\alpha$ of the Cook et al. (2016) line at
$\Sigma_{\mathrm{SFR}}$$=-2.6$ and the mean $\alpha$ of the points at
$\Sigma_{\mathrm{SFR}}$$>-2.6$ (see Sec. 6.1 for further details).

To look at what may drive galaxy-to-galaxy variations of the LF slope in our
sample, we verified that $\alpha$ does not have any dependence on galaxy
inclination reported in Table 3 and then estimate a few global galaxy
parameters, excluding the areas defined as galaxy centers (see Sec. 2.3), to
be in line with the H ii regions probed by our LFs. More specifically, we
estimate:

  * •

The gas-phase metallicity at a representative radius, taken to be the median
galactocentric radius of the H ii regions detected in each galaxy. The
metallicity at this radius is calculated from a linear fit to the radial
metallicity gradient (shown in Fig. 30 of Appendix B).

  * •

The total stellar mass ($M_{*}$) extracted by integrating the PHANGS–MUSE
stellar mass maps produced by the DAP.

  * •

The total star formation rate (SFR) extracted by integrating the PHANGS–MUSE
H$\alpha$ maps (not corrected for dust extinction) and applying the
calibration reported in Kennicutt & Evans (2012b) taken from the work of Hao
et al. (2011) and Murphy et al. (2011), who adopted a Kroupa & Weidner (2003)
IMF, with a Salpeter slope ($\alpha_{*}=-2.35$) from $1$ to
$100\leavevmode\nobreak\ M_{\odot}$ and $\alpha_{*}=-1.3$ for
$0.1{-}1\leavevmode\nobreak\ M_{\odot}$, and solar metallicity.

  * •

The specific SFR ($\mathrm{sSFR}=\mathrm{SFR}/M_{*}$) and the SFR surface
density ($\Sigma_{\mathrm{SFR}}$) obtained by dividing the SFR by $M_{*}$ and
the area covered by the MUSE FoV, respectively.

Taking advantage of the study on the resolved star formation scaling relations
at ${\sim}100$ pc scales carried out by Pessa et al. (2021), we also consider
the offset $\Delta$ of individual galaxies with respect to the trend of the
full sample for the resolved star formation main sequence (rSFMS;
$\Sigma_{\mathrm{M_{*}}}$ vs. $\Sigma_{\mathrm{SFR}}$), the Kennicutt–Schmidt
relation (rKS; $\Sigma_{\mathrm{H_{2}}}$ vs. $\Sigma_{\mathrm{SFR}}$), and the
molecular gas main sequence (rMGMS; $\Sigma_{\mathrm{M_{*}}}$ vs.
$\Sigma_{\mathrm{H_{2}}}$). The offsets $\Delta$ rSFMS, $\Delta$ rKS, and
$\Delta$ rMGMS have been estimated taking as a reference the modeled trends
for the star formation scaling relations of the whole sample and fitting the
same trends for each galaxy. The slope of the individual galaxy trends, for
each scaling relation, was set to be the same as the one of the entire sample,
and the offset of each galaxy from the global relation was estimated as the
difference between the normalization factors of the linear fits (see Pessa et
al., 2021, for additional details). It should be noted that these measurements
are not available for NGC0628 due to the quality of the MUSE data while, due
the lack of significant CO(2–1) emission, IC5332 does not have measurements of
$\Delta$ rKS and $\Delta$ rMGMS. The $\Delta$ rSFMS, $\Delta$ rKS, and
$\Delta$ rMGMS offsets are interesting quantities because, by construction of
the star formation scaling relations, they measure the relative variations of
the sSFR ($\Sigma_{\mathrm{SFR}}$/$\Sigma_{\mathrm{M_{*}}}$), the gas
depletion time ($t_{\mathrm{dep}}=$
$\Sigma_{\mathrm{H_{2}}}$/$\Sigma_{\mathrm{SFR}}$), and the molecular gas
fraction ($\Sigma_{\mathrm{H_{2}}}$/$\Sigma_{\mathrm{M_{*}}}$) within our
sample. We verified that estimating such offsets by excluding galaxy centers
has no significant effect on our findings.

The trends between the LF slope $\alpha$ and the aforementioned parameters,
with the addition of the morphological T type, are shown in Fig. 4 together
with their Spearman correlation coefficient $\rho$ and $p$-value, indicating
the probability that the two sets of data are uncorrelated. We summarize the
properties for which we look for a correlation in Table 7. We define a
correlation to be negligible when $\lvert\rho\rvert=[0{-}0.2]$, weak when
$\lvert\rho\rvert=[0.2{-}0.4]$, moderate when $\lvert\rho\rvert=[0.4{-}0.6]$,
strong when $\lvert\rho\rvert=[0.6{-}0.8]$, and very strong when
$\lvert\rho\rvert=[0.8{-}1]$; using the $p$-value to evaluate the probability
that, despite showing a correlation, two variables may be uncorrelated. It
should be noted that only a handful of studies so far have looked at the
correlation between $\alpha$ and global galaxy properties: Kennicutt et al.
(1989); Elmegreen & Salzer (1999); Youngblood & Hunter (1999); van Zee (2000),
and Thilker et al. (2002) investigated nebular LFs as in this paper, while Liu
et al. (2013) identified H ii regions via Pa$\alpha$, and Cook et al. (2016)
studied the GALEX far-ultraviolet (FUV) LFs of H ii regions. While the sample
of Cook et al. (2016) includes a few hundred galaxies, the other studies are
based on samples ranging from $10$ to $35$ galaxies, similar to our study. In
this section and in Sec. 6.1, we compare our results to those studies which,
like in our case, applied a uniform analysis methodology on galaxy samples. It
should be noted that using different tracers means probing different source
ages and, as reported by Oey & Clarke (1998a), older H ii regions tend to have
steeper LF slopes, mainly due to the short main-sequence lifetimes of the more
massive stars constituting the brighter H ii regions. This is the reason why
e.g. FUV observations, probing H ii regions with ages less than $100$ Myr, are
expected to deliver steeper LF compared to H$\alpha$ observations, typically
probing H ii regions younger than 10 Myr, and our comparison remains
qualitative.

We find a weak correlation between $\alpha$ and the galaxy morphology (Fig.
4a), and a negligible correlation with the gas-phase metallicity (Fig. 4b) and
the total stellar mass of our galaxies (Fig. 4c). Kennicutt et al. (1989) and
Elmegreen & Salzer (1999) found that the LF slope flattens for later-type
galaxies where the relative number of giant H ii regions is noticeably higher.
On the other hand, Cook et al. (2016) found no trend between $\alpha$ and
morphology when studying a sample more representative for irregulars, arguing
that the trends reported in the literature were relying on few data points in
this regime. Due to the selection of the PHANGS–MUSE sample we are only
probing typical star-forming galaxies with morphologies in the range Sa–Sb
(corresponding roughly to T types between 1 and 6). As noted also by Kennicutt
et al. (1989) the scatter between the $\alpha$ of spirals of the same type in
their sample is comparable to the magnitude of the trend they report to exist
between $\alpha$ and morphology. The lack of a correlation between $\alpha$
and morphology in our sample is thus fully compatible with the results in the
literature. The lack of a correlation with the gas-phase metallicity and the
total stellar mass of our galaxies, despite the differences extracting these
quantities with respect to previous studies, can also be explained by sample
selection effects. Our sample mainly comprises galaxies with high
metallicities (i.e. $12+\log(\mathrm{O/H})>8.3$) and stellar masses (i.e.
$\log M_{*}\,[M_{\odot}]>9$). Cook et al. (2016) found no trend between
$\alpha$ and gas-phase metallicity in the range
$7.2<12+\log(\mathrm{O/H})<9.2$ and only a moderate trend between $\alpha$ and
$M_{*}$ in the range $6<\log M_{*}\,[M_{\odot}]<11$. The parameter space
occupied by our galaxies is compatible with what has been shown by Cook et al.
(2016) and we do not expect to find strong correlations.

On the other hand, SFR-related quantities show a better correlations with
$\alpha$ (Fig. 4d-e-f). SFR, sSFR, and $\Sigma_{\mathrm{SFR}}$ have a
moderate, strong, and very strong correlation with $\alpha$, respectively. The
presence of a tighter correlation with $\Sigma_{\mathrm{SFR}}$ is in line with
the findings of Cook et al. (2016), whose best-fitting relation is also shown
in Fig. 4f. The negative Spearman correlation coefficient indicates that the
quantities anti-correlate with $\alpha$, meaning that the LF has a shallower
slope (i.e. higher relative number of bright H ii regions) for galaxies with
higher SFR, sSFR, and $\Sigma_{\mathrm{SFR}}$.

Looking at where our galaxies lie in the plane of the three classical star
formation scaling relations (Fig. 4g-h-i), we find that $\alpha$ anti-
correlates strongly with $\Delta$ rSFMS and $\Delta$ rMGMS, and moderately
with $\Delta$ rKS. The weaker correlation in the latter case can be somehow
explained by the fact that, between the three, the rKS relation is the one
that shows the lowest galaxy-to-galaxy scatter in our sample (Pessa et al.,
2021). While the trend with $\Delta$ rSFMS confirms the trend we find with
sSFR, $\Delta$ rMGMS adds information on the molecular gas fraction and
indicates that galaxies with a lower fraction of molecular gas tend to have
steeper LFs.

In this paper we choose to not perform any fit for the strongest correlations
shown in Fig. 4 because there is no strong evidence for a particular
functional form from the data (e.g. considering the scatter with respect to
the errors) nor is there one expected from theoretical studies. We note that
future work could explore local correlations (e.g. using H ii region samples
from different areas of the same galaxy) which would give more insights about
this matter but this is behind the scope of the current paper.

Table 3: Global properties of the PHANGS–MUSE galaxies

Galaxy T type Slope $\alpha$ $12+\log(\mathrm{O/H})_{\langle
r_{\mathrm{gal}}\rangle}$ $\log M_{*}$ $\log\mathrm{SFR}$ $\log\mathrm{sSFR}$
$\log\Sigma_{\mathrm{SFR}}$ $\Delta\mathrm{rSFMS}$ $\Delta\mathrm{rKS}$
$\Delta\mathrm{rMGMS}$ $M_{\odot}$ $M_{\odot}\leavevmode\nobreak\
\mathrm{yr}^{-1}$ $\mathrm{yr}^{-1}$ $M_{\odot}\leavevmode\nobreak\
\mathrm{yr}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-2}$ NGC5068 6.0
1.58$\pm$0.09 8.31 9.33 0.14 -10.17 -2.20 0.10 -0.31 0.03 IC5332 6.8
1.82$\pm$0.13 8.37 9.49 0.05 -10.76 -2.82 -0.49 – – NGC1087 5.2 1.71$\pm$0.10
8.40 9.96 0.70 -10.11 -2.25 0.36 0.09 0.39 NGC1385 5.9 1.52$\pm$0.08 8.42 9.94
0.98 -9.95 -2.10 0.52 0.24 0.45 NGC7496 3.2 1.69$\pm$0.10 8.49 10.01 0.29
-10.52 -2.49 -0.09 -0.03 0.02 NGC2835 5.0 1.76$\pm$0.11 8.38 9.87 0.37 -10.28
-2.36 0.11 -0.00 -0.22 NGC0628 5.2 1.71$\pm$0.10 8.48 10.15 0.37 -10.50 -2.36
– – – NGC3351 3.1 1.98$\pm$0.16 8.59 10.41 0.11 -11.26 -2.81 -0.86 -0.25 -0.45
NGC4254 5.2 1.61$\pm$0.08 8.53 10.37 1.29 -10.19 -2.11 0.38 0.06 0.43 NGC4303
4.0 1.56$\pm$0.09 8.56 10.71 2.02 -10.36 -2.02 0.16 0.18 0.11 NGC4535 5.0
1.63$\pm$0.09 8.54 10.44 0.42 -10.81 -2.47 -0.21 -0.18 0.04 NGC1300 4.0
1.77$\pm$0.12 8.52 10.61 0.40 -10.96 -2.94 -0.45 -0.06 -0.42 NGC1512 1.2
2.04$\pm$0.17 8.55 10.66 0.23 -11.21 -3.05 -0.73 -0.17 -0.58 NGC1672 3.3
1.67$\pm$0.10 8.52 10.73 1.23 -10.56 -2.30 0.04 0.02 0.11 NGC4321 4.0
1.95$\pm$0.15 8.56 10.62 0.80 -10.64 -2.42 -0.01 -0.13 0.15 NGC1566 4.0
1.57$\pm$0.09 8.55 10.73 1.46 -10.51 -2.15 -0.04 0.03 0.06 NGC3627 3.1
1.57$\pm$0.09 8.55 10.69 0.84 -10.71 -2.00 0.06 -0.07 0.18 NGC1433 1.5
1.96$\pm$0.16 8.54 10.85 0.38 -11.20 -3.04 -0.85 -0.23 -0.67 NGC1365 3.2
1.74$\pm$0.11 8.52 10.99 1.29 -10.77 -2.49 0.03 -0.02 0.07

666The table reports the galaxy name (col 1), morphological T type taken from
HyperLEDA (Makarov et al., 2014) (col 2), the LF slope (col 3), the
metallicity at the median H ii regions galactocentric radius (col 4), the
stellar mass (col 5), the star formation rate (col 6), the specific star
formation rate (col 7), the star formation rate density (col 8), the offset
from the resolved star formation main sequence (col 9), the offset from the
resolved Kennicutt–Schmidt relation (col 10), and the offset from resolved the
molecular gas main sequence (col 11). All the properties reported in col 3-11
have been extracted using the PHANGS-MUSE data.

###  4.3 LF variations within galaxies

After having investigated LF slope variations across our sample, in this
section, we focus our attention on possible variations within a given galaxy
by looking at the LF of H ii region sub-samples. In the following, we use
galactic environment (spiral arm and inter-arm areas), galactocentric radius,
and ionization parameter to draw, for each galaxy, H ii region sub-samples and
fit their LFs while keeping $L_{\mathrm{min}}$ fixed to the value found in
Sec. 4.1 for each parent sample. All the details of the LFs fits performed in
this section are reported in Table 7, while the actual fits are shown in
Appendix B.

####  4.3.1 Spiral arm vs. inter-arm areas

Differences in the LF slope between spiral arm and inter-arm areas have been
reported in the past (Rand, 1992; Banfi et al., 1993; Thilker et al., 2000),
even though not always significant, with steeper LFs in inter-arm areas.
However, this does not appear to be a universal property of spiral galaxies,
as numerous cases are known where such variations have not been found (Knapen,
1998; Rozas et al., 1996; Azimlu et al., 2011; Gutiérrez et al., 2011). This
has to do with the fact that a simple spatial classification of H ii regions
as spiral or inter-arm regions does not guarantee to probe the conditions
under which star formation happens in the two environments. Part of the H ii
region population in the inter-arm areas, for example, can consist of aged H
ii regions previously formed within spiral arms and/or as part of a recent
star formation burst, depending on the star formation history (SFH) of a
galaxy.

Figure 5 summarizes the result of fitting the LFs of the H ii regions
belonging to spiral arm and inter-arm areas (see Sec. 2.3 for our definition
of the environments and Fig. 31 for the actual LFs fits). Six of our 19
galaxies do not show evident spiral arms, for these galaxies we report the LF
slope for the H ii regions across the entire galaxy disk but excluding bar
regions (i.e. this is the reason for slight variations with respect to the
slopes reported in Table 5).

We find that the LFs in spiral arm areas are either comparable to or shallower
than in inter-arm areas. We find that six out of the 13 galaxies with spiral
arms, preferentially at higher stellar masses, show variations between spiral
arm and inter-arm LF $\alpha$ that go beyond the uncertainties. By comparing
the LF shapes we see that the fraction of bright to faint H ii regions tends
to be higher in spiral arms and lower in inter-arms areas (especially for
NGC1566 and NGC1300), however, in general, the two LFs span quite similar
ranges of luminosities and have similar turnover points. The interpretation of
these results is further discussed in Sec. 6.2.

Figure 5: Slope of the LF for H ii regions located in spiral arms (red points)
and inter-arm (blue points) areas. Galaxies are ordered by increasing stellar
mass from left to right and their names are indicated along the abscissa.
Galaxies with no evident spiral arms are marked by black points and their LF
slope refers to the disk area (excluding bars).

####  4.3.2 Inner vs. outer disk

There exists evidence that the conditions under which SF happens in the outer
disk of spiral galaxies are different from the inner disk, for example, as
suggested by an observed steepening of the GMC mass function (e.g. Rosolowsky
et al., 2007; Colombo et al., 2014; Rice et al., 2016; Faesi et al., 2018;
Schruba et al., 2019). In line with this, radial variations of the LF are
expected and have been reported in the literature suggesting that the LF is
steeper in the outer disk (i.e. extending beyond $R_{25}$; see e.g. Lelièvre &
Roy, 2000) where H ii regions tend to be smaller and fainter (see, e.g.
Ferguson et al., 1998; Helmboldt et al., 2005).

To assess the presence of radial trends (which may also be connected with
parameters varying radially such as the gas-phase metallicity), we study the
LF at small and large galactocentric radii. To mitigate effects due to
different sub-sample sizes, for each galaxy we select the median H ii region
galactocentric radius as a demarcation between the inner and the outer star-
forming disk and extract two H ii region sub-samples of equal sizes. It should
be noted that the area covered by our MUSE observations is largely contained
within ${\sim}1$$R_{25}$ of our targets, as can be seen looking at the
$r_{\mathrm{max}}$ values reported in Table 3. What we define as the outer
star-forming disk is still covering the main (i.e. molecular gas dominated)
star-forming disk of our spiral galaxies and should not be confused with what
is conventionally called outer (i.e. atomic gas dominated) disk in the
literature (typically extending far beyond $R_{25}$).

Fig. 6 summarizes the results of our fits of the LF slope in the inner and
outer star-forming disk (shown in Fig. 32). The median H ii region
galactocentric radii, used to draw the H ii region sub-samples, varies between
${\sim}0.2$$R_{25}$ and $0.5$$R_{25}$ across our galaxies. Most of our
galaxies do not show significant variations between the LF slope of H ii
regions located at small and large galactocentric radii. Our results do not
change if we adopt a fixed value of $0.3$$R_{25}$ to draw H ii region sub-
samples (not shown here). Only two out of 19 galaxies (i.e. NGC1385 and
NGC1566) show a significant variation, indicating a steeper LF in the outer
disk. Our results indicate that within the star-forming disk of spiral
galaxies there are no significant radial trends for the LF slope. This
suggests that the radial metallicity trends in our galaxies (no more than
$0.2$dex variations; see Fig. 30) do not affect the LF slope.

Figure 6: Slope of the LF for H ii regions located in the inner (red points)
and outer (blue points) star-forming disk. Galaxies are ordered by increasing
value of their H ii regions median galactocentric radius from left to right
and their names are indicated along the abscissa.

####  4.3.3 High vs. low ionization parameter

Lastly, we investigate variations of the H ii region LF slope connected with
the ionization parameter $q$ of the H ii regions. As for the radial bins, for
each galaxy we use the median H ii region $q$ to draw two sub-samples of H ii
regions of equal size. The median values that we find are in the range $\log
q=6.4{-}7$. We visually inspected the spatial location of the two sub-samples
of H ii regions for each galaxy and they both appear similarly distributed
across the galaxy disks.

Fig. 7 summarizes the results of the LF fit for the population of H ii regions
with high and low ionization parameter (shown in Fig. 33). It is evident that
for 15 out of 19 galaxies the population of H ii regions with high $q$ has,
within the uncertainties, a shallower LF compared to the H ii regions with low
$q$. Also, looking at the shapes of the LFs we see that for a number of
targets (e.g. IC5332, NGC0628, NGC2835, NGC3627, and NGC5068) the H ii regions
with high/low $q$ constitute the bulk of the bright/faint H ii region
populations. The interpretation of these results is further discussed in Sec.
6.3.

Figure 7: Slope of the LF for H ii regions with high (blue) and low (red)
ionization parameter. Galaxies are ordered by increasing stellar mass from
left to right and their names are indicated along the abscissa.  Table 4: LF
properties of H ii regions sub-samples in the PHANGS–MUSE galaxies.

Galaxy | $\alpha_{a}$ | $n_{a}$ | $\alpha_{i}$ | $n_{i}$ | $\alpha_{r<}$ | $n_{r<}$ | $\alpha_{r>}$ | $n_{r>}$ | $\alpha_{q<}$ | $n_{q<}$ | $\alpha_{q>}$ | $n_{q>}$
---|---|---|---|---|---|---|---|---|---|---|---|---
NGC5068 | – | – | – | – | 1.57$\pm$0.09 | 607 | 1.60$\pm$0.08 | 598 | 2.01$\pm$0.17 | 525 | 1.39$\pm$0.00 | 587
IC5332 | – | – | – | – | 1.87$\pm$0.14 | 199 | 1.76$\pm$0.11 | 158 | 2.39$\pm$0.25 | 134 | 1.58$\pm$0.09 | 189
NGC1087 | – | – | – | – | 1.66$\pm$0.09 | 303 | 1.84$\pm$0.13 | 135 | 2.31$\pm$0.23 | 170 | 1.55$\pm$0.08 | 268
NGC1385 | 1.59$\pm$0.08 | 197 | 1.49$\pm$0.08 | 359 | 1.46$\pm$0.08 | 389 | 1.80$\pm$0.12 | 167 | 1.85$\pm$0.13 | 224 | 1.42$\pm$0.00 | 331
NGC7496 | – | – | – | – | 1.73$\pm$0.11 | 174 | 1.65$\pm$0.09 | 127 | 2.14$\pm$0.20 | 136 | 1.52$\pm$0.08 | 164
NGC2835 | 1.65$\pm$0.09 | 186 | 1.87$\pm$0.14 | 229 | 1.70$\pm$0.10 | 238 | 1.83$\pm$0.13 | 194 | 2.86$\pm$0.31 | 119 | 1.62$\pm$0.08 | 311
NGC0628 | 1.66$\pm$0.09 | 883 | 1.85$\pm$0.13 | 389 | 1.76$\pm$0.11 | 572 | 1.66$\pm$0.09 | 700 | 2.17$\pm$0.20 | 510 | 1.48$\pm$0.09 | 633
NGC3351 | – | – | – | – | 2.03$\pm$0.17 | 204 | 1.93$\pm$0.15 | 165 | 2.07$\pm$0.18 | 158 | 1.74$\pm$0.11 | 152
NGC4254 | 1.58$\pm$0.09 | 792 | 1.65$\pm$0.09 | 637 | 1.57$\pm$0.09 | 921 | 1.72$\pm$0.10 | 508 | 1.74$\pm$0.11 | 583 | 1.55$\pm$0.08 | 846
NGC4303 | 1.49$\pm$0.08 | 864 | 1.68$\pm$0.10 | 668 | 1.52$\pm$0.08 | 937 | 1.65$\pm$0.09 | 630 | 1.83$\pm$0.13 | 731 | 1.44$\pm$0.00 | 823
NGC4535 | 1.54$\pm$0.08 | 341 | 1.71$\pm$0.10 | 545 | 1.64$\pm$0.09 | 455 | 1.62$\pm$0.09 | 479 | 1.70$\pm$0.10 | 459 | 1.50$\pm$0.08 | 393
NGC1300 | 1.67$\pm$0.10 | 397 | 2.13$\pm$0.19 | 184 | 1.79$\pm$0.12 | 344 | 1.75$\pm$0.11 | 289 | 1.87$\pm$0.14 | 337 | 1.67$\pm$0.10 | 289
NGC1512 | 1.96$\pm$0.16 | 125 | 2.18$\pm$0.20 | 77 | 2.18$\pm$0.21 | 108 | 1.92$\pm$0.15 | 101 | 2.83$\pm$0.31 | 78 | 1.81$\pm$0.13 | 128
NGC1672 | 1.57$\pm$0.09 | 145 | 1.79$\pm$0.12 | 265 | 1.67$\pm$0.10 | 307 | 1.67$\pm$0.10 | 263 | 1.97$\pm$0.16 | 228 | 1.55$\pm$0.08 | 337
NGC4321 | 1.84$\pm$0.13 | 217 | 2.15$\pm$0.20 | 134 | 1.91$\pm$0.14 | 233 | 2.02$\pm$0.17 | 150 | 2.00$\pm$0.17 | 213 | 1.89$\pm$0.14 | 168
NGC1566 | 1.47$\pm$0.09 | 586 | 1.86$\pm$0.14 | 407 | 1.52$\pm$0.08 | 644 | 1.71$\pm$0.10 | 383 | 1.70$\pm$0.10 | 506 | 1.48$\pm$0.09 | 513
NGC3627 | 1.48$\pm$0.09 | 320 | 1.70$\pm$0.10 | 282 | 1.51$\pm$0.08 | 350 | 1.65$\pm$0.09 | 307 | 1.80$\pm$0.12 | 259 | 1.45$\pm$0.08 | 371
NGC1433 | – | – | – | – | 2.01$\pm$0.17 | 203 | 1.92$\pm$0.15 | 204 | 2.40$\pm$0.25 | 182 | 1.74$\pm$0.11 | 217
NGC1365 | 1.60$\pm$0.08 | 111 | 1.89$\pm$0.14 | 182 | 1.77$\pm$0.12 | 146 | 1.71$\pm$0.10 | 200 | 2.17$\pm$0.20 | 119 | 1.61$\pm$0.08 | 226

777The table reports the galaxy name (col 1) and the LF slope and the number
of H ii regions involved in the LF fit for the H ii regions sub-samples in
spiral arms areas (col 2–3), in inter-arms areas (col 4–5), in the inner star-
forming disk (col 6–7), in the outer star-forming disk (col 8–9), with low
ionization parameter (col 10–11), and with high ionization parameter (col
12–13).

##  5 Completeness, blending, and selection effects

Figure 8: Left panel. LF slope plotted as a function of the
$\mathrm{FWHM}_{\mathrm{PSF}}\times\sqrt{|\Sigma_{\mathrm{SFR}}|}$, used as a
proxy to quantify the effect of blending. Right panel. LF $L_{\mathrm{min}}$
plotted as a function of the observations spatial resolution. For both panels
we report the Spearman correlation coefficient $\rho$ and $p$-value of the
plotted quantities.

Observational studies of LFs suffer from two main limitations. On the one
hand, faint H ii regions are more difficult to detect against the diffuse
H$\alpha$ background. On the other hand, the brighter H ii regions, often
located in crowded areas such as spiral arms, are more subject to blending
effects. This can artificially lower the number of detected faint regions and
increase the luminosity of the brightest regions, respectively, with
repercussions on the slope of the LF.

The first limitation is related to the sensitivity and the completeness limit
of the observations. We estimate the completeness limit of our LFs (shown in
Fig. 2) in an empirical way by looking at the distribution of the H$\alpha$
flux outside the ionized nebula footprints in the H$\alpha$ emission line
maps, which we refer to as the H$\alpha$ diffuse emission for simplicity.
Given the PSF of the MUSE observations for a given galaxy, we estimate the
completeness limit as the luminosity of a mock circular H ii region with the
size of the PSF and a uniform $\Sigma_{\mathrm{H\alpha}}$ equal to the 90th
percentile (i.e. slightly larger than $1\sigma$ of the distribution) of the
H$\alpha$ diffuse emission distribution. In this way, we take into account the
bright tail of the diffuse H$\alpha$ background against which H ii regions can
remain undetected. It is worth noting that tests injecting mock H ii regions
are a good way forward to assess more robustly the completeness limit in the
future.

As shown in Fig. 2, the completeness limit estimated in this way lies below
the best-fitting $L_{\mathrm{min}}$ (i.e. the LF turnover luminosity)
suggesting that the presence of a turnover point is not a completeness problem
but rather an intrinsic feature of the LFs of our galaxies (see Youngblood &
Hunter, 1999, and references therein for a discussion). Models have shown that
a turnover point is to be expected when considering a constant SFR over time,
while different SFH (e.g. a star formation burst in the past) can affect the
shape of the LF in a more complex way (see e.g. M31’s double-peaked LF is
attributed to a recent starburst; Azimlu et al., 2011). We note that, by
taking into account the bright tail of the diffuse H$\alpha$ emission, we
estimate a quite conservative completeness limit especially for the targets
with stronger diffuse background (i.e. the galaxies for which the completeness
limit is located close to the peak of the LF in Fig. 2). However, it should be
noted that higher sensitivity and resolution are expected to move the turnover
point to lower luminosities and, as shown in the right panel of Fig. 8,
$L_{\mathrm{min}}$ is indeed driven by resolution of the observations in our
survey. This constitutes a validation of our fitting approach and ensures
that, at a given spatial resolution, the LFs $\alpha$ is extracted from a
well-sampled luminosity regime.

On the other hand, the issue of blending is related to the spatial resolution
of the observations and the number density (or filling factor) of star-forming
regions. Intuitively, observations with coarser spatial resolution or a higher
filling factor will cause nearby H ii regions to blend and be detected as a
single H ii region, leading to a flattening of the LF. To assess the effect of
blending to a first order, for each target, we multiply the
$\mathrm{FWHM}_{\mathrm{PSF}}$, that takes into account the spatial resolution
of the data, by $\sqrt{|\Sigma_{\mathrm{SFR}}|}$, that gives us a handle on
the mean separation between H ii regions. We use this parameter as a proxy for
blending and check its relation with $\alpha$; as shown in the left panel of
Fig. 8, we find a moderate correlation. Different studies (Kennicutt et al.,
1989; Bastian et al., 2007; Cook et al., 2016) showed that spatial resolution
does not significantly change the slope of the LF at resolutions finer than
few hundred parsecs like in the case of our data (i.e. the coarser physical
resolution in our sample is $\sim 110$ pc for NGC1512 and NGC1365). This
indicates that, although present, the effect of blending is not the main
driver of the correlations we report in Sec. 4.2. It is worth noting that our
checks for the blending effect are either indirect or made a posteriori. In
the future, progressively degrading the MUSE data cubes at a common resolution
and analyzing how the LF parameters vary will give us a more comprehensive
understanding of this matter.

In Appendix C, we also discuss how the selection criteria, we use to select
our H ii regions (see Sec. 3.2), have no significant effect on the
measurements of the LF slope presented in Sec. 4.1.

##  6 Discussion: what sets the LF slope?

###  6.1 The SFR surface density

In Sec. 4.2, we find that in general $\alpha$ correlates better with the
global star formation properties of our galaxies and especially with
$\Sigma_{\mathrm{SFR}}$. We do not find any clear trends within the T-type,
$\mathrm{[O/H]}$, and $M_{*}$ ranges probed by our sample. This indicates that
the properties of star-forming regions are more closely connected to the
global star formation properties in our sample.

One caveat to keep in mind is that the correlation between the LF and the SFR-
related properties may be driven by the so-called size-of-sample effect
(Larsen, 2002; Weidner et al., 2004; Bastian, 2008; Cook et al., 2012;
Whitmore et al., 2014). This is a statistical effect (i.e. stochastic
sampling) due to the fact that galaxies with higher SFRs tend to have more H
ii regions and, in turn, the chances of observing a higher number of bright H
ii regions increase in these galaxies. This is expected to drive a tight
correlation between the total number of detected H ii regions (and
subsequently the SFR) and the luminosity of the brightest star-forming region.
Intuitively, the resulting effect is a flattening of the LF for galaxies with
higher SFRs. To check for this effect in our sample, in Fig. 9, we plot the
number of H ii regions (or similarly the number of ionized nebulae) against
the maximum H ii region luminosity $L(\mathrm{H}\alpha)_{\mathrm{max}}$ and
the global SFR. We find only a weak and a moderate trend, respectively. This
indicates that, thanks to the resolution of our observations and the nature of
our sample, the flatter slopes we observe for the more star-forming galaxies
are not simply due to a statistical effect. In addition, Cook et al. (2016)
tested the size-of-sample effect using simulated LFs with $\alpha=-2$ and
showed that the stochastic scatter of the LF slope starts to increase
symmetrically when the number of detected H ii regions drops below $100$. This
implies that low number statistics increases the scatter in the measured LF
slope but does not drive systematic changes of $\alpha$ in a specific
direction. Moreover, our LFs are built with H ii region samples containing at
least ${\sim}500$ objects, and the LF fits presented in Sec. 4.1 rely on at
least $200$ objects. For this reason, we expect the stochastic scatter of
$\alpha$ to be lower in our sample compared to previous studies that often
relied on smaller H ii region samples.

Figure 9: Maximum luminosity among the detected H ii regions (left panel) and
the total SFR (right panel) as a function of the total number of detected H ii
regions for the galaxies of our sample. For each panel we report the Spearman
correlation coefficient $\rho$ between the plotted quantities.

In agreement with Cook et al. (2016), who studied the LF in a sample of $258$
nearby galaxies using GALEX FUV data, we find that the LF slope correlates
best with $\Sigma_{\mathrm{SFR}}$. In Fig. 4f, we show the best-fitting
relation found by Cook et al. (2016) for a qualitative comparison. We note
that our points mostly lie above this relation. A direct comparison is,
however, not straightforward because of the different type of observations. In
fact, FUV and H$\alpha$ LFs probe different ages of H ii regions,
${\lesssim}100$ Myr and ${\lesssim}10$ Myr, respectively (see e.g. Haydon et
al., 2020). Oey & Clarke (1998b) predicted that older H ii region populations
(e.g. derived from FUV) should have steeper slopes than younger populations
(e.g. derived from H$\alpha$). Overall, the average slope of the nebular LF
($-2\pm 0.5$; Kennicutt et al., 1989; Elmegreen & Salzer, 1999), including our
study ($-1.73\pm 0.15$), is compatible with what has been found in the FUV
($-1.76\pm 0.3$; Cook et al., 2016).

The observed trend indicates that a shallower LF and, thus, the formation of
higher relative number of massive, luminous H ii regions is connected with a
higher $\Sigma_{\mathrm{SFR}}$. We know that the star formation process is
fueled by cold gas. We thus expect the SFR-related properties we measure to be
connected with the way how cold gas reservoir fuels star formation, which in
turn is connected with the sSFR, the $t_{\mathrm{dep}}$, and the molecular gas
fraction.

As discussed in Sec. 4.2, the offsets from the resolved star formation scaling
relations can be used as a proxy for these parameters. We find that the
stronger correlation is between the LF slope and $\Delta$ rSFMS (i.e. a
measurement of the sSFR) indicating that, when the mass fraction between young
and old stars is higher, galaxies tend to form a higher relative number of
luminous H ii regions, in agreement with what we find when estimating the
global sSFR directly from the MUSE data (Fig. 4e). A slightly weaker but still
significant trend is found between the LF slope and $\Delta$ rMGMS (i.e. a
measurement of the molecular gas fraction) suggesting that, when in a galaxy
the fraction of cold gas that will eventually be converted into stars is
higher, the relative number of bright and massive H ii regions produced by
ongoing star formation is also higher.

Kruijssen (2012) formulated a theoretical framework (invoked also by Cook et
al., 2016 to explain their results), in which higher SFR and cold gas
densities can result in a higher pc-scale star formation efficiency that, in
turn, enhances the cluster formation efficiency, defined as the ratio between
the cluster formation rate (CFR) and the SFR
($\Gamma\equiv\mathrm{CFR}/\mathrm{SFR}$). This means that, at higher SFR and
gas densities, more stars form in (gravitationally bound) star clusters.
Observational evidence has indeed been found that $\Gamma$ correlates well
with the SFR surface density (Goddard et al., 2010; Adamo et al., 2011; Cook
et al., 2012; Adamo et al., 2015; Johnson et al., 2016; Adamo et al., 2020b,
a, see e.g.). The trend that we find in Fig. 4e between $\alpha$ and
$\Sigma_{\mathrm{SFR}}$ goes in the direction predicted by this scenario. It
is worth noting that, with the available PHANGS–HST data (Lee et al., 2021),
it will be possible to directly measure $\Gamma$ and compare it to our
measured H ii region LF slopes.

Quantitatively, the slope of the H ii region LF might be expected to be
similar to the slope of the GMC mass function (historically estimated to lie
between $1.6$ and $2.0$, e.g. Kennicutt & Evans, 2012b, although more recent
works find a wider range of slopes, see e.g. Colombo et al. 2014; Rosolowsky
et al. 2021). If star formation within GMCs is self-similar (Efremov &
Elmegreen, 1998) and the star formation efficiency is independent of the GMC
mass at the observable, high-mass end of the mass range, then the slope of the
H ii region LF is also expected to be $\alpha=1.6{-}2.0$. Steeper slopes may
then be the result of incomplete sampling of the stellar initial mass function
(IMF), which occurs if the stellar populations born within the parent GMCs
have a low mass (${\la}10^{4}\leavevmode\nobreak\ M_{\odot}$; Krumholz et al.
2015). Simple analytical models, in which the maximum GMC and stellar
population mass scales are regulated by gravitational instability and stellar
feedback, predict that the maximum mass increases approximately with the total
gas surface density as $\Sigma_{\textrm{gas}}^{1.4}$, and therefore
approximately linearly with $\Sigma_{\mathrm{SFR}}$ (Kruijssen & Longmore,
2014; Reina-Campos & Kruijssen, 2017). This linear relation is indeed observed
by Johnson et al. (2017), who found that a maximum mass of
$10^{4}\leavevmode\nobreak\ \textrm{M}_{\odot}$ is reached at
$\log{(\Sigma_{\textrm{SFR}}\ [\textrm{M}_{\odot}\leavevmode\nobreak\
\textrm{yr}^{-1}\leavevmode\nobreak\ \textrm{kpc}^{-2}])}\approx-2.6$. Above
this SFR surface density, we expect the H ii region LF to be close to the
slope of the GMC mass function (which in environments of high gas and SFR
surface densities is about $1.6$, e.g. Colombo et al. 2014; Hughes et al.
2016). Toward lower SFR surface densities, we expect the H ii region LF to
steepen. The same behavior, in which $\alpha$ declines with
$\Sigma_{\mathrm{SFR}}$ for $\log{(\Sigma_{\textrm{SFR}}\
[\textrm{M}_{\odot}\leavevmode\nobreak\ \textrm{yr}^{-1}\leavevmode\nobreak\
\textrm{kpc}^{-2}])}\la-2.6$ and stays approximately constant (or declines
less steeply) for $\log{(\Sigma_{\textrm{SFR}}\
[\textrm{M}_{\odot}\leavevmode\nobreak\ \textrm{yr}^{-1}\leavevmode\nobreak\
\textrm{kpc}^{-2}])}\ga-2.6$ was also found by Cook et al. (2016) and is
illustrated in (Fig. 4f). Here, as a reference to illustrate this flattening,
we marked with horizontal lines the $\alpha$ corresponding to the Cook et al.
(2016) best-fitting line at $\Sigma_{\mathrm{SFR}}\approx 2.6$ and the mean
$\alpha$ of our samples at $\Sigma_{\mathrm{SFR}}>2.6$. Both values are
remarkably close to $\alpha=1.6$, the slope of the GMC mass function in the
high gas and $\Sigma_{\mathrm{SFR}}$regime. Based on the discussion in this
section, we expect this behavior to arise from a combination of increased
stellar clustering at birth toward higher gas densities (and therefore
pressures; Kruijssen 2012), as well as of decreased IMF sampling toward low
SFR surface densities. Higher-resolution observations may be necessary to
assess how blending of adjacent H ii regions may affect these trends.

###  6.2 The strength of spiral arms

In Sec. 4.3.1, we showed that a flattening of the LF slope in spiral arm areas
compared to inter-arm areas is evident for only about half of the galaxies
with spiral arms in our sample (i.e. six out of 13 galaxies, namely NGC1300,
NGC1365, NGC1566, NGC1672, NGC3627, and NGC4303). The Monte Carlo simulations
by Oey & Clarke (1998a) suggested that such variations can be expected when
the H ii regions populating the spiral arm areas trace a current burst of
coeval star formation while the inter-arm population is an aged version
thereof (i.e. aging effect). However, the same authors noted that, at any
given location in a galaxy disk, the timescale between the passage of spiral
density waves is of the order of $40$ Myr (see e.g. Rand, 1993), much longer
than the observed lifetime of H ii regions (in the interval $5{-}10$ Myr;
Chevance et al., 2020b). This implies that H ii regions detected in inter-arm
areas are, potentially, tracing both a population of aging H ii regions and
recent star formation that is genuinely happening in inter-arm areas.

Large-scale dynamical processes shaping the distribution of cold gas that
ultimately fuels the star formation can also influence the variation (or not)
in the LF slopes between spiral arm and inter-arm areas. We might expect this
to depend on the nature of the dynamical perturbation present in the disk
(material spiral arms vs. density waves) and its strength. Rand (1992), for
example, found a significant difference between the LF slopes of spiral arm
and inter-arm H ii regions in the grand-design interacting galaxy M51 and
attributed it to a shallower GMCs MF in the spiral arms. In NGC 6814, on the
other hand, Knapen et al. (1993) argued that the spiral arms are not strong
enough to build up the large cloud masses needed to produce the number of
giant H ii regions that make the spiral arm LF shallower than the inter-arm
LF.

Strong spiral arms are not only capable of building high gas densities and
growing high mass molecular clouds, they can also impact the likelihood of in
situ inter-arm star formation, thus influencing the LFs of H ii regions.
Indeed, the presence of spiral arms imposes a strong organization on the gas
distribution, which acts to concentrate star formation into the spiral arms.
Comparing a high resolution simulation of an M51-like galaxy with a comparable
isolated galaxy, Tress et al. (2020, 2021) concluded that strong spiral arms
gather molecular gas and clouds without dramatically affecting their
properties, yielding similar overall star formation rates.

In our case, we might thus expect that variations in cloud mass functions and
H ii region LF slopes will be weak in systems with more uniformly distributed
H ii regions across the disk (weaker spirals) and stronger in systems with
strong spiral arms. In line with this, we find that the galaxies showing a
significant environmental change in $\alpha$ are those in which the
distribution of the H ii regions traces quite clearly the spiral arms and
fewer H ii regions are detected in between spiral arms (i.e. less ongoing star
formation in the inter-arm areas). The galaxies with more uniformly
distributed H ii regions across the disk (e.g. NGC0628, NGC4254, and NGC4321),
on the other hand, exhibit comparable $\alpha$ in the two environments.

We quantify this further using the azimuthal contrast in CO brightness across
our targets, which serves as a proxy for spiral arm strength (Meidt et al.,
2021). For our targets, we measure the CO contrasts in a set of 150-pc wide
radial bins, calculating the ratio between the brightness of the $84$th
percentile of the CO distribution in a given ring and a reference level
defined to capture the level of low brightness inter-arm emission in Meidt et
al. (2021). We then take the average of the CO contrasts measured across the
full area of interest.

As shown in Fig. 10, for the six galaxies showing environmental changes of the
LF slope, the magnitude of such variations appears to correlate with the
spiral arm contrast. The remaining galaxies with evidently weaker spiral arms,
more uniformly distributed H ii regions, and no significant variation in LF
slope from arm to inter-arm, do not exhibit a similar trend.

Overall, these findings suggest that environmental variations in $\alpha$ are
sensitive to how spiral arms shape the distribution of the more massive
molecular clouds that serve as the birthplace of the brightest H ii regions. A
study of environmental variations in GMC mass functions in our sample (e.g.
Rosolowsky et al., 2021, Hughes et al. in preparation) will be an important
tool for recognizing how much of the change in LF slope from arm to inter-arm
is due to the evolution in the distribution of progenitor clouds and how much
is due to, i.e. aging, as considered in the next section. The morphology of
inter-arm gas (in the form of spurs/feathers and non-spiral clump features)
may also yield important clues on the nature of inter-arm star formation and
the change in H ii region LF slope.

Figure 10: Change in the LF slope between inter-arm and spiral arm
environments $\Delta\alpha_{ENV}$ (derived as $\alpha_{i}-\alpha_{a}$ as
reported in Table 7) as a function of the CO(2–1) contrast between the $84$th
percentile and the reference level defined in Meidt et al. (2021) for six
galaxies in our sample.

Given the resolution of our observations, we cannot fully exclude that H ii
region blending, which affects spiral arms to a greater extent, has an effect
on our findings. An additional factor that we are not considering and that can
play a role in this context is the SFH of the galaxies. Feinstein (1997)
predicted that a past starburst can manifest as an additional peak below the
LF turnover luminosity and result in a double-peaked LF as the one observed
for M31 by Azimlu et al. (2011). We would need observations that are at least
one order of magnitude deeper in terms of $L_{\mathrm{H\alpha}}$ to search for
such features and, for this reason, we cannot deduce or fold in information on
the SFH in our study of the LF.

###  6.3 The aging effect and the ionization parameter

Figure 11: Change in LF slope between H ii regions with low and high
ionization parameter (derived as $\alpha_{q<}-\alpha_{q>}$ as reported in
Table 7) as a function of the feedback timescales (upper panel) and the GMC
lifetime (lower panel). For each panel we report the Spearman correlation
coefficient $\rho$ between the plotted quantities.

According to Charlot & Longhetti (2001, eq. (10)), by assuming that H ii
regions are well approximated by Strömgren spheres, the ionization parameter
changes with time as $U(t)\propto(Q(t)n_{\mathrm{H}})^{1/3}$, where $Q$ is the
rate of ionizing photons produced by the central star/cluster and
$n_{\mathrm{H}}$ is the hydrogen density. Considering the typical evolution of
an H ii region, we can expect the ionization parameter to drop as time passes:
as the central ionizing star/cluster evolves its ionizing flux drops while the
gas density decreases due to stellar feedback dispersing the H ii region’s gas
(see e.g. Rahner et al., 2017, 2019; Pellegrini et al., 2020).

The LF slope can, in principle, be used to test this expectation by using the
ionization parameter to draw sub-samples of H ii regions in the way we did in
Sec. 4.3.3. If the ionization parameter $q$ is related to the evolutionary
stage of an H ii region, it would separate less and more evolved H ii regions
and, ultimately, we would expect the aging effect to drive changes in the LF
slope (Oey & Clarke, 1998a). This is the same effect discussed for
environmental variations in the previous section. What we find in Sec. 4.3.3
is indeed a steeper LF for the population of H ii regions with lower $q$ that,
for a number of galaxies, constitutes the bulk of the faint H ii regions. This
corroborates the expectations of Charlot & Longhetti (2001) from an empirical
point of view, indicating that the ionization parameter can be used as a good
proxy for the evolutionary stage of H ii regions.

However, we know that H ii regions are not always well represented by
Strömgren spheres, e.g. they can be asymmetric or even broken shells (see e.g.
Pellegrini et al., 2012). Also, assuming equal physical conditions for an H ii
region, the same ionizing flux (and thus ionization parameter) can be produced
by a single very young star or a cluster of slightly older stars. Another
parameter that is relevant for the interpretation of our results is the
metallicity, and its long-debated degeneracy with age. In fact, an increase in
stellar metallicity is also expected to lower the ionization parameter: as
stellar atmospheres of O stars become cooler (Massey et al., 2005) and the
mechanical energy due to photon scattering increases, dispersing gas more
efficiently (Dopita et al., 2006).

Interestingly, using PHANGS–MUSE data, Kreckel et al. (2019) showed that H ii
regions with higher [S iii]/[S ii]line ratios (and thus higher ionization
parameter) also tend to have higher H$\alpha$-to-FUV flux ratios, indicating
younger star cluster ages. We note that the PHANGS–HST data set will be an
ideal test-bed to actually connect our H ii regions (and their ionization
parameters) to the ages of their ionizing star(s)/cluster and shed a light on
this topic. The [S iii]/[S ii]line ratio we use to compute the ionization
parameter has been found to have only a weak dependence on the gas-phase
metallicity for H ii regions located in spiral galaxies. (Dors & Copetti,
2005; Garnett et al., 1997; Kennicutt & Garnett, 1996; Dors et al., 2011).
However, the topic is still debated as other works presented contrasting
results (see e.g. Bresolin et al., 1999 for disk galaxies and Dopita et al.,
2014 for Luminous Infrared Galaxies) and, above all, in our sample a clear
correlation has been found between [S iii]/[S ii]and the gas-phase metallicity
(see fig. 4 in Kreckel et al., 2019). We confirm this by looking at the
distribution of the H ii region metallicities for the sub-samples with high
and low $q$. However, if the gas-phase metallicity plays a role in setting the
slope of the LF, we would see it by defining H ii region sub-samples based on
their metallicity and looking for variations in $\alpha$. We conducted this
test and do not find such variations. We thus conclude that the variations in
the LF slope we observe in Fig. 33 are mainly driven by aging of the H ii
regions rather than their metal content.

By comparing the spatial distribution of the warm ionized gas (i.e. H ii
regions) and the cold molecular gas (i.e. GMCs), it is now possible to infer
e.g. the lifetime of GMCs $t_{\mathrm{CO}}$ and the timescale over which
stellar feedback acts $t_{\mathrm{fb}}$ (see e.g. Kruijssen et al., 2019;
Chevance et al., 2020a; Kim et al., 2021). Assuming that the aging effect is
the main driver behind the LF slope variations we observe between the H ii
region sub-samples with high and low $q$, we can reasonably expect the
magnitude of such variations $\Delta\alpha_{q}$ to correlate with the
timescales that regulate the star formation process (e.g. the lifetime of GMCs
serving as the stellar nurseries, the timescales of stellar feedback) in a
given galaxy. Following the “uncertainty principle for star formation” method
of Kruijssen et al. (2018) and expanding the results of Chevance et al.
(2020a), Kim et al. (in preparation) derived $t_{\mathrm{CO}}$ and
$t_{\mathrm{fb}}$ for most of the galaxies in our sample. In that work, H ii
regions are traced via the H$\alpha$ emission in the PHANGS–H$\alpha$ data
(Razza et al. in preparation), while GMCs are traced via the CO(2–1) emission
in the PHANGS–ALMA data (Leroy et al., 2021). It should be noted that for
IC5332 there is no significant CO(2–1) emission detected to apply this method
and the $t_{\mathrm{CO}}$ and $t_{\mathrm{fb}}$ measurements for NGC0628 and
NGC3627 are taken from Chevance et al. (2020a). As shown in Fig. 11, we find a
weak-to-moderate/moderate anti-correlation between $\Delta\alpha_{q}$ and both
$t_{\mathrm{fb}}$ and $t_{\mathrm{CO}}$. These trends suggest that when
stellar feedback occurs on shorter timescales and the GMC lifetime is shorter
(i.e. they get dispersed faster due to stellar feedback) the difference
between the LF slope of H ii regions with high and low $q$ increases, as it
would be expected if those populations are intrinsically different in terms of
their age.

##  7 Conclusions

In this paper, we have studied the nebular luminosity functions (LFs; i.e.
built using the H$\alpha$ luminosity of H ii regions) in the star-forming
disks of 19 nearby galaxies making up the PHANGS–MUSE sample. Thanks to the
exquisite spatial resolution (mean $\mathrm{PSF}_{\mathrm{FWHM}}=67$ pc) and
sensitivity of the data, we were able to build a catalog of about $31\,400$
ionized nebulae from which we extracted an unprecedented sample of about
$23\,000$ H ii regions. With MUSE covering a large part of the optical
spectrum we were able to derive global properties (e.g. $M_{*}$, SFR, sSFR)
for our galaxies and characterize our H ii regions in terms of their dust
attenuation (via the H$\alpha$/H$\beta$ Balmer decrement), gas-phase
metallicity $\mathrm{O/H}$, and ionization parameter $q$. The average number
of H ii regions detected per galaxy is about $1200$, marking a significant
improvement with respect to previous spectroscopic studies at a comparable
spatial resolution.

We fit the LFs via maximum likelihood estimation (MLE; Alstott et al., 2014).
This method does not require us to bin the data, is especially suited for
heavily tailed distributions, and overcomes typical pitfalls of more commonly
used methods (e.g. histograms combined with linear regression). We find an LF
slope $\alpha=-1.73\pm 0.15$ in agreement with what has been found in the
literature for Sa–Sc type galaxies (see e.g. Kennicutt et al., 1989; Elmegreen
& Salzer, 1999; Whitmore et al., 2014). Given the sensitivity of the MUSE
observations, our LFs are mainly probing massive star formation happening in
star clusters and associations. We report evidence that for a number of
galaxies their LFs steepen at luminosities above $\log(L_{\mathrm{H\alpha}}\
[\mathrm{erg\leavevmode\nobreak\ s^{-1}}])=38.6$ (i.e. type II LF). However,
this appears as a subtle feature compared to what has been shown in past
studies (e.g. Beckman et al., 2000). In general, a power-law with a single
slope offers a good representation of the LFs and allows us to analyze
uniformly the LFs across our sample.

We find that the LF slope remains largely unchanged for the inner and outer
parts of the star-forming galaxy disks (Sec. 4.3.2). For the galaxies with H
ii regions more uniformly distributed across the star-forming disk, we do not
find any significant difference between the spiral arm and inter-arm LFs
slope. On the other hand, we find environmental (i.e. arm vs. inter-arm) LF
variations for galaxies showing H ii regions preferentially located along
spiral arms. We attribute these variations to the spiral arms increasing the
molecular clouds arm–inter-arm mass contrast and find suggestive evidence that
they are more pronounced in galaxies with spiral arms driving stronger
dynamical perturbations.

In line with the results of Cook et al. (2016), who studied the LF of star-
forming regions as traced by GALEX FUV imaging, we find that galaxies with a
higher $\Sigma_{\mathrm{SFR}}$ have a flatter LF, meaning that their relative
number of bright H ii regions is higher compared to galaxies with lower
$\Sigma_{\mathrm{SFR}}$. This potentially connects to fundamental changes in
the physics regulating star formation in galaxy disks. The trend we observe
between $\alpha$ and $\Sigma_{\mathrm{SFR}}$ is not driven by low-number
statistics and brings further evidence to the prediction that the clustering
of young stars is enhanced at high gas surface densities and star formation
efficiencies (Kruijssen, 2012). We propose that $\alpha$ may increase even
further at $\log{(\Sigma_{\textrm{SFR}}\ [M_{\odot}\leavevmode\nobreak\
\textrm{yr}^{-1}\leavevmode\nobreak\ \textrm{kpc}^{-2}])}\la-2.6$, where
observations and models find that the maximum masses of the new-born stellar
populations should become so low (e.g. Johnson et al., 2017; Reina-Campos &
Kruijssen, 2017) that they are affected by stochastic sampling of the stellar
IMF (e.g. Krumholz et al., 2015).

Finally, we find that, within each of our galaxies, H ii regions with high
ionization parameter have a shallower LF compared to H ii regions with low
ionization parameter. Such variations are compatible with an aging effect
suggesting that H ii regions with high $q$ are the youngest star-forming
regions while the ones with low $q$ are more evolved. The lack of $\alpha$
trends/variations related to the H ii regions’ gas-phase metallicity persuades
us that the main parameter regulating changes in $q$ is age. We also bring
some tentative evidence suggesting that $\alpha$ variations between these H ii
region populations are more evident when stellar feedback acts on shorter
timescales and GMCs lifetimes are shorter.

##  8 Future prospects

We note that this work can be expanded in a number of promising directions in
the near future, especially thanks to the availability of HST and ALMA data
for the entire PHANGS–MUSE sample. The PHANGS–ALMA data will allow us to
characterize GMCs and study their mass function (MF). Comparing GMC MF slopes
with the H ii regions LF slopes reported in this paper can bring new insights
into if/how physical properties of the molecular gas influence the physics of
star formation. On the other side, the PHANGS–HST data can be used to estimate
the ages and metallicities of the star(s)/cluster ionizing the H ii regions
and help disentangling what is the main driver behind changes in the H ii
regions’ ionization parameter. Furthermore, by providing a direct estimate of
the clustered star formation efficiency (defined as the ratio between the
cluster formation rate and the SFR; $\Gamma\equiv\mathrm{CFR}/\mathrm{SFR}$)
it will be possible to compare $\Gamma$ to the H ii region LF slopes measured
in this paper and test whether this is the mechanism behind the flattening of
the LF that we observe with increasing $\Sigma_{\mathrm{SFR}}$.

###### Acknowledgements.

This work was carried out as part of the PHANGS collaboration. FS, ES, RME,
TS, and TGW acknowledge funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 694343). ATB would like to acknowledge funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No.726384/Empire). JMDK and MC
gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through an Emmy Noether Research Group (grant
number KR4801/1-1), as well as from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme via the
ERC Starting Grant MUSTANG (grant agreement number 714907). JMDK, MC, and JK
gratefully acknowledge funding from the DFG through the DFG Sachbeihilfe
(grant number KR4801/2-1). ER acknowledges the support of the Natural Sciences
and Engineering Research Council of Canada (NSERC), funding reference number
RGPIN-2017-03987. KK gratefully acknowledges funding from the German Research
Foundation (DFG) in the form of an Emmy Noether Research Group (grant number
KR4598/2-1, PI Kreckel). EJW acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
138713538 – SFB 881 (“The Milky Way System”, subproject P2). Based on
observations collected at the European Southern Observatory under ESO
programmes 1100.B-0651, 095.C-0473, and 094.C-0623 (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).

## References

  * Adamo et al. (2020a) Adamo, A., Hollyhead, K., Messa, M., et al. 2020a, MNRAS, 499, 3267
  * Adamo et al. (2015) Adamo, A., Kruijssen, J. M. D., Bastian, N., Silva-Villa, E., & Ryon, J. 2015, MNRAS, 452, 246
  * Adamo et al. (2011) Adamo, A., Östlin, G., & Zackrisson, E. 2011, MNRAS, 417, 1904
  * Adamo et al. (2020b) Adamo, A., Zeidler, P., Kruijssen, J. M. D., et al. 2020b, Space Sci. Rev., 216, 69
  * Alstott et al. (2014) Alstott, J., Bullmore, E., & Plenz, D. 2014, PLoS ONE, 9, e85777
  * Anand et al. (2021) Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621
  * Azimlu et al. (2011) Azimlu, M., Marciniak, R., & Barmby, P. 2011, AJ, 142, 139
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 773508
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Banfi et al. (1993) Banfi, M., Rampazzo, R., Chincarini, G., & Henry, R. B. C. 1993, A&A, 280, 373
  * Barnes et al. (2020) Barnes, A. T., Longmore, S. N., Dale, J. E., et al. 2020, MNRAS, 498, 4906
  * Bastian (2008) Bastian, N. 2008, MNRAS, 390, 759
  * Bastian et al. (2007) Bastian, N., Ercolano, B., Gieles, M., et al. 2007, MNRAS, 379, 1302
  * Bastian et al. (2009) Bastian, N., Gieles, M., Ercolano, B., & Gutermuth, R. 2009, MNRAS, 392, 868
  * Beckman et al. (2000) Beckman, J. E., Rozas, M., Zurita, A., Watson, R. A., & Knapen, J. H. 2000, AJ, 119, 2728
  * Belfiore et al. (2021) Belfiore, F., Santoro, F., & Groves, B. 2021, A&A, submitted
  * Bittner et al. (2019) Bittner, A., Falcón-Barroso, J., Nedelchev, B., et al. 2019, A&A, 628, A117
  * Blanc et al. (2009) Blanc, G. A., Heiderman, A., Gebhardt, K., Evans, Neal J., I., & Adams, J. 2009, ApJ, 704, 842
  * Bresolin et al. (1999) Bresolin, F., Kennicutt, Robert C., J., & Garnett, D. R. 1999, ApJ, 510, 104
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cepa & Beckman (1990) Cepa, J. & Beckman, J. E. 1990, A&AS, 83, 211
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot & Longhetti (2001) Charlot, S. & Longhetti, M. 2001, MNRAS, 323, 887
  * Chevance et al. (2020a) Chevance, M., Kruijssen, J. M. D., Hygate, A. P. S., et al. 2020a, MNRAS, 493, 2872
  * Chevance et al. (2020b) Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E., et al. 2020b, Space Sci. Rev., 216, 50
  * Clauset et al. (2009) Clauset, A., Shalizi, C. R., & Newman, M. E. J. 2009, SIAM Review, 51, 661
  * Colombo et al. (2014) Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ, 784, 3
  * Cook et al. (2016) Cook, D. O., Dale, D. A., Lee, J. C., et al. 2016, MNRAS, 462, 3766
  * Cook et al. (2012) Cook, D. O., Seth, A. C., Dale, D. A., et al. 2012, ApJ, 751, 100
  * Diaz et al. (1991) Diaz, A. I., Terlevich, E., Vilchez, J. M., Pagel, B. E. J., & Edmunds, M. G. 1991, MNRAS, 253, 245
  * Dopita et al. (2006) Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJS, 167, 177
  * Dopita et al. (2014) Dopita, M. A., Rich, J., Vogt, F. P. A., et al. 2014, Ap&SS, 350, 741
  * Dors & Copetti (2005) Dors, O. L., J. & Copetti, M. V. F. 2005, A&A, 437, 837
  * Dors et al. (2011) Dors, O. L., J., Krabbe, A., Hägele, G. F., & Pérez-Montero, E. 2011, MNRAS, 415, 3616
  * Efremov & Elmegreen (1998) Efremov, Y. N. & Elmegreen, B. G. 1998, MNRAS, 299, 588
  * Elmegreen et al. (1996) Elmegreen, B. G., Elmegreen, D. M., Salzer, J. J., & Mann, H. 1996, ApJ, 467, 579
  * Elmegreen & Falgarone (1996) Elmegreen, B. G. & Falgarone, E. 1996, ApJ, 471, 816
  * Elmegreen & Salzer (1999) Elmegreen, D. M. & Salzer, J. J. 1999, AJ, 117, 764
  * Emsellem et al. (2021) Emsellem, E., Schinnerer, E., Santoro, F., et al. 2021, arXiv e-prints, arXiv:2110.03708
  * Espinosa-Ponce et al. (2020) Espinosa-Ponce, C., Sánchez, S. F., Morisset, C., et al. 2020, MNRAS, 494, 1622
  * Faesi et al. (2018) Faesi, C. M., Lada, C. J., & Forbrich, J. 2018, ApJ, 857, 19
  * Feinstein (1997) Feinstein, C. 1997, ApJS, 112, 29
  * Ferguson et al. (1998) Ferguson, A. M. N., Wyse, R. F. G., Gallagher, J. S., & Hunter, D. A. 1998, ApJ, 506, L19
  * Garnett et al. (1997) Garnett, D. R., Shields, G. A., Skillman, E. D., Sagan, S. P., & Dufour, R. J. 1997, ApJ, 489, 63
  * Goddard et al. (2010) Goddard, Q. E., Bastian, N., & Kennicutt, R. C. 2010, MNRAS, 405, 857
  * Gutiérrez et al. (2011) Gutiérrez, L., Beckman, J. E., & Buenrostro, V. 2011, AJ, 141, 113
  * Haffner et al. (2009) Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Reviews of Modern Physics, 81, 969
  * Hao et al. (2011) Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124
  * Haydon et al. (2020) Haydon, D. T., Kruijssen, J. M. D., Chevance, M., et al. 2020, MNRAS, 498, 235
  * Helmboldt et al. (2005) Helmboldt, J. F., Walterbos, R. A. M., Bothun, G. D., & O’Neil, K. 2005, ApJ, 630, 824
  * Herrera-Endoqui et al. (2015) Herrera-Endoqui, M., Díaz-García, S., Laurikainen, E., & Salo, H. 2015, A&A, 582, A86
  * Ho et al. (2019) Ho, I. T., Kreckel, K., Meidt, S. E., et al. 2019, ApJ, 885, L31
  * Hoare (2005) Hoare, M. G. 2005, Ap&SS, 295, 203
  * Hopkins (2013) Hopkins, P. F. 2013, MNRAS, 430, 1653
  * Hughes et al. (2016) Hughes, A., Meidt, S., Colombo, D., et al. 2016, in From Interstellar Clouds to Star-Forming Galaxies: Universal Processes?, ed. P. Jablonka, P. André, & F. van der Tak, Vol. 315, 30–37
  * Johnson et al. (2016) Johnson, L. C., Seth, A. C., Dalcanton, J. J., et al. 2016, ApJ, 827, 33
  * Johnson et al. (2017) Johnson, L. C., Seth, A. C., Dalcanton, J. J., et al. 2017, ApJ, 839, 78
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kennicutt et al. (1989) Kennicutt, Robert C., J., Edgar, B. K., & Hodge, P. W. 1989, ApJ, 337, 761
  * Kennicutt & Garnett (1996) Kennicutt, Robert C., J. & Garnett, D. R. 1996, ApJ, 456, 504
  * Kennicutt et al. (2008) Kennicutt, Robert C., J., Lee, J. C., Funes, J. G., et al. 2008, ApJS, 178, 247
  * Kennicutt & Evans (2012a) Kennicutt, R. C. & Evans, N. J. 2012a, ARA&A, 50, 531
  * Kennicutt & Evans (2012b) Kennicutt, R. C. & Evans, N. J. 2012b, ARA&A, 50, 531
  * Kewley & Dopita (2002) Kewley, L. J. & Dopita, M. A. 2002, ApJS, 142, 35
  * Kewley et al. (2001) Kewley, L. J., Heisler, C. A., Dopita, M. A., & Lumsden, S. 2001, ApJS, 132, 37
  * Kim et al. (2021) Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2021, MNRAS, 504, 487
  * Knapen (1998) Knapen, J. H. 1998, MNRAS, 297, 255
  * Knapen et al. (1993) Knapen, J. H., Arnth-Jensen, N., Cepa, J., & Beckman, J. E. 1993, AJ, 106, 56
  * Kreckel et al. (2020) Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2020, MNRAS, 499, 193
  * Kreckel et al. (2019) Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2019, ApJ, 887, 80
  * Kroupa & Weidner (2003) Kroupa, P. & Weidner, C. 2003, ApJ, 598, 1076
  * Kruijssen (2012) Kruijssen, J. M. D. 2012, MNRAS, 426, 3008
  * Kruijssen & Longmore (2014) Kruijssen, J. M. D. & Longmore, S. N. 2014, MNRAS, 439, 3239
  * Kruijssen et al. (2019) Kruijssen, J. M. D., Schruba, A., Chevance, M., et al. 2019, Nature, 569, 519
  * Kruijssen et al. (2018) Kruijssen, J. M. D., Schruba, A., Hygate, A. P. S., et al. 2018, MNRAS, 479, 1866
  * Krumholz (2014) Krumholz, M. R. 2014, Phys. Rep, 539, 49
  * Krumholz et al. (2015) Krumholz, M. R., Adamo, A., Fumagalli, M., et al. 2015, ApJ, 812, 147
  * Lang et al. (2020) Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122
  * Larsen (2002) Larsen, S. S. 2002, AJ, 124, 1393
  * Lawton et al. (2010) Lawton, B., Gordon, K. D., Babler, B., et al. 2010, ApJ, 716, 453
  * Lee et al. (2021) Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2021, arXiv e-prints, arXiv:2101.02855
  * Lee et al. (2011) Lee, J. H., Hwang, N., & Lee, M. G. 2011, ApJ, 735, 75
  * Lelièvre & Roy (2000) Lelièvre, M. & Roy, J.-R. 2000, AJ, 120, 1306
  * Leroy et al. (2021) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, arXiv e-prints, arXiv:2104.07739
  * Liu et al. (2013) Liu, G., Calzetti, D., Kennicutt, Robert C., J., et al. 2013, ApJ, 772, 27
  * Luridiana et al. (2015) Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42
  * Maíz Apellániz & Úbeda (2005) Maíz Apellániz, J. & Úbeda, L. 2005, ApJ, 629, 873
  * Makarov et al. (2014) Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13
  * Mascoop et al. (2021) Mascoop, J. L., Anderson, L. D., Wenger, T. V., et al. 2021, ApJ, 910, 159
  * Massey et al. (2005) Massey, P., Puls, J., Pauldrach, A. W. A., et al. 2005, ApJ, 627, 477
  * Meidt et al. (2021) Meidt, S. E., Leroy, A. K., Querejeta, M., et al. 2021, ApJ, 913, 113
  * Murphy et al. (2011) Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67
  * O’Donnell (1994) O’Donnell, J. E. 1994, ApJ, 422, 158
  * Oey (1996) Oey, M. S. 1996, ApJ, 465, 231
  * Oey & Clarke (1998a) Oey, M. S. & Clarke, C. J. 1998a, AJ, 115, 1543
  * Oey & Clarke (1998b) Oey, M. S. & Clarke, C. J. 1998b, AJ, 115, 1543
  * Oey et al. (2007) Oey, M. S., Meurer, G. R., Yelda, S., et al. 2007, ApJ, 661, 801
  * Pellegrini et al. (2010) Pellegrini, E. W., Baldwin, J. A., & Ferland, G. J. 2010, ApJS, 191, 160
  * Pellegrini et al. (2012) Pellegrini, E. W., Oey, M. S., Winkler, P. F., et al. 2012, ApJ, 755, 40
  * Pellegrini et al. (2020) Pellegrini, E. W., Rahner, D., Reissl, S., et al. 2020, MNRAS, 496, 339
  * Pessa et al. (2021) Pessa, I., Schinnerer, E., Belfiore, F., et al. 2021, A&A, 650, A134
  * Pietrinferni et al. (2004) Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168
  * Pilyugin & Grebel (2016) Pilyugin, L. S. & Grebel, E. K. 2016, MNRAS, 457, 3678
  * Querejeta et al. (2021) Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, arXiv e-prints, arXiv:2109.04491
  * Rahner et al. (2017) Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 4453
  * Rahner et al. (2019) Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2019, MNRAS, 483, 2547
  * Rand (1992) Rand, R. J. 1992, AJ, 103, 815
  * Rand (1993) Rand, R. J. 1993, ApJ, 410, 68
  * Reina-Campos & Kruijssen (2017) Reina-Campos, M. & Kruijssen, J. M. D. 2017, MNRAS, 469, 1282
  * Rice et al. (2016) Rice, T. S., Goodman, A. A., Bergin, E. A., Beaumont, C., & Dame, T. M. 2016, ApJ, 822, 52
  * Rosolowsky et al. (2021) Rosolowsky, E., Hughes, A., Leroy, A. K., et al. 2021, MNRAS, 502, 1218
  * Rosolowsky et al. (2007) Rosolowsky, E., Keto, E., Matsushita, S., & Willner, S. P. 2007, ApJ, 661, 830
  * Rousseau-Nepton et al. (2018) Rousseau-Nepton, L., Robert, C., Martin, R. P., Drissen, L., & Martin, T. 2018, MNRAS, 477, 4152
  * Rozas et al. (1996) Rozas, M., Beckman, J. E., & Knapen, J. H. 1996, A&A, 307, 735
  * Sánchez et al. (2012) Sánchez, S. F., Rosales-Ortega, F. F., Marino, R. A., et al. 2012, A&A, 546, A2
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schruba et al. (2019) Schruba, A., Kruijssen, J. M. D., & Leroy, A. K. 2019, ApJ, 883, 2
  * Scoville et al. (2001) Scoville, N. Z., Polletta, M., Ewald, S., et al. 2001, AJ, 122, 3017
  * Thilker et al. (2000) Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2000, The Astronomical Journal, 120, 3070
  * Thilker et al. (2002) Thilker, D. A., Walterbos, R. A. M., Braun, R., & Hoopes, C. G. 2002, AJ, 124, 3118
  * Thilker et al. (2021) Thilker, D. A., Whitmore, B. C., Lee, J. C., et al. 2021, arXiv e-prints, arXiv:2106.13366
  * Tress et al. (2020) Tress, R. G., Smith, R. J., Sormani, M. C., et al. 2020, MNRAS, 492, 2973
  * Tress et al. (2021) Tress, R. G., Sormani, M. C., Smith, R. J., et al. 2021, MNRAS, 505, 5438
  * Vacca (1994) Vacca, W. D. 1994, ApJ, 421, 140
  * van Zee (2000) van Zee, L. 2000, AJ, 119, 2757
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409
  * Wei et al. (2020) Wei, P., Zou, H., Kong, X., et al. 2020, PASP, 132, 094101
  * Weidner et al. (2004) Weidner, C., Kroupa, P., & Larsen, S. S. 2004, MNRAS, 350, 1503
  * Whitmore et al. (2014) Whitmore, B. C., Chandar, R., Bowers, A. S., et al. 2014, AJ, 147, 78
  * Whitmore et al. (2021) Whitmore, B. C., Lee, J. C., Chandar, R., et al. 2021, MNRAS
  * Youngblood & Hunter (1999) Youngblood, A. J. & Hunter, D. A. 1999, ApJ, 519, 55
  * Zhang et al. (2017) Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217

##  Appendix A Nebulae spatial masks and metallicity gradients

In this section, we present the spatial masks of the H ii regions for 18
galaxies of the PHANGS–MUSE sample (from Fig. 12 to Fig. 29) highlighting the
H ii region footprints and the different galaxy environments. At the end of
this section in Fig. 30, we also show the best-fitting radial metallicity
gradients for all the galaxies in our sample. They provide the measurement of
the metallicity at the median galactocentric radius covered by the MUSE
observations used in Sec. 4.2.

Figure 12: H ii regions and environments for IC5332. The figure shows the
H$\alpha$ emission in the background, color coded according to the color
scheme on the right, with overlaid the borders of the H ii regions in our
catalog. The centers of the nebulae which have been discarded by our selection
criteria are marked with crosses. In the lower left corner, the black circle
indicates the PSF of the MUSE observations while the black line marks a
physical scale corresponding to $1$ kpc. Both the H ii regions and the
discarded nebulae are color coded according to our definition of environments
as outlined by the color scheme at the bottom.  Figure 13: H ii regions and
environments for NGC0628. As in Fig. 12.  Figure 14: H ii regions and
environments for NGC1087. As in Fig. 12.  Figure 15: H ii regions and
environments for NGC1300. As in Fig. 12.  Figure 16: H ii regions and
environments for NGC1365. As in Fig. 12.  Figure 17: H ii regions and
environments for NGC1385. As in Fig. 12.  Figure 18: H ii regions and
environments for NGC1433. As in Fig. 12.  Figure 19: H ii regions and
environments for NGC1512. As in Fig. 12.  Figure 20: H ii regions and
environments for NGC1566. As in Fig. 12.  Figure 21: H ii regions and
environments for NGC1672. As in Fig. 12.  Figure 22: H ii regions and
environments for NGC2835. As in Fig. 12.  Figure 23: H ii regions and
environments for NGC3351. As in Fig. 12.  Figure 24: H ii regions and
environments for NGC3627. As in Fig. 12.  Figure 25: H ii regions and
environments for NGC4254. As in Fig. 12.  Figure 26: H ii regions and
environments for NGC4303. As in Fig. 12.  Figure 27: H ii regions and
environments for NGC4535. As in Fig. 12.  Figure 28: H ii regions and
environments for NGC5068. As in Fig. 12.  Figure 29: H ii regions and
environments for NGC7496. As in Fig. 12.  Figure 30: Fit of the radial