© The Authors 2024

1. Introduction

Non-star-forming, passively evolving galaxies are often referred to as quiescent galaxies (QGs). These galaxies are characterized by their red optical colors and are dominated by evolved stellar populations. Originally, QGs were thought to be poor in cool and warm gas. Nevertheless, due to improvements in observations, it is found that diffuse warm ionized gas is present in such galaxies (Knapp et al. 1985; Phillips et al. 1986; Kim 1989; Buson et al. 1993; Goudfrooij et al. 1994; Knapp & Rupen 1996; Macchetto et al. 1996; Zeilinger et al. 1996; Sarzi et al. 2006; Davis et al. 2011; Kehrig et al. 2012). As revealed by Phillips et al. (1986), the spectrum of such gas exhibits strong emission lines from low-ionization ions, similar to low-ionization nuclear emission-line regions (LINERs) (Heckman 1980). Ionization of LINERs was widely believed to be related to the central active galactic nuclei (AGNs) in the galaxy (Ferland & Netzer 1983; Halpern & Steiner 1983; Ho 1999, 2009; Kauffmann et al. 2003; Groves et al. 2004; Kewley et al. 2006) and the ionized interstellar medium (ISM) in QGs were therefore thought to have the same ionization mechanism.

In recent years, observations with enhanced spatial and spectral resolution have led to discoveries of similar emission regions extending outside the nuclear region of galaxies (Sarzi et al. 2010; Yan & Blanton 2012; Singh et al. 2013; Belfiore et al. 2016), raising a request for a different ionization mechanism. Following Belfiore et al. (2016), these regions are referred to as low-ionization emission regions (LIERs). While the dominant ionizing source remains a mystery, two main candidates have been proposed, which include photoionization by post-asymptotic giant branch (pAGB) stars (Binette et al. 1994; Stasińska et al. 2008; Sarzi et al. 2010; Cid Fernandes et al. 2011; Yan & Blanton 2012; Singh et al. 2013; Belfiore et al. 2016; Yan 2018a; Byler et al. 2019) and collisional ionization by fast shocks (Heckman 1980; Dopita & Sutherland 1995; Groves et al. 2004; Allen et al. 2008).

The spectral emission line ratio is one of the tools that has been used to investigate the physical condition of the ISM. One of the most typical examples is the Baldwin-Phillips-Terlevich (BPT) diagram (Baldwin et al. 1981; Veilleux & Osterbrock 1987). With different line ratios plotted against each other, we can characterize galaxies by their ionization mechanisms. Nonetheless, shocks and post-AGB stars can produce similar ratios in emission lines used by the BPT diagrams. Models of the two mechanisms cover similar areas on these diagrams, and thus do not provide any leverage to distinguish the dominant ionization source. In addition to the BPT diagrams, different emission line ratios can be used to measure the electron density, metallicity, ionization parameters, or temperature of specific ion species. The physical conditions of the ionized ISM are expected to depend on its ionization sources. Given that collisional ionization by interstellar shocks and photoionization by pAGB stars predict different temperatures in the gas, measuring the electron temperatures provides a way to distinguish between the two mechanisms. In general, gas ionized by interstellar shocks is expected to have a higher temperature than photoionized gas. The temperature of different ions can be measured by the ratio of auroral lines to strong lines from the gas.

For this study we made use of spectroscopic data from the Sloan Digital Sky Survey IV (SDSS-IV) Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey to detect temperature-sensitive line ratios: [O III] λ4363/[O III] λ5007, [N II] λ5755/[N II] λ6583, [S II] λ4068, 4076/[S II] λ6716, 6731, and [O II] λ7320, 7330/[O II] λ3727, 3729. Following procedures similar to those used in Yan (2018a), we made use of stacking techniques to increase the S/N of the observed spectrum. Because LIERs have weak emission lines in general and the temperature-sensitive auroral lines are even weaker, we needed very accurate subtraction of the stellar continua, in addition to high S/N, to detect these weak lines. Full-spectrum fitting using stellar population models turned out to be insufficient. Empirical spectra from the same dataset produced a much better result. Thus, we selected those quiescent spaxels without emission lines, matching them to those spaxels with emission lines in terms of stellar population. From these line-free spaxels, we stacked their spectra together to get modeled stellar continua. After subtracting the stellar continua, the fluxes of the emission lines were measured and the temperature-sensitive line ratios were compared to the prediction of models.

This paper is organized as follows. In Sect. 2 we briefly describe the data. In Sect. 3 we present the methodology of our spaxel selection, stacking, and line flux measurement. In Sect. 4 different physical properties of the data are computed. The results are then compared to the models assuming different ionization sources in Sect. 5. In Sect. 6 we summarize our main findings.

All absolute magnitudes in this paper are calculated with H₀ = 100h km s⁻¹ Mpc⁻¹, so they should be interpreted as M – 5 log₁₀h, in the AB magnitude system.

2. Data

We used data from MaNGA survey (Bundy et al. 2014; Yan et al. 2016), which is part of the SDSS-IV (Blanton et al. 2017). The MaNGA instrument contains 17 integral field units (IFUs) which are composed of 19–127 fibers (Drory et al. 2015). These IFUs cover diameters between 12 and 32 arcsec in the sky. Using the 2.5 m Sloan telescope (Gunn et al. 2006) at the Apache Point Observatory, the light from the galaxies is guided through the fibers and fed into the BOSS spectrograph (Smee et al. 2013), where they are split into blue and red channels. In each channel, the light is dispersed producing spectra with a resolution of R ∼ 2000 and a wavelength range of 3622–10 354 Å.

The SDSS Data Release 17 (Abdurro’uf et al. 2022) contains the complete release of the MaNGA survey, consisting of over 10 000 nearby galaxies. Most of the targets are covered with radial coverage between 1.5 to 2.5 effective radii (R_e) (Wake et al. 2017). Since all these galaxies have a redshift of 0.01 < z < 0.15, the wavelength coverage of MaNGA covers all the emission lines of interest in our project. The spectroscopic data we used is processed by the MaNGA Data Reduction Pipeline (DRP) (Law et al. 2016, 2021). The DRP converts raw data into sky-subtracted, flux-calibrated data cubes. In this work, we also use the spectral indices, emission line properties and kinematics derived by the MaNGA Data Analysis Pipeline (DAP) (Westfall et al. 2019; Belfiore et al. 2019). The DAP analyzes the data from DRP to produce a range of physical properties of the spaxels. Each spaxel is 0.5″ in size. The typical point spread function of the MaNGA data is around 2.5″ to 3″. MaNGA employs dithered observation to produce near-critical sampling of the PSF (Law et al. 2015).

In the following section, Gaussian flux, Equivalent Width (EW) of emission lines, and spectral indices provided by the DAP are used.

3. Methods

This section consists of four subsections. Section 3.1 describes the LIERs spaxel selection process. The construction of the corresponding stellar continuum for the selected spaxels is discussed in Sect. 3.2. Section 3.3 talks about the de-redshifting and stacking process. The method used to measure emission line fluxes and estimate uncertainties is described in Sect. 3.4.

3.1. Spaxel selection

To select quiescent galaxies with negligible ongoing star formation, we constructed a plot with D_n4000 versus absolute magnitude in the r band, which is shown in Fig. 1. The absolute magnitudes are from the NASA-Sloan Atlas (NSA) catalogue. D_n4000 is a dust-insensitive stellar age indicator, which can help us identify quiescent galaxies without including dusty star-forming galaxies. We adopted the definition of D_n4000 from Balogh et al. (1999), who define this index as the ratio of average flux density between 4000–4100 Å and 3850–3950 Å. In Fig. 1, the galaxies are divided into two groups, the blue cloud and the red sequence. Located at the bottom right, the blue cloud represents active galaxies hosting young stellar populations. They have lower D_n4000 values and on average fainter in the r band. On the contrary, the quiescent galaxies with older stellar populations are at the top left region.

Fig. 1. Dn4000 vs. absolute magnitude in the r band. The values of Dn4000 are from the MaNGA DAP and the absolute magnitudes in the r band are from the NSA catalogue. The gray region shows the density distribution of all galaxies, while the orange contours are drawn at 10%, 30%, 68%, 95%, and 99% of all galaxies. Galaxies between the red dashed lines are selected as our quiescent galaxy sample. Galaxies with invalid Mr and Dn4000 are not plotted.

To select our targeted galaxies, we first drew a line following the ridge-line of the quiescent galaxies distribution. Then, we shifted it upward and downward such that it covered the whole quiescent sequence. The formulae of the boundary lines are given by

(1)

and

(2)

We only included spaxels with mean g-band weighted signal-to-noise ratio per pixel greater than 15 in the selected galaxies, yielding a total of 492 275 spaxels. We then applied a zero point correction as described in Sect. 3.3 of Yan (2018a). The rationale is to remove any systematic error in the measurements of small EW by DAP. We defined EW to be

(3)

where F_s is the flux of the observed spectrum and F_c is the flux of the continuum. For emission lines, the EW is positive.

To measure the zero-point of EW, we adopted the assumption that EWs of strong emission lines correlate with each other. This assumption applies because LIERs have fairly uniform line ratios and also very uniform stellar continuum shapes. Here we considered [O II] λλ3727, 3729, [N II] λ6583, [S II] λλ6718, 6731, [O III] λ5007, Hβ, and [O I] λ6300. If an emission line EW is consistent with zero, the EW value of each of the other lines of this spaxel should also be among the weakest. To evaluate the real zero-point of a certain emission line, we looked at other lines. For each of the other lines, spaxels are divided into different percentiles according to their EW. Spaxels that lies below a certain percentile in all emission line except the targeted one were selected. The median EW value of the targeted line among the selected spaxels was then plotted against percentiles in Fig. 2. The median EW of the 35th percentile was chosen as the zero-point since it represents the flattest part in most of the plots, showing a stable zero-point without much disruption by noise. Figure 3 shows the distribution of EWs after the correction. One can see from these figures that, after the correction, the most densely populated regions are around (0,0) in all EW versus EW diagrams. Besides correcting the EWs for the zero-point offsets, we also correct the Gaussian emission-line fluxes accordingly.

Fig. 2. Median EW value of the selected line among spaxels with EW lower than different percentile in other strong lines. The horizontal lines denote the EWs corresponding to 35 percentiles, which are the zero points we chose.

Fig. 3. Different strong-line EWs against Hα EW after zero point correction. The blue cross denotes the origin. The grayscale represents the density of the spaxels, while the orange contours enclose 10%, 68%, 95%, and 99% of the spaxels in each panel. The 10% contour is used to show that most of the data are centered at the zero value. The line in the first subplots represents the cut we applied to exclude contamination from transition objects and dusty star-forming galaxies. All spaxels on the right side of the line are removed from the next step.

As seen in the [O II] EW versus Hα EW panel in Fig. 3, there is another group of spaxels to the right of the main group. This might be due to contamination from other transition objects and star-forming galaxies which appear red due to dust. According to Yan et al. (2006), these objects should be excluded using an [O II]-Hα EW cut. Hence, we excluded these objects by the cut

(4)

which is shown in Fig. 3. This removed 25 526 spaxels with the remaining 466 749.

To prevent possible contamination of spaxels coming from Seyfert galaxies, we applied another cut to the sample using the definition of Seyfert galaxies in Heckman (1980). We selected spaxels with a fractional error¹ of [O III]/[O II] flux lower ratio than 0.3 and excluded those with [O III]/[O II] flux ratio greater than 1. This removed 2676 spaxels and 464 073 spaxels were remaining. The distributions of EWs after removing these dusty and Seyfert spaxels are displayed in Fig. 4.

Fig. 4. Different strong-line EWs against Hα EW after excluding dusty and Seyfert spaxels. The gray regions and contours are generated in the same way as those in Fig. 3.

Next, the spaxels are separated into two classes, strong-line samples and zero-line samples. The former contains spaxels with strong emission from the gas and will be regarded as the data to be stacked and measured. Gas emission is very weak in the latter class, making these spaxels suitable candidates for constructing stellar continua. To select strong-line samples, we constructed a parameter named total EW using all EWs from the targeted lines. Here we considered Hα, [O II] λλ3727, 3729, [N II] λ6583, [S II] λλ6718, 6731, and [O III] λ5007 they are the strongest lines in the optical spectrum. The Total EW is defined as follows:

(5)

The coefficients were chosen to match the slope of the distribution in Fig. 3 so that the Total EW increases along the sequence in each diagram. Spaxels with the highest 20% of total EW index are often classified as the strong-line sample, and the spaxels are regarded as strong-line spaxels. In the selection, we only included the spaxels with unmasked EW values on all of the above emission lines. There are a total of 90 359 strong-line spaxels.

For a strong-line spaxel to be stacked, we need to know its stellar and gas kinematics, which can differ from each other. We describe how we deal with this in Sect. 3.3. Therefore, we only included strong-line spaxels with valid line-of-sight stellar velocity and Hα velocity from the DAP. We also required the difference between the two velocities to be less than 300 kms⁻¹, to remove the most extreme velocity offsets. In total, these selections excluded 2390 strong-line spaxels with 87 969 remaining.

Since metals are effective coolants in the ISM, variations in metallicity among different spaxels can cause large discrepancies in measured temperatures. Therefore, we divided the strong-line sample into three bins, namely the high, mid, and low [N II]/Hα bins. We plotted two metallicity indicators, [N II]/[O II] and [N II]/Hα, against each other in Fig. 5. Due to the large wavelength difference of [N II] and [O II] strong lines, dust extinction can significantly affect the [N II]/[O II] values. Hence, we divided the sample only based on dust insensitive [N II]/Hα values. We only selected spaxels with a fractional error for [N II]/Hα flux ratio lower than 0.3. The selected spaxels were then allocated into three equal-number bins, with 28 098 spaxels in the high [N II]/Hα bin and 28 097 spaxels in each of the remaining bins. The [N II]/Hα values separating the three bins are –0.0247 and 0.1374, respectively. Although our rationale was to use [N II]/Hα as a metallicity indicator, it actually may be sensitive to other parameters, such as the ionization parameter in photonionization or shock velocity in the case of shocks. Therefore, the bins are referred to as [N II]/Hα bins instead of metallicity bins throughout the paper.

Fig. 5. [N II]/[O II] vs. [N II]/Hα for all spaxels in the strong-line sample. Instead of a narrow sequence, there is noticeable scatter in this relation. The two blue dashed lines separate the spaxels into high, mid, and low [N II]/Hα bins. The data is meshed into 100 × 100 bins. The grayscale represents the density of all selected spaxels in each bin, while the orange contours enclose 68%, 95%, and 99% of the spaxels. The blue error bars in the top left corner denote the typical 1σ uncertainty of the line ratios. The extinction vector is given at the top right corner.

Spaxels with negligible gas emission are selected as the zero-line sample which provides the stellar continuum to be subtracted from the spectra with strong gas emission. Here we considered the same emission lines as the strong-line sample selection, which include Hα, [O II] λλ3727, 3729, [N II] λ6583, [S II] λλ6718, 6731, and [O III] λ5007. Following Yan (2018a), we constructed a multi-dimensional ellipsoid around 0 EW. To select points that are consistent with 0, we need to know the typical noise in the measurement. Since we are measuring the EWs of emission lines, all negative EW values should be due to systematics and random errors. We used these negative EW to estimate the typical uncertainty for each emission line by calculating their root-mean-squared (RMS) value from 0. The EW of each line was normalized by this typical uncertainty. We then squared them and summed them up, then took the square root to calculate an uncertainty-normalized distance away from the origin. For zero-line sample selection, we required this distance in the multidimensional space to be less than 2.5. In total, there are 185 604 zero-line spaxels.

3.2. Matching of strong- and zero-line samples

In order to measure the weak auroral lines, we needed an accurate subtraction of the stellar continuum. We selected and stacked spectra from the zero-line sample to construct a stellar continuum matching that of the strong-line stack. In order to get a good matching, we tried to select spaxels from the zero-line sample matching the stellar population of each spaxel in the strong-line sample. Each strong-line spaxel contains a spectrum composed of a stellar continuum and gas emission lines. We used four proxy parameters, D_n4000, stellar velocity dispersion, surface luminosity in the r-band, and total spectral broadening to match the stellar continuum shapes.

1. D_n4000: The D_n4000 values are from the DAP, following the definition mentioned in Sect. 3.1. D_n4000 is a stellar age and metallicity indicator. Both of them largely affect the shape of the stellar continuum.

2. Stellar velocity dispersion (σ_*): DAP provides the measurement of the line-of-sight stellar velocity dispersion. Because σ_* correlates with the stellar metallicity [Fe/H], and especially [Mg/Fe] (Graves et al. 2009; Bernardi et al. 2003), matching it can help match the stellar continuum.

3. Surface luminosity in the r-band: The r-band flux per pixel is from the DRP. We corrected for the galactic reddening effect using the E(B-V) value given by the DAP. As all spaxels in our sample have similar colors, matching the surface luminosity in r-band can help match the stellar mass surface density between spaxels, which can facilitate the matching of properties affecting the stellar continuum.

4. Total spectral broadening: Total spectral broadening is the quadratic sum of stellar velocity dispersion σ_* and the instrumental line spread function (LSF). The instrumental LSF values are the post-pixelized LSF given by the DRP. We first convert the post-pixelized LSF value into velocity space, then take the median across the whole spectrum before combining with σ_*. After that, we convert back the median value into wavelength space and quadratically summed it with σ_*. Matching the total spectral broadening matches the absorption line width between the strong and zero-line spaxels.

We excluded strong and zero-line spaxels with any of the above parameters masked by the DRP or DAP. There are 183 736 zero-line spaxels. The high, mid, and low [N II]/Hα bins of the strong-line sample have 27 715, 27 883, and 27 924 spaxels, respectively. Next, we normalized these four parameters using their 5th and 95th percentiles as 0 and 1, respectively. We used a K-Dimensional tree for matching the spaxels. A four-dimensional space was created using the four normalized parameters. For each strong-line spaxel, we found the 2 non-occupied zero-line spaxels that are closest to it in the 4-D space. In each iteration, each strong-line spaxel is matched with the two closest zero-line spaxels. It is possible that a zero-line spaxel will be matched with more than one strong-line spaxel. In such a case, only the matching between the zero-line spaxel and its closest strong-line spaxel will be kept. Strong-line spaxels with fewer than two matches will be kept for the next iteration. The maximum iteration is 25, and the remaining strong-line spaxels will be abandoned. After all the iterations, the matching distance is required to be less than 0.05 in the four-dimensional space, or the trio will be excluded. There are 25 139, 26 435 and 25 229 spaxel trios in the high, mid, and low [N II]/Hα bins respectively.

3.3. Drizzling and stacking

The position of a certain emission line in the spectrum depends on the redshift of the source spaxel. Within the region covered by the spaxel, ionized gas velocity can be different from the stellar kinematics. Therefore, the stellar absorption lines and the gas emission lines are shifted by different extents. To detect emission lines in the stacked spectrum of each [N II]/Hα bin, the spectra of the spaxels must be shifted according to their gas kinematics. In order to get purely gas emission spectra, the stellar continuum of the strong-line spaxel has to be subtracted. Using the zero-line stacked spectrum as a stellar continuum template, the zero-line spectrum has to be shifted to the same position as the strong-line spaxels’ stellar continua.

We defined a velocity offset as the difference between the stellar and ionized gas velocities. The redshift value, stellar velocity, and gas velocity are all from the DAP. The Hα velocities are used as the gas velocities of the spaxels. We divided all the strong-line spaxels into 25 bins of velocity offset ranging from –300 km s⁻¹ to 300 km s⁻¹. The distributions of strong-line spaxels in the velocity offset bins are given in Fig. 6. To align with the strong-line spaxel, the matched zero-line spectra are first shifted by their redshifts and stellar kinematics. Then they are shifted by the center value of the velocity offset of the velocity-offset bin in which its associated strong-line spaxel is located. In practice, the shifts were done together to minimize errors introduced by the shifts. The overall redshift value is given by Eq. (6),

(6)

Fig. 6. Distribution of velocity offsets of all strong-line spaxels in each [N II]/Hα bin. The velocity offset is defined as the difference between the stellar velocity and the gas velocity. There are a total of 25 velocity-offset bins covering the range from –300 km s−1 to 300 km s−1.

where z_p is the overall redshift, v_* is the stellar velocity, v_off is the center value of the velocity offset bin and z is the redshift of the galaxy. The term v_* − v_off means that the spaxels are aligned by the gas velocity, which corresponds to the positions of the emission lines in the spectrum.

All spectra are first corrected for galactic extinction using the extinction curve of Fitzpatrick (1999) with R_V = 3.1 before the drizzling. MaNGA galaxies are nearby with redshift values smaller than 0.15. Since the covered wavelength is up to 10 354 Å, the maximum wavelength after deredshifting is set to be 8500 Å which allows us to measure all the targeted emission lines. We stored the shifted spectra in a common wavelength grid, evenly spaced logarithmically with the same spacing as the MaNGA logcube. The process of assigning flux among the output pixel is referred to as drizzling.

In the following context, a pixel refers to the region between two wavelengths in the wavelength array of the MaNGA logcube. The re-distribution of flux among the output pixels should conserve flux locally. In the drizzling process, the original wavelength array should shrink by a factor of (1 + z), and the flux in each pixel is multiplied by this amount to conserve the total flux. The shifted pixel could fall across the boundary of two pixels in the output array. We allocated the flux of each shifted pixel according to its overlapping fractions with the final pixel grid. Since both the input and output array are logarithmic, the proportion depends only on the z_p value of the spaxel and it remains the same across the whole wavelength array.

All the spectra after de-redshifting are normalized by their median flux values of the region between 6000–6100 Å. Within each [N II]/Hα bin, strong-line spaxels and zero-line spaxels in each velocity offset bin are stacked together respectively. Similar to Yan (2018a), there are small discrepancies between the overall shapes of the stacked strong-line and zero-line spectra. We followed the same procedure to solve the issue. We first divided the strong-line spectrum by its corresponding zero-line spectrum and applied a B-spline fitting on the resulting curve with a spacing of 400 pixels between the breakpoints. Emission lines are masked in the fitting as the goal is just to get the overall shape difference. The stacked zero-line spectrum is then multiplied by the B-spline curve to match the shape of the stacked strong-line spectrum. The gas emission spectra are obtained by subtracting the shape-corrected zero-line spectrum from the strong-line spectrum. Finally, the gas emission spectra of all velocity offset bins are stacked together, weighted by the number of spectra stacked in that velocity offset bin. Finally, the gas emission spectrum of each of the three [N II]/Hα bins is obtained. Figure 7 shows the strong-line and zero-line spectra of the three bins. The gas emission spectra are displayed in Fig. 8.

Fig. 7. Stacked strong-line spectrum (blue) and zero-line spectrum (orange) of the three [N II]/Hα bins. The zero-line spectrum has been corrected for overall shape using B-spline and is vertically shifted downward for visibility.

Fig. 8. Continuum-subtracted residual spectrum of the three [N II]/Hα bins. All the lines in the spectrum should originate from gas emissions.

3.4. Emission line measurement

We followed the same method in Yan (2018a) and Yan et al. (2006) for flux measurement. Two sidebands and a central region are selected for each line. A linear fit of the sidebands is used as an estimation for the local continuum shape and is subtracted before integrating the flux inside the central region. We divided the emission lines into two groups, strong lines and weak lines. The definition of strong and weak lines with the corresponding regions defined can be found in Tables 1 and 2. Assuming all the lines have a symmetric profile, the central region is defined to be centered at the vacuum wavelength of the line in most of the singlet cases. A special case is Hε, which is treated as a doublet with the coinciding [Ne III] λ3968. For the doublets, the region extends to the left of the bluer line with the same amount to the right of the redder line. For [N II] λλ6548, 6583 and Hα, summed flux is not used and the regions are only defined to estimate the noise and fit the Voigt profile as described in the next paragraph. By visual inspection, the central region can cover the whole emission line profile in all of the lines.

Since [N II] λλ6548, 6583 are overlapping with Hα, we estimate the fluxes of the three lines with pseudo-voigt fits. The pseudo-voigt profile is a linear combination of a Gaussian and a Lorentzian profile which has one more degree of freedom in the fitting process than the Gaussian and Lorentzian profiles. We tried to fit all three profiles and the pseudo-voigt has the best estimation for the shape of the lines as it has one more degree of freedom than the other two profiles. Before the fitting process, the spectrum is first subtracted with the linear fit estimated from the sidebands. Singlets are fitted only by providing the vacuum wavelength as the initial guess. For the doublets, the stronger line position is provided as the initial guess, and the relative position of the other line is fixed. For [N II] and Hα, the three line profiles are fitted simultaneously as a combination of a singlet and a doublet. The relative strength of the two [N II] lines has a fixed ratio as required by quantum mechanics, where I₆₅₈₃/I₆₅₄₈ = 3.071 (Storey & Zeippen 2000) and I is the flux intensity of the line. In our following discussion, only the flux from [N II] λ6583 is considered. The flux of Hα and [N II] λ6583 used are both estimated by integrating the pseudo-Voigt profile, which is referred to as the Voigt flux. The measured Voigt fluxes of other strong lines with comparison to the summed fluxes are provided in Table 3. All the fluxes are normalized relative to the summed flux of Hβ. It can be seen that the pseudo-Voigt profiles can capture the fluxes of the lines well to better than 3% assuming the summed flux is the true line flux.

The uncertainties of the measurements are estimated by sliding boxes in the sidebands. After the correction of the linear fitline, the spectrum is locally flat with tiny fluctuations. We assumed the fluctuations in the central region are the same as in the sidebands. We slid a box, which has the same size as the central region, across the sidebands and measured the summed flux inside the box. The position of the box is shifted by one pixel each time and the RMS flux value of all the boxes is taken as the uncertainty of the measurement. For some of the lines, narrower sidebands are chosen to capture the local continuum shape. For lines with narrow sidebands in which fewer than 20 boxes can be fitted, the uncertainty is calculated by first computing the RMS value of all pixels in the sidebands, then multiply it with , where N is the number of pixels in the central region. The uncertainties and the fluxes of all the lines are given in Table 4. Among them, the fluxes of [N II] and Hα are given by the pseudo-Voigt fit.

In the stacking process, the spectra are normalized using the region of 6000–6100 Å. Therefore, the flux ratio between two lines A and B in the stacked spectrum becomes a weighted average of the per-spaxel line ratios with the weight proportional to the ratio of the line flux in the denominator to the median flux of 6000–6100 Å. This effect can be illustrated using the following equation:

(7)

Here F_A and F_B are the flux from lines A and B, and C_i is the continuum value used for normalization. We refer to Sect. 4.3 of Yan (2018a) for further details.

Figure 9 shows the continuum and emission line of the four targeted auroral lines. Since the auroral lines are weak, in this study a 2σ detection threshold is used. Shown in the bottom left panel, [O II] λλ7320, 7330 can be well detected in all the [N II]/Hα bins. Moreover, on the right end of the plotting range, a nickel line [Ni II] λ7378 is observed, especially in the low [N II]/Hα spectrum. A little bump is also observed in two other bins but their levels are close to the noise level in that region. Shown in the bottom right panel, [S II] λλ4068, 4076 is also 2σ detected in all the bins. The sidebands of this doublet are chosen to be narrower to capture the local shape of the continuum. The widths of the sidebands of the two lines mentioned are limited by [Ni II] λ7378 and Hδ respectively. In the top left panel, the flux of [O III] λ4363 increases with increasing metallicity. Although in the low [N II]/Hα plot, the continuum seems to be lower than 0 which leads to underestimation of flux, a narrower band was used and the result remains the same. As a result, [O III] λ4363 is detected in the high and mid [N II]/Hα bins but not the remaining bin. [N II] λ5755 is not detected in all of the bins despite showing a little bump. The continuum appears to be noisier in that region because the wavelength of the line, when redshifted, lies in the transition region of the blue and red spectrographs in MaNGA. The throughput curves can be found in Yan et al. (2015).

Fig. 9. Spectra around the four targeted auroral lines corrected by the fit line derived from the sidebands. The dark gray region is the central region, the light gray regions are the sidebands, and the orange line denotes the zero value. [O III] λ4363 is on the top left, [N II] λ5755 is on the top right, [O II] λλ7320, 7330 is at the bottom left, and [S II] λλ4068, 4076 is at the bottom right.

4. Data analysis

In this section, we discuss several physical properties deduced from the measured line ratios. In Sect. 4.1 we compute the dust extinction coefficient using the Balmer decrement. In Sect. 4.2 we use [S II] doublets to estimate the electron density. Last but not least, temperatures of different ions are computed in Sect. 4.3.

4.1. Balmer decrement and dust extinction

Under case B approximation, the intrinsic ratios of Balmer lines are fixed at certain values at a given electron temperature and density (Osterbrock & Ferland 2006). The observed ratios between the Balmer lines can be used to estimate the degree of dust extinction by comparing it with the theoretical value. We used the extinction curve from Fitzpatrick (1999) with R_V = 3.1. In Fig. 10, we plotted the ratios of the observed Balmer ratios to the theoretical ratios as a function of wavelength for the first five Balmer lines. Using Hα/Hβ ratio under case B approximation and at low-density limit at 10 000 K, the total extinction (A_v) were measured to be 0.73, 0.96, 0.73 magnitudes for high, mid, and low [N II]/Hα bins, respectively. The black curves in the figure show the predicted values for other wavelengths given these measured extinctions. The red points with error bars are the measured Balmer line fluxes relative to the Case-B values. Hε overlaps with [Ne III] λ3967. Since the flux ratio of [Ne III] λ3967 and [Ne III] λ3869 is locked by quantum mechanics, we subtracted 31% of the flux of [Ne III] λ3869 from the measured Hε flux. In the same way, part of Hζ overlaps with He Iλ3889 and He Iλ3889 is fixed to be 88% of He Iλ5876 theoretically. However, since He Iλ5876 and He Iλ3889 are far separated which makes them suffer from different dust extinctions, we had to correct the flux using the predicted dust extinction value before subtracting the expected He Iλ3889 flux from Hζ.

Fig. 10. Dust extinction value predicted by the Balmer decrement under case B approximation and low-density limit at 10 000 K. While Hα, Hβ, Hγ, and Hδ are more consistent, Hε is consistent with the prediction only in the low [N II]/Hα bin. Hζ values are not shown due to strong contamination and stellar continuum systematics.

The Hγ flux detected is pretty consistent with the predictions in all three [N II]/Hα bins. Hδ is on the prediction curve in the high [N II]/Hα bin, but below it in the mid and low [N II]/Hα bin. Hε fluxes are above the prediction curves. One of the reasons causing the offsets is that the right sideband (3898.86–3912.86 Å) of [Ne III] λ3869 is significantly higher than the left sideband due to imperfect stellar continuum subtraction. The fitline then has a positive slope making an underestimation of the flux. As a result, Hε flux is overmeasured. Since Hζ shares the similar sidebands with [Ne III] λ3869, the same effect largely affects the measurement of Hζ and therefore we do not trust the measurement on the weak Hζ line which is not included in the figure.

In conclusion, the well-measured high-order Balmer lines are consistent with the prediction of A_V by Hα and Hβ and the Fitzpatrick (1999) law. In Yan (2018a), higher-order Balmer lines show much less extinction than that derived from Hα/Hβ and that may be caused by noisier measurement and the difference of [N II]/Hα binning method. In spite of the consistency between the Balmer lines, the measured degree of dust extinction may not be applicable to the emission lines of other species. This will be discussed in more detail in Sect. 4.3. It is possible that different emission lines have different effective extinctions in these low spatial resolution observations, just like what Ji et al. (2023) found in star-forming galaxies. This finding is consistent with Yan (2018a).

4.2. Density derivation

We derived the electron density using the ratio of [S II] λ6731/[S II] λ6716. This doublet is optimal for measuring density as it has a very tiny dependence on temperature and dust extinction. We computed the density with the PYNEB package (Luridiana et al. 2015), which is the Python version of the IRAF package NEBULAR. Assuming the temperature to be 10 000 K, the measured electron densities are 45.20 cm⁻³, 57.52 cm⁻³, and 42.33 cm⁻³ for the three [N II]/Hα bins, respectively. All the measured densities lie in the low-density regime. Collisional de-excitation is not likely to happen in such diffuse regions and it is reliable to measure the temperature using the deduced densities. The measured densities are shown in Table 5.

The error of the measured electron density in Table 5 was estimated by repeating the line ratio measurement after adding uncertainties. For pixels in both the sidebands and the central region, we assumed the flux errors follow Gaussian distribution with the width of the Gaussian profile set to be 1σ standard deviation propagated from the stacking process. In each measurement, we randomly picked a number in the Gaussian distribution of each pixel as the flux. For each [N II]/Hα bin, the error was given by the highest and lowest derived density in 10 000 simulations.

4.3. Temperature derivation

For the same ion species, the temperature can be derived from the auroral-to-strong lines ratio. From our emission line measurement, we deduced temperatures of O⁺, O²⁺, S⁺, and N⁺. These ions populate different parts of the interstellar gas in different temperatures. Among the four, [O III] lines have the highest ionization potential and therefore represent more ionized regions. The other three ions, O⁺, S⁺, and N⁺, have similar ionization energy, and therefore emission lines from these originate from similar regions with similar temperatures.

Same as electron densities, temperatures were calculated using PYNEB. For undetected lines, the upper limit of temperature is given by assuming the flux value of the line to be 2σ. The upper and lower bounds of the measurement were calculated by adding 1σ uncertainty derived from the sliding box as described in Sect. 3.4. The uncertainty in electron density was also considered. Using the fact that the effect of electron density on temperature measurement is monotonic, we computed the temperatures using the maximum and minimum electron density given in Sect. 4.2.

All the temperatures derived without and with extinction correction are given in Tables 5 and 6, respectively. However, we believe that those without extinction correction are closer to the truth than those corrected using the extinction measured by Balmer decrement. The reason is explained below. Here we quote the temperatures using those in Table 5, without extinction correction. We discuss the results of different ions following the order of high, mid, and low [N II]/Hα bins:

1. O²⁺ zone: [O III] λ4363 was only detected in the high and mid [N II]/Hα bins as the flux of [O III] λ4363 significantly decreases with metallicity and becomes undetectable in the low [N II]/Hα bin. We note that some of the flux from the measured [O III] λ4363 flux may be contributed by [Fe II] λ4359, which can be estimated by another forbidden line [Fe II] λ4287. However, no visible bump was found around 4288 Å in all stacked spectra. Therefore, we did not attempt any correction on the strength of [O III] λ4363. The temperatures of the O²⁺ zone in the three [N II]/Hα bins were measured to be 3.17 K, 2.19 K and < 1.54 × 10⁴ K. Consistent with the findings of Yan (2018a), the temperature of the O²⁺ zone increases with metallicity, contrary to the common wisdom. With a higher abundance of metal, the ionized gas should be more efficiently cooled, resulting in a lower temperature. This measurement raises the question of whether metallicity is affecting ionic temperature through other mechanisms in LIERs, or whether our binning of [N II]/Hα might be tracing other properties.

2. N⁺ zone: Since [N II] λ5755 was not detected in all [N II]/Hα bins, we gave an upper limit for its temperature. Using 2σ value as the flux of [N II] λ5755, the derived temperature limits are < 0.9 × 10⁴ K, < 0.99 × 10⁴ K and < 1.18 × 10⁴ K in the three bins, respectively.

3. S⁺ zone: We had very robust detections on [S II] λλ4068, 4076. The temperatures derived are 0.94 K, 0.98 K and 1.02 K. The S⁺ zone temperatures do not show large variations among the three bins. It is not surprising that the temperature decreases when metallicity increases. Given that neutral nitrogen atoms have similar ionization energy to neutral sulfur atoms, the temperatures of the N⁺ zones may be also close to the S⁺ zones.

4. O⁺ zone: [O II] λλ7320, 7330 were significantly detected in all three bins. [O II] λλ7320, 7330 overlaps with another emission line [Ca II] λ7324 which can be estimated through [Ca II] λ7291. Nevertheless, [Ca II] λ7291 was not detected in all stacked spectra. Therefore, we assumed the contribution of [Ca II] λ7324 to the [O II] auroral lines is negligible. From the measurement, we derived the temperature of O⁺ zone to be 1.38 K, 1.33 K and 1.33 K in the three bins. This is consistent with the expectation that the O⁺ zone temperatures are higher than the S⁺ zone because neutral oxygen atoms have a higher ionization energy than neutral sulfur atoms.

The reason that we prefer the temperatures derived without extinction correction is that applying the dust extinction derived in Sect. 4.1 will lead to an unphysical scenario similar to that encountered by Yan (2018a). After correction, the temperatures of S⁺ zone would be much higher than the O⁺ zone, by 37–90% in different bins. Based on the first ionization energies of neutral oxygen atoms and neutral sulfur atoms, the S⁺ zone spans deeper toward the neutral gas compared to the O⁺ zone. The line in Fig. 11 shows the [S II] and [O II] line ratios when the temperatures of the two zones are the same. The shock models predict S⁺ zone to have a lower temperature due to the reason stated above. The photoionization model grid is located slightly above the line which means that the S⁺ zone has slightly higher temperatures than the O⁺ zone because fewer coolants are available in deeper regions and ionizing radiation becomes harder. Nevertheless, none of the models predict a temperature of S⁺ zone to be a lot higher than the O⁺ zone. Therefore the extinction measured from the Balmer decrement is likely inapplicable for the correction of O⁺ and S⁺ lines. To get rid of the effect of dust, we constructed a dust-insensitive temperature indicator combining the auroral-to-strong line ratios of [S II] and [O II]. This will be discussed in Sect. 5.

Fig. 11. log ([S II] λλ4068, 4076/ [S II] λλ6716, 6731) vs. log ([O II] λλ7320, 7330/ [O II] λλ3727, 3729). All the data points are plotted with error bars indicating 1σ uncertainty. The blue grid is the photoionization model and the other three grids are the shock models. The indigo grid has a pre-shock density of n = 10 cm−3 and solar metallicity. The red grid has a pre-shock density of n = 1 cm−3 and solar metallicity. The green grid has a pre-shock density of n = 1 cm−3 and sub-solar LMC metallicity. The black arrow represents the effect of 1 mag extinction in AV. Extinction correction using the coefficient predicted by the Balmer decrement will move the data points to the locations indicated by the white data points, which cannot be explained by any model and current understanding of ISM. Therefore, the correction may not be physical. The gray line represents the line ratios when the temperature measured by [S II] is the same as [O II].

5. Discussion

In this section, we compare the measured temperatures to the predictions of models by the two possible ionization mechanisms of LIERs, photoionization from hot, evolved stars and fast radiative shocks. In Sect. 5.1 we describe the models we are using in this study. The measured temperature-sensitive line ratios will be compared to the theoretical values predicted by the models in Sect. 5.2.

5.1. Shocks and photoionization models

1. The photoionization models: The models were computed with CLOUDY c17.02 (Ferland et al. 2017). The ionizing SED was generated by Python-FSPS (Flexible Stellar Population Synthesis) (Conroy et al. 2009; Conroy & Gunn 2010), assuming a simple stellar population (SSP) from a starburst of 5 Gyr ago. At this age, ionizing photons are primarily contributed by post-asymptotic giant branch (pAGB) stars. We used a Kroupa IMF (Kroupa 2001; Kroupa & Weidner 2003) and MIST isochrones (Choi et al. 2016; Dotter 2016) for the SSP. Each set of SSPs contains 6 different stellar metallicities: log(Z/Z_⊙) = –0.6, –0.4, –0.2, 0.0, 0.2, 0.4. The solar abundance set we used is taken from Grevesse et al. (2010), where 12 + log(O/H) = 8.69. We matched the gas-phase metallicities with the stellar metallicities when computing photoionization models. In reality, the stellar metallicity does not need to be close to the gas-phase metallicity in these early-type galaxies. Regardless, the dependence of the shape of the ionizing spectrum on the stellar metallicity is small if pAGB stars are the dominant population (Byler et al. 2019). For each photoionization model, there are 4 different ionization parameters: log(U) = –4.5, –3.5, –2.5, –1.5. The ionization parameter U is defined as the ionizing flux divided by the gas density. For the elements having secondary origins including carbon and nitrogen, we assumed their relative abundances to oxygen change with O/H following the prescriptions given by Groves et al. (2004). We chose not to use the relation given by Dopita et al. (2013) as such photoionization models are not able to cover the data points in Fig. 5. This is consistent with Ji & Yan (2020) that only the prescription of Groves et al. (2004) is able to cover the whole AGN/LIER region in the [N II]/Hα BPT diagram. Finally, we used the default dust depletion factors of CLOUDY v17.02 to remove part of the elements in the gas phase that condense onto dust grains (Cowie & Songaila 1986; Jenkins 1987).

2. The shock models: We used the shock models from the 3MdBs database by Alarie & Morisset (2019), which contains models generated using version 5.1.13 of MAPPINGS V (Sutherland et al. 2018; Sutherland & Dopita 2017). The fast radiative shock models in the database cover a large range of pre-shock densities n, magnetic fields B and shock velocities v. Using different combinations of B and n, different values of magnetic field parameter (B/n^1/2) can be covered. We only included shocks without precursors grid since the amount of gas beyond the shock is insufficient to form an ionized precursor in LIERs (Dopita & Sutherland 1995). The typical compression factor after the shock passed through the gas is 4 and the magnetic field might reduce the number by providing pressure support. When the gas started to cool, the gas was further compressed and the density increases. It is hard to calculate the exact compression factor because it depends on the magnetic field, the shock velocity, the pre-shock density, and the temperature. However, in the shock models we consider, the compression factor is between 100 and 1000 after the cooling phase is finished. Compared to our measured density, which is about 50 cm⁻³, the model grids of n = 1 cm⁻³ and n = 10 cm⁻³ might be able to represent the gas in different stages of cooling and we computed the ratio of [S II] doublets for checking. We included both of the model grids with solar metallicities. Variation in metallicities may also affect the line ratio, and therefore we also plotted a grid with sub-solar LMC metallicity and n = 1 cm⁻³ for reference. For each of these grids, there are 4 different shock velocities: V = 200, 300, 400 and 500 km s⁻¹, and 5 magnetic parameters: B/n^1/2 = 0, 1, 2, 4, 5 μG cm^3/2. Since there are no B/n^1/2 = 0 μG cm^3/2 models, we chose the weakest magnetic field model available. For the two n = 1 cm⁻³ models, the weakest magnetic field available gives B/n^1/2 = 10⁻⁴ μG cm^3/2. For the n = 10 cm⁻³ model, B/n^1/2 = 3.16 × 10⁻⁴ μG cm^3/2 is the weakest option. Dust depletion is not considered in the shock models since fast shocks will destroy dust by grain-grain collisions (Allen et al. 2008). Even if there is remaining dust depletion, our conclusion will not be largely affected since we are looking at ratios of lines from the same ions, except [N II]/Hα. To check the presence of dust depletion, we compare the flux of [Ca II] λ7291 with the predictions from the shock models. The shock models indicate the strength of [Ca II] λ7291 should be at least 10–50% of Hβ depending on the shock velocities, magnetic parameter and metallicity. Nevertheless, we did not detect the line in our spectrum. The 2σ upper limits are 9.5%, 7.5%, and 3.8% of their respective Hβ flux in the high, mid, and low [N II]/Hα bins. With the fact that the strength of [Ca II] λ7291 is predicted to be stronger in lower metallicity, we concluded that dust depletion is present in the gas.

5.2. Data-model comparison

In Figs. 11–13, we included the model grids mentioned in the last subsection, the measured line ratios from this study and those from Yan (2018a). We note that the data from Yan (2018a) is based on the SDSS Legacy. They are not observed using IFUs, but are done using single fiber centered on the centers of galaxies. The model used in this study is also not the same as that used by Yan (2018a). In Yan (2018a), the shock models applied are from MAPPINGS III and in this paper, models from MAPPINGS V are used. Two main differences exist between the photoionization models in this paper and those used by Yan (2018a). First, the ionizing spectra are slightly different. In Yan (2018a), the photoionization model was generated using a 13-Gyr-old SSP from Bruzual & Charlot (2003). In such a model, the isochrone is computed by the Padova 1994 stellar evolutionary tracks (Alongi et al. 1993; Bressan et al. 1993; Fagotto et al. 1994a,b; Girardi et al. 1996), the STELIB/BaSeL 3.1 spectral library, and the Chabrier Initial Mass Function. The STELIB/BaSeL 3.1 uses STELIB spectra (Le Borgne et al. 2003), which cover the wavelength range of 3200–9500 Å and are extended into infrared and ultraviolet using the color–temperature calibrations of BaSeL 3.1 ‘WLBC 99’ SED library (Westera et al. 2002). In this study, we assume a 5-Gyr-old SSP from FSPS using the MIST isochrone (Choi et al. 2016; Dotter 2016). Second, the N/O ratio in Yan (2018a) follows the relation from Vila-Costas & Edmunds (1993), and this study uses the relations in Groves et al. (2004) instead, which better describes the AGN/LIER spaxels in MaNGA (Ji & Yan 2022).

Fig. 12. log ([O III] λ4363/ [O III] λ5007) vs. SOT. The error bars, the grids, and the black arrow are the same as those in Fig. 11. The dust extinction corrected data points are in white. These corrections could be significant overcorrections for the low-ionization lines (see text). With a given value of SOT, the photoionization model prefers a lower O2+ zone temperature than the shock models. For undetected lines, the upper limits are denoted by arrows and are calculated using the 2σ flux values.

Fig. 13. log ([N II] λ5755/ [N II] λ6583) vs. SOT. The error bars, the grids, and the black arrow are the same as those in Fig. 11. The dust extinction corrected data points are in white. These corrections could be significant overcorrections for the low-ionization lines (see text). For undetected lines, the upper limits are denoted by arrows and are calculated using the 2σ flux values.

Figure 11 shows the temperature-sensitive line ratios for [S II] and [O II]. As mentioned in Sect. 4.3, the partially ionized zone of singly ionized sulfur extends deeper toward the neutral zone and should have a similar, or lower temperature than the O⁺ zone. This is predicted by both the shocks and photoionization models, which only cover the central and the lower right regions of the figure. The extinction coefficient predicted by the Balmer decrement is around A_V ∼ 0.7–0.9. If we apply extinction correction using these A_V values, all data points measured would shift to the ‘unphysical’ upper left region, as shown by the data points in white. The gray line shows the line ratios where the temperature measured by [S II] is the same as [O II]. We computed the line by PYNEB assuming the density to be 50 cm⁻³. Varying the density within 1–100 cm⁻³ will induce a small shift to the line and will not change the conclusion.

All of our data are located near the n = 1 cm⁻³ shock model with a shock velocity of around 500 km s⁻¹ and magnetic parameter between 0 and 1 μG cm^3/2. Looking at the position of other shock model grids, both an increase in pre-shock density and a decrease in metallicity will cause the models to move upward. Therefore, it is also possible that our measurements do not indicate such high-velocity shocks, but an environment with sub-solar metallicity or higher density, or a combination of both. Our data are generally consistent with the increasing S⁺ temperature with decreasing metallicity of shock models. Looking at the SDSS-I data, the high [N II]/Hα data is on the n = 10 cm⁻³ grid with a pre-shock velocity between 200–300 km s⁻¹ and magnetic parameter between 2 and 3 μG cm^3/2. It is reasonable as the post-shock density measured for that bin was 136 ± 14 cm⁻³, which is significantly higher than any other bins. The mid [N II]/Hα SDSS-I data is not on any of the grids. That set of data has lower metallicity and density compared to the high [N II]/Hα data. While a decrease in metallicity may cause the point to move upward, it is not expected for these galaxies to have metallicity lower than the LMC. In contrast, the low [N II]/Hα Sloan data is consistent with the photoionization model. It is important to note that we do not account for any dust extinction effect in the above discussion. Since it is expected that the line ratios are affected by a certain degree of dust extinction, we cannot rule out the photoionization model from this plot alone.

While both [S II] and [O II] line ratios are sensitive to dust, a dust-insensitive temperature indicator, SOT, can be constructed using the combination of the two line ratios. This was proposed by Yan (2018b). Since the auroral [O II] lines are on the red side and its strong lines are on the blue side, introducing dust extinction will increase the ratio. On the contrary, the [S II] auroral lines are on the blue side and strong lines are on the red side. As a result, we can combine the two line ratios to make a dust-insensitive temperature indicator. The definition of SOT in this paper is

(8)

We note that the definition is slightly different from the original paper. We put the weak lines in the numerator to be consistent with the signs in other plots. The constant 1.3 is to account for the difference in wavelength separation and has negligible dependence on the choice of extinction curve.

A plot of log ([O III] λ4363/[O III] λ5007) versus SOT is given in Fig. 12. From the extinction vector, it can be seen that SOT is insensitive to dust, any dust extinction will only move the data points in the vertical direction. The high and mid-[N II]/Hα data are consistent with the LMC shock models with magnetic parameters between 2 and 4 μG cm^3/2. The high-metallicity data has a pre-shock velocity of 300–400 km s⁻¹ and the mid metallicity data has one between 400 to 500 km s⁻¹. However, quiescent galaxies are not likely to have metallicities similar to the LMC. We note that SOT shows a degeneracy between pre-shock density and metallicity. Thus, we suspect the positions of the data could be consistent with a pre-shock density higher than 1 cm⁻³ and a metallicity close to solar value. The SOT is sensitive to both metallicity and pre-shock density, and the actual pre-shock density values of these data might be between 1 cm⁻³ and 10 cm⁻³ which cause the data points to fall between the two models. The data of these two bins provide strong evidence for the interstellar shocks because the photoionization model cannot explain such high temperature in the O²⁺ zone. The [O III] measurement of low metallicity gives a 2σ upper limit. It is hard to distinguish whether it is consistent with the shock, or photoionization model.

Another point to note is how the locations of the data vary with metallicity. Both the shock models and photoionization models predict the low [N II]/Hα bin will be shifted to the upper right with respect to higher [N II]/Hα bins. But our data show the opposite trend. None of the models can explain the positive correlation between [O III] temperature and metallicity as we mentioned in Sect. 4.3. On the other hand, for the shock models, increasing the magnetic parameters and shock velocities matches the tendency of the variation. It is possible that these variables might be correlated with the metallicity in the shock mechanism, or related with the increase of [N II]/Hα ratio.

Considering the SDSS-I data points, the high-metallicity data lies on the 10 cm⁻³ data grid which has solar metallicity. It has a B/n^1/2 ≈ 4 μG cm^3/2 and pre-shock velocity between 300 and 400 km s⁻¹. [O III] λ4363 is not detected in mid and low [N II]/Hα bins. Interestingly, the low-metallicity data is on the left of the other two bins. This contradicts what both the shock and photoionization models predict: when metallicity decreases, the data will go to the right. Nevertheless, this is suspected to be also the effect of pre-shock density. In the measurement of Yan (2018a), the post-shock density of the low [N II]/Hα bin is significantly lower than the remaining bins, which is 34 ± 15 cm⁻³ compared to 136 ± 14 cm⁻³ and 71 ± 12 cm⁻³ for the high and mid bins, respectively. As a result, the position of the low [N II]/Hα bin data may be affected by density.

Figure 13 shows the data in log ([N II] λ5755/[N II] λ6583) versus SOT plot. Although none of our data can detect [N II] λ5755, the upper limits of the three bins all show that the photoionization model predicts a higher [N II] temperature than the measured values. The trend of the upper limit shows a similar relationship to the model predictions. Compared to our data, the SDSS-I data show a larger variation. Since the x-axis is the same as Fig. 12, we believe the horizontal position difference between the low [N II]/Hα bin data and the rest of the bins is caused by the difference in pre-shock density. While the high [N II]/Hα data favors the shock models, SDSS-I data of the mid and low [N II]/Hα bins are not covered by any of the grids. It is possible that the position of SDSS-I mid-metallicity data contains both the effects of pre-shock density and metallicity. We do not have a model grid with n = 10 cm⁻³ and sub-solar metallicity which might be able to cover the data point. However, it is important to note again that the SDSS quiescent galaxies are not likely to have a metallicity as low as the LMC. Last but not least, the SDSS-I low [N II]/Hα data shows a higher [N II] ratio than all model predictions and the reason remains unknown.

Finally, SOT versus log([N II] λ6583/Hα) is given in Fig. 14. Both of the line ratios are dust insensitive and the dust extinction vector is too small to be shown. Since both [N II] λ6583 and Hα are robustly detected, we did not include the error bars in the x direction which are smaller than the marker. This figure provides a shred of evidence to rule out the photoionization model which predicts a smaller SOT than the measured values. Our measured data are consistent with the solar metallicity shock model with n = 1 cm⁻³. The high, mid, and low [N II]/Hα points have B/n^1/2 between 0 to 1 μG cm^3/2 and pre-shock velocity of 500, 400, and 300 km s⁻¹ respectively. However, the position of the data points on the grid depends not only on these two parameters but also on metallicity. By comparing the LMC metallicity grid and the solar metallicity grid, a smaller metallicity will shift the model to the left as predicted by the [N II]/Hα. It is worth mentioning that from the model grids in the figure, [N II]/Hα does not solely trace metallicity. In shock models, [N II]/Hα increases with shock velocity while in photoionization models, it is sensitive to the ionization parameter. Therefore, the binning of [N II]/Hα might not efficiently divide the data based on metallicity. Without accurate measurement of metallicity, we cannot constrain the model grid to be used, and the actual shock velocity and magnetic parameter cannot be tightly constrained.

Fig. 14. log ([N II] λ6583/ Hα) vs. SOT. Both values are insensitive to dust, and therefore the extinction vector is not drawn. The error bars on the x-axis are also neglected as they are smaller than the marker and the error bars indicate 1σ uncertainties. All the grids are the same as those in Fig. 11. The blue photoionization grid predicts a lower SOT value compared to all of the data points. All of the data support the interstellar shock models. However, the actual shock parameters cannot be tightly constrained due to the uncertainty in metallicity and pre-shock density.

From the shock models, SOT is not only sensitive to temperature but also the pre-shock density. Again, the difference in the SOT for the SDSS-I data should be attributed to density variations. The high and mid metallicity SDSS-I data are on the n = 10 cm⁻³ shock model. The high [N II]/Hα Sloan data has a magnetic parameter between 2 to 4 μG cm^3/2 and a shock velocity of 300 km s⁻¹. The mid [N II]/Hα data has a similar magnetic parameter and a slightly lower shock velocity. However, it may also be the effect of metallicity difference and it is not possible to constrain shock parameters of the data without a more precise metallicity and pre-shock density value. The low-metallicity SDSS-I data point shows consistency with the n = 1 solar metallicity model, having a shock velocity of 300 km s⁻¹ and a magnetic parameter between 0 and 1 μG cm^3/2. We emphasize that this data point has the lowest post-shock density measured, which can lower the measured SOT value. Overall, all of the data points favor the shock models.

6. Conclusions

In this paper we constrained the ionization source of LIERs by measuring the temperatures of different ions inside these regions in red non-star-forming galaxies. We selected quiescent galaxies from the spatially resolved MaNGA survey and further selected the spaxels that contain LIERs. We divided the sample into three bins based on [N II]/Hα and stacked them together. For each [N II]/Ha bin, we constructed a stellar continuum based on the gas-poor spaxels. The stellar continuum was then subtracted from the stacked spectra of spaxels with strong gas emissions to obtain spectra with purely gas emission lines. Finally, we measured the temperatures using the auroral-to-strong line flux ratios of different ions and compared them to the models of photoionization and interstellar shocks. A detailed summary of our results follows:

1. The Balmer decrement overestimated the dust extinction effect experienced by low-ionization ion species. Since the first ionization energy of neutral sulfur is lower than neutral oxygen, we expect the S⁺ zone to have a temperature that is similar to or lower than the O⁺ zone. Nevertheless, from the measurements of [S II] λ4068, 4076/[S II] λ6716, 6731, and [O II] λ7320, 7330/[O II] λ3727, 3729, correcting the extinction effect with the extinction coefficient predicted by the Balmer line ratio will lead to a S⁺ zone that is much hotter than the O⁺ zone, which is not explained by any of the models and physical pictures. As a result, we suspect that the dust extinction derived from the Balmer lines should not be applied to other emission lines.

2. The O²⁺ zone temperature increases with metallicity, which contradicts the general expectation. In general, metals are effective coolants and temperatures are expected to have an inversely proportional relationship with metallicity. In this study we divided the sample into three bins based on [N II]/Hα, which is a dust-insensitive metallicity indicator when other parameters are fixed. The temperature derived from [O III] λ4363/[O III] λ5007 in the high, mid, and low [N II]/Hα bins are found to be 3.17 K, 2.19 K, and < 1.54 × 10⁴ K, respectively. Therefore, there could be other unconsidered processes affecting the O²⁺ zone temperature.

3. The measured data are mostly consistent with the interstellar shock models. Most of the derived line ratios are predicted by the interstellar shock model with n = 1 cm ⁻³ and solar metallicity. On the other hand, the high and mid [N II]/Hα data from Yan (2018a) are matched by the shock model with n = 10 cm ⁻³ in some of the cases. The intrinsic shock parameters cannot be tightly constrained since the predicted line ratios depend not only on shock velocity and magnetic parameters, but also on pre-shock density and metallicity. Uncertainties in pre-shock density and metallicity estimation make it more difficult to determine the shock parameters.

We note here that our conclusion is different from that reached by Yan (2018a). Yan (2018a) concluded that neither photoionization nor shock could explain the measured auroral-to-strong line ratios. That conclusion was based on an older shock model. With the new shock model adopted in this paper, the Yan (2018a) measurements are generally consistent with shock models as well.

This work provides strong evidence for interstellar shocks as the ionization source of LIERs. It also shows the power of spatially resolved IFU data. Using such data, smaller-scale structures inside a galaxy can be investigated, giving conclusions different from previous investigations. In the future, continuous development of integral field spectroscopy with deeper observations will help detect weak auroral lines in LIERs. Such measurements can provide even more robust constraints on the current mystery.

Acknowledgments

M.-Y.L.L. and R.Y. acknowledge the support of two grants from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No: CUHK 14303123, CUHK 14302522], and the Direct Grant of CUHK Faculty of Science. R.Y. also acknowledges support from the Hong Kong Global STEM Scholar scheme, by the Hong Kong Jockey Club Charities Trust through the JC STEM Lab of Astronomical Instrumentation. Z.S.L. acknowledges the support from Hong Kong Innovation and Technology Fund through the Research Talent Hub program (GSP028). M.-Y.L.L thanks Ricky Wai Kiu Wong and Ziming Peng for their helpful discussions. This research makes use of data from the SDSS-IV. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

All Tables

All Figures

Fig. 1. Dn4000 vs. absolute magnitude in the r band. The values of Dn4000 are from the MaNGA DAP and the absolute magnitudes in the r band are from the NSA catalogue. The gray region shows the density distribution of all galaxies, while the orange contours are drawn at 10%, 30%, 68%, 95%, and 99% of all galaxies. Galaxies between the red dashed lines are selected as our quiescent galaxy sample. Galaxies with invalid Mr and Dn4000 are not plotted.

In the text

	Fig. 2. Median EW value of the selected line among spaxels with EW lower than different percentile in other strong lines. The horizontal lines denote the EWs corresponding to 35 percentiles, which are the zero points we chose.
In the text

Fig. 3. Different strong-line EWs against Hα EW after zero point correction. The blue cross denotes the origin. The grayscale represents the density of the spaxels, while the orange contours enclose 10%, 68%, 95%, and 99% of the spaxels in each panel. The 10% contour is used to show that most of the data are centered at the zero value. The line in the first subplots represents the cut we applied to exclude contamination from transition objects and dusty star-forming galaxies. All spaxels on the right side of the line are removed from the next step.

In the text

	Fig. 4. Different strong-line EWs against Hα EW after excluding dusty and Seyfert spaxels. The gray regions and contours are generated in the same way as those in Fig. 3.
In the text

Fig. 5. [N II]/[O II] vs. [N II]/Hα for all spaxels in the strong-line sample. Instead of a narrow sequence, there is noticeable scatter in this relation. The two blue dashed lines separate the spaxels into high, mid, and low [N II]/Hα bins. The data is meshed into 100 × 100 bins. The grayscale represents the density of all selected spaxels in each bin, while the orange contours enclose 68%, 95%, and 99% of the spaxels. The blue error bars in the top left corner denote the typical 1σ uncertainty of the line ratios. The extinction vector is given at the top right corner.

In the text

	Fig. 6. Distribution of velocity offsets of all strong-line spaxels in each [N II]/Hα bin. The velocity offset is defined as the difference between the stellar velocity and the gas velocity. There are a total of 25 velocity-offset bins covering the range from –300 km s−1 to 300 km s−1.
In the text

	Fig. 7. Stacked strong-line spectrum (blue) and zero-line spectrum (orange) of the three [N II]/Hα bins. The zero-line spectrum has been corrected for overall shape using B-spline and is vertically shifted downward for visibility.
In the text

	Fig. 8. Continuum-subtracted residual spectrum of the three [N II]/Hα bins. All the lines in the spectrum should originate from gas emissions.
In the text

	Fig. 9. Spectra around the four targeted auroral lines corrected by the fit line derived from the sidebands. The dark gray region is the central region, the light gray regions are the sidebands, and the orange line denotes the zero value. [O III] λ4363 is on the top left, [N II] λ5755 is on the top right, [O II] λλ7320, 7330 is at the bottom left, and [S II] λλ4068, 4076 is at the bottom right.
In the text

	Fig. 10. Dust extinction value predicted by the Balmer decrement under case B approximation and low-density limit at 10 000 K. While Hα, Hβ, Hγ, and Hδ are more consistent, Hε is consistent with the prediction only in the low [N II]/Hα bin. Hζ values are not shown due to strong contamination and stellar continuum systematics.
In the text

Fig. 11. log ([S II] λλ4068, 4076/ [S II] λλ6716, 6731) vs. log ([O II] λλ7320, 7330/ [O II] λλ3727, 3729). All the data points are plotted with error bars indicating 1σ uncertainty. The blue grid is the photoionization model and the other three grids are the shock models. The indigo grid has a pre-shock density of n = 10 cm−3 and solar metallicity. The red grid has a pre-shock density of n = 1 cm−3 and solar metallicity. The green grid has a pre-shock density of n = 1 cm−3 and sub-solar LMC metallicity. The black arrow represents the effect of 1 mag extinction in AV. Extinction correction using the coefficient predicted by the Balmer decrement will move the data points to the locations indicated by the white data points, which cannot be explained by any model and current understanding of ISM. Therefore, the correction may not be physical. The gray line represents the line ratios when the temperature measured by [S II] is the same as [O II].

In the text

Fig. 12. log ([O III] λ4363/ [O III] λ5007) vs. SOT. The error bars, the grids, and the black arrow are the same as those in Fig. 11. The dust extinction corrected data points are in white. These corrections could be significant overcorrections for the low-ionization lines (see text). With a given value of SOT, the photoionization model prefers a lower O2+ zone temperature than the shock models. For undetected lines, the upper limits are denoted by arrows and are calculated using the 2σ flux values.

In the text

	Fig. 13. log ([N II] λ5755/ [N II] λ6583) vs. SOT. The error bars, the grids, and the black arrow are the same as those in Fig. 11. The dust extinction corrected data points are in white. These corrections could be significant overcorrections for the low-ionization lines (see text). For undetected lines, the upper limits are denoted by arrows and are calculated using the 2σ flux values.
In the text

Fig. 14. log ([N II] λ6583/ Hα) vs. SOT. Both values are insensitive to dust, and therefore the extinction vector is not drawn. The error bars on the x-axis are also neglected as they are smaller than the marker and the error bars indicate 1σ uncertainties. All the grids are the same as those in Fig. 11. The blue photoionization grid predicts a lower SOT value compared to all of the data points. All of the data support the interstellar shock models. However, the actual shock parameters cannot be tightly constrained due to the uncertainty in metallicity and pre-shock density.

In the text