© ESO 2022

1. Introduction

H II regions are key tracers of the rate of star formation and the degree of chemical enrichment of the interstellar medium (ISM) in galaxies across the entire range of cosmological distances accessible to current instrumentation (Maiolino & Mannucci 2019; Kewley et al. 2019). Detailed studies of H II regions in nearby galaxies have quantified the effect of feedback from massive stars on their environment (e.g. Lopez et al. 2011, 2014; Pellegrini et al. 2012; McLeod et al. 2019, 2020; Barnes et al. 2020; Olivier et al. 2021; Della Bruna et al. 2021) and allowed the development of empirical relations that describe the conversion of cold gas into stars (e.g. Kennicutt 1989; Schinnerer et al. 2019; Chevance et al. 2020). Feedback from massive stars, in the form of radiation, winds, and supernovae, is a crucial ingredient in determining the state of the ISM in which new stars form (Hopkins et al. 2014; Klessen & Glover 2016; Gatto et al. 2017).

A large fraction of the ionising radiation from massive stars may escape from H II regions in their immediate vicinity and contribute to powering the diffuse ionised gas (DIG): a tenuous plasma with low electron densities (n_e ∼ 10⁻¹ − 10⁻² cm⁻³) and temperatures slightly higher than those of typical H II regions (∼10⁴ K). Understanding what fraction of the DIG is powered by leaking radiation from H II regions, and therefore the ionising photon escape fraction from regions of star formation, is of fundamental importance for drawing a full budget of the feedback mechanisms acting to disrupt giant molecular clouds, and therefore driving the timescales of the cycle of star formation in galaxies (Chevance et al. 2022). The distribution of ionised hydrogen in the DIG is traced by the Hα recombination line. In general, the DIG is found to contain roughly half of the total Hα emission in galaxies (Oey et al. 2007; Kreckel et al. 2016; Chevance et al. 2020). This emission is often observed to be spatially associated with nearby H II regions (Ferguson et al. 1996; Zurita et al. 2000, 2002; Madsen et al. 2006; Seon 2009; Howard et al. 2018; Pellegrini et al. 2020a), supporting the idea that the DIG may be ionised by leaking radiation. Detailed assessments of the energetics involved support this argument for the DIG in the Milky Way (Reynolds 1990; McKee & Williams 1997; Haffner et al. 2009).

In external galaxies, collisionally excited optical line ratios constitute our main window into the physical conditions of the DIG and can shed light onto its ionisation source. In nearby edge-on galaxies (Otte et al. 2001; Jones et al. 2017; Levy et al. 2019), spectroscopic observations show an increase in emission of low-ionisation species in the DIG with respect to H II regions. In particular, the line ratios [S II]λλ6717,31/Hα, [N II]λ6584/Hα, and [O I]λ6300/Hα are found to increase with scale height from the midplane, or as a function of decreasing Hα surface brightness (Σ_Hα; Hill et al. 2014). Models for gas photoionised by leaking radiation from H II regions in the low-density regime are able to explain the increase in low-ionisation line ratios such as [S II]/Hα, [N II]/Hα, and [O I]/Hα. Furthermore, intervening absorption from neutral gas preferentially absorbs photons near the hydrogen ionisation potential, therefore hardening the ionising spectrum, which can also lead to an increase in low-ionisation line ratios (Wood & Mathis 2004).

However, leaking H II region models are generally unable to reproduce the relative strengths of high-ionisation lines (e.g. [O III]λ5007/Hβ) in the extraplanar DIG without additional sources of heating (Rand 1998; Reynolds et al. 1999). Several physical processes have been suggested as sources for an additional heating term (e.g. shocks: Collins & Rand 2001, turbulence: Binette et al. 2009, or magnetic reconnection: Lazarian et al. 2020). An alternative solution, provided by Flores-Fajardo et al. (2011) in a study of the prototypical edge-on galaxy NGC 891, is that the DIG emission can be modelled as a combination of photoionisation from leaky H II regions and hot low-mass evolved stars (HOLMES). HOLMES are low- to intermediate-mass stars (0.8 − 8.0 M_⊙) in all the stages of stellar evolution subsequent to the asymptotic giant branch (post-AGB), including the hydrogen-burning, white-dwarf-cooling, and intermediate phases (Stanghellini & Renzini 2000; Stasińska et al. 2008). The term ‘post-AGB’ has been used in other studies (Binette et al. 1994; Belfiore et al. 2016; Byler et al. 2019) to refer to the same population of sources. HOLMES have a harder ionising spectrum than H II regions and are therefore capable of producing higher [O III]/Hβ ratios, even with relatively weak radiation fields. HOLMES are a prime candidate to explain the [O III]/Hβ line ratios in the DIG of star-forming galaxies because, unlike the other potential candidates, their presence in galaxies is undisputed.

Data from kiloparsec-resolution integral field spectroscopy (IFS) surveys of nearby galaxies, such as CALIFA (Calar Alto Legacy Integral Field spectroscopy Area survey, Sánchez et al. 2012), SAMI (Sydney-Australian-Astronomical-Observatory Multi-object Integral-Field Spectrograph, Croom et al. 2012), and MaNGA (Mapping Nearby Galaxies at APO, Bundy et al. 2015), has demonstrated the importance of photoionisation by HOLMES in early-type galaxies (Sarzi et al. 2010; Singh et al. 2013) and central regions of massive late-type galaxies (Belfiore et al. 2016), which are devoid of recent star formation. These systems show line ratios that fall in the low-ionisation (nuclear) emission-line region (LI[N]ER) corner of classical Baldwin–Phillips–Terlevich (BPT; Baldwin et al. 1981; Kewley et al. 2001, 2006) diagrams, as expected from ionisation sources other than young stars.

There is also increasing evidence with regard to the importance of HOLMES in explaining the observed line ratios in the DIG of star-forming galaxies. Zhang et al. (2017) invoke a contribution from HOLMES to explain the high [O III]/Hβ line ratios observed at low Σ_Hα in star-forming galaxies from the MaNGA survey. Other authors (Lacerda et al. 2018; Vale Asari et al. 2019; Espinosa-Ponce et al. 2020) use data at kiloparsec resolution to effectively redefine the DIG as regions with a low equivalent width of Hα emission (e.g. EW(Hα) < 3 Å), which models suggest may be consistent with ionisation from HOLMES. However, this redefinition is problematic for studies interested in determining the ionisation sources of the DIG since it defines the DIG a priori as the ionised gas component consistent with ionisation from HOLMES. Unfortunately, at the typical resolution provided by CALIFA, SAMI, or MaNGA (∼0.8 for CALIFA to 2 kpc for MaNGA), regions classified as star-forming consist of an agglomeration of several bright H II regions (Mast et al. 2014; Poetrodjojo et al. 2019), and it is not possible to spatially resolve these regions from any nearby, potentially associated, DIG.

This uncertainty in identifying the DIG leads to disagreement over the importance of the DIG in determining chemical abundances. Abundance measurements rely on photoionisation models, which do not generally include DIG contamination, or on empirical relations calibrated on bright H II regions (themselves generally expected to suffer only negligible DIG contamination). The [N II]/Hα line ratio, for example, is sometimes used as a metallicity tracer, especially at high redshift (Pettini & Pagel 2004). Considering that roughly half of the Hα emission may originate outside H II regions, and that the [N II]/Hα ratio is enhanced in the DIG, the abundance derived without taking the DIG into account inevitably results in a biased measurement (e.g. Zhang et al. 2017; Poetrodjojo et al. 2019). A reassessment of the relative importance of leaking radiation and HOLMES in powering the DIG in star-forming disc galaxies is therefore urgently needed to make progress in quantifying the H II regions’ ionising photon escape fraction and to measure chemical abundances accurately.

Making such advances would require sensitive, high-spatial-resolution spectral mapping of nearby star-forming galaxies. The ideal dataset would have sufficient spatial resolution to resolve the average separation between H II regions (∼100 − 200 pc; Chevance et al. 2020) in order to allow H II region emission to be spatially separated from the DIG.

A dataset suitable for a high-spatial-resolution study of the DIG is now available thanks to the Physics at High Angular resolution in Nearby GalaxieS (PHANGS)¹ project, which includes mapping of 19 nearby galaxies at < 100 pc resolution obtained via a Large Programme (PI: E. Schinnerer) with the MUSE (Multi Unit Spectroscopic Explorer) integral field spectrograph on the ESO/Very Large Telescope (VLT; Emsellem et al. 2022). In this paper we make use of the PHANGS–MUSE dataset to study the spatial structure of the DIG and its line ratios. In Sect. 2 we describe the data and our analysis tools and present our methodology for separating H II regions from DIG based on Hα morphology. We present the general characteristics of the DIG in our galaxy sample in Sect. 3. We then test the canonical assumption that the DIG is powered by leaked radiation from H II regions by comparing our observations with simple predictions of the spatial distribution of the DIG under this scenario (Sect. 4). We describe the change in DIG line ratios as a function of proxies for metallicity and the ionisation parameter in Sect. 5 and explain the trends via photoionisation models. In Sect. 6 we present our main take-away points with regards to the powering source of the DIG and what role HOLMES play in its ionisation state.

For the rest of the paper we make use of the following shorthand notation: [N II] ≡ [N II]λ6584; [S II] ≡ [S II]λ6717 + [S II]λ6731; [O III] ≡ [O III]λ5007; and [O I] ≡ [O I]λ6300.

2. Sample and data analysis

2.1. Observations and data reduction

PHANGS–MUSE consists of 19 galaxies, selected to be nearby (D < 20 Mpc), close to the star-formation main sequence, and moderately inclined. The main properties of the sample, which spans the stellar mass range log(M_⋆/M_⊙) = 9.4 − 11.0, are summarised in Table 1. A description of the PHANGS–MUSE sample selection, observations, and associated data is presented in Emsellem et al. (2022). In the rest of this section we offer a brief summary.

The MUSE IFS instrument (Bacon et al. 2010) at the ESO/VLT was used to cover the star-forming part of the disc. For each galaxy we obtain a contiguous mosaic of several 1′×1′ MUSE pointings (between 3 and 15 pointings per galaxy). The majority of the data was obtained via a dedicated ESO Large Programme (PI: E. Schinnerer), running from ESO periods 100–106 (October 2017–March 2021). Three extra MUSE pointings (one each in NGC 1566, NGC 1385 and NGC 4321) are not included here because their reduction was pending during the development of this work, but they are included in the PHANGS public data release. Eight galaxies in the sample were observed by making use of the MUSE wide-field ground layer adaptive optics mode, while the remaining ones were observed in natural seeing. All the data were reduced via the PYMUSEPIPE python package², which wraps the ESOREX MUSE reduction recipes (Weilbacher et al. 2020) and additionally provides routines for astrometric alignment, mosaicking and point spread function (PSF) homogenisation.

Astrometry and flux calibration were performed by comparing the MUSE data with ESO/MPG 2.2 m Wide Field Imager and LCO/DuPont Direct CCD (for NGC 7496 only) broadband images (Razza et al., in prep.), whose astrometry was matched to that of Gaia data release 2 (Gaia Collaboration 2018). For each galaxy we produced a sky-subtracted, flux-calibrated mosaicked datacube, with a spaxel size of 0.2″. Masks for foreground stars were generated based on a combination of the Gaia point-source positions and identification of the rest-frame Ca II triplet absorption features at 8498, 8542, and 8662 Å. The Ca II triplet absorption helped to avoid masking compact H II regions and galactic nuclei, which may appear as point sources in Gaia.

The effective R-band PSF of the MUSE data, as measured in our data as detailed in Emsellem et al. (2022), varies from 0.4″ (only achieved in wide-field adaptive optics mode) to 1.2″. The median PSF for each mosaic is provided in Table 1, together with the offsets to the largest (max) and smallest (min) PSF measured among the MUSE pointings that constitute the final mosaic. Our median PSF corresponds to an average resolution across the sample of ∼50 pc.

Throughout the paper, radial distances are deprojected assuming the inclinations and positions angles derived by Lang et al. (2020) from kinematic analysis of molecular gas data, or analysis of near-IR imaging when not available or when the fit to the velocity field was deemed unreliable³ (see Table 1).

2.2. Spectral fitting

The mosaicked cubes are processed by a custom-built data analysis pipeline (Emsellem et al. 2022), based on the public GIST (Galaxy IFU Spectroscopy Tool; Bittner et al. 2019) software package, and modified by F. Belfiore and I. Pessa for obtaining maps of stellar kinematics, stellar population properties (including stellar mass surface density), and emission-line fluxes and kinematics. The pipeline makes use of the penalised pixel fitting python package (PPXF; Cappellari & Emsellem 2004; Cappellari 2017), in three subsequent steps: (i) stellar kinematics extraction, (ii) stellar population analysis, and (iii) determination of emission-line parameters.

The stellar kinematics and stellar population steps were performed on binned spectra to increase the continuum signal-to-noise ratio (S/N) over that of single spaxels. In particular, we made use of Voronoi binning, as implemented in the python package VORBIN (Cappellari & Copin 2003), to reach a continuum S/N of 35 in the wavelength range 5300−5500 Å. To fit the stellar continuum we used E-MILES simple stellar population (SSP) models (Vazdekis et al. 2016) generated with a Chabrier (2003) initial mass function, BaSTI isochrones (Pietrinferni et al. 2004), eight ages (0.15 − 14 Gyr), and four metallicities ([Z/H] = [ − 1.5, −0.35, 0.06, 0.4]) to determine the stellar kinematics. A slightly larger grid of 78 templates was used for the stellar population analysis. We fitted the wavelength range 4850 − 7000 Å in order to avoid strong sky residuals in the redder part of the MUSE wavelength range. The templates were convolved to the MUSE spectral resolution, as derived by Bacon et al. (2017).

For the emission-line analysis, we performed an additional, simultaneous fit of stellar continuum and line emission. We also performed this fitting stage with PPXF, exploiting the extensive support for emission-line fitting available within that package, by treating emission lines as additional Gaussian templates. The emission-line analysis was performed on individual spaxels to maximise the resolution of the emission-line maps. During this step, the kinematics of the stellar continuum was fixed to that determined from the binned spectra analysis. All fluxes and maps presented in this work were corrected for foreground Galactic extinction using the Cardelli et al. (1989) extinction law and the E(B − V) values from Schlafly & Finkbeiner (2011). We reach a typical 3σ Hα surface brightness (Σ_Hα) sensitivity of 1.5 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² in 0.2″ spaxels. More details on the data analysis and quality assessment can be found in Emsellem et al. (2022).

2.3. Morphological definitions of DIG and H II regions

H II regions present a complex challenge for photometry. Nebulae show irregular morphologies, including shells and arcs, and do not have well-defined boundaries. Sophisticated algorithms have been developed for the task of isolating bright regions from the diffuse background and defining their boundaries. The spatial resolution of the PHANGS–MUSE data (median 50 pc) is in general not sufficient to resolve the internal structure and size of typical H II regions (∼1 − 100 pc), except for the most extended, most luminous giant nebulae with size > 100 pc. For comparison, 30 Doradus in the LMC, the most luminous giant H II in the Local Group, is ∼400 pc in size (Kennicutt 1984).

In this work we employed a pragmatic definition of H II regions as bright regions with relatively sharp boundaries at the ∼50 pc resolution of our dataset. We defined H II region masks by running the IDL tool HIIPHOT (Thilker et al. 2000) in conjunction with a point source finder on PSF-homogenised Hα line maps. HIIPHOT is optimised for photometry of extended regions, and separates nebulae from the surrounding DIG by generating seed regions around bright peaks, and then growing them until a termination criterion determined by the spatial gradient of the Σ_Hα profile is reached. By visually inspecting the resulting H II region masks we determined that HIIPHOT was systematically missing a large population of compact point sources, corresponding to unresolved H II regions. In order to obtain a more complete catalogue, we generated masks for these regions using the point-source finder DAOSTARFINDER, provided by the python PHOTUTILS module. Our final segmentation maps consists of the union of the HIIPHOT and DAOSTARFINDER nebulae masks, and contains ∼1200 nebulae per galaxy on average. The segmentation maps are shown for all the galaxies in the survey in Figs. A.1 and A.2.

The DIG is defined as line emission originating from outside the area covered by the masks used to identify ionised nebulae. Vice versa, the H II region fluxes are taken as simply the flux within the H II region masks, without any correction for the local DIG background. The final catalogue of nebulae is dominated by number by bona fide H II regions, but includes contamination from supernovae remnants, and planetary nebulae (PNe). A detailed description of the catalogue, and further discussion of algorithmic choices made to define H II regions will be presented in a future paper. In the rest of the section we address some of the main limitations, related to blending of nebulae near large complexes, faint nebulae merging with the background, and the difference in kinematics between DIG and H II regions.

Blending of H II regions in crowded regions may affect our ability to determine number counts of H II regions in luminosity bins (i.e. determine the H II region luminosity function, Santoro et al. 2022), but should have a small impact on separating H II regions and DIG, since in such regions we expect most of the area to be dominated by H II regions. It should also be noted that defining borders of H II regions in the case of large complexes is an intrinsically ill-defined problem, as the gas may not be easily associated with a unique source of ionising photons. In our analysis, the size of H II regions is set by the termination gradient value we employed in HIIPHOT, which is driven by the noise level of the data and visual expert assessment of the performance of the algorithm, rather than by physical considerations. Choosing a different termination gradient would extend the H II region masks into the surrounding area, which we have classified as DIG. We experimented with different values of the termination gradient, and find that the line ratios in H II regions are unaffected, while some of the brighter DIG regions get included into the H II region masks.

A more crucial limitation is the blending of faint H II regions into the diffuse background. The completeness limits of our H II region catalogue are estimated by Santoro et al. (2022) to be in the range log(L_Hα/erg s⁻¹) ∼ 36 − 37, depending on the distance of the galaxy. H II regions fainter than log(L_Hα/erg s⁻¹) ∼ 36 have not been studied extensively in external galaxies. In Local Group galaxies, such faint H II regions are found to only contribute marginally to the total Hα flux (< 10%, Kennicutt et al. 1989 for M 33 and the Magellanic Clouds, although M 31 may be an exception; see Azimlu et al. 2011). Higher-resolution data (provided e.g. by Hubble Space Telescope narrow-band imaging) would be necessary to estimate the contribution from faint, undetected H II regions to the DIG.

Planetary nebulae are bright in [O III] but comparatively faint in Hα, and are thus not efficiently detected by our algorithm. In Scheuermann et al. (2022) the authors perform a detailed study of the planetary nebula luminosity function (PNLF) in the PHANGS-MUSE galaxy sample, applying a point source detection algorithms to an [O III] flux map generated from our MUSE data to detect PNe. Comparing with their catalog of PNe, we find that only 22% are detected in Hα and are therefore identified by our algorithm as nebulae. The remaining PNe identified in [O III] are merged into to our DIG mask. We can compute the fractional contribution of PNe to the total DIG flux by integrating the [O III] PNLF within suitable limits and normalising it according to the total number of detected PNe in our galaxies, taken from Scheuermann et al. (2022). If we adopt the PNLF parametrisation of Ciardullo et al. (1989) and integrate it from the bright magnitude cutoff down to sources 8 mag fainter, we find that the total [O III] flux contributed by PNe is on average 1% of the [O III] flux observed in our DIG masks. The contribution of PNe to the flux in other emission lines considered in this work, including Hα, will be even smaller. We therefore conclude that, even though our DIG fluxes certainly include a contribution from PNe, this contribution is negligible. Our estimate is in agreement with Yan & Blanton (2012), who perform a similar computation, and find that PNe fall short by about 1.5 dex in explaining the observed [O III] emission in early-type galaxies.

Detailed decomposition of line profiles may offer a promising tool to distinguish DIG from H II regions, as the DIG is generally found to lie in a thicker, and more slowly rotating structure than the star-forming disc (Fraternali et al. 2004; Heald et al. 2006; Boettcher et al. 2017; Levy et al. 2018; den Brok et al. 2020). However, the rotational velocity lag between DIG and H II regions is of the order of a few tens of km s⁻¹, which is not large enough to be easily detectable at the average spectral resolution of the MUSE data (spectral resolution of ∼115 km s⁻¹ at the wavelength of the Hα line). In this paper, therefore, we do not use kinematic information to define the DIG.

2.4. Spatial binning optimised for low-surface-brightness gas

At the depth of our data, Hα is detected throughout a large fraction of the surveyed area (Hα is detected at the 3σ level in 95% of the 0.2″ spaxels in the MUSE observations of our galaxies, Emsellem et al. 2022), but even the strongest metal lines, such as [S II] and [O III], fall below the detection limit in low-surface-brightness regions. This problem can be better appreciated by inspecting Fig. 1 (panels a–c), where we show maps of Hα flux, log([S II]/Hα) and log([O III]/Hβ) for the galaxy NGC 4535. The black contour corresponds to Σ_Hα = 8 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² or about a 16σ detection for Hα in 0.2″ spaxels. This surface-brightness level was chosen to correspond roughly to the lower 10th percentile of Σ_Hα for H II regions in this galaxy, and therefore constitutes a useful visual delimiter between bright areas associated with catalogued nebulae and the DIG. In Fig. 1, areas with S/N < 3 in the relevant emission lines are masked out. The figure readily demonstrates the impossibility of deriving representative line ratios in the DIG using individual spaxels.

Fig. 1. Maps of the Hα flux, binning scheme, and key line ratios for individual spaxels and after Voronoi binning in NGC 4535. a) Map of log(Hα/10−17 erg s−1 cm−2 arcsec−2), at the native 0.2″ single-spaxel resolution. The circular areas that appear white in this and other panels are locations of masked foreground stars. b) Map of log([S II]/Hα) at single-spaxel resolution. White areas have S/N < 3 for the relevant emission lines. The black contour corresponds to ΣHα = 8 × 10−17 erg s−1 cm−2 arcsec−2 or about a 16σ detection for Hα. c) Same as b) but for the log([O III]/Hβ) ratio. d) Representation of the Voronoi binning scheme adopted to recover the low-surface-brightness line emission in the DIG, as described in the text. Each bin is assigned a random colour. e) Map of log([S II]/Hα) using the binned data. f) Same as e) but for the log([O III]/Hβ) ratio.

Our binning choice in the DIG is motivated by the need to obtain a minimum average detection for the metal emission lines of interest. Assuming line ratios typical of the DIG of log([S II]/Hα) ∼ −0.5 (see also Sect. 5.2.2), and log ([O III]/Hβ) = −0.5 (see Sect. 5.3), we conclude that in order to achieve a ∼5σ detection for [O III] (and consequently ∼10σ for [S II]λ6717), one needs an Hα S/N of 60. The median Hα S/N of single spaxels in the DIG regime is around 9, implying we need to increase the S/N by a mean factor of seven. We achieve this goal by performing the following procedure. First, we make use of the segmentation map defining ionised nebulae, and treat each nebula as one bin. These regions are not binned further, in order to avoid mixing of flux between H II regions and DIG. Moreover, H II regions are intrinsically bright and need no further binning to ensure detection of the lines of interest. The DIG area outside the nebulae masks is re-gridded to spaxels of 1.4″ (i.e. 7 by 7 native spaxels), in order to increase the S/N by a factor of rougly seven. The Hα S/N in these larger pixels is calculated based on error propagation from the original spaxel-level Hα map, since the error vectors in each spaxel are nearly statistically independent (Weilbacher et al. 2020; Emsellem et al. 2022).

These square 1.4″ × 1.4″ bins are used as a starting point for a Voronoi binning procedure, following Cappellari & Copin (2003), which groups together the 1.4″ × 1.4″ pixels to reach a target Hα S/N of 60. This step is necessary to generate larger bins in regions of lower than average surface brightness in the DIG, and is especially important for galaxies that show large regions of faint line emission (e.g. NGC 1300, NGC 1433, NGC 1512). We could have applied the Voronoi binning to the original spaxels, but we found that S/N levels for individual spaxels are not always reliable, leading to bins with unnecessarily complicated small-scale (smaller than the PSF) morphologies.

We perform a new run of the spectral fitting algorithm to derive the stellar kinematics and emission-line properties on the new bins. The resulting maps for log([S II]/Hα) and log([O III]/Hβ) are shown in Fig. 1 (panels e and f). The figure demonstrates the power of our approach in recovering emission-line ratios in the low-surface-brightness regime. For the remainder of this paper we use these binned emission-line maps as a basis for further analysis of line emission from the DIG. This effectively lowers the spatial resolution of the DIG maps to the size of the average Voronoi bin (∼60 times the original 0.2″ MUSE spaxels), which corresponds to an average scale of 120 pc across our sample. We note, however, that the higher resolution of the original dataset is crucial to the identification of H II regions and their separation from the DIG.

3. General characteristics of the DIG

In this section, we discuss the general characteristics of DIG and H II regions obtained according to our masks derived in Sect. 2.3. In particular, we discuss key observables that have been used in the literature to distinguish DIG from H II regions: Σ_Hα (Zhang et al. 2017), EW(Hα) (e.g. Lacerda et al. 2018), and the location in the BPT diagram. The quantitative conclusions of this section are somewhat dependent on resolution, and will necessarily differ when galaxies are observed at lower spatial resolution (e.g. Poetrodjojo et al. 2019).

3.1. Hα emission in the DIG and H II regions

In our sample of galaxies, the DIG defined according to our masks shows variations of ∼1.2 dex in Hα surface brightness, log(Σ_Hα DIG/erg s⁻¹ kpc⁻²) = 38.2 ± 0.6 (where the range shows the 16th and 84th percentiles). In Fig. 2, we show histograms of the Σ_Hα distribution in the DIG and H II regions in our sample of galaxies, ordered by increasing stellar mass from left to right. We observe a wide variety in the resulting distributions. Because of our source extraction algorithm, the DIG is effectively defined in terms of its contrast in Σ_Hα with respect to clumpy emission, identified as H II regions. As a consequence, while the average Σ_Hα of DIG and H II regions varies across galaxies, the average contrast in Σ_Hα between DIG and H II region is rather stable at ∼1 dex. A few galaxies show tails of fainter DIG, in particular NGC 1087 and NGC 1385 show a tail of faint diffuse gas corresponding to the extreme galactic outskirts (galactocentric radii > 1 R₂₅, only covered in these two galaxies and in NGC 4254). For the high-mass galaxies NGC 1300, NGC 1512, and NGC 1433, on the other hand, the faintest regions in Hα correspond to the central regions of high mass surface density.

Fig. 2. ΣHα in H II regions (in blue) and DIG (in orange) in the PHANGS–MUSE sample ordered by increasing stellar mass from left to right. DIG and H II regions are defined according to a morphological criterion (using HIIPHOT; see Sect. 2.3). The 16th, 50th, and 84th percentiles of the distributions are shown with dashed lines on each violin. The blue and orange horizontal lines (and associated colour-shaded areas) represent the 50th (and the 16th and 84th) percentiles of the ΣHα distribution for H II regions and DIG across all galaxies.

We measure the fraction of the Hα flux associated with the DIG (f_DIG) by summing the Hα flux within the DIG masks and correcting this flux for dust extinction using the Balmer decrement. The resulting f_DIG is listed in Table 1, and ranges between 20% and 55%, with a median value of 37%. These numbers are in general agreement with previous estimates obtained from Hα narrow-band imaging (Oey et al. 2007), and also with the DIG fractions measured independently by Chevance et al. (2020) (median diffuse Hα fraction of 34%) for a subset of eight PHANGS–MUSE galaxies, following the Fourier-filtering methodology of Hygate et al. (2019). Our estimate depends on the choice of termination criterion for defining H II region masks. We investigated the effect of changing the termination gradient from our fiducial value to either steeper (resulting in smaller regions) or shallower gradients (resulting in larger regions). These cases are meant to reproduce the range of termination gradient values considered by the authors of HIIPHOT (Thilker et al. 2000) and bracket the range employed in the study of nearby galaxies by Oey et al. (2007)⁴. We find that the median value of f_DIG varies from 43% in the case of steep termination gradients to 22% for the shallowest one.

The brightest DIG emission is located near bright H II regions in all our galaxies, as can be qualitatively appreciated from Figs. A.1 and A.2. In particular, for galaxies with well-defined spiral arms, the DIG appears brighter in Hα near the arms than in the inter-arm regions. Hα emission falls off in intensity with distance from H II regions but otherwise shows a smooth morphology, except in the two galaxies with the lowest masses and smaller distances, NGC 5068 and IC 5332, where prominent shells and filaments are visible in diffuse Hα. Overall, the spatial correspondence of DIG and H II region strongly suggests that leaking radiation powers the DIG, an idea that we develop further in Sect. 4.

3.2. The DIG and EW(Hα)

Several studies have discussed the usefulness of EW(Hα) in distinguishing ionisation due to HOLMES from that resulting from star formation or active galactic nuclei (AGN; e.g. Stasińska et al. 2008; Cid Fernandes et al. 2011). The value of the EW(Hα) expected from HOLMES remains uncertain, even in the simple case where all the radiation is absorbed ‘locally’. Cid Fernandes et al. (2011) compared various population-synthesis codes and found overall agreement to within a factor of two for the number of ionising photons produced by HOLMES. This flux results in EW(Hα) = 1.5 − 2.5 Å for different assumptions regarding metallicity, initial mass function and stellar evolution tracks.

The number of ionising photons produced by HOLMES declines slowly as a function of age (by a factor of two from 100 Myr to 10 Gyr). The EW(Hα) they produce, however, increases with age (by a factor of five from 100 Myr to 10 Gyr) because of the decreasing brightness of the stellar continuum, according to the predictions from FSPS code (Conroy et al. 2009; Byler et al. 2019). In particular, a 10 Gyr old, solar-metallicity SSP generated with FSPS emits 5 × 10⁴⁰ ionising photons ${\mathrm{s}}^{-1}{M}_{\odot }^{-1}$, corresponding to an EW(Hα) ∼ 0.5 Å⁵. In contrast, the PEGASE population synthesis predicts a flux of 7 × 10⁴⁰ ionising photons ${\mathrm{s}}^{-1}{M}_{\odot }^{-1}$ for HOLMES, corresponding to EW(Hα) ∼ 1.1 Å. Given these uncertainties, we take the latter as our fiducial value in this work, as it lies somewhat in the middle of the range between the prediction from FSPS and the largest values quoted in Cid Fernandes et al. (2011).

In Fig. 3 we show histograms of log(EW(Hα)) of both DIG and H II regions for all the galaxies in the PHANGS–MUSE sample, ordered by stellar mass from left to right. The value of EW(Hα) = 3 Å, employed in the literature to distinguish H II regions from DIG, does not provide a clean division at the spatial resolution of our data. In fact, the median value of EW(Hα) in the DIG in our sample is $5.3_{-3.0}^{+12.1}$ Å. The value of 14 Å, suggested by Lacerda et al. (2018) to select pure H II regions, approximately corresponds to the 84th percentiles of the EW(Hα) distribution of the DIG, and the 16th percentile of the H II regions distribution, therefore providing a reasonable dividing value for our sample. The EW(Hα) distributions are smoother than those of Σ_Hα, simply as a result of the spatially smooth continuum.

Fig. 3. EW(Hα) in H II regions (blue) and DIG (orange) in the PHANGS–MUSE sample. DIG and H II regions are defined according to a morphological criterion (using HIIPHOT; see Sect. 2.3). The 16th, 50th, and 84th percentiles of the distributions are shown with dashed lines on each violin. Dashed horizontal black lines represent demarcations employed in the literature to distinguish H II regions from DIG (3 and 14 Å) and the value of EW(Hα) from HOLMES using the PEGASE models (1.1 Å). The blue and orange horizontal lines (and associated colour-shaded areas) represent the 50th (16th and 84th) percentiles of the EW(Hα) distribution for H II regions and DIG across all galaxies.

In most of the area we classify as DIG, the EW(Hα) is greater than our fiducial value of 1.1 Å predicted for pure HOLMES. In the central regions of some of the more massive galaxies in our sample, however, we encounter areas of high stellar mass surface density and EW(Hα) < 2 Å, where HOLMES are predicted to play a major role in keeping the gas ionised. This is particularly evident in, for example, NGC 3351, NGC 1433, NGC 1300, NGC 1512, and NGC 3627. In some of these galaxies (e.g. NGC 3627) there is evidence for an increase in EW(Hα) towards the centre, potentially associated with the AGN.

We show EW(Hα) maps for the full sample in Fig. 4, again ordered by stellar mass along rows from top to bottom. The maps reveal notable differences in the morphology and pervasiveness of regions with EW(Hα) < 3 (grey contour) and < 1.1 Å (white contour and white hashed for emphasis), respectively. Regions with EW(Hα) < 3 Å start to appear in inter-arm regions (e.g. NGC 7496, NGC 0628) of low-mass galaxies, but become particularly evident in the central regions of higher-mass galaxies swept by a bar (e.g. NGC 3351, NGC 1512, NGC 1433). These are also the only places where we observe extremely low EW(Hα) ∼ 1 Å, compatible with the EW(Hα) predictions for HOLMES at old ages.

Fig. 4. Maps of EW(Hα) for all the galaxies in the PHANGS–MUSE sample. Galaxies are ordered by increasing stellar mass from the top left to the bottom right. The contours correspond to EW(Hα) of 1.1 (white) and 3 Å (light grey). For added clarity, regions with EW(Hα) < 1.1 Å are shown with white hatches.

3.3. The DIG in the BPT diagram

In Fig. 5 we show maps of the galaxies in our sample colour-coded by the position of each region in the [S II] BPT diagram ([S II]/Hα versus [O III]/Hβ). Similar trends are seen in the [N II] ([N II]/Hα versus [O III]/Hβ) and [O I] ([O I]/Hα versus [O III]/Hβ) BPT diagrams. We restrict our discussion to the [S II] here because of its added ability to discriminate between AGN and LI(N)ER regions with respect to the [N II] diagram (Kewley et al. 2001; Belfiore et al. 2016; Law et al. 2021), and because the [S II] emission lines are stronger than [O I], therefore allowing classifications for a larger fraction of regions than the [O I] BPT. We perform only a simple classification into star-forming (blue), LI(N)ER (orange) and AGN (red) regions according to the Kewley et al. (2001, 2006) lines. Here we summarise four key points that are evident from these maps.

Fig. 5. Maps colour-coded by the position of each region in the [S II] BPT diagram (blue: star formation; orange: LI(N)ER; red: AGN). The classification is performed according to the demarcation lines of Kewley et al. (2001) (to classify star-forming regions) and the AGN-LI(N)ER demarcation line of Kewley et al. (2006). The black contours correspond to EW(Hα) of 1.1 and 3 Å (same as in Fig. 4).

First, in low-mass galaxies (e.g. NGC 5068, IC 5332, NGC 2835, NGC 1087, NGC 1672, NGC 628), the LI(N)ER-like DIG is found in inter-arm regions and galactic outskirts. These regions correspond to areas of low Σ_Hα, but not necessarily of low EW(Hα).

Second, four galaxies in our sample have AGN signatures in their nuclear regions. NGC 1672 and NGC 1512 show Seyfert-like ionisation in their central regions according to the [S II] BPT diagram (also confirmed in the [N II]/Hα versus [O III]/Hβ plane, and for NGC 1672 by independent evidence from X-ray observations; Jenkins et al. 2011). NGC 1566 and NGC 1365 host Type 1 (or intermediate type; see Oknyansky et al. 2020) AGN, with broad hydrogen Balmer lines. In NGC 1672, a small region to the south-west of the nucleus shows Seyfert-like line ratios and can be associated with an AGN outflow. NGC 1365 hosts a well-studied biconical AGN ionisation cone (Storchi-Bergmann & Bonatto 1991; Kristen et al. 1997; Veilleux et al. 2003; Venturi et al. 2018) to the south-east and north-west of the nucleus, evident as a region of Seyfert-like ionisation in the [S II] (and all other) BPT diagrams.

Third, NGC 3351, NGC 1300, NGC 1512, and NGC 1433 host a circum-nuclear star-forming ring, surrounded by a ‘star-formation desert’, with LI(N)ER and Seyfert-like line ratios. These star-formation deserts have very low EW(Hα) ∼ 1 − 2 Å, supporting the idea that their ionisation is largely powered by HOLMES.

Finally, NGC 1566 and NGC 3627 show the signature of kiloparsec-scale, spatially extended LI(N)ER emission in their central regions (equivalent to the typical ‘central LIER’ galaxies observed in kiloparsec-resolution data of nearby galaxies; Belfiore et al. 2016). These central regions also show low EW(Hα) < 3 Å. This combined evidence supports the idea that line emission in the central regions of these galaxies appears LI(N)ER-like because of the substantial contribution of HOLMES to the local ionisation budget, despite the presence of a Seyfert 1 nucleus in NGC 1566⁶.

In summary, six of our 19 galaxies show extended LI(N)ER regions, four of them in combination with circum-nuclear rings. While the maps in Figs. 4 and 5 show large areas with EW(Hα) < 3 Å and LI(N)ER-like ionisation, respectively, these regions contribute a negligible amount of the total Hα flux in the galaxy. In particular, on average across our galaxies 88% of the flux in our DIG masks is contributed by regions that fall in the star-forming section of the [S II] BPT, while these same regions only occupy 58% of the DIG on-sky area. While this is clearly different from galaxy to galaxy, it highlights the fact that the DIG with extreme LI(N)ER-like ratio contributes only a small amount of the total Hα flux, while occupying a significant fraction of the area.

4. Spatial model for the propagation of leaking radiation from H II regions

4.1. Thin slab model

Having determined that HOLMES are subdominant in powering the DIG on galaxy-wide scales, we proceed here to test the association between DIG Hα emission and H II regions by building a toy model for the propagation of ionising radiation in the ISM.

We consider a thin slab model for the propagation of ionising radiation in the DIG, as originally suggested by Miller & Cox (1993), and subsequently explored by Zurita et al. (2002) and Seon (2009). The model described by Seon (2009) is identical to the one presented in this work. We summarise its derivation in Appendix B. In short, the model considers propagation of ionising photons leaking from H II regions due to both geometric dilution and absorption by a population of DIG clouds. The latter effect is parametrised via an effective absorption coefficient, k₀. The model predicts that the observed Hα surface brightness Σ_Hα, i of pixel i is given by

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{H}\alpha ,i}={\displaystyle \sum _j{f}_{\mathrm{scale}}\mathrm{\Delta }h{F}_{\mathrm{H}\alpha ,j}^{\mathrm{H}\mathrm{II}}\frac{{\mathrm{e}}^{-k_{0}r}}{4\pi r^2}\frac{1-{\mathrm{e}}^{-k_0\mathrm{\Delta }r}}{\mathrm{\Delta }r},}\end{array}$$(1)

where $F_{\mathrm{H}\alpha ,j}^{\mathrm{H}\mathrm{I}\mathrm{I}}$ is the Hα flux of catalogued H II regions, and the scaling factor, f_scale = f_esc/(1 − f_esc), where f_esc is the escape fraction of ionising photons from individual H II regions. Δr is the linear size of the pixel considered and Δh the height of the thin slab, and k₀ is the effective absorption coefficient (such that the mean free path of ionisation radiation λ₀ ≡ 1/k₀).

In order to compare the data with this model, we first subtracted the Hα flux due to HOLMES, assuming our fiducial value of the HOLMES ionising photon flux (7 × 10⁴⁰ ionising photons s${}^{-1}{M}_{\odot }^{-1}$) discussed in Sect. 3.2, and calculating the ionising flux expected at each location based on the stellar mass surface density of each pixel. This is only an approximate estimate, since it does not account for changes in the ionising flux from HOLMES as a function of age of the stellar population. Nor does it model the propagation of ionising photons from HOLMES. Since we performed a detailed spectral fitting of the MUSE data, it is in principle possible to calculate the HOLMES flux for the specific age and metallicity of the underlying stellar population at each position, but the degeneracies involved in the spectral fitting and uncertainties in model predictions for the post-AGB phase do not currently warrant this level of sophistication. We checked the effect of decreasing or increasing the HOLMES photon flux by a factor of two, and found that this had negligible impact on the conclusions presented in this section. In fact, not performing the subtraction at all only had a small overall impact, causing only marginal shifts in the value of the best-fit parameters.

For each galaxy we assumed k₀ to be constant across the disc and generated the model Σ_Hα via the following sequence of steps. First, we considered the Hα flux corrected for extinction (via the Balmer decrement) from each H II region in our catalogue as $F_{\mathrm{H}\alpha ,j}^{\mathrm{H}\mathrm{I}\mathrm{I}}$. The flux for each region was taken to be simply the sum of the flux within the region mask. Sources in the H II region catalogue that do not fall within the star-forming demarcation lines of Kewley et al. (2001) in the BPT diagrams of [N II]/Hα versus [O III]/Hβ, [S II]/Hα versus [O III]/Hβ and [O I]/Hα versus [O III]/Hβ were excluded from this computation since such sources are likely not to be bona fide H II regions, but rather PNe, supernova remnants, or clouds photoionised by the AGN. These sources will clearly also produce ionising radiation, but we excluded them here for the sake of simplicity. Their inclusion does not change any of the results since they contribute only a small fraction (generally a few percent) of the total Hα flux in each galaxy. If Hβ is not detected within the H II region (this only happens for 0.2% of regions), we assumed no dust extinction. If we do not detect all the metal lines necessary for the BPT classification, we assumed the source is a faint H II region and retain it in the sample (this only affects 2% of the regions). Different choices for these edge cases are inconsequential to the subsequent results.

For each H II region, we considered its Hα luminosity-weighted centre (using the geometric centre does not change any of our conclusions) and computed the deprojected distance to other pixels in the galaxy, taking the position angle and inclination into account. We only considered radiation propagating along the midplane of the galaxy. We used the deprojected distance, r_ij, to compute the Hα emission in the DIG according to Eq. (1). In particular, since f_scale and Δh are degenerate with each other, we computed the radiation field described in Eq. (1) up to a multiplicative factor.

The multiplicative scaling factor was subsequently determined by comparing the model with the data (Sect. 4.2).

We did not correct the DIG Hα emission for dust extinction, therefore effectively assuming no extinction in the diffuse medium. This assumption is not strictly correct as we measure an average E(B − V) = 0.15 mag in the DIG, based on the Balmer decrement. We also did not treat the effect of dust in the DIG on the propagation of ionising photons. Radiative transfer models from Pellegrini et al. (2020a) demonstrate that this is a reasonable assumption for low-inclination galaxies.

Finally, we convolved the model with the median PSF of the MUSE mosaic and spatially binned it using the same bins as in the observed Hα map.

4.2. Comparing the observed Hα with a leakage model

Models were computed for nine logarithmically spaced values of k₀ = [0.1, 10.0] kpc⁻¹. We also computed a model with k₀ = 10⁻⁴ ∼ 0 kpc⁻¹, which corresponds to the case of spherical dilution, since in such a model the mean free path of ionisation radiation λ₀ = 10⁴ kpc is much larger than the physical size of the system (∼10 kpc).

For each galaxy in our sample and each value of k₀ we compare the model predictions with the observed Hα map for the areas classified as DIG (H II regions are masked). Since we only computed the model Hα distribution up to a multiplicative scale factor, we need to rescale the model to the observations. We do so by fitting a straight line of the form

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{H}\alpha ,\mathrm{observed},i}=a_0+a_1{\displaystyle \sum _j{F}_{\mathrm{H}\alpha ,j}^{\mathrm{H}\mathrm{II}}\frac{{\mathrm{e}}^{-k_0{r}_{\mathit{ij}}}}{4\pi r_{\mathit{ij}}^2}·}\end{array}$$(2)

According to Eq. (1), the slope parameter a₁ is proportional to f_scale Δh, while a₀ should be zero if the model can correctly predict the DIG observations without any additional background.

In Figs. 6 and 7 we illustrate the results of this procedure for the two example galaxies, NGC 4535 and NGC 4254. For each galaxy we present the observed Σ_Hα map of the DIG (H II regions are masked and appear white) and the model predictions for different values of k₀ (top row), maps of the logarithmic difference between model and data (middle row) and examples of the median relations between log(best-fit model/data) and Σ_{Hα,  observed}. The ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ for each model is also indicated, which is computed by considering an intrinsic scatter of 10%, added in quadrature to the observed flux errors.

Fig. 6. Example of our DIG modelling and its comparison with the observed ΣHα for NGC 4535. Top row: observed Hα flux map (H II regions are masked and appear white in all the maps) and the model predictions for different values of k0, generated according to Eq. (1). Middle row: a spatial comparison of the model and observations, in the form of log(Model/Data). Bottom row: median relation between log(Model/Data) and ΣHα,  observed with associated scatter (the shaded region corresponds to the range between the 16th and 84th percentiles). The ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ for each model is also listed above each panel. The model with k0 ∼ 1 kpc−1 offers a better representation of the data with respect to models with higher (k0 = 10 kpc−1) or lower values (k0 = 0 kpc−1, i.e. only spherical dilution).

Fig. 7. Same as Fig. 6 but for NGC 4254.

Both galaxies manifest similar trends in this comparison. Considering the model with k₀ = 10 kpc⁻¹, the radiation field does not travel far enough from the H II regions. The observed Σ_Hα is higher than that predicted by the model in low Σ_Hα regions, explaining the large areas of negative log(model/data). At the other extreme, in the case of spherical dilution (k₀ = 0), the model predicts too much diffuse flux with respect to the data. For both example galaxies, the spherical dilution model overestimates the Hα flux at low Σ_Hα and underestimates it in the vicinity of H II regions (high Σ_Hα).

For k₀ = 1 kpc⁻¹, both NGC 4535 and NGC 4254 exhibit a nearly flat relation between log(model/data) and Σ_{Hα,  observed}. The model predicts roughly the correct amount of Hα both in the vicinity of H II regions and in the diffuse component. While such models still show a relatively high ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$, we consider the fit satisfactory given the simplistic assumptions we have employed.

4.3. Best-fit model parameters

We present the variations in the ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ of the fit as a function of k₀ for each galaxy in Fig. 8a. A star symbol denotes the position of the model with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$. The figure demonstrates that for all galaxies we observe a well-defined ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ minimum. For the example galaxies in Figs. 6 and 7, NGC 4535 and NGC 4254, the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ corresponds to log(k₀/kpc⁻¹) = 0.0 and log(k₀/kpc⁻¹) = − 0.25, respectively, or mean free paths of 1.0 and 1.8 kpc. Across the full sample, the average mean free path of the best-fit model is 1.9 kpc (k₀ = 0.52 kpc⁻¹). This scale can be compared with the mean nearest neighbour distance between H II regions in our catalogue (65 − 300 pc, median 170 pc), and that measured by Chevance et al. (2020, 80 − 190 pc, median 180 pc) for a subsample of eight of our targets, showing that the DIG emission is truly diffuse and occurs on scales much larger then the separation scales of the units of the star-formation cycle⁷. We do not find any significant correlations between the best-fit value of k₀ and the stellar mass, gas fraction, distance or inclination of the galaxies in our small sample.

Fig. 8. Summary of best-fit model parameters for the DIG leakage model applied to our entire sample. a) Variation in ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ for the DIG spatial models described in the text as a function of the absorption coefficient, k0, for all the galaxies in our sample. Low values of k0 correspond to the case of spherical dilution, while high values correspond to exponential attenuation with scale length 1/k0. For each galaxy, a coloured star denotes the model with minimum ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$. b) Value of the additive offset, a0, defined in Eq. (2), as a function of k0. According to the simple model described in the text, a0 ∼ 0. Star symbols correspond to the k0 values with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$, which are clustered around a0 ∼ 0. c) fesc, the escape fraction of ionising radiation from H II regions as a function of k0 for each galaxy, assuming Δh = 1 kpc. Star symbols correspond to the k0 values with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$.

Figure 8b presents the variation in the additive offset a₀ as defined in Eq. (2). This parameter should be equal to zero, if the DIG is well described by our modelling framework. We find in Fig. 8b that a₀ is close to zero for the models with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ (coloured stars). The value of a₀ increases from negative to positive with increasing k₀, implying that there is a value of k₀ for which a₀ = 0 exactly (i.e. it is suitable to reproduce the DIG without any constant background). For NGC 3351 and NGC 1365, such a change in sign of a₀ is not observed. Both galaxies have extended regions (the central regions in NGC 3351 and an AGN ionisation cone in NGC 1365; see Sect. 5.6) that are not ionised by leaking radiation from H II regions.

Finally, the multiplicative offset a₁, also defined in Eq. (2), is related to the product f_scale Δh. Measurements of the DIG scale height in our own galaxy or in edge-on spiral galaxies span a range of 0.5 − 1.5 kpc (Reynolds 1991; Ferrière 2001; Gaensler et al. 2008; Levy et al. 2019). Here we take the median value of 1 kpc obtained by Levy et al. (2019) using a compilation of literature measurements. The thickness Δh is given by twice the scale height, but it is reasonable to assume that the Hα emission from the far side of the galaxy with respect to the disc midplane is obscured, and therefore not observed. Assuming Δh = 1 kpc and taking into account a correction for galaxy inclination, we can derive an estimate of the escape fraction f_esc. We show the derived value of f_esc for each galaxy as a function of k₀ in Fig. 8c. For each galaxy the model with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ is shown with a coloured star. As expected, for most galaxies a smaller k₀ (and therefore a larger mean free path) requires a larger escape fraction to fit the data. For the best-fit models, the median escape fraction is 40%, with IC 5332 being an outlier with f_esc > 0.6. These values of escape fraction are consistent with those obtained in Sect. 3.1, if one assumes f_DIG = f_esc.

Considering the assumed thickness of the DIG layer, the toy model for a thin disc presented in Sect. 4 may not be a good approximation since the emitting region is not particularly thin, even though the galaxy as a whole is. Seon (2009) compared the thin-slab model also used in this work with a 3D Monte Carlo radiative transfer code, confirming that they lead to similar results for Δh/2 < 1/k₀.

In summary, the simple model for the propagation of ionisation radiation that we have developed in this section is capable of reproducing reasonably well the observed Σ_Hα in the DIG for a median value of the mean free path of ionising radiation of 1.9 kpc. While our models do not reproduce the small-scale variations in Σ_Hα observed in the data, as reflected by the high best-fit ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$, they are nonetheless capable of reproducing the average observed Σ_Hα in the DIG from the lowest to the highest observed Hα surface brightness. Our conclusions are comparable to those drawn by Zurita et al. (2002) and Seon (2009) for NGC 157 and M 51, respectively.

The spatial model for the DIG we have developed potentially allows us to calculate the fraction of the Hα flux within each H II region spatial mask due to radiation leaked from all other regions in the galaxy. Subtracting this DIG contamination from the Hα flux provided in our H II region catalogue would allow for a more accurate estimate of the intrinsic emission of ionising photons for each H II region (however, the contribution from leaking radiation along the line of sight from the region itself would have to be estimated independently). A change in the intrinsic H II region flux would, in turn, result in a change in the DIG model, requiring a recursive approach to achieve convergence. Given the uncertainties in the modelling framework, we do not attempt to follow this approach here. We note, however, that the main effect of these corrections would be to lower the luminosity of H II regions in areas of high spatial clustering (e.g. spiral arms) and therefore larger DIG surface brightness. This may potentially reduce the discrepancies observed in Figs. 6 and 7, where the best-fit models show too much flux in the vicinity of the spiral arms.

The conclusion from this modelling effort is that the DIG Σ_Hα is consistent with photoionisation from LyC radiation leaking from H II regions. Since successful models can be built with k₀ being constant throughout the disc, the clumping (or volume filling) factor of photoionised clouds in the DIG is roughly constant within galaxies, and does not change much across our sample. These findings motivate the study of line ratios in the DIG as a function of Σ_Hα presented in the next section.

5. Line ratios in the DIG

5.1. Dependence of line ratios on Σ_Hα and metallicity

According to the analysis in the previous section, in the DIG the distance between the source of the ionising photons and the absorbing medium is much larger than the typical radii of H II regions. In order to predict the line ratios in a photoionised medium one needs to consider the ratio between the photon number density and the electron density, referred to as the ionisation parameter, U₀, defined as

$$\begin{array}{c}\hfill U_0=\frac{{\mathrm{\Phi }}_{\mathrm{i}}}{n_{\mathrm{e}}c},\end{array}$$(3)

where Φ_i is the number of ionising photons reaching the incident surface of the ionised cloud per unit time and unit area⁸. In the low electron density limit, which always applies to the DIG, the thermal and ionisation structure of a nebula depends only on U₀, and not on Φ_i and n_e individually.

The electron density in the DIG cannot be directly measured with our data. Typical density-sensitive line ratios in the optical wavelength range, such as the [S II]λ6717/[S II]λ6731 ratio, saturate at their low-density limit for n_e less than a few tens of cm⁻³ and are therefore unsuitable to measure the density of the DIG. Studies of the diffuse emission in the Milky Way have determined that the typical electron density of DIG clouds is n_e ∼ 0.1 cm⁻³ (Reynolds 1991; Ferrière 2001), about 2 − 3 dex lower than typical values for H II regions (Kennicutt 1984; Dale et al. 2009; Kewley et al. 2019). Since the typical sizes of H II regions (10 − 100 pc) are 1 − 2 dex smaller than the mean free path of ionising photons in the DIG as determined in Sect. 4.3, we predict the ionisation parameter in the DIG to be up to ∼2 dex lower than in H II regions. Typical values for the ionisation parameter in H II regions lie in the range log(U₀) = [−3.5, −2.5] (Kewley et al. 2019; Mingozzi et al. 2020). We therefore expect the DIG to have log(U₀) = [−4.5, −3.5]. If the electron density is roughly constant within the DIG layer, Σ_Hα is proportional to the ionisation parameter, since we have demonstrated in Sect. 4.1 that Σ_Hα traces the diluted ionising radiation field.

The interpretation of line ratios in the DIG, especially in external galaxies, is complicated by the difficulty of isolating the effects of changes in ionisation parameter from variations in both metallicity and the shape of the ionising continuum. In light of the discussion in the previous section, we assume here that the DIG is mostly ionised by radiation from young stars leaking from H II regions. We discuss in Sect. 5.4 the effect of the harder radiation field from HOLMES.

Having assumed a fixed input spectrum, we marginalise over the effects of metallicity on the line ratios in the DIG by looking for trends in bins of galactocentric radius. Since galaxies present well-defined metallicity gradients (Moustakas et al. 2010; Sánchez et al. 2014; Ho et al. 2015; Belfiore et al. 2017; Mingozzi et al. 2020), and azimuthal variations in chemical abundances in our galaxies are small (< 0.05 dex; Kreckel et al. 2019, 2020), we assume that within our radial bins (described in Sect. 5.2.2) the DIG is roughly chemically homogeneous.

We study the changes of line ratios in the DIG following the approach of Zhang et al. (2017), who studied the line ratios in the DIG as a function of both Σ_Hα and galactocentric distance using data from the MaNGA survey. Such an approach allows some separation of the effect of ionisation parameters, via Σ_Hα as its proxy, from metallicity, assumed to decrease with galactocentric distance. Unlike Zhang et al. (2017), however, we resolve the mean separation scale between H II regions and can therefore measure the DIG line ratios without being contaminated by flux from H II regions.

5.2. Low-ionisation line ratios in the DIG

5.2.1. Model predictions for low-ionisation lines

Before describing the data, it is helpful to consider the expected trends in the scenario where the DIG consists of low-U₀ radiation-bounded clouds, subject to the diluted ionising spectrum of young stellar populations. Under these conditions the DIG is effectively described by the same model ingredients as canonical H II regions, albeit with lower log U₀. In Fig. 9, we show the predictions for the [S II]/Hα, [N II]/Hα and [O I]/Hα line ratios from gas photoionised by a 2 Myr SSP generated with the flexible stellar population synthesis (FSPS) code. As discussed in Byler et al. (2017), the line ratios obtained from a 2 Myr old ‘burst’ model are roughly equivalent to those expected from a model with a constant star-formation history, once it has reached a steady state. Moreover, for the adopted MIST isochrones the line ratios do not change substantially for ages between 1 and 4 Myr (see Fig. 24 of Byler et al. 2017).

Fig. 9. Model predictions for changes in low-ionisation line ratios with ionisation parameter in the case of leaking radiation. Models are computed with CLOUDY v17.02 and a 2 Myr SSP spectrum generated with FSPS as input. The line ratios, from left to right log([S II]/Hα), log([N II]/Hα), and log([O I]/Hα), are plotted as a function of the ionisation parameter at the incident surface of the cloud (log(U0)). Different colours correspond to different gas-phase metallicities, as indicated in the legend. The light grey band corresponds to the 10th and the 90th percentiles of the observed line ratios in the DIG across all the galaxies in our sample. We further highlight in dark grey the area of log(U0) = [−4.5, −3.5], corresponding to the expected range of ionisation parameter values in the DIG. The models span the range of observed line ratios in the DIG for the expected values of log(U0) in the case of log([S II]/Hα) and log([O I]/Hα), while they fall somewhat short of reproducing the highest line ratios for log([N II]/Hα).

The photoionisation model calculations are run using the CLOUDY photoionisation code v17.02, last described by Ferland et al. (2017), and more details on the input parameters are given in Appendix C. For each line ratio we show its dependence on the ionisation parameter log U₀, for different values of oxygen abundance, 12 + log(O/H) = [8.09, 9.09] in intervals of 0.2 dex, where the solar oxygen abundance is taken to be 12 + log(O/H) = 8.69 (Asplund et al. 2009). The absolute scale of ISM metallicity is notoriously difficult to determine accurately from H II region line ratios (see e.g. Blanc et al. 2015; Kewley et al. 2019). Making use of the Pilyugin & Grebel (2016) metallicity calibration, which is found to be in good agreement with abundances measured via the auroral line temperatures, the H II regions in PHANGS–MUSE galaxies span the metallicity range 12 + log(O/H) = [8.3, 8.6] (Kreckel et al. 2019). Calibrations that directly compare observed line ratios with photoionisation models, however, predict considerably higher abundances (by ∼0.2 − 0.4 dex), especially at the high-metallicity end. Considering these systematic uncertainties, we computed models for a range in metallicity wide enough to span the entire range of plausible abundances present in our galaxies.

[S II]/Hα, [N II]/Hα and [O I]/Hα are predicted to generally increase with decreasing log U₀ (i.e. moving towards lower Σ_Hα), although the trends can be seen to flatten or mildly invert for the lowest ionisation parameters and highest-metallicity models. For both [S II]/Hα and [O I]/Hα the highest values are associated with near-solar metallicity, while higher and lower metallicities correspond to lower line ratios. [N II]/Hα, in contrast, increases monotonically with metallicity, as expected given its use as a metallicity diagnostic. In Fig. 9 we show as a light grey horizontal band the range of observed line ratios (from the 10th to the 90th percentile of the distribution) in the DIG of our sample. The dark grey box corresponds to the intersection with the expected log U₀ in the DIG (log U₀ = [ − 4.5, −3.5]).

The absolute values of the line ratios observed in our data also agree, to first order, with those predicted by the models (as seen by comparing the models with the grey boxes in Fig. 9). The match is not perfect, however: our models do not reach the highest values of [N II]/Hα observed in the DIG at the lowest Σ_Hα and in the range of ionisation parameters expected in the DIG, the [N II]/Hα ratio does not increase, but seems roughly flat. We do not consider these discrepancies large enough to discredit the assumption that the DIG is ionised by leaked radiation from H II regions, although such a discrepancy has been used in the past to motivate the need for an additional shock component to the DIG (Haffner et al. 2009).

The models presented here do not include changes in the N/O abundance ratio, age of the ionising clusters, or hardening of the ionisation spectrum due to intervening absorption, effects that would increase the values of the considered line ratios and potentially lead to a better match with the data. Hardening of the spectrum due to intervening absorption will cause an increase in [N II]/Hα and [S II]/Hα with respect to the models shown in Fig. 9. A change in N/O will manifest itself mostly in an increase in the flux of the [N II] line and a higher [N II]/Hα for all values of the ionisation parameter. We discuss the effect of adding a radiation source with a harder spectrum in Sect. 5.4.

5.2.2. Observed trends for low-ionisation lines in the DIG

In Fig. 10 we show log([S II]/Hα) as a function of Σ_Hα for regions included within the DIG mask for all the galaxies in our sample. Each galaxy is sub-divided into five radial bins, going from the galaxy centre (red) to the outskirts (purple) in steps of 0.2 R_max, where R_max is a measure of the maximum radial coverage for each galaxy. On average R_max ∼ 0.6 R₂₅ or ∼8 kpc for our sample, meaning each radial bin is ∼1.5 kpc in width. Each galaxy has a different metallicity gradient (Kreckel et al. 2019), so stepping in relative radius aims to qualitatively capture radial variations in metallicity. The line ratios are not corrected for extinction in the DIG. Such a correction would have a small impact since we only consider ratios of lines close in wavelength.

Fig. 10. Log([S II]/Hα) in the DIG as a function of log(${\mathrm{\Sigma }}_{\mathrm{H}\alpha }/{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{kpc}}^{-2}$) for different radial bins in individual galaxies in the PHANGS–MUSE sample. Galaxies are ordered by stellar mass from lowest (top left) to highest (bottom right). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, and NGC 1433) have higher values of this (and other) line ratios in their central radial bin (in red). Aside from this feature, the galaxies in the sample follow a remarkably tight relation between log([S II]/Hα) and ΣHα.

The relation between [S II]/Hα and Σ_Hα is remarkably similar among different radial bins and across different galaxies. At low Σ_Hα (i.e. log(Σ_Hα/erg s⁻¹ kpc⁻²) < 38) we find log([S II]/Hα) ∼ 0.0, while at the highest Σ_Hα [log(Σ_Hα/erg s⁻¹ kpc⁻²) ∼ 40] we find log([S II]/Hα) ∼ −0.5. Variations in [S II]/Hα with radius within individual galaxies are small and are most evident at high Σ_Hα, especially for NGC 628 and NGC 4254. These variations go in the direction of higher-metallicity inner regions having lower [S II]/Hα at fixed Σ_Hα, which is the trend expected from the photoionisation models, assuming these massive galaxies have super-solar metallicities in their central regions.

The most evident deviations from the general trend are seen in a set of five galaxies, NGC 3351, NGC 1300, NGC 1512, NGC 3627, and NGC 1433, which show enhanced line ratios at given Σ_Hα in their central regions. This behaviour is seen also in [N II]/Hα and [O I]/Hα and even more clearly in the high-ionisation [O III]/Hβ line ratio discussed in the next section. We therefore postpone a discussion of the origin of this trend to Sect. 5.4.

In Figs. 11 and 12 we show the changes in [N II]/Hα and [O I]/Hα with Σ_Hα, presented in the same format as in Fig. 10. [N II]/Hα reveals its metallicity dependence through the vertical offset between radial bins, with the innermost radius being highest in [N II]/Hα (particularly evident in e.g. NGC 2835). Such patterns demonstrate the existence of a metallicity gradient, such that the metal-rich inner regions have higher [N II]/Hα in their DIG. In the case of [O I]/Hα all radial bins, except for the inner regions of massive galaxies, follow a similar trend and no strong metallicity dependence is observable. This is also in accordance with model predictions, which show that [O I]/Hα is only marginally affected by metallicity except for 12 + log(O/H) > 8.9.

Fig. 11. Same as Fig. 10 but for log([N II]/Hα). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) have exceptionally high values of this (and other) line ratios in their central radial bin (in red). Aside from this feature, the galaxies in the sample follow an average relation between log([N II]/Hα) and ΣHα, with a clear secondary dependence on metallicity. Inner regions have higher metallicity and show elevated [N II]/Hα at fixed ΣHα.

Fig. 12. Same as Fig. 10 but for log([O I]/Hα). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) show exceptionally high line ratio values in their central regions (in red).

5.3. High-ionisation line ratios in the DIG: [O III]/Hβ

Literature studies of the change of [O III]/Hβ with Σ_Hα in the DIG present conflicting results, with some authors finding little or no change (for the Milky Way, Haffner et al. 2009; for M 31, Greenawalt et al. 1998) and others finding an increase in the DIG (with scale height in edge-on galaxies, e.g. Rand 1998; Otte et al. 2001; Collins & Rand 2001, or with decreasing Σ_Hα towards the inner regions of the Milky Way, e.g. Madsen & Reynolds 2005). Either trend contradicts models for clouds photoionised by massive stars, since such models predict [O III]/Hβ to decrease with decreasing log U₀. This discrepancy has led to proposals of additional ionisation sources in the DIG, for example shocks (Collins & Rand 2001). An important piece of evidence was uncovered by the analysis of Zhang et al. (2017) based on MaNGA data. They found that the change of [O III]/Hβ with Σ_Hα depends on stellar mass. Low-mass galaxies show a marginal decrease in [O III]/Hβ, while massive galaxies show a clear increase in [O III]/Hβ with decreasing Σ_Hα.

In Fig. 13 we reassess these trends using the PHANGS–MUSE sample. Following the same formalism as in Figs. 10–12, we show the change of [O III]/Hβ with Σ_Hα for different radial bins in individual galaxies. As previously noted for low-ionisation line ratios in Sect. 5.2.2, the inner regions of several massive galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) have noticeably higher [O III]/Hβ at fixed Σ_Hα than all other data. In addition, we find systematic changes of the observed [O III]/Hβ versus Σ_Hα trends with stellar mass. Lower-mass galaxies show a roughly constant [O III]/Hβ as a function of Σ_Hα (e.g. NGC 5068, IC 5332, NGC 1087, and NGC 1385). For higher-mass galaxies, in contrast, we find a decrease in [O III]/Hβ as a function of Σ_Hα of up to 1 dex from log(Σ_Hα/erg s⁻¹ kpc⁻²) = 38 − 40. These trends agree with the results of Zhang et al. (2017), obtained using kiloparsec-resolution MaNGA data.

Fig. 13. Log([O III]/Hβ) in the DIG as a function of log(ΣHα/erg s−1 kpc−2) for different radial bins for individual galaxies in the PHANGS–MUSE sample. Galaxies are ordered by stellar mass from lowest (top left) to highest (bottom right). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) show exceptionally high line ratio values in their central regions (in red). Galaxies present a systematic change in the slope of the log([O III]/Hβ) versus ΣHα relation, with low-mass galaxies showing comparable line ratios at low and high ΣHα and high-mass galaxies showing a clear decrease in log([O III]/Hβ) with ΣHα. Both of these trends are at odds with predictions from photoionisation models.

5.4. Model predictions for HOLMES

In the left panel of Fig. 14 we show the predicted changes in [O III]/Hβ for clouds photoionised by a 2 Myr SSP (solid lines) as a function of log U₀ (same model parameters as Fig. 9). [O III]/Hβ is predicted to increase with log U₀, at odds with our observations where this ratio remains constant or decreases with Σ_Hα. The light grey horizontal band corresponds to the range between the 10th and the 90th percentile of the observed [O III]/Hβ in the DIG, and only overlaps with the models at log U₀ ∼ −3.5. The models, however, do not cover the data for the range of expected values of log U₀ in the DIG (dark grey box). In fact, the only similarity between models and data can be found at the high Σ_Hα end, where more metal-rich, higher-mass galaxies show lower [O III]/Hβ, as predicted by the models. An additional ionisation source, which would contribute little to the Balmer line emission, but generate a higher [O III]/Hβ ratio, is therefore needed to explain the data.

Fig. 14. Model predictions for the [O III]/Hβ line ratio considering both young stars and HOLMES. a) Model predictions for photoionised clouds, computed with CLOUDY and a 2 Myr (solid lines) and a 10 Gyr (dashed lines) SSP spectrum. The line ratio log([O III]/Hβ) is plotted as a function of the ionisation parameter at the incident face of the cloud (log(U0)). The vertical extent of the grey area highlights the range of line ratios observed in the DIG in the PHANGS–MUSE data. Different colours correspond to different gas-phase metallicities, with the same colour-coding in both panels, as indicated in the legend. All models predict log([O III]/Hβ) to increase with ΣHα, in contrast to what is observed. b) Model predictions for the log([O III]/Hβ) line ratio as a function of the fraction of the ionising radiation originating from old stars, fion = ΦOB/ΦHOLMES. We show models for a fixed value of log(U0) = − 3.5 and different metallicities. The horizontal light grey band corresponds to the range of values of log([O III]/Hβ) observed in the DIG of our sample, while the dark grey box corresponds to the range of fion observed in the PHANGS data. fion is estimated by computing the ionising photons expected to be produced by HOLMES given the local stellar mass surface density.

Several authors have pointed out that there is a natural candidate for a pervasive ionisation source with the necessary hard spectrum: HOLMES (Binette et al. 1994; Stasińska et al. 2008; Flores-Fajardo et al. 2011; Sarzi et al. 2010; Yan & Blanton 2012; Belfiore et al. 2016; Byler et al. 2019). Byler et al. (2019) presented photoionisation models for an old stellar population based on MIST stellar evolution tracks, which provide a consistent framework for the evolution of low- and intermediate-mass stars from the main sequence to the white dwarf cooling sequence and are therefore very well-suited to the study of the ionising spectrum of HOLMES. We build on the work of Byler et al. (2019) and generate a 10 Gyr spectrum using FSPS with MIST isochrones. This spectrum is used to generate a grid of CLOUDY models, covering the same values of gas phase metallicity and ionisation parameter as our 2 Myr models. A full description of the models is presented in Appendix C.

The [O III]/Hβ sequence produced by the 10 Gyr SSP models are shown as dashed lines in Fig. 14a. At solar metallicity and above, the increase in [O III]/Hβ with respect to the 2 Myr models is > 0.5 dex, and gets as large as 1.5 dex at 12 + log(O/H) = 9.09. The old star models are able to produce sufficiently high [O III]/Hβ at low log U₀, as evidenced by the fact that their position now overlaps with the dark grey box (expected position of the data) in Fig. 14a. Since the [O III]/Hβ ratios predicted by young star models are very low, even a small fraction (a few percent) of ionising flux from HOLMES can dominate the [O III] emission at low Σ_Hα.

5.5. DIG as a mixing sequence: From leaking radiation to HOLMES

We next explore the possibility that the DIG is ionised by a combination of ionising photons from both young stars (which we represent with the 2 Myr SSP model) and HOLMES (represented by the 10 Gyr model). We build CLOUDY models with variable ratios of the ionising flux from OB stars with respect to that from HOLMES. We parametrise the models via f_ion = Φ_OB/Φ_HOLMES, where we take Φ_OB and Φ_HOLMES are the ionising fluxes from the 2 Myr and 10 Gyr models, respectively. We compute models with log(f_ion) = [−3, −2, −1, 0, 1, 2]. In order to demonstrate the effect of this mixing sequence on [O III]/Hβ, we pick models with log U₀ = −3.5 (which we assume typical of the DIG, as discussed in Sect. 5.1) and present the change of the line ratio as a function of f_ion for different metallicities in Fig. 14b.

For log(f_ion) > 2, one recovers the line ratios predicted by the pure 2 Myr models. A contribution of HOLMES larger than 1%, however, leads to noticeable changes in the line ratios. Namely, for −1 < log(f_ion) < 2, the predicted line ratios are intermediate between those of the 2 Myr and 10 Gyr populations, while for lower and higher values they are indistinguishable from those generated by the two extremes. [O III]/Hβ increases with an increasing contribution of HOLMES to the ionising flux, and the changes in the line ratio are largest for the high-metallicity models. At 12 + log(O/H) = 9.09, the difference in [O III]/Hβ between the OB stars and HOLMES models reaches 1.5 dex. However, at sub-solar metallicity the change in [O III]/Hβ between the 2 Myr and 10 Gyr models is much more modest (∼0.3 dex for 12 + log(O/H) = 8.09).

In order to compare the data with these model predictions, we estimated f_ion by computing Φ_HOLMES in each spaxel of the MUSE maps from the stellar mass surface density and our fiducial value of the ionising flux of HOLMES (Sect. 3.2). The Φ_OB was then derived from the observed Σ_Hα after subtracting the contribution from HOLMES. According to this estimate, the DIG shows a median log(f_ion) = 1.2, and the 10th and 90th percentiles of its distribution across the surveyed area are log(f_ion) = 0.5 − 2.0. This range is shown as the dark grey area in Fig. 14b, and the intersection with the observed range of [O III]/Hβ, shown as a horizontal light grey band, corresponds quite closely to the location of the models. Remarkably, the range of observed [O III]/Hβ can be reproduced assuming a contribution from HOLMES consistent with the expectations from population synthesis modelling of old stellar populations. We therefore confirm that HOLMES are a ‘natural’ solution to explain the DIG line ratios in star-forming galaxies.

5.6. Mixing sequence explains the locus of the DIG in the BPT diagram

In Fig. 15 we show the position of various models along the mixing sequence in the classical BPT diagrams: log([O III]/Hβ) versus log([N II]/Hα) ([N II] BPT, left), log([O III]/Hβ) versus log([S II]/Hα) ([S II] BPT, middle), and log([O III]/Hβ) versus log([O I]/Hα) ([O I] BPT, right). We present the 2 Myr (top row) and 10 Gyr (bottom row) models, in addition to models along the mixing sequence with log(f_ion) = 1 and log(f_ion) = 0. We limit the grid to the values of log U₀ = [ − 4.0, −3.5, −3.0], because these models are the ones that more closely match the observed DIG line ratios. Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003) and solid black lines (Law et al. 2021). The latter is a recent empirical determination of the BPT demarcation lines based on the MaNGA survey.

Fig. 15. Position in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) of photoionisation models with varying fion = ΦOB/ΦHOLMES. We show models for different metallicity and log(U0) = [−4.0, −3.5, −3.0], with colours given in the legend. Top row: models where the input spectrum includes only OB stars (2 Myr SSP), while in the bottom row the input spectrum only includes old stars (10 Gyr model). Two middle rows: illustrate the effect of a decreasing contribution from OB stars to the line ratios in the BPT diagram. Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021). The grey contours represent the locus of the DIG regions from the galaxies in the PHANGS–MUSE sample (the lowest contour encloses 10% of the data, and successive contours are in steps of 20%). The comparison shows that adding a varying contribution of HOLMES allows all the line ratios observed in the DIG to be broadly covered.

In the [N II] BPT, the original line suggested by Kewley et al. (2001) was not a good representation of the data. Even though a new demarcation line was suggested by Kauffmann et al. (2003) to provide a closer representation of the sequence of star-forming galaxies in the legacy Sloan Digital Sky Survey (SDSS-I) data, researchers have continued to use the Kewley et al. (2001) demarcation line in the [N II] BPT diagram because it provides a useful ‘intermediate’ zone between the star-forming sequence and the AGN locus. The empirical Law et al. (2021) line, derived from MaNGA data, matches the Kauffmann et al. (2003) line very well.

We show the positions of all DIG regions in the three BPT diagrams in Fig. 15 as grey contours (the lowest contour encloses 10% of the data, and successive contours are shown in steps of 20%). For this comparison we exclude data from NGC 1365 to avoid its prominent AGN ionisation cone.

Figure 15 shows that with decreasing f_ion (i.e. larger contribution from HOLMES), models move ‘diagonally up’, towards the top-right corner in all diagrams. log(f_ion) = 1 (HOLMES contribute 1/10 of the ionising flux), a small number of models cross the relevant demarcation lines for the star-forming sequence into the area of the BPT diagram associated with LI(N)ERs and AGN. The shift into the LI(N)ER region is almost complete for log(f_ion) = 0 (equal contributions of HOLMES and OB stars to the ionising flux). Coincidentally, models with log(f_ion) = 0 in the [N II] BPT provide an excellent match to the old Kewley et al. (2001) line, highlighting the usefulness of this line as a probe for the DIG mixing sequence.

While the models extend into the top left area of the BPT diagrams, generally associated with Seyfert AGN, that area is occupied by models with log U₀ ≥ −3. Qualitatively, the best match between models and data is achieved with log U₀ = [ − 4.0, −3.5] and requires the full range of the mixing sequence between OB stars and HOLMES to explain the line ratios in the top-most right corner of the diagram, especially in the [N II] and [O I] BPT diagrams, while models only containing HOLMES at high metallicity over-predict [S II]/Hα. Remarkably, considering all the models along the mixing sequence, we can span the entire range of line ratios found in the data.

In Fig. 16 we show the location of the DIG regions for the entire sample (again excluding NGC 1365) in the three BPT diagrams colour-coded by the median EW(Hα) in each bin, as a proxy for the relative importance of HOLMES. As a reminder, we expect EW(Hα) ∼ 1.1 Å for pure HOLMES under our fiducial assumptions. We observe a strong correlation between EW(Hα) and the position of DIG regions in the three BPT diagrams, with EW(Hα) decreasing diagonally up in excellent agreement with the expectations of an increased contribution of HOLMES moving from the star-forming locus to the LI(N)ER regime. Lacerda et al. (2018) and Law et al. (2021) also showed similar trends between EW(Hα) and the position of regions in the BPT diagram, but this work demonstrates unambiguously that this trend is also present within the DIG in the disc of star-forming galaxies.

Fig. 16. Position of the DIG regions in our sample in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) colour-coded by the median EW(Hα). We excluded NGC 1365 because of its prominent AGN cone. Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021). The colour bar for EW(Hα) highlights the transition region between the 16th and 84th percentiles of the EW(Hα) distribution of all regions.

6. Discussion

6.1. Mean free path for ionising radiation

A key challenge for models that consider leaking radiation as the ionising source for the DIG is to explain the propagation of ionising photons within the multi-phase ISM. Early Monte Carlo 3D photoionisation models were based on the assumption of a uniform medium and needed to assume a lower midplane gas density than observed in order to fully ionise the extra-planar DIG (Wood & Mathis 2004). More recently, works coupling 3D photoionisation with more realistic ISM structures (e.g. from hydrodynamical simulations) have demonstrated that low-density chimneys naturally emerge, and allow relatively unimpeded transport of ionising photons away from the midplane (Wood et al. 2010; Barnes et al. 2015; Kado-Fong et al. 2020).

In particular, a mean free path for ionising photons of the order of ∼1 kpc implies a small filling factor of neutral gas in the galaxy disc. For 13.6 eV photons, an optical depth τ ∼ 1 is reached after a neutral column density of ∼1.6 × 10¹⁷ cm⁻². For our best-fit value for the mean free path of 1.9 kpc, we obtain a mean atomic hydrogen number density n_H of ∼3 × 10⁻⁵ cm⁻³. More energetic photons are less affected, but even for 30 eV photons, the estimated mean free path implies n_H < 10⁻³ cm⁻³ on average. Our model, therefore, requires most of the volume of the disc to be filled with either very ionised or very diffuse gas. This is in qualitative agreement with the fact that most of the DIG emission does not originate from the region where the cold gas disc lies (scale height of ∼100 pc; Ferrière 2001; Kalberla & Kerp 2009; Heyer & Dame 2015), but from a thicker ionised gas layer, where the filling factor of neutral material is expected to be significantly lower. It has also been suggested that young regions the Milky Way could be situated above or below the disc midplane (Alves et al. 2020), meaning that photons escaping such regions may be less hampered by neutral gas.

6.2. Spectral hardening of leaking radiation and alternative heating sources

Considering an inhomogeneous H II region, one can assume two different scenarios for the leakage of ionising radiation. Either the ionised gas shell around a massive cluster presents fully transparent lines of sight, allowing a fraction of the ionising radiation to leak into the ISM unprocessed (the ‘picket fence’ scenario), or radiation leaks if the ionised cloud is density-bounded, that is, if the cloud does not terminate at an ionisation front with the neutral medium but lacks the outer zone. In this case, the ionising spectrum propagating into the ISM will have been partially absorbed by neutral hydrogen and helium in the nebula (the ‘filtered radiation’ scenario). Since the absorption cross-section of hydrogen decreases with frequency, ν, roughly as ν⁻³, the escaping radiation will be harder (at least in the energy interval between the ionisation potentials of hydrogen and helium). Real H II regions may lie somewhere in between these limiting cases. The models we have discussed in this work correspond to the picket fence scenario, while in the filtered radiation scenario the spectrum ionising the DIG will be harder than the unprocessed 2 Myr SSP we have used. A harder spectrum will lead to stronger emission from low-ionisation lines but also higher [O III]/Hβ (Sokolowski & Bland-Hawthorn 1991; Hoopes & Walterbos 2003; Flores-Fajardo et al. 2011; Zhang et al. 2017; Weber et al. 2019; Pellegrini et al. 2020b).

To demonstrate the effect of hardening of the radiation field, we computed a set of density-bounded photoionisation models with the 2 Myr SSP as input, and terminated the models at different values of hydrogen column density, resulting in different escape fractions of ionising radiation. The filtered (transmitted) radiation escaping from these density-bounded regions was then used as input for a second set of CLOUDY models, aiming to reproduce DIG clouds illuminated by a filtered radiation field. In Fig. 17, we show a representative set of such models with log U₀ = −3.5, and varying metallicity and escape fraction (as shown in the legend).

Fig. 17. Photoionisation models of clouds illuminated by a hardened radiation field leaking from H II regions in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) diagrams. The input spectrum is a 2 Myr SSP, which has been transmitted from a density-bounded H II region. We consider H II regions with a range of hydrogen column densities (NH, reported in the legend). The legend also shows the escape fraction from the H II region corresponding to each column density value. All models shown have log U0 = −3.5 and varying metallicity (three representative metallicities are labelled in the left panel to demonstrate that metallicity generally increases along the x axis). Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021).

The dependence of line ratios on escape fraction is non-monotonic. While in general the largest enhancements for both low (e.g. [S II]/Hα) and high ([O III]/Hβ) ionisation lines occur for the smallest escape fractions, fluorescent production⁹ of hydrogen Balmer lines leads to a decrease in these line ratios for very small (< 1%) escape fractions (Flores-Fajardo et al. 2011; Zhang et al. 2017). Our results confirm the conclusions of Zhang et al. (2017): enhancements in line ratios in hardened radiation models are modest (< 0.2 dex) and do not allow leaky radiation models to expand sufficiently far into the LI(N)ER region of the BPT diagram. Moreover, the largest enhancements are found for escape fractions of ∼10%, which are not sufficient to power the DIG.

In order to compute the more complex process of radiation hardening in the DIG itself as a function of distance from H II regions one needs to take into account the 3D radiative transfer problem in a clumpy medium, as considered by, for example, Pellegrini et al. (2020a). In this framework, the DIG is assumed to be ionised by (filtered) radiation from H II regions, and the clumpy ISM will allow propagation of ionising photons along low-density paths. In this case, one may think of the whole galactic disc as a ‘Strömgren volume’, and the line ratios changing within this volume as they would within an H II region. In analogy to the structure of H II regions, low-ionisation line ratios ([S II]/Hα, [N II]/Hα, and [O I]/Hα) increase towards the edge of the Strömgren volume, while [O III]/Hβ decreases. Weber et al. (2019) computed 3D photoionisation models in a clumpy medium and quantitatively demonstrate this point. This effect, however, cannot explain the increase in [O III]/Hβ as a function of decreasing Σ_Hα that we observe in the DIG.

In summary, hardening of leaking radiation is not sufficient to explain the observed properties of the DIG in our galaxy sample. Predicted line ratios from 1D models do not extend sufficiently into the LI(N)ER regime to explain the more extreme observed line ratios. Moreover, 3D models predict a decrease in [O III]/Hβ with distance from the ionising sources, therefore, not matching the observed decreasing or flat trend of [O III]/Hβ versus Σ_Hα.

Additional heating processes and/or alternative sources for a hard spectral component may be able to explain the observations presented in this paper. For example, shocks models (Allen et al. 2008; Rich et al. 2011) largely overlap with the LI(N)ER section of the BPT diagram. However, we find a trend of the BPT line ratios with EW(Hα) across the entire sample – meaning that the extra heating mechanism needs to be associated with stellar mass density, supporting HOLMES.

The presence of LI(N)ER emission in barred galaxies deserves a more specific discussion. As a result of their underlying orbital structure, bars are expected to drive gas flows and shocks in the gaseous component. Hydrodynamical models predict the existence of shocks along the leading edge of the bar, encompassing the thin gas (or dust) lanes often observed in barred galaxies (Athanassoula 1992; Sormani et al. 2015). The gas within the bar is not expected to be shocked, but often shows enhanced velocity dispersions (Sun et al. 2020) due to high shear (Emsellem et al. 2015). Bars are also predicted to sweep the regions of the disc inside of the bar radius (the ‘interbar’ region) clear of cold gas, as is indeed observed by inspecting the molecular gas distribution of PHANGS galaxies (Leroy et al. 2021). Absence of cold gas and star formation within the interbar region leads to the appearance of so-called star-formation deserts (James et al. 2009; Neumann et al. 2020). Within the PHANGS-MUSE sample, 12 of our targets are classified as barred galaxies by Querejeta et al. (2021) based on Spitzer imaging (NGC 1300, NGC 1365, NGC 1433, NGC 1512, NGC 1566, NGC 1672, NGC 3351, NGC 3627, NGC 4303, NGC 4321, NGC 4535, and NGC 7496). Of these, at least six show an extensive star-formation desert in the interbar region, showing low EW(Hα) and LI(N)ER classification (NGC 1300, NGC 1433, NGC 1512, NGC 1566, NGC 3351, and NGC 3627). While somewhat less conspicuously, the interbar region is associated with lower EW(Hα) and higher fraction of LI(N)ER pixels also in the other barred galaxies.

Overall, the bar exerts an important influence on the distribution of star-formation and line emission in both the bar and interbar region, but HOLMES, rather than shocks, seem to be the likely cause of the observed line ratios in the interbar region. Moreover, our proposed solution relying on HOLMES has the advantage of not requiring fine-tuning, since the observed line ratios are reproduced by adding the predicted amount of ionising photons from HOLMES.

6.3. Consequences for SFR measurements

We have demonstrated that the DIG is powered primarily by leaking radiation from H II regions, and only in small part by a spatially diffuse, harder spectral component originating from HOLMES. This conclusion argues in favour of including the Hα contribution from the DIG in the total computation of the integrated star-formation rate (SFR) in galaxies. As a consequence, the inability to spatially map the leaking photons from H II regions does not prevent accurate SFR measurements from observations at kiloparsec resolution (as in e.g. MaNGA). However, the significant escape fraction from H II regions does caution against applying standard SFR calibrations, which assume no photon loss, to individual H II regions.

If we calculate the expected contribution from HOLMES to the Hα flux directly from the stellar mass surface density, the average contribution of HOLMES to the total Hα flux within area mapped by our MUSE data is 2% (0.9% if one uses the ionising flux predicted by FSPS for a 10 Gyr SSP). NGC 1433 shows the highest fractional contribution of HOLMES to the total Hα emission in our sample, with a value of 9%. The Hα flux here is taken to include both H II regions and DIG, but not correcting for extinction. For comparison, the fraction of Hα flux in the LI(N)ER region of the BPT diagram is on average 4%. These numbers may be interpreted remembering that our analysis of line ratios in the BPT diagram demonstrates that a contribution from HOLMES of > 10% to the total number of ionising photons can push the line ratios in the DIG into the LI(N)ER region of the BPT diagram, as shown in Fig. 15. Spaxels falling the in the LI(N)ER section of the BPT diagram therefore populate the mixing sequence between ionisation from OB stars and HOLMES, and may have contributions from both ionising sources. Taking the numbers derived above at face value, we conclude that in our sample of galaxies about half of the Hα emission in the LI(N)ER part of the BPT diagram is powered by ionising photons originating from OB stars.

The contribution of HOLMES to the Hα flux on integrated scales is small in our sample, and is therefore not likely to represent a significant correction to the SFR estimation of nearby galaxies in the blue cloud. To confirm this, we consider a galaxy on the star-formation main sequence as parametrised by Leroy et al. (2019). In this parametrisation the main sequence for local galaxies has a sub-linear slope of 0.68, leading to an estimated larger contribution of HOLMES at the high stellar masses with respect to the low-mass end. The predicted contribution from HOLMES is nonetheless very small, 0.5% for galaxies of log(M_⋆/M_⊙) = 9 going up to 2% for galaxies with log(M_⋆/M_⊙) = 11. As galaxies move below the main sequence, the contribution of HOLMES will become more important (Belfiore et al. 2018), and eventually become the dominant, and potentially only source of Hα line emission on the red sequence. We conclude that the SFR measured from Hα by kiloparsec-resolution surveys (e.g. CALIFA, MaNGA, SAMI) are likely to be robust to the small correction due to HOLMES, especially if they exclude low EW(Hα) regions from the computation, and, thanks to their low resolution, already include the contribution from radiation leaking from H II regions.

6.4. Line ratios and DIG correction

The fact that line ratios in the DIG change systematically with Σ_Hα was already evident in both Galactic and extra-galactic work (Madsen et al. 2006; Zhang et al. 2017). Here we confirmed and discussed these trends for galaxies of different masses (and therefore metallicities) and at a spatial resolution sufficient to distinguish DIG and H II regions. These trends are of particular significance in the context of determinations of the chemical abundances of H II regions. Given the lack of a comprehensive modelling framework for the DIG, several authors have attempted to empirically correct the line ratios of H II regions for DIG contamination (Sanders et al. 2017; Vale Asari et al. 2019; Espinosa-Ponce et al. 2020). Despite a large body of detailed work on individual galaxies (e.g. NGC 891; Rand 1998), IFS for a representative sample of galaxies and their DIG component is currently limited to the kiloparsec-resolution surveys (e.g. CALIFA, SAMI, MaNGA). Such studies can only access uncontaminated DIG line ratios at large distances from H II regions, corresponding to the lowest Σ_Hα levels probed in this work.

Vale Asari et al. (2019), for example, define the DIG in terms of EW(Hα) (EW(Hα) < 6 Å) using kiloparsec-resolution MaNGA data and argue that the inclusion of the DIG contamination makes little to no difference to the derived abundances. This conclusion is a direct consequence of the fact that the fraction of Hα line emission due to HOLMES is small in typical main sequence galaxies (Sect. 6.3), but it does not address the impact of the change in lines ratios in dominant component of the DIG due to leaking radiation. The impact of leaking radiation on abundance measurement on kiloparsec scales is likely to be subtle. Focusing on the DIG component due to leaking radiation, this gas is likely to have lower log U₀ than H II regions. However, metallicity calibrations implicitly assume a relation between metallicity and ionisation parameter for H II regions, which will not apply to the DIG. Moreover, tools such as IZI (Blanc et al. 2015; Mingozzi et al. 2020), which directly compare model grids with data, do not allow for the possibility of mixing models with different ionisation parameters. The corrections to integrated metallicities presented in Vale Asari et al. (2019), which only take into account the negligible fraction of the DIG powered by HOLMES, are therefore likely in need of revision.

Sanders et al. (2017), in contrast, assume the line ratios observed in the lowest Σ_Hα regions of MaNGA galaxies to be representative of the DIG component in its entirety (i.e. at all Σ_Hα levels). They further assume that the DIG contributes a substantial fraction of the total Hα emission (on average 55%), based on the SINGS Hα narrow band survey of nearby galaxies (Oey et al. 2007). While the latter assumption is in agreement with our data for DIG powered by leaking radiation, the former one associates its line emission with extreme line ratios, generally found only in the part of the DIG ionised by HOLMES. These are easy to identify observationally, since they occupy a large area on sky, but contribute a negligible fraction of the Hα flux. Using these assumptions, Sanders et al. (2017) compute corrections for the mass–metallicity relation derived from integrated spectra, which can be as large as 0.3 dex, but could potentially be overestimated.

Based on the results presented here, the PHANGS dataset will allow the validity of these DIG corrections to be tested by taking the fraction of the Hα flux and the variations in the line ratios as a function of Σ_Hα into account.

7. Summary and conclusions

This work presents a systematic study of the spatial and spectral characteristics of the DIG in a sample of 19 nearby (D < 20 Mpc) star-forming galaxies spanning a range of ∼1.5 dex in stellar mass (9.4 < log(M_⋆/M_⊙) < 11) at an average spatial resolution of ∼50 pc. Our data consist of optical IFS obtained in the context of the PHANGS–MUSE Large Programme. Unlike large IFS surveys of nearby galaxies (e.g. CALIFA, SAMI, and MaNGA), which have a resolution of approximately a kiloparsec, we were able to resolve the average distance between H II regions (∼100 − 200 pc) and therefore spatially isolate H II regions from the surrounding DIG.

We defined H II regions from Hα line emission maps by using the HIIPHOT algorithm (Thilker et al. 2000), complemented with a point-source finder. We defined the DIG as the ionised gas component found outside the H II region masks. In order to study the line ratios of this low-surface-brightness component, we performed a spatial binning based on the observed Σ_Hα, which allowed us to recover fluxes for [N II], [S II], [O III], and Hβ above the 3σ level over almost the entirety of the observed area. We summarise the main results obtained from these data below.

• The distributions of Σ_Hα and EW(Hα) in H II regions and in the DIG for the full sample are nearly non-overlapping to within 1σ, with log(Σ_Hα/erg s⁻¹ kpc⁻²) = 38.8 and EW(Hα) ∼ 16 Å representing the best boundaries between the two. Individual galaxies, however, show a wide variety in their distributions, and in several cases the tail of the DIG distribution extends to higher Σ_Hα and EW(Hα).

• Six of our 19 galaxies show extended LI(N)ER emission in their central regions, generally associated with EW(Hα) < 3 Å. These regions are consistent with being predominantly ionised by HOLMES. On integrated scales, however, the fraction of the total Hα emission due to HOLMES is small within our sample (2% on average, going up to 9% in the most extreme galaxy, NGC 1433).

• We model the propagation of leaking radiation from H II regions with a simple thin-slab model. The model takes geometric dilution of the radiation field into account and uses an effective absorption coefficient (k₀) to parametrise the mean free path (λ₀ = 1/k₀) of the ionising photons. Fitting this model to the DIG Hα emission in our sample, we obtain an average best-fit mean free path of 1.9 kpc. Despite relatively large ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$, we find that in all galaxies the model has a well-defined ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ minimum as a function of k₀ and provides a good qualitative representation of the data at all Σ_Hα levels. Our approach demonstrates that the observed Hα emission in the DIG is consistent with being powered by leaking radiation, after subtracting the small contribution to the Hα flux due to HOLMES.

• If we assume that the DIG disc has a thickness comparable to the ionised gas scale heights measured from the literature (1 kpc), then we can estimate the escape fraction from H II regions in the context of the leakage model we have developed. We find a median value of the escape fraction of 40% for our sample. This value is in good agreement with the more naive estimate of the escape fraction obtained by directly comparing the Hα flux inside the H II region masks with that outside of them, which results in a median escape fraction of 37%.

• In the leaking radiation model, the ionisation parameter is predicted to decrease following the dilution of the ionising photon flux. Assuming a roughly constant electron density, this is reflected in the decrease in Σ_Hα. Metallicity gradients complicate the analysis of line ratios as a function of Σ_Hα. In order to minimise the effect of metallicity, we consider changes in line ratios within radial bins, on average 1.5 kpc in width. In accordance with previous work, we find that low-ionisation line ratios ([S II]/Hα, [N II]/Hα, and [O I]/Hα) increase with decreasing Σ_Hα. Photoionisation models with a 2 Myr SSP as input are able to reproduce these trends and the observed values of the line ratios assuming that decreasing Σ_Hα is linked to a decrease in log U₀. Different radial bins do not show large differences in [S II]/Hα and [O I]/Hα, because the effect of metallicity on these line ratios is relatively small. [N II]/Hα, on the other hand, is found to be systematically higher in the inner regions with respect to the outer regions at fixed Σ_Hα, demonstrating the role of metallicity in setting the value of this line ratio in the DIG.

• Models of leaking radiation predict that an increase in the ionisation parameter leads to an increase in [O III]/Hβ, in clear tension with the observation of the flat or decreasing trend of [O III]/Hβ with increasing Σ_Hα. A second class of ionising sources with a harder spectrum is needed in order to explain the high observed value of [O III]/Hβ in the DIG. Moreover, the relative importance of this second component must increase with distance from H II regions in order to reproduce the observed trends with Σ_Hα. We demonstrate that the predicted ionising flux and spectral hardness of HOLMES make them the natural candidate to explain the observations. We find that the line ratios in the DIG can be explained considering a mixing sequence between HOLMES and leaking radiation.

• For a fractional contribution of HOLMES to the total ionising photon budget larger than about 10%, our models extend above the typical demarcation lines used to define the star-forming sequence in the BPT diagram. The mixing sequence spans the full range of line ratios observed in the DIG in our sample and also explains the trend between the position in the BPT diagram and the EW(Hα).

• Central regions from five of the more massive galaxies in our sample deviate from the average line ratio versus Σ_Hα relations, showing systematically higher line ratios at fixed Σ_Hα. These regions correspond to areas of low EW(Hα) and LI(N)ER-like emission in the BPT diagram. These observed properties are naturally explained by a locally dominant contribution of HOLMES to the ionisation budget.

Based on the understanding of the DIG as a mixing sequence between leaking radiation from H II regions and HOLMES, we aim to develop a detailed photoionisation modelling framework for the DIG of star-forming galaxies. Such models, which need to be tested against high-spatial-resolution data, would be a powerful tool for interpreting emission-line ratios observed on kiloparsec scales at both low and high redshift.

Acknowledgments

This work has been carried out as part of the PHANGS collaboration. This work is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme IDs 1100.B-0651, 095.C-0473, and 094.C-0623. FB acknowledges funding from the ESO-Garching fellowship during the earlier stages of development of the project. FB also thanks ASTRO3D for supporting a research stay at ICRAR, which was fundamental to the development of this project. We thank the referee for their insightful suggestions, which improved the quality of this work. ES, FS, HAP, and TGW acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). RSK and SCOG acknowledge financial support from the German Research Foundation (DFG) via the Collaborative Research Center (SFB 881, Project-ID 138713538) ‘The Milky Way System’ (subprojects A1, B1, B2, and B8). They also acknowledge funding from the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany’s Excellence Strategy (grant EXC-2181/1 – 390900948) and from the European Research Council via the ERC Synergy Grant ECOGAL (grant 855130). EC acknowledges support from ANID project Basal AFB-170002. KK gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI Kreckel). MB gratefully acknowledges support by the ANID BASAL project FB210003. ATB would like to acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 726384/Empire). MC and JMDK gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through an Emmy Noether Research Group (grant number KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1), as well as from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). The work of AKL was partially supported by the National Science Foundation (NSF) under Grants No. 1615105, and 1653300. Science-level MUSE mosaicked datacubes and high-level analysis products (e.g. emission line fluxes) are provided via the ESO archive phase 3 interface (https://archive.eso.org/scienceportal/home?data_collection=PHANGS). A full description of the first PHANGS data release is presented in Emsellem et al. (2022). The H II region catalogue used in this work will be made available in a forthcoming publication.

References

Appendix A: Atlas of Hα and H II region mask images for the full sample

In this Appendix, we show Hα maps for all PHANGS–MUSE galaxies at their native (spaxel) resolution, with the H II region masks superimposed in blue (Fig. A.1–A.2).

Fig. A.1. Ionised nebula masks (blue) superimposed on the ΣHα maps (grey scale). Galaxies are ordered by increasing stellar mass from top left to bottom right.

Fig. A.2. Same as Fig. A.1.

Appendix B: Derivation of the thin-slab model solution for the propagation of ionising radiation in the DIG

One may model the propagation of ionising radiation in a thin disc filled with a population of opaque clouds with an effective attenuation law, described by the usual radiative transfer equation,

$$\begin{array}{c}\hfill \frac{\mathit{dI}}{\mathit{dr}}=-I{k}_0,\end{array}$$(B.1)

where I(r) is the ionising photon – Lyman continuum (LyC) – surface brightness at distance r from the ionising source, and k₀ an effective absorption coefficient, which we take as a free parameter. The ionising luminosity at position r is then given by I(r) = I₀ e^−τ(r), where $I_0=I(r=0)=F^{\mathrm{H}\mathrm{I}\mathrm{I}}/4\pi$ and $F^{\mathrm{H}\mathrm{I}\mathrm{I}}$ is the H II region ionising photon flux and τ(r) = ∫k₀ dr is the optical depth.

Considering a pixel with linear size Δr, the amount of ionising radiation absorbed within it, which will be proportional to the Hα recombination line luminosity, is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }I=I(r)-I(r+\mathrm{\Delta }r)=I_0{\mathrm{e}}^{-\tau (r)}(1-{\mathrm{e}}^{-k_0\mathrm{\Delta }r}).\end{array}$$(B.2)

In a thin slab of height Δh, the flux F(r) absorbed in a pixel of size Δr is given by

$$F(r)=\mathrm{\Delta }I\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Delta }I\frac{\mathrm{\Delta }r\mathrm{\Delta }h}{r^2}=F^{\mathrm{H}\text{I}\mathrm{I}}\frac{\mathrm{\Delta }r\mathrm{\Delta }h{\mathrm{e}}^{-\tau (r)}}{4\pi r^2}(1-{\mathrm{e}}^{-k_0\mathrm{\Delta }r}).$$(B.3)

Summing over all the sources of ionising photons (H II regions) and dividing by the pixel area (Δr)², one arrives at Eq. 1, noting that conversion factors between ionising photon and Hα flux are applied to both sides of Eq. 1, and therefore cancel out.

Appendix C: Photoionisation models: Young and old star models with Cloudy v17.02 and FSPS spectra with MIST isochrones

In this paper we have presented photoionisation models based on the FSPS code (Conroy et al. 2009) and MIST isochrones (Choi et al. 2016; Dotter 2016), using CLOUDY v17.02 (Ferland et al. 2017). These models update and expand the works of Byler et al. (2017, 2018, 2019), who used CLOUDY v13.03. CLOUDY v> 17.00 include a new treatment of the S²⁺ → S⁺ dielectronic recombination coefficient. In short, dielectronic recombination (S²⁺ + e⁺ → S^+ * * → S⁺ + hν) is the main recombination pathway for doubly ionised sulphur under nebular conditions. Unfortunately the necessary atomic data to accurately calculate the dielectronic recombination rate for S²⁺, or other third row elements are not available. CLOUDY v< 17.00 therefore used the charge-normalised mean dielectronic recombination rates for the C, N and O as a guess for the dielectronic recombination for sulphur (Ali et al. 1991). From v17.0 onwards, a new empirically derived formalism for the dielectronic recombination of S²⁺ is included in CLOUDY (Badnell et al. 2015). The net effect of the change is an increase in the predicted flux of the [S II] lines, and, in particular, an increase in the [S II]/[S III] ratio of ∼20 − 50% at fixed [O II]/[O III] (see Badnell et al. 2015, their Fig. 8). We considered this update significant enough in the context of this work, and we therefore proceeded to compute models based on almost identical input physics to Byler et al. (2018) but using CLOUDY v17.02.

In particular, we generate input spectra for a 2 Myr and a 10 Gyr SSP using FSPS v3.1 and its python wrappings, PYTHON-FSPS¹⁰. The 2 Myr SSP is chosen because the slope of the ionising UV spectrum at this age is roughly consistent with that of a continuous star-formation model once it reaches equilibrium (Byler et al. 2017). The input spectra are converted into the format required for input into CLOUDY using the ad hoc script provided in CLOUDY-FSPS¹¹. The SSPs are generated using the default settings in FSPS, MIST isochrones and a Kroupa (2001) initial mass function. FSPS makes use of O and B star spectra generated with WMBasic (Pauldrach et al. 2001), Wolf–Rayet stars are taken from the spectral library of Smith et al. (2002), and post-AGB stars are modelled with non-LTE model spectra from Rauch (2003). The MILES stellar library is used for other stellar types. We note that, while the choice of the MILES stellar library sets us apart from the choices made in the Byler papers (BaSEL was used in Byler et al. 2017 and a new theoretical library still in preparation in Byler et al. 2018, 2019), we do not expect this to have a significant effect on the predicted line fluxes, which are determined by the choice of library for the hot stars.

The gas-phase abundance model is identical to Byler et al. (2018). The solar abundance model of Asplund et al. (2009) is used together with the depletion factors by Dopita et al. (2013). We scale the abundances of carbon and nitrogen with metallicity according to the fitting function of Byler et al. (2018) (their Eqs. 1 and 2; see also the discussion in their Appendix B), which are based on existing measurements of these abundance ratios in H II regions. The other heavy elements are scaled linearly with metallicity. Unlike the fiducial Byler models, we include ISM grains in our CLOUDY runs. They are found to have a small impact on line ratios, but with the largest impact on the high-metallicity high-log U₀ models (see discussion in Byler et al. 2017, their Sect. 6.3).

In the case of young stars, the ionising spectrum has a significant dependence on the metallicity of the input spectrum, since lower-metallicity stars have harder spectra. Our 2 Myr models are computed assuming the same metallicity for stars and gas. For the 10 Gyr model there is no physical reason to assume the gas clouds to have the same metallicity as the stars. Byler et al. (2019) find that, however, post-AGB ionising spectra have nearly identical shapes at all metallicities. We therefore arbitrarily fix the metallicity of the input SSP to solar for all the 10 Gyr models.

Models were run in each case for six values of gas-phase metallicity ([Z/H] = [ − 0.6, −0.4, −0.2, 0.,0.2, 0.4]) and eight values of ionisation parameter (log U = [ − 5, −4.5, −4, −3.5, −3, −2.5, −2, −1]) using PYCLOUDY. The CLOUDY computation was stopped when the temperature drops below 1000 K or the ionised fraction below 0.01. We verified that model runs with CLOUDY 13.03, no grains, and the old abundance model described in Byler et al. (2017) reproduce their model predictions based on the MIST isochrones well.

Both the 2 Myr and the 10 Gyr model grids are made publicly available¹².

All Tables

All Figures

Fig. 1. Maps of the Hα flux, binning scheme, and key line ratios for individual spaxels and after Voronoi binning in NGC 4535. a) Map of log(Hα/10−17 erg s−1 cm−2 arcsec−2), at the native 0.2″ single-spaxel resolution. The circular areas that appear white in this and other panels are locations of masked foreground stars. b) Map of log([S II]/Hα) at single-spaxel resolution. White areas have S/N < 3 for the relevant emission lines. The black contour corresponds to ΣHα = 8 × 10−17 erg s−1 cm−2 arcsec−2 or about a 16σ detection for Hα. c) Same as b) but for the log([O III]/Hβ) ratio. d) Representation of the Voronoi binning scheme adopted to recover the low-surface-brightness line emission in the DIG, as described in the text. Each bin is assigned a random colour. e) Map of log([S II]/Hα) using the binned data. f) Same as e) but for the log([O III]/Hβ) ratio.

In the text

Fig. 2. ΣHα in H II regions (in blue) and DIG (in orange) in the PHANGS–MUSE sample ordered by increasing stellar mass from left to right. DIG and H II regions are defined according to a morphological criterion (using HIIPHOT; see Sect. 2.3). The 16th, 50th, and 84th percentiles of the distributions are shown with dashed lines on each violin. The blue and orange horizontal lines (and associated colour-shaded areas) represent the 50th (and the 16th and 84th) percentiles of the ΣHα distribution for H II regions and DIG across all galaxies.

In the text

Fig. 3. EW(Hα) in H II regions (blue) and DIG (orange) in the PHANGS–MUSE sample. DIG and H II regions are defined according to a morphological criterion (using HIIPHOT; see Sect. 2.3). The 16th, 50th, and 84th percentiles of the distributions are shown with dashed lines on each violin. Dashed horizontal black lines represent demarcations employed in the literature to distinguish H II regions from DIG (3 and 14 Å) and the value of EW(Hα) from HOLMES using the PEGASE models (1.1 Å). The blue and orange horizontal lines (and associated colour-shaded areas) represent the 50th (16th and 84th) percentiles of the EW(Hα) distribution for H II regions and DIG across all galaxies.

In the text

	Fig. 4. Maps of EW(Hα) for all the galaxies in the PHANGS–MUSE sample. Galaxies are ordered by increasing stellar mass from the top left to the bottom right. The contours correspond to EW(Hα) of 1.1 (white) and 3 Å (light grey). For added clarity, regions with EW(Hα) < 1.1 Å are shown with white hatches.
In the text

	Fig. 5. Maps colour-coded by the position of each region in the [S II] BPT diagram (blue: star formation; orange: LI(N)ER; red: AGN). The classification is performed according to the demarcation lines of Kewley et al. (2001) (to classify star-forming regions) and the AGN-LI(N)ER demarcation line of Kewley et al. (2006). The black contours correspond to EW(Hα) of 1.1 and 3 Å (same as in Fig. 4).
In the text

Fig. 6. Example of our DIG modelling and its comparison with the observed ΣHα for NGC 4535. Top row: observed Hα flux map (H II regions are masked and appear white in all the maps) and the model predictions for different values of k0, generated according to Eq. (1). Middle row: a spatial comparison of the model and observations, in the form of log(Model/Data). Bottom row: median relation between log(Model/Data) and ΣHα,  observed with associated scatter (the shaded region corresponds to the range between the 16th and 84th percentiles). The ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ for each model is also listed above each panel. The model with k0 ∼ 1 kpc−1 offers a better representation of the data with respect to models with higher (k0 = 10 kpc−1) or lower values (k0 = 0 kpc−1, i.e. only spherical dilution).

In the text

	Fig. 7. Same as Fig. 6 but for NGC 4254.
In the text

Fig. 8. Summary of best-fit model parameters for the DIG leakage model applied to our entire sample. a) Variation in ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$ for the DIG spatial models described in the text as a function of the absorption coefficient, k0, for all the galaxies in our sample. Low values of k0 correspond to the case of spherical dilution, while high values correspond to exponential attenuation with scale length 1/k0. For each galaxy, a coloured star denotes the model with minimum ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$. b) Value of the additive offset, a0, defined in Eq. (2), as a function of k0. According to the simple model described in the text, a0 ∼ 0. Star symbols correspond to the k0 values with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$, which are clustered around a0 ∼ 0. c) fesc, the escape fraction of ionising radiation from H II regions as a function of k0 for each galaxy, assuming Δh = 1 kpc. Star symbols correspond to the k0 values with the lowest ${\chi }_{\text{d}.\text{o}.\text{f}.}^2$.

In the text

Fig. 9. Model predictions for changes in low-ionisation line ratios with ionisation parameter in the case of leaking radiation. Models are computed with CLOUDY v17.02 and a 2 Myr SSP spectrum generated with FSPS as input. The line ratios, from left to right log([S II]/Hα), log([N II]/Hα), and log([O I]/Hα), are plotted as a function of the ionisation parameter at the incident surface of the cloud (log(U0)). Different colours correspond to different gas-phase metallicities, as indicated in the legend. The light grey band corresponds to the 10th and the 90th percentiles of the observed line ratios in the DIG across all the galaxies in our sample. We further highlight in dark grey the area of log(U0) = [−4.5, −3.5], corresponding to the expected range of ionisation parameter values in the DIG. The models span the range of observed line ratios in the DIG for the expected values of log(U0) in the case of log([S II]/Hα) and log([O I]/Hα), while they fall somewhat short of reproducing the highest line ratios for log([N II]/Hα).

In the text

Fig. 10. Log([S II]/Hα) in the DIG as a function of log(${\mathrm{\Sigma }}_{\mathrm{H}\alpha }/{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{kpc}}^{-2}$) for different radial bins in individual galaxies in the PHANGS–MUSE sample. Galaxies are ordered by stellar mass from lowest (top left) to highest (bottom right). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, and NGC 1433) have higher values of this (and other) line ratios in their central radial bin (in red). Aside from this feature, the galaxies in the sample follow a remarkably tight relation between log([S II]/Hα) and ΣHα.

In the text

Fig. 11. Same as Fig. 10 but for log([N II]/Hα). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) have exceptionally high values of this (and other) line ratios in their central radial bin (in red). Aside from this feature, the galaxies in the sample follow an average relation between log([N II]/Hα) and ΣHα, with a clear secondary dependence on metallicity. Inner regions have higher metallicity and show elevated [N II]/Hα at fixed ΣHα.

In the text

	Fig. 12. Same as Fig. 10 but for log([O I]/Hα). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) show exceptionally high line ratio values in their central regions (in red).
In the text

Fig. 13. Log([O III]/Hβ) in the DIG as a function of log(ΣHα/erg s−1 kpc−2) for different radial bins for individual galaxies in the PHANGS–MUSE sample. Galaxies are ordered by stellar mass from lowest (top left) to highest (bottom right). Several galaxies (NGC 3351, NGC 1300, NGC 1512, NGC 3627, NGC 1433, and NGC 1365) show exceptionally high line ratio values in their central regions (in red). Galaxies present a systematic change in the slope of the log([O III]/Hβ) versus ΣHα relation, with low-mass galaxies showing comparable line ratios at low and high ΣHα and high-mass galaxies showing a clear decrease in log([O III]/Hβ) with ΣHα. Both of these trends are at odds with predictions from photoionisation models.

In the text

Fig. 14. Model predictions for the [O III]/Hβ line ratio considering both young stars and HOLMES. a) Model predictions for photoionised clouds, computed with CLOUDY and a 2 Myr (solid lines) and a 10 Gyr (dashed lines) SSP spectrum. The line ratio log([O III]/Hβ) is plotted as a function of the ionisation parameter at the incident face of the cloud (log(U0)). The vertical extent of the grey area highlights the range of line ratios observed in the DIG in the PHANGS–MUSE data. Different colours correspond to different gas-phase metallicities, with the same colour-coding in both panels, as indicated in the legend. All models predict log([O III]/Hβ) to increase with ΣHα, in contrast to what is observed. b) Model predictions for the log([O III]/Hβ) line ratio as a function of the fraction of the ionising radiation originating from old stars, fion = ΦOB/ΦHOLMES. We show models for a fixed value of log(U0) = − 3.5 and different metallicities. The horizontal light grey band corresponds to the range of values of log([O III]/Hβ) observed in the DIG of our sample, while the dark grey box corresponds to the range of fion observed in the PHANGS data. fion is estimated by computing the ionising photons expected to be produced by HOLMES given the local stellar mass surface density.

In the text

Fig. 15. Position in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) of photoionisation models with varying fion = ΦOB/ΦHOLMES. We show models for different metallicity and log(U0) = [−4.0, −3.5, −3.0], with colours given in the legend. Top row: models where the input spectrum includes only OB stars (2 Myr SSP), while in the bottom row the input spectrum only includes old stars (10 Gyr model). Two middle rows: illustrate the effect of a decreasing contribution from OB stars to the line ratios in the BPT diagram. Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021). The grey contours represent the locus of the DIG regions from the galaxies in the PHANGS–MUSE sample (the lowest contour encloses 10% of the data, and successive contours are in steps of 20%). The comparison shows that adding a varying contribution of HOLMES allows all the line ratios observed in the DIG to be broadly covered.

In the text

Fig. 16. Position of the DIG regions in our sample in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) colour-coded by the median EW(Hα). We excluded NGC 1365 because of its prominent AGN cone. Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021). The colour bar for EW(Hα) highlights the transition region between the 16th and 84th percentiles of the EW(Hα) distribution of all regions.

In the text

Fig. 17. Photoionisation models of clouds illuminated by a hardened radiation field leaking from H II regions in the [N II] BPT (left), [S II] BPT (middle), and [O I] BPT (right) diagrams. The input spectrum is a 2 Myr SSP, which has been transmitted from a density-bounded H II region. We consider H II regions with a range of hydrogen column densities (NH, reported in the legend). The legend also shows the escape fraction from the H II region corresponding to each column density value. All models shown have log U0 = −3.5 and varying metallicity (three representative metallicities are labelled in the left panel to demonstrate that metallicity generally increases along the x axis). Demarcation lines from the literature are shown as dashed (Kewley et al. 2001, 2006), dotted (Kauffmann et al. 2003), and solid black lines (Law et al. 2021).

In the text

	Fig. A.1. Ionised nebula masks (blue) superimposed on the ΣHα maps (grey scale). Galaxies are ordered by increasing stellar mass from top left to bottom right.
In the text

	Fig. A.2. Same as Fig. A.1.
In the text