© The Authors 2024

1. Introduction

Understanding the relationship between stellar formation processes and the interstellar medium (ISM) is a key step in discerning the complexity in the evolutionary history of galaxies. One problem in this regard has been understanding the nature and importance of feedback processes in which massive stars deposit energy into the interstellar medium through photoionisation, stellar winds, and supernovae. This feedback mechanism affects the physical and dynamic state of the ISM and, therefore, it has an influence on the rate and distribution of stellar formation in galaxies (Ceverino & Klypin 2008, 2009; Hopkins et al. 2014; Klessen et al. 2016; Gatto et al. 2017; Grisdale 2017).

In this context, the existence of a warm and ionised component of the ISM that is ubiquitously distributed in our galaxy has been known for decades. This component of the Milky Way has been referred to as the warm ionised medium (WIM; Hoyle & Ellis 1963; Hewish et al. 1968; Reynolds 1971, 1985, 1989; Reynolds et al. 1973; Guélin 1974; Kennicutt et al. 1989; Finkbeiner 2003). With the warm temperatures and low densities found in this medium (0.6 − 1 × 10⁴ K and ∼10⁻¹ cm⁻³ respectively, with an emission measurement of ∼35 cm⁻⁶ pc), the WIM represents the 20% of the ISM in volume and the 90% of the ionising hydrogen in mass, being one of the most important components of our galaxy (Reynolds 1971; Monnet 1971). The intensity of the forbidden [S II] λ6716 and [N II] λ6584 lines is found to be higher with respect to the Hα intensity in the WIM in comparison with the typical intensities of the H II regions; in addition, the measurement of these lines reveals that the ionisation state of the WIM and its temperature vary significantly along the galactic plane, increasing towards the galactic centre (Reynolds 1985; Madsen & Reynolds 2005). Moreover, the filling fraction increases as the height with respect to the galactic plane increases, being around 0.1 at the midplane and reaching 0.3–0.4 at 1 kpc from the galactic plane (Kulkarni et al. 1987; Reynolds 1991; Berkhuijsen et al. 2006; Gaensler et al. 2008).

The first detection of this medium in other galaxies, now referred as diffuse ionised gas (DIG) in the extragalactic context, were carried out through narrowband Hα images of edge-on star-forming galaxies (Rand et al. 1990; Dettmar 1990). This was followed up by studies carried by small samples (< 5) of face-on and edge-on galaxies using narrowband Hα and [S II] images (Ferguson et al. 1996a; Greenawalt et al. 1997; Zurita et al. 2000) and long-slit spectroscopy with limited spatial coverage (Reynolds 1985; Wang et al. 1997; Otte et al. 2001; Hoopes & Walterbos 2003; Haffner et al. 2009).

All these studies have shown that 30%−50% of the total Hα luminosity of their galaxies correspond to the DIG, being the 80 ± 10% of the projected area of the disks occupied by DIG (Oey et al. 2007), thereby proving the ubiquitousness and omnipresence of this component in the star-forming galaxies. Moreover, the spectroscopic studies shown that the emission of low-ionisation lines, such as [S II]λλ6717, 31, [N II]λ6584, or [O I]λ6300, are enhanced in DIG and WIM regimes in comparison with the emission of the H II regions (Wang et al. 1997). Specifically, the ratios of those lines with respect the Hα emission in the DIG are larger than in H II regions (found firstly in edge-on galaxies; Otte et al. 2001). The leakage of the Lyman continuum (Lyc) photons, produced by the OB stars within H II regions, is commonly regarded as the primary source of DIG and WIM ionisation (Haffner et al. 2009; Seon 2009; Belfiore et al. 2022). However, it has been observed that in some cases, the total energy output of these stars does not appear to provide the necessary power to ionise the DIG uniformly across galactic discs (Ferguson et al. 1996a,b). The substantial power demand, combined with the presence of high-ionisation species like [O III] (enhanced at more than 1 kpc above the galactic plane) and He I (Reynolds et al. 1998; Wood & Mathis 2004; Haffner et al. 2009) underscores the insufficiency of Lyc photons leaking from H II regions in explaining the existence of the DIG. The question of the ionisation source for this component continues to be a topic of ongoing debate, with numerous proposals suggesting alternative heating sources for the DIG. Some of these sources include photoelectric heating from interstellar dust particles or large molecules (Reynolds & Cox 1992; Weingartner & Draine 2001), fast shocks provided by Wolf-Rayet stars or supernovae (Reynolds et al. 1998; Collins & Rand 2001; Hidalgo-Gámez 2005), turbulent mixing layers and dissipation of turbulence (Slavin et al. 1993; Minter & Spangler 1997; Reynolds et al. 1998; Binette et al. 2009), cosmic ray heating (Wiener et al. 2013), or microflares and magnetic field reconnections (Raymond 1992; Birk et al. 1998). Another alternative source of ionisation, as proposed by Flores-Fajardo et al. (2011) and Cid et al. (2011), is the photoionisation from hot low-mass evolved stars (HOLMES), that could provide an explanation for the high [O III]/Hα ratio observed in the extraplanar DIG (Rand et al. 1990).

The extent to which HOLMES contribute to the ionising radiation responsible for the DIG remains uncertain. While they may play a role in ionising specific regions within the sparsely populated interstellar medium, their overall significance is not yet established. Initial calculations conducted by Hills (1974) proposed that the ionising radiation emitted by these hot pre-white dwarf stars could have a substantial influence on the interstellar medium.

The exploration of the combination of Lyc photon leakage and HOLMES as potential primary sources of ionisation for the DIG has been extensively investigated using integral field spectrographs (IFS). The IFS offers simultaneous spatial and spectral information, enabling the examination of these ionisation mechanisms across a wide range of spatial resolutions and high spectral resolution. Some IFS surveys provides large samples of nearby galaxies at kiloparsec resolution such as the Calar Alto Legacy Integral Field spectroscopy Area survey (CALIFA; Sánchez et al. 2012; Husemann et al. 2013); 391 galaxies at ∼0.8 kpc resolution; Lacerda et al. 2018, hereafter Lac18) or the Mapping Nearby Galaxies at APO (MaNGA; Bundy et al. 2015; 356 galaxies at 2 kpc resolution; Zhang et al. 2017, hereafter Z17).

Z17 found an enhancement in the DIG of the ratios [S II]/Hα, [N II]/Hα, [O II]/Hα, and [O I]/Hα, increasing in function of the galactocentric radius, as well as a variable trend of the [O III]/Hβ ratio in comparison with H II regions. Besides, Lac18 defines a criterion of DIG definition in function of the Hα equivalent width (EW_Hα), based on if the region if dominated by photoionisation by HOLMES, corresponding to a EW_Hα < 3 Å. This would then be the primary regime in E-S0 galaxies, as the old stellar populations conform to these galaxies.

On the other hand, IFS studies based on the MUSE instrument (Multi Unit Spectroscopic Explorer; Bacon et al. 2010) or wide-field spectroscopic coverage such as TYPHOON (Grasha 2022) offer datasets with a spatial resolution ranging from ∼50 pc to the resolution of MaNGA-like galaxies, which allows for a spatially resolved study of the DIG and H II regions. Nevertheless, the studies of the DIG using resolutions of ∼50 pc have predominantly focused on the best galaxies within these surveys in terms of spatial resolution. These studies typically involve individual galaxies, such as M83 (Poetrodjojo et al. 2019; Della et al. 2022a,b), or the 19 galaxies encompassed by the Physics at High Angular resolution in Nearby GalaxieS project (PHANGS; Emsellem et al. 2022; Belfiore et al. 2022, hereafter Bel22), capable of resolving the ISM structure at ∼50 pc of resolution. Studies conducted within the PHANGS galaxies reveals that photoionisation by HOLMES is more prevalent in central regions where EW_Hα < 3 Å. In these regions, approximately 2% of the total Hα emission is powered by HOLMES. The remaining emission arises from photon leakage from H II regions. The mean free path of ionised photons is estimated to be < 1.9 kpc, based on a straightforward thin-slab model for photon propagation through the ISM (Zurita et al. 2002; Seon 2009). This model also predicts an increase in the [O III]/Hβ ratio and the ionisation parameter, which contradicts the decrease in the ratio as Σ_Hα increases. However, the contribution of the spectral hardness of HOLMES to the ionisation budget explains this discrepancy and the trend of increasing [S II]/Hα, [N II]/Hα, and [O I]/Hα line ratios as Σ_Hα decreases.

In order to make progress in quantifying the escape fraction of ionising photons from H II regions, it is imperative to reassess the importance of leaked radiation and HOLMES in contributing to the ionisation of the DIG in star-forming disk galaxies. However, previous studies aimed at this task have relied on large surveys with low spatial resolution or small samples with high spatial resolution. In most cases, these studies have employed methodologies that primarily rely on Hα emissions for the definition and exploration of the DIG.

In this paper, we introduce the Bidimensional Exploration of the warm-Temperature Ionised gaS (BETIS) project. The main goals of BETIS are the exploration of potential ionisation mechanisms, assessment of the DIG’s influence on calculating chemical abundances and star formation rates, and examination of the correlations between the DIG and various factors such as galaxy morphology, star formation rate, neutral hydrogen abundance, and other physical parameters. For this first paper of the BETIS series, we present the methodology employed in this project, designed to be applicable to galaxies with varying linear resolutions and morphologies, based not only on Hα, but on other emission lines. In addition, we discuss the problematical aspects of performing a general and global study of the ionisation mechanisms of the DIG. We validate this methodology through testing on a representative sample of seven galaxies with diverse characteristics. The goal is to determine the implications of incorporating galaxies exhibiting different physical processes into the same diagnostic of the ionisation mechanisms of the DIG.

This paper is structured as follows. Section 2 provides a definition of the showcase sample selection from different MUSE observations. Section 3 outlines the methodology for distinguishing between DIG emission and H II regions emission. We also introduce a novel adaptive binning method that will be applied to the observed datacubes, designed to enhance the signal-to-noise ratios (S/Ns) of our data. This methodology is tested on a showcase sample in Sect. 4, showing preliminary results, as well as a discussion of the challenges encountered when using EW_Hα as a proxy for characterising the DIG. Finally, Sect. 5 provides a overview of the key findings and outcomes presented in this paper. Additionally, we outline our plans for future research, specifically focusing on the forthcoming parts of the DIG analysis with a full BETIS sample.

The following notation is used throughout the paper: [N II] ≡ [N II]λ6584; [S II] ≡ [S II]λ6717+[S II]λ6731; [O III] ≡ [O III]λ5007; and [O I] ≡ [O I]λ6300. We adopt the standard ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_λ = 0.7, Ω_M = 0.3.

2. The BETIS showcase sample

The task of characterising and understanding the DIG in terms of star formation requires spatially resolved information, so that we may discern the DIG emission from the H II regions and explore the possible ionisation mechanisms. For this reason, the IFS surveys are the ideal data sets for studying the DIG, for instance, MaNGA (Z17) and CALIFA (Lac18), since they bring the necessary optical coverage to undertake this challenge. However, the limited spatial resolution has been the main obstacle in such studies. In terms of resolution and spectral coverage, MUSE is one of the most suitable instruments to study the DIG at extragalactic level.

MUSE is an IFS located at the ESO-VLT 8.1 m telescope at Cerro Paranal, Chile. The instrument brings 1 arcmin² of field of view (FoV), with a sampling of 0.2 × 0.2 arcsec, a spectral range of 4650 − 9300 Å, with a resolution of 1750 at 4650 Å to 3750 at 9300 Å, and a spectral sampling of 1.25 Å. The data observed by this instrument has already been used to study the DIG (PHANGS-MUSE; Bel22; Della et al. 2022a,b; Congiu et al. 2023) at a spatial resolution of ∼50 pc, but constrained to only < 20 objects.

In order to develop a methodology to study the DIG at different physical resolutions (from PHANGS-like to CALIFA-like) and present initial findings of this analysis, a representative subset of seven galaxies of different morphologies, characteristics and resolutions has been selected from different MUSE projects (which will be included in the full BETIS sample for the forthcoming papers). This showcase sample has been carefully chosen to provide insight into the broader study and serves as a demonstration of the methodology employed. In addition, this sample will also serve to prove the impact of using galaxies of different characteristics and physics involved in the same analysis of the DIG.

The showcase sample was selected from the All-weather MUse Supernova Integral-field of Nearby Galaxies (AMUSING; Galbany et al. 2016) and AMUSING+ (Galbany et al. in prep.) project samples. These surveys are ongoing projects aimed at studying the environments of supernovae by means of the analysis of a large number of nearby supernova host galaxies (0.005 < z < 0.1). The AMUSING survey¹ currently comprises 571 supernova hosts observed during the semesters P95-P104 (April 2015-March 2020) and is composed by a wide variety of galaxy types with the common characteristic of having hosted a known supernova. The AMUSING+ sample is aimed to increase the number of supernova host, adding 143 objects to the previous AMUSING sample. The criteria followed for selecting the showcase sample galaxies from the AMUSING and AMUSING+ samples are as follows.

First, to test the methodology for different linear resolutions, we chose galaxies with resolutions between the full width at half maximum for PHANGS and CALIFA (41 pc ≲ FWHM ≲ 1 kpc). Next, we aimed to analyse the radial behaviour of the DIG. However, since the physical properties of H II regions can vary between the inner and outer parts of galaxies (Rodríguez-Baras et al. 2019), to avoid limiting ours study to the inner regions, we assumed that at least 0.5R₂₅² of the galaxy must be in the datacube FoV if the resolution is high (FWHM ≲ 100 pc). Then, to test the methodology for different galaxy morphologies, we assumed every galaxy must have a different Hubble type, while also exhibiting star-forming regions. Finally, we did not consider galaxies with low depth or data-quality. In particular, if the average spectroscopic S/N of the Hα line in the DIG regions was lower than 3, the galaxy was disregarded.

Seven galaxies were chosen based on the previous criteria, belonging to the ESO periods P95, P98, P99, P101, and P103. The respective datacubes were already reduced by the AMUSING collaboration using the ESO reduction pipeline (Galbany et al. 2016; Weilbacher et al. 2020). Table 1 shows the general characteristics of the BETIS showcase sample and Table 2 shows the characteristics of the observations. The FWHM measured of the galaxies were derived assuming the seeing reported by the telescope via the TEL.IA.FWHMOBS parameter. We caution that this value corresponds to the PSF measured from the differential image motion monitor (DIMM) at Paranal and propagated to the data headers. This value may be unreliable, especially in cases of high wind, however, we can assume this value as the seeing since we are not working with the native resolution (see Sect. 3). The integrated stellar mass of the galaxies, obtained from a single stellar population (SSP) fitting with STARLIGHT (see Sect. 3), spans between 3.31 × 10¹⁰ M_⊙ < M_⋆ < 1.04 × 10¹² M_⊙, while the integrated star formation rate (SFR) spans between 1.31 M_⊙/yr < SFR < 4.92 M_⊙/yr. The SFR was obtained using the Kennicutt & Evans (2012) relation: log(SFR) = log(L_Hα) – log(C_Hα), with L_Hα as the Hα luminosity, corrected for interstellar extinction assuming the Cardelli extinction law, assuming R_V = 3.1 (Cardelli et al. 1989) and a Balmer decrement Hα/Hβ = 2.85. Then, log(C_Hα) is the conversion factor between SFR and L_Hα, corresponding to 41.27 (Kennicutt & Evans 2012). Figure 1 shows false-colour images of the sample constructed as a composition of the [S II] (green), Hα (red), and [O III] (blue) emission line maps obtained following the methodology described in Sect. 3.2, but on a spaxel-by-spaxel basis. The native spatial resolutions vary from 49 pc (IC3476, i.e. PHANGS-like) to 960 pc (ESO325-43, CALIFA-like), with a median resolution of 395 pc. The sample includes three known active galactic nuclei (AGNs): NGC863 and NGC3393 are Seyfert 2 (Weedman 1977; Lipovetsky et al. 1988), while NGC692 is a low-luminosity AGN (LLAGN; Zaw et al. 2019). We deliberately incorporated these galaxies presenting AGN emission. Although the AGN emission of the centre of the galaxy are masked during the methodology and analysis of the results, our aim is to investigate the impact of including galaxies exhibiting different physical processes via a global analysis of the ionisation mechanisms of the DIG in galaxies.

Fig. 1. RGB synthetic images of the galaxies selected for our sample. The images are constructed as a composition of the [S II]λ6716 (green), Hα (red), and [O III]λ5007 (blue) spaxel-by-spaxel emission maps obtained following the methodology explained in Sect. 3.2.

3. Methodology

The study of extragalactic DIG has conventionally relied on analyses of high-resolution narrowband Hα and images of nearby galaxies (z < 0.01, Rand et al. 1990; Dettmar 1990; Ferguson et al. 1996a; Zurita et al. 2000) or long-slit spectroscopy with limited spatial coverage (Reynolds 1985; Wang et al. 1997; Otte et al. 2001; Hoopes & Walterbos 2003; Haffner et al. 2009). In these studies, the methodology for subtracting the DIG from the galaxies and distinguish their emission from the H II regions emission is typically based on using cut-offs in the surface brightness (Σ_Hα) of the Hα line. In studies with images at a high resolution, a morphological definition of the H II regions using automatised tools is performed (Walterbos & Braun 1994; Zurita et al. 2000, 2002; Oey et al. 2007). Nevertheless, these studies have been constrained only to less than five face-on and edge-on galaxies in their sample, based only in Hα in those studies with high resolution and with low resolution in those with spectroscopic information. However, the use of the IFS in recent years has allowed for the adoption of methodologies for DIG subtraction and analysis that are based on spatially-resolved spectroscopic information with broad samples, (e.g. MaNGA, CALIFA, or PHANGS).

For example, DIG studies using data from the MaNGA survey (Z17; Jones et al. 2017) made use of 365 nearly face-on star-forming galaxies with 2.5 arcsec PSF FWHM, which is not enough to resolve individual H II regions which typical sizes are ∼100 pc. In this particular case, Z17 performed a cut out on Σ_Hα > 10³⁹ erg s⁻¹ kpc⁻² to select those spaxels dominated by H II regions emission. On the other hand, the studies based on CALIFA data (Lac18), made use of 391 galaxies with a median PSF FWHM of ∼0.8 kpc and proposed an alternative method to separate the DIG and the H II regions using the Hα equivalent width (EW_Hα), since (according to the authors) this parameter is a more appropriate proxy to distinguish the fundamental differences between the DIG and the star forming regions in comparison with the Σ_Hα. In addition, the authors argue that the usage of Σ_Hα could lead to misclassifications of low-surface-brightness H II regions such as DIG.

Other authors have used the criteria to define and classify the DIG regime based on EW_Hα continued to be used to define the DIG (Vale et al. 2019; Espinosa-Ponce et al. 2020). However, redefining the DIG a priori as the ionised gas consistent with ionisation from HOLMES could be problematic if the main goal is to discern the different ionisation mechanisms (Bel22). This is a challenging problem, as the standard resolution of CALIFA or MaNGA can not discern individual H II regions in order to characterise the DIG morphologically. In addition, the high spatial resolution of some galaxies observed by IFS (e.g. M83-TYPHOON/PrISM; Poetrodjojo et al. 2019, M83-MUSE; Della et al. 2022a, PHANGS-MUSE; Bel22, Congiu et al. 2023), enables the use of automated tools in order to detect and remove individual H II regions in IFU datacubes with resolutions of ∼50 pc, but constrained to a few objects.

In our work, we outline the approach for determining the DIG by utilising the spectroscopic and morphological data generated by the MUSE-IFS, using a broad dataset with high resolution and with spectroscopic information, not only Hα. A summary of the methodology, designed for galaxies of inclinations below 45°, is provided below and it is subsequently explained in further detail in the following sub-sections:

First, to increase the S/N of the weak, low-surface brightness emission lines involved in the DIG study, a modified version of the adaptive binning technique from Li et al. (2023) is performed to the S/N map of the [S II] line, obtaining then a segmentation map of the galaxy. The segmentation map is then applied to the datacube, to obtain a binned datacube, in which each bin corresponds to the integrated spectra as the sum of the spaxels contained in the bin.

Then, a spectral fitting is performed on the binned spectra in order to derive emission line maps of the most important species for the DIG study. We considered nine species: the hydrogen Hα and Hβ Balmer recombination lines, the He Iλ5876 recombination line, and the collisionally excited forbidden lines: [O III] λ5007, [O I] λ6300, and the doublets [N II] λ6548, 6584 and [S II] λ6717, 6731. The Hα and Hβ equivalent widths (EWs) are also calculated in this step.

Finally, the DIG is separated from the emission of the H II regions using a combination of an automated tool to detect and subtract the H II regions from the binned Hα maps and a cut-off in the Hα surface brightness (Σ_Hα) to subtract bright, irregular H II regions not detected by the automated tools. This results in a mask that corresponds to the lower-limit of the DIG. This mask is applied to all the binned emission line maps, creating a set of the lower-limit DIG emission. The upper-limit DIG emission is derived from the lower-limit binned maps considering a constant DIG emission column above the projected areas of the H II regions. We will generalise the methodology for inclinations above 45° and edge-on galaxies in future papers.

3.1. Adaptive binning

Typically, the distinction between SF regions and DIG has been made using Hα, but if we aim to explore the DIG in using all available spectroscopic information, we must consider key emission lines in the study of the DIG, such as [O I], [N II], [O III], and [S II]. For a reliable analysis of the DIG, is crucial to take into account S/N limitations of the data, specially in those lines of lower surface brightness, such as those mentioned above.

The signal of an emission line feature of a spectrum can be defined as the difference between the maximum flux value of an emission line f(λ_em) centred on λ_em of width w_em (in Å), and the mean of the fluxes in the two adjacent pseudo-continuum bands f(λ_c1) and f(λ_c1) of widths, w_c1, w_c2, measured on the spectrum of a given spaxel; and the noise corresponds to the mean of the flux within the two adjacent pseudo-continuum bands, f(λ_c1) and f(λ_c1) (Rosales-Ortega et al. 2012).

For the sake of clarity, in this work, the S/N of an emission line is the ratio of the amplitude of the emission line, defined as the peak of the line minus the mean flux of the pseudo-continuum, over the standard deviation of the adjacent pseudo-continuum on the spectrum:

$$\begin{array}{c}\hfill S/N({\lambda }_{\mathrm{em}})=\frac{\mu }{\sigma }=\frac{f({\lambda }_{\mathrm{em}})-⟨f({\lambda }_{c1}),f({\lambda }_{c2})⟩}{\sqrt{{\sigma }^2(f({\lambda }_{c1}),f({\lambda }_{c2}))}}.\end{array}$$(1)

The most commonly employed approach to enhance the S/N ratio of the data is to conduct adaptive binning, where individual pixels are combined into larger entities known as "bins" until the desired target S/N is achieved, However, this comes at the expense of diminished spatial resolution. The Voronoi binning method (Okabe et al. 2000; Cappellari & Copin 2003) generally resolves the issue of retaining the highest spatial resolution of the images while adhering to the minimum S/N limitation. This is achieved by tessellating the image and adjusting the bin size to ensure that every bin attains the desired S/N. However, this technique is particularly well suited for analysing elliptical and featureless galaxies as it relies on the continuum S/N. It may not be as effective for our purpose, as our maps show irregular structures, such as those depicting nebular gas emission lines. Consequently, we require an alternative method to avoid losing the primary morphology of spiral galaxies.

To solve this, we used a modified version of the adaptive binning method, introduced by Li et al. (2023) in order to enhance the S/N of weak emission lines. The code takes the flux of the line, the noise, and target S/N as its input to create a series of maps that cover the same area as the input data (map₁, ..., map_N, ..., map_Nmax), as seen in Fig. 2. For each map_N the set of N × N pixels are averaged, and, as a consequence, the S/N will be increased. If the S/N of the bin is not the target S/N, the code takes the map_N + 1 value instead. In our algorithm, we modified this code to get a segmentation map for binning the observed cube and getting a set of binned emission line maps (one per each considered line). The new binning technique was carried out as follows:

Fig. 2. Example of an adaptive binning of NGC863 using S/N([S II]) = 10 as target. The first figure represent the map numbers (mapN) from Li et al. (2023) algorithm. The second figure is the segmentation map obtained from our algorithm, where the bin index goes from k = 1 to Nbins. The third figure is a close-up of the segmentation map plotting the borders of the bins over the Hα map. It is noticeable that the H II regions, whose S/N is higher, are not binned, maintaining the structure pixel by pixel, and the bins are getting bigger as we move out of the H II regions.

First, we take as input the signal and noise defined in Eq. (1) and a target S/N. The code then creates the series of map_N; however, instead of recovering the new averaged flux, we save an index, k, in the coordinates [i₀, j₀], [i₀, j₁], [i₁, j₀], [i₁, j₁], ..., [i_N, j_N], whose flux can then be averaged. The next subset of coordinates whose flux will be averaged will have an index of k + 1. Then, we perform the previous step for all subsets of coordinates for k = 1, ..., N_bins. This results in a segmentation map, with the lowest bin indices k corresponding to the higher bins, namely, the 1 × 1 bins that correspond to the pixels with S/N higher than the S/N target and the lower indices corresponding to the bigger bins, namely, those needing more pixels to reach the target S/N (as seen in Fig. 2). Finally, for each bin of the segmentation map, we save the integrated spectra of the observed cube as the sum of the spaxels contained in the bin. This result in a binned observed datacube.

Binning the observed cube offers the advantage of enhancing the S/N of the spectra, rather than merely reducing the relative error associated to the Gaussian fitting of the lines, as is the case with direct binning the emission line maps. Moreover, if our goal includes generating binned EW maps, it is not methodologically correct to bin the EW map obtained through spaxel-by-spaxel fitting. This is because EW is not an additive quantity; it varies when there are changes in the underlying stellar continuum. Hence, it is important to calculate the new EW from the integrated spectra. If we want to increase the S/N of all the lines of interest, we need to choose a target feature as a basis to construct the segmentation map for the observed cube, ensuring to recover the maximum spatial information and signal in the DIG regimes. We employed S/N([S II]) as the target for adaptive binning because it is a low-ionisation collisional excitation line that remains unaffected by the correction for stellar population after the spectral fitting process.

Considering the most important low-ionisation species that are found to be enhanced in the DIG regimes in SF galaxies, such as [N II] and [S II] (Z17, Lac18), the [S II] line exhibits the lowest S/N (mean of 5, median of 3) among the low-ionisation species. Furthermore, the [O III] line exhibits a mean S/N of 2.5, therefore, if a ∼5σ detection is required, the S/N needs to be increased by a factor of 7 (Bel22). To achieve this, we performed our adaptive binning method using S/N([S II]) = 10 as the target S/N. Figure 3 (right panel) shows the distribution of S/N([S II]) obtained in the integrated spectra, where we can find that we get the average S/N of 10 that we expected.

Fig. 3. NGC863. Left panel: [S II] emission line map obtained from the Gaussian fitting of the [S II] line from the nebular gas cube spaxel-by-spaxel. The nebular cube was obtained as the difference between the observed cube and stellar cube as described in Appendix A. Central panel: binned [S II] emission line map, from the binned nebular gas cube, after performing an adaptive binning with target S/N([S II]) = 10 to the observed cube, as described. Both panels are close-up of the same region than Fig. 2. Right panel: distributions of S/N of the [S II] line for NGC863, calculated with the Eq. (1). Red distribution corresponds to the S/N([S II]) measured in the observed cube, namely, the S/N([S II]) of the spaxels. Green distribution corresponds to the S/N([S II]) measured in the binned cube, namely, the S/N([S II]) of the bins. The vertical dark green line marks the mean value of the S/N([[S II]]) of the bins, showing that the adaptive binning technique fulfils the goal of reaching a target S/N([S II]) = 10, on average.

3.2. Spectral fitting

For every new spectrum, we performed a SSP fitting with STARLIGHT (Cid Fernandes et al. 2005), with the obtained nebula emission spectra being the difference between the integrated observed spectra and the synthetic stellar spectra obtained after the fitting. We used the CB07 base spectra for the fittings (Bruzual & Charlot 2003; Bruzual 2007b). The number of SSPs (N_* = 100) from CB07 comprises three metallicities, Z = 0.2 Z_⊙, Z_⊙ and 2.5 Z_⊙, and 15 ages from t = 0.001 to t = 13 Gyr. All SSPs are normalised to 1 M_⊙ at t = 0. Their spectra were computed with Padova-2004 evolutionary tracks models, and Chabrier (2003) IMF (0.1 M_⊙ < M < 100 M_⊙). Then we computed a Gaussian fitting around the nine lines of interest on every nebular spectra, using the python library MPFIT, resulting in a set of nine binned emission line maps. We refer to Appendix A for further details about the SSP fitting and emission line fits. Figure 3 shows the difference between the [S II] emission line map obtained after performing the SSP fitting on the observed cube spaxel-by-spaxel (left panel) and on the binned cube applying the exposed method (central panel).

The Hα and Hβ EWs are also calculated during this step using the ratio between the integrated observed spectra and the fitted stellar model, resulting in a normalised stellar absorption spectra, where the EW is fitted for every bin using again MPFIT.

3.3. Morphological definition of the DIG limits

The separation between SF regions and the DIG has been always the first topic of discussion when studying and analysing the physics of the DIG. Historically, the most common method has been to separate the DIG from SF regions based on a Hα surface brightness (Σ_Hα) cut-off (Zurita et al. 2000; Oey et al. 2007; Z17). However, this method presents the problem of misclassifying low-surface-brightness H II regions as DIG, in addition to the possibility of classifying the emission of two overlapped DIG regions as H II regions (Lac18).

A suitable alternative for this problem is using automatised tools for the detection and subtraction of the H II regions individually. Several tools have been developed for this task, including SEXTRACTOR (Bertin & Arnouts 1996), HIIPHOT (Thilker et al. 2000), HIIEXPLORER (Sánchez et al. 2012), PYHIIEXPLORER (Espinosa-Ponce et al. 2020), ASTRODENDRO (Della et al. 2022a), PYHIIEXTRACTOR (Lugo-Aranda et al. 2022), and CLUMPFIND (Congiu et al. 2023).

The spatial resolution is an important key to detect and define morphologically individual H II regions. For the specific range of resolutions of the BETIS sample, we used PYHIIEXTRACTOR (Lugo-Aranda et al. 2022) in order to detect H II regions candidates. The code detects the H II regions candidates and assigns them a radius and a centroid, in the image coordinates. It is highly efficient for the detection of circular H II regions, especially those with low Σ_Hα, solving the aforementioned problem. Nevertheless, the limiting factor of this algorithms is in the complexity of the regions and in the variety of shapes and sizes of the brighter ones; hence, a circular extraction is not sufficient for those complex regions, since it would not be adapted to their morphology.

In the pursuit of the optimal methodology for segregating DIG emissions from H II regions, we conducted an extensive series of tests employing a variety of galaxies, techniques, and algorithms, considering the different observing conditions and depth of the BETIS sample. Our findings have led us to the conclusion that to establish a reliable morphological definition of the DIG, a combination of two methods is necessary: (i) the detection and masking of bona-fide H II regions exhibiting regular (circular) morphology based on the Hα emission line map; and (ii) the implementation of a cut-off in Hα surface brightness to account for irregular and highly luminous H II regions that may elude automated masking tools.

For the first step, we made use of the code PYHIIEXTRACTOR for a first detection of the H II regions. The algorithm takes a Hα map and detects the positions and radii of all the candidates H II regions. We used those positions and radii to perform a mask of H II regions, applied to the binned Hα map.

For the second step, we performed a Σ_Hα cut-off to the binned map. This cut-off is defined as three times the standard deviation of the surface brightness distribution of the masked map, as defined in the initial step (3σ_Σ(Hα)). This cut-out corresponds to three times the average Hα background level of the showcase sample. Performing this cut-off will result in a second DIG mask, that in combination with the first step DIG mask will give us the lower limit of DIG, as seen in Fig. 4, since we are assuming that all the emission coming from the H II positions is due to the star formation and the DIG contribution can be neglected.

Fig. 4. Steps to get the lower and upper DIG limit for NGC863 galaxy. The upper-left part corresponds to the DIG mask obtained by masking the H II candidates given by the centroids and radii obtained from PYHIIEXTRACTOR to the binned Hα map. The upper-right part corresponds to the binned Hα map after applying a 3σΣ(Hα) cut-off. The lower-left part is the mask of the lower limit of the DIG, obtained as the combination of the upper-left and upper-right masks. The lower-right part corresponds to the upper limit, assuming a constant non-zero flux level above the H II regions, interpolating flux around these regions.

We can also estimate an upper limit for the DIG, rejecting the previous assumption and considering a constant DIG emission column above the projected areas of the H II regions (Zurita et al. 2000; Congiu et al. 2023). This limit is estimated in two steps. Firstly, for every H II region of centroid (‘X’, ‘Y’) and radius ‘R’ detected by PYHIIEXTRACTOR, we defined an annulus centred in that detection, with inner radius equal to the radius obtained and outer radius 1.4 times the inner radius. We filled the H II region with the mean surface brightness measured inside the annulus. After performing this ‘filling’ with every H II region detection, part of the galaxy is still masked due to the Σ_Hα cut-off. To estimate the DIG contribution in these remaining regions, we filled them with the mean value of surface brightness of the border of the region. This result in a upper limit Hα DIG map, as seen in Fig. 4. The lower-limit DIG mask is then applied to the rest of the binned emission line and EW maps, which is the basis of our data.

4. BETIS: first results

After performing the binning methodology outlined above, we found an average bin size of the sample of 627 pc (see Table 3 for a overview of the increment of the S/Ns after binning). The average bin size in our sample is larger than the typical bin sizes in studies of narrowband Hα images (e.g. ∼250 pc, Ferguson et al. 1996a; ∼35 pc (native), Zurita et al. 2000; ∼475 pc, Oey et al. 2007). Moreover, the bin dimensions exceed those of the mean Voronoi bins of the PHANGS-MUSE sample (∼200 pc; Bel22). Nonetheless, they are smaller than the average PSF FWHM of the CALIFA galaxies (∼800 pc; Lac18).

To ensure the robustness of our analysis, we exclusively considered DIG bins with a relative error of 40% or lower and a S/N greater than 3 for all the lines involved in the subsequent analyses of this section. We also exclude the central part of those exhibiting an AGN: NGC863, NGC692 and NGC3393. Those constraints leave us with 80% of the DIG bins.

4.1. The DIG fraction

The visual inspection and variation of the DIG, both within individual galaxies and across them, provide valuable insights into its origins. In the binned Hα maps of our showcase sample galaxies, a diverse range of structures within the ionised gas regions becomes evident. The variations in Hα emissions in different directions yield crucial information about the evenness on the intensity of the DIG distribution and its spatial association with prominent H II regions. Likewise, investigating changes in the diffuse fraction and Hα surface distribution within a particular galaxy is of significant interest as it aids in pinpointing the source of the DIG.

Figure 5 shows the distribution of the DIG fraction, defined as the ratio of the DIG flux to the total Hα flux, f(Hα)_DIG/f(Hα)_total, for each galaxy of the sample, in both the radial and cumulative distributions. The radial distribution is obtained by measuring the DIG fraction from a series of de-projected rings centred at the nucleus of each galaxy with a width of 0.05R₂₅, while the cumulative distribution is obtained by performing aperture photometry of the DIG fraction in the same series of de-projected rings that are used for deriving the radial profiles, for both the lower and upper DIG limits. Moreover, using the definition of radial distribution, we present in Fig. 6 the Hα surface-brightness distributions for the lower DIG limit (green) after applying the mask (defined in Sect. 3.3) for the galaxy with only H II regions (red) and the total galaxy (blue). All the fluxes from this section and hereafter have been corrected for interstellar extinction assuming the Cardelli extinction law, assuming R_V = 3.1 (Cardelli et al. 1989) and a Balmer decrement Hα/Hβ = 2.85.

Fig. 5. Distribution of the DIG fraction (f(Hα)DIG/f(Hα)total) for each galaxy of the sample. Red area represents the cumulative distribution between the upper and lower DIG limits, with a step of 0.05R25, being the solid red lined the mean value. Blue area is the radial distribution of the DIG obtained integrating the flux between rings of width 0.05R25 and with a step of 0.05R25, being the solid blue lined the mean value. R25 is the isophote at the blue brightness of 25 mag arcsec−2. This parameter is obtained from Hyperleda database.

Overall, the radial distribution of the DIG fraction tends to increase, while the cumulative distributions remain constant, but it is more affected by the galaxy morphology than the cumulative distribution. The presence of a barred structure in NGC692, NGC6627, and ESO325-43 is evident through a decline in both distributions between 0.2–0.4R₂₅, which decreases up to 20%. This effect is also reflected on the radial distributions of the Hα surface brightness (Fig. 6), which show an increment of the Σ_Hα for the H II regions distribution, reaching up to 10³⁹ erg s⁻¹ kpc⁻² between 0.2–0.4R₂₅ for those galaxies. NGC3393 exhibits the highest DIG fraction for both lower and upper limits (0.69–0.87), presenting a surface brightness that surpasses the H II regions, throughout the entire galaxy, resulting in a maximum difference of around 1 dex. Furthermore, NGC863 also exhibits higher Σ_Hα in the DIG than in the H II regions, notably towards greater galactocentric radii.

Fig. 6. Radial distribution of ΣHα for each galaxy, obtained by measuring the surface brightness of a series of de-projected annulus centred at the nucleus of each galaxy with a width of 0.05R25. The distributions are performed for the entire galaxy (blue), the lower DIG limit (green) and H II regions (red). The last panel represent the mean distributions for the sample.

The Sc galaxy ESO584-7 and the dwarf IBm galaxy IC3476 present the lower DIG fraction, 0.37–0.66 and 0.4–0.63, respectively. Additionally, they have the highest surface brightness for the H II regions in the sample, being > 10³⁹ erg s⁻¹ kpc⁻² in the inner parts of the galaxy and decreasing to ∼5 × 10³⁸ at 0.5R₂₅. The average radial distribution of the DIG fraction is not showing any tendency in particular, remaining constant probably due to the mixture of morphologies of our showcase sample. A full analysis considering a general distribution by morphological type will be carried out as part of BETIS in a forthcoming paper.

Table 4 summarises the results obtained from the DIG fraction for the showcase sample. The general findings for this showcase sample show that the DIG fraction ranges from 0.4 to 0.7, which coincides with the results of previous research, in terms of both the lower limit (e.g. Ferguson et al. 1996a,b; Zurita et al. 2000; Thilker et al. 2002; Oey et al. 2007; Bel22), and upper limit (e.g. Zurita et al. 2000; Congiu et al. 2023). In addition, the average Σ_Hα cut-out in also consistent with the average upper DIG level of previous studies (Poetrodjojo et al. 2019, Bel22, Congiu et al. 2023). In addition, the tendency of the cumulative distribution is to be constant and indicating approximately a 60% of DIG in these galaxies. The surface brightness of the DIG vary from ∼6 × 10³⁸ in the inner parts of the galaxies and decreasing monotonically to ∼5 × 10³⁷ erg s⁻¹ kpc⁻² in the outer regions, with a notably high integrated Hα luminosity between 10⁴⁰ and 5 × 10⁴¹ erg s⁻¹. This is consistent with previous studies (e.g. Ferguson et al. 1996a,b; Zurita et al. 2000).

The similarity in the Σ_Hα radial distributions between the H II and DIG regimes and its impact in the DIG fraction supports the correlation between H II regions emission and DIG photoionisation (Ferguson et al. 1996a,b; Zurita et al. 2000). However, performing this exploration individually shows that there are instances where DIG Σ_Hα values exceed those of the H II regions, in particular, in the two Seyfert galaxies, NGC863 and NGC3393. For instance, for NGC863, a total DIG Σ_Hα of 6.6 × 10³⁸ erg s⁻¹ kpc⁻² will require a power per unit area of 1.5 × 10⁻³ erg s⁻¹ cm⁻² to keep the DIG ionised, while the total Σ_Hα of the H II regions of 5.02 × 10³⁸ erg s⁻¹ kpc⁻² provides a power of 1.1 × 10⁻³ erg s⁻¹ cm⁻², which is insufficient to ionise the entire DIG³.

Hence, it is imperative to contemplate alternative ionisation sources that can provide additional energy supply to the ISM apart from Lyman continuum photons escaping from H II regions to keep the DIG ionised for those particular cases. The incorporation of collisional, low-excitation lines such as [N II], [S II], or [O I], as well as high-excitation lines such as [O III] and He I in the analysis is essential for comprehensively investigating the diverse array of ionisation mechanisms that exist within the (ISM).

4.2. BPT diagnosis of the DIG

We can explore the location of the DIG bins in the classical [N II] BPT diagram, as shown in Fig. 7. In the same figure we have included a diagnostic from Kopsacheili et al. (2020) for the separation of shock excited from photoionised regions. The BPT of the global sample shows the DIG falls mostly below the Kewley et al. 2001 demarcation, showing a photoionisation feature due to H II regions, but with high-excitation regions above the demarcation corresponding to AGN-like emission, as noted in previous studies and usually explained as photoionisation due to HOLMES (Lac18, Bel22). However, when performing this diagnosis for individual galaxies, only NGC863 and NGC3393, both Seyfert-2 AGNs, exhibit line ratios characteristic of AGN emission. All the DIG for the rest of the galaxies show photoionisation feature due to H II regions. The central region of NGC863 and NGC3393 is not considered in this analysis, so the AGN outflows could be the source of the gas ionisation with high-excitation lines found in the BPTs. We compare our results with the theoretical models of gas that is highly excited via fast-shocks of Allen et al. (2008), assuming solar metallicity and a pre-shock density of 1 cm⁻³. We plotted the predicted line ratios including (red lines) or not (blue lines) a photoionisation by precursor, and given a range of shock velocities (v_s) and magnetic field intensities (B). The line ratios observed in NGC863 are consisted with fast-shock without precursor with v_s between 200 and 500 km s⁻¹ and B between 0.0001 and 5. In the case of NGC3393, the ratios observed corresponds to fast-shock without a precursor, with v_s between 200 and 1000 km s⁻¹ and B between 0.0001 and 10. Furthermore, the fact that the ionisation bi-cone of NGC3393 and its continuum emission are uncoupled (Maksym et al. 2016), along with the emission line ratios in the BPT diagram indicative of ionisation by fast shocks, suggests that we are tracing gas outflows rather than DIG emission (López-Cobá et al. 2020).

Fig. 7. BPT diagnosis of the sample for the DIG bins that verifies that the relative errors of [O III] line bins are below 40% and its S/N above 3. The last panel displays the seven prior panels together. The central parts of NGC863, NGC3393, and NGC692 were excluded in all panels due to their strong AGN emissions. Each contour encloses the 10% of the points, with every point a single bin. Black lines are given by Kewley et al. (2001) (solid) for the classic H II regions photoionisation and AGN demarcations. The dashed-dot line depicts one of the 2D diagnostics developed by Kopsacheili et al. (2020) for the separation of shock excited (e.g. supernova remnants) from photoionised regions (e.g. H II regions). Coloured lines represent the models of photoionisation by fast shocks from Allen et al. (2008). The blue and red models illustrate photoionisation where only front shocks occur and when pre-ionisation by a precursor is taken into account. The solid model curves plotted represents shocks winds of 200, 400, 500 and 1000 km s−1, and dashed model curves represents magnetic field intensities of 0.0001, 1.0, 5.0 and 10. Red and blue arrows represent the direction of increasing wind velocity in each model.

The presence of AGN outflows affects the overall BPT diagnosis, revealing that a section of the DIG is ionised by sources separate from the photon leakage originating from H II regions, but what we are introducing is the ionisation cone of the AGNs, mimicking the DIG emission. Therefore, when conducting a comprehensive diagnostic assessment of the DIG across a global sample, it is crucial to completely exclude galaxies that exhibit AGN emission. This occurs even if the whole DIG is ionised by star-forming regions, as observed in the remaining galaxies in the sample. This effect also explains the high Σ_Hα in the DIG found in the same two galaxies; NGC863 and NGC3393, as we are incorporating the AGN emissions.

4.3. DIG line ratios

Further evidence of the connection between the DIG and H II regions lies in the behaviour of the [S II]/Hα, [N II]/Hα, [O I]/Hα, and [O III]/Hβ ratios. Historically, studying the radial dependence of the DIG has been important, as its trend along different galactocentric distances can constrain its origin (Ferguson et al. 1996a). Furthermore, recent investigations have highlighted the importance of considering possible contamination and/or contribution of the DIG to physical parameters derived from emission lines, such as metallicity and star-formation (Bel22, Lugo-Aranda et al. 2024). Therefore, aperture (or, correspondingly, redshift) biases should be taken into account, for which the study of the DIG in a radial basis is a key element. An observational fact is that the low-ionisation [N II]/Hα and [S II]/Hα line ratios increase with decreasing Σ_Hα, the best example being the increase in these line ratios with increasing distance from the mid-plane, both in the Milky Way (Haffner et al. 1999) and other galaxies (Rand 1999; Tüllmann et al. 2000). Although the values for [N II]/Hα and [S II]/Hα could vary considerably in the DIG, they are correlated, often with a nearly constant [S II]/[N II] ratio over large regions (Haffner et al. 2009).

Figure 8, shows the average radial distribution of the [S II]/Hα, [N II]/Hα, [O III]/Hβ, and [O I]/Hα line ratios for the lower DIG limit. The distributions have been generated using the same method as for Σ_Hα (as shown in Fig. 6) and f_DIG (blue curves of the Fig. 5) radial distributions. Each data point is computed as the average value from the seven galaxies in the sample at a specific radius. We created these distributions for both the DIG (in green) and the area corresponding to H II regions (red). In all instances, we exclusively consider the bins with a relative error below 40% for each line. In all cases, the DIG line ratios show higher values than the line ratios corresponding to the H II regions.

Fig. 8. Radial distribution of the [S II]/Hα ratio of the DIG bins (green), and H II regions according to our definition (red), obtained as the ratio of the mean values of the respective fluxes within a ring of width 0.05R25 with a step of 0.05R25 from 0.2R25 to R25. We are only considering those bins whose relative error in [S II] and Hα lines are lower than 40%. To avoid the strong AGN emissions, we start all the distributions from 0.2R25. In general we see a tendency to increase the ratio as we move to outer regions of the galaxy. Besides, the ratio is in all cases higher in DIG regions, as we expected from the literature (Haffner et al. 2009).

The [N II]/Hα radial distributions (upper-left panel), for both the DIG and the H II regions, decrease with galactocentric distance. The [N II]/Hα DIG distribution ranges from 0.67 to 0.51, while for the H II regions from 0.44 to 0.37 between 0.2R₂₅ and 1R₂₅, being the [N II]/Hα line ratio of the DIG 0.15 dex higher in average. The [N II]/Hα line ratio is a well-know metallicity indicator (Storchi-Bergmann et al. 1994; Denicoló et al. 2002; Ho et al. 2015); thus, the decreasing value as a function of radius is primarily reflecting a change in metallicity of the ionised gas, in this case of the metallicity gradient of the spiral galaxies. Nevertheless, given that N⁺/N and H⁺/H vary little within the DIG, for a given metallicity, variations in the [N II]/Hα line ratio essentially trace variations in T_e (Haffner et al. 2009). However, calculating absolute temperatures is uncertain because of the requisite assumptions about the precise ionic fractions and elemental abundances.

Figure 8 suggest a flattening in the [N II]/Hα index in both the H II regions and the DIG distributions for r < 0.6R₂₅, indicating a potential contribution from the DIG to the observed radial metallicity gradients within galactic planes (Z17).

On the other hand, the [S II]/Hα line ratio increases radially in both H II regions and DIG regimes (lower-left panel of Fig. 8), with higher values observed in the DIG; from 0.28 to 0.63 for the DIG and from 0.19 to 0.12 for the H II regions between 0.2R₂₅ and 1R₂₅, with the [S II]/Hα DIG line ratio higher by 0.15 dex, on average. The shapes of both distributions display comparable patterns, reaching a maximum at 0.75–0.8R₂₅.

The [O I] line is produced by collisions of neutral oxygen with thermal electrons, while its intensity is a measure of the neutral hydrogen content within the DIG. The first ionisation potential of the oxygen is close to that of the hydrogen, and the large H⁺+O⁰↔ H⁰+O⁺ charge-exchange cross-section keeps O⁺/O nearly equal to H⁺/H. Thus, the [O I]/Hα ratio is related to the amount of H⁰ relative to H⁺, and it is a sensitive probe of the ionisation state of the emitting gas (Reynolds et al. 1998; Haffner et al. 2009). The volume photon emissivity e of the O¹D λ6300 transition relative to Hα is related to the hydrogen ionisation ratio n(H⁺)/n(H⁰), with a linear dependence of the gas-phase abundance of oxygen n(O)/n(H) and a weak dependence that tracks changes in T_e (Haffner et al. 2009). Observations of [O I] in the Milky Way and objects at z ∼ 0 are difficult because of the [O I] λ6300 air glow line, which is of order 100 times brighter than the interstellar line. Voges & Walterbos (2006) made the first DIG detection of [O I] λ6300 in any non-edge-on spiral other than the Milky Way near the H II region NGC 604 in M33, with observed [O I]/Hα ratios in the range 0.038 − 0.097. In the upper-right panel of Fig. 8, we observe the azimuthally averaged values of the [O I]/Hα ratio with a significantly steeper slope in the DIG compared to the emission in the H II regions, increasing from 0.12 to 0.79 between 0.25R₂₅ to 0.75R₂₅ (−0.9 to −0.1 in log), while the H II distribution increases from 0.06 to 0.18 (−1.2 to −0.7 in log), which are higher than the values reported by Voges & Walterbos (2006) for M33, and higher than the values found by Bel22 with an upper limit ∼0.3. The high [O I]/Hα line ratio in the DIG found in the BETIS showcase sample is challenging, and will be matter of study in forthcoming studies.

The radial distribution of the [O III]/Hβ line ratio is shown in the lower-right panel of Fig. 8. Generally, the variation in this ratio between DIG and H II regions depends on the specific physical characteristics of the ISM (Z17). Normal spiral galaxies show an increasing value of the [O III]/Hβ line ratio with increasing radius, mainly due to secondary dependence on metallicity, with a high dispersion in the central regions (∼1 dex for r < 0.4R₂₅). Previous studies also show a [O III]/Hβ ratio both higher in H II regions (Greenawalt et al. 1997; Galarza et al. 1999) and lower than the DIG regions (Collins & Rand 2001; Otte et al. 2002). In our sample, this ratio remains elevated in DIG regions, but the difference between the DIG and H II distributions is less pronounced compared to the other line ratio distributions. This trend suggests a radial decrease in this ratio; from 1.58 to 1.0 in the DIG (0.2 to 0.0 in log) and from 1.07 to 0.79 in the H II regions (0.05 to −0.1 in log). Our results indicate that these line ratios are typically higher in the DIG compared with the H II regions. This trend aligns with what is commonly reported in the existing literature (Haffner et al. 2009).

The similarity in the trends observed in the radial distribution of Σ_Hα and line ratios between both DIG and H II regions, especially in the case of [N II]/Hα and [S II]/Hα ratios, along with the majority of DIG bins falling within the photoionisation regime on the BPT diagram, suggests that the explanation for DIG behaviour can be attributed to photoionisation processes within the galactic plane, without the need for alternative sources of ionisation.

However, another phenomenon needs to be explained, the low values of EW_Hα found in the DIG, specially in regions with photoionisation regimes not corresponding to SF regions (Lac18, Bel22).

4.4. EW_Hα in the DIG

The EW_Hα has been used by many authors to differentiate between ionisation caused by star formation and AGNs and ionisation caused by a smooth background of hot evolved stars. Cid et al. (2011) used SDSS data to demonstrate that the emission-line galaxy population exhibits a bimodal distribution in EW_Hα, and that 3 Å serves as an empirical demarcation between these two. Later, Belfiore et al. (2016) using MaNGA data, showed the presence of extended (kpc scale) low-ionisation emission-line regions (LIERs) in both star-forming and quiescent galaxies, associated with low EW_Hα (< 3 Å). In SF galaxies, the LIER emission was associated with diffuse ionised gas, most evident as extraplanar emission in edge-on systems.

Lac18 proposed a separation of DIG ionisation regimes based on the EW_Hα and applied over all types of galaxies, including elliptical and S0. The regions where EW_Hα > 14 Å traces SF regime, 3 < EW_Hα < 14 Å reflects a mixed regime, and regions where EW_Hα < 3 Å define the component of the DIG where photoionisation is dominated by hot, low-mass, evolved stars (HOLMES; Flores-Fajardo et al. 2011; Cid et al. 2011). Those stars were proposed as an additional ionisation source of the DIG, in order to explain the high [O III]/Hβ ratio found in the extraplanar DIG in edge-on galaxies (Reynolds et al. 1998; Rand 1999; Flores-Fajardo et al. 2011).

The significance of employing the EW_Hα to distinguish between various ionisation sources, while considering the contribution of HOLMES to the energy budget, has been a topic of frequent discussion among several authors, even in the context of face-on galaxies (e.g. Cid et al. 2011; Lac18,Bel22); however it is not yet clear whether it reveals a true division (if any) between ionisation carried by HOLMES, star formation, or shocks. Moreover, Bel22 also noticed an unclear division for ionisation by HOLMES at EW_Hα = 3 Å and by star formation at EW_Hα = 14 Å suggested by Lac18, since their EW_Hα for DIG and H II regions tends to overlap in lower EW_Hα regimes.

However, the general methodology employed in IFS studies in order to measure both the line intensities and the EW of Balmer recombination lines poses a problem for spectra with low S/N values for the emission lines. This methodology calls for measuring the line intensities and EWs from a nebular spectrum obtained by subtracting a fitted stellar model to the observed spectrum. In the case of the Balmer recombination lines, this procedure typically makes a small correction to the emission line intensity due to the underlying stellar absorption. For EWHα > 3 Å, this correction is usually of a few percent. Nevertheless, when the S/N of the Hα line is ∼1 (Hα emission embedded in the stellar continuum), the resulting Hα emission, and therefore EW_Hα, would depend on the Balmer stellar absorption feature due to the SSP fitting. In this case, the underlying stellar absorption would not represent just a correction of the emission line, but the resulting flux would be fully originated from the fitted stellar spectrum, generating low EW_Hα values ≲3 Å.

This is shown in Fig. 9, where we plot the observed (green), fitted stellar model (black) and nebular (blue) spectra of regions with EW_Hα of 0.6, 1, and 3 Å. In the left and central panels, the Hα emission in the observed spectra is totally embedded by the stellar continuum. However, due to the fitted stellar model, the nebular spectra includes a Hα emission line of the same magnitude as the stellar absorption feature of the model, from which an EW_Hα of 0.6 Å (left) and EW_Hα = 1 Å (right) are measured. On the other hand, the right panel of Fig. 9 shows that for a region with EW_Hα = 3 Å, the underlying stellar absorption of the model makes just a small correction to the total flux of the emission line.

Fig. 9. Three examples of observed (green), model (black) and nebular (blue) spectra of different EWHα = 0.6, 1 and 3 Å. All observed spectra are integrated spectra from NGC863 bins, the model spectra are obtained from STARLIGHT following our methodology and the nebular spectra are obtained subtracting the model to the observed (see text). Gold lines represent the [N II] doublet.

Consequently, when the S/N of the Hα is ∼1 and/or when the emission is embedded in the stellar continuum (EW_Hα ≲ 3 Å), the resulting EW_Hα is dependent on the selection of SSPs and the fitting methodology. This sharp differentiation between HOLMES dominated regions, as well as the differentiation between LIERs and passive galaxies (Cid et al. 2011; Belfiore et al. 2016) may not be reliable, as it could potentially originate from a methodological artefact. We show a test of this effect for an integrated region of NGC863 in Appendix B.

Further evidence that the low EW_Hα regime (≲3 Å) imposes a methodological challenge is manifested when verifying the validity of physical parameters derived from the embedded emission lines, such as the Hα/Hβ Balmer decrement. Figure 10 shows the distribution of the Hα/Hβ ratio for different EW_Hα regimes, while for EW_Hα > 3 Å the distributions peaks at the theoretical Hα/Hβ ratio ≳2.87 (Osterbrock & Ferland 2006), for the regions with EW_Hα < 3 Å, the ratio shows non-physical values.

Fig. 10. Distribution of Hα/Hβ ratio for three different EWHα regimes. Each distribution represent the Hα/Hβ flux ratio for all bins between 6 Å < EWHα ≤ 10 Å (green), 3 Å < EWHα ≤6 Å (blue), and EWHα ≤3 Å (red), for the seven binned galaxies listed in Table 1. Black vertical dashed line presents the theoretical ratio of Hα/Hβ = 2.87 (Osterbrock & Ferland 2006). We can see that a substantial part of the bins with EWHα ≤ 3 Å shows a non-physical ratio (< 2.87).

In addition, the values of EW_Hα < 14 Å and < 3 Å are not exclusive to the DIG emission, so the demarcation with H II regions using the EW_Hα is not clear (Bel22). This is also evident in Fig. 11, which shows the distribution of EW_Hα for both DIG and H II bins within our sample. We found that the galaxies with lower median EW_Hα in the DIG regions (∼9 Å) are NGC863 and NGC3393, both Sa type, followed by NGC692 (Sbc), NGC6627 (Sb) and ESO325-43 (Sc), with EW_Hα ∼ 18 − 20 Å. ESO584-7 (Sc) and IC3476 (Im) have exceptionally high EW_Hα in the DIG regions, with 34 and 36 Å respectively. ESO584-7 is a H II galaxy (Contini et al. 1998), thus, the high luminosity and star formation rate may be the cause of this high median EW_Hα. The case of IC3476 is also special. It is a galaxy suffering the effects of the ram pressure stripping due to the Virgo Cluster intergalactic environment (Boselli et al. 2021). These authors show that the effects of this perturbations reach scales of individual H II regions r_eq ∼ 50 pc. Furthermore, the compression of the gas along its stellar disc may be the cause of the increase of the star formation activity. This increment could explain the exceptionally high EW_Hα in the DIG regions found in this galaxy, as happened in the case of ESO584-7.

Fig. 11. EWHα distributions for all the bins of the sample. Green histograms represents the distribution for the DIG emission bins, and red for the H II regions emission bins. Last histogram represents the distribution for the seven galaxies.

In general, this results are in concordance with previous authors (Lac18, Bel22), being the early type galaxies those with lower EW_Hα, due to the older stellar populations of their bulges, and late-types those with higher EW_Hα. The median EW_Hα of the DIG regions for all sample is ∼25 Å, substantially higher in comparison with previous studies (∼5 Å) due to the bias given by ESO584-7 and IC3476, and due to the small sample selected.

4.5. Dust reddening in the DIG regime

In Fig. 12, we show the reddening in the DIG bins, obtained as the radial distribution of the Hα/Hβ ratio for the seven galaxies of our sample. This ratio is always lower in the DIG regime (between 3.85 and 2.5), following the same tendency to decrease radially as in the H II regions, with ratio between 4 and 3.25. The decline of the Hα/Hβ ratio in DIG regimes shown are in concordance to the expectations. Since the Hα/Hβ ratio reflects the attenuation of young stars by dust both in H II regions and in the ISM (Chevallard et al. 2013), this ratio is expected to be lower in the DIG, due to the lower optical depth and the increased scatter in the dust attenuation-line luminosity relation (Vale et al. 2020). The fact that we find a Hα/Hβ ratio lower in the DIG than in the H II regions suggests that the extinction is consistently lower in the DIG. There could also be an effect of the temperature, as it has been found in the literature that the DIG temperature is ∼2000 K warmer than in H II regions (Madsen & Reynolds 2005).

Fig. 12. Mean radial Hα/Hβ ratio distribution of the sample. The distributions are obtained as in Fig. 8. Red represent the distribution for the H II regions bins, and green for the DIG bins. Both distributions decrease radially, on average, but the extinction is always lower in the DIG regime.

5. Summary and conclusions

In this work, we present the Bidimensional Exploration of the warm-Temperature Ionised gaS (BETIS) project, designed for the spatially resolved and spectral study of the diffuse ionised gas (DIG) in a selection of nearby galaxies. We present a methodology for characterising and studying, both spatially and spectroscopically, the DIG optimised for galaxies with different linear resolutions and physical characteristics. To validate our methodology, we selected a showcase sample consisting of seven galaxies with diverse morphological and characteristic traits. This methodology involves the following steps:

• An adaptive binning was performed to the observed datacube in order to increase the S/N of the fainter lines such as [O III], [O I], and [S II]. Our technique is based on the spectroscopic S/N of the [S II] line, with a target S/N of 10.

• We conducted a SSP synthesis using the STARLIGHT code for each integrated spectrum within the binned datacube. Subsequently, we employed Gaussian fitting to the residuals of each SSP fitting to derive the binned emission line maps for the nine lines of interest. Additionally, we generated a binned EW_Hα map.

• The DIG was separated from the H II regions using a combination of an automated tool to detect and subtract the H II regions from the binned Hα maps, together with a cut-off in Hα surface brightness to subtract bright, irregular H II regions not detected by the automated tools.

We found an average DIG faction of 40%–70% in the showcase sample, with NGC3393 the one with the higher DIG fraction (69%–87%), followed by NGC863, ESO325-43, NGC692, NGC6627, ESO584-7, and IC3476. Those with higher DIG fractions are the two Seyferts of the sample: NGC3393 and NGC863. This is further exemplified when analysing the radial distributions of the Σ_Hα in the DIG and H II regions. In these two galaxies, the DIG exhibits higher surface brightness compared to the H II regions, with the disparity between these two regimes reaching up to 1 dex. The overall radial distributions of Σ_Hα, as depicted in Fig. 6, generally reveal similar trends for both H II regions and DIG, with a radial decrease. However, there is an increase in Σ_Hα within the bars of NGC692 and NGC6627 for the H II regions.

On average, we observe in Fig. 8 higher [S II]/Hα, [N II]/Hα, [O III]/Hβ, and [O I]/Hα ratios in the DIG compared to the H II regions. Additionally, the radial trends of the DIG and H II distributions are similar in all cases, indicating a correlation between the ionisation of these species in both the DIG and the H II regions.

Computing the [N II] BPT diagram also highlights a significant distinction between the two Seyfert galaxies and the rest of the sample. It is evident that the DIG is predominantly photoionised by H II regions in all galaxies, except for NGC3393 and NGC863. In these two cases, the ionisation source of the DIG appears to be accounted for by the fast shock models proposed by Allen et al. (2008). Nevertheless, it is worth noting that these two galaxies host prominent AGNs, which can mimic the emission of the DIG when assessing a global BPT for the entire sample. In particular, NGC3393 presents a strong ionisation cone due to galactic outflows.

We also addressed the challenge of employing the EW_Hα as a proxy for delineating DIG regions and different ionisation regimes. This issue arises because at low EW_Hα, typically used to identify HOLMES or AGNs regimes, the Hα line derived from synthetic spectra after conducting a SSP fitting can be an artefact of the model. This results in an artificial Hα emission line when correcting the observed Hα emission with the stellar model absorption, after subtracting the model from an observed Hα line with a spectroscopic S/N ∼ 1. Therefore, the low EW_Hα value may be a result of the characteristics of the stellar models, making it dependent on these models – if we do not consider the spectroscopic S/N of the observed Hα.

These results suggest that conducting a global analysis of the DIG using a sample of galaxies with diverse characteristics may lead to misleading conclusions about the ionisation mechanisms. This is because each galaxy can present distinct physical processes, for example large-scale AGNs, which can mimic the emission of high-excitation DIG, making the large-scale AGN emission and high-excitation DIG indistinguishable in the same diagnosis. For this reason, and due to the lack of reliability of the low EW_Hα regimes, every galaxy needs to be considered individually when performing a DIG diagnosis.

The distributions of EW_Hα for the DIG regions in our sample exhibit a morphological pattern. Sa-type galaxies, with prominent bulges of older stellar populations, have lower EW_Hα, followed by Sb and Sc galaxies. Notably, ESO584-7 and IC3476 show significantly higher EWs. ESO584-7 is an H II galaxy, while IC3476 experiences elevated EWs due to ram pressure stripping in the Virgo Cluster’s intergalactic environment.

Lastly, we examined the impact of dust reddening in the DIG by assessing the radial distribution of Hα/Hβ ratios within both the DIG and H II regions across our entire sample. Our results suggest that extinction is consistently lower in DIG regions.

In the forthcoming papers of this series, we will apply the methodology and expand the analysis to a full BETIS sample selected from the AMUSING, AMUSING+, and AMUSING++ project samples (López-Cobá et al. 2020). We will explore how the results and methodology may vary based on galaxy morphology and inclination (de-projected face-on galaxies i < 45° vs edge-on galaxies i > 75°), we will investigate the influence of the DIG on the determination of various parameters, including chemical abundances and star formation rates. Additionally, we will explore other spectroscopic lines of interest for the DIG study, such the high-excitation He I, aiming to uncover potential new ionisation mechanisms of the DIG. Furthermore, we will explore the extragalactic DIG for edge-on galaxies, shedding light on the high-excitation regimes found at more than 1 kpc above the galactic plane, reported by previous authors.

Acknowledgments

Based on the observations at ESO with program IDs: 095.D-0172, 098.B-0551, 099.B-0294, 099.D-0022, 0101.D-0748 and 0103.D-0440. R.G.D. acknowledges the CONAHCyT scholarship No. 1088965 and INAOE for the PhD program. The authors also acknowledge Laboratorio Nacional de Supercómputo del Sureste de México (LNS) for allowing the usage of their cluster with the project No. 202201027C in collaboration with INAOE, and Itziar Aretxaga for allowing us the usage of the Toltec/GTM cluster. L.G. acknowledges financial support from the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de Investigación (AEI) 10.13039/501100011033, and the European Social Fund (ESF) “Investing in your future” under the 2019 Ramón y Cajal program RYC2019-027683-I and the PID2020-115253GA-I00 HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC) under the PIE project 20215AT016, and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. The images of Fig. 1 were created with the help of the NOIRLab/IPAC/ESA/STScI/CfA FITS Liberator and the free web-based image editor https://photopea.com.

References

Appendix A: Emission line measurement

A.1. Cube preprocessing

Each individual spectrum of the cube, after being brought to the rest-frame, must be corrected for Milky Way extinction. The correction is carried out by multiplying each spectrum by a factor 10^{0.4 ⋅ a_λ}, where a_λ is the extinction function of Fitzpatrick (1999), using the Python library extinction.fitzpatrick99⁴, that reads the R_V = A_V/E(B − V) ratio (fixed at 3.1) and the total V-band extinction in magnitudes for each galaxy, obtained from Hyperleda. Values can be found in the table of BETIS characteristics online.

In this work we use CB07 (Bruzual 2007a,b) base spectra for the SSP models. The model spectra have a resolution of R ∼ 2000 and 1λ of spectral sampling. Our observed spectra have R ∼ 3000 and 1.25λ of spectral sampling, so the observed spectra must be resampled at 1λ. Performing a linear interpolation of every single spectra of the cubes resolves the conflict between the samplings. After these steps, we obtain a new observed cube, rest-framed, resampled to Δλ = 1 Å and corrected to Milky Way extinction. This new resampled cubes are those using to perform the SSP fittings.

A.2. SSP synthesis

Performing a SSP synthesis (Tinsley 1968) requires an estimation of the type of stellar populations, namely, the masses, ages, and metallicities found in a galaxy, star cluster, or region of a galaxy based upon its spectra. The SSP synthesis is carried by the STARLIGHT software (Cid Fernandes et al. 2005), which fits an observed spectrum, O_λm with a model, M_λ, which adds up N* spectral components from a pre-defined set of base spectra. The synthetic model spectra that the program generates take the form:

$$\begin{array}{c}\hfill M_{\lambda }=M_{{\lambda }_0}\left(\sum _{j=1}^{N_{\ast }}{x}_j{b}_{j,\lambda }{r}_{\lambda }\right)\otimes G(v_{\ast },{\sigma }_{\ast }),\end{array}$$(A.1)

where b_j, λ is the spectrum of the jth SSP normalised; r_λ is the reddening term; M_λ0 is the synthetic flux at the normalisation wavelength; ⊗ is the convolution operator; G(v_*, σ_*) is a Gaussian distribution centred at velocity, v_* and with a dispersion, σ_*; and x is the population vector. Each component x_j (j = 1, ..., N_*) represents the fractional contribution of the SSP with age, t_j, and metallicity, Z_j, to the model flux.

As previously mentioned, we use the CB07 base spectra for the fittings. The N_* = 100 SSPs comprises three metallicities: Z = 0.2, 1, and 2.5 Z_⊙, and 15 ages, from t = 0.001 to t = 13 Gyr. All SSPs are normalised to 1 M_⊙ at t = 0. Their spectra were computed with Padova-2004 evolutionary tracks models, and Chabrier (2003) IMF (0.1 M_⊙< M< 100M⊙).

The CB07 base comprises SSPs with the same metallicities and age range than the BC03 base (Bruzual & Charlot 2003). The difference between them are the TP-AGB spectra that CB07 adds in the base and the use of MILES-2007 evolutionary tracks models instead of Padova-2004. For the preliminary results in this work, we performed the SSP fitting on the CB07 base spectra.

STARLIGHT takes as its input: the wavelengths (λ), the observed spectrum (O_λ, resampled, rest-framed and corrected by MW extinction), the errors (e_λ), a base-master file with the SSPs, and a mask file. The mask file contains the regions of the spectra that we do not want to model with STARLIGHT, such as emission lines, artefacts, and holes in O_λ. Once STARLIGHT is running, it goes on to build the synthetic spectra M_λ of the form of Eq. (A.1) and to find the one that best fits the observed O_λ. The outputs files contains the population mixture of the best fit; the x_j, Z_j, t_j, (L/M)_j, stellar masses and the percentage of each component, among the synthetic spectrum M_λ.

Once we run STARLIGHT with all spaxels of our MUSE datacubes (∼100.000 spectra per cube), we can build a ’nebular gas’ or ’residue’ datacube of the form O_λ − M_λ, containing the nebular emission of the object (see Fig. A.1).

Fig. A.1. Example of a SSP fitting for a simple spaxel of the galaxy NGC 863. The observed, rest-framed and resampled spectrum is coloured in green. The black spectrum is the best fit synthetic spectrum we get from the STARLIGHT SSP fitting, and the blue spectrum is the nebular emission of the spaxel, obtained from the residue of the fit, subtracting the synthetic spectrum to the observed.

Then, using the Python MPFIT⁵ module, we perform a Gaussian fitting around the emission lines of interest for all the spaxels of the nebular gas cube. This module performs a multipeak Gaussian fit for each individual line, with $\sigma =FWH{M}_i/(2\sqrt{2ln2})$, with FWHM_i as the MUSE initial full width at half maximum, corresponding to 3 Å. The module search the peak of the line in a 10 Å (200-300 km/s) range centred in a central wavelength λ_c and performs the Gaussian fit in a range between 0.5σ − 2σ. The flux of the i, j spaxel for the λ line is then defined as $f_{i,j}(\lambda )=\sqrt{2\pi }{I}_{\mathit{peak}}$, and I_peak the peak of the Gaussian fit (see Fig. A.2).

Fig. A.2. Example of the multipeak fitting of Hα and [N II] lines (top). The purple curve is the Gaussian curve fitted. Example for Hβ (bottom). Blue curve is the Gaussian curve fitted. In both cases, black dashed line is the input spaxel spectrum.

Appendix B: EW_Hα for different SSPs

Figure B.1 shows the different values of EW(Hα) obtained in function of the type of SSP fitting. We used an alternative version of STARLIGHT, called PHOTOMETRICSTARLIGHT (Werle et al. 2019), that combines spectroscopic and photometric constraints to perform SSP synthesis using photometric points to extrapolate the model spectra to the bluer part, solving the problem of the lack of blue constraint in the MUSE spectra due to its spectral coverage. We used the spectrum of a bin with EW(Hα)≈1 Å of the NGC863 observed cube binned performing the methodology previously exposed. We selected three different photometric constraints to perform this fittings; near-ultraviolet (NUV, λ_eff = 2310 Å) and far-ultraviolet (FUV, λ_eff = 1528 Å) from All-Sky Survey of the Galaxy Evolution Explorer (AIS-GALEX; Bianchi et al. 2017) and u-SDSS (λ_eff = 3543 Å) from Sloan Digital Sky Survey (SDSS) DR16 (Ahumada 2020). Then, PHOTOMETRICSTARLIGHT reads the spectrum and the AB magnitudes obtained from the GALEX and SDSS images and performs the SSP fitting setting a certain base spectra (see Sect. 3.2). We make this fittings selecting both BC03 and CB07 bases and selecting all different combinations of photometric constraints with the three photometric points, obtaining then 12 different model spectra. Considering this effect as a cause of systematic error, we can get an error as high as 80% due to the selection of stellar populations considered in the fitting and the type of SSP fitting performed, at low EW(Hα) regimes.

Fig. B.1. EW(Hα) measurement for a NGC863 bin of EW(Hα)≈1 Å using 12 different SSP fittings with the Photometric STARLIGHT code (Werle et al. 2019). Red values correspond to the those obtained from the synthetic spectrum after setting the BC03 base. Idem for the green, using the CB07 base. The x-axis represent the six different photometric constraints; 1) No photometric constraints; 2) GALEX-NUV and GALEX-FUV photometric points added to the fitting; 3) only u-SDSS point added; 4) u-SDSS and GALEX-FUV points added; 5) u-SDSS and GALEX-NUV points added; 6) u-SDSS, GALEX-FUV, and GALEX-NUV points added to the fitting. We can see that the type of fitting affects the EW(Hα) measurement by as much as 80%.

All Tables

All Figures

	Fig. 1. RGB synthetic images of the galaxies selected for our sample. The images are constructed as a composition of the [S II]λ6716 (green), Hα (red), and [O III]λ5007 (blue) spaxel-by-spaxel emission maps obtained following the methodology explained in Sect. 3.2.
In the text

Fig. 2. Example of an adaptive binning of NGC863 using S/N([S II]) = 10 as target. The first figure represent the map numbers (mapN) from Li et al. (2023) algorithm. The second figure is the segmentation map obtained from our algorithm, where the bin index goes from k = 1 to Nbins. The third figure is a close-up of the segmentation map plotting the borders of the bins over the Hα map. It is noticeable that the H II regions, whose S/N is higher, are not binned, maintaining the structure pixel by pixel, and the bins are getting bigger as we move out of the H II regions.

In the text

Fig. 3. NGC863. Left panel: [S II] emission line map obtained from the Gaussian fitting of the [S II] line from the nebular gas cube spaxel-by-spaxel. The nebular cube was obtained as the difference between the observed cube and stellar cube as described in Appendix A. Central panel: binned [S II] emission line map, from the binned nebular gas cube, after performing an adaptive binning with target S/N([S II]) = 10 to the observed cube, as described. Both panels are close-up of the same region than Fig. 2. Right panel: distributions of S/N of the [S II] line for NGC863, calculated with the Eq. (1). Red distribution corresponds to the S/N([S II]) measured in the observed cube, namely, the S/N([S II]) of the spaxels. Green distribution corresponds to the S/N([S II]) measured in the binned cube, namely, the S/N([S II]) of the bins. The vertical dark green line marks the mean value of the S/N([[S II]]) of the bins, showing that the adaptive binning technique fulfils the goal of reaching a target S/N([S II]) = 10, on average.

In the text

Fig. 4. Steps to get the lower and upper DIG limit for NGC863 galaxy. The upper-left part corresponds to the DIG mask obtained by masking the H II candidates given by the centroids and radii obtained from PYHIIEXTRACTOR to the binned Hα map. The upper-right part corresponds to the binned Hα map after applying a 3σΣ(Hα) cut-off. The lower-left part is the mask of the lower limit of the DIG, obtained as the combination of the upper-left and upper-right masks. The lower-right part corresponds to the upper limit, assuming a constant non-zero flux level above the H II regions, interpolating flux around these regions.

In the text

Fig. 5. Distribution of the DIG fraction (f(Hα)DIG/f(Hα)total) for each galaxy of the sample. Red area represents the cumulative distribution between the upper and lower DIG limits, with a step of 0.05R25, being the solid red lined the mean value. Blue area is the radial distribution of the DIG obtained integrating the flux between rings of width 0.05R25 and with a step of 0.05R25, being the solid blue lined the mean value. R25 is the isophote at the blue brightness of 25 mag arcsec−2. This parameter is obtained from Hyperleda database.

In the text

	Fig. 6. Radial distribution of ΣHα for each galaxy, obtained by measuring the surface brightness of a series of de-projected annulus centred at the nucleus of each galaxy with a width of 0.05R25. The distributions are performed for the entire galaxy (blue), the lower DIG limit (green) and H II regions (red). The last panel represent the mean distributions for the sample.
In the text

Fig. 7. BPT diagnosis of the sample for the DIG bins that verifies that the relative errors of [O III] line bins are below 40% and its S/N above 3. The last panel displays the seven prior panels together. The central parts of NGC863, NGC3393, and NGC692 were excluded in all panels due to their strong AGN emissions. Each contour encloses the 10% of the points, with every point a single bin. Black lines are given by Kewley et al. (2001) (solid) for the classic H II regions photoionisation and AGN demarcations. The dashed-dot line depicts one of the 2D diagnostics developed by Kopsacheili et al. (2020) for the separation of shock excited (e.g. supernova remnants) from photoionised regions (e.g. H II regions). Coloured lines represent the models of photoionisation by fast shocks from Allen et al. (2008). The blue and red models illustrate photoionisation where only front shocks occur and when pre-ionisation by a precursor is taken into account. The solid model curves plotted represents shocks winds of 200, 400, 500 and 1000 km s−1, and dashed model curves represents magnetic field intensities of 0.0001, 1.0, 5.0 and 10. Red and blue arrows represent the direction of increasing wind velocity in each model.

In the text

Fig. 8. Radial distribution of the [S II]/Hα ratio of the DIG bins (green), and H II regions according to our definition (red), obtained as the ratio of the mean values of the respective fluxes within a ring of width 0.05R25 with a step of 0.05R25 from 0.2R25 to R25. We are only considering those bins whose relative error in [S II] and Hα lines are lower than 40%. To avoid the strong AGN emissions, we start all the distributions from 0.2R25. In general we see a tendency to increase the ratio as we move to outer regions of the galaxy. Besides, the ratio is in all cases higher in DIG regions, as we expected from the literature (Haffner et al. 2009).

In the text

	Fig. 9. Three examples of observed (green), model (black) and nebular (blue) spectra of different EWHα = 0.6, 1 and 3 Å. All observed spectra are integrated spectra from NGC863 bins, the model spectra are obtained from STARLIGHT following our methodology and the nebular spectra are obtained subtracting the model to the observed (see text). Gold lines represent the [N II] doublet.
In the text

Fig. 10. Distribution of Hα/Hβ ratio for three different EWHα regimes. Each distribution represent the Hα/Hβ flux ratio for all bins between 6 Å < EWHα ≤ 10 Å (green), 3 Å < EWHα ≤6 Å (blue), and EWHα ≤3 Å (red), for the seven binned galaxies listed in Table 1. Black vertical dashed line presents the theoretical ratio of Hα/Hβ = 2.87 (Osterbrock & Ferland 2006). We can see that a substantial part of the bins with EWHα ≤ 3 Å shows a non-physical ratio (< 2.87).

In the text

	Fig. 11. EWHα distributions for all the bins of the sample. Green histograms represents the distribution for the DIG emission bins, and red for the H II regions emission bins. Last histogram represents the distribution for the seven galaxies.
In the text

	Fig. 12. Mean radial Hα/Hβ ratio distribution of the sample. The distributions are obtained as in Fig. 8. Red represent the distribution for the H II regions bins, and green for the DIG bins. Both distributions decrease radially, on average, but the extinction is always lower in the DIG regime.
In the text

	Fig. A.1. Example of a SSP fitting for a simple spaxel of the galaxy NGC 863. The observed, rest-framed and resampled spectrum is coloured in green. The black spectrum is the best fit synthetic spectrum we get from the STARLIGHT SSP fitting, and the blue spectrum is the nebular emission of the spaxel, obtained from the residue of the fit, subtracting the synthetic spectrum to the observed.
In the text

	Fig. A.2. Example of the multipeak fitting of Hα and [N II] lines (top). The purple curve is the Gaussian curve fitted. Example for Hβ (bottom). Blue curve is the Gaussian curve fitted. In both cases, black dashed line is the input spaxel spectrum.
In the text

Fig. B.1. EW(Hα) measurement for a NGC863 bin of EW(Hα)≈1 Å using 12 different SSP fittings with the Photometric STARLIGHT code (Werle et al. 2019). Red values correspond to the those obtained from the synthetic spectrum after setting the BC03 base. Idem for the green, using the CB07 base. The x-axis represent the six different photometric constraints; 1) No photometric constraints; 2) GALEX-NUV and GALEX-FUV photometric points added to the fitting; 3) only u-SDSS point added; 4) u-SDSS and GALEX-FUV points added; 5) u-SDSS and GALEX-NUV points added; 6) u-SDSS, GALEX-FUV, and GALEX-NUV points added to the fitting. We can see that the type of fitting affects the EW(Hα) measurement by as much as 80%.

In the text