EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 708 (April 2026) A&A, 708 (2026) A202 Full HTML
Open Access
Issue
A&A
Volume 708, April 2026
Article Number A202
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202556811
Published online 08 April 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Stellar feedback plays a crucial role in the evolution of the interstellar medium (ISM; Lopez et al. 2014; Klessen & Glover 2016; Krumholz et al. 2019). It shapes the surrounding gas, halts or triggers new star formation, and enriches the gas with metals. However, not all energy and momentum from stellar feedback is deposited in the H II regions (Hanish et al. 2010). Turbulence and radiative feedback can fragment the cloud (Kakiichi & Gronke 2021), opening channels through which some radiation escapes, while gas and dust absorb part of the stellar feedback. This contributes substantially to the diffuse ionised gas (DIG; Haffner et al. 2009).

A common way to constrain the leakage is to compare the number of predicted ionising photons emitted by individual stars or clusters, Q(H0), with the ionising photon flux inferred from the observed H α luminosity of an associated H II region, QHαMathematical equation: ${Q_{{\rm{H}}\alpha }}$ (Niederhofer et al. 2016). This yields the escape fraction of ionising photons, which is usually defined as fesc=Q(H0)−QHαQ(H0).Mathematical equation: ${f_{esc}} = {{Q\left( {{{\rm{H}}^0}} \right) - {Q_{{\rm{H}}\alpha }}} \over {Q\left( {{{\rm{H}}^0}} \right)}}.$(1)

The escape fraction is expected to evolve over time (Trebitsch et al. 2017) and depends on various local ISM conditions, such as gas density, metallicity, and star-gas geometry (Ramambason et al. 2022). Therefore, a large sample is necessary to marginalise over these variations or explore trends.

For H II regions located within a massive galaxy, any escaping radiation primarily impacts the surrounding galactic ISM, and it is clear that within galaxies there are kiloparsec-scale consequences of stellar feedback (Zurita et al. 2000; Belfiore et al. 2022). Connecting this feedback with the galactic escape fraction, which is vital to understanding cosmic reionisation (Paardekooper et al. 2011; Mitra et al. 2013; Japelj et al. 2017; Kado-Fong et al. 2020; Ma et al. 2020; Ramambason et al. 2020; Wang et al. 2021; Chisholm et al. 2022; Saldana-Lopez et al. 2022; Flury et al. 2022a,b), is challenging.

Different measurement techniques produce a wide range of estimates for fesc. One of the first attempts to apply a method similar to ours, i.e. to analyse individual ionising sources in nearby galaxies, was conducted by McLeod et al. (2020) for NGC 300. Because the galaxy is nearby (2.09 Mpc; Jacobs et al. 2009), they were able to use the Multi Unit Spectroscopic Explorer (MUSE) to identify individual O-stars as the source of the ionising radiation. By matching them to two H II regions, they measured escape fractions of 28 ± 6% and 51 ± 1%, respectively. Della Bruna et al. (2021) analysed eight H II regions in NGC 7793 and report a value of fesc=67−12+8%Mathematical equation: ${f_{{\rm{esc}}}} = 67_{ - 12}^{ + 8}\% $. Della Bruna et al. (2022), applied their technique to M 83 and identify 541 matches between H II regions and young star clusters. However, only 13 percent of them return a positive fesc, i.e. for most objects, the photons emitted by the stars were insufficient to explain the observed H α emission. Recently, Teh et al. (2023) studied 42 H II regions in NGC 628 and find an escape fraction of 9−6+6%Mathematical equation: $9_{ - 6}^{ + 6}\% $.

Studies of even closer galaxies do not yield more consistent results, although individual ionising stars can be characterised and rich multi-wavelength constraints are available. Pellegrini et al. (2012) analysed 13 H II regions in the Large Magellanic Cloud with known ionising stars and infer typical escape fractions ranging from 30 to 60%. However, in the extreme starburst region of the Tarantula Nebula, Doran et al. (2013) measured extremely low escape fractions of 5 to 10%. The dwarf galaxies Holmberg I (Egorov et al. 2018) and Sextans A (Gerasimov et al. 2022) have well-characterised O star populations and escape fractions of ∼60%, but larger dwarf galaxy samples, in which modelling included infrared constraints, typically return much lower 10 to 20% escape fractions.

In recent years, simulations have made considerable progress to include more realistic treatments of stellar feedback physics, but models still predict a range of escape fractions as a function of both local conditions and time. Cloud-scale simulations with simplified feedback prescriptions (e.g. Dale et al. 2012, 2013) indicate that the escape fraction can vary widely over time as photoionisation reshapes the cloud structure. Radiation-hydrodynamic models of turbulent clouds, some including magnetic fields (Howard et al. 2017, 2018; Kim et al. 2019, 2021b), also find strongly time-dependent escape fractions that can briefly approach 100%, but average only a few per cent, up to ten per cent, over a cloud’s lifetime. Semi-analytic feedback models (Rahner et al. 2017) predict similar trends, with escape fractions rising as the feedback clears the gas, while galaxy-scale simulations (Tacchella et al. 2022) infer global values of the order of a few per cent.

Large field-of-view integral field unit (IFU) spectroscopy instruments, such as MUSE or the Spectromètre Imageur à Transformée de Fourier pour l’Etude en Long et en Large de raies d’Emission (SITELLE), now enable studies of the ionised gas using multiple optical emission-line diagnostics, isolating individual H II regions while imaging the full star-forming disc for a large sample of nearby galaxies. In combination with space-based imaging, for example, the Hubble Space Telescope (HST), we can obtain an accurate census of the ionising sources in galaxies beyond the Local Group. The Physics at High Angular resolution in Nearby GalaxieS (PHANGS)1 survey, with its combination of MUSE and HST imaging, is ideally suited to systematically study star formation in nearby galaxies (Leroy et al. 2021a,b; Emsellem et al. 2022; Lee et al. 2022, 2023). In this work, we used the large sample of H II regions and stellar associations compiled in Scheuermann et al. (2023) to calculate fesc and search for systematic variations with local physical conditions and star cluster properties.

This paper is organised as follows. In Sect. 2, we present the catalogues and models used in this work. In Sect. 3, we use these data to compute the escape fractions of the H II regions and discuss trends with physical properties. We conclude in Sect. 4.

2 Data and models

Our sample comprises 19 PHANGS galaxies observed with both MUSE and HST and is listed in Table 1. Scheuermann et al. (2023) describe the matching of nebulae and ionising sources in detail; hereafter, we refer to this work as Paper I. In the following, we provide a short summary of the individual data sets, how we combined them to form matched catalogues, and the stellar models used during the analysis.

Table 1

Properties of the PHANGS galaxy sample.

2.1 MUSE HII regions

The MUSE instrument is a powerful optical IFU spectrograph mounted on the Very Large Telescope in Chile (Bacon et al. 2010). It has a 1 × 1 arcmin2 field of view and covers the wavelength range from 4800 to 9300 Å in the nominal mode. The PHANGS–MUSE survey (PI: Schinnerer) provides wide-field mosaics of 19 nearby spiral galaxies. Emsellem et al. (2022) describe the analysis and the main data products. Among these are continuum-subtracted emission-line maps, constructed by fitting single Gaussians to a set of bright emission lines (e.g. H β, [O III]λ5007, H α, [N II]λ6583, [S II]λλ6717, 6731, [S III]λ9531), which form the basis for our analysis.

Based on the H α morphologies, Santoro et al. (2022) and Groves et al. (2023) used this dataset to construct an H II region catalogue. In the first step, we used HIIPHOT (Thilker et al. 2000) to draw the boundaries of the emission-line nebulae based on their H α brightness. This initial catalogue is composed of 31 497 nebulae, but includes other H α bright objects such as planetary nebulae (PNe) or supernova remnants (SNRs). We removed these by analysing diagnostic line ratios in the Bald-win–Phillips–Terlevich (BPT) diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987). We required regions to fall below the Kauffmann et al. (2003) diagnostic curve in the [O III]/H β versus [N II]/H α diagrams and below the Kewley et al. (2006) diagnostic curve in the [O III]/H β versus [S II]/H α diagrams, with S/N > 5 in all lines. Where [O I] is detected with S/N > 5, we further required that regions fall below the Kewley et al. (2006) diagnostic curve in the [O I]/H β versus [S II]/H α diagram.

For this work, we used only the objects that fell inside the HST coverage, which yielded 20 341 H II regions in the final catalogue. In addition, we derived properties such as metallicity, ionisation parameter, and density from the measured emission lines. Paper I provides a full summary of the catalogue and the derived properties, including detailed statistics.

We were primarily interested in the H II region H α luminosity, which we compared to the predicted ionising photon rate from the stars. The flux measurements for the nebulae are very precise, with a typical uncertainty of less than 1%, but we did not account for uncertainties arising from the determination of the boundaries. We corrected the measured fluxes for extinction. For the Milky Way, we used the extinction curve from O’Donnell (1994), E(B − V) from Schlafly & Finkbeiner (2011), and RV = 3.1. For the internal extinction, we adopted the same extinction curve and determined E(B − V) from the Balmer decrement, assuming H α/H β = 2.86. To convert the fluxes to luminosities, we used the distances listed in Table 1.

To compare these values with the emission from the stars, we converted the H α luminosity to the number of ionising photons. This was performed using a conversion factor such as Hummer & Storey (1987). Assuming Case B recombination, an electron temperature of 104 K, and a density of 102 cm−3 (well matched to our sample; Barnes et al. 2021), we obtain QHα/s−1=7.31⋅1011L(Hα)/ergs−1.Mathematical equation: ${Q_{{\rm{H}}\alpha }}/{{\rm{s}}^{ - 1}} = 7.31 \cdot {10^{11}}L({\rm{H}}\alpha )/erg\,{{\rm{s}}^{ - 1}}.$(2)

The observed temperatures in our H II regions can vary between 5000 and 10 000 K (Kreckel et al. 2022). We therefore ran INTRAT by Storey & Hummer (1995) to compute the conversion factor for a number of temperatures and densities. A lower temperature of 5000 K decreases the conversion factor by about 8%, and implies that we slightly underestimated the number of ionising photons (and hence overestimated the escape fraction) for some objects. An increase of the density to 103 cm−3 does not change the result by more than 0.3%.

2.2 HST stellar associations and star clusters

All 19 galaxies are included in the PHANGS–HST survey (Lee et al. 2022), providing a sharp ∼4 pc view of ionising stellar sources with coverage in five broadband filters: F275W (NUV), F336W (U), F438W (B), F555W (V), and F814 (I). From these images, we identify and distinguish both compact, gravitationally bound star clusters and physically larger, loosely bound stellar associations (Elmegreen 2008; Gouliermis 2018). In Paper I, we found that stellar associations provide a better one-to-one match to individual H II regions; however, compact star clusters are easier to identify and model, and they have been the focus of most previous studies constraining H II region escape fractions (e.g. Della Bruna et al. 2021; Teh et al. 2023). Therefore, in this work we include both star clusters and stellar associations in our analysis, flexibly allowing H II regions to contain multiple ionising sources.

The stellar association catalogue created by Larson et al. (2023) employs a watershed algorithm to join nearby young stars into a single object. Starting from either near-UV (NUV) or V tracer images, we used DOLPHOT to identify an initial set of point sources. We then used these to construct density maps, which we Gaussian smoothed to physical scales ranging from 8 to 64 pc. A watershed algorithm then groups multiple peaks to form a single larger structure or stellar association. We obtained five-band photometry of each catalogued object using the integrated fluxes of all stellar peaks within the association. As demonstrated in Paper I, our analysis yields similar results when we use different scales and tracer bands. We chose to use the 32 pc NUV catalogue for our analysis, as it has the largest overlap with our H II region catalogue. This catalogue contains 16 273 associations across 19 galaxies that fall within the MUSE coverage. This number differs slightly from that used in Paper I. The previous association catalogue was an internal release, and some associations in NGC 0628 lay outside the overlapping footprint of all five HST bands. This caused some issues with the spectral energy distribution (SED) fit; hence, those few objects are now removed from the catalogue.

For the star clusters, we used the catalogue created by Thilker et al. (2022) and Maschmann et al. (2024), which identifies marginally resolved sources in HST imaging and characterises them based on their morphology. The sources are detected with DOLPHOT and DAOSTARFINDER (v2.0, Dolphin 2000; Stetson 1987) and classified into Classes 1, 2, 3, and 4 using machine learning models (Wei et al. 2020; Hannon et al. 2023). Class 1 and 2 objects represent the compact clusters, while Class 3 objects are often multi-peaked and are typically referred to as ‘compact associations’ of stars. Class 4 objects are not considered to be star clusters; instead, they are classified as artefacts, such as background galaxies or single stars, and are therefore excluded. As with the stellar associations, we obtained integrated fluxes for each object to provide five-band photometry. The catalogue we used contains 24 529 star clusters that fall within the MUSE coverage (6593 of Class 1 or 2 and 17 936 of Class 3).

For both the stellar association and compact star cluster catalogues, age, mass, and reddening were derived by fitting the HST five-band SED as described in Turner et al. (2021). To do this, we used Code Investigating GALaxy Emission (CIGALE; Boquien et al. 2019) with theoretical single stellar population models from Bruzual & Charlot (2003), a metallicity of Z = 0.02, and a fully sampled Chabrier (2003) initial mass function (IMF). Since the publication of Paper I, Thilker et al. (2025) has updated the cluster catalogue. However, the catalogues are broadly consistent and the changes do not alter our results; for reproducibility, we carry out our analysis on the catalogue presented in Paper I.

We note a few caveats on the use of these catalogues to study the youngest stellar populations. We used a limited colour-space to derive the stellar properties, with five HST bands used to fit three stellar properties (mass, age, and reddening). Although this will be refined in the near future with the inclusion of additional James Webb Space Telescope (JWST) bands (Lee et al. 2023), the lack of suitable FUV coverage still presents a problem, especially for the youngest clusters. Although galaxies were observed with GALEX (Martin et al. 2005) and AstroSat (Hassani et al. 2024), the resolution of these instruments is insufficient for this analysis. Another issue that we cannot address with our current data is the contribution of Wolf-Rayet stars (Crowther 2007) or other massive single stars (Zastrow et al. 2013) to the ionising photon output.

2.3 Calculating the number of ionising photons

The published catalogues do not include the number of ionising photons, Q(H0), emitted by the stars. This number can depend on a variety of assumptions, such as the adopted stellar atmosphere model (Simón-Díaz & Stasińska 2008), stellar rotation (Levesque et al. 2012), and binarity (Lecroq et al. 2024). In this section, we use several codes to compute Q(H0) in order to investigate the systematics affecting this calculation.

Generally, the models adopt a given IMF to create a grid of stars and then compute the ionising photon flux as a function of age, based on different stellar models and atmospheres. If the IMF is fully sampled, the output scales linearly with mass. This scaling should be treated with caution, as it does not hold for lower-mass clusters (Cerviño & Luridiana 2006). We mitigate this problem by applying certain mass cuts, as discussed in Sect. 2.4. To propagate uncertainties in the ages, we assumed a Gaussian age distribution that is clipped at 1 Myr and drew a random sample for which we then computed Q(H0). Since the flux is expected to scale linearly with mass (for a fully sampled IMF), error propagation is straightforward. The uncertainty of the distance affects both the H α luminosity and the derived cluster mass in the same way and is therefore not included.

For the comparison, we considered results from CIGALE, STARBURST99 (the models without rotation, SB99 v00, and those with rotation, SB99 v40), and the binarity model from the Binary Population and Spectral Synthesis code (BPASS). Appendix A provides an overview of the details and the age dependence of the different models. Another popular model is Stochastically Lighting Up Galaxies (SLUG, da Silva et al. 2012), which treats individual properties as probability distributions rather than discrete values. This makes it difficult to include SLUG in the main analysis, so we refer to Appendix B for a simple comparison of our results with the SLUG models.

2.3.1 BC03

As previously stated, we used CIGALE (Boquien et al. 2019) to fit the SED and derive stellar properties. We also applied a self-consistent approach by using the ionising photon flux computed by this code. The catalogues published by Thilker et al. (2022) and Larson et al. (2023) do not include this quantity. We therefore ran CIGALE with the same input parameters in savefluxes mode to compute the ionising photon flux as a function of age. We then interpolated the value of Q(H0) for each stellar cluster and association using the previously fitted ages and masses. We used stellar models from Bruzual & Charlot (2003), together with an IMF from Chabrier (2003).

2.3.2 SB99

Another widely used stellar population synthesis model is STARBURST99 (Leitherer et al. 1999, 2010, 2014; Vázquez & Leitherer 2005). We selected a single burst of star formation as input, adopting a fully sampled Kroupa (2001) IMF, Geneva stellar models at solar metallicity (Ekström et al. 2012), and stellar atmospheres from Lejeune et al. (1997). We left the remaining parameters at their default values. We then used the stellar ages and masses to interpolate the outputs for our sample. We also consider a model with rotating stars, spinning at 40% of the break-up velocity, as an extreme limiting case.

2.3.3 BPASS

Most models track single-star evolution, but many stars form in binaries, which impacts their evolution (Sana et al. 2012). The difference in ionising photons becomes significant only after 8 Myr (Lecroq et al. 2024) and is therefore negligible for our purposes. To estimate the impact of binarity on Q(H0), we included BPASS by Eldridge & Stanway (2009) in our comparison. We used the fiducial binary model from version 2.1 (Eldridge et al. 2017), adopting a Kroupa IMF (Kroupa 2001) and an upper mass limit of 300 M⊙ at solar metallicity.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Comparison of different population synthesis models. We use a sample with random ages and masses and compute the ionising photon flux Q(H0) with BC03, SB99, and BPASS. The grey points represent the full sample, while the red points indicate a young and massive sub-sample (age ≤ 8 Myr and mass > 104 M⊙). The Spearman correlation coefficient ρ indicates good agreement between the different models across the robust subsamples, and all correlations are statistically significant (p-values < 0.05).

2.3.4 Comparison between the different models

To quantify the variations between the different models introduced in the previous sections, we compared the fluxes they predicted for a random sample. For this, we generated 10 000 clusters with random ages and masses, selected to be representative of the objects in our real sample. This provides order-of-magnitude understanding of how different processes impact our results.

Figure 1 presents a corner plot comparing the predicted ionising photon flux Q(H0). Overall, we find decent agreement between most models, but there are a few noticeable and systematic differences. The models with rotation predict higher ionising photon fluxes, as expected from Levesque et al. (2012). For BPASS, there are two distinct branches: during the first few million years the agreement is good, but around 6 Myr the difference between single stars and binaries becomes relevant.

Although the models differ significantly, these differences mostly affect specific parts of the sample, such as older or low-mass clusters. Young, massive objects show much smaller differences, with all models returning similar values. For the remainder of the paper, we use the STARBURST99 model without rotation as our fiducial model.

2.4 Matched catalogues

From the H II region catalogue and the two stellar catalogues (clusters and associations), we define different matched catalogues. Here we describe the samples used in this paper. Figure 2 illustrates examples of possible overlaps, and Table 2 provides an overview of the different samples and their sizes.

2.4.1 one-to-one sample

The ideal scenario involves a single H II region that is powered by a single stellar population, as this enables us to study trends of the escape fraction with different nebula and stellar properties. This is the case for the catalogue published in Paper I. It consists of H II regions and stellar associations that have a clear one-to-one match and contains 4169 objects in the full sample. This includes 20.5% of all H II regions and 25.6% of the associations.

For lower-mass associations, the presence or absence of individual massive stars can have a significant impact on the measured fluxes (Fouesneau & Lançon 2010; Fouesneau et al. 2012; Hannon et al. 2019). To alleviate stochastic sampling effects of the IMF, we also defined a robust one-to-one sub-sample, consisting of objects that are more massive than 104 M⊙, with ages of 8 Myr or less, and fully contained within the nebulae. This yields 474 objects and constitutes our fiducial sample, which we use as a reference throughout the paper. They are spread across all galaxies except NGC 5068 (the lowest-mass galaxy in our sample). Panel a) in Fig. 2 shows an example of such an object.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Examples of the overlap between H II regions and their ionising sources. Each cutout shows a three-colour composite image based on the five available HST bands, overlaid with the H α line emission from MUSE in red. (a) H II region with a fully contained association; (b) H II region with fully and partially contained associations plus two compact clusters; (c) H II region with two fully contained associations and one compact cluster; (d) Multiple H II regions that form a single H II region complex with multiple associations and clusters; (e) Association overlapping with two H II regions.
Table 2

Matched samples defined in this study.

2.4.2 extended sample

As reported in Paper I, a considerable number of our detected H II regions (14.7%) overlap with more than one association. We identified 2981 H II regions that contain 8871 stellar associations, but some of the associations overlap with another H II region, making it difficult to map their ionising photon flux to the correct nebula. We therefore limited our analysis to the 966 H II regions in which the associations do not overlap with another H II region (a total of 2248 associations). More than three-quarters of them contain two associations, and only a single region contains more than six associations. We refer to this sample as the extended sample, and panels b) and c) in Fig. 2 show postage-stamp cutouts for this sample.

We again defined a robust subsample; however, in this case, the definition was not straightforward. We required that at least one association be young and massive (>104 M⊙ and ≤8 Myr), which yielded 707 such objects in the initial catalogue. For the overlap, we tolerated old or light associations, but removed objects with partially overlapping young, massive associations. This yields 291 objects for our robust extended subsample. Compared to the one-to-one sample, the H II regions in the extended catalogue are brighter in H α by a factor of five.

2.4.3 complexes sample

As alluded to in Sect. 2.1, difficulties in drawing the boundaries of the H II regions can lead to large uncertainties for the measured fluxes. This is particularly problematic for the 34% of our H II regions that touch another ionised nebula.

For this sample, we grouped the bordering nebulae together and labelled them a complex. We find 9434 nebulae in 3019 complexes, but many contain ionised nebulae inconsistent with photoionisation (e.g. SNRs or PNe), indicating that other processes may contribute to the H α emission. Some complexes do not contain an ionising source, while for others, the ionising source is shared with another complex (two complexes that do not touch, but share an association that overlaps with both, as shown in panel e) of Fig. 2). Excluding these leaves 1124 complexes in our sample that contain 3722 H II regions and 3923 associations. Panel d) of Fig. 2 shows an example of such an object.

For the robust subsample, we applied the same criteria as for the robust extended sample, i.e. at least one fully contained young and massive association and no partially overlapping young and massive associations. This results in 317 objects in the robust complexes subsample.

2.4.4 cluster sample

To directly compare our results with existing studies that focus on star clusters rather than stellar associations (e.g. Della Bruna et al. 2021; Teh et al. 2023), we also constructed a sample focused on the overlap of H II regions and star clusters. We used Class 1, 2, and 3 clusters from Thilker et al. (2022), identifying 22 632 clusters that overlap with 8905 H II regions. Similar to the extended sample, multiple clusters inside a single H II region pose no issue, and this constitutes our cluster sample. We also applied the same mass and age restrictions for the robust subsample (the overlap criterion does not apply, as the clusters are treated as point sources), yielding 2223 objects in the robust cluster subsample.

3 Escape fractions

The previous section shows that there are many assumptions involved in the calculation, but the choice of the stellar population model should only slightly alter the value of Q(H0) for the robust subsamples. Here, we examine the different defined samples. We estimated the percentage of leaking radiation by combining the values from the H II regions with those from the stellar catalogues.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Comparison between the predicted ionising photon flux Q(H0) and the observed values from the H II region QHα. Full samples are shown in grey and robust subsamples in red. Contours indicate the distribution of the robust sample and diagonal lines correspond to different escape fractions.

3.1 Ionising photon fluxes for different samples

We took the samples defined in Sect. 2.4 and plotted their observed ionising photon flux from the H II regions, QHα, against the predicted fluxes from the stars, Q(H0). Figure 3 shows the results, with robust subsamples marked in red.

In the one-to-one sample, we find that for more than half of the full sample, the ionising photon production rate from the associations is insufficient to explain the observed H α flux. However, most of those objects are either low-mass or only partially overlap. In the robust one-to-one subsample, we find that 80% of the H II regions can be ionised by their associated stars. The stellar association catalogues created at different physical scales (16 pc or 64 pc; Paper I) provide similar results.

For the other samples, we find similar trends, emphasising the importance of selecting robust subsamples. It also became clear that the different catalogues probe different regimes in terms of the intensity of the ionising photon flux. Unsurprisingly, the extended sample and the complexes sample have ionising photon fluxes that are, on average, more than four times as large as those of the one-to-one sample.

3.2 Distribution of the escape fraction

Using Q(H0) from the stars and QHα from the H II regions, we applied Eq. (1) to compute the percentage of escaping ionising radiation. As already obvious from Fig. 3, for a substantial number of objects, the ionising photon flux from the stars is insufficient to ionise their matched H II region. As a result, some objects have a negative escape fraction, creating ambiguity in determining a representative escape fraction for the sample as a whole. Excluding all regions with negative fesc leads to a higher average, while including clearly unphysical values also does not make sense. Our approach is a compromise to account for objects that have a positive escape fraction within their uncertainties. We assumed that the escape fraction for each object followed a normal distribution and summed the individual distributions: p(x)=∑i=1n12πσfex,iexp((μfexc,i−x)22σfecc,i2).Mathematical equation: $p(x) = \mathop \sum \limits_{i = 1}^n {1 \over {\sqrt {2\pi } {\sigma _{{f_{ex}},i}}}}\exp \left( {{{{{\left( {{\mu _{{f_{exc}},i}} - x} \right)}^2}} \over {2\sigma _{{f_{ecc}},i}^2}}} \right).$(3)

We then calculated the median (50th percentile) and the 1σ interval (16th and 84th percentiles) within the range of valid values (0 to 100%). Figure 4 shows the results as histograms; most regions have an escape fraction greater than 50%. For robust subsamples, independent of the source of ionisation, we measured similar medians, ranging from 76 to 82%. The full sample medians are consistently within a few per cent of the robust sample values. The number of objects with negative escape fractions varies between 12 and 19%, and these are grouped together in the left bar.

3.3 Negative escape fractions

There are a considerable number of objects with a negative escape fraction, so here we take a closer look at what may cause this. The objects with a negative escape fraction in the robust one-to-one subsample are consistent with fesc = 0 within their uncertainties. Considering the contributors to the absolute uncertainty, the dominant factor is the age, followed closely by the mass, while the error on the H α flux is mostly negligible. Most objects with a negative escape fraction are therefore 4 Myr or older and have a median uncertainty of 2 Myr, making their predicted fluxes particularly uncertain. An overestimation of the age could explain an underestimation of Q(H0) and therefore resolve most of the negative escape fractions.

Another possibility is that we missed some ionising sources altogether. To explore this, we visually inspected the NUV images of the 92 objects with negative escape fraction and looked for missed peaks inside the H II regions. Only about 10% of them show peaks that were not included by the source detection of the associations.

Stars outside of the H II region might also ionise the cloud. To investigate this, we compiled a list of associations within 10 arcsec (250 to 950 pc, measured from boundary to boundary) of each H II region. We then took the flux of each of the nearby associations (excluding associations that are fully contained in another H II region) and scaled the estimated ionising photons that could be absorbed by the H II region, assuming that the flux is emitted isotropically. For each H II region, we compared the external-source ionising photon flux to overlapping-source flux. Figure 5 shows the results. We find that in the subsample with a negative escape fraction, the estimated contribution from external ionising photons is approximately a factor of 10 greater than the internal contribution, which is more than a factor 100 larger than that found in the subsample with a positive escape fraction. However, even after accounting for these external contributions, the majority of objects with a negative escape fraction remain negative. Nevertheless, this result shows that the escape fraction is underestimated if external sources are not taken into account.

So far, we have mainly focused on the stellar side, but the QHα estimated from the H II regions is also subject to uncertainties. Given that we directly measured the H α fluxes and they have a high S/N, the stated errors are significantly smaller for QHα than for Q(H0). However, they do not include uncertainties arising from choices such as where we define the region boundaries. Even a small change in the size can significantly alter the measured flux. This is especially true for closely clustered and unresolved H II regions, where we might attribute the measured flux to a neighbouring object. Increasing or decreasing the boundary by just one pixel radially (corresponding to 5 to 19 pc, depending on the distance to the galaxy) alters the measured flux on average by ±40% for the full H II regions sample, ±30% for the matched objects, and ±18% for the robust subsample. However, the flux measured from the H II regions is only roughly one-fifth of that predicted from the stars, so although the uncertainty itself is high, it typically changes the escape fraction by less than 5%.

In addition, particularly for faint H II regions, the DIG might contribute the majority of H α emission, although our brighter robust sample should only be affected at the ∼20% level (Belfiore et al. 2022; Congiu et al. 2023).

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Histograms of the observed escape fractions for the different samples. Regions with negative escape fractions are grouped together in the shaded left bar. The full sample is shown in grey, and the robust subsample and its median appear in red.
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Contribution of external ionising sources. For H II regions with negative escape fractions, the average contribution from nearby associations, not overlapping with the nebula, is far greater than for those with positive escape fractions.

3.4 Photoionisation budget of entire galaxies

As we have previously seen, it is challenging to correctly assign each H II region to its source of ionising radiation. It requires a careful selection to include all real sources while excluding false ones. Another approach is to consider all objects together to avoid problems with matching. As shown in Fig. 5, sources outside the H II region can contribute significantly to the ionising photon budget. We therefore performed this analysis on all H II regions and associations within the overlapping MUSE and HST coverage. Because we cannot exclude H II regions ionised by older or low-mass associations, this analysis must include all associations and not only the robust subsample. One potentially problematic area is the galactic centre, as some galaxies possess active nuclei. We therefore used the environment masks from Querejeta et al. (2021) to mask out the galaxy centres (eliminating 286 H II regions and 690 associations from our sample).

We summed the observed ionising photon flux across all H II regions that fall inside the HST coverage and measured QHα = 2.3 · 1054 s−1 across all galaxies. We then summed the predicted ionising photon flux from all associations within the MUSE coverage and find Q(H0) = 7.2 · 1054 s−1, yielding an escape fraction of fesc = 68.9 ± 1.3%. Similarly, summing the flux from all star clusters yields Q(H0) = 7.4 · 1054 s−1 and an escape fraction of fesc = 69.6 ± 5.0%. We list the values for individual galaxies in Table 3 and show them in Fig. 6.

With this approach, the upper end of the luminosity function should be fully sampled for both the associations and H II regions, but the sensitivity differs at the lower end. The cumulative luminosity function (see Appendix C for more details) is dominated by the brightest objects. The H II regions fainter than 1038 erg s−1 contain less than 10% of the total ionising photon flux; this is already ten times brighter than the completeness limit estimated by Santoro et al. (2022). Therefore, completeness should not be an issue for the H II regions. For the stars, the upper end of the luminosity function dominates even more.

Other minor adjustments to our approach have a negligible impact on the global result. If some real H II regions were excluded by the BPT cuts, artificially reducing the observed ionising photon flux and increasing the escape fraction, we estimate that they contribute only around ∼10% to the total H α flux in all nebulae. This would only decrease the resulting escape fraction by a few per cent. Some compact clusters are also not contained within an association. Including them in the analysis only slightly increases the predicted ionising photon flux, and the escape fraction increases by only ∼3%.

Overall, both the global and resolved approaches yield comparatively high values of 70% and above, providing confidence in our overall result for this large statistical sample of stellar sources and H II regions.

Table 3

Ionising photon budget for entire galaxies.

Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Same as Fig. 3 but with the ionising photon budget for entire galaxies. The values for the galaxies are listed in Table 3 and include all H II regions and stellar associations in the overlapping coverage. For reference, the values of individual objects from the robust one-to-one and robust extended sample are plotted as grey dots.

3.5 Trends with nebulae properties

The escape fraction is expected to change as the cloud evolves (Dayal & Ferrara 2018) and as a function of local conditions. For example, one might expect higher fesc in lower density environments or older H II regions, as the gas is already partially cleared. Figure 7 presents trends with selected properties. This analysis uses the robust one-to-one subsample. In general, the large scatter in fesc makes it challenging to identify strong correlations.

Contrary to our expectations, we find that the H II regions containing younger associations are leakier (similar to the results reported by Teh et al. 2023). As a function of the association stellar mass, the scatter is too large to identify any trends. In terms of the dust reddening, traced by the Balmer decrement, we observe an anti-correlation, consistent with the scenario that a deeply embedded cloud leaks less radiation than a more permeable one. However, when using stellar E(B − V), we do not find a strong correlation. For the ionisation parameter, log q, and the local metallicity offset, ∆(O/H), we find no correlation. Lastly, for the individual H II regions, we examined the H α luminosity and find an anti-correlation with the escape fraction, similar to Teh et al. (2023).

Based on our global measurements of fesc, we also searched for galaxy-wide trends. The escaping photons are thought to ionise the DIG and hence we expect a positive correlation. However, we see no trend between our observed fesc and the DIG fraction, fDIG, measured by Belfiore et al. (2022); our values are generally larger. Another property is the global metallicity 12 + log O/H (taken from Kreckel et al. 2020); here too, we see no dependency.

Finally, we examined the variations between different galactic environments within each galaxy, separately considering the bar, interbar, interarm, and spiral arm environments. Overall, they are very similar, with medians ranging from 75 to 83%; however, the H II regions in the bar have slightly lower escape fractions.

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Trends between different nebula properties and the escape fraction. Shown are the age and mass of the stellar associations, the reddening E(B − V)Balmer, the ionisation parameter log q, the local metallicity offset ∆(O/H) of the H II region, and the luminosity L(H α). The last row shows the global escape fraction (see Sect. 3.4) versus the fraction of DIG emission, the global escape fraction versus the average metallicity of the galaxy, and the escape fraction for different environments. The black lines represent the median values in each bin and the gray shaded areas indicate the 68 and 98 percentiles.

3.6 Comparison with existing studies

The large sample of H II regions and stellar sources used to measure the escape fraction results in a wide range of values. By examining their distribution, we see that quoting a single value is deceptive. The values published in the literature reinforce this, as they also cover a wide range depending on the conditions of the studied objects and the techniques used for the analysis.

We therefore focused our comparison on similar studies that also match observed H II regions to their stellar counterparts in nearby galaxies. Most of our escape fraction measurements exceed 60%, which is itself at the upper end of the values present in the literature (9 to 67%) for comparable studies such as McLeod et al. (2020) or Della Bruna et al. (2021). Even for independent calculations in the same galaxy (NGC 0628 by Teh et al. 2023; see Appendix B), our approach yields much larger values of fesc. Although specific modelling and measurement approaches can introduce large scatter, it is not evident that our approach systematically differs from previous studies. We aim to highlight here that these modelling uncertainties remain dominant when attempting to quantitively constrain the escape fraction, and below we emphasise several issues that stand out.

The large size of our sample prompted us to introduce certain cuts to define a robust subsample, in which we removed objects that we deemed unreliable. The stellar component must be fully contained within the H II region, its mass must be larger than 104 M⊙ to ensure full sampling of the IMF, and its age must be ≤ 8 Myr because the ionised cloud should not exist for much longer than a few million years (Hollyhead et al. 2015; Kim et al. 2021a; Hannon et al. 2022). This procedure is more objective than selecting the objects by hand, but it might still be biased towards certain physical parameters. Partly because of this selection, the regions in our sample are relatively young, with a median age of 2 Myr. Nevertheless, the age distribution appears physically reasonable, as Kim et al. (2023) find that 50% of the LyC photons are already emitted by this age. Around this age, the ionising photon flux is strongly age-dependent, so even a small uncertainty of 1 Myr can reduce Q(H0) by half. This is also reflected in the propagated error of fesc, for which it is the largest contributor.

Other issues also have uncertainties that are difficult to incorporate. As mentioned in Sect. 2.1, a small change in the size of a region can alter the flux considerably, and this is not accounted for in the stated error. Manual inspection of some of the outliers reveals that, in particular, the larger H II regions are often complex structures and difficult to segment into physically meaningful pieces. Furthermore, we did not consider external sources, which can have a major impact on small H II regions as shown in Sect. 3.3.

One aspect that we have not yet addressed is the impact of dust. This is a complex issue, and while some simulations account for the ionising photons that are directly absorbed by the dust, we can only correct the H α fluxes for reddening. According to Kim et al. (2019), a double-digit percentage may be lost in this way. This cannot explain the difference with other observations, which have generally followed the same approach as we have, but it may help explain the higher values that we observe when compared to simulations.

4 Conclusions

We combined catalogues of H II regions with catalogues of young stellar associations and star clusters to measure the escape fraction across an unprecedentedly large sample of thousands of regions across 19 nearby galaxies. This dataset provides a statistical view of how much ionising radiation escapes from the local H II region environment into the galaxy disc.

To calculate fesc, we compared QHα, derived from the extinction-corrected H α luminosity, and Q(H0), derived from the associated stellar counterpart. After comparing different methods of matching stars and gas, we focused our analysis on regions with a one-to-one match between a single H II region and a single stellar association. To mitigate the impact of stochastic sampling, we established a robust subsample of 474 massive (>104 M⊙) and young (≤8 Myr) stellar associations that are fully contained within the nebula.

  • (i)

    From the robust one-to-one subsample, we measure a comparatively high value of fesc=82−24+12%Mathematical equation: ${f_{{\rm{esc}}}} = 82_{ - 24}^{ + 12}\% $;

  • (ii)

    We derive similar escape fractions with different cross-matching methods or when using compact clusters instead of the loosely bound associations;

  • (iii)

    The escape fraction depends heavily on the derived stellar parameters and the boundaries of the H II regions;

  • (iv)

    For many H II regions, a significant fraction of the ionising photon budget comes from stars that reside outside their boundaries;

  • (v)

    If we consider all stellar associations and H II regions together, to account for the entire ionising photon budget in the galaxies, we measure slightly lower but consistent values;

  • (vi)

    We do not identify any pronounced trends between fesc and local or global galaxy properties.

This work demonstrates that a large fraction of ionising photons can escape their immediate H II region bubble and become available to power the surrounding pervasive DIG observed in local galaxies (Haffner et al. 2009; Belfiore et al. 2022).

Overall, more work is needed on the precise determination of the stellar ionising photon output before it is possible to establish more concrete connections between changes in fesc and local conditions in the ISM.

We expect that significant progress in parametrising the properties of the stellar clusters will come from the inclusion of additional infrared JWST bands and the use of priors based on the presence of H α emission (Whitmore et al. 2025). In addition, the recent availability of new HST H α high-resolution maps (Chandar et al. 2025) enables more detailed cross-matching between ionised gas and stellar populations at matched resolution for the brightest, most compact regions (Barnes et al. 2026).

These maps reveal how clustering and hierarchical structuring of young stars make simple one-to-one matching of ionised gas and stars extremely challenging.

Acknowledgements

We thank the anonymous referee for the helpful comments that improved this work. This work has been carried out as part of the PHANGS collaboration. Based on observations from the PHANGS-MUSE program, collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), 1100.B-0651 (PHANGS-MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). In addition, this research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. Support for Program number 15654 was provided through a grant from the STScI under NASA contract NAS5-26555. FS and KK gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI Kreckel) and the European Research Council’s starting grant ERC StG-101077573 “ISM-METALS”). MB acknowledges support from the ANID BASAL project FB210003. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Univer-sité Côte d’Azur under reference number ANR-15-IDEX-01. OE acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project-ID 541068876. SCOG and RSK acknowledge financial support from the European Research Council via the ERC Synergy Grant “ECO-GAL” (project ID 855130) and from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”. JEMD acknowledges support from project UNAM DGAPA-PAPIIT IG 101025, Mexico. MC and LR gratefully acknowledge funding from the DFG through an Emmy Noether Research Group (grant number CH2137/1-1). COOL Research DAO is a Decentralised Autonomous Organisation supporting research in astrophysics aimed at uncovering our cosmic origins (Chevance et al. 2025). This research made use of ASTROPY (Astropy Collaboration 2013, 2018, 2022), NUMPY (Harris et al. 2020), and MATPLOTLIB (Hunter 2007).

References

  1. Adamo, A., Ryon, J. E., Messa, M., et al. 2017, ApJ, 841, 131 [Google Scholar]
  2. Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621 [Google Scholar]
  3. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bacon, R., Accardo, M., Adjali, L., et al. 2010, Proc. SPIE, 7735, 773508 [Google Scholar]
  7. Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 [Google Scholar]
  8. Barnes, A. T., Glover, S. C. O., Kreckel, K., et al. 2021, MNRAS, 508, 5362 [NASA ADS] [CrossRef] [Google Scholar]
  9. Barnes, A. T., Chandar, R., Kreckel, K., et al. 2026, A&A, 706, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [Google Scholar]
  14. Cerviño, M., & Luridiana, V. 2006, A&A, 451, 475 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  16. Chandar, R., Barnes, A. T., Thilker, D. A., et al. 2025, AJ, 169, 150 [Google Scholar]
  17. Chevance, M., Kruijssen, J. M. D., & Longmore, S. N. 2025, arXiv e-prints [arXiv:2501.13160] [Google Scholar]
  18. Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104 [CrossRef] [Google Scholar]
  19. Congiu, E., Blanc, G. A., Belfiore, F., et al. 2023, A&A, 672, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Crowther, P. A. 2007, ARA&A, 45, 177 [Google Scholar]
  21. da Silva, R. L., Fumagalli, M., & Krumholz, M. 2012, ApJ, 745, 145 [CrossRef] [Google Scholar]
  22. da Silva, R. L., Fumagalli, M., & Krumholz, M. R. 2014, MNRAS, 444, 3275 [Google Scholar]
  23. Dale, J. E., Ercolano, B., & Bonnell, I. A. 2012, MNRAS, 424, 377 [NASA ADS] [CrossRef] [Google Scholar]
  24. Dale, J. E., Ercolano, B., & Bonnell, I. A. 2013, MNRAS, 430, 234 [Google Scholar]
  25. Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1 [Google Scholar]
  26. Della Bruna, L., Adamo, A., Lee, J. C., et al. 2021, A&A, 650, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Della Bruna, L., Adamo, A., McLeod, A. F., et al. 2022, A&A, 666, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Dolphin, A. E. 2000, PASP, 112, 1383 [Google Scholar]
  29. Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Smirnov-Pinchukov, G. V. 2018, MNRAS, 478, 3386 [NASA ADS] [CrossRef] [Google Scholar]
  31. Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]
  32. Eldridge, J. J., & Stanway, E. R. 2009, MNRAS, 400, 1019 [NASA ADS] [CrossRef] [Google Scholar]
  33. Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058 [Google Scholar]
  34. Elmegreen, B. G. 2008, ApJ, 672, 1006 [Google Scholar]
  35. Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJS, 260, 1 [NASA ADS] [CrossRef] [Google Scholar]
  37. Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022b, ApJ, 930, 126 [NASA ADS] [CrossRef] [Google Scholar]
  38. Fouesneau, M., & Lançon, A. 2010, A&A, 521, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Fouesneau, M., Lançon, A., Chandar, R., & Whitmore, B. C. 2012, ApJ, 750, 60 [NASA ADS] [CrossRef] [Google Scholar]
  40. Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Gerasimov, I. S., Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Oparin, D. V. 2022, MNRAS, 517, 4968 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gouliermis, D. A. 2018, PASP, 130, 072001 [NASA ADS] [CrossRef] [Google Scholar]
  44. Grasha, K., Calzetti, D., Adamo, A., et al. 2015, ApJ, 815, 93 [NASA ADS] [CrossRef] [Google Scholar]
  45. Groves, B., Kreckel, K., Santoro, F., et al. 2023, MNRAS, 520, 4902 [Google Scholar]
  46. Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Rev. Mod. Phys., 81, 969 [CrossRef] [Google Scholar]
  47. Hanish, D. J., Oey, M. S., Rigby, J. R., de Mello, D. F., & Lee, J. C. 2010, ApJ, 725, 2029 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hannon, S., Lee, J. C., Whitmore, B. C., et al. 2019, MNRAS, 490, 4648 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hannon, S., Lee, J. C., Whitmore, B. C., et al. 2022, MNRAS, 512, 1294 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hannon, S., Whitmore, B. C., Lee, J. C., et al. 2023, MNRAS, 526, 2991 [NASA ADS] [CrossRef] [Google Scholar]
  51. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hassani, H., Rosolowsky, E., Koch, E. W., et al. 2024, ApJS, 271, 2 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hollyhead, K., Bastian, N., Adamo, A., et al. 2015, MNRAS, 449, 1106 [NASA ADS] [CrossRef] [Google Scholar]
  54. Howard, C., Pudritz, R., & Klessen, R. 2017, ApJ, 834, 40 [Google Scholar]
  55. Howard, C. S., Pudritz, R. E., Harris, W. E., & Klessen, R. S. 2018, MNRAS, 475, 3121 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hummer, D. G., & Storey, P. J. 1987, MNRAS, 224, 801 [NASA ADS] [CrossRef] [Google Scholar]
  57. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  58. Jacobs, B. A., Rizzi, L., Tully, R. B., et al. 2009, AJ, 138, 332 [NASA ADS] [CrossRef] [Google Scholar]
  59. Japelj, J., Vanzella, E., Fontanot, F., et al. 2017, MNRAS, 468, 389 [Google Scholar]
  60. Kado-Fong, E., Kim, J.-G., Ostriker, E. C., & Kim, C.-G. 2020, ApJ, 897, 143 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kakiichi, K., & Gronke, M. 2021, ApJ, 908, 30 [CrossRef] [Google Scholar]
  62. Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055 [Google Scholar]
  63. Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961 [Google Scholar]
  64. Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2019, ApJ, 883, 102 [Google Scholar]
  65. Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2021a, MNRAS, 504, 487 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021b, ApJ, 911, 128 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kim, J.-G., Gong, M., Kim, C.-G., & Ostriker, E. C. 2023, ApJS, 264, 10 [NASA ADS] [CrossRef] [Google Scholar]
  68. Klessen, R. S., & Glover, S. C. O. 2016, Saas-Fee Adv. Course, 43, 85 [NASA ADS] [CrossRef] [Google Scholar]
  69. Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2020, MNRAS, 499, 193 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kreckel, K., Egorov, O. V., Belfiore, F., et al. 2022, A&A, 667, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  72. Krumholz, M. R., Adamo, A., Fumagalli, M., et al. 2015, ApJ, 812, 147 [NASA ADS] [CrossRef] [Google Scholar]
  73. Krumholz, M. R., McKee, C. F., & Bland -Hawthorn, J. 2019, ARA&A, 57, 227 [Google Scholar]
  74. Larson, K. L., Lee, J. C., Thilker, D. A., et al. 2023, MNRAS, 523, 6061 [NASA ADS] [CrossRef] [Google Scholar]
  75. Lecroq, M., Charlot, S., Bressan, A., et al. 2024, MNRAS, 527, 9480 [Google Scholar]
  76. Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2022, ApJS, 258, 10 [NASA ADS] [CrossRef] [Google Scholar]
  77. Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17 [NASA ADS] [CrossRef] [Google Scholar]
  78. Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3 [Google Scholar]
  79. Leitherer, C., Ortiz Otálvaro, P. A., Bresolin, F., et al. 2010, ApJS, 189, 309 [Google Scholar]
  80. Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14 [Google Scholar]
  81. Lejeune, T., Cuisinier, F., & Buser, R. 1997, A&AS, 125, 229 [Google Scholar]
  82. Leroy, A. K., Hughes, A., Liu, D., et al. 2021a, ApJS, 255, 19 [NASA ADS] [CrossRef] [Google Scholar]
  83. Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021b, ApJS, 257, 43 [NASA ADS] [CrossRef] [Google Scholar]
  84. Levesque, E. M., Leitherer, C., Ekstrom, S., Meynet, G., & Schaerer, D. 2012, ApJ, 751, 67 [NASA ADS] [CrossRef] [Google Scholar]
  85. Lopez, L. A., Krumholz, M. R., Bolatto, A. D., et al. 2014, ApJ, 795, 121 [NASA ADS] [CrossRef] [Google Scholar]
  86. Ma, X., Quataert, E., Wetzel, A., et al. 2020, MNRAS, 498, 2001 [Google Scholar]
  87. Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1 [Google Scholar]
  88. Maschmann, D., Lee, J. C., Thilker, D. A., et al. 2024, ApJS, 273, 14 [Google Scholar]
  89. McLeod, A. F., Kruijssen, J. M. D., Weisz, D. R., et al. 2020, ApJ, 891, 25 [NASA ADS] [CrossRef] [Google Scholar]
  90. Mitra, S., Ferrara, A., & Choudhury, T. R. 2013, MNRAS, 428, L1 [Google Scholar]
  91. Niederhofer, F., Hilker, M., Bastian, N., & Ercolano, B. 2016, A&A, 592, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. O’Donnell, J. E. 1994, ApJ, 422, 158 [Google Scholar]
  93. Paardekooper, J. P., Pelupessy, F. I., Altay, G., & Kruip, C. J. H. 2011, A&A, 530, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Pellegrini, E. W., Oey, M. S., Winkler, P. F., et al. 2012, ApJ, 755, 40 [NASA ADS] [CrossRef] [Google Scholar]
  95. Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 4453 [NASA ADS] [CrossRef] [Google Scholar]
  97. Ramambason, L., Schaerer, D., Stasińska, G., et al. 2020, A&A, 644, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Ramambason, L., Lebouteiller, V., Bik, A., et al. 2022, A&A, 667, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Rousseau-Nepton, L., Robert, C., Martin, R. P., Drissen, L., & Martin, T. 2018, MNRAS, 477, 4152 [NASA ADS] [Google Scholar]
  100. Rousseau-Nepton, L., Martin, R. P., Robert, C., et al. 2019, MNRAS, 489, 5530 [NASA ADS] [CrossRef] [Google Scholar]
  101. Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2022, A&A, 663, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  103. Santoro, F., Kreckel, K., Belfiore, F., et al. 2022, A&A, 658, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Scheuermann, F., Kreckel, K., Barnes, A. T., et al. 2023, MNRAS, 522, 2369 [NASA ADS] [CrossRef] [Google Scholar]
  105. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [Google Scholar]
  106. Simón-Díaz, S., & Stasińska, G. 2008, MNRAS, 389, 1009 [CrossRef] [Google Scholar]
  107. Stetson, P. B. 1987, PASP, 99, 191 [Google Scholar]
  108. Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41 [NASA ADS] [CrossRef] [Google Scholar]
  109. Tacchella, S., Smith, A., Kannan, R., et al. 2022, MNRAS, 513, 2904 [NASA ADS] [CrossRef] [Google Scholar]
  110. Teh, J. W., Grasha, K., Krumholz, M. R., et al. 2023, MNRAS, 524, 1191 [CrossRef] [Google Scholar]
  111. Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2000, AJ, 120, 3070 [Google Scholar]
  112. Thilker, D. A., Whitmore, B. C., Lee, J. C., et al. 2022, MNRAS, 509, 4094 [Google Scholar]
  113. Thilker, D. A., Lee, J. C., Whitmore, B. C., et al. 2025, ApJS, 280, 1 [Google Scholar]
  114. Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224 [Google Scholar]
  115. Turner, J. A., Dale, D. A., Lee, J. C., et al. 2021, MNRAS, 502, 1366 [NASA ADS] [CrossRef] [Google Scholar]
  116. Vázquez, G. A., & Leitherer, C. 2005, ApJ, 621, 695 [Google Scholar]
  117. Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295 [Google Scholar]
  118. Wang, B., Heckman, T. M., Amorín, R., et al. 2021, ApJ, 916, 3 [NASA ADS] [CrossRef] [Google Scholar]
  119. Wei, W., Huerta, E. A., Whitmore, B. C., et al. 2020, MNRAS, 493, 3178 [NASA ADS] [CrossRef] [Google Scholar]
  120. Whitmore, B. C., Chandar, R., Lee, J. C., et al. 2025, ApJ, 982, 50 [Google Scholar]
  121. Zastrow, J., Oey, M. S., & Pellegrini, E. W. 2013, ApJ, 769, 94 [NASA ADS] [CrossRef] [Google Scholar]
  122. Zurita, A., Rozas, M., & Beckman, J. E. 2000, A&A, 363, 9 [Google Scholar]

Appendix A Ionising photon flux from different models

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Ionising photon flux Q(H0) of a star cluster with a mass of 106 M⊙ as a function of age. The top panel shows the values predicted by different models, all at solar metallicity, while the lower panel shows the SB99 v00 model at various metallicities.

In Fig. A.1 we compare the ionising photon flux Q(H0) for different models (top panel) and metallicities (lower panel). A detailed listing of the models can be found in Table A.1. These plots demonstrate that, at constant metallicity, all models except the one with rotation are similar for the first few million years. Around 6 Myr, however, the effects of binarity becomes relevant as its flux decreases much more slowly than the other models. For the impact of metallicity, we include the low-metallicity models from STARBURST99, based on the stellar model by Georgy et al. (2013). The variation between models that are close to solar metallicity is practically negligible, while the predicted flux is slightly higher at lower metallicities. Given that the flux decreases by an order of magnitude every few million years and the uncertainty of the age is comparatively high, the choice of model is likely to have only a minor influence.

Appendix B Comparison with Teh et al. (2023)

The stellar catalogues presented in Sect. 2.2 rely on the assumption that the clusters and associations in our sample are massive enough such that a deterministic relationship between the physical and photometric properties of the stellar populations exists. While we try to ensure this by focusing our analysis on high mass objects from the robust subsamples defined in Sect. 2.4, this means that we are omitting a large part of the sample. Another approach is to employ a stochastic stellar population synthesis code like SLUG (da Silva et al. 2012, 2014) to infer a distribution of possible values for Q(H0).

Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Comparisons between the H II regions from PHANGS (Groves et al. 2023, observed with MUSE) and the ones from SIGNALS (Rousseau-Nepton et al. 2019, observed with SITELLE).

A particularly relevant literature result is the study by Teh et al. (2023), who measured the escape fraction of NGC 0628 (which is also in our sample) with the help of SLUG. They combined H II regions from the Star formation, ionised gas and nebular abundances legacy (SIGNALS) survey (Rousseau-Nepton et al. 2018) with star clusters from LEGUS (Grasha et al. 2015; Adamo et al. 2017). A large fraction of the 139 objects in their matched catalogue have negative escape fractions and based on a subsample of 42 objects, they measured a comparatively low value of fesc = 9 ± 6%. The striking difference to our results asks for a closer comparison.

We start by comparing the H II regions from MUSE to those from SITELLE and 643 of them match within 0.8 arcsec. In Fig. B.1 we compare the measured fluxes, and while the ones identified by SIGNALS are more luminous at the fainter end, the ones from PHANGS are brighter for the most luminous objects.

To compare the stellar component, we take the five-band filters from our association catalogue and redo the SED fit with SLUG, following the procedure described in Teh et al. (2023): First, a library of 107 stellar populations is generated to account for the effects of stochasticity in the IMF. For this, we adopt the following physical parameters: Solar metallicity, the IMF by Kroupa (2001), Geneva evolutionary tracks (Ekström et al. 2012; Georgy et al. 2012), STARBURST99 stellar atmosphere models (Leitherer et al. 1999; Vázquez & Leitherer 2005), and a Milky Way extinction curve (Cardelli et al. 1989). Then we use a modified version of cluster_slug (Della Bruna et al. 2021, 2022; Teh et al. 2023) to compute the posterior probability density function (PDF) of the age, mass, and ionising luminosity (see Krumholz et al. 2015, for more details). We represent the Q(H0) values with the median (50th percentile) and the 1σ uncertainty with the 16th–84th percentile of the PDF distribution.

The result is shown in Fig. B.2. As far as the masses are concerned, there is a large scatter, but not systematic differences. There are significant deviations in the derived ages, with those from SLUG generally being older. This is particularly problematic in the case of very young stars. While the lower limit of the age distribution is often consistent with being 3 Myr or younger, the median rarely is. Since the majority of the ionising photons are emitted during this period, the predicted fluxes are generally lower.

Table A.1

Codes used to compute the ionising photon flux Q(H0).

When matched to the H II region catalogue from MUSE, we can compute the resulting escape fractions as shown in Fig. B.3. Only 776 of them (∼19%) have positive escape fractions, and their escape fractions are uniformly distributed across the entire range from 0 to 100%. If we focus on the young, massive, and contained objects that form a robust subsample, we find 228 objects that fulfil our criteria and only 50 of them (∼22%) have positive escape fractions. It was to be expected that SLUG would deliver different values at the faint end, due to stochastic sampling. Overall, it predicts systematically lower Q(H0), even at the bright end, where stochastic sampling should not be an issue. The lower Q(H0) is responsible for significantly lower escape fractions, which are mostly uniformly distributed between 0 to 45%, and are in better agreement with results found by previous studies that are based on SLUG (Della Bruna et al. 2022; Teh et al. 2023).

Appendix C Completeness issues

In Sect. 3.4 we consider the entire ionising photon budget of each galaxy. In this scenario, different completeness limits for the H II regions and stellar populations could bias those results. We therefore estimate the completeness for both to assess the contribution of missed objects. In Fig. C.1 we show the normalised cumulative ionising photon flux to estimate the contribution of faint H II regions and low-mass associations.

In the case of IC 5332, NGC 3351 and NGC 5068, faint H II regions (below the completeness limit from Santoro et al. 2022) contribute more than 10% of the total ionising photon flux. The potentially missing H α flux could increase the escape fraction, but while the values of the first two galaxies are above the average (see Table 3), both also appear on the stellar side, weakening this conclusion.

For IC 5332, NGC 1300, NGC 1433, NGC 1512, NGC 2835 and NGC 3351, the contribution of low-mass associations (< 5000 M⊙) is larger than 10%. In contrast to the case of the H II regions, the potentially overlooked clusters would lead to a decrease of the escape fraction. While half of galaxies have lower escape fractions, the other half, including the two galaxies mentioned above, have higher values.

For the 12 other galaxies, the completeness limit should not be an issue, and even for the seven discussed above, there is no noticeable impact.

Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

We ran SLUG on the stellar associations in the matched catalogue. For the comparison we define a robust subsample, where either CIGALE or SLUG predict an age younger than 8 Myr and a mass above 104 M⊙ (859 objects). The masses show a wide scatter in both directions, but SLUG generally predicts older ages. As a result, the ionising photon flux is systematically lower than the one from CIGALE.
Thumbnail: Fig. B.3 Refer to the following caption and surrounding text. Fig. B.3

Histograms of the escape fractions based on stellar parameters inferred from SLUG. The full sample is shown in gray and the robust subsample of young and massive, fully contained associations is shown in red. All regions with negative escape fractions are grouped together in the shaded left bar.
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Normalised cumulative Q(H0) as a function of Q(H0) for H II regions and stellar associations. The individual galaxies are shown in pastel and the sum of all galaxies is shown in bold colours. We only show the H II regions and associations in the overlapping coverage that do not fall in the centre (corresponding to the sample used in Sect. 3.4). The H II region completeness limit from Santoro et al. (2022), corresponding to a luminosity of 1037 erg s−1 is indicated by the left red dashed line. The two blue lines corresponds to a 103 M⊙ and 104 M⊙ cluster that is 1 Myr old.

All Tables

Table 1

Properties of the PHANGS galaxy sample.

In the text
Table 2

Matched samples defined in this study.

In the text
Table 3

Ionising photon budget for entire galaxies.

In the text
Table A.1

Codes used to compute the ionising photon flux Q(H0).

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Comparison of different population synthesis models. We use a sample with random ages and masses and compute the ionising photon flux Q(H0) with BC03, SB99, and BPASS. The grey points represent the full sample, while the red points indicate a young and massive sub-sample (age ≤ 8 Myr and mass > 104 M⊙). The Spearman correlation coefficient ρ indicates good agreement between the different models across the robust subsamples, and all correlations are statistically significant (p-values < 0.05).
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Examples of the overlap between H II regions and their ionising sources. Each cutout shows a three-colour composite image based on the five available HST bands, overlaid with the H α line emission from MUSE in red. (a) H II region with a fully contained association; (b) H II region with fully and partially contained associations plus two compact clusters; (c) H II region with two fully contained associations and one compact cluster; (d) Multiple H II regions that form a single H II region complex with multiple associations and clusters; (e) Association overlapping with two H II regions.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Comparison between the predicted ionising photon flux Q(H0) and the observed values from the H II region QHα. Full samples are shown in grey and robust subsamples in red. Contours indicate the distribution of the robust sample and diagonal lines correspond to different escape fractions.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Histograms of the observed escape fractions for the different samples. Regions with negative escape fractions are grouped together in the shaded left bar. The full sample is shown in grey, and the robust subsample and its median appear in red.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Contribution of external ionising sources. For H II regions with negative escape fractions, the average contribution from nearby associations, not overlapping with the nebula, is far greater than for those with positive escape fractions.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Same as Fig. 3 but with the ionising photon budget for entire galaxies. The values for the galaxies are listed in Table 3 and include all H II regions and stellar associations in the overlapping coverage. For reference, the values of individual objects from the robust one-to-one and robust extended sample are plotted as grey dots.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Trends between different nebula properties and the escape fraction. Shown are the age and mass of the stellar associations, the reddening E(B − V)Balmer, the ionisation parameter log q, the local metallicity offset ∆(O/H) of the H II region, and the luminosity L(H α). The last row shows the global escape fraction (see Sect. 3.4) versus the fraction of DIG emission, the global escape fraction versus the average metallicity of the galaxy, and the escape fraction for different environments. The black lines represent the median values in each bin and the gray shaded areas indicate the 68 and 98 percentiles.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Ionising photon flux Q(H0) of a star cluster with a mass of 106 M⊙ as a function of age. The top panel shows the values predicted by different models, all at solar metallicity, while the lower panel shows the SB99 v00 model at various metallicities.
In the text
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Comparisons between the H II regions from PHANGS (Groves et al. 2023, observed with MUSE) and the ones from SIGNALS (Rousseau-Nepton et al. 2019, observed with SITELLE).
In the text
Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

We ran SLUG on the stellar associations in the matched catalogue. For the comparison we define a robust subsample, where either CIGALE or SLUG predict an age younger than 8 Myr and a mass above 104 M⊙ (859 objects). The masses show a wide scatter in both directions, but SLUG generally predicts older ages. As a result, the ionising photon flux is systematically lower than the one from CIGALE.
In the text
Thumbnail: Fig. B.3 Refer to the following caption and surrounding text. Fig. B.3

Histograms of the escape fractions based on stellar parameters inferred from SLUG. The full sample is shown in gray and the robust subsample of young and massive, fully contained associations is shown in red. All regions with negative escape fractions are grouped together in the shaded left bar.
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Normalised cumulative Q(H0) as a function of Q(H0) for H II regions and stellar associations. The individual galaxies are shown in pastel and the sum of all galaxies is shown in bold colours. We only show the H II regions and associations in the overlapping coverage that do not fall in the centre (corresponding to the sample used in Sect. 3.4). The H II region completeness limit from Santoro et al. (2022), corresponding to a luminosity of 1037 erg s−1 is indicated by the left red dashed line. The two blue lines corresponds to a 103 M⊙ and 104 M⊙ cluster that is 1 Myr old.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.