© G. Micheva et al. 2022

1 Introduction

H II regions are tracers of recent star formation in the interstellar medium (ISM) and typically display well-detected nebular emission lines such as the Balmer lines, [O iii]λλ4959, 5007, [S II]λλ6716, 6731, and [N II]λλ6548, 6583 (e.g., Shields 1990). The diffuse ionized gas (DIG) is another component of the ISM, contributing a significant part (~50%) of the total Hα luminosity in late-type spirals (e.g., Haffner et al. 2009). This component is often also referred to as the warm ionized medium (WIM) in the literature (e.g., Haffner et al. 2009; Weilbacher et al. 2018). The properties of the DIG have been shown to be significantly different not only from those of classical H II regions (Madsen et al. 2006), but also from one location to the next within the same galaxy, with a range of measured electron densities and temperatures, and a great variation in the [O iii]λλ4959, 5007/Hα, [S II]λλ6716, 6731/Hα, and [N II]λλ6548, 6583/Ha emission line ratios (e.g., Mathis 2000; Haffner et al. 2009). This led to the conclusion that there is no single emission line value or threshold in the equivalent width of the lines that allows for a separation of DIG-dominated spaxels from those with non-DIG emission (e.g., Tomičić et al. 2021). The DIG has been shown to be ionized by Lyman continuum radiation escaping from the nearby H II regions at least partially (e.g. Reynolds 1984; Ferguson et al. 1996; Oey et al. 2007; Chevance et al. 2020; den Brok et al. 2020; Della Bruna et al. 2020; Belfiore et al. 2022). This is supported by the fact that the DIG is usually observed in the vicinity of H II regions (e.g., Ferguson et al. 1996; Zurita et al. 2000, 2002; Madsen et al. 2006; Howard et al. 2018). In order to study DIG-dominated regions, therefore, a selection criterion based on physical proximity to nearby H II regions is often used (e.g., Weilbacher et al. 2018; Della Bruna et al. 2020; den Brok et al. 2020). There are also alternative scenarios for the DIG ionization in which the DIG emission is attributed to turbulence (Binette et al. 2009), shocks (Collins & Rand 2001), and contributions by hot evolved low-mass stars (HOLMES) as defined by Flores-Fajardo et al. (2011), and supported by Belfiore et al. (2022). None of these, however, offer a robust and systematic solution to determining the location of the DIG-dominated regions. Understanding the characteristics of H II and DIG regions is fundamentally important for understanding the star formation process in galaxies. Most of the H II regions studied in the literature are rarely fainter than 10³⁶ erg s⁻¹ in Hα luminosity (e.g., Youngblood & Hunter 1999; Bradley et al. 2006; Hakobyan et al. 2007; Weilbacher et al. 2018; Della Bruna et al. 2020), with the DIG being similarly bright. To study fainter samples of both H II and DIG regions, we must therefore look to nearby galaxies.

NGC 300 is a flocculent spiral galaxy in the Sculptor Group and one of the closest galaxies to the Local Group (Sérsic 1966; Gieren et al. 2005). Due to its proximity and therefore the higher spatial resolution it offers, it is the subject of numerous studies on its own merit and as a case study (e.g., Sérsic 1966; Deharveng et al. 1988; Galvin et al. 2012; Stasińska et al. 2013; Faesi et al. 2014, 2016; Gazak et al. 2015; Niederhofer et al. 2016; Hillis et al. 2016; Rodríguez et al. 2016; Kang et al. 2016; Riener et al. 2018; Roth et al. 2018; Mondal et al. 2019; Gross et al. 2019; McLeod et al. 2019). The MUSE observations of NGC 300 are presented at length in the pilot paper of Roth et al. (2018), which showed the use of MUSE for crowded field 3D spectroscopy. The observed fields in NGC 300 are purposefully chosen to cover much of the central regions of this galaxy while avoiding the brightest H II regions. Our work is a continuation of this project, but concentrates exclusively on the study of the H II regions and the DIG. In this paper we study extremely faint H II regions and even fainter DIG. We aim to construct a sample of the H II regions and localize the DIG. We make use of the high spatial resolution that MUSE affords in NGC 300 to characterize both samples in terms of their physical properties, study their spatial distribution, and highlight any detected differences between these ISM components.

Similar to Roth et al. (2018), the first paper in this series, the adopted distance to NGC 300 is 1.88 Mpc (Gieren et al. 2005; Bresolin et al. 2005). The sample of H II and DIG regions we use in our analysis is available at the CDS.

The paper is organized as follows: in Sect. 2 we briefly present the data and give details of the line extraction. The classification of emission line regions is explained in Sect. 3. In Sect. 4 and its subsections, we derive physical properties for the H II and DIG regions such as metallicity, extinction, and electron density, as well as kinematic properties such as velocity dispersion and shear. To increase the signal-to-noise ratio of the data, we stack the spectra per field in Sect. 5. We discuss our findings in Sect. 6 and summarize our conclusions in Sect. 7.

2 Data

Several fields of NGC 300 were observed as part of the Guaranteed Time Observations with the MUSE instrument (Bacon et al. 2010) in 2014–2016¹. All observations used the wide-field mode, which covers 1′ × 1′, and individual exposures of 900 s were made, with 90° rotations in between exposures and with an offset sky exposure taken after two science exposures. The initial observations including all exposures taken in 2014 and 2015 are described by Roth et al. (2018). They included the fields A, B, C, D, E, J, and I. All of these observations were carried out without adaptive optics (AO) support (extended Wide Field Mode with natural seeing, WFM-NOAO-E), with a contiguous wavelength coverage of 4600–9350 Å. While all of the fields were targeted to be 1.5 h deep (6 × 900 s), fields B, E, and J had only partial coverage. After combination, fields B and C had an effective image quality of 1″.0 or lower. In August and September 2016, observations of fields D (two new exposures) and E (four new exposures) were completed using the same mode. Two new fields (P and Q) at larger radii from the galaxy center were added, both receiving the full coverage 6 × 900 s in good atmospheric conditions. In October and November 2018, the innermost three fields (A, B, and C) were reobserved, this time with AO support (Ströbele et al. 2012) in good conditions with an average external seeing of 0″.75. They were again observed in extended mode (Wide Field Mode with Adaptive Optics - Extended, WFM-AO-E, 4600–9350 Å), except for two exposures of field A, which by mistake were taken in nominal mode (Wide Field Mode with Adaptive Optics - Nominal, WFM-AO-N, 4700-9350 Å). In AO-supported observations, a wide gap around the NaD line blocks out the laser light (about 5750–6010 Å). Again, 120 s offset sky observations followed each on-target exposure. Illumination flat-fields were taken before or after the observing blocks, and a standard star was observed in the same instrument mode in the same night.

2.1 Reductions

The data were reduced with the standard MUSE pipeline (Weilbacher et al. 2020), but we used the version integrated into the MUSE-Wise (Vriend 2015) environment. MUSE-Wise automatically associates the appropriate calibrations for each science exposure. Due to the long time span of the observations, the data were processed with different software versions. The basic CCD-level processing of the data of 2014 and 2015 is the same as used by Roth et al. (2018) and was performed with v1.0 and v1.2; the 2016 data use v1.6. A stack of all exposures was created with v2.4. The AO data of 2018 then used a development version (v2.5.4) that soon after was publicly released as v2.6². As usual for most MUSE data, all science exposures were corrected for the CCD bias, but did not correct for the negligible dark current. Flat-fields and arc exposures taken in the morning after the science data were used to correct for pixel-to-pixel response, to locate the illuminated areas on the CCDs, and to assign wavelengths to each pixel. Spatial positions in the field of view were taken from the geometric calibration specific to the corresponding observing run. Further flat-fielding was applied to minimize spatial slice-to-slice variations using the corresponding illumination-flat exposure and to correct global gradients across the field with sky-flats taken in the same mode within a few days of the science exposure.

The data were then corrected for atmospheric refraction and were flux calibrated. Because the observations were obtained during clear or photometric nights, a single standard star in the same night and instrument mode was used to determine the response curve and also to apply a correction of the telluric absorption features. The sky subtraction first employed the creation of a sky continuum spectrum from the offset sky field taken closest in time. To do this, the sky spectrum was averaged over 80% of the sky field, and the emission lines were fitted and subtracted. The fitted sky emission line fluxes were then taken as first-guess inputs to the fit of the sky lines over a 5–10% fraction of the science field, from which the sky continuum was directly subtracted. The telluric Hα emission is only imperfectly removed by this procedure; all other lines typically leave 1–2% residuals in the science data. The data were further corrected for distortions using an astrometric calibration of the corresponding observing run, and were shifted to barycentric velocity. One AO exposure of field B (with header DATE-OBS=2018-11-10T00:16:43) was contaminated by a satellite. After finding the trail in visual inspection of the cube, the pixels corresponding to ±5 pixels from the peak were masked and the data were processed again. Relative offsets were corrected before combining all exposures of each field. In case of the AO fields, all offsets were computed relative to the Gaia DR2 (Lindegren et al. 2018) sources. For the non-AO fields, offsets were computed relative to the first exposure, and the world coordinates of the final cubes were manually adjusted to have absolute astrometry in agreement with the Gaia DR2 positions. The absolute positional accuracy is on the order of one spatial MUSE pixel or better, that is, <0″.2.

All cubes were constructed with the default 0″.2 × 0″.2 spatial sampling and have approximately 60″ × 60″ spatial coverage. The wavelength coverage is about 4600–9350 Å and sampled at 1.25 Å per pixel. The exact wavelength value of the first pixel depends on the velocity correction applied to the data.

2.2 Definition of regions

Before the emission line flux could be extracted, we prepared the data in the following way. A first-order approximation of the Hα spatial distribution was obtained for all fields for the AO and non-AO data using the Hα narrowband filter range defined in Roth et al. (2018). To reduce noise, the Hα maps were smoothed with a 2D Gaussian kernel with a standard deviation equal to 1.25, using astropy.convolve. Unfortunately, the signal-to-noise ratio (S/N) in individual spaxels is low, and hence reliable spectral information on the scales of individual pixels often cannot be obtained. Therefore, we binned the preliminary Hα maps in dendrograms, using the astrodendro³ package (Robitaille et al. 2019). Through trial and error, we determined that the best parameters for dendrogram extraction, resulting in a reduced number of spurious detections and yet capturing all main structures in the fields, are min_delta= 0.005 and min_npix = 10. For the purpose of illustration, in Fig. 1a we show an example region in field P and some typical dendrogram leaf masks. We manually examined the dendrogram tree, and found that more often than not, only the compact peaks of otherwise large filaments and bubbles were captured by the automatic algorithm. We merged down and trimmed leaves and branches until the resulting shapes traced the entirety of the visible filaments and bubbles. This procedure is subjective in the sense that there are cases in which the morphological appearance suggested to merge individual leaves that were not reported as connected by the automatic procedure. The spectra in each dendrogram leaf were summed.

Next, as an adaptation of classical aperture photometry techniques (e.g., Stetson 1987), the diffuse background contribution was estimated locally around each dendrogram leaf in a ring with a width of 7 pixels, with a gap between the dendrogram and the background region of 5 pixels. The ring shapes were determined via binary dilation of the leaf masks. These ring regions are illustrated in Fig. 1b. Nearby dendrogram leaves overlapping the rings in crowded areas were excluded from the ring masks. For each background ring, all spectra were ordered in increasing continuum levels, as measured in the wavelength range of 6320 Å ≤ λ ≤ 6480 Å, where no strong lines could be seen. This ordering was used to identify the bottom 25^th percentile of the background spectra. The median of these was taken as the final representation of the local background and was subtracted from the corresponding dendrogram leaf. The underlying assumption here is that the diffuse background in each ring is representative of the background inside its corresponding dendrogram leaf, and that the convolution of the true image of a nebula with the point spread function (PSF) has become negligible at the location of the ring.

The ring masks serve the additional purpose of representing the DIG regions. In a separate process, the spectra in each ring mask were summed to represent the local DIG. This method resulted in large gaps with sizes on the order of hundreds of spaxels between leaves and DIG rings, which remain unassigned to any of the binned regions and therefore outside of our analysis. The diffuse gas signal there is very low, and we binned them into bins that were as large as possible, resulting in an additional 113 DIG regions with average equivalent radii of 27.8 pc. These regions were manually added to the list of DIG regions. An example is shown in Fig. 1c, where these regions are shaded pink and fill the available space between ring regions.

Fig. 1 Region definitions for an example region in field P. Overplotted are (a) dendrogram leaves, (b) DIG ring regions, which also serve as sky background regions, and (c) dendrogram leaves and DIG rings both, as well as the manually added DIG regions, shaded in pink.

2.3 Line extraction

We measured emission line fluxes of all sources using the pPXF Python program (v6.7.17 Cappellari & Emsellem 2004; Cappellari 2017). The advantage of this program is that it can model any residual stellar contribution to the extracted spectra, and fit emission lines in NGC 300 as well as the residual telluric Hα at the same time. After testing various spectral libraries of individual stars and simple stellar populations (SSPs), we settled on the GALEV (Kotulla et al. 2009) SSPs computed with the Munari et al. (2005) spectral library as created by Weilbacher et al. (2018). As a second component, we input 32 typical lines of hydrogen and helium as well as significant forbidden lines as listed in Table 1 in the wavelength range of the MUSE data. We masked the wavelength ranges around strong sky emission lines, and for AO, we also masked the region in the NaD gap. We excluded the [O i] lines from the fit because they are still in the masked area at the redshift of NGC 300. All object emission lines were fit as one component with consistent velocity, and hence all contribute to the fit of the ionized gas kinematics simultaneously. As a third component to the fit, we included an Hα line at near-zero redshift, which effectively models the strong residual of the telluric Hα emission. Because our main aim was to extract the emission line fluxes, we included a multiplicative polynomial of degree 15 in the fit and iteratively fit each spectrum. The first iteration determines the best fit and allows us to compute the residual noise. We used the residual noise in the wavelength range 6050–6250 Å to scale the variance of each spectrum. This is useful because binning the MUSE cubes causes the actual variance to be underestimated because of correlation across spaxels (see Weilbacher et al. 2020 and, e.g., Herenz et al. 2020). Improving the variance allowed us to obtain realistic error estimates for the kinematics in a second iteration of the fit. We also used the scaled residual noise to build a 6-pixel sliding standard deviation spectrum, which served as input to 100 Monte Carlo iterations of the emission line fit. We therefore expect to obtain realistic flux error estimates for all fitted emission lines. As output of pPXF, we then obtained kinematics (V, σ) for both stars and ionized gas, error estimates for them, and emission line fluxes and corresponding errors for every spectrum.

2.4 Sample sizes

The total number of dendrogram leaves is 661. We discarded 68 leaves whose spectra looked like those from emission line stars (compact regions, broad Hα, and none or extremely weak other emission lines), and an additional 45 leaves that were completely noise dominated after the background subtraction. In Sect. 3, we describe the further cleaning of these leaves, which decreased the total number to 499 dendrogram leaves.

The total number of DIG regions are 774. Of these, 661 were created by binary dilation of the leave masks, and 113 were the manually added, large, and extremely faint DIG regions described in Sect. 2.2. The further cleaning of the DIG regions, described in the next section, resulted in a total of 696 DIG regions.

3 Classification of emission line objects via a BPT diagram

To distinguish between H II regions, supernova remnants (SNR), and planetary nebulae (PNe), the traditionally used diagnostics are based on the relative strengths of the Hα, [N II]λλ6583, and the sum of [S II]λλ6716+[S II]λλ6731 (hereafter [S II]λλ6716 + 31) emission lines. These diagnostics are defined in greater detail in Fesen et al. (1985), Frew & Parker (2010), and Sabin et al. (2013), for example, and were used for MUSE datacubes to identify PNe in Roth et al. (2021).

In Fig. 2, we show the Baldwin et al. (1981) diagram (BPT) with respect to [S II]λ6716 + 31 and [N II]λ6583 and the theoretical limit of star-forming regions of Kewley et al. (2006) and the PNe/SNR separation line (Sabin et al. 2013). These separation lines are typically used in the literature to distinguish between H II, PNe, and SNR (e.g., Frew & Parker 2010; Sabin et al. 2013). The Kewley et al. (2006) line is the revised extreme starburst line of Kauffmann et al. (2003), obtained via a semiempirical fit to the outer bound of galaxy spectra. To alleviate the concern that H II regions may not share the same properties as star-forming galaxies, and hence that the Kewley line is not an appropriate separator, we also plot the line from Dopita et al. (2000), which is based on instantaneous burst modeling of H II regions. The figure shows that the Kewley and Dopita lines are very similar in the [N II]λ6583/Hα diagram. We therefore performed the formal separation of the objects into H II, PNe, and SNR in the [N II]λ6583-BPT diagram, but we kept in mind that their location with respect to [S II]λ6716 + 31/Hα may vary. This is shown in Fig. 2, where objects classified based on [N II]λ6583/Hα (top panel) cross over the separation lines in the [S II]λ6716 + 31/Hα diagram (bottom panel), highlighting the continuous nature of the sample properties. While our goal is to obtain a sample of H II regions, we are aware that the emission line properties comprise a continuous distribution, and hence we cannot claim that the blue data points in Fig. 2 constitute a pure H II region sample.

Without cleaning the sample, we obtain 429 H II-regions, 96 SNR, and 23 PNe over all nine fields with the [N II]λ6583/Hα classification. With the [S II]λ6716 + 31/Hα classification, the numbers instead are 334, 203, and 17, respectively. The numbers do not add up to the same total because some dendrogram leaves have a low S/N in [Ν II]λ6583 and some in [S II]λ6716 + 31, and so the number of points in the top and bottom panels of Fig. 2 is different. Next, we adjusted the labels for five extended bubble and filament-like structures in fields C, D, and P, which otherwise fall in the PNe region in the [N II]λ6583-based BPT diagram. Because PNe are compact objects, we assigned either an H II or SNR label to these structures, depending on their label in the [S II]λ6716 + 31-based BPT diagram. The assigned SNR labels are consistent with the log₁₀([S II]λ6716 + 31/Hα) > −0.5 criterion used to separate PNe from SNR in Roth et al. (2021). We compared our resulting labels to the list of identified PNe in Roth et al. (2018) and find that while most of the 22 objects in common have matching PNe labels, 6 of the objects are not identified as PNe by the [N II]λ6583-based BPT diagram. We visually examined these dendrogram leaves to confirm that they are small, compact, and somewhat symmetrical regions, and manually relabeled these regions as PNe. The number of PNe detections is consistent with the more targeted search for such objects in Paper IV of this series, which is concerned with the planetary nebula luminosity function in NGC 300 (Soemitro et al. 2022). We filtered the remaining sample to discard all spectra with an S/N < 5 in Hα and the [S II]λλ6716, 6731 (hereafter, [S II]λ6716, 31) doublet. Because the MUSE pointings purposefully avoid bright H II regions (Roth et al. 2018), the population of H II regions that we observe is likely older, fainter, and not very strongly ionized. We therefore did not add a selection criterion based on [O iii]λ5007 emission. Finally, we visually examined the spectrum of each region and the pPXF fit to the emission lines and discarded any remaining clearly failed fits. The final numbers are 390 H II-regions, 82 SNR, and 27 PNe over all nine fields, based on the [N II]λ6583-BPT diagram.

The DIG regions were cleaned separately by removing spectra that were too noisy, resulting in either a failed pPXF fit, or in large uncertainties consistent with non-detections in the brightest forbidden lines. This manual inspection removed 78 regions, leaving a total of 696 DIG regions.

The spatial distribution of objects is shown in Fig. 3. The H II regions can be found in every pointing, with a noticeable dearth of regions in the inter-arm field J. The DIG is also ubiquitous, but again, field J is the exception. The SNR are mostly found in the central field A, with a possible second number density peak in the inter-arm field J, trailing the spiral arm that goes through fields C, D, and I. We note that field E only partially covers the spiral arm, and partially covers another inter-arm region, which likely accounts for the lack of detected regions in the upper part of E, which overlaps with the inter-arm region.

In Appendix A, we attempt to classify the emission line objects into H II, SNR, and PNe via machine learning. We conclude, however, that an unsupervised learning algorithm such as UMAP and a labeling algorithm such as HDBscan are unable to separate the objects into the desired groups when the feature sets listed in Table A.1 are used.

Fig. 2 BPT diagram in [N II]λ6584/Hα (top) and [S II]λ6716, 31/Hα (bottom). The objects are classified in the [N II]λ6584/Hα diagram via the Kewley et al. (2006; solid line) star formation limit and the PNe-SNR (dashed line) separation lines. The Dopita et al. (2000; dotted line) is also shown for reference. The color-coding in the [S II]λ6716, 31/Hα diagram is based on the [N II]λ6584/Hα classification.

4 Physical properties

4.1 Brightness and size

The Hα surface brightness Σ_Hα of the H II dendrograms ranges between 3 × 10⁻²¹ and 1.6 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻², with a mean and median of 8.0 × 10⁻¹⁹ and 2.5 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², respectively, as illustrated by the histograms in Fig. 4a. We obtained their equivalent radius via $r_{\text{equiv}}=\sqrt{N_{\text{pix}}/\pi }$, shown in Fig. 4b. The radius ranges between r_equiv ∈ [3.3pc, 91pc], with a mean and median of 14.3, and 11.7pc, respectively. The diffuse regions are about eight times fainter on average, with an Hα surface brightness ∈ [1.4 × 10⁻²², 2.2 × 10⁻¹⁸] erg s⁻¹ cm⁻² arcsec⁻², with a mean and median of 1.0 × 10⁻¹⁹ and 5.2 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², respectively. They are also larger on average by design because they are constructed via dilated masks, r_equiv ∈ [4pc, 197pc], with a mean and median 20 and 18.8 pc, respectively.

In terms of Hα luminosity, the H II regions range between log₁₀ L_Hα [erg s⁻¹] ∈ (33.2, 37.3), and have a mean and median of 34.8 and 34.7, respectively. This makes our H II sample one of the faintest in the literature to date. The DIG, represented by the ring masks in Sect. 2.3, is similarly faint, with log₁₀ L_Hα [erg/s] ∈ (33.6, 36.5), and a mean and median of 35.0 and 34.9, respectively. This is visualized in Fig. 4c. In Fig. 5, we display the luminosity-size relation. We note that the median luminosity of 34.9 corresponds to regions that are ≳ 100 spaxels in size (≳ 10 pc equivalent radius), and hence we conclude that the extreme faintness of the sample cannot be the result of random peaks in the noise and must instead be due to the physically small sizes.

Although we do not analyze SNR and PNe regions in this paper, we mention the range of their Hα luminosities for completeness. For PNe, log₁₀ L_Hα.[erg s⁻¹] ∈ (33.1, 36.1), with a mean equal to the median of 34.6. For SNR, log₁₀ L_Hα.[erg s⁻¹] ∈ (33.1, 36.7), with amean and median 34.5 and 34.3, respectively.

Fig. 3 2D map of the BPT classification.

Fig. 4 Normalized histograms of (a) the Hα surface brightness ΣHα (10−20 erg s−1 cm−2 arcsec−2) of the H II and DIG dendrogram regions, (b) the equivalent radius (pc), and (c) the Hα luminosity (erg s−1).

Fig. 5 Hα luminosity (erg s−1) as a function of region size in spaxels (bottom x-axis) or equivalent radius in pc (top x-axis).

4.2 Electron temperature and density

Unfortunately, the auroral [S iii] λ6312 and [N II] λ5755 lines in our data are so weak as to be completely in the noise of individual regions, and hence we cannot directly map an [S iii]-or [N II]-based estimate of the electron temperature T_e. This is likely due to the near-solar metallicity of NGC 300 (e.g., Stasińska et al. 2013). The [O iii] λ4363 line is outside of the wavelength coverage, and in addition the [O II] doublet at λλ7320, 7330 is also extremely noisy. The ratio [S II]λ6731/16 is an electron density (N_e) tracer and is well detected in the data.

For calculations of N_e, we used the PyNeb (Luridiana et al. 2015; Morisset et al. 2020) software package, with sulfur atomic data from Podobedova et al. (2009) and collision strengths from Tayal & Zatsarinny (2010). In Fig. 6, we show the [S II]λ6731/16 distribution of the H II and DIG regions, as well as the theoretical [S II]λ6731/16 emissivity grid as a function of T_e and N_e. The theoretical limits of the [S II]λ 6731/16 ratio are 0.6806 ≲ [S II]λ6731/16 ≲ 2.351. Because the observed distribution peaks at the theoretical minimum, implying very low densities for most of the regions (Fig. 6, upper panel), we would expect some regions to fall left of this minimum due to statistical scatter. In Appendix B we motivate our choice to include [S II]λ6731/16 down to 0.548 for H II regions and 0.638 for the DIG.

The theoretical model in the lower panel of Fig. 6 is degenerate in T_e for [S II]λ6731/16 between approximately 0.706 to 2.27, but the [S II]λ 6731/16 contours are fairly straight, and hence we can give an acceptably narrow range of possible N_e values from the [S II]λ 6731/16 ratio alone. The peak of the [S II]λ 6731/16 histogram falls in the small range of 0.6806 ≲ [S II]λ6731/16 ≲ 0.70, where the contours curve with T_e. This range is highly sensitive to the accuracy of the atomic data, and hence here we can only say that N_e ≲ 30 cm⁻³, but are unable to make a statement about T_e.

Figure 7 shows the estimated range in N_e as a function of galactocentric distance for H II and DIG regions. Toribio San Cipriano et al. (2016) reported electron temperatures in the range of 7–10 kK for a handful of H II regions in NGC 300, but they are much larger and brighter than our sample. We therefore did not constrain the temperature range based on their results and instead assumed that all of our H II regions lie within a temperature range of T_e= 4–20 kK, which is hardly constraining. This resulted in a range of possible N_e, represented by a blue-shaded area between the low and high N_e limits for each data point. The same was done for the DIG regions, but shaded orange in the figure. We also computed the average N_e between the extremes for each object and from these values the medians of the full samples, with N_e ≈ 100 cm⁻³ for H II (solid blue line) and N_e ≈23 cm⁻³ for DIG (solid orange line).

4.3 Extinction

The wavelength range covers both Hα and Hβ, allowing us to estimate the dust attenuation via the Balmer decrement. Roussel et al. (2005) investigated extinction variations in NGC 300 and reported that compact and diffuse H II regions may be associated with very different extinction laws. However, 80% of the H II regions in their sample are larger than 90 pc in diameter, and 100% are larger than 55 pc, while 93% of our sample are smaller than 30 pc in projected diameter, and 98% are smaller than 55 pc. It is therefore not immediately obvious whether we can apply their findings in the choice of extinction law. Furthermore, Tomičić et al. (2017) found that a weakly or nonattenuated DIG affects the global attenuation. Global attenuation is not included in most dust-gas geometry models. We therefore chose the Cardelli et al. (1989) extinction law, which is applicable to the dense and to the diffuse ISM, with the default value of R_V = 3.1.

Figure 8 shows E(B − V) as a function of galactocentric distance D for individual H II and DIG regions. By propagating the uncertainties of the objects, we also constructed the uncertainty-weighted average E(B − V) in bins containing 20 objects each for both DIG and H II, which show a similar fluctuation in the data, but with a narrower range. While individual regions show highly variable E(B − V), the general trend is for a flat radial profile for both H II and DIG, with the former being ~0.1 mag redder than the latter.

Because the 1D plot with galactocentric distance erases any possible azimuthal trends in the fields, we also show in the bottom panel of Fig. 8 the 2D map across all pointings and for all objects with derivable E(B − V), including PNe and SNR. We see very little measurable variation across all pointings, but individual H II or SNR regions occasionally reach higher extinction values. The Balmer decrements of most of the DIG are consistent with the theoretical value, and hence an E(B − V) = 0. These values are omitted. Apparently, only DIG regions in close proximity to H II regions have E(B − V) > 0. The predominantly low DIG extinction values support the typical assumption in the literature that the DIG suffers no extinction (e.g., Belfiore et al. 2022).

On a final note, the range of extinction values for our H II regions is consistent with the value of Roussel et al. (2005), who obtained an A(Hα) = 0.6 mag average extinction with RMS = 0.64 mag for 23 H II regions in NGC 300, which compares well to our A(Hα) = 0.53 mag average with RMS = 0.67.

Fig. 6 Sample electron densities. Top: histogram of the [S II]λ6731/16 ratio for H II (solid blue lines) and DIG (dashed orange lines). The area hatched with backslashes represents non-physical values, and the forward slash covers the range of Ne that is fairly invariant with Te. Bottom: [S II]λ6731/16 emissivity grid for 3kK ≤ Te ≤ 20 kK and 1 ≤ log10 Ne(cm−3) ≤ 7. Contours of [S II]λ6731/16 are overplotted in white, and the corresponding value is shown inline.

Fig. 7 Electron density Ne as a function of galactocentric distance D for H II (blue) and DIG regions (orange). At each D, the width of the Ne range is shown, obtained from Ne(Te= 20 kK) and Ne(Te= 4 kK). The medians for H II and DIG are also shown with solid lines.

Fig. 8 Sample nebular extinction. Top: E(B − V) as a function of galactocentric distance D for individual H II (black dots) and DIG regions (blue diamonds). Uncertainty-weighted averages in bins of 20 objects are also shown in red circles for H II and in orange diamonds for DIG. The propagated error bars are often smaller than the markers. The x-axis position of each uncertainty-weighted average is the average galactocentric distance for the dendrograms inside each bin. Bottom: 2D map of E(B − V).

Fig. 9 Sample metallicity. Top: 12 + log10 O/H as a function of galactocentric distance D, using the Marino et al. (2013) O3N2 calibration. Same legend as in Fig. 8, except for the data by Bresolin et al. (2009; open gray circles), which were obtained via the direct-Te method for NGC 300. Bottom: same as top, but using the Pilyugin & Grebel (2016) S-calibration.

4.4 Metallicity

The lack of detections in the auroral lines of [S iii]λ6312 and [N II]λ5755 implies that we cannot use the so-called direct method, which is based on the electron temperature T_e, to estimate the metallicity across our H II and DIG samples. Instead, we have to resort to the strong-line method based on [O iii] and [N II] that was first introduced by Alloin et al. (1979), which was later recalibrated by Pettini & Pagel (2004). We used the more recent recalibration by Marino et al. (2013). The O3N2 index is defined as O3N2 ≡ log₁₀ {([O iii]λ5007/Hβ)/([N II]λ6583/Hα)}, and is valid for the range −1 < O3N2 < 1.9. This is a suitable index for our data because NGC 300 is fairly metal rich in the low electron density regime, and all of these lines are well detected in our H II and DIG samples.

Figure 9 shows the resulting 12 + log₁₀ O/H as a function of galactocentric distance D for individual H II and DIG regions, as well as the uncertainty-weighted averages in bins containing 20 objects each. The O3N2 calibration by Marino et al. (2013) was obtained for H II regions and hence may not automatically apply for the DIG. However, by showing that chemically homogeneous H II-DIG pairs show minimum offset (0.01−0.04 dex) and little dispersion in the metallicity differences (0.05 dex) when using the O3N2 diagnostic, Kumari et al. (2019) concluded that the O3N2 metallicity calibration for H II regions can be used for DIG regions as well. A comparison of the H II and DIG region metallicity is therefore still meaningful, and we overplot the DIG results in Fig. 9.

We observe a flat metallicity gradient for D ≲ 245 arcsec ≈2.2 kpc. We have propagated the individual uncertainties of each line to the O3N2 index and then to the metallicity. Bresolin et al. (2009) suggested that strong-line diagnostics such as O3N2 (and N2) consistently underestimate the slope of abundance gradients because they ignore the effects of the ionization parameter, as extensively demonstrated by Pilyugin & Grebel (2016). Therefore, as a sanity check, we overplot the NGC 300 12 + log₁₀ O/H measurements by Bresolin et al. (2009) in the figure. These authors used auroral lines to obtain the H II electron temperature and then the metallicity via the so-called direct method. We note that our data show a flat gradient only in the central region, while the observed metallicity beyond D ≳ 400 arcsec ≈3.6 kpc seems consistent with the gradient suggested by the Bresolin data points at D ≳ 500 arcsec. Via single linear fit to their data, these authors found a y-intercept of 12 + log₁₀ O/H= 8.57 ± 0.02. A linear regression to our data gives a similar y-intercept of 12 + log₁₀O/H= 8.58 ± 0.01, which is consistent with the Bresolin value. Furthermore, we note that within the uncertainties shown in Fig. 9, our measurements are fully consistent with those of Bresolin at all galactocentric distances, not just in the y-intercept or at large distances. When it is superimposed on our data, the Bresolin et al. (2009) direct method metallicity measurements themselves appear to be also consistent with a flat inner disk gradient. Similarly, the results in Toribio San Cipriano et al. (2016) are consistent with those of Bresolin et al. (2009) in terms of oxygen abundance measurements of giant H II regions in NGC 300. We saw no indication of a steepening of the metallicity gradient in the inner disk, as suggested by Vila-Costas & Edmunds (1992). In both Vila-Costas & Edmunds (1992) and Bresolin et al. (2009), the number of data points near the galactic center (D < 1 kpc ≈110 arcsec) is much lower than in our data, which could explain the different conclusions reached by these authors. A flatter metallicity gradient for NGC 300 was also reported by Stasińska et al. (2013), who used the direct T_e method and 37 giant and compact H II regions, most of which are at comparable galactocentric distances as our sample (≤0.5R₂₅ in their Fig. 8 corresponds to ≤292.5 arcsec in our Fig. 9).

However, we must also point out that the intrinsic uncertainty in the O3N2 index is on the order of 0.02 dex (Marino et al. 2013). Such a large uncertainty range would easily be able to accommodate a metallicity gradient. The inherent insensitivity of the O3N2 index to variations in what Pilyugin & Grebel (2016) called the “excitation parameter P” is perhaps reason enough to dismiss the observed flat metallicity gradient in Fig. 9 as an artifact of the methodology. We refer to this parameter as P_O, where the subscript indicates that it is based on oxygen. It is possible that O3N2 consistently underestimates the slope only in the central disk regions (D ≲ 245 arcsec ≈2.2 kpc) where P_O might vary rapidly, and not at larger distances where it is fairly constant. We therefore also examined the resulting metallicity when the so-called S-calibration (Pilyugin & Grebel 2016) was used, based on R3 ≡ (I_{[O iii](λ4959+λ5007)}/I_Hβ), S2 ≡ (I_{[S II](λ6717+λ6731)}/I_Hβ), and N2 ≡ (I_{[N II](λ6548+λ6584)}/I_Hβ). This is shown in the lower panel of Fig. 9. Consistent with our previous findings, the H II regions continue to display a flat metallicity gradient in the central fields, albeit at ~0.2 dex lower metallicity. When the direct-T_e metallicity of Bresolin et al. (2009) is considered as the true metallicity, then the S-calibration systematically underestimates the true metallicity for H II regions. For the DIG, we now instead observe a steepening of the gradient at very small D ≲ 25 arcsec. In the outer fields at D ≳400 arcsec, the results from the S-calibration method suggest that the metallicity for both H II and DIG regions has dipped lower than for the inner regions, and steadily increases until it reaches the Bresolin et al. (2009) direct-T_e values for H II-regions. The A and B star metallicity measurements of Bresolin et al. (2009) suggest no such rapid increase in the outer fields, however. We summarize the obtained metallicity values per field in Table 2, where we list the formal uncertainties obtained through error propagation.

As another sanity check and in an attempt to understand why the H II metallicities from both methods show a similarly flat behavior although one method accounts for variations in P_O (the S-calibration) and the other does not (the O3N2 calibration), we constructed and examined the parameter P_O itself. Pilyugin & Grebel (2016) defined it as P_O ≡ R₃/(R₂ + R₃), where R2 ≡ (I_{[O II](λ3727+λ3729)}/I_Hβ), and R₃ as above. The MUSE data cover the [O II]λ7320, 30 doublet, but the detections are insignificant and the propagated uncertainties on P are too large, with relative errors ≳100%. The definition of P_O as given above traces the ionization more than the excitation, because the difference between ionization potentials⁴ of R2 and R3 is ΔE_ion = 21.5 eV, while the difference between collisional excitation potentials⁵ is only ΔE_exc = 0.8 eV. We can therefore attempt an approximation of P_O by using other ionic species of fairly large ionization potential differences, but only marginally different excitation potentials. The only option with our data is the [S iii]λ9068 and [S II]λ6716, 31 doublet, with ΔE_ion = 13 eV and ΔE_exc = 0.5 eV. Our modified parameter is then P_S ≡ S₃/(S₂ + S₃), where S3 ≡ (I_{[S iii](λ9068)}/I_Hβ), and S2 as above, and the subscript indicates that it is based on sulfur, not oxygen. The result is shown in Fig. 10. The P_s parameter for H II regions shows a flat monotone behavior with galactocentric distance D, which explains why the S calibration metallicity profile shows a flat gradient similar to that of the 03N2 method. The uncertainty-weighted average P_s parameter for all H II regions from all fields is 〈P_S〉 = 0.028 ± 0.002. At this low level, there is little to no difference between the P_O of Pilyugin & Grebel (2016) and our P_S variant. However, 90% of the data points in Fig. 10 lie below a value of P_S ≤ 0.15. We therefore proceed with this higher limit on P_S. To translate P_S values into P_O values, we used the Mexican million models database (3MdB; Morisset et al. 2015) to obtain a grid of Cloudy (v 17.02, Ferland et al. 2017) models. We selected only a subset of the 3MdB_17 model grid under the reference BOND_2 and based on the selection criteria of Amayo et al. (2021, their Eqs. (1)−(3)), describing the relations between the ionization parameter U and the ratio N/O to the oxygen abundance O/H. For the remaining set of models, we defined P_O ≡ [O iii] λ15007/([O iii] λ25007+ [O II]λ3726, 29) and P_S ≡[S iii]λ9068/([S iii]λ9068+ [S II]λ6716, 31) and plot the resulting models in Fig. 12. The figure shows that for P_S ≤ 0.15, the value of P_O can vary up to P_O ≲ 0.4. With this value, we can finally check the consistency of our metallicity results. A P_O ≲ 0.4 is indeed consistent with metallicities of 12 + log₁₀ O/H ~ 8.4−8.7 (Pilyugin & Grebel 2016, their Fig. 2), which comfortably brackets our results in Fig. 9.

An enlightening comparison could have been made for the DIG regions, but most have insignificantly detected [S iii] λ9068. The few available DIG data points at galactocentric distances D > 425 arcsec (orange diamonds in Fig. 10) are not enough to hazard a guess about radial variations in the DIG physical conditions.

In conclusion, we can note that a linear regression on the H II data has y-intercepts 12 + log₁₀ O/H = 8.58 ± 0.19 and 8.33 ± 0.06 with the O3N2 and S calibrations, respectively, where we have estimated the error via bootstrapping. Similarly, the DIG data has y-intercepts 12 + log₁₀ O/H = 8.51 ± 0.02 and 8.57 ± 0.15 with the O3N2 and S calibrations, respectively. A final comparison of our H II metallicity to the literature shows that the O3N2 calibration results are consistent with Stasińska et al. (2013), who obtained a y-intercept of 12 + log₁₀ O/H = 8.48 ± 0.03 for 37 compact and giant H II regions in NGC 300, based on the direct T_e method. Because two independent results from the literature are consistent with the O3N2 results in our sample, we consider only the O3N2-based metallicity for the remainder of this paper.

For completeness, in Fig. 11 we also show the 2D map of the O3N2-based metallicity. Taken at face value, the apparent decrease in metallicity throughout the outer fields p and q is consistent with expectations of lower metal abundances at larger radii in spiral galaxies (e.g., Holwerda et al. 2005).

Fig. 10 Excitation parameter P ≡ S3/(S2 + S3) as a function of galactocentric distance D.

Fig. 11 2D map of 12 + log10O/H from the 03N2 calibration for H II and DIG regions.

Fig. 12 3MdB models showing the relation between the PO parameter and our PS substitute. The limit of PS ≤ 0.15 and the corresponding PO ≲ 0.4 are indicated with solid lines. See Sect. 4.4 for more details.

4.5 Kinematics

Table 3 shows the kinematic properties of the fields for H II and DIG. For completeness, we also show the results of the stellar continuum. At the risk of stating the obvious, the local stellar background is best represented by itself and not by the residual of the local stellar background and the nearby stellar background, which is our method to obtain the H II region spectra before the pPXF line extraction. Therefore, we ran a separate pPXF extraction on spectra from H II region dendrogram leaves, but without subtraction of the diffuse background approximation. We stress that this pPXF run was solely used for the purpose of obtaining precise stellar continuum velocities V_cont and velocity dispersions σ_cont.

The focal point of our analysis, however, are the H II and DIG regions. We note that the values of σ_HII and σ_DIG fall below the instrumental resolution (47 km s⁻¹ at Hα), which means that it is unclear whether they are reliable. Similar values were obtained by den Brok et al. (2020) for 41 star-forming spirals in the MUSE Atlas of Discs, and these authors performed idealized simulations to examine the intrinsic compared to the measured dispersion of the gas. They concluded that although the absolute values of the dispersion are biased below σ = 25 km s⁻¹, the relative values are robust. Additionally, it has been demonstrated that pPXF robustly measures dispersion below the instrumental dispersion (Cappellari 2017). This enables us to perform a comparative analysis of the H II and DIG regions.

If the velocities V and velocity dispersions σ were weighted by the Hα flux, then the handful of bright regions in our sample would dominate the results. This is not representative of the sample, and therefore we take the uncertainty-weighted averages instead in Table 3. The lowest uncertainty points have the highest weights. The shear velocity is a measure of the large-scale motion of the gas, and we measured it as V_shear = (V_max – V_min)/2. Similar to Herenz et al. (2016) and Micheva et al. (2019), we took the median of the upper and lower fifth percentile of the velocities to obtain V_min and V_max, and propagated the full width of each percentile as a conservative estimate of the V_shear uncertainties. We note that the systemic stellar velocity we obtain (field A in Table 3) is $V_{\text{sys}}^A$ = 144.3 ± 0.4 km s⁻¹, which agrees well with the literature, for instance, Rogstad et al. (1979, V_sys = 145 ± 2 km s⁻¹)and Westmeier et al. (2011, V_sys = 144 ± 0.2 from radio data).

The table is visualized in Fig. 13, showing the velocity, velocity dispersion, the V_shear/σ ratio, and the residual between H II and DIG gas kinematics in panels a, b, c, and d, respectively. Figure 13a shows that the H II and DIG velocities are similar, implying that the two components move together. This is perhaps better visible in panel d, where ΔV is fairly flat for the inner fields, with the DIG becoming marginally faster than the H II gas in the outer fields P and Q. In contrast, Fig. 13b shows an increased velocity dispersion for the DIG in the central regions, with a decreasing tendency for larger galactocentric distances, while the H II gas apparently maintains a constants σ. This results in an increase in Δσ with galactocentric distance in panel d. Figure 13c shows that V_shear/〈σ〉 has a flat behavior that remains fairly constant within the uncertainties for both H II and DIG. Taking the V_shear/〈σ〉_H II ~ 1.0 average at face value this implies that H II gas is rotationally supported, while the DIG is predominantly dispersion supported, with Vshear/〈σ〉_DIG ~ 0.4. When the propagated uncertainties are considered as well, the separation between rotationally and dispersion supported kinematics of the H II and DIG regions is less clear, but still systematic at all distances. The exception is field J, which covers the inter-arm region and features high extinction from a prominent dust lane (see Fig. 1 in Roth et al. 2018). It is possible that the inter-arm region has different DIG kinematics, which we sample only with field J and partially with field E, which apparently covers a fraction of the spiral arm and also another inter-arm region with a visible dust lane (Roth et al. 2018).

In Table 3, the median σ_HII = 21.2 ± 0.2 km s⁻¹, and σ_DIG ~ 25.1 ± 0.1 km s⁻¹. The DIG velocity dispersion is higher than for H II gas, which is consistent with the DIG being hotter than H II gas in spiral galaxies, including the Milky Way (e.g., Madsen et al. 2006; Haffner et al. 2009; Della Bruna et al. 2020). The expected temperature difference in the literature is ΔT_e ~ 2–4 kK at most. We measured a difference in average velocity dispersion of Δσ ~ 3.9km s⁻¹. When we assume that all of this difference is due to the thermal component of the velocity dispersion, σ_th, we can calculate the corresponding electron temperature as $T_{\text{e}}\left[K\right]={\sigma }_{\text{th}}^2{m}_{\text{H}}/k_{\text{B}}$, where m_H and k_B are the mass of a hydrogen atom and the Boltzmann constant, respectively. With Δσ ~ 3.9 km s⁻¹, we obtain T_e ~ 1.8 kK, which is well within the range observed in the literature to date. The difference per field is shown in Fig. 13e, where we have again assumed that all of the velocity dispersion is purely thermal, and show the corresponding difference in temperature between H II and DIG. It is clear from the figure that at galactocentric distance of D ≳ 75 arcsec all of the observed velocity dispersion difference between H II and DIG could be accounted for simply by temperature differences. We note that this does not mean that turbulence is negligible in these fields because the DIG is heterogeneous (e.g., Haffner et al. 2009; Madsen et al. 2006) and may be hotter than or of equal temperature as the H II gas. The implied temperature differences are consistent with the velocity dispersion being dominated by its thermal component. For the innermost regions A and B, however, the figure suggests that the velocity dispersion cannot be exclusively thermal because that would correspond to a ΔT_e = 6–10 kK, that is, the DIG would have to be 6–10 kK hotter than the H II gas, which is much higher than the expected temperature differences found in the literature. We therefore conclude that in these central regions, the DIG must be significantly more turbulent than the H II gas and the DIG in the rest of the fields. Alternatively, similar to the spirals in den Brok et al. (2020), the DIG emission might partially originate from a somewhat thicker layer than the H II gas.

In Appendix C, we show the complementary figure for the stellar kinematics. We do not dwell on the stars, other than to note that the stellar continuum is always slower than both the H II and DIG regions, and always with a higher velocity dispersion. This observation is consistent with the early-to-intermediate type spirals in Vega Beltrán et al. (2001), for example. As these authors noted, these differences in kinematics between the gas and the stars can be explained by the gas being confined to the disk and supported by rotation, as implied by the measured V_shear/〈σ〉_H II ~ 1.0, while the observed stellar continuum mostly belongs to the bulge and is dispersion supported, with V_shear/〈σ〉_stars ~ 0.3.

Two-dimensional velocity maps of the stellar continuum, H II, and DIG regions are presented in Fig. 14, and the corresponding velocity dispersion maps are shown in Fig. 15. The data trace the velocity field of the galaxy. Most of the pointings show a velocity redshift from the central field a, with the inter-arm field J showing a mix of red- and blueshifted regions. The line of midpoint velocity, defined on stars in field A, is clearly visible for the stellar component (white in Fig. 14, top panel), becomes diluted for the H II regions (Fig. 14, middle panel), and appears to be offset in the DIG (Fig. 14, bottom panel). In Fig. 15, the velocity dispersion shows an increase toward the outskirts (fields P, Q) for stars (top panel) and a reversed trend for the DIG (bottom panel). We note that the velocity dispersion in the inter-arm field J is lower for both stars and DIG. The H II regions (middle panel) show a more complex velocity dispersion behavior that is apparently uncoupled from the distance of the pointing to the galactic center in field A.

Fig. 13 Sample kinematics. (a) Velocity. (b) Velocity dispersion. (c) Vshear/σ. (d) ΔV and Δσ between Hii and DIG. (e) ΔTe between H II and DIG, assuming Δσ is solely due to the thermal component. The pink shaded region is the ΔTe range found in the literature. Panels (b) and (c) have the same legend as in (a). All panels share the galactocentric distance in arcsec as x-axis. The twin axis shows the corresponding field labels.

Fig. 14 Velocity in km s−1 for stars (top), H II regions (middle), and DIG (bottom). We note that the color map is a divergent two-slope color map, linear in both slopes around the centered midpoint (white color) at 144.34 km s−1 in all panels, which is the uncertainty-weighted average velocity of the stars in field A.

5 Stacked spectra

To increase the S/N in the spectra, we average-stacked all H ii-region spaxels per field and ran pPXF on the stacks to extract nebular emission lines and kinematics. To obtain a galaxy stack, we also stacked all spaxels from all fields. The same method was used to obtain a per field and galaxy stack for the DIG regions. These spectra are shown in Fig. D.1.

In the stacked spectra of non-AO fields, we identify the faint auroral [N II]λ5755 line for both H II and DIG stacks. However, the propagated uncertainties in the temperature-sensitive [N II]λ5755/[N II]λ6548 ratio are too large for meaningful results of the electron temperature and density. For example, for field J, the only field in which the combination of [N II] and [S II] ratios is strictly within the theoretical limit, we obtain $\left(T_{\text{e}},N_{\text{e}}\right)=({9178}_{-3223}^{+2329}$ K, ${1.9}_{-1.9}^{+14.8}$ cm⁻³) for H II regions and $\left(T_{\text{e}},N_{\text{e}}\right)=({15611}_{-8788}^{+9555}$ K, ${1.2}_{-12.7}^{+42.3}$ cm⁻³) for DIG regions. This implies that T_e for the DIG could be 1 kK to over 10 kK hotter than for the H II regions. This result is hardly constraining. For the rest of the non-AO fields, we were able to obtain upper limits on T_e within the uncertainties at best. We list them in Table 4. Within the uncertainties, the H II region electron densities in field J are indistinguishable from those of the DIG. Taken at face value, even within the uncertainties, the N_e for H II regions in field J are at the low end of the expected densities (e.g., Osterbrock & Ferland 2006, N_e ∈ (10, 1000) cm⁻³), which is consistent with our sample containing extremely faint H II regions.

From these stacked spectra, we calculated the DIG/Hα fractions per field. They are listed in Table 5. In order to test the robustness of these fractions, we checked how the sizes of the H II regions and the physical locations of the DIG regions affect the measured DIG contribution. For this purpose, we expanded the H II region masks by the gap annulus of five spaxels and simultaneously shifted the DIG annulus by the same amount away from the H II regions. This resulted in an increasing number of H II spaxels by a factor of 2.0, 2.0, 2.0, 2.15, 2.66, 2.48, 2.38, 1.97, and 2.68 for fields A, B, C, D, E, I, J, P, and Q, respectively. Simultaneously, the number of DIG spaxels increased by a factor of 1.24, 1.24, 1.24, 1.16, 1.22, 1.34, 1.20, 1.22, and 1.13, respectively. This affects the DIG/H II fractions by <10%, namely, by +4.1%, +4.1%, +4.1%, +9.3%, −4.3%, +5.9%, +2.4%, +5.7%, and +0.7%, respectively. We therefore conclude that our definition of H II and DIG regions is robust, and that the DIG fractions are not overly sensitive to the exact shape of our H II regions. Doubling the H II region size does not influence the conclusions of our analysis.

Fig. 15 Velocity dispersion in km s−1 for stars (top), H II regions (middle), and DIG (bottom).

6 Discussion

In Sect. 4, we presented a number of properties that can be extracted from our data for the H II and DIG regions. We obtained no convincing differences between the two types of objects in terms of metallicity, velocity, and velocity dispersion, and marginal differences in electron density and extinction. The DIG apparently has a lower density and E(B−V) values on average. The DIG fraction we measure per field appears to agree in general with expectations from the literature (e.g., Oey et al. 2007; Chevance et al. 2020; Belfiore et al. 2022; Della Bruna et al. 2020), as shown in Table 5.

Belfiore et al. (2022) reported that the main differences between the DIG and H II regions in spiral galaxies are enhanced low-ionization line ratios toward the centers of spiral galaxies, an offset of some DIG regions from the locus of H II regions in the BPT diagram extending beyond the star-formation region and into the Low-ionization nuclear emission-line/Active galactic nuclei (LINER/AGN) domain, and a decreasing [O iii]λ5007/Hβ ratio with the Hα surface brightness of the DIG. Their H II and DIG regions, however, are brighter than ours, with log L(Hα[erg s⁻¹]) ~ 36. Our regions are much fainter, with a median log L(Hα[erg s⁻¹]) = 34.7, and thus fall in a brightness range that has been studied much less extensively. It is therefore prudent to verify whether these trends are still present at these extremely faint levels. In what follows, we have found no indication that the extremely faint, manually added large DIG regions described in Sect. 2 and shown in Fig. 1 behave differently than normal DIG regions created by binary dilation of the H II region masks.

In Fig. 16 we show the BPT diagram for DIG and H II field stacks, as well as the Kewley et al. (2006) limit for star-forming regions. In this parameter space, the H II regions are separated from the corresponding DIG in the field in all fields. Similar to the findings in Belfiore et al. (2022), the DIG is found predominantly in the top right corner of the diagram, crossing into the region associated mainly with LINERs and AGN. This is likely due to the increasing contribution to DIG ionization from HOLMES with a harder spectrum than H II regions (Flores-Fajardo et al. 2011; Belfiore et al. 2022). In addition, the DIG ionization state is expected to be lower than that of the H II regions (Madsen et al. 2006), which we also observe in the middle panel of Fig. 16, using the [S iii]λ9068/[S II]λ6716 + 31 ratio as a tracer of the ionization parameter. The same figure demonstrates that the DIG high-to-low ionization line ratio is flat within the uncertainties (fields A-J, Q) or decreasing (fields P and, e.g., B) with increasing Hα surface brightness. Either of these trends requires the added emission from HOLMES (Flores-Fajardo et al. 2011; Belfiore et al. 2022). The final argument in support of the HOLMES contribution to DIG ionization is the behavior of low-ionization line ratios such as [S II]λ6716/Hα and [N II]λ6583/Hα with galactocentric distance. For the DIG, an increase in these ratios is observed with increasing distance to the galactic plane in an edge-on spiral (Flores-Fajardo et al. 2011) and with galactocentric distance in 19 low-inclination spirals from the PHANGS-MUSE survey (Belfiore et al. 2022). The bottom panel of Fig. 16 shows the expected systematic enhancement toward small distances in both of these line ratios in the DIG, but not in the H II regions.

Despite any contribution of HOLMES to the ionization of the DIG, it is likely that the DIG line ratios are further enhanced by contribution from shocks due to supernova explosions. There is a positive correlation between the velocity dispersion and the shock velocity in a scenario in which multiple shocks propagate in random directions (Ho et al. 2014). In a BPT diagram, this correlation translates into an increased velocity dispersion with the [N II]λ6583/Hα or [S II]λ6716/Hα line ratios. Similar to Oparin & Moiseev (2018) and Della Bruna et al. (2020), we construct in Fig. 17 the BPT diagram for the DIG and color-code the data points with their gas velocity dispersion. We observe an increase in velocity dispersion toward the top right corner of the diagram, in the LINER region. Overplotted are the MAPPINGS IV shock models of Ho et al. (2014) for different shock fractions (20–100%). The line segments connect model grids of different shock velocities. All of the measured DIG velocity dispersions in any field (Table 3) are well below the shock velocities of the models (100, 200, and 300 km s⁻¹), but we note that we did not account for multiple components in the line extraction with pPXF. To investigate the effect of fitting a single line to a narrow and broad line composite velocity dispersion, we performed the following test. We created a narrow line with a velocity dispersion σ = 22 km s⁻¹ and three broad lines of 100, 200, and 300 km s⁻¹ each. Then we created composite lines by varying the contribution of the broad component from 5 to 100% of the strength of the narrow line. Finally, we extracted the lines with pPXF and examined the recovered velocity dispersion for each composite line as a function of the fractional contribution of the broad-line component. This is visualized in Fig. 18, where we plot the residual of the velocity dispersion of the composite line and the narrow-line baseline (22 km s⁻¹). For a 100 km s⁻¹ broad component, up to 20% of the contribution remains completely undetected by pPXF within the uncertainties, and it only increases the recovered velocity dispersion by ~10 km s⁻¹ up to ~90% of the contribution by the broad component. This means that even if significantly strong broad components with a velocity dispersion of 100 km s⁻¹ exist, we would not detect them, and would instead measure a single line of a velocity dispersion that is marginally broadened by ~10 km s⁻¹. Similar arguments can be made for the 200 and 300 km s⁻¹ shocks, but with a much smaller contribution range of up to ~20% of the line strength. We therefore cannot exclude the possibility that shocks of these velocities exist in our data.

By comparing Fig. 17 to the upper panel of Fig. 16, we see that the shock models clearly predict a substantial contribution of shock ionization (0.2–0.4) for the DIG regions in fields C, B, D, and I, and a dominant contribution (~0.6) in fields A and J. Field A is the center of the galaxy and contains several morphological features reminiscent of ancient SNR arcs; this has been reported in Roth et al. (2018). This field also contains several X-ray emission sources including a black hole X-ray binary (Read & Pietsch 2001) and a bipolar microquasar jet (McLeod et al. 2019), and hence we intuitively expect field A to be kinematically active. Field J is the inter-arm region, in which a shock can form along the trailing edges of a spiral arm (e.g., Foyle et al. 2010). Fields P and Q are at the largest galactocentric distances in our data and are apparently consistent with a zero shock fraction contribution.

Finally, in Appendix E, we determine whether machine learning can detect the subtle differences between H II and DIG regions. However, the two object types completely overlap in the UMAP 2D projection, implying that the H II (DIG) regions show larger differences within themselves than when compared to the DIG (H II) regions.

Fig. 16 Evidence of HOLMES. Top: BPT diagram for average stacked spectra per field. Middle: high-to-low sulfur ionization ratio as a function of Hα surface brightness for average stacked spectra. Bottom: low-ionization ratios of sulfur and nitrogen with galactocentric distance for individual H II and DIG regions. The total spectrum of stacking all spaxels from all fields is marked with a star.

Fig. 17 BPT-σ diagram for DIG regions. The shock models of Ho et al. (2014) are overplotted for different shock velocities. The line segments are color-coded by shock fractions, as indicated in the legend.

Fig. 18 pPXF recovered velocity dispersion of composite lines and a single narrow line. The x-axis represents the per cent increase of the strength of the broad component relative to the narrow line. The gray-shaded horizontal area demarcates the average uncertainty in the DIG velocity dispersion measurement (single narrow line results).

7 Conclusions

We presented MUSE data of nine pointings in the flocculent spiral galaxy NGC 300. For all fields, we identified emission line regions via a dendrogram algorithm on preliminary Hα images. We extracted nebular emission lines such as the Banner and Paschen lines, and the strongest forbidden lines such as [S II]λ6716, [S iii]λ9068, [O iii]λ5007, and [N II]λ6583 via pPXF. We proceeded to separate the detections into H II, SNR, and PNe classes. All extracted properties are provided in a catalog of emission line objects. We also extracted DIG regions via dilated H II-region masks. The H II sample is extremely faint, with log₁₀ L_Hα[erg/s] ∈ (33.2, 37.3), and a median of 34.7. This makes our sample one of the faintest H II samples to date. The DIG is similarly faint, with log₁₀ L_Hα[erg/s] ∈ (33.6, 36.5), and a median of 34.9. We examined the properties of the H II and DIG regions in terms of kinematics, abundances, density, and extinction. Our conclusions from this analysis are listed below.

1. The distribution of [S iii]λ6731/16 ratios for the H II and DIG regions peak at the extreme low-density limit. For a handful of objects, we calculated a range of possible electron densities and found averages of 100 and ~23 cm⁻³ for the H II and DIG regions, respectively.

2. Most of the DIG is consistent with no extinction, E(B − V) = 0. The H II average is E(B − V) = 0.1. For the handful of DIG objects with nonzero extinction, we find 〈E(B − V)〉 = 0.05 on average.

3. The 03N2-calibration of the strong line method gives oxygen abundances 12 + log O/H consistent with those of the direct-T_e method, while the S-calibration on average underestimates the oxygen metallicity. In the inner fields, the average metallicity shows a completely flat profile with galactocentric distance, without evidence of a metallicity gradient or increase in the central region. The metallicity of the DIG, 8.48 on average, is consistent with that of the H II regions, 8.53 on average, at any galactocentric distance.

4. The H II and DIG regions move together with similar velocities. The DIG becomes marginally faster than the H II gas in the outer fields.

5. The average velocity dispersion is 21.2 ± 0.2km s⁻¹ for H II gas and 25.1 ± 0.1 km s⁻¹ for DIG. This difference can be fully attributed to the thermal velocity dispersion component, which would correspond to the DIG being 1.8kK hotter than the H II regions.

6. The DIG has an increased velocity dispersion in the central galactic region, consistent with models of a dominant (~60%) contribution of shocks to the DIG ionization.

7. The DIG fraction per field varies between 42–77% of Hα. The inter-arm region field J shows a much lower DIG fraction of 15%. These numbers are consistent with the literature.

8. The DIG has a lower ionization state than H II gas, as traced by the high-to-low ionization line ratio [S iii]λ9068/[S II]λ16716 + 31.

9. Signs of a contribution to DIG ionization by hot low-mass evolved stars are detected:

   • (i)

     a completely flat trend of the DIG [S iii]λ9068/[S II]λ6716 + 31 ratio with Hα surface brightness, in contrast to a positive correlation for H II regions,

   • (ii)

     low ionization line ratios such as [S II]λ6716/Hα and [N II]λ6583/Hα show a systematic enhancement toward small galactocentric distances, in contrast to a flat trend for H II regions.

10. Unsupervised machine-learning algorithms such as UMAP/HDBscan are unable to distinguish between DIG and H II regions, implying that both the DIG and the H II regions are so heterogeneous that the differences within them are larger than between them.

11. The differences between extremely faint H II and DIG regions follow the same trends as their brighter counterparts.

On a final note, we acknowledge that we have completely omitted any comparison of the number of available ionizing photons and the necessary photons that the DIG regions imply. This is beyond the scope of this paper and will be addressed in a separate paper that is currently in preparation.

Acknowledgements

The authors would like to thank the anonymous referee for the careful reading of the paper. P.M.W. received support from BMBF Verbundforschung (ErUM, project VLT-BlueMUSE, grant 05A20BAB). N.C. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) – CA 2551/1-1. CM. acknowledges support from grant UNAM / PAPIIT - IN101220.

Appendix A Attempt of classifying emission line objects via machine learning

Machine learning has been developed to facilitate the analysis of large data, and we wish to test whether it can help classify the objects into H II, PNe, and SNR. An unsupervised machine-learning classification algorithm is apparently most suitable for this task; t-SNE⁶ (van der Maaten & Hinton 2008) and UMAP⁷ (McInnes et al. 2018) are used widely. They reduce the dimensionality of the data and can project the data in 2D, clustering points of significantly similar features. We selected UMAP because it appears to have advantages over t-SNE in the selection of its cost function, which better preserves global structure in the projection.

The choice of features for the machine-learning algorithm is more important than the selection of the algorithm. The features have to be relevant to the desired classification outcome. We tested several sets of features and list them in Table A.1. It is possible to add more features and explore subgroups of each object class (H II, PNe, and SNR), but the more features are added, the harder the interpretation of the resulting clustering projection. Another important step is the normalization of the features by removing the mean and scaling by the variance. This is necessary in order to ensure that all features are within the same numerical range and hence have an equal probability to influence the learning algorithm. A problem we immediately ran into is that emission line fluxes are not independent features, and hence scaling along columns does not seem suitable. We have explored both options, scaling per row, that is, scaling each individual object by some feature, and scaling per column, that is, scaling one feature across all objects.

The most important UMAP parameter is n_neighbours, which determines the number of neighboring points that influence the clustering of dendrograms in the 2D projection. We experimented with n_neighbours= 5, 6, 7, 8, 9, 10, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23, 35, 70, 100, and 300. Increasing the value of n_neighbours resulted in an ever more tightly clustered projection. For display purposes, we chose n_neighbours= 23, but we verified that all 20 projections give very similar results in terms of cluster identification and classification. In addition, UMAP is a stochastic algorithm in the way it optimizes and speeds up the various approximation steps. Hence, different runs of UMAP can produce slightly different results, depending on the random state that is selected. To ensure reproducibility, we set the random seed to 42 in this section. For all 20 values of n_neighbours, we verified that changing the random state to one, for example, gives very similar results and does not change the outcome. Finally, we set min_dist to zero, meaning that the points were allowed to cluster so tightly in the 2D projection that they completely overlapped, and spread to 5, which increases the distance between clusters but does not dilute the clustering itself.

The UMAP results for n_neighbours = 20 and random_seed = 1 are shown in Figure A.1. No clear separation of H II, SNR, and PNe is obtained for any feature set. At best, the run with two features pushes the SNR and PNe towards the edges of the 2D projection, as shown in the middle and right panels in the first row of Figure A.1. This is not reflected in the clustering found by HDBscan in the first panel, however. Instead, all three classes of objects in the different subclusters have identical labels. The six-feature run shows two well-defined clusters in the HDBscan panel, which is the effect of adding the galactocentric distance as a feature in this run. The smaller cluster contains fields P and Q, and the larger cluster contains all remaining fields. We conclude that unsupervised machine-learning algorithms such as UMAP and t-SNE are not suitable for classifying emission line objects into H II, SNR, and PNe.

Fig. A.1 Sample classification via machine learning. (Left) UMAP 2D projection with HDBscan cluster labels. (Middle) UMAP 2D projection with [N II]λ6583-based BPT diagram labels. (Right) UMAP 2D projection with [S II]λ6716 + 31-based BPT diagram labels. The inset text shows the number of features in each run. All displayed runs are with n_neighbours= 20.

Appendix B Justification of the [S II]λ6731/16 cutoff

The theoretical limits of [S II]λ6731/16 are 0.6806 ≤ [S II]λ6731/16≤ 2.351. Simply due to statistics, we would expect there to be a distribution of [S II]λ6731/16 measurements around some characteristic (typical) value. In NGC 300, the characteristic value is at 0.675 ± 0.036. This poses a problem because it is located practically at the theoretical limit, and we now have to decide how many of the data points falling to the left of the limit need to be included in the sample.

When the line fluxes of the [S II]λ6716, 31 individually follow normal distributions, their ratio is not necessarily a Gaussian, although it can be approximated by a Gaussian if the coefficient of variation CV= σ_gauss/µ_gauss is sufficiently small (Díaz-Francés & Rubio 2013, CV< 0.1). The top panel of Figure B.1 shows a histogram of the [S II]λ6731/16 with propagated uncertainties of the data points falling inside each bin. We used astropy to fit a 1D Gaussian to the data points and obtained σ_gauss = 0.06 and a peak location µ_gauss = 0.66. This gives CV= 0.09, and hence we can safely continue to approximate the distribution of [S II]λ6731/16 data points with a Gaussian.

To the left of the Gaussian peak, up to 34.1%, 13.6%, and 2.1% of the data points are expected to be found in the three intervals (−1σ, peak), (−2σ, −1σ), and (−3σ, −2σ), respectively. For the background-subtracted spectra, we calculated the fractions of data points in these intervals as 14.8%, 7.3%, and 3.7%, respectively. The first two intervals are consistent with the data being normally distributed, and hence we kept these points in the final sample. In contrast, we expect only 2.1% of the data to fall between −3σ and −2σ, but instead, we measure 3.7%, implying that in addition to statistical scatter there is also some noisy data with low S/N. In summary, we discarded all background-subtracted spectra that were associated with [S II]λ6731/16 more than −2σ away from the peak. This corresponds to a limit of [S II]λ6731/16≤ 0.548. In a similar fashion, the bottom panel of Figure B.1 shows the histogram and fit for the DIG data. We discarded all DIG spectra with [S II]λ6731/16=≤ 0.638.

We note that we cannot apply a cutoff on the right side of the Gaussian peak because in addition to statistics, the physical structures present in the NGC 300 fields affect the distribution.

Fig. B.1 Justification of the [S II]λ6731/16 cutoff. (Top) Histogram of the [S II]λ6731/16 for background-subtracted regions (gray markers with error bars), with the best-fit Gaussian model (solid blue line) and vertical lines indicating the location of −1σ (dotted line), −2σ (dash-dotted), and −3σ (dashed) plotted for convenience. The uncertainties are propagated from the points in each bin. (Bottom) Same as top panel, but for DIG regions.

Appendix C Stellar kinematics

For completeness, we present the stellar kinematics in Figure C.1.

Fig. C.1 Sample stellar kinematics, (a) Velocity, (b) Velocity dispersion, (c) Vshear/σ for the stellar continuum and the DIG. Panels b) and c) have the same legend as panel a). All panels share the galactocentric distance in arcsec as x-axis. The twin axis shows the corresponding field labels.

Appendix D Stacked spectra

The average stacked spectra per field are shown in Figure D.1. For both H II and DIG spectra, only the residuals of the original spectra minus the pPXF best-fit stellar continuum are shown. For all fields, the figures are capped at half of the maximum Hα peak value in order to provide more detail for the fainter lines. We note that the DIG contribution and the initial stellar background are removed from the H II regions already before the pPXF fit, as explained in Section 2.3.

Fig. D.1 Residual average stacked spectra for H II (blue) and DIG (orange) for each field. The field is indicated in the top right corner of each panel. The main nebular emission lines that are of interest in this work are annotated in the top panel. Bad pixels that were omitted in the pPXF fit are grayed out. The y-axis has units of 10−20erg/s/cm2/Å and is truncated at half the Hα flux.

Appendix E Attempt to separate DIG and H II regions via machine learning

With the same method as in Section A, but for the feature sets in Table E.1, we determined whether the UMAP projection separates the DIG and H II regions. The result is shown in Figure E.1.

The two types of objects overlap in all runs, and moreover, the labeling of the HDBscan algorithm implies larger differences within different subgroups in the DIG or the H II regions than between the DIG and H II regions as a whole.

Fig. E.1 Sample classification in to H II and DIG via machine learning. (Left) UMAP 2D projection with HDBscan cluster labels. (Middle) UMAP 2D projection with DIG/H II labels. (Right) UMAP 2D projection with field labels. The inset text shows the number of features in each run. All displayed runs are with n_neighbours= 20.

References

All Tables

All Figures

	Fig. 1 Region definitions for an example region in field P. Overplotted are (a) dendrogram leaves, (b) DIG ring regions, which also serve as sky background regions, and (c) dendrogram leaves and DIG rings both, as well as the manually added DIG regions, shaded in pink.
In the text

Fig. 2 BPT diagram in [N II]λ6584/Hα (top) and [S II]λ6716, 31/Hα (bottom). The objects are classified in the [N II]λ6584/Hα diagram via the Kewley et al. (2006; solid line) star formation limit and the PNe-SNR (dashed line) separation lines. The Dopita et al. (2000; dotted line) is also shown for reference. The color-coding in the [S II]λ6716, 31/Hα diagram is based on the [N II]λ6584/Hα classification.

In the text

	Fig. 3 2D map of the BPT classification.
In the text

	Fig. 4 Normalized histograms of (a) the Hα surface brightness ΣHα (10−20 erg s−1 cm−2 arcsec−2) of the H II and DIG dendrogram regions, (b) the equivalent radius (pc), and (c) the Hα luminosity (erg s−1).
In the text

	Fig. 5 Hα luminosity (erg s−1) as a function of region size in spaxels (bottom x-axis) or equivalent radius in pc (top x-axis).
In the text

Fig. 6 Sample electron densities. Top: histogram of the [S II]λ6731/16 ratio for H II (solid blue lines) and DIG (dashed orange lines). The area hatched with backslashes represents non-physical values, and the forward slash covers the range of Ne that is fairly invariant with Te. Bottom: [S II]λ6731/16 emissivity grid for 3kK ≤ Te ≤ 20 kK and 1 ≤ log10 Ne(cm−3) ≤ 7. Contours of [S II]λ6731/16 are overplotted in white, and the corresponding value is shown inline.

In the text

	Fig. 7 Electron density Ne as a function of galactocentric distance D for H II (blue) and DIG regions (orange). At each D, the width of the Ne range is shown, obtained from Ne(Te= 20 kK) and Ne(Te= 4 kK). The medians for H II and DIG are also shown with solid lines.
In the text

Fig. 8 Sample nebular extinction. Top: E(B − V) as a function of galactocentric distance D for individual H II (black dots) and DIG regions (blue diamonds). Uncertainty-weighted averages in bins of 20 objects are also shown in red circles for H II and in orange diamonds for DIG. The propagated error bars are often smaller than the markers. The x-axis position of each uncertainty-weighted average is the average galactocentric distance for the dendrograms inside each bin. Bottom: 2D map of E(B − V).

In the text

	Fig. 9 Sample metallicity. Top: 12 + log10 O/H as a function of galactocentric distance D, using the Marino et al. (2013) O3N2 calibration. Same legend as in Fig. 8, except for the data by Bresolin et al. (2009; open gray circles), which were obtained via the direct-Te method for NGC 300. Bottom: same as top, but using the Pilyugin & Grebel (2016) S-calibration.
In the text

	Fig. 10 Excitation parameter P ≡ S3/(S2 + S3) as a function of galactocentric distance D.
In the text

	Fig. 11 2D map of 12 + log10O/H from the 03N2 calibration for H II and DIG regions.
In the text

	Fig. 12 3MdB models showing the relation between the PO parameter and our PS substitute. The limit of PS ≤ 0.15 and the corresponding PO ≲ 0.4 are indicated with solid lines. See Sect. 4.4 for more details.
In the text

Fig. 13 Sample kinematics. (a) Velocity. (b) Velocity dispersion. (c) Vshear/σ. (d) ΔV and Δσ between Hii and DIG. (e) ΔTe between H II and DIG, assuming Δσ is solely due to the thermal component. The pink shaded region is the ΔTe range found in the literature. Panels (b) and (c) have the same legend as in (a). All panels share the galactocentric distance in arcsec as x-axis. The twin axis shows the corresponding field labels.

In the text

	Fig. 14 Velocity in km s−1 for stars (top), H II regions (middle), and DIG (bottom). We note that the color map is a divergent two-slope color map, linear in both slopes around the centered midpoint (white color) at 144.34 km s−1 in all panels, which is the uncertainty-weighted average velocity of the stars in field A.
In the text

	Fig. 15 Velocity dispersion in km s−1 for stars (top), H II regions (middle), and DIG (bottom).
In the text

	Fig. 16 Evidence of HOLMES. Top: BPT diagram for average stacked spectra per field. Middle: high-to-low sulfur ionization ratio as a function of Hα surface brightness for average stacked spectra. Bottom: low-ionization ratios of sulfur and nitrogen with galactocentric distance for individual H II and DIG regions. The total spectrum of stacking all spaxels from all fields is marked with a star.
In the text

	Fig. 17 BPT-σ diagram for DIG regions. The shock models of Ho et al. (2014) are overplotted for different shock velocities. The line segments are color-coded by shock fractions, as indicated in the legend.
In the text

	Fig. 18 pPXF recovered velocity dispersion of composite lines and a single narrow line. The x-axis represents the per cent increase of the strength of the broad component relative to the narrow line. The gray-shaded horizontal area demarcates the average uncertainty in the DIG velocity dispersion measurement (single narrow line results).
In the text

	Fig. A.1 Sample classification via machine learning. (Left) UMAP 2D projection with HDBscan cluster labels. (Middle) UMAP 2D projection with [N II]λ6583-based BPT diagram labels. (Right) UMAP 2D projection with [S II]λ6716 + 31-based BPT diagram labels. The inset text shows the number of features in each run. All displayed runs are with n_neighbours= 20.
In the text

Fig. B.1 Justification of the [S II]λ6731/16 cutoff. (Top) Histogram of the [S II]λ6731/16 for background-subtracted regions (gray markers with error bars), with the best-fit Gaussian model (solid blue line) and vertical lines indicating the location of −1σ (dotted line), −2σ (dash-dotted), and −3σ (dashed) plotted for convenience. The uncertainties are propagated from the points in each bin. (Bottom) Same as top panel, but for DIG regions.

In the text

	Fig. C.1 Sample stellar kinematics, (a) Velocity, (b) Velocity dispersion, (c) Vshear/σ for the stellar continuum and the DIG. Panels b) and c) have the same legend as panel a). All panels share the galactocentric distance in arcsec as x-axis. The twin axis shows the corresponding field labels.
In the text

	Fig. D.1 Residual average stacked spectra for H II (blue) and DIG (orange) for each field. The field is indicated in the top right corner of each panel. The main nebular emission lines that are of interest in this work are annotated in the top panel. Bad pixels that were omitted in the pPXF fit are grayed out. The y-axis has units of 10−20erg/s/cm2/Å and is truncated at half the Hα flux.
In the text

	Fig. E.1 Sample classification in to H II and DIG via machine learning. (Left) UMAP 2D projection with HDBscan cluster labels. (Middle) UMAP 2D projection with DIG/H II labels. (Right) UMAP 2D projection with field labels. The inset text shows the number of features in each run. All displayed runs are with n_neighbours= 20.
In the text