© ESO 2021

1. Introduction

As a result of hierarchically-structured, accreting molecular clouds, stars preferentially form in clustered environments. More than 70% of massive stars are located in clusters and OB associations (Lada & Lada 2003), making young star clusters (YSCs) preferential sites from which radiative and mechanical feedback originates. Typically, only a minority of these stellar aggregates are gravitationally bound and will, therefore, survive as star clusters (Kruijssen 2012; Grudić et al. 2020). These are the clusters that sit in the densest regions of the hierarchy and are therefore able to survive the removal of gas long enough to reach a high star formation efficiency and become dynamically relaxed and well-mixed (Kruijssen et al. 2019). The remaining loosely-bound clusters will disperse in short timescales and blend into the field stellar population.

A detailed understanding of the process of stellar feedback is important as this mechanism leads to the suppression of global star formation in galaxies and to a multiphase and metal-enriched interstellar medium (ISM). In particular, it will shed light on the origin of the diffuse ionised gas (DIG) that has been observed to constitute up to 50% of the total Hα luminosity in local spiral galaxies (Ferguson et al. 1996; Hoopes et al. 1996; Zurita et al. 2000; Thilker et al. 2002; Oey et al. 2007). The origin of this gas is unclear (see e.g. Mathis 2000; Haffner et al. 2009, for a review), and it has been linked to various processes including leaking H II regions (Zurita et al. 2002; Weilbacher et al. 2018), evolved field stars (Hoopes & Walterbos 2000; Zhang et al. 2017), shocks (Collins & Rand 2001), and cosmic rays (Vandenbroucke et al. 2018). Both theory and observation have shown that stellar feedback can facilitate the escape of ionising radiation from the star forming regions (see e.g. Bik et al. 2015; Menacho et al. 2019) if the youngest star clusters have had time to clear out channels through the dense leftover gas from the parent giant molecular clouds (GMCs) from which they were born (Dale 2015; Howard et al. 2018). In order to study how the radiation from massive stars affects the ISM, it is therefore necessary to investigate the timescale for photons to find a channel through which they can leave the region. Various studies have found that small scale physics plays an important role in the escape of ionising radiation. In particular, the effect of temporal and spatial clustering of supernovae (SN) can have important consequences (Gentry et al. 2017; Kim et al. 2017; Fielding et al. 2018), resulting in a transferred momentum up to an order of magnitude larger than for isolated SN. Simulations of individual galaxies do indeed indicate that the escape of ionising radiation is largely driven by the distribution of dense gas (see e.g. Paardekooper et al. 2011), and high-resolution GMC simulations have proven that photons are more likely to escape from less massive GMCs (Howard et al. 2018). Consequently, late-emitted ionising photons are more likely to escape the star forming region as the gas surrounding the region is dispersed by ionising radiation and supernova feedback (Kim et al. 2013). The same trend has been recently observed in high-resolution cosmological simulations, showing that photons are preferentially leaking from star forming regions containing feedback-driven superbubbles (Ma et al. 2020). Finally, stars are also responsible of yielding metals that will enrich the ISM, and a detailed understanding of abundance variations in star forming regions is needed in order to constrain chemical evolution models. However, until recently, it has been unclear whether star forming galaxies feature, in general, a homogeneous metal enrichment or whether the enrichment varies in different star forming regions and within the regions themselves. Recent high-resolution studies (e.g. James et al. 2016; McLeod et al. 2019) seem to indicate the latter, implying that feedback mechanisms might play an important role in local metallicity variations.

High resolution studies of star forming regions are therefore an essential piece of information in order to bridge the spatial scales that critically connect the sources of feedback (stars and clusters, i.e. a few parsec scales) and their immediate surroundings (10 s parsec scale) to galactic scale dynamics (kpc scale). Such studies, initially limited to narrow-band photometry (e.g. Pellegrini et al. 2012), were revolutionised with the advent of instruments with high spatial resolution integral field spectroscopy capabilities, such as the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) instrument at ESO Very Large Telescope (VLT). Recently, McLeod et al. (2019, 2020) exploited the MUSE capability for a detailed study of H II regions in the Large Magellanic Cloud and in the dwarf spiral NGC 300, mapping the ionisation structure of the regions and constraining their escape of ionising radiation. In this work, we perform a similar study in the local spiral galaxy NGC 7793 at a spatial resolution of ∼10 pc with MUSE. In Della Bruna et al. (2020, henceforth Paper I), we presented the data, constructed a sample of H II regions, determined the fraction of DIG in our field of view (FoV) and studied the properties of the ionised gas in emission line diagrams. Here, we combine the MUSE results on the ionised gas with the stellar and cluster population studied with the Hubble Space Telescope (HST) and the molecular gas mapped with the Atacama Large Millimeter/submillimeter Array (ALMA), construct a photoionisation budget for each H II region and study how it relates to its stellar and GMC content and to its ionisation structure.

This work is structured as follows. In Sect. 2 we give a brief overview of the data. In Sect. 3 we present the stellar census in our FoV and analyse the distribution of YSC properties. In Sect. 4 we compare extinction maps derived from individual stellar extinction and from the ionised gas and in Sect. 5 we derive the oxygen abundance of the H II regions. In Sect. 6 we analyse the ionisation structure of the gas, and in Sect. 7 we compile a photoionisation budget for the entire FoV and in more detail for a subset of the H II regions. Finally, in Sect. 8 we discuss our findings and in Sect. 9 we give a brief summary and present our conclusions.

2. Data description

The data and data reduction steps are described in detail in Paper I: we give here a brief summary. We observed the nearby flocculent spiral galaxy NGC 7793 with MUSE¹. The dataset consists of two pointings in the Wide Field Mode (WFM) adaptive optics assisted (AO) configuration, with the extended wavelength setting. On the final datacube, we measure a PSF of 0.62″ FWHM (at 7000 Å), corresponding to a spatial resolution of ∼10 pc at the distance of NGC 7793.

NGC 7793 has also been observed in the framework of the HST Treasury programme LEGUS² (Calzetti et al. 2015) and in ALMA CO(J = 2 − 1). The catalogues of YSCs derived from the LEGUS data are publicly available³: in this work, we make use of the catalogue generated with the Geneva stellar evolution models and Milky Way extinction. We are moreover selecting only clusters labelled as ‘Class 1’ (compact and symmetric), ‘Class 2’ (concentrated but with some degree of asymmetry), and ‘Class 3’ (multiple peak systems). For more details on the cluster catalogues (e.g. cluster identification, classification and photometry) see Adamo et al. (2017). Lee et al. (in prep.) have identified O stars based on HST photometry. The catalogue is described in Wofford et al. (2020), and consists of main sequence O stars of M > 20 M_⊙, selected in a colour-magnitude and colour-Q diagram⁴. The resulting completeness of this stellar catalogue is estimated to be about 75% at the lower mass limit of 20 M_⊙ and rapidly growing to 90% at stellar masses above 25 M_⊙⁵. The authors report also an average contamination rate ∼30%. We use the latter rate to account for uncertainties in the contribution of O stars to the ionisation budget.

The catalogue of GMCs built on the ALMA dataset is presented in Grasha et al. (2018). The ALMA data have an angular resolution of 0.85″ and a velocity resolution of 1.2 km s⁻¹. The data are sensitive to emission up to 11″ and are mapping GMC down to masses ∼10⁴ M_⊙.

3. Properties of the stellar population

3.1. Stellar census

We present here the full census of the stellar population in our FoV. Figure 1 shows a map of the Hα emission line obtained from the integration of the MUSE stellar continuum subtracted datacube in the rest-frame wavelength range 6559 − 6568 Å. Overlaid on the Hα emission, we show the position of candidate O stars from Lee et al. (in prep., white circles), YSCs from LEGUS (green triangles) and of Wolf Rayet (WR) stars that we spectroscopically identified in Paper I (listed in Table 1; blue crosses). We also show the contours of H II regions identified in Paper I (in black, corresponding to a flux brightness cut in Hα = 6.7 × 10⁻¹⁸ erg s⁻¹ cm⁻² spaxel⁻¹) and of GMCs (light blue filled contours). The total stellar census in our FoV consists of 651 O stars, 93 YSC, and 9 WR stars.

Fig. 1. Stellar census in the MUSE FoV. The location of candidate O stars (white circles), YSCs (green triangles) and WR stars (blue crosses) is indicated, overlaid on a map of the Hα emission. Black contours indicate the H II regions identified in Paper I, light blue filled contours denote the position of GMCs.

From Fig. 1, we observe that both YSCs and candidate O stars are abundant over the entire FoV. However, we note that several of the smaller and medium-bright H II regions do not seem to host any cluster. Overall, we find ∼49% field⁶ clusters (35% if considering only clusters with median age < 10 Myr, as motivated in Sect. 7.5) and ∼38% field candidate O stars. We note that the numbers of O-stars quoted here are not corrected for completeness, so that the latter fraction could in part be due to the high contamination rate of catalogue (see Sect. 2). We observe that WR stars are populating the H II regions, with the exception of one object (located in the upper left corner and catalogued as WR 3 in Table 1). Rate et al. (2020) recently found that isolated WR stars are not an exception, and amount to at least 64% in the Milky Way (as observed in Gaia Collaboration 2018). The authors explore the possible origin of field WR stars in simulations, finding that the most frequent mechanisms are the formation in a low density association that then expands during the star lifetime and the ejection from a cluster through dynamical interaction or binary disruption. This has been observed in galactic star clusters (e.g. Westerlund 2 and NGC 3603 Drew et al. 2018, 2019). The fact that WR 3 is located ∼100 pc away from a nearby HII-region cluster (see Fig. 1) could support the latter formation scenario. In order to investigate this hypothesis, we have computed the minimum velocity needed for the star to be ejected from the nearby cluster. We assume an age range of (4 ± 1) Myr for the WR star⁷. Assuming that the star has been ejected immediately after birth, we find ${\mathit{v}}_{\mathrm{proj},\min }={25}_{-8}^{+5}$ km s⁻¹, which meets the standard runaway criterium of v > 25 km s⁻¹ (e.g. Drew et al. 2018; Portegies Zwart 2000).

3.2. YSCs age, mass, and ionising flux distribution

In this subsection, we investigate the physical properties of the YSCs, namely their age, mass, and ionising photon flux. The latter, Q(H⁰) [s⁻¹], is a measure of the number of Lyman continuum (LyC) hydrogen-ionising photons (hν > 13.6 eV) emitted per unit time.

For this purpose, we use the stochastic stellar population synthesis code SLUG (da Silva et al. 2012; Krumholz et al. 2015a, v2). SLUG implements a Monte Carlo technique by randomly drawing stars from an input initial mass function (IMF). This approach allows one to account for the effect of stochastic sampling of the IMF, which is especially important when populating low mass clusters (see e.g. Hannon et al. 2019). This is of a particular relevance in our target: NGC 7793 has a very low SFR ∼ 0.1 M_⊙ yr⁻¹ (Calzetti et al. 2015), resulting in H II regions hosting low mass (∼a few 1000 s M_⊙) clusters. The SLUG software package includes the cluster_slug tool (Krumholz et al. 2015a), that implements Bayesian inference techniques to calculate posterior probability distribution functions (PDFs) of physical parameters of star clusters based on observed cluster photometry, as described in Krumholz et al. (2015b). This is achieved by comparing the data to a large library of simulated clusters. We modified cluster_slug to also return the ionising photon luminosity Q(H⁰) as a further output parameter and ran the code on HST broadband photometry of the YSCs. We consider age t, mass M and visual extinction A_V as free parameters, and assume a flat prior in A_V and in log t, and a log(M)∼1/M prior on the mass. We use the library of mock star clusters described in Ashworth et al. (2018), a Milky Way extinction law by Fitzpatrick (1999) and the non-rotating solar metallicity stellar population models from Ekström et al. (2012). The abundance map derived in Sect. 5 confirms that NGC 7793 has a metallicity close to solar (see also Pilyugin et al. 2014); the same metallicity was also assumed in the SLUG models of Krumholz et al. (2015b). We discuss the impact of including binary stars and the effect of stellar rotation on the resulting physical properties in Sect. 8.

Figure 2 shows the spatial distribution across the FoV of YSC ages (left panel) and masses (right panel). We are showing here as best estimate the median of each cluster PDF. As remarked in Krumholz et al. (2015b), reducing a full PDF to a single point inevitably leads to some imprecision, especially if the PDF is multi-peaked. We follow the method tested in Krumholz et al. (2015b) and we refer the interested reader to this work for evaluations of the methodology and effects on the recovered cluster physical parameters. In Fig. 3, we show the age, mass, and ionising flux PDF for YSCs located inside H II regions (in blue) and for field clusters (in orange). Thin and thick lines in the plot indicate, respectively, the PDFs of single clusters and the combined PDFs obtained by summing the fractional probability of each cluster in every age, mass, and flux bin, and re-normalising the resulting distribution. In the case of the age, both the single and combined PDFs are (most often) multi-peaked; the main peak of the combined PDF is located at t ∼ 5 Myr, both for field and H II region clusters. However, we see that the PDF of field clusters flattens towards older ages, while the probability distribution of H II region clusters is more peaked at ages below 10 Myr, that is clusters located in H II regions have a higher probability to be young and potentially ionising sources. The mass distribution is single-peaked for both categories of clusters, with a tail towards higher masses. The peak of the distribution is located at log(M) = 3.0 and 3.4 M_⊙ for H II region clusters and field clusters respectively. The Q(H⁰) distribution is, as the age, multi-peaked (both the single and combined PDFs) with two main modes (both for H II region and field clusters) at log(Q(H⁰)) = 44.1 and 48.8 s⁻¹. Nevertheless, we notice that clusters in H II regions have a higher probability to be associated with the high Q(H⁰) value, while the opposite is true for field clusters. In general however, the inferred masses and Q(H⁰) are on the low range for YSCs: for a reference, the peak at log(Q(H⁰)) = 48.8 s⁻¹ corresponds to the emission of a single O7.5 star (Martins et al. 2005). We comment more on the mass distribution below (Fig. 4 and associated text).

Fig. 2. Spatial distribution of YSC ages (left panel) and masses (right panel) in the MUSE FoV, overlaid on the Hα map shown in Fig. 1. Circles and diamonds indicate, respectively, clusters younger and older than 10 Myr. Ages and masses are best estimates computed from the median of the corresponding PDF (see Sect. 3.2). Black contours indicate the H II regions identified in Paper I.

Fig. 3. Posterior probability distributions of YSC age (left) mass (centre) and ionising flux (right) for clusters located in H II regions (in blue) and for field clusters (in orange). Thick lines indicate the combined PDFs, thin lines correspond to the PDFs of single clusters. The grey dashed lines indicate the location of the main modes of the combined distributions.

Fig. 4. Distribution of YSC mass (clusters with t < 10 Myr, in blue), GMC mass (in green) and H II regions luminosity (in black, reddening corrected) in our FoV. YSC masses are best estimates computed from the median of the corresponding PDF (see Sect. 3.2). The histograms are normalised to a peak value of 1.

In synthesis, we observe that YSCs in H II regions have a higher probability to be younger, slightly lower mass, but a more significant source of ionising photons than clusters in the field. These trend can easily be visualised in the plots of Fig. 2, where we colour code the positions of the clusters in the FoV according to their median ages and masses.

The age range observed for the clusters within the H II regions is in agreement with studies of YSCs in local galaxies (Whitmore et al. 2011; Hollyhead et al. 2015; Hannon et al. 2019; Grasha et al. 2019), that find that clusters emerge from their natal gas in < 4 − 5 Myr, but that the process can start as early as 2 − 3 Myr. Also simulations of cluster formation including solely the effect of photoionisation and stellar winds (Dale et al. 2014; Bending et al. 2020), have shown that within 3 Myr, and before the onset of SN, several clusters have vacated some channels with low gas density. These channels are easily ionised; therefore ionising radiation can escape from the H II regions even if star formation is still taking place. In Fig. 1 we observe that a significant amount of molecular gas is still detectable in the H II regions hosting YSCs, in agreement with the findings of Dale et al. (2014).

Finally, in Fig. 4 we show the total demographics of GMCs (in green), YSCs with age < 10 Myr (in blue) and H II regions (in black) within our FoV, highlighting the various phases of the star formation cycle. We observe that GMCs, hosting the molecular gas that will give rise to the next generation of stars, span a rather low range of masses compared to what observed in other local galaxies (log M_GMC ∼ 4.5 − 8 M_⊙, Hughes et al. 2016).

The collapse of the clouds results in turn in low-mass YSCs. While there is not a one-to-one correlation between the mass spectrum of GMCs and YSCs, we observe that overall the two distributions would suggest an integrated SFE of a few percent of what is typically observed in local galaxies (see e.g. the review by Krumholz et al. 2019). The bulk of the YSCs have masses around 1000 M_⊙, and the distribution reaches up to ∼10^3.5 M_⊙. These mass distributions are typical of galaxies with low SFR, such as M 31 (Larsen 2009; Johnson et al. 2017; Adamo et al. 2020). If we consider that the ionising feedback from a cluster of 1000 M_⊙ gives rise to an Hα luminosity⁸ of log(L)∼37.3 erg s⁻¹, then the observed luminosity distribution is consistent with the mass distributions of the YSCs in the FoV. If compared to studies of H II regions luminosity functions in the local universe, our sample occupies the low-L end of the luminosity function of local spirals (log(L)∼37 − 40 erg s⁻¹, Kennicutt et al. 1989). Overall, the low masses and luminosity observed are in agreement with the flocculent morphology and low SFR of our target.

4. Comparison between stellar and ionised gas extinction with maps of cold dense gas and hot dust emission

The HST, MUSE, ALMA, and Spitzer/MIPS⁹ coverage gives us the rare opportunity of comparing extinctions as traced by stellar reddening and the ionised gas phase to the distribution of the cold dense gas and the hot dust emission. In Paper I, we have derived an extinction map from the MUSE Hα/Hβ ratio using PYNEB (Luridiana et al. 2015). Hereby, we assumed a theoretical Hα/Hβ_int = 2.863 and we obtained the observed ratio by binning the stellar continuum subtracted cube with the Voronoi technique of Cappellari & Copin (2003) to a S/N = 20 in Hβ and fitting both lines with a single Gaussian profile. Here, we compare the gas extinction map derived in Paper I with the map derived by Kahre et al. (2018) from the individual stellar extinctions obtained from HST, and correlate the two maps with the position of GMCs traced by ALMA. Figure 5 (left panel) shows overlaid on the gas extinction map the contours of the weighted average map with adaptive resolution from Kahre et al. (2018). The latter was derived by spatially binning the HST data in the smallest possible regions (of varying size between 1 and 10 arcsec²) hosting at least 10 O stars. In order to compare the two maps, we have converted the colour excess E(V − I) to E(B − V) assuming a Milky Way extinction law with R_V = 3.1, as in Kahre et al. (2018). We observe a clear correspondence between the two maps, despite the difference in binning and in the method used to derive the extinction. In general, regions in which we measure a high extinction also feature a comparable stellar extinction, although the position of the extinction peaks is slightly shifted in the two maps. Conversely, the ionised gas extinction seems to miss some of the regions of enhanced stellar extinction. We notice that these regions lie outside the bright ionised gas, where Hβ is poorly detected and the ionised gas might not be a sensitive enough tracer of extinction. However, we observe that the gas extinction seems to better correlate with the position of the GMCs (orange contours in Fig. 5). We also note that there seems to be a better agreement between gas extinction peaks and GMC position within H II regions.

Fig. 5. Extinction map derived from the ionised gas observed with MUSE assuming an intrinsic Hα/Hβ ratio (from Paper I) compared to: Left panel: contours of the extinction map derived by Kahre et al. (2018) from the individual stellar extinctions obtained from HST data. The contours correspond to values of E(B − V) = 0.2 − 0.5, in steps of 0.033. Right panel: 24 μm emission from Spitzer/MIPS. The contours correspond to a flux of [1.5, 2, 3, 4, 6, 9, 11] MJy sr−1. The thick orange contours indicate the position of GMCs.

In the right panel of Fig. 5, we furthermore show overlaid on the gas extinction map contours of 24 μm emission from Spitzer/MIPS. The 24 μm emission is tracing hot dust in the galaxy, and inside the H II regions it is indicative of where star formation is taking place. We observe a good correlation between gas extinction, GMC position and 24 μm emission, indicative of the fact that, in HII regions, the dust is well mixed with gas. We comment more on this fact in relation with the effect of removal of LyC photons by dust prior to their absorption by neutral hydrogen in Sect. 8.

5. Mapping the oxygen abundance within the H II regions

We estimate an oxygen abundance in the H II regions using the strong-line method of Pilyugin & Grebel (2016). In Paper I, we derived electron temperature and density from the [S III]6312/9069 and [S II]6716/6731 line ratios¹⁰; however, due to the weakness of the [S III]6312 line and to the fact that the [S II] ratio is largely insensitive to n_e < 30 cm⁻³, we were only able to obtain a coarse T_e map and marginally constrained n_e, preventing us from measuring direct abundances. The Pilyugin & Grebel (2016) method bases on ratios of strong lines calibrated on a sample of ∼300 H II regions. Here we make use of the S calibration, based on the following three line ratios:

$$\begin{array}{cc}& N_2=(I_{[\text{N}\text{II}]}\lambda 6548+I_{[\text{N}\text{II}]}\lambda 6584)/I_{\mathrm{H}\beta },\hfill \\ \hfill & S_2=(I_{[\text{S}\text{II}]}\lambda 6716+I_{[\text{S}\text{II}]}\lambda 6731)/I_{\mathrm{H}\beta },\hfill \\ \hfill & R_3=(I_{[\text{O}\text{III}]}\lambda 4959+I_{[\text{O}\text{III}]}\lambda 5007)/I_{\mathrm{H}\beta }.\hfill \end{array}$$

We use the upper branch of the calibration, valid for log N₂ ≳ −0.6, as our data lie well above this limit throughout the FoV. The line ratios are determined from the reddening corrected fluxes, which have been obtained by binning the data with the Voronoi technique to a S/N = 20 in Hβ as described in Sect. 4. Figure 6 shows the resulting metallicity variation in a subset of the H II region sample, for which we construct a detailed photoionisation budget in Sect. 7. The exact location of these regions in the FoV is shown in Fig. 8. We observe that the abundance ranges from 12 + log(O/H) ∼ 8.25 to 8.50, with a median value of 8.37. This range of values is in agreement with the strong-line estimate of Pilyugin et al. (2014) for NGC 7793 and with the metallicity range inferred in Paper I when comparing our data with the models of star forming galaxies from Levesque et al. (2010) in ‘BPT’ diagrams (Baldwin et al. 1981).

Fig. 6. Oxygen abundance for the subset of H II regions labelled in Fig. 8, determined with the S-strong line calibration of Pilyugin & Grebel (2016). The location of YSCs is indicated with filled circles, colour-coded based on their age (best estimates derived from the median of each PDF). Red crosses and purple stars correspond to the position of WR stars and PNe, respectively. Dark red contours indicate a ratio of [S II]/[O III] = 0.5 (see Fig. 8).

In small and isolated H II regions, such as regions 4, 5, 6, and 7, we notice that the oxygen abundance varies of about 0.1 dex at most. The largest H II region complexes, such as regions 1, 2, 3, and 8 show larger gradients and pockets of enriched gas coincident or arranged in shell-like structures around clusters and massive stars¹¹. This is particularly visible in region 4, where the YSC is surrounded by an enriched shell, and in regions 6 and 8 (centre left and bottom right) where we estimate a local enhancement in the oxygen abundance at the location of YSCs. We furthermore observe that some of the WR stars are associated with enriched shells or pocket of ionised gas: this is the case in region 1 (WR2) and 8 (WR4), likely due to nitrogen-enriched gas expelled by their stellar winds. In other cases, we do not observe any enrichment surrounding the stars: this is especially visible in region 3 and in the north-west, and south of region 1. Finally, we observe a large enhancement of the abundance at the location of one of the planetary nebulae identified in Paper I (PN2).

We would like to point out that, although such indirect methods are generally calibrated for unresolved H II regions, they have recently been applied also to study abundance variations within resolved H II regions (e.g. James et al. 2016; McLeod et al. 2019). Using narrow-band HST maps, James et al. (2016) show chemical inhomogeneity (∼0.1 dex) on scales smaller than 50 pc within the star forming region Mrk 71. McLeod et al. (2019) also report abundance variations and structures similar to what we observe here in two star-forming regions of the Large Magellanic Cloud, resolved at scales of a few parsec. However, as remarked, for example, by McLeod et al. (2019), abundances derived with a strong line method can have a secondary non-trivial dependence on temperature and density, resulting for instance in falsely low abundances in regions with a high ionisation state (Ercolano et al. 2012; McLeod et al. 2019, 2016). In Fig. 6, we show contours of [S II]/[O III] = 0.5 (dark red), delimiting strongly ionised regions in the nebulae (see Sect. 6). We observe indeed a close correspondence between the contours and areas of low metallicity in all the regions, pointing to the fact that a degeneracy between oxygen abundance and ionisation conditions is likely playing a role closest to the young massive stars. In the future, IFU having a spectral and spatial capability comparable to MUSE, but covering a bluer wavelength range (e.g. BlueMUSE, Richard et al. 2019) will help us to better understand abundances in resolved H II regions by addressing all such degeneracies.

6. Ionisation structure of the H II regions

The ionisation structure of a nebula can be studied through the relative emission of two ions with different ionisation potential. The structure of an ideal ionisation bounded H II region is shown in Fig. 7. The figure illustrates a model generated with CLOUDY for a nebula at the metallicity of our target (Z ∼ 0.008, see Sect. 5), surrounding a 2 Myr YSC with a mass of 10⁴ M_⊙. The simulation is run until the gas interface reaches an optical depth τ = 2. Two zones can be distinguished in the nebula: a central zone with a higher ionisation state, dominated by S++ and O++ and an external zone with a lower ionisation state, where these ions are singly ionised. The ratio of S+/O++ (thick grey) can then be used as a proxy for the optical depth of the nebula, that progresses from optically thin in the inner zone to optically thick in the outer zone, with a sharp cutoff at the Strömgren radius. In reality, the inner zone is often not, or only partially embedded in the optically thick outer envelope, so that a fraction of hydrogen-ionising photons are not absorbed and can escape into the ISM. Studying the spatial extent of the two zones through emission from their most abundant ions can provide us with information about the ionisation structure of an H II region. Pellegrini et al. (2012) exploited this dependency to develop the ionisation parameter mapping (IPM) method and applied it to study the optical depth of H II regions in the Large and Small Magellanic Clouds (LMC and SMC) with narrow-band photometry in [S II], [O III] and Hα. Regions were classified based on the fraction of the central high-ionisation region being surrounded by a low-ionisation envelope in an [S II]/[O III] emission map. Regions completely embedded in the envelope are labelled as optically thick, regions with a low [S II]/[O III] throughout as optically thin and those with an in between morphology as ‘blister nebulae’. Each morphological class was then assigned an escape fraction; we discuss the results of this approach in the light of our work in Sect. 8.

Fig. 7. CLOUDY model illustrating the ionisation structure of an ideal ionisation-bound nebula surrounding a 2 Myr YSC of M = 104 M⊙ at the metallicity of NGC 7793. The radial variation in ionic fraction of H I (dashed black), O III (dashed purple), S II (solid green) and the S II/O III ratio (solid grey) are shown. The nebula consists of two main zones: an inner zone dominated by S III and O III, and an outer zone dominated by S II and O II.

Figure 8 shows a map of the [S II]6716,31/[O III]4959,5007 ratio in our FoV. The map has been reddening corrected assuming an intrinsic Hα/Hβ ratio of 2.863, as described in Paper I. The data have been Voronoi binned to a S/N = 20 in Hβ. The numbered black contours in Fig. 8 indicate the subset of the H II regions – out of the full sample identified in Paper I and shown in Fig. 1 – that we study in detail in Sect. 7. The bar shows increasing ratios of [S II]/[O III], from black to yellow, corresponding to an increasing optical depth of the medium (or, a decrease in its transparency to ionising photons). A ratio of 1 (purple) indicates where the abundance of S+ prevails over O++, as shown in the model in Fig. 7¹².

Fig. 8. [S II]6716,31/[O III]4959,5007 ratio map (reddening corrected), tracing the ionisation structure of the gas. The cyan contours indicate a ratio of 0.5, which we use as fiducial limit for the transition to optically thick gas. Numbered black contours correspond to the H II regions analysed in Sect. 7. The position of WR stars (yellow crosses) and planetary nebulae (white plus symbols) is also indicated.

We observe that the DIG has a lower ionisation state than the H II regions (higher ratios of [S II]/[O III]), as we also observed in Paper I by analysing the ratios of [S II] and [N II] to Hα and in agreement with measurements in nearby galaxies (Haffner et al. 2009). In about half of the regions we observe zones of optically thin gas extending until the regions edge, which could potentially indicate escape channels cleared by stellar feedback. In Sect. 7 we explore the link between the fraction of photons escaping the regions and the ionisation structure traced by [S II]/[O III].

In Table 2 we summarise the visual appearance of the regions in [S II]/[O III], with a similar approach as Pellegrini et al. (2012). Regions are classified into a ‘classical’ ionisation bounded nebula (‘IB’), with an ionisation structure as the model in Fig. 7, or a nebula featuring one (or more) optically thin channels (‘CH’) cutting through the low ionisation zone and reaching the regions edge. We base our classification on the value of [S II]/[O III] along the region contour, as shown in Fig. 9. In the figure, the solid points indicate the median value of the ratio along the contour, and the grey bars span the entire range of observed values. In order to determine which regions provide a potential escape path, we consider [S II]/[O III] = 0.5 as fiducial limit for the transition to optically thick gas: regions having a ratio below this threshold along their contours are classified as CH; all other regions are classified as IB. We test the effect of spatial binning on this classification by using both a finer and a coarser binning pattern, corresponding to S/N = 15 and 25 in Hβ. The classification remains unchanged, suggesting that spatial tessellation does not play a decisive role. We stress that this classification does not take into account the geometry of the nebula and can only be conclusive when having a statistically significant sample of regions that allows to eliminate selection effects.

Fig. 9. Values of [S II]6716,31/[O III]4959,5007 along the contour of the H II regions labelled in Fig. 8. The solid points and error bars show, respectively, the median and the range in values spanned by the ratio. The cyan line indicates a ratio of 0.5, which we use as fiducial limit for the transition to optically thick gas. Regions for which the range in [S II]/[O III] extends below this limit are classified as featuring optically thin channels (CH); all other regions are classified as ionisation bounded (IB).

7. Ionisation budget

We construct a ionisation budget for the entire FoV and inspect the H II regions labelled in Fig. 8 in detail. The latter regions were selected to lie for their most part in the FoV, and not at the edge of the field. We note that we consider here the outermost contours of the H II regions sample, based on a flux brightness cut in Hα = 6.7 × 10⁻¹⁸ erg s⁻¹ cm⁻² spaxel⁻¹ as defined in Paper I. We therefore consider the largest H II region complexes (such as regions 1 and 8 in Fig. 8) in their entirety. This choice is made because dividing the regions within smaller subregions would imply arbitrary choices of their outermost border and division of their fluxes.

For each region, we model the expected Q(H⁰) from the stellar content as described below, and compare it to the observed ionising flux computed from the (reddening corrected) Hα luminosity, following the calibration of Kennicutt (1998):

$$\begin{array}{c}\hfill Q{(H^0)}_{\mathrm{obs}}=7.31\times {10}^{11}L(\mathrm{H}\alpha )[{\mathrm{erg}\mathrm{s}}^{-1}],\end{array}$$

based on the assumption of case B recombination and an electron temperature T_e ∼ 10 000 K.

7.1. O stars

Due to the selection method for the O stars catalogue (see Sect. 2), we lack information such as mass, age or spectral class of the O stars. In order to compute what Q(H⁰) is expected from the O stars population, we therefore have to assume a mass distribution. We consider a Salpeter mass function (Salpeter 1955):

$$\begin{array}{c}\hfill p(m)=A{m}^{-\alpha },\end{array}$$

where p(m) denotes the probability of finding a star of mass m, and α = 2.35. We sample the distribution from m_min = 20 M_⊙, the lowest stellar mass probed by the catalogue, to m_max = 60 M_⊙ with 1000 Montecarlo realisations. This upper limit corresponds to bright HII regions in the Milky Way (Feigelson et al. 2013) and is the highest stellar mass tabulated in the work of Martins et al. (2005) (see below). We also obtain a lower limit for Q(H⁰) of field stars by sampling the Salpeter distribution up to a maximum mass m_max = 30 M_⊙, corresponding to the most massive stars observed in typical galactic star forming regions (Bik et al. 2010; Gennaro et al. 2012; Bik et al. 2012; Ellerbroek et al. 2013). Moreover, stars of mass M > 20 M_⊙ are generally expected to emit enough ionising radiation to be embedded in an H II region, and if observed in the field are most likely runaway objects. For each Montecarlo realisation, we compute the total number of O stars by taking into account the 30% contamination rate and ∼70% completeness of the stellar catalogue (see Sect. 2):

$$\begin{array}{c}\hfill n_{\mathrm{Ostars},\mathrm{tot}}=\mathcal{N}(n_{\mathrm{Ostars}},0.3),\end{array}$$

where 𝒩 denotes a normal distribution. Hereby, we also excluded O stars coinciding with the position of WR stars and YSCs (r < 0.4″). We then draw n_{Ostars, tot} times a mass

$$\begin{array}{cc}& m_i=p(m<m\prime )={\displaystyle {\int }_{m_{\min }}^{m\prime }p(m)\mathrm{d}m=m_{\min }{[C·\mathrm{unif}(0,1)+1]}^{1/(1-\alpha )}}\hfill \\ \hfill & C={\left(\frac{m_{\max }}{m_{\min }}\right)}^{1-\alpha }-1,\hfill \end{array}$$

where unif denotes a uniform distribution. The resulting masses are then converted into Q(H⁰) values by interpolation of Table 1 from Martins et al. (2005), and the total flux of the region Q(H⁰)_tot is stored for each iteration. We finally compute the best Q(H⁰) value and its uncertainty from, respectively, the median m and m ± 1σ of the resulting Q(H⁰)_tot distribution.

7.2. YSC

We determine Q(H⁰) for the YSC using the SLUG code, as described in Sect. 3.2. We obtain a distribution of Q(H⁰)_tot from 1000 Montecarlo realisations: for each realisation, we sample one value from the PDF of each cluster in the region, and sum the resulting quantities. As for the O stars, we then compute the best Q(H⁰) value and its uncertainty from, respectively, the median m and m ± 1σ of the resulting Q(H⁰)_tot distribution. Hereby, we also check that no YSCs are coincident with the position of WR stars, within a radius of 0.4″.

7.3. WR stars

Lastly, we compute the contribution of WR stars. We classify each star as carbon-dominated (WC type) or He II-dominated (WN-type), based on its spectral features as listed in Table 1. We then compute Q(H⁰)_exp from Tables 3 and 4 in Smith et al. (2002), assuming a temperature T = 60 000 K for WN-type stars and T = 120 000 K for WC-type stars. We compute the uncertainty on Q(H⁰) with 1000 Montecarlo realisations of T in the range [40 − 80] × 10³ K for WN-type stars and [100 − 140] × 10³ K for WC-type. Also in this case, we consider the median m and m ± 1σ of the resulting Q(H⁰)_tot distribution.

7.4. Resulting budget

Tables 2 and 3 summarise, respectively, the observed properties and the photoionisation budget for the entire FoV and in detail for the eight regions labelled in Fig. 8. In Table 2 we list the stellar content, the age of the oldest star or cluster in each region, the total GMC mass, median E(B − V) (estimated from the Balmer decrement) and Hα luminosity (reddening corrected), as well as the morphological classification derived in Sect. 6. In Table 3 we indicate the total flux Q(H⁰)_exp modelled from the stellar content, as well as the relative contribution from O stars, YSC, and WR stars, and the observed ionising photon flux Q(H⁰)_obs derived from the reddening corrected Hα luminosity. We also list the ratio of expected to observed flux and the corresponding escape fraction f_esc = 1 − Q(H⁰)_obs/Q(H⁰)_exp. We note that the extremely large uncertainties on regions 5 and 7 are driven by the fact that these regions host (almost) exclusively O stars, for which Q(H⁰) is more loosely constrained. We also note that in region 6 the uncertainty ranges to unphysical values f_esc < 0, indicating that either the models are underestimating the photon flux, or the observed luminosity is being overestimated, possibly due to the reddening correction or to the exact location of the region boundaries.

From Table 3, we see that WR and O stars dominate the contribution to the $Q_{\text{exp},\text{tot}}^0$ value of the regions. Overall, we find an $f_{\mathrm{esc},\mathrm{H}\mathrm{I}\mathrm{I}}=0.{67}_{-0.12}^{+0.08}$ for the entire population of H II regions (black contours in Fig. 1). We also observe that the stellar population in the DIG produces a more than sufficient amount of ionising photons (Q(H⁰)_exp > Q(H⁰)_obs), and that the DIG is therefore consistent with being self-ionised, with $f_{\mathrm{esc},\mathrm{DIG}}=0.{87}_{-0.05}^{+0.03}$. This holds also if considering a maximum mass of 30 M_⊙ for field O stars (last row in Table 3), in which case we find $f_{\mathrm{esc},\mathrm{DIG}}=0.{57}_{-0.17}^{+0.11}$. In our FoV, we observe hence that the sources of ionising photons produce a photon flux that is more than sufficient to explain the emission of the ionised ISM, both within and outside the H II regions.

7.5. Escape fraction from individual H II regions

In this section, we focus in more detail on individual H II regions. In order to better visualise the link between f_esc and the visual appearance of the regions, in Fig. 10 we plot the ratio of Q(H⁰)_exp/Q(H⁰)_obs, and label each region according to the morphological criteria described in Sect. 6 (IB: ionisation bounded, CH: optically thin channels). We do not see any clear trend with visual appearance; we discuss this result in the light of other studies in Sect. 8.

Fig. 10. Ratio of expected to observed ionising luminosity Q(H0) for the regions labelled in Fig. 8 versus their visual appearance: ionisation bounded regions (IB, green stars) vs. regions featuring optically thin channels (CH, blue stars). The dashed horizontal line indicates an fesc = 0.

In Fig. 11, we furthermore investigate the dependence of the Q(H⁰)_exp/Q(H⁰)_obs ratio on the age of the oldest cluster in each region. We consider the latter as a proxy for the age range spanned by the stellar population inside each H II region, assuming that all regions are currently forming stars. We remark that this approximation might not be accurate for regions 6 and 7, in which no GMCs are detected with ALMA and the extinction traced by the ionised gas is very low (see Table 2); however, both regions are still hosting young candidate O stars. When determining the age of the oldest cluster, we use as best estimate the median of each cluster PDF, and we exclude clusters with t > 10 Myr, as such objects are not anymore associated with H II regions, and are most likely line-of-sight objects. For region 5, which is hosting exclusively O stars, we assume an upper limit of 3 Myr for the age. We also colour-code the regions according to their median E(B − V) (left panel) and total GMC mass (right panel), to indicate which regions are more affected by reddening and are still actively forming stars. Dust can indeed absorb part of the LyC photons before they have the possibility to ionise hydrogen, and re-emit them at longer wavelengths, mimicking a larger f_esc. We see that for regions 3 and 8, this could in part explain the f_esc > 0 observed, whereas – given their low dust and GMC content – regions 4 and 5 could be a true leaking region. We comment more on the effect of photons removal by dust in Sect. 8. In our small sample of regions and due to the large uncertainties, we do not find any trend between f_esc and the age of the stellar population in the region. This is discussed in more detail in Sect. 8.

Fig. 11. Ratio of expected to observed ionising photon flux Q(H0) for the regions labelled in Fig. 8, versus the age of the oldest cluster in each region (best-value estimate obtained from the median of the age PDF; masking cluster with ages > 10 Myr). For region 5, which is hosting exclusively O stars, we assume an upper limit of 3 Myr for the age. The dashed horizontal line indicates an fesc = 0. Left panel: points are colour-coded according to the median E(B − V) in each region (Balmer decrement estimate). Right panel: regions are colour-coded according to their total GMC mass. Empty black circles indicate regions hosting no GMCs.

8. Discussion

In Paper I, we found a DIG fraction f_DIG, obs = 0.15 from the Hα luminosity in our FoV¹³. This is in good agreement with the fraction of diffuse gas expected by modelling the stellar and cluster population:

$$\begin{array}{c}\hfill f_{\mathrm{DIG},\exp }=\frac{Q_{\exp ,\mathrm{DIG}}^0}{Q_{\exp ,\mathrm{DIG}}^0+Q_{\exp ,\mathrm{HII}}^0}\gtrsim 0.17.\end{array}$$

Moreover, we assessed from our BPT diagram analysis in Paper I that the bulk of the DIG emission in our FoV was consistent with being photoionised. Here we confirm that the DIG is consistent with being self-ionised by field stars and clusters, with an overabundance of ionising photons $f_{\mathrm{esc}}\sim 0.{87}_{-0.05}^{+0.03}$; we confirm this trend also when considering our lower-limit estimate. This result is perhaps surprising considering the fact that we observe in Fig. 8 that the DIG is not optically thin to Lyman continuum emission (high [S II]/[O III] ratios). However, if the DIG has a clumpy structure, the emission line ratio could be tracing the denser, optically thick clumps; the apparent lack of Hα emission can originate from photons escaping through the dilute, optically thin interclump medium, which through its low density will not contribute significantly to the observed emission line intensities. We also observe that overall the H II regions are leaking ionising photons at a rate $f_{\mathrm{esc}}\sim 0.{67}_{-0.12}^{+0.08}$, and that seven out of eight have an f_esc ≳ 0.3.

We note that the uncertainties on f_esc are rather large, as we try to account for a variety of effects, including: the contamination rate ∼30% in the O stars catalogue, the large range in Q(H⁰) obtained by Montecarlo sampling O stars of different mass, the range in values spanned by the PDFs of the clusters and the assumption on the temperature of the WR stars. The large uncertainties recovered are consistent with the work of Niederhofer et al. (2016): by varying parameters of synthetic YSCs, the authors conclude that escape fractions derived from broadband photometry data are typically dominated by uncertainties in the spectral types of the stars. However, even when considering the uncertainties, we find f_esc > 0 for the DIG overall and for six out of the eight H II regions inspected (see Table 3 and Fig. 10).

An additional factor that is not taken into account in the computation of f_esc is the effect of absorption of LyC photons by dust within the H II regions before these photons have the chance of ionising hydrogen atoms. In Sect. 4, we have compared the extinction map derived from the Hα/Hβ ratio with 24 μm emission from Spitzer/MIPS (Fig. 5, right panel), and found a good correspondence, indicative that dust and gas are well mixed (see e.g. Choi et al. 2020). The resolution of mid-IR data is currently insufficient to confirm potential dust absorption in proximity of the sources of Lyman continuum photons at the distance of our target, and will only become possible with new generation IR telescopes such as JWST. However, in local H II regions, it is typically observed that a large fraction of LyC photons contributes to hydrogen ionisation. This fraction anti-correlates with metallicity (e.g. Inoue et al. 2001), and at solar metallicity, an increasing fraction of Lyman continuum photons are absorbed by dust, instead. In the metallicity range spanned by our target (12 + log(O/H) ∼ 8.25 − 8.5, see Sect. 5), Inoue et al. (2001) estimate that ≳80% Lyman continuum photons are absorbed by neutral hydrogen. Therefore, we expect that the effect of dust absorption in the H II regions in our FoV can reduce Q(H⁰)_obs of up to ∼20%, which lies well within the given uncertainties.

Other studies targeting H II regions in nearby galaxies at high resolution found similar escape rates. Doran et al. (2013) have estimated an f_esc > 0.5 in the 30 Dor star forming region in the LMC, based on a complete spectroscopic census of hot, luminous stars. In a sample of H II regions in the LMC and SMC studied with narrowband photometry, Pellegrini et al. (2012) have found f_esc > 0.4 for all the regions. McLeod et al. (2019) have studied in more detail two of the largest H II region complexes in the LMC with MUSE, finding f_esc > 0.2 for the complexes and the respective sub-regions. Finally, in two regions recently observed with MUSE and HST in the nearby dwarf galaxy NGC 300, McLeod et al. (2020) have estimated an f_esc ≳ 0.3.

Following the non-negligible escape of Lyman continuum photons from H II regions, several studies of local and nearby galaxies have proven the DIG consistent with being photoionised. For example, Pellegrini et al. (2012) have estimated that both in the LMC and SMC the total galactic escape fraction would be sufficient to account for the observed f_DIG, and that, when considering the additional contribution of field stars, the galaxies could leak a substantial amount of ionising radiation into the circumgalactic medium. Similarly, in a recent MUSE+HST study of the Antennae merger system, Weilbacher et al. (2018) have estimated a sufficient fraction of Lyman-continuum leakage in order to explain the amount of DIG observed.

We do not find any correlation between ionising photon leakage in the HII regions and their visual appearance (‘classic’ ionisation bound regions vs. regions featuring optically thin channels). This is somewhat in contrast with the results of Pellegrini et al. (2012) that, by comparing the fiducial f_esc assigned to H II regions in the LMC and SMC based on their appearance in [S II]/[O III] (see Sect. 6) with previous f_esc estimated based on the stellar population, found a good agreement overall. On the other hand, a systematic trend is not apparent from the MUSE results of McLeod et al. (2019) on the LMC, that used the [O II]7320,7330/[O III]4959,5007 ratio as optical depth tracer. We interpret this apparent tension as likely being a simple effect of 3D geometry in number limited samples.

In our small sample of regions, we do not find evidence for a correlation between age range of the stellar population and f_esc. Increasing the sample size will be essential to investigate any such correlation, that could be expected in light of the results of galactic and cosmological simulations. Kim et al. (2013) simulated ionising radiation and supernova feedback in a low-redshift galactic disk, at a high spatial resolution ∼4 pc. The study indicated that photons emitted at later ages are more likely to escape the star forming region, as the gas surrounding the region is dispersed by ionising radiation and supernova feedback. A similar result was obtained in the recent high-resolution cosmological simulations of Ma et al. (2020), that studied ∼30 zoom-in simulations of galaxies at z ≳ 5. The simulations suggest that ionising photons are preferably leaking from star forming regions containing a kpc-size superbubble, likely created by clustered SNe set off by stars of age > 3 Myr. As the bubble expands, new stars are formed at its edge, inside a dense shell of compressed gas. The shell keeps expanding while forming stars and, as a consequence, the young stars formed in it end up inside the superbubble and are able to fully ionise channels of low-column density pre-cleared by feedback from the previous population of stars. Therefore, regions hosting stars spanning a large range in age seem to be advantaged in leaking photons. Also observations have already pointed to the importance of the local star formation history in shaping the ionisation conditions in H II regions. By analysing a sample of ∼5000 H II regions from CALIFA (Calar Alto Large Integral Field Area survey, Sánchez et al. 2012), Sánchez et al. (2015) found for example a correlation between the position occupied by the regions in BPT diagrams and the age and metallicity of their stellar population.

Finally, we comment on the impact of using different stellar population model in our cluster_slug analysis on the cluster physical properties. The use of non rotating and single star models produces a lower limit to the total estimated Q(H⁰) and therefore to f_esc, in particular at ages larger than 3−4 Myr. Binaries and rotating stars have the same effect, that is they increase the production of ionising photons at ages older than 3 Myr with respect to predictions from non-rotating single star models (see Leitherer et al. 2014; Götberg et al. 2019). Any existing trend between f_esc and age range of the stellar population would therefore be reinforced by different model assumptions, as also reported in the simulations of Ma et al. (2020).

9. Conclusions

We have studied the ionised gas in the nearby galaxy NGC 7793 with MUSE, and complemented our observations with HST and ALMA data tracing the stellar content and the molecular gas.

We have constructed a census of YSCs, O stars, and WR stars in the MUSE FoV, and studied the properties of the stellar population. We modelled the age, mass, and ionising flux of YSCs with the stochastic stellar population synthesis code SLUG, finding that clusters located in H II regions have a higher probability to be younger, less massive, and to emit a higher number of ionising photons than clusters in the field.

We have investigated the link between the stellar population and the dense and ionised gas. We have contrasted the reddening map derived from individual stellar extinctions with the one constructed from the Hα/Hβ ratio and found that they compare well, but that the latter correlates better with the position of GMCs traced by ALMA. We have estimated an oxygen abundance for the H II regions from the MUSE data, using the S-strong line method from Pilyugin & Grebel (2016). We found a median abundance of 12 + log(O/H) ∼ 8.37 with a scatter of 0.25 dex, in agreement with the previous estimate by Pilyugin et al. (2014). The abundance map appears to be rich in substructures, especially surrounding YSCs and WR stars. We caution however against possible degeneracies with for example the temperature and density, and we do indeed observe a correlation with the ionisation state of the regions.

We have studied the ionisation structure of the H II regions using the [S II]6716,31/[O III]4959,5007 ratio as a proxy for the optical depth. We have focused on a subset of regions and classified them based on the value of [S II]/[O III] along their border into ionisation bounded or as featuring channels of optically thin gas. Finally, we have compiled a photoionisation budget for the entire FoV and for the subset of H II regions.

Overall, we find an escape fraction $f_{\mathrm{esc}}=0.{67}_{-0.12}^{+0.08}$ for the population of H II regions, and that the DIG in our FoV is more than consistent with being self-ionised, with an $f_{\mathrm{esc}}=0.{87}_{-0.05}^{+0.03}$. This holds even when considering a lower-limit estimate for the DIG flux, derived by assuming a maximum mass of 30 M_⊙ for the field O stars. We furthermore find that the f_DIG, exp ≳ 0.17 obtained by modelling the DIG stellar population is in good agreement with the DIG fraction derived from the observed Hα luminosity in Paper I, f_DIG, obs = 0.15.

We observe an f_esc ≳ 0.3 in seven out of the eight studied regions, and investigate how the f_esc is linked to the regions properties. Among others, we do not find any trend between f_esc and the visual appearance of the regions; we read this as the effect of 3D geometry in our number limited sample.

Acknowledgments

We would like to thank the anonymous referee for the constructive feedback on the second draft of the manuscript. This work is based on observations collected at the European Southern Observatory under ESO programme 60.A-9188(A). A.A. acknowledges the support of the Swedish Research Council, Vetenskapsrådet, and the Swedish National Space Agency (SNSA). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 757535). This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

All Tables

All Figures

	Fig. 1. Stellar census in the MUSE FoV. The location of candidate O stars (white circles), YSCs (green triangles) and WR stars (blue crosses) is indicated, overlaid on a map of the Hα emission. Black contours indicate the H II regions identified in Paper I, light blue filled contours denote the position of GMCs.
In the text

	Fig. 2. Spatial distribution of YSC ages (left panel) and masses (right panel) in the MUSE FoV, overlaid on the Hα map shown in Fig. 1. Circles and diamonds indicate, respectively, clusters younger and older than 10 Myr. Ages and masses are best estimates computed from the median of the corresponding PDF (see Sect. 3.2). Black contours indicate the H II regions identified in Paper I.
In the text

	Fig. 3. Posterior probability distributions of YSC age (left) mass (centre) and ionising flux (right) for clusters located in H II regions (in blue) and for field clusters (in orange). Thick lines indicate the combined PDFs, thin lines correspond to the PDFs of single clusters. The grey dashed lines indicate the location of the main modes of the combined distributions.
In the text

	Fig. 4. Distribution of YSC mass (clusters with t < 10 Myr, in blue), GMC mass (in green) and H II regions luminosity (in black, reddening corrected) in our FoV. YSC masses are best estimates computed from the median of the corresponding PDF (see Sect. 3.2). The histograms are normalised to a peak value of 1.
In the text

Fig. 5. Extinction map derived from the ionised gas observed with MUSE assuming an intrinsic Hα/Hβ ratio (from Paper I) compared to: Left panel: contours of the extinction map derived by Kahre et al. (2018) from the individual stellar extinctions obtained from HST data. The contours correspond to values of E(B − V) = 0.2 − 0.5, in steps of 0.033. Right panel: 24 μm emission from Spitzer/MIPS. The contours correspond to a flux of [1.5, 2, 3, 4, 6, 9, 11] MJy sr−1. The thick orange contours indicate the position of GMCs.

In the text

Fig. 6. Oxygen abundance for the subset of H II regions labelled in Fig. 8, determined with the S-strong line calibration of Pilyugin & Grebel (2016). The location of YSCs is indicated with filled circles, colour-coded based on their age (best estimates derived from the median of each PDF). Red crosses and purple stars correspond to the position of WR stars and PNe, respectively. Dark red contours indicate a ratio of [S II]/[O III] = 0.5 (see Fig. 8).

In the text

Fig. 7. CLOUDY model illustrating the ionisation structure of an ideal ionisation-bound nebula surrounding a 2 Myr YSC of M = 104 M⊙ at the metallicity of NGC 7793. The radial variation in ionic fraction of H I (dashed black), O III (dashed purple), S II (solid green) and the S II/O III ratio (solid grey) are shown. The nebula consists of two main zones: an inner zone dominated by S III and O III, and an outer zone dominated by S II and O II.

In the text

Fig. 8. [S II]6716,31/[O III]4959,5007 ratio map (reddening corrected), tracing the ionisation structure of the gas. The cyan contours indicate a ratio of 0.5, which we use as fiducial limit for the transition to optically thick gas. Numbered black contours correspond to the H II regions analysed in Sect. 7. The position of WR stars (yellow crosses) and planetary nebulae (white plus symbols) is also indicated.

In the text

Fig. 9. Values of [S II]6716,31/[O III]4959,5007 along the contour of the H II regions labelled in Fig. 8. The solid points and error bars show, respectively, the median and the range in values spanned by the ratio. The cyan line indicates a ratio of 0.5, which we use as fiducial limit for the transition to optically thick gas. Regions for which the range in [S II]/[O III] extends below this limit are classified as featuring optically thin channels (CH); all other regions are classified as ionisation bounded (IB).

In the text

	Fig. 10. Ratio of expected to observed ionising luminosity Q(H0) for the regions labelled in Fig. 8 versus their visual appearance: ionisation bounded regions (IB, green stars) vs. regions featuring optically thin channels (CH, blue stars). The dashed horizontal line indicates an fesc = 0.
In the text

Fig. 11. Ratio of expected to observed ionising photon flux Q(H0) for the regions labelled in Fig. 8, versus the age of the oldest cluster in each region (best-value estimate obtained from the median of the age PDF; masking cluster with ages > 10 Myr). For region 5, which is hosting exclusively O stars, we assume an upper limit of 3 Myr for the age. The dashed horizontal line indicates an fesc = 0. Left panel: points are colour-coded according to the median E(B − V) in each region (Balmer decrement estimate). Right panel: regions are colour-coded according to their total GMC mass. Empty black circles indicate regions hosting no GMCs.

In the text