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

We use an unprecedented sample of about 23000 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 α 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, α , in our sample is = 1 . 73 with a σ of 0 . 15. We ﬁnd that α decreases with the galaxy’s star formation rate surface density, Σ 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 ﬁnd that no signiﬁcant 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 ﬁnd suggestive evidence that they are more evident for galaxies with stronger spiral arms. Furthermore, we ﬁnd systematic variations in α between samples of H ii regions with a high and low ionization parameter, q , and argue that they are driven by the aging of H ii regions.


Introduction
Hα 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 whose 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), though we should not forget that star formation currently ongoing within dense, dust-enshrouded cores may remain hidden for about 3 Myr (Kim et al. 2021). Since the Hα 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 function and the luminosity function of starforming regions are 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 The catalog of HII regions is only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http: //cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/658/A188 this expectation, measuring LF slopes close to −2 with minor variations of ±0.2−0.5 among galaxies (e.g., UV: Cook et al. 2016, Hα: Kennicutt et al. 1989, Paα: 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α emission line to ground-based observations favored the study of the so-called nebular LF (i.e., built using the Hα luminosity of H ii regions) in numerous extra-galactic Hα surveys (e.g., Knapen 1998;Thilker et al. 2002;Kennicutt et al. 2008). Thanks to studies in the Milky Way (MW) and the 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 ∼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, but 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, and 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 in 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 in the LF slope with the galaxy's star formation rate (SFR), morphology, or local environmental properties can elucidate whether, for example, environment affects the demographics of young stellar populations and may provide new clues that help us better understand 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-MUSE 1 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. 2019Kreckel et al. , 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 for studying the H ii region LF, marking a turning point in terms of statistics and ability to resolve H ii regions. Furthermore, the availability of PHANGS-ALMA ) and PHANGS-HST ) observations opens up the possibility of comparing 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. 2022;Whitmore et al. 2021) and investigating 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, SFR, stellar mass, and gas-phase metallicity) and galactic environment (e.g., spiral arm versus inter-arm) across the entire PHANGS-MUSE sample. The paper is structured as follows: Sect. 2 briefly describes the PHANGS-MUSE data set, the data reduction strategy, and the ancillary data that have been used. In Sect. 3 we describe the source identification technique and present our H ii region catalogs. In Sect. 4 we present the nebular LF of our galaxies and the fitting technique (Sect. 4.1), and we investigate variations in the LF slope as a function of global properties across the sample (Sect. 4.2) or for different populations of H ii regions within individual galaxies (Sect. 4.3).
In Sect. 5 we assess the effects of completeness and crowding on our results. In Sect. 6 we discuss our results and what drives variations in the LF slope. We finally summarize our main findings in Sect. 7 and highlight future prospects in Sect. 8.

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.

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 55 degrees). We summarize their main properties in Table 1. The galaxies are located close to the so-called star formation main sequence and span a stellar mass range of log(M /M ) = 9.4−11. The galactocentric radii used in this paper have been deprojected according to the inclinations and position angles reported in Table 1, 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 1.

PHANGS-MUSE data reduction and analysis
The PHANGS-MUSE data set (Emsellem et al. 2021) is the result of an extensive (∼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 Very Large Telescope (VLT), with the addition of archival observations for NGC 0628 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 ×1 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 2 , 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.2 m Wide Field Imager and with the Du Pont Direct CCD camera (Razza et al., in prep.) were used as a reference to astrometrically align and flux calibrate the MUSE data. The complete data reduction provides a sky-subtracted, fluxcalibrated, mosaicked data cube for each galaxy with a spaxel size of 0.2 . The point spread function (PSF) of the observations was determined for each pointing, and for each galaxy the largest value (reported in Table 1) 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 was 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).  Col. 11), angular resolution of the copt data (Col. 12), the parsec per arcsecond scaling factor (Col. 13), and the maximum deprojected galactocentric radius in units of R 25 (Col. 14). References.
(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) 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)  . 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 make use of the copt Hα maps to localize H ii regions and the native MUSE mosaicked data cubes to extract their integrated spectra (see Sect. 3).

Environmental and foreground stars masks
Morphological masks delimiting stellar structures by using Spitzer near-infrared (NIR) 3.6 µm imaging are available for the PHANGS-MUSE galaxies from Querejeta et al. (2021). We used 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 "disk" for the galaxies without spiral arms. These five environments can be visualized for NGC 4321 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 that 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.

Ionized nebulae identification
To isolate H ii regions across our sample, we make use of the PHANGS-MUSE copt Hα emission line maps provided by the DAP procedure summarized in Sect. 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α narrowband 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α maps to the mean Hα flux of the pixels with surface brightness below Σ Hα < 1 × 10 −17 erg s −1 arcsec −2 and the detection threshold to 3σ above the background, with σ 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α map. We consider the total integrated flux of the seed region and impose a S/N cut of 50 (using the Hα emission line error maps to assess the noise) and a Σ Hα cut above 3σ 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α equivalent width, Hα 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α 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 (FWHM PSF ≤ 1.2 ), after visual inspection we select a termination gradient of 2.43 × 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 Sect. 4.1 of this paper. Excluding the regions A188, page 4 of 41 overlapping with foreground stars (see Sect. 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 Sect. 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]λ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α 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α emission outside the ionized nebulae footprints can represents 20−50% of the total Hα emission (see Belfiore et al. 2021).

The final H ii region catalog
In order to obtain our fiducial H ii regions catalog, we excluded the following five types of nebulae. The first are nebulae located in galaxy centers (as defined in Sect. 2.3); this is meant to avoid areas with a high SFR surface density (e.g., starburst rings in NGC 1300, NGC 1512, NGC 1672, where the deblending of H ii regions is not possible at our resolution (i.e., typically a few parsec in radius; see, e.g., Barnes et al. 2020), and contamination by active galactic nuclei (AGN; i.e., NGC 1365, NGC 1512, NGC 1566, and NGC 1672 show evidence of hosting low-luminosity AGN when looking at emission line maps). This cut comprises about 1.5% of the nebulae.
Second are 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)  ii]λ6731Å emission lines; this is meant to ensure a robust estimate of both the dust correction based on the Hα/Hβ 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.
Third are nebulae with velocity dispersion σ > 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 candidates. This cut comprises about 1% of the nebulae.
Fourth are nebulae whose geometric center matches within 0.5 the positions of planetary nebulae identified using [O iii] emission line maps (Scheuermann et al., in prep.). This cut comprises about 0.6% of the nebulae. Fifth are nebulae where the distance of the geometric center from the edges of the MUSE mosaic FoV is less than one FWHM PSF . This cut comprises about 2% of the nebulae.
We note that the S/N of the Hα and Hβ lines is usually high, more specifically, 99% of the region spectra have a S /N 32 for the Hα line and 12 for the Hβ 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 NGC 4321, highlighting also the different galactic environments defined in Sect. 2.3. Analogous figures for the remaining galaxies in our sample are presented in Appendix A.
We assumed a screen geometry and used pyneb 3 (Luridiana et al. 2015) to correct line fluxes for dust extinction via the Hα/Hβ ratio, adopting the O'Donnell (1994) reddening law with R V = 3.1 and a theoretical Hα/Hβ = 2.86.
The extinction-corrected emission line luminosities of the H ii regions were then computed using the distances reported in Table 1. For every H ii region in our catalog, we also estimated the gas-phase metallicity and the gas ionization parameter by using extinction-corrected emission line fluxes. The gas-phase metallicity O/H was calculated using the Pilyugin & Grebel (2016)  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 leastsquare linear fitting of the trend between 12 + log(O/H) and the deprojected galactocentric radius (see Fig. A.19). 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 × c = Q(H 0 )/4π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) (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]λ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α and Hβ observed fluxes, the Hα extinction-corrected flux, the ionization parameter, and an environmental flag for all the H ii regions discussed in the current paper is only available at the CDS.

The H ii region luminosity function
We make use of the extinction-corrected Hα luminosities of our H ii regions to obtain the nebular 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, α. Several different methods for deriving the slope 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, but the noise of the logarithm of a histogram is not Gaussian. Furthermore, the choice of the binning scheme (e.g., equalluminosity bins versus equal-number bins) and the definition of the bin center can affect the LF modeling and add uncertainties to α (see, e.g., Maíz Apellániz & Úbeda 2005; Cook et al. 2016).
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 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 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 min and α of a given LF in a statistically robust way (see Sect. 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 The algorithm performs an MLE fit by recursively fixing L min equal to each empirical data point. For a given L 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 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 min to an interval of ±1σ around the median of the distribution (i.e., where, visually, the LFs start to flatten) and the LF slope to the interval of α = 1−3. The error on L min is taken to be the standard deviation of the L 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 α, we consider the trend of the likelihood estimator as a function of α for the models with the best L min . This trend can be well represented by a Gaussian that peaks at the best α and we compute the error on α 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 2. The slopes we find are in the interval of α = 1.5−2 with a mean value of 1.73 and a σ 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 NGC 7496) and that our fits can be considered robust against biases due to low-number statistics (see Sect. 5 for further discussion). In Fig. 2, we also mark the completeness limit of our LFs that will be discussed in Sect. 5. In Appendix C, we discuss how the results would change if Hα luminosities were not corrected for dust extinction (e.g., in narrowband 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 In this regime, the Hα 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 Large Magellanic Cloud) and are typical of late-type galaxies (Sc and Irr Hubble types; see, e.g., Kennicutt et al. 1989;. The ionizing luminosities corresponding to L min in each galaxy of our sample are reported in Table 2 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 is log(Q O7 V 0 [photons 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 min it is in the range of 0.05−6 in our sample. This testifies that our LF fits, extending ∼2−3 orders of magnitude above L min , are mainly probing H ii regions ionized by star clusters and associations, and only for a few targets (i.e., NGC 5068, IC 3352 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 Hα [erg s −1 ]) = 38.6 ± 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 Hα [erg 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.
A188, page 6 of 41 The number of putative density bounded H ii regions in our galaxies (i.e., L Hα ≥ 10 38.6 erg s −1 ) varies from a few up to about 350 and is reported in Table 2. 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).
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 Hα ≥ 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., NGC 1365 and NGC 4321) the variations are within the errors. However, looking at the LFs of our sample in Fig. 2, the steepening at log(L Hα [erg s −1 ]) > 38.6±0.1 and the glitch at log(L Hα [erg s −1 ]) = 38.6 ± 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 take the advantage of this simpler parametriza-tion of the LF to directly compare the LF slopes of different galaxies and of H ii region subsamples within individual galaxies.

LF variations among galaxies
To look at what may drive galaxy-to-galaxy variations in the LF slope in our sample, we verified that α does not have any dependence on galaxy inclination reported in Table 1 and then estimate a few global galaxy parameters, excluding the areas defined as galaxy centers (see Sect. 2.3), to be in line with the H ii regions probed by our LFs. More specifically, we estimate the following four properties. The first is 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. A.19). The second is the total stellar mass (M * ) extracted by integrating the PHANGS-MUSE stellar mass maps produced by the DAP.
The third is the total SFR extracted by integrating the PHANGS-MUSE Hα 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) A188, page 7 of 41 A&A 658, A188 (2022)

Galaxy
Slope IMF, with a Salpeter slope (α * = −2.35) from 1 to 100 M and α * = −1.3 for 0.1−1 M , and solar metallicity. The fourth is the specific star formation rate (sSFR; sSFR = SFR/M * ) and the SFR surface density (Σ 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 ∼100 pc scales carried out by Pessa et al. (2021), we also consider the offset ∆ of individual galaxies with respect to the trend of the full sample for the resolved star formation main sequence (rSFMS; Σ M * versus Σ SFR ), the Kennicutt-Schmidt relation (rKS; Σ H 2 versus Σ SFR ), and the molecular gas main sequence (rMGMS; Σ M * versus Σ H 2 ). The offsets ∆ rSFMS, ∆ rKS, and ∆ 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 NGC 0628 due to the quality of the MUSE data while, due the lack of significant CO(2-1) emission, IC 5332 does not have measurements of ∆ rKS and ∆ rMGMS. The ∆ rSFMS, ∆ rKS, and ∆ rMGMS offsets are interesting quantities because, by construction of the star formation scaling relations, they measure the relative variations in the sSFR (Σ SFR /Σ M * ), the gas depletion time (t dep = Σ H 2 /Σ SFR ), and the molecular gas fraction (Σ H 2 /Σ 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 α and the aforementioned parameters, with the addition of the morphological T type, are shown in Fig. 4 together with their Spearman correlation coefficient (ρ) and their 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 4. 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 Sect. 6.1, we compare our results to those studies that, as 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 (1998), 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, for example, FUV observations, probing H ii regions with ages less than 100 Myr, are expected to deliver a steeper LF compared to Hα observations, typically probing H ii regions younger than 10 Myr, and our comparison remains qualitative.
We find a weak correlation between α 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 A188, page 8 of 41 ( 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 α 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 α of spirals of the same type in their sample is comparable to the magnitude of the trend they report to exist between α and morphology. The lack of a correlation between α 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(O/H) > 8.3) and stellar masses (i.e., log M * [M ] > 9). Cook et al. (2016) found no trend between α and gas-phase metallicity in the range 7.2 < 12 + log(O/H) < 9.2 and only a moderate trend between α and M * in the range 6 < log M * [M ] < 11. The parame-ter 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 α (Figs. 4d-f). SFR, sSFR, and Σ SFR have a moderate, strong, and very strong correlation with α, respectively. The presence of a tighter correlation with Σ 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 anticorrelate with α, 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 Looking at where our galaxies lie in the plane of the three classical star formation scaling relations (Fig. 4g-i), we find that α anticorrelates strongly with ∆ rSFMS and ∆ rMGMS, and moderately with ∆ 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-togalaxy scatter in our sample (Pessa et al. 2021). While the trend with ∆ rSFMS confirms the trend we find with sSFR, ∆ rMGMS adds information on the molecular gas fraction and indicates that galaxies with a lower fraction of molecular gas tend to have steeper LFs. A188, page 9 of 41   Notes. The 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 region subsamples in spiral arms areas (Cols. [2][3], in inter-arms areas (Cols. [4][5], in the inner star-forming disk (Cols. 6-7), in the outer star-forming disk (Cols. [8][9], with low ionization parameter (Cols. 10-11), and with high ionization parameter (Cols. [12][13]. 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 beyond the scope of the current paper.

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 subsamples. 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 subsamples and fit their LFs while keeping L min fixed to the value found in Sect. 4.1 for each parent sample. All the details of the LFs fits performed in this section are reported in Table 4, while the actual fits are shown in Appendix B.

Spiral arm versus 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. 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 2).
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 α 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 NGC 1566 and NGC 1300); 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 Sect. 6.2.

Inner versus outer disk
There exists evidence that the conditions under which SF happens in the outer disk of spiral galaxies are different from those of the inner disk, as suggested, for example, by an observed steepening of the GMC MF (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 in 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 A188, page 11 of 41 A&A 658, A188 (2022) Fig. 5. Slope of the LF for H ii regions located in spiral arms (red points) and inter-arm areas (blue points). 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).
gas-phase metallicity), we study the LF at small and large galactocentric radii. To mitigate effects due to different subsample 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 subsamples of equal sizes. It should be noted that the area covered by our MUSE observations is largely contained within ∼1R 25 of our targets, as can be seen looking at the r max values reported in Table 1. 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 ). Figure 6 summarizes the results of our fits of the LF slope in the inner and outer star-forming disk (shown in Fig. B.2). The median H ii region galactocentric radii, used to draw the H ii region subsamples, varies between ∼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 subsamples (not shown here). Only two out of 19 galaxies (i.e., NGC 1385 and NGC 1566) 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. A.19) do not affect the LF slope.

High versus low ionization parameter
Lastly, we investigate variations in 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 subsamples 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 subsamples of H ii regions for each galaxy and they both appear similarly distributed across the galaxy disks. Figure 7 summarizes the results of the LF fit for the population of H ii regions with high and low ionization parameter (shown in Fig. B.3). 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 num-  ber of targets (e.g., IC 5332, NGC 0628, NGC 2835, NGC 3627, and NGC 5068) the H ii regions with high and low q constitute, respectively, the bulk of the bright and faint H ii region populations. The interpretation of these results is further discussed in Sect. 6.3.

Completeness, blending, and selection effects
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α 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α flux outside the ionized nebula footprints in the Hα emission line maps, which we refer to as the Hα 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 Σ Hα equal to the 90th percentile (i.e., slightly larger than 1σ of the distribution) of the Hα diffuse emission distribution. In this way, we take into account the bright tail of the diffuse Hα background against which H ii regions can remain undetected. It is worth noting that tests A188, page 12 of 41 injecting mock H ii regions are a good way forward to more robustly assess the completeness limit in the future.
As shown in Fig. 2, the completeness limit estimated in this way lies below the best fitting L 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α 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 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 α 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 multiplied the FWHM PSF , which takes into account the spatial resolution of the data, by √ |Σ SFR |; this gives us a handle on the mean separation between H ii regions. We used this parameter as a proxy for blending and check its relation with α; 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) have shown that spatial resolution does not significantly change the slope of the LF at resolutions finer than few hundred parsecs, as is the case of our data (i.e., the coarser physical resolution in our sample is ∼110 pc for NGC 1512 and NGC 1365). This indicates that, although present, the effect of blending is not the main driver of the correlations we report in Sect. 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 Sect. 3.2), have no significant effect on the measurements of the LF slope presented in Sect. 4.1.

Discussion: What sets the LF slope?
6.1. The SFR surface density In Sect. 4.2 we find that in general α correlates better with the global star formation properties of our galaxies and especially with Σ SFR . We do not find any clear trends within the T-type, [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 starforming 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(Hα) 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 α = −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 lownumber statistics increases the scatter in the measured LF slope but does not drive systematic changes of α in a specific direction.
Moreover, our LFs are built with H ii region samples containing at least ∼500 objects, and the LF fits presented in Sect. 4.1 rely on at least 200 objects. For this reason, we expect the stochastic scatter of α to be lower in our sample compared to previous studies that often relied on smaller H ii region samples.
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 Σ SFR . In Fig. 4f,   A188, page 13 of 41 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α LFs probe different ages of H ii regions, 100 Myr and 10 Myr, respectively (see, e.g., Haydon et al. 2020). Oey & Clarke (1998) predicted that older H ii region populations (e.g., derived from FUV) should have steeper slopes than younger populations (e.g., derived from Hα). Overall, the average slope of the nebular LF (−2 ± 0.5; Kennicutt et al. 1989;Elmegreen & Salzer 1999), including our study (−1.73 ± 0.15), is compatible with what has been found in the FUV (−1.76±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 Σ 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 dep , and the molecular gas fraction.
As discussed in Sect. 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 ∆ 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 ∆ 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 parsec-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 (Γ ≡ CFR/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 Γ correlates well with the SFR surface density (see, e.g., Goddard et al. 2010;Adamo et al. 2011Adamo et al. , 2015Adamo et al. , 2020aCook et al. 2012;Johnson et al. 2016). The trend that we find in Fig. 4e between α and Σ SFR goes in the direction predicted by this scenario. It is worth noting that, with the available PHANGS-HST data , it will be possible to directly measure Γ 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 MF (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 α = 1.6−2.0. Steeper slopes may then be the result of incomplete sampling of the stellar IMF, which occurs if the stellar populations born within the parent GMCs have a low mass ( 10 4 M ; 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 Σ 1.4 gas , and therefore approximately linearly with Σ 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 M is reached at log (Σ SFR [M yr −1 kpc −2 ]) ≈ −2.6. Above this SFR surface density, we expect the H ii region LF to be close to the slope of the GMC MF (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 α declines with Σ SFR for log (Σ SFR [M yr −1 kpc −2 ]) −2.6 and stays approximately constant (or declines less steeply) for log (Σ SFR [M yr −1 kpc −2 ]) −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 α corresponding to the Cook et al. (2016) best fitting line at Σ SFR ≈ 2.6 and the mean α of our samples at Σ SFR > 2.6. Both values are remarkably close to α = 1.6, the slope of the GMC MF in the high gas and Σ 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.

The strength of spiral arms
In Sect. 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 NGC 1300, NGC 1365, NGC 1566, NGC 1672. The Monte Carlo simulations by Oey & Clarke (1998) 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 on 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. 2020a). 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 versus 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 granddesign 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 A188, page 14 of 41 Fig. 10. Change in the LF slope between inter-arm and spiral arm environments, ∆α ENV (derived as α i −α a as reported in Table 4), as a function of the CO(2-1) contrast between the 84th percentile and the reference level defined in Meidt et al. (2021) for six galaxies in our sample.
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. (2020Tress et al. ( , 2021 concluded that strong spiral arms gather molecular gas and clouds without dramatically affecting their properties, yielding similar overall SFRs.
In our case, we might thus expect that variations in cloud MFs 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 that show a significant environmental change in α are those in which the distribution of the H ii regions quite clearly traces 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., NGC 0628, NGC 4254, and NGC 4321), on the other hand, exhibit comparable α 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 . 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 84th 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 α 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 MFs in our sample (e.g., Rosolowsky et al. 2021;Hughes et al., in prep.) 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, for example, 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.
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 Hα 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.

The aging effect and the ionization parameter
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) ∝ (Q(t)n H ) 1/3 , where Q is the rate of ionizing photons produced by the central star or cluster and n 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. 2017Rahner et al. , 2019Pellegrini et al. 2020).
The LF slope can, in principle, be used to test this expectation by using the ionization parameter to draw subsamples of H ii regions in the way we did in Sect. 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 1998). This is the same effect discussed for environmental variations in the previous section. What we find in Sect. 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; for example, 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 it is now possible to infer, for example, the lifetime of GMCs t CO and the timescale over which stellar feedback acts t fb (see, e.g., Kruijssen et al. 2019;Chevance et al. 2020b;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 subsamples with high and low q, we can reasonably expect the magnitude of such variations ∆α 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. (2020b), Kim et al. (in prep.) derived t CO and t fb for most of the galaxies in our sample. In that work, H ii regions are traced via the Hα emission in the PHANGS-Hα data (Razza et al., in prep.), while GMCs are traced via the CO(2-1) emission in the PHANGS-ALMA data . It should be noted that for IC 5332 there is no significant CO(2-1) emission detected to apply this method and the t CO and t fb measurements for NGC 0628 and NGC 3627 are taken from Chevance et al. (2020b). As shown in Fig. 11, we find a weak-to-moderate and a moderate anticorrelation between ∆α q and t fb and t CO , respectively. 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.

Conclusions
In this paper we have studied the nebular LFs (i.e., built using the Hα luminosity of H ii regions) in the star-forming disks of the 19 nearby galaxies that make up the PHANGS-MUSE sample. Thanks to the exquisite spatial resolution (mean PSF 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α/Hβ Balmer decrement), gas-phase metallicity 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 MLE (Alstott et al. 2014). This method does not require us to bin the data, is especially suited for heav- Fig. 11. Change in LF slope between H ii regions with low and high ionization parameters (derived as α q< − α q> as reported in Table 4) as a function of the feedback timescales (upper panel) and the GMC lifetime (lower panel). For each panel we report the Spearman correlation coefficient, ρ, between the plotted quantities.
ily tailed distributions, and overcomes typical pitfalls of more commonly used methods (e.g., histograms combined with linear regression). We find an LF slope of α = −1.73 ± 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 mainly probe massive star formation happening in star clusters and associations. We report evidence that the LFs of a number of galaxies steepen at luminosities above log(L Hα [erg 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 uniformly analyze 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 (Sect. 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 slopes of the spiral arm and interarm LFs. On the other hand, we find environmental (i.e., arm versus inter-arm) LF variations for galaxies that show H ii regions preferentially located along spiral arms. We attribute these variations to the spiral arms increasing the molecular clouds' arminter-arm mass contrast and find suggestive evidence that they are more pronounced in galaxies with spiral arms that drive 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 Σ SFR have a flatter LF, meaning that their relative number of bright H ii regions is higher compared to galaxies with lower Σ SFR . This potentially connects to fundamental changes in the physics regulating star formation in galaxy disks. The trend we observe between α and Σ SFR is not driven by low-number statistics and provides further evidence for the prediction that the clustering of young stars is enhanced at high gas surface densities and star formation efficiencies (Kruijssen 2012 and models find that the maximum masses of the newborn 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 a high ionization parameter have a shallower LF compared to H ii regions with a low ionization parameter. Such variations are compatible with an aging effect, suggesting that H ii regions with a high q are the youngest star-forming regions, while the ones with a low q are more evolved. The lack of α trends or 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 that suggests that α variations between these H ii region populations are more evident when stellar feedback acts on shorter timescales and GMC lifetimes are shorter.

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 MF. Comparing GMC MF slopes with the H ii regions' LF slopes reported in this paper can bring new insights into if and 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) or cluster ionizing the H ii regions and help in disentangling what the main driver behind changes in the H ii regions' ionization parameter is. Furthermore, by providing a direct estimate of the clustered star formation efficiency (defined as the ratio between the CFR and the SFR; Γ ≡ CFR/SFR), it will be possible to compare Γ 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 Σ SFR .

Appendix A: Nebula 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. A.1 to Fig. A.18) highlighting the H ii region footprints and the different galaxy environments. At the end of this section in Fig. A.19, 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 Sect. 4.2.

Appendix B: Variations in the LF within single galaxies
In this section we show the LFs and their fits for the H ii region subsamples described in the main text in Sect. 4.3. Figure B.1 shows the LFs of the spiral arm and inter-arm areas, six of our 19 galaxies do not show evident spiral arms, for these galaxies we show the LF and best fitting models of the H ii regions across the entire disk but excluding bars (i.e., this is the reason for slight variations with respect to the fits shown in Fig. 2). Figure

Appendix C: The effect of dust correction and BPT cuts on the LF slopes
In this section we discuss in more detail the effect that dust extinction and selection criteria might play on the LF slopes presented in our paper. In Sect. 4.1, we discuss LFs built from extinction-corrected Hα luminosities of H ii regions. The spectroscopic nature of the MUSE data allow us to perform extinction correction via the Hα/Hβ Balmer decrement; however, in some cases (i.e., when the nebular LF is obtained from narrowband Hα imaging) this is not possible. To understand the effect of dust extinction, we fit the LF built with observed Hα luminosities and in Fig. C.1  to the fact that we detect a higher number of faint H ii regions.
The most significant drop in the number of ionized nebulae entering our H ii region catalog is due to the cuts on the line ratios applied using semi-empirical demarcation lines in the three classical BPT diagrams ("BPT cuts" hereafter; see Sect. 3.2). It is known that, depending on their gas-phase metallicity and ionization parameter, H ii regions can lie outside the boundaries that are commonly adopted to have a clean H ii region sample and avoid nebulae whose gas is ionized by mechanisms other than photoionization from young OB stars (e.g., AGN photoionization/shocks). As shown by Fig. C.3 and Fig. C.4, if we drop the BPT cuts when cleaning our ionized nebula catalogs and rerun the LF fit, we find no significant variation in the measured LF slope. Remarkably, this is true also for NGC 1365 for which there is clear evidence of an AGN ionization cone extending further out the central region . As can be seen by comparing the empirical LFs, the regions that are flagged by the BPT cuts tend to populate the fainter end of the LF and thus have negligible effects on the LF shape at the bright end. In addition, we note that the LF L min remains substantially unchanged aside from the case of NGC 3627, the only strongly interacting galaxy in our sample.
In conclusion, we find that both the dust correction and selection effects do not have a significant impact on the slopes of the LFs that we measure in Sect. 4. Galaxies are ordered by increasing stellar mass from left to right, and their names are indicated along the abscissa.   Fig. 2, red colors) and not applying the BPT cut (blue colors).
Galaxies are ordered according to increasing stellar mass from left to right, and their names are indicated along the abscissa.